A new type of self-similarity is found in the problem of a plane-parallel, ultra-relativistic blastwave, propagating in a power-law density profile of the form ρ z k. Self-similar solutions of the first kind can be found for k < 7∕4 using dimensional considerations. For steeper density gradients with k > 2, second type solutions are obtained by eliminating a singularity from the equations. However, for intermediate power-law indices 7 / 4 < k < 2, the flow does not obey any of the known types of self-similarity. Instead, the solutions belong to a new class in which the self-similar dynamics are dictated by the non-self-similar part of the flow. We obtain an exact solution to the ultra-relativistic fluid equations and find that the non-self-similar flow is described by a relativistic expansion into vacuum, composed of (1) an accelerating piston that contains most of the energy and (2) a leading edge of a fast material that coincides with the interiors of the blastwave and terminates at the shock. The dynamics of the piston itself are self-similar and universal and do not depend on the external medium. The exact solution of the non-self-similar flow is used to solve for the shock in the new class of solutions.

Self-similar solutions offer great mathematical simplicity by reducing a set of partial differential equations (PDEs) to a set of ordinary differential equations (ODEs). They are used to characterize the asymptotic behavior of a physical system when its initial scales become unimportant and can sometimes be described by analytical solutions. Self-similar solutions are traditionally classified into one of two types: when the self-similar properties of the system can be found from dimensional considerations, the solutions are of first type. Second type solutions are obtained when one is required to solve an eigenvalue problem.1,2

A canonical example for self-similarity in hydrodynamics is the strong explosion problem. Consider a powerlaw density profile of the form ρ r k, where r is the distance from the origin and k is a constant. The release of a large amount of energy at r = 0 produces a strong shock wave that scales with time like R t α, where α is a function of k. Dimensional analysis can be used when the energy is a relevant parameter of the self-similar flow, such that α is found by enforcing energy conservation. First type solutions for Newtonian spherically symmetric explosions were originally obtained by Taylor,3 von Neumann,4 and Sedov.5 Waxman and Shvarts6 showed that when the density falls sufficiently fast (k > 3) dimensional considerations give a wrong temporal scaling for the shock as it accelerates and loses causal contact with the bulk of the flow. Second type solutions are then obtained for k > 3.26 by eliminating a singularity from the equations, corresponding to a smooth crossing of the sonic line. Waxman and Shvarts6 concluded that when 3 < k < 3.26 the solutions cannot be described by either of the known types of self-similarity. Gruzinov7 subsequently showed that the solutions in that parameter space belong to a new class (see also Kushnir and Waxman8). The uniqueness of the new class will be explained when discussing its ultra-relativistic analogue.

Blandford and McKee9 obtained self-similar solutions of the first type for ultra-relativistic blastwaves, produced when the explosion energy is much greater than the rest mass of the exploding material. The blastwave is engulfed by a shock wave with a Lorentz factor Γ 1, assumed to follow a temporal scaling characterized by the self-similar index
m = d log Γ 2 d log t .
(1)
The pioneering work of Blandford and McKee9 was followed by many other studies of first and second type solutions in various geometries (e.g., Best and Sari,10 Perna and Vietri,11 Hidalgo and Mendoza,12 and Sari13).
This work is focused on ultra-relativistic blastwaves in plane-parallel geometry. We consider an ultra-relativistic shock traveling in the powerlaw density profile given by
ρ z k ,
(2)
where z is the planar spatial coordinate. By studying the solutions at large distances behind the shock, we find that the pressure and Lorentz factor asymptotically approach power-laws in time and position. Three different types of self-similar solutions are then inferred from the asymptotic flow. When the energy in the asymptotic solution converges, the solution is of first type (k < 7∕4). If the asymptotic solution, when applied to the entire space, has diverging energy, the dynamics cannot be inferred from energy conservation considerations. We will show that two different classes of solutions can describe such flows; the first is the well-known second type solutions, in which the self-similar part of the flow is not in causal contact with its non-self-similar part due to the existence of a sonic point (k > 2). The similarity exponent is found by requiring that the solution smoothly crosses the sonic line. In the absence of a singularity, the solution is neither of first nor of second type, defining a solution gap ( 7 / 4 < k < 2). In this region, we find a new class of self-similar solutions, in which the dynamics of the shock are dictated by the rear region of the flow that is not part of the same self-similar solution. A solution for the non-self-similar flow, which behaves as expansion into vacuum, is required in order to solve for the similarity index m. As already mentioned, the new class of solutions was first identified in Newtonian blastwaves by Gruzinov7 and was termed “self-similarity of the third type.” However, the Newtonian case is a trivial manifestation of third type solutions; the shock is driven by a piston of infinite mass and momentum, and simple considerations show that both the piston and the shock must travel at a constant velocity. In this work, we introduce a non-trivial example of the new type of self-similarity in the ultra-relativistic regime.

The paper is organized as follows: after deriving the self-similar equations in Sec. II, we give an overview of first and second type solutions in the ultra-relativistic blastwave problem in Sec. III. We then study the asymptotic behavior of the flow at large distances behind the shock in Sec. IV and obtain exact solutions for expansion into vacuum in Sec. V. The self-similar scaling of the shock in new class of solutions is found in Sec. VI. We compare to numerical simulations in Sec. VII and summarize in Sec. VIII.

This section broadly follows the formalism of Blandford and McKee9 in plane-parallel geometry. We consider a fluid with an ultra-relativistic equation of state, such that its internal energy and pressure satisfy the relation p = e / 3. As a consequence, the velocity and pressure of the flow are decoupled from the material density, simplifying the problem to the conservation of energy and momentum alone
(3)
t γ 2 ( e + β 2 p ) + z γ 2 β ( e + p ) = 0 ,
(3a)
t γ 2 β ( e + p ) + z γ 2 β 2 ( e + p ) + p z = 0 ,
(3b)
where γ is the Lorentz factor of the flow and β = 1 γ 2 is its associated reduced velocity. Equation set (3) can be rewritten in a more compact form by making the transformation to a new variable,
x t z .
(4)
Keeping only leading orders of 1 / γ 2 and defining
q γ 2 ,
(5)
the conservation equations become
(6)
4 q p t + p x = 0 ,
(6a)
p t + p / q x = 0 .
(6b)
Consider a flow with a characteristic Lorentz factor Γ and an associated characteristic position R, related by R ̇ = 1 1 / Γ 2 1 1 / ( 2 Γ 2 ), where the second equality assumes the ultra-relativistic limit. In the context of ultra-relativistic blastwaves, R and Γ are taken to be the position and Lorentz factor of a shock wave, respectively, and the typical scale height of the flow in the immediate downstream of the shock is of order R / Γ 2. It is then natural to define the following similarity variable:
χ = 1 + 2 ( m + 1 ) x R / Γ 2 ,
(7)
which places the shock at χ = 1. The shock is assumed to propagate in a density profile defined by Eq (2). The profiles of q and p behind the shock are given by the self-similar functions g ( χ ) and f ( χ ), defined by
(8)
q ( x , t ) = Γ 2 ( t ) 2 g ( χ ) ,
(8a)
p ( x , t ) = P ( t ) f ( χ ) .
(8b)

The temporal scalings of Γ ( t ) and P(t) are given by the ultra-relativistic jump conditions across the shock, such that taking g ( 1 ) = f ( 1 ) = h ( 1 ) = 1 implies d log Γ 2 / d log t = m and d log P / d log t = m k. Substituting the definitions in Eqs. (7)–(8b) into equation set (6), we obtain the following self-similar equations:

(9)
1 g χ d log g d log χ = ( 7 m + 3 k ) m g χ ( m + 1 ) ( g 2 χ 2 8 g χ + 4 ) η g ( g χ ) ,
(9a)
1 g χ d log f d log χ = 4 ( 2 m + k ) ( m + k ) g χ ( m + 1 ) ( g 2 χ 2 8 g χ + 4 ) η f ( g χ ) .
(9b)
In Sec. III, we derive first and second type solutions to equation set (9).

The results of this section reproduce Sari13's solutions to Eq. (9).

Dimensional considerations imply that the explosion energy is a relevant parameter of the solution and the similarity index m is found by requiring energy conservation. The energy of the shocked fluid has the following scaling:
E Γ 2 ρ R t 1 k m .
(10)
Therefore, for adiabatic blastwaves, m must satisfy
m I = 1 k .
(11)
Solving equations (9) with m = m I, one obtains an analytic solution for the self-similar profiles
g = χ 1 ,
(12a)
f = χ 4 k 7 3 ( 2 k ) .
(12b)
This solution is physical only if the energy in it converges. Integrating the energy density inside the blastwave gives
E = 0 R sh 4 p γ 2 d z 1 f g d χ χ 4 k 7 3 ( 2 k ) | 1 .
(13)
One immediately finds that first type solutions are restricted to k < 7∕4. The pressure gradient behind the shock is decreasing and becomes flat at k = 7/4.
When the density profile becomes sufficiently steep, the shock accelerates rapidly and loses causal contact with the flow behind it. Loss of causal contact happens when the solution crosses a sonic line, defined by the locus of points on which C + characteristics, moving at the speed of sound with respect to the local fluid velocity, stay at a fixed self-similar coordinate (in the same way, a sonic line can be defined for C characteristics that move in the negative z direction with respect to the fluid). Energy conservation can no longer be imposed on the self-similar flow in order to solve for the similarity index m: the dynamics of the shock wave are independent of the part of the flow that contains most of the energy, which lies on the opposite side of the sonic point. Instead, m is found by eliminating the singularity from the equations. Equation set (9) has two singular points χ ± that satisfy
g χ ± = 4 ± 2 3 .
(14)
The positive branch is associated with the sonic line of the negative C characteristic, while the negative branch is the relevant one and defines the sonic line for a C + characteristic. Along g χ + , C + characteristics stay at the same coordinate χ. Barring the formation of discontinuities, a well-behaved solution is required to pass smoothly through the sonic line. This constitutes the eigenvalue problem for m. Eliminating the singularity in equation set (9) requires having
m II = ( 3 2 3 ) k .
(15)

The solutions in this case are implicit in g χ,

(16)
g = [ g χ + 4 + 2 3 2 3 k 3 + 2 3 2 3 k ] ( 2 3 3 ) k ,
(16a)
f = [ g χ + 4 + 2 3 2 3 k 3 + 2 3 2 3 k ] ( 2 3 4 ) k .
(16b)
In contrast to first type solutions, here, the energy is allowed to diverge since most of it is carried by the non-self-similar flow, which is out of causal contact with the shock. One should keep in mind that, however, the eigenvalue problem is defined only if the solution crosses the sonic line. In the absence of a singularity, the solution must be consistent with energy convergence and conservation. In Sec. IV, we will show that second type solutions can be found only for k > 2. One can immediately see that a problem arises when 7 / 4 < k < 2, where the flow does not obey either of the known types of self-similarity. If the energy in this regime diverges, one must first solve for the non-self-similar flow as it dictates the dynamics of the shock.

In order to study energy convergence, we obtain in Sec. IV the asymptotic properties of the flow in the far downstream of the shock.

In this section, we investigate the asymptotic behavior of the blastwave at χ 1 ( x 0) and find that p and q attain power-law profiles in x and t. To acquire some intuition, let us consider a fluid element that is passed by the shock at time t0. Assuming that its Lorentz factor changes like γ fe 2 t m fe at t , the distance between the shock and the fluid parcel at time t is
Δ z = t 0 t ( β sh β fe ) d t [ t 2 ( m fe + 1 ) γ fe 2 t 0 2 ( m + 1 ) Γ 0 2 ] ,
(17)
where Γ 0 is the Lorentz factor of the shock at t = t0. Since fluid elements do not overtake the shock ( m fe > 1), the first term dominates and we find that γ fe 2 t 2 ( m fe + 1 ) x at t t 0, where we approximate Δ z x. Thus, away from the immediate vicinity of the shock, γ is a function only of t and x and is independent of Γ 0 (see also Faran and Sari14) Using Eqs. (7) and (8a) and the approximation for γ fe 2, we can expect to find that g χ constant in the limit χ 1. The asymptotic value of g χ can be found from Eq. (9a) by rewriting it in the following form:
d log g χ d log χ = 1 + η g ( g χ ) · g χ = ( g χ g χ 1 ) ( g χ g χ 2 ) ( g χ g χ + ) ( g χ g χ ) ,
(18)
where
g χ 1 , 2 = 8 3 k + m ± D 2
(19)
are the two possible asymptotic values of g χ, and we define the discriminant
D 16 ( m + 1 ) + ( 8 3 k + m ) 2 .
(20)
Generally, there are two asymptotic values for g χ and the solution settles on the one it encounters first, which depends on the boundary conditions at χ = 1. A single asymptotic value exists if one of g χ 1 , 2 coincides with g χ ±, or if D = 0. The former case corresponds to second-type solutions, where one of the asymptotic values is equal to g χ +. Substituting Eq. (11) or Eq. (15) for m into Eq. (19) returns the following asymptotic solutions:
g χ as = { 9 4 k ± | 7 4 k | 2 , m = m I , 4 3 k 3 | 2 k | , m = m II .
(21)
The limiting behavior of first type solutions can be alternatively derived in the following way: let us write the temporal change in energy through a surface of constant χ:
F ( χ ) = ( w γ 2 p ) d z d t | χ w γ 2 β | z = p ( 1 g χ ) ,
(22)
where w = 4 p is the enthalpy. Integrating the energy conservation equation over a constant self-similar volume implies d F ( χ ) / d χ = d E ( χ ) / d t, where E ( χ ) is the energy contained in the interval d χ. This equality translates to
( 1 g χ ) η f = 1 m k m + 1 .
(23)
Since Eq. (11) enforces ( 1 g χ ) η f = 0, the two possible values for g χ are either g χ = 1, giving zero energy flux, or g χ = 4 ( 2 k ), satisfying constant pressure ( η f = 0). Equation (21) reduces to these two values when taking k 7 / 4. Imposing zero flux at the shock forces the solution to settle on g χ = 1. An analogous interpretation for second type solutions is less straightforward.

Second type solutions exist only when the solution crosses the sonic line, i.e., g χ as < g χ +. Equations (14) and (21) show that this happens only for k > 2, such that g χ as = 4 + 2 3 ( 1 k ).

The asymptotic hydrodynamical profiles are found by substituting Eq. (19) into equation set (9). We obtain power-law profiles with exponents that depend on m and k,
q = Q t x , p t α x δ ,
(24)
where
α = 1 3 ( 4 Q 1 ) , δ = 4 3 ( Q 1 )
(25)
and
Q 1 8 ( m + 1 ) [ m + 8 3 k + D ] .
(26)
It is clear from Eqs. (24) and (25) that the energy diverges at x if δ 0 (Q > 1). We have already shown that the energy in first type solutions converges for k < 7∕4. Second type solutions have Q > 1 and, therefore, the energy in them diverges asymptotically. Nevertheless, this is not a problem since most of the energy is carried by the non-self-similar part of the flow [not described by Eq. (16)], which is causally disconnected from the shock. In reality, the solution beyond the sonic point deviates from that of the self-similar flow such that the energy in the blastwave converges. One can, therefore, remain oblivious to the properties of the inner flow when solving for the dynamics of the shock.

Although m and Q are unknown when 7 / 4 < k < 2, upper and lower limits can be set on m. The asymptotic solution necessarily satisfies g χ as = max { g χ 1 , g χ 2 } and also g χ as > g χ +. These two conditions give m ( 3 2 3 ) k. Another restriction comes from the requirement that D 0, which is satisfied for m 4 3 ( k 1 ) + 3 k or m 4 3 ( k 1 ) + 3 k. Only the first inequality is consistent with a continuous solution for m as a function of k. Therefore, m must lie in within ( 3 2 3 ) k m 4 3 ( k 1 ) + 3 k. These limits imply that inside the gap Q > 1 and the energy diverges asymptotically at x . One, therefore, faces a new challenge: the energy diverges, while a sonic point is not part of the solution. This necessarily means that the self-similar dynamics are dictated by the non-self-similar part of the flow, which contains most of the energy. In order to study the dynamics of the shock, one must first solve for the non-self-similar flow.

A clue to how one should address the problem comes from the fact that when the pressure increases behind the shock, the asymptotic solution satisfies p ( 0 ) = 0 and q ( 0 ) . This kind of behavior is expected in the leading edge of a relativistically hot gas that expands into vacuum. We, therefore, expect the non-self-similar flow to follow expansion into vacuum and look for exact solutions that coincide with the power law asymptotic of Eq. (24) at x 0. We impose initial conditions in which the energy in the flow converges, and then check what kind of shocks can be driven by these flows. The similarity index m will then be found from the condition for the coexistence of the shock with the fluid that expands into vacuum behind it. The associated self-similar solutions will naturally be consistent with energy conservation and convergence.

In this section, we obtain exact solutions describing the expansion of an ultra-relativistic gas into vacuum. Our goal is to find solutions in which the flow approaches the powerlaw asymptotic [(Eq. (24)] expected to appear at x 0, while deviating from it at some large x in order to comply with energy convergence.

In the first stage, we reduce the two fluid equations (6) to a single PDE by performing a hodograph transformation. We proceed by finding solutions that approach the powerlaw asymptotic at x 0, and in the final stage, construct a solution that satisfies specific boundary conditions imposed by initial profiles with finite energy.

The fluid equations (6) can be further simplified via a hodograph transformation. We change independent variables from {t, x} to {u, v}
u = 3 12 ln [ q p 3 / 2 ] ,
(27a)
v = 3 12 ln [ q p 3 / 2 ] ,
(27b)
where u and v are the Riemann invariants of the plane-parallel fluid equations in the ultra-relativistic limit (see, e.g., Johnson and McKee15) conserved along C + and C characteristics, respectively. The numerical prefactor 3 / 12 is introduced to simplify the equations obtained below. The pressure and Lorentz factor can be represented in terms of the new variables by
p = e 4 ( u + v ) ,
(28)
q = e 2 3 ( u v ) .
(29)
The parameters u and v are defined in the entire space < u , v < and correspond to 0 < p , q < . We make an additional transformation to a new dependent variable ψ ( u , v ),
ψ = ( q p ) 1 / 2 Φ = e u ( 2 3 ) v ( 2 + 3 ) Φ ,
(30)
where Φ ( u , v ) is defined by its partial derivatives
t = p · Φ p = 1 8 ( u + v ) Φ
(31)
and
x = p q Φ p + 3 4 Φ q = 1 16 q [ ( 2 + 3 ) u + ( 2 3 ) v ] Φ .
(32)
The above-mentioned variable transformation brings equation set (6) into the form of a remarkably simple Klein–Gordon equation
ψ u v = ψ .
(33)
Next, we find solutions to Eq. (33) that approach the power law asymptotic of Eq. (24).
In this section, we find exact solutions to Eq. (33) that will subsequently be used to describe expansion into vacuum. It is easy to confirm that
ψ pl ( u λ , v / λ ) = A · e λ u + v / λ
(34)
is a special solution of Eq. (33) that returns the power law asymptotic of Eq. (24), where the characteristic parameter λ determines the powerlaw indices α and δ through its relation to Q,
λ = 1 + 2 Q ( 2 + 3 ) 2 + 3 2 Q ,
(35)
and A is a constant that only affects the pressure amplitude.
Equation (33) is also solved by the function ψ = F ( u λ , v / λ ), defined as
F ( u , v ) = n = 1 u n n ! k = 0 n 1 v k k ! .
(36)
In the limit u 1 and v 1 , F ( u λ , v / λ ) simplifies to
F ( u λ , v / λ ) { e u λ + v / λ + I 0 ( 2 u v ) v / ( u λ 2 ) 1 , u v / λ 2 , I 0 ( 2 u v ) v / ( u λ 2 ) 1 , u v / λ 2 ,
(37)
where I 0 ( 2 u v ) e 2 u v 4 π ( u v ) 1 / 4 is the modified Bessel function I n ( 2 u v ) with n = 0, in the limit u v 1. Since e 2 u v e u λ + v / λ when u v / λ 2 , F ( u λ , v / λ ) coincides with the powerlaw asymptotic solution [Eq. (34)] in that limit.

The problem defined by Eq. (33) with boundary conditions on C ± characteristics (corresponding to fixed values of u and v) is a well posed problem, which can be translated to an initial condition problem in p and q. The solution to Eq. (33) is then a superposition of ψ pl and F ( u λ , v / λ ) [and any other additional solution of Eq. (33)] that satisfies the boundary conditions. We obtain the solution that develops from specific initial profiles in Subsection V C.

The results obtained in Secs. V A and V B will now be used to achieve our stated goal: obtaining a solution that contains finite energy and asymptotes to the powerlaw profiles of Eq. (24) at x 0. For that, we impose initial conditions in which the energy converges: in the region 0 < x < x 0, the pressure is an increasing power-law of x, represented by ψ pl ( u λ , v / λ ) with λ > 2 + 3. At x > x 0, the pressure is a decreasing power-law of x, also represented by ψ pl ( u μ , v / μ ) with a different index μ < 2 + 3. This profile ensures energy convergence at x . The two power-law solutions intersect at the point x0, at which the Riemann invariants are chosen (without loss of generality) to be u = v = 0, such that q ( x 0 ) = p ( x 0 ) = 1. The initial profiles at time t0 are, therefore,
p ( t 0 ) = { ( x / x 0 ) δ ( λ ) , x < x 0 ( x / x 0 ) δ ( μ ) , x > x 0 , q ( t 0 ) = x 0 / x ,
(38)
where δ ( λ ) > 0 and δ ( μ ) < 0. In terms of ψ ( u , v ), the initial conditions are given by Eq. (34),
ψ 0 = { A λ · e u λ + v / λ + Δ ψ λ , u < 0 , v > 0 , A μ · e u μ + v / μ + Δ ψ μ , u > 0 , v < 0 ,
(39)
where the terms Δ ψ λ and Δ ψ μ are introduced to account for the offsets x Δ x and t Δ t, which together with A λ and A μ are required for constructing the initial profiles of p and q (see the  Appendix).

The onset of perturbations to the initial profiles occurs upon the arrival of limiting C ± characteristics that emerged from x0 at the initial time t0. The C + characteristic along which u = 0 propagates within the small x power-law asymptotic, while the C characteristic carrying the value v = 0 propagates within the larger-x asymptotic. As illustrated in Fig. 1, the perturbed flow lies in the region between the two characteristics, while outside of it the fluid retains its original profiles. The boundary conditions of Eq. (33) are, therefore, defined on the limiting C ± characteristics that emerge from x0. Equation (33) needs to be solved only in the first quadrant, where u , v > 0; everywhere else the flow is unperturbed and the solution is given by Eq. (39).

FIG. 1.

A schematic characteristic plot illustrating the different regions that develop from the initial conditions of Eq. (38) at times t > t 0. The flow is perturbed between the limiting C + and C characteristics that emerge from x0, while outside this region it retains its exact initial profiles.

FIG. 1.

A schematic characteristic plot illustrating the different regions that develop from the initial conditions of Eq. (38) at times t > t 0. The flow is perturbed between the limiting C + and C characteristics that emerge from x0, while outside this region it retains its exact initial profiles.

Close modal
The problem, thus, admits the boundary conditions ψ ( 0 , v ) = A λ · e v / λ + Δ ψ λ ( 0 , v ) and ψ ( u , 0 ) = A μ · e u μ + Δ ψ μ ( u , 0 ). Although one can find a solution that satisfies the exact boundary conditions, doing so requires adding multiple terms to ψ that are only important for small values of u and v. Recalling that we are interested in the asymptotic evolution of the flow at t t 0 where u , v 1, the boundary conditions can be simplified by taking the limit of u , v , thus neglecting the contribution of the subleading terms Δ ψ λ and Δ ψ μ,
ψ ( u = 0 , v ) A λ · e v / λ , ψ ( u , v = 0 ) A μ · e u μ .
(40)
We, therefore, solve the problem defined by Eqs. (33) and (40), which provide an accurate description of the problem in the limit t . We find the following solution for u , v > 0:
ψ = A λ [ e u λ + v / λ F ( u λ , v / λ ) ] + A μ F ( u μ , v / μ ) , { A λ e u λ + v / λ , u v / λ 2 , A μ e u μ + v / μ , u v / μ 2 .
(41)
Interestingly, the limits taken in the second equality of Eq. (41) show that the same powerlaw asymptotics that were imposed on the initial profiles form independently in parts of the perturbed flow. The boundaries of the powerlaw asymptotics, thus, do not lie along the axes but rather along the lines u = v / λ 2 and u = v / μ 2, whose respective characteristic widths are Δ u λ = 2 v / λ 3 and Δ u μ = 2 v / μ 3. Despite the fact that perturbations from the large x asymptotic are carried into 0 < u < v / λ 2 along C + characteristics, the flow there remains effectively undisturbed. The same is true for the large x asymptotic, where the region 0 < v < u μ 2 remains effectively unchanged to perturbations arriving from the small x asymptotic along C characteristics.
Using Eq. (37) in the limit u v / λ 2, it can be shown that between the two powerlaw asymptotics, the solution is
ψ ( v / μ 2 Δ u μ < u < v / λ 2 + Δ u λ ) e 2 u v 4 π ( u v ) 1 / 4 ( A λ v / ( u λ 2 ) 1 A μ v / ( u μ 2 ) 1 ) .
(42)
We point out that the flow in this regime has negligible dependence on λ and μ, as they do not appear in the exponent. In what follows, we will show that this part of the flow constitutes a piston that contains most of the energy in the solution.

We use Eqs. (31) and (41) with x 0 = 1, λ = 5, and μ = 2.8 to solve for u(v) at time t = 10 80. Figure 2 shows that in the limits u v / λ 2 and u v / μ 2 the exact solution coincides with those of the small and large x asymptotics, respectively. Corresponding profiles of p(x) and q(x) are shown in Fig. 3 using Eqs. (28), (29), (31), and (32).

FIG. 2.

Exact solution obtained from Eq. (41) (black solid line) at time t = 10 80 using λ = 5 and μ = 2.8. The asymptotic power-law solutions of the small and large x asymptotic are the red dashed and green dashed lines, respectively. The line tracing the pressure maximum is also plotted (gray dotted–dashed line).

FIG. 2.

Exact solution obtained from Eq. (41) (black solid line) at time t = 10 80 using λ = 5 and μ = 2.8. The asymptotic power-law solutions of the small and large x asymptotic are the red dashed and green dashed lines, respectively. The line tracing the pressure maximum is also plotted (gray dotted–dashed line).

Close modal
FIG. 3.

Exact profiles of p and q obtained from Eq. (41) (black solid lines) and asymptotic power-law solution of the small x (red dashed line) and large x (green dashed line) asymptotics at t = 10 80. The solution at small x terminates at u = 0 and is extended into the fourth quadrant, where u < 0 and v > 0 using the unperturbed profiles set by the initial conditions.

FIG. 3.

Exact profiles of p and q obtained from Eq. (41) (black solid lines) and asymptotic power-law solution of the small x (red dashed line) and large x (green dashed line) asymptotics at t = 10 80. The solution at small x terminates at u = 0 and is extended into the fourth quadrant, where u < 0 and v > 0 using the unperturbed profiles set by the initial conditions.

Close modal
We are interested in the properties of the flow in the intermediate region given by Eq. (42). Using Eqs. (30), (31), and (42) in the limits u , v , one obtains a solution for u(v, t) that is independent of λ and μ,
u = ( 2 + 3 ) [ ln ( t ) 1 / 2 v ( 2 + 3 ) ] 2 .
(43)
Since the solution lies between the powerlaw asymptotics of increasing and decreasing pressure, a maximum must be described by Eq. (43). It can be seen from Eq. (28) that a maximum in pressure is described by the requirement d u / d v = 1. Taking the full derivative of Eq. (43) with respect to v at constant t, we have
d u d v [ ( 2 3 ) + v u ] + u v + 2 + 3 = 0 .
(44)
Setting d u / d v = 1 and solving the quadratic equation (keeping in mind that u and v are real and, therefore, only one root is relevant), we find that the line that traces the pressure maximum on the {u, v} plane is
u p = ( 7 4 3 ) v .
(45)
The temporal scaling of v at the peak is found by substituting Eq. (45) into Eq. (43)
v p = ln ( t ) 16 ( 2 3 ) .
(46)
Finally, we find the scaling of q and p at the peak by substituting Eqs. (45) and (46) into Eqs. (28) and (29),
p p = e 16 ( 2 3 ) v = t 1 ,
(47)
q p = e 12 ( 2 3 ) v = t 3 / 4 .
(48)
The flow around the peak contains most of the energy and constitutes a hot, accelerating piston. It is interesting to note that at the lower boundary of the solution gap (k = 7/4), the boundary of the powerlaw asymptotic satisfies v / λ 2 = u p ( v ), which implies that powerlaw asymptotic extends all the way to the peak of the pressure. This result is reassured by the fact that the temporal scaling relations of Eqs. (47) and (48) are identical to those of the shock's for k = 7/4. When k > 7∕4, the profiles transitions from the power-law asymptotic to the intermediate solution before the peak.

The width of the pressure maximum on the {u, v} plane can be estimated by Δ v p = ( p d 2 p / d v 2 ) 1 / 2 | v p, yielding a narrowing relative width of Δ v p / v p = 2 ( 2 3 ) ln ( t ) 1 / 2 and Δ u p / u p = 2 ( 2 + 3 ) ln ( t ) 1 / 2. The two characteristic lines u = v / λ 2 and u p ( v ) are shown in Fig. 4.

FIG. 4.

The two components of expansion into vacuum. The line u = v / λ 2 is the boundary between the power-law asymptotic and the rear solution and the respective width over which the solution changes is shaded in red. The peak of the pressure is traced by the line u = ( 7 4 3 ) v, and the shaded gray region designates the width of the peak.

FIG. 4.

The two components of expansion into vacuum. The line u = v / λ 2 is the boundary between the power-law asymptotic and the rear solution and the respective width over which the solution changes is shaded in red. The peak of the pressure is traced by the line u = ( 7 4 3 ) v, and the shaded gray region designates the width of the peak.

Close modal
The Lorentz factor of the flow can be written in a form similar to q = Q t / x of the power-law solution, while the equivalent of Q is now a function of u and v,
q = t x [ 1 3 16 ( u v ) Φ t ] = t x [ u ( 2 3 ) + v ( 2 + 3 ) + 2 u v 2 ( u + v + 4 u v ) ] .
(49)
We confirm that substituting u = v / λ 2 returns q = Q ( λ ) t / x, as expected. Equation (49) also suggests that the profile of q around the piston should not show considerable deviation from those in the powerlaw asymptotics, as observed in Fig. 3; while t/x depends exponentially on u and v, its prefactor has only linear dependence.
The position of the pressure maximum is found by substituting Eq. (45) into Eq. (49), thus obtaining
x p ( t ) = t q = t 0 q p , i ( t t i ) 1 / 4 ,
(50)
with t i being some reference time at which q p ( t i ) = q p , i. While the relative widths of the pressure maximum in terms of Δ v p / v p and Δ u p / u p decrease with time, the spatial relative width of the peak increases logarithmically
Δ x p x p = 3 / 2 2 ln ( t ) 1 / 2 .
(51)
We confirm that the energy in the solution is conserved by recalling that most of it is carried by the piston, whose energy can be estimated by E q p p p Δ x q t 1 t 3 / 4 t 1 / 4 = const, where the typical scale height of the energy density, Δ x q x p t 1 / 4, is dictated by the profile of q.

Further intuition about the nature of the piston can be gained by a heuristic description of a hot, relativistic blob that undergoes free expansion, while exchanging internal energy with bulk kinetic energy. Let us associate Γ b with the bulk motion of the blob and γ th with the random motion of its particles (thermal energy). Since E Γ b ( γ th m c 2 ) = constant due to energy conservation, we have γ th Γ b 1. Adiabatic expansion implies γ th m c 2 ( Γ b Δ ) 1 / 3, where Δ t / Γ b 2 is the volume of the expanding blob in the lab frame. This gives Γ b 2 t 1 / 2. The scaling of the pressure can be found from energy conservation: E Γ b 2 p ( t / Γ b 2 ) = p · t = const p t 1. The evolution of the pressure is inferred solely from energy conservation and, therefore, agrees with our solution for the piston. However, this analysis returns a different scaling for q, which immediately implies that the expansion of the piston is not adiabatic. We conclude that the piston's entropy is not conserved, which suggests that it is not occupied by the same material throughout its evolution.

In Sec. V, we obtained an exact solution for relativistic gas expanding into vacuum, whose leading front approaches the powerlaw asymptotic expected to form in the far downstream of a strong shock. Our objective is to smoothly connect the rear part of the self-similar blastwave solution with the leading front of expansion into vacuum. In order for both solutions to coexist with one another, one must find the appropriate values of the temporal index m for which the powerlaw asymptotic that forms by the piston at u , v > 0 is not destroyed by the shock.

The above-mentioned condition is more easily interpreted by considering the trajectories of the shock and the boundary of the powerlaw asymptotic on the {u, v} plane. A shock can form self-consistently from expansion into vacuum if d u / d v | sh < 1 / λ 2, so that it never crosses the line u = v / λ 2. Using the ultra-relativistic shock jump conditions, this requirement translates to the following inequality:
d u d v | sh = 2 m + 3 ( m + k ) 2 m + 3 ( m + k ) 1 λ 2 .
(52)
At the same time, as Q (and, therefore, λ) depends on the discriminant D [see Eq. (26)], an additional constraint comes from the requirement D 0. Assuming λ is finite, we find that both conditions are satisfied only if d u / d v | sh = 1 / λ 2, which then gives the desired value of m:
m III = 4 3 ( k 1 ) + 3 k ,
(53)
thereby closing the gap between first and second type solutions. We confirm that both m(k) and m ( k ) are continuous at k = 7/4 and k = 2 and note that m III returns D = 0. A self-similar solution for the blastwave can now be found by solving the self-similar ODEs [Eq. (9)] with m = m III.

One might expect that the interpretation of a shock driven by expansion into vacuum continues to be correct also in the regime of second type solutions, as the external density gradient becomes even steeper. Indeed, in solving the inequality in Eq. (52), we have made the assumption that λ is finite. However, in second type solutions, all C + characteristics originate from the immediate vicinity of the sonic point, and therefore, u = constant in space. This corresponds to taking the limit λ for which Eq. (53) returns m = m II = ( 3 2 3 ) k. The analysis in this paper, thus, applies to all k > 7∕4.

A characteristic plot summarizing the main features of the flow is shown in Fig. 5. A fluid element that crossed the shock at an arbitrary time t0 joins the asymptotic self-similar flow designated by the red region between the lines x sh and x pl. As it is advected away from the shock, it leaves the self-similar flow after crossing x pl and eventually joins the bulk of the fluid at the piston when reaching x p. A pair of C ± characteristics are also plotted to demonstrate that the flow is in causal contact. The position of the shock and the boundary of the powerlaw asymptotic have the same scaling and, therefore, the distance between them increases with time.

FIG. 5.

An illustration of the blastwave's dynamics. The power-law asymptotic region that develops behind the shock is shown in red, and the region between the power-law boundary and the peak of the pressure is in gray. We plot the trajectories of C ± characteristics emerging from the power-law asymptotic behind the shock that demonstrate that the flow is in causal contact. A fluid element, whose position is designated by x fe starts at the shock and joins the piston at later times. The dashed lines assume scalings of the asymptotic power-law solution.

FIG. 5.

An illustration of the blastwave's dynamics. The power-law asymptotic region that develops behind the shock is shown in red, and the region between the power-law boundary and the peak of the pressure is in gray. We plot the trajectories of C ± characteristics emerging from the power-law asymptotic behind the shock that demonstrate that the flow is in causal contact. A fluid element, whose position is designated by x fe starts at the shock and joins the piston at later times. The dashed lines assume scalings of the asymptotic power-law solution.

Close modal
The problem introduced by equation set (6) at x > x s ( t ), with boundary conditions
q ( t , x s ( t ) ) = 1 4 x ̇ s , p ( t , x s ( t ) ) = t k 3 x ̇ s ,
(54)
where x s ( t ) is the shock position, can be solved numerically as follows: introducing the Riemann invariants:
V = q 1 2 p 3 4 , U = q 1 2 p 3 4 .
(55)
Equation (6) can be written as
V ̇ + ( 1 3 2 ) V U V = 0 ,
(56)
U ̇ + ( 1 + 3 2 ) V U U = 0.
(57)
where the overdot and prime represent partial derivative with respect to time and space, respectively. Making the following transformation of the independent variable:
x x x s ( t ) ,
(58)
we have for x > 0,
V ̇ + [ ( 1 3 2 ) V U Z ] V = 0 ,
(59)
U ̇ + [ ( 1 + 3 2 ) V U Z ] U = 0 ,
(60)
with
Z 1 4 V ( t , 0 ) U ( t , 0 ) .
(61)
Solving Eqs. (59) and (60) requires only one boundary condition at x = 0. This is because 1 3 2 < 1 4 < 1 + 3 2 and V propagates into the shock, while U propagates out of the shock. It follows that only the U boundary condition is needed. From the original boundary conditions [Eq. (54)],
U ( t , 0 ) = [ 4 3 t k V ( t , 0 ) 2 3 1 ] 2 3 3 .
(62)
The numerical algorithm is now straightforward. We solve Eqs. (59) and (60) at x > 0, propagating the Riemann invariants to the left or to the right, as dictated by the corresponding velocities. At the shock, at x = 0, we compute V ( t , 0 ) by propagating V from the right, and then, we compute U ( t , 0 ) from the boundary condition (62).

We run the simulation for k = 1.98 with initial conditions at t = 1 where a shock is placed at x s = 1 and the pressure behind it decreases as p x 0.0231. In Fig. 6, we show the profiles of p vs q at t = 10 15 , 10 25 , 10 35, together with the theoretical curve (thick line). The numerical profiles approach the theoretical curve as time increases. At t = 10 35, the numerical shock velocity satisfies m = d log x s ̇ / d log t = 0.9138, whereas the theoretical value given by Eq. (53) is m = 0.9186. While our numerical results seem to agree with theory, the simulation converges very slowly and we are, therefore, not able to claim perfect agreement. We note that the choice of k = 1.98 for this example is arbitrary and that numerical results obtained for other values of k indicate a qualitatively similar behavior.

FIG. 6.

Numerical simulations of a shock propagating in a medium with k = 1.98. The profiles of p vs q normalized to their values at the shock are plotted at times t = 10 15 , 10 25 , 10 35 and the theoretical curve is in black. The numerical profiles approach the theoretical curve, yet exact self-similarity is not reached by t = 10 35 due to the slow convergence of the simulation.

FIG. 6.

Numerical simulations of a shock propagating in a medium with k = 1.98. The profiles of p vs q normalized to their values at the shock are plotted at times t = 10 15 , 10 25 , 10 35 and the theoretical curve is in black. The numerical profiles approach the theoretical curve, yet exact self-similarity is not reached by t = 10 35 due to the slow convergence of the simulation.

Close modal

We study the propagation of a plane-parallel, ultra-relativistic shock down a powerlaw density profile of the form ρ z k. In the parameter space 7 / 4 < k < 2, the blastwave cannot be described by any of the known types of self-similar solutions; global conservation laws do not apply in the self-similar domain, and at the same time, no eigenvalue problem can be defined. Instead, the solution gap is populated by a new class of self-similar solutions. Following, we summarize our main results:

  • We find that within the solution gap the flow obeys self-similarity of the third type. The unique characteristic of this new class of solutions is that the similarity index m is determined by the part of the flow that does not participate in the self-similar solution. This is a consequence of the fact that most of the energy in the blastwave resides in the non-self-similar flow, which is in causal contact with the shock. Solving for the shock, therefore, requires a physical understanding of the solution in the entire space. In this work, we first obtained an exact solution for the non-self-similar flow, described by relativistic expansion into vacuum, and then solved for the similarity index m by requiring that the leading edge the expanding flow is not destroyed by the shock. We find that m = 4 3 ( k 1 ) + 3 k for 7 / 4 < k < 2.

  • The flow that expands into vacuum has two general components: (1) an accelerating piston that contains most of the energy and (2) a leading edge described by a powerlaw asymptotic in t and x. The properties at the pressure maximum within the piston are characterized by universal exponents that are independent of k: p p t 1 , q p t 3 / 4, and x p t 1 / 4. Despite containing most of the energy, the piston does not go through free expansion, so that its entropy is not conserved.

  • Expansion into vacuum drives shock waves that obey second type self-similarity when k > 2. Therefore, the analysis in this paper generally applies whenever k > 7∕4. Nonetheless, if k > 2 solving for the non-self-similar flow is not required to describe the shock and its near downstream; a loss of causal contact between the shock and the bulk of the blastwave means that the solution beyond the sonic point is unable to affect the shock and is, therefore, not important for that purpose.

  • Self-similar solutions are important because they describe the asymptotic behavior of a physical system for a large class of initial conditions. In addition, they also offer the technical benefit of having to solve a set of ODEs instead of a set of PDEs, which makes them a very useful tool. This simplifying feature was not used in this work, where an exact solution to the PDEs was obtained in order to derive the self-similar asymptotic thereafter. Nonetheless, the exact solution we find for expansion into vacuum does not probe the flow directly behind the shock, before it settles on the powerlaw asymptotic. In order to study the small scales behind the shock, one must indeed solve the self-similar equations (9) with the corresponding index m(k).

  • Numerical simulations seem to agree with our analytic results for the shock. However, due to the very slow numerical convergence we are not able to confirm perfect agreement with theory.

  • We note that the original classification to first and second type solutions,1,2 only makes the distinction between problems that can be solved by dimensional considerations versus the need to solve an eigenvalue problem. The definition of second type self-similarity should be refined to include the requirement that the eigenvalue problem is imposed exclusively by the self-similar flow.

  • Plane parallel, ultra-relativistic shock waves of the kind discussed in this work do not have obvious physical applications. The spherically symmetric analogue of this problem, on the other hand, is often discussed in the astrophysical context of gamma-ray burst afterglows. In addition, plane-parallel implosions in which the shock wave is converging instead of diverging (explosions) describe the dynamics of ultra-relativistic shock waves that break out from a stellar edge.16,17 According to Sari,13 both first and second type solutions can be found for spherical explosions when 5 3 / 4 < k < 17 / 4 and for plane parallel implosions when 7 / 4 < k < 1 + 3 / 4, such that there is an overlap rather than a gap. While numerical simulations can be used to break those degeneracies, exact solutions can be found by generalizing this work to spherical symmetry or to plane-parallel implosions. We note that another solution gap exists in spherically symmetric implosions13 and is likely analogous to the plane parallel case.

R.S. was partially supported by ISF, NSF/BSF, and MOS grants.

The authors have no conflicts to disclose.

Tamar Faran: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrei Gruzinov: Formal analysis (equal); Writing – original draft (supporting). Re'em Sari: Conceptualization (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

A temporal offset t t + Δ t and a spatial translation x x + Δ x can be introduced by adding the following term to ψ ( u , v ):
Δ ψ = 4 3 [ u ( 2 3 ) v ( 2 + 3 ) ] e u ( 2 3 ) v ( 2 + 3 ) Δ t + 4 3 e u ( 2 + 3 ) v ( 2 3 ) Δ x .
(A1)
The initial conditions defined by Eq. (38) impose Δ t λ = 0 and Δ x λ = Δ x μ = 0. The parameters Δ t μ , A λ, and A μ in this case are
Δ t μ = ( 1 Q ( λ ) Q ( μ ) ) t 0 ,
(A2a)
A λ = 16 λ x 0 ( 2 + 3 ) λ 2 + 2 λ + 2 3 ,
(A2b)
A μ = 16 μ x 0 ( 2 + 3 ) μ 2 + 2 μ + 2 3 ,
(A2c)
where t 0 = A λ ( 4 + λ + 1 / λ ) / 8 and Q ( λ ) , Q ( μ ) are given by Eq. (35). The exact boundary conditions are then ψ ( 0 , v ) = A λ · e v / λ and ψ ( u , 0 ) = A μ · e u μ + 4 Δ t μ 3 u ( 2 3 ) e u ( 2 3 ), where the subleading term that appears due to Δ t μ vanishes in the limit u .
1.
Y. B.
Zel'dovich
and
Y. P.
Raizer
,
Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena
(Academic Press,
1967
).
2.
G. I.
Barenblatt
, “
Self-similarity—Similarity and intermediate asymptotics
,”
Radiofizika
19
,
902
931
(
1976
).
3.
G.
Taylor
, “
The formation of a blast wave by a very intense explosion. II. the atomic explosion of 1945
,”
Proc. R. Soc. London, Ser. A
201
,
175
186
(
1950
).
4.
J.
von Neumann
,
Blast Waves
(
Los Alamos National Laboratory
,
Los Alamos
,
1947
), Vol.
7
.
5.
L.
Sedov
,
Similarity and Dimensional Methods in Mechanics
(
Academic
,
New York
,
1969
).
6.
E.
Waxman
and
D.
Shvarts
, “
Second-type self-similar solutions to the strong explosion problem
,”
Phys. Fluids A
5
,
1035
1046
(
1993
).
7.
A.
Gruzinov
, “
Self-similarity of the third type in the strong explosion problem
,” arXiv:astro-ph/0303242 (
2003
).
8.
D.
Kushnir
and
E.
Waxman
, “
Closing the gap in the solutions of the strong explosion problem: An expansion of the family of second-type self-similar solutions
,”
Astrophys. J.
723
,
10
19
(
2010
).
9.
R. D.
Blandford
and
C. F.
McKee
, “
Fluid dynamics of relativistic blast waves
,”
Phys. Fluids
19
,
1130
1138
(
1976
).
10.
P.
Best
and
R.
Sari
, “
Second-type self-similar solutions to the ultrarelativistic strong explosion problem
,”
Phys. Fluids
12
,
3029
3035
(
2000
).
11.
R.
Perna
and
M.
Vietri
, “
A self-similar solution for the propagation of a relativistic shock in an exponential atmosphere
,”
Astrophys. J. Lett.
569
,
L47
L50
(
2002
).
12.
J. C.
Hidalgo
and
S.
Mendoza
, “
Self-similar imploding relativistic shock waves
,”
Phys. Fluids
17
,
096101
(
2005
).
13.
R.
Sari
, “
First and second type self-similar solutions of implosions and explosions containing ultrarelativistic shocks
,”
Phys. Fluids
18
,
027106
(
2006
).
14.
T.
Faran
and
R.
Sari
, “
The non-relativistic interiors of ultra-relativistic explosions: Extension to the Blandford-McKee solutions
,”
Phys. Fluids
33
,
026105
(
2021
).
15.
M. H.
Johnson
and
C. F.
McKee
, “
Relativistic hydrodynamics in one dimension
,”
Phys. Rev. D
3
,
858
863
(
1971
).
16.
E.
Nakar
and
R.
Sari
, “
Relativistic shock breakouts—A variety of gamma-ray flares: From low-luminosity gamma-ray bursts to type Ia supernovae
,”
Astrophys. J.
747
,
88
(
2012
)..
17.
T.
Faran
and
R.
Sari
, “
Shock breakout from stellar envelopes: The relativistic limit
,”
Astrophys. J.
943
,
97
(
2023
).