Relativistic Spherical Shocks in Expanding Media

We investigate the propagation of spherically symmetric shocks in relativistic homologously expanding media with density distributions following a power-law profile in their Lorentz factor. That is, $\rho_{ej} \propto t^{-3}\gamma_{e}(R,t)^{-\alpha}$, where $\rho_{ej}$ is the medium proper density, $\gamma_{e}$ is its Lorentz factor, $\alpha>0$ is constant and $t$, $R$ are the time and radius from the center. We find that the shocks behavior can be characterized by their proper velocity, $U'=\Gamma_s'\beta_s'$, where $\Gamma_s'$ is the shock Lorentz factor as measured in the immediate upstream frame and $\beta_s'$ is the corresponding 3-velocity. While generally, we do not expect the shock evolution to be self-similar, for every $\alpha>0$ we find a critical value $U'_c$ for which a self-similar solution with constant $U'$ exists. We then use numerical simulations to investigate the behavior of general shocks. We find that shocks with $U'>U'_c$ have a monotonously growing $U'$, while those with $U'<U'_c$ have a decreasing $U'$ and will eventually die out. Finally, we present an analytic approximation, based on our numerical results, for the evolution of general shocks in the regime where $U'$ is ultra-relativistic.


INTRODUCTION
The strong explosion problem, consisting of a sudden release of energy that drives a blast wave into the surrounding medium has been studied extensively.Special attention has been given to shocks propagating in media with power-law density profiles, where self-similar solutions can be identified.While in most astrophysical scenarios blast waves propagate in stationary media, there are cases in which we expect to find blast waves in expanding ones.These types of shocks can be formed when a source of energy launches the shock into a medium that was energized by an earlier explosion, setting it in motion.For example, such a scenario may arise in a binary neutron star merger, where mass is ejected during the coalescence phase, known as the dynamical ejecta and the compact merger-product injects an additional energy into the ejecta as it accretes the remaining material.Similarly, in various types of supernovae a central engine may form following the initial explosion injecting its energy into the expanding envelope.Indication for the existences of such engines can be found, for example, in cases where the energy measured in supernovae of type IC exceeds the amount of energy that can be formed by Nickel decay (Afsariardchi et al. 2021).
In neutron star mergers, the dynamical ejecta is likely to include a fast precursor that reaches mildly or even ultrarelativistic velocities.Mildly relativistic components were observed in numerical simulations (Bauswein et al. 2013;Hotokezaka et al. 2018;Radice et al. 2018;Ishii et al. 2018;Hotokezaka et al. 2018) though the maximal velocities, masses, and density profiles were hard to determine due to limited resolution.Analytic considerations suggest the possibility of an ultra-relativistic precursor (Kyutoku et al. 2014;Beloborodov et al. 2020).For example, in a model explored by Beloborodov et al. (2020) the precursor can reach high Lorentz factors with a power-law mass distribution ρ ∝ γ −α , where 2 ≲ α ≲ 4. Finally, a mildly relativistic precursor has been suggested as an explanation for the gamma-rays observed in GW170817.
The gamma-rays in this model originate from a mildely relativistic shock, driven by a relativistic jet that is launched following the merger, and breaks out from the ejecta (Kasliwal et al. 2017;Gottlieb et al. 2018;Nakar 2020;Beloborodov et al. 2020).While in the case of a merger, the density profile of the fast outflows is hard to constrain, we can consider instead the density profile created as a result of a shock passing through the edge of a star.Barniol Duran et al. (2015) found (based on shock propagation theory of Johnson &Mc-Kee 1971 andPan &Sari 2006) that the density profile of a product ejecta formed by the passage of a mildly relativistic shock through the stellar edge is ρ ∝ γ −α , with α ≃ 1.6, where ρ is the proper density, and γ is the ejecta Lorentz factor.
Motivated by the case of neutron-stars merger, where both the jet driven shock as well as the relativistic ejecta have angular dependence profiles, we consider here the simpler problem of the propagation of a spherical shock in an expanding medium with a power-law density profile in the Lorentz factor.Namely, the outflow is cold and ballistic with a proper density distribution ρ ∝ γ −α e , α > 0 where γ e is the ejecta Lorentz factor and α > 0 is chosen so that the ejecta energy is convergent.Though our work focuses on spherically symmetric systems, it may be applicable also to systems where over a limited angular range the symmetry is nearly spherical, such as relativistic uncollimated jets, and quasi-spherical shocks (as in the case of a shock driven by the cocoon of a choked jet).This is because relativistic blast waves are causally connected over an angular scale of ∼ 1 γ , so flows that are approximately spherical over a larger angular scale will evolve roughly as part of a spherical blast wave.
Most previous works on relativistic shocks focused on shocks propagating in a static medium, in which case a selfsimilar solution can be obtained.These include Blandford & McKee (1976); Sari (2006) to list a few.One exception is Lyutikov (2017), who considered a double shock system, i.e.; an initial shock wave propagating in a power-law density profile followed by an additional shock or wind, and found an approximate self-similar solution for the second shock structure.While the problem bears some resemblance, it is different in the setup of the ambient medium.Lyutikov (2017) considered the second shock to be very close to the first one, thus it propagates in a downstream density profile that is described by the Blandford & McKee (1976) solution.We, on the other hand, consider a shock propagating in a cold, homologously expanding medium with a power-law density profile.Govreen-Segal et al. (2021) studied the Newtonian analogous to the setup considered here.In their setup the shock was propagating in a homologously expanding medium with a power-law density profile: ρ ∝ v −α .They found that for profiles with α ≲ 8, all shocks decay, i.e. the ratio between the shock velocity and the immediate upstream velocity decreases with time.In steeper density profiles, for every density profile, they found a critical ratio between the shock velocity and immediate upstream velocity, such that in shocks with a ratio larger than the critical value, the ratio grows monotonically with time, while shocks with a smaller ratio, monotonically decay.Separating these two regimes is a self-similar solution, describing a shock with a constant ratio between the shock and the immediate upstream velocities, equal to the critical value.
Similar to the Newtonian case (Govreen-Segal et al. 2021), since there are two velocity scales, we generally do not expect to find a self-similar solution.However, we find that for every density profile with α > 0, there exists a critical ratio between the shock Lorentz factor and the immediate upstream Lorentz factor for which a self-similar solution exists.We then use numerical simulations to study the evolution of general shocks, which are not self-similar.
We proceed as follows: in §2 we define the setup we're considering and give an overview of the solutions and go on to derive the self-similar solutions in §2.1.We then use numerical simulations to study general solutions in §3.Finally, in §4 we conclude.

SHOCK PROPAGATION IN A HOMOLOGOUSLY EXPANDING MEDIUM
Consider a spherically symmetric, relativistic shock wave that propagates through a cold, expanding medium.The expansion of the medium is assumed to be homologous, with a velocity v ej related to the radius r and time t by the equation v ej = r/t for r < ct, where c is the speed of light.We model the density profile of the medium as a power law where ρ ej is the proper density, a ej is a constant, and γ ej = γ ej (r, t) is the Lorentz factor of the ejecta, which we assume to be ultra-relativistic.We focus on cases with α > 0, which correspond to a scenario where most of the ejecta energy is stored in slow material.We assume a relativistic ideal equation-of-state with an adiabatic index γ = 4/3.Our goal is to characterize the evolution of a spherical blast wave in such a medium.Specifically, we aim to find under which conditions such a shock decays and ultimately dies out and under which conditions it grows and crosses an infinite number of mass shells.
We denote by capital letters shock-related properties and by sub-index e the medium properties at the location of the shock front, i.e. at the immediate upstream, measured in the lab frame.Namely, R, V, U, Γ represent the shock radius, absolute values of the 3-velocity and 4-velocity and the Lorentz factor respectively, while v e ≡ R t and γ e represent the medium 3-velocity and Lorentz factor at the immediate upstream respectively.Untagged quantities are measured in the lab, or in the proper frame, according to the regular convention.In addition we denote by capital-tagged letters shock properties measured in the frame of the immediate upstream1 .
Let us examine which parameters may affect the shock evolution.Once the initial conditions are forgotten, in addition to the adiabatic index γ and the density power-law index α, the shock evolution can only depend on the shock radius R, the shock velocity V , the medium density and velocity at the immediate upstream ρ e , v e and on time.Since all parameters evolve as powerlaws, a e cannot be a relevant parameter.In addition, since R and t are connected via the ratio R t , which by definition is v e , we are left with two parameters {v e , V }.Finally, since all velocities are relativistic we replace the 3-velocities with the corresponding Lorentz factors and obtain that the shock evolution depends only on γ, α and on the ratio Γ γ e , or equivalently on , the shock 4-velocity in the immediate upstream frame, appropriate for a case where Γ, γ e ≫ 1.Note that while we assume that the shock and the immediate upstream are ultra-relativistic in the lab frame, the shock may be mildly relativistic or even Newtonian in the upstream frame.
As stated above, the fact that the flow has two velocity scales implies that we generally do not expect the solution to be self-similar.An exception is if U ′ is constant throughout the shock evolution.In such a case, there may be a self-similar solution, i.e.; shocks with a critical value of U ′ c (α) = Const.For density profiles in which such a solution exists, U ′ will either monotonically increase or monotonically decrease, depending on whether the value of U ′ is larger or smaller than the self-similar value, thus defining two qualitatively different types of regimes.The first, termed growing shocks, for which U ′ increases monotonically with time asymptotically approaching infinity, and a second type of decaying shocks, consisting of shocks with U ′ that decreases with time.In all decaying shock the shock Lorentz factor ultimately approaches the local Lorentz factor of the moving medium, viz., Γ → γ e , U ′ → 0, where the shock dies out.This condition differs from the static case, in which a decaying shock decelerates, and its velocity approaches zero in the lab frame.It is worth noting that as viewed in the lab frame, both decaying and growing shocks accelerate.
The shock behaviour in the two regimes, i.e. shocks with U ′ > U ′ c grow while shocks with U ′ < U ′ c decay, renders the self similar solution a repelling one (a bifurcation point of shock solutions).In addition, as the energy in the ejecta increases for smaller α, becoming infinite at α < 0, we expect a minimal value of α below which no self-similar solution exists.i.e.U ′ c → α→α min ∞.Below, we derive a self-similar solu-tion with a constant U ′ , and find that such a solution exists for every α > 0.

Self-similar solutions for U
We seek a self-similar solution where for a given value of α, U ′ = U ′ c , or equivalently Γ/γ e = Const.The requirement that Γ/γ e = Const obeys the scaling Γ 2 (t) = At m , where m is a free parameter and A is a scaling constant.This can be shown through the following argument.Denoting the shock velocity as V , the shock trajectory is given to an order O(Γ −2 ) by (2) Using Eq. ( 2) we can obtain the medium velocity at the shock immediate upstream, and its Lorentz factor (to an order O(Γ −2 )): The density at the shock location is given by From ( 4) it is seen that γ e < Γ implies 0 < m < 1.A convenient choice of a similarity parameter is such that the shock is located at σ(R) = 1 + 1 2γ 2 e → 1 to the order we are working with here.Using the similarity parameter, we may re-write the shock downstream parameters γ, p, ρ in terms of the self-similar parameters and reduce the equations of relativistic fluid dynamics to ordinary differential equations, which can be easily solved.The solution requires boundary conditions, attained by the shock jump conditions.

Characteristics
The velocities of the C± characteristics in a relativistic flow are given by: Now, define σ+(t) to be the value of the self-similar coordinate of the C+ characteristic.We have Substituting Eq. ( 25) we obtain It is seen that for the particular characteristic for which ∆ = 0, the term in the parentheses vanishes and dσ+/dt = 0.That is, the sonic point is located at a fixed self-similar coordinate, as in the Newtonian case.

The Sonic Point
Since at the sonic point ∆ = 0, a smooth crossing of this point requires that ∆1 = ∆2 = 0 at that point as well.This, in turn, fixes m for every value of α.In order to determine the eigenvalue m for a given setup with known values of α and γ, we integrate Eqs. ( 17)-( 19) numerically together with the boundary conditions on the shock front, and search for a value of "m" that allows for a smooth transition through the sonic point.
Figure 1 displays the relationship between U ′ c and α for the family of self-similar solutions.We find that a solution exists for cases with α > 0. As expected, U ′ c approaches infinity (m → 1) as α → 0, and converges to zero as α becomes large.We find that to a very good approximation, U ′ c = 4 α in the range tested above.The drop in U ′ c is notably abrupt, with a value of Γ γ e ≈ 1.5 at α = 5.

NUMERICAL RESULTS
To study the shock evolution in the general, non-self-similar cases, we use numerical simulations.

Simulation Setup
We use the publicly available code GAMMA2 (Ayache et al. 2022) to carry out 1D relativistic hydrodynamic (RHD), spherically symmetric simulations.We set up an initial blastwave that propagates in an expanding medium and follow its evolution to times when it is no longer affected by the initial conditions.Using GAMMA, we are able to properly resolve shocks up to a Lorentz factor of 400.
The initial conditions are set such that the upstream has a density profile of ρ e ∝ (γ e v e ) −α and a velocity profile v ∝ r.As the simulation starts, the shock, which was present in the computational domain begins to propagate in the medium.Once the initial conditions are forgotten, the shock location, time, immediate downstream Lorentz factor, and immediate upstream Lorentz factor are collected from each simulation.The simulations are then grouped according to the value of α.Note that for every value of α, there are several simulations ranging in different values of Γ/γ e (R).The fact that in Fig 3 (which will be discussed in the following subsection) the different simulations form a continuous curve in the d log γ e (R) space, shows that the initial conditions are indeed forgotten and that the simulation is at a high enough resolution to accurately resolve the shock Lorentz factor.The simulation details and initial setup are discussed in more detail in appendix A.

Simulation Results
Fig. 2 shows U ′ (t) for several simulations in a density profile with α = 2.The simulations are plotted from the time where the shock reaches a self-consistent structure, independent of the initial conditions.Each simulation is marked with a different color, and the dashed line mark the value of U ′ c .As expected, shocks with an initial U ′ < U ′ c (U ′ > U ′ c ) have a monotonically increasing (decreasing) U ′ (t) throughout the entire simulation.
Fig. 3 shows the different regimes of shock evolution for different values of α.Each point represent a snapshot from a simulation, where all simulations with the same α are given the same color, according to the color legend.The self-similar solution is marked with an x colored according to the corresponding value of α.It sits on the dU ′ d log γe = 0 line (shown with a black dotted line), which divides the parameter plane into two regimes: The domain below the line corresponds to shocks that decay with time.In this case dU ′ d log γe < 0 and the trajectory evolution of the shock is downwards and to the left (shown with a black arrow).In the region above the self-similar line, dU ′ d log γe > 0 and shocks move upwards and to the right with time as they grow.For α = 1 (red dots), we find that the simulations are in slight disagreement with the self-similar solution, and do not pass through it exactly, but rather are slightly below it.The reason is likely numerical.For smaller values of α, it takes the simulation longer to forget the initial conditions, and within the limited dynamical range allowed by our computational resources it is likely that the shock is still affected by the initial conditions.

Analytic Approximation
Examining Fig. 1, we notice that at least for U ′ ≫ 1, dU ′ d log γ e appears to be linear with respect to U ′ .As this approximation must pass through (U ′ c , 0), our ansatz takes the form: where −a is the intercept.Measuring a from the numerical simulations in this manner is not robust, as the numerical differentiation introduces a lot of noise.We therefore first integrate the relation to find: where c is an arbitrary constant that depends on the initial conditions.We can use this expression to fit a linear relation and find a.In our simulations, we find that a ≃ 0.75 − 0.8 for U ′ ≃ 0 − 2, and a ≃ 0.9 − 0.95 for U ′ ≫ 1.In Fig. 1, we adopt a = 0.92.The analytic approximation is shown in Fig. 3 with thin solid lines.For more general use, it may be useful to re-write (29) as: where U ′ 0 and γ e,0 are the initial values taken at a time when the simulation has evolved to a point where the shock evolution becomes independent of the initial conditions.Note that for U ′ 0 > U ′ c (U ′ 0 < U ′ c ) the shock is indeed growing (decaying).

SUMMARY
In this paper, we study the propagation of a spherically symmetric shock in a relativistic homologously expanding medium with a power-law density gradient (ρej ∝ t −3 γ e (R, t) −α , α > 0).The medium is assumed to consist of an ideal cold gas where the adiabatic index of the shocked gas is γ = 4/3.This index is applicalbe for all relativistic shocks and for Newtonian shocks where the downstream internal energy is dominated by radiation.Note that while the setup is spherical, our solutions are applicable to jets and quasi-spherical shocks if they are relativistic enough so the jet core is not causally connected with the edges.We find that while the shock always accelerates in the lab frame, the shock behavior can be characterized according to the shock four-velocity as measured in the immediate upstream frame -U ′ .Once the initial conditions are forgotten, U ′ may either increase monotonously, corresponding to a growing shock, or decrease, meaning the shock is decaying throughout the evolution.Separating these two regimes there exists an unstable self-similar solution, for which U ′ is constant.That is, for every α > 0 there is a critical value U ′ = U ′ c above which the shock grows and below which it decays.
The self-similar value diverges as α → 0, and decreases with α, so that for α ≳ 5, in the self-similar case, U ′ is mildlyrelativistic or even Newtonian.We present an analytical approximation that can be used to describe general shocks, and seems to be robust as long as U ′ is ultra-relativistic.

Figure 1 .
Figure 1.The phase space of the analytic solution.The black line shows the dependence of the self-similar value of U′ (U ′ c = m 2 √ 1 − m) on α, found by requiring a smooth transition through the sonic point.The white line at α = 0 marks the asymptote at which the self-similar solution diverges and terminates.The black dotted line depicts the relation U ′ = 4/α, which is in good agreement with the semi-analytic solution for U ′ c .The self-similar solution divides the parameter space into two distinct regions; growing shocks, and decaying shocks.

Figure 2 .
Figure 2. The evolution of the shock 4-velocity measured in the immediate upstream frame, U ′ (t), in several simulations with a medium having a powerlaw density profile with α = 2.The date from each simulation is marked with a different color.The black dashed line marks self-similar value U ′ c .Time is measured in arbitrary units.

Figure 3 .
Figure 3.The evolution of the shock strength defined by the parameter dU ′ d log γe , which measures the change in U ′ as it crosses the upstream as a function of U ′ , is plotted for different values of the density profile.The figure shows several simulations for every value of α.The values from the numerical simulations are marked by dots, and the self-similar solutions all lie on the dotted black line at (U ′ c (α) , 0) and are marked with 'x'.Regions above the dotted line correspond to growing shocks, while the ones below it are decaying shocks.The solid lines denote the analytical approximation (see §3.3).
Figure 3.The evolution of the shock strength defined by the parameter dU ′ d log γe , which measures the change in U ′ as it crosses the upstream as a function of U ′ , is plotted for different values of the density profile.The figure shows several simulations for every value of α.The values from the numerical simulations are marked by dots, and the self-similar solutions all lie on the dotted black line at (U ′ c (α) , 0) and are marked with 'x'.Regions above the dotted line correspond to growing shocks, while the ones below it are decaying shocks.The solid lines denote the analytical approximation (see §3.3).