Dynamics of Relativistic Shock Waves Subject to a Strong Radiation Drag: Similarity Solutions and Numerical Simulations
Ilia Leitus Amir Levinson
astro-ph.HE


AIP/123-QED

]Dynamics of Relativistic Shock Waves Subject to a Strong Radiation Drag: Similarity Solutions and Numerical Simulations





 Levinson@wise.tau.ac.il
 Raymond and Beverly Sackler School of Physics & Astronomy, Tel Aviv University, Tel Aviv 69978, Israel



 





             


We examine the effect of Compton drag on the dynamics of a relativistic shock wave. 
Similarity solutions describing a radiation-supported shock are obtained for certain 
profiles of the external radiation intensity and the density of the unshocked ejecta, and
are compared with 1D numerical simulations of a blast wave expanding into an ambient 
medium containing isotropic seed radiation.  Both the analytic model and the simulations indicate 
that under realistic conditions the radiation drag should strongly affect the dynamics and structure of the shocked layer.
 In particular, our calculations show significant strengthening of the reverse shock and weakening of the forward shock 
 over time of the order of the inverse Compton cooling time.   We conclude that the effect of radiation drag on the 
 evolution of the emitting plasma should affect the resultant light curves and, conceivably,
spectra of the observed emission during strong blazar flares.





                             

                              
[
    Amir Levinson
    December 30, 2023
=====================







§  INTRODUCTION


High-energy emission is a defining property of compact relativistic astrophysical sources.   There is ample 
evidence that the emission is produced in relativistic outflows launched by  an accreting
black hole or a magnetar.    The interaction of these outflows with their environments results in formation 
of expanding cocoons and shocks, that decelerate the flow and lead to lower energy emissions over a  large range of scales,
e.g., afterglow emission in GRBs and radio lobes in AGNs and micro-quasars.    In certain circumstances the relativistic jets interact 
also with ambient radiation fields emitted by the surrounding gas.   This interaction can considerably alter the dynamics 
of dissipative fluid shells at relatively small scales,  and is likely to be the origin of the variable gamma-energy emission observed
in certain classes of high-energy sources. 


Direct evidence for such interactions is found in blazars, where the properties of the ambient radiation field can be measured.
This radiation is contributed primarily by  the accretion disk around the black hole,  by gaseous clouds in the broad line 
region (BLR) and, at larger scales, by a dusty molecular torus <cit.>.   The conventional wisdom has
been that  the high-energy emission seen in powerful  blazars 
results from inverse Compton scattering of ambient seed photons by non-thermal electrons accelerated 
in dissipative regions in the jet <cit.>. In most early works the dynamics of the radiating fronts is not computed in a self-consistent manner.  More recent works 
incorporated realistic shock models to compute light curves of flared emission <cit.>.
However, the effect of the radiation drag on the structure and dynamics of the shock has been neglected. 

In this paper we construct a model for the dynamics of a relativistic shock wave in the presence of intense radiation field. 
In section <ref> we derive the governing equations of a radiation-supported shock.
In section  <ref>  we present self-similar solutions, obtained under several simplifying 
assumptions regarding the properties of the ejecta and the intensity profile of the external radiation field, 
and study the dependency of the evolved structure on the magnitude of the radiation drag.  In section  <ref> 
we present results of numerical simulations of a uniform, spherical ejecta propagating through an ambient medium 
containing isotropic seed radiation.  The simulations enable analysis of the temporal evolution from the 
initial stage under more realistic conditions.  While the analysis is motivated by the application to blazars, mainly 
because there we have sufficient information to asses the conditions in the source, it might also be relevant to 
other high-energy astrophysical systems.




§ RADIATION-SUPPORTED SHOCK MODEL


Quite generally, the collision of relativistic fluid shells leads to the formation 
of a double-shocked layer.  Some examples include the interaction of expanding ejecta with 
an ambient medium and collision of fast and slow shells (internal shocks).  
In the absence of external forces, the dynamics of the double-shock
system is dictated solely by the properties of the unshocked media.  In the presence of
external radiation field, the shocked layers are subject to radiation drag that 
can alter the dynamics significantly.  When the drag is strong enough it should
lead to a weakening of the forward shock and strengthening of the reverse shock.  In 
essence, the compression of the shocked ejecta (or the shocked plasma of the fast shell 
in case of internal shocks) is supported by the radiation pressure imposed on it.  
Thus, as a simple approximation one can ignore the entire region beyond the contact discontinuity, and
obtain solutions describing a radiation-supported shock.  This is the treatment adopted in section <ref> below.
The validity of this approximation is checked in section <ref> using numerical simulations.


Below we adopt the following notation:
The velocities of the reverse shock and the contact surface are denoted by V_r and V_c, respectively,
and the corresponding Lorentz factors by Γ_r and Γ_c.  Plasma quantities are denoted by 
lower case letters, and subscripts u and s refer to the unshocked and shocked ejecta, respectively. 





 §.§ Flow equations



The dynamics of a relativistic fluid is governed by the continuity equation,  

    ∂_μ(ρ u^μ)=0,

and energy-momentum equations,

    ∂/∂ x^αT^αμ=S^μ,

here expressed in terms of the energy-momentum tensor

    T^μν=w  u^μ u^ν + g^μν p,

where ρ, p and w are, respectively, the proper density, pressure and specific enthalpy, 
u^μ=γ(1,v_1,v_2,v_3) 
is the fluid 4-velocity, and g^μν is the metric tensor, here taken to be the Minikowsky metric
g_μν= diag(-1,1,1,1).  The source term on the right hand side of Eq. (<ref>) accounts
for energy and momentum exchange between the fluid and the ambient radiation field, and is derived explicitly below.
The above set of equations must be augmented by an equation of state.  
We adopt the relation 

    w=ρ+γ̂/γ̂-1p≡ρ+ap,

with a=4 for a relativistically hot plasma and a=5/2 for sub-relativisric temperatures.

We shall restrict our analysis to spherical (or conical) flows, and adopt a spherical coordinate
system (r,θ,ϕ).  Equations (<ref>) and (<ref>) 
for the shocked fluid then reduce to:

    ∂_t(ργ)+1/r^2∂_r(r^2ργ v)=0,
    ∂_t(wγ^2-p)+1/r^2∂_r(r^2 wγ^2 v)=S^0,
    
    wdlnγ/dt+dp/dt-γ^-2∂_tp=v(S^r-v S^0),

here d/dt=∂_r +v∂_r denotes the convective derivative.

A common approach is to seek solutions of Eqs. (<ref>)-(<ref>) inside the 
shocked layer, given the properties 
of the unshocked flows upstream of the forward and reverse shocks as input.  The jump 
conditions at the forward and reverse shocks then serve as boundary conditions for
the shocked flow equations.  In addition, the fluid velocity and total momentum flux must 
be continuous across the contact discontinuity surface.  In the approximate treatment adopted here, 
whereby the reverse shock is fully supported by radiation and the layer enclosed between the contact and the 
forward shock can be ignored, the location of the contact is fixed by the requirement that the total energy is 
conserved, given formally by Eq. (<ref>) below.



 §.§ Jump conditions at the reverse shock

For the range of conditions envisaged here the shocked plasma is optically thin.  The shock, in this case, is collisionless, i.e., 
mediated by collective plasma processes (as opposed to radiation mediated shocks), and its width, roughly the skin depth, is much smaller than 
the Thomson length.   Thus, the shock transition layer can be treated as a discontinuity in the flow parameters. 
The local jump conditions at the reverse shock can be expressed in terms of the local shock velocity V_r=dR_r/dt, where R_r(t) is the 
shock radius at time t, as <cit.>

    ρ_uγ_u(v_u-V_r)=ρ_sγ_s(v_s-V_r), 
            w_uγ_u^2v_u(v_u-V_r)+p_u=w_sγ_s^2v_s(v_s-V_r)+p_s, 
           w_uγ_u^2(v_u-V_r)+p_uV_r=w_sγ_s^2(v_s-V_r)+p_sV_r.
   
We shall henceforth assume that the unshocked medium is cold and set p_u=0.
The jump conditions at the forward shock can, in principle, be derived in a similar manner, but are redundant under our radiation-supported shock approximation.  




 §.§ Compton drag terms

Behind the shock electrons are heated and accelerated to relativistic energies.  
We denote the energy distribution of electrons there by F_e(γ_e), where 
γ_e=ϵ/m_ec^2 is the dimensionless electron energy, as measured in the fluid frame.  The total electron 
proper density is then given by n_e=∫ F_e(γ_e) dγ_e, and the average energy and second energy moment by 

    <γ_e>=1/n_e∫γ_e  F_e(γ_e) dγ_e,   <γ_e^2>=1/n_e∫γ_e^2  F_e(γ_e) dγ_e.

If only the thermal electrons are taken into account one expects  m_e<γ_e>≃ m_pγ^'_u/2 for plasma 
in rough equipartition, where 
γ^'_u=γ_uΓ_r(1-v_uV_r) is the Lorentz factor of the upstream flow, as measured in the shock frame, and it is assumed that γ^'_u>1.  If the contribution of non-thermal electrons is substantial, then <γ_e> may be larger. 
Since the shocked electron plasma is relativistic we have m_ec^2n_e<γ_e^2>=w<γ_e>.

The derivation of the source terms is given in Ref GL15.   In terms of the energy density of the 
radiation intercepted by the flow, u_rad(r), and the quantities defined above they can be expressed as 

    S^0=-8/3γ_s^3<γ_e^2>u_radσ_T n_e,

and 

    S^r=v_s S^0+S^0/3γ_s^2.

From Eqs. (<ref>) and (<ref>) it is readily seen that the deceleration time 
of the shocked fluid is given by 
t_dec=-3wγ_s^2/v_s S^0=2 γ_s t_IC^', where t^'_IC≃ 3m_ec/(4σ_T u^'_rad<γ_e>)
is the inverse Compton cooling time of an electron having an energy m_ec^2<γ_e>, as measured in the fluid rest frame, and 
u_rad^'=4γ_s^2 u_rad/3.  Substantial deceleration occurs if t_dec is much shorter than the outflow time, t_f=r/c.
For convenience, we define the drag coefficient as

    α=3t_f/t_dec=8/3 m_e c^2γ_s<γ_e> σ_T u_rad r.

As noted above, for thermal electrons γ_s<γ_e>≃ (m_p/2m_e)γ^'_uγ_s
=(m_p/4m_e)γ_u.  If only the contribution of thermal electrons is accounted for, then

    α=8 m_p σ_T/12 m^2_e c^2γ_u u_rad  r.

The energy loss term can now be expressed as:

    S^0=-αγ_s^2 w r^-1.


To get an estimate of the value of α anticipated in blazars, suppose that a fraction η of the nuclear luminosity, L=10^45 L_45
erg/s, is scattered across the jet.  The corresponding energy density of the radiation intercepted by the jet is roughly 
u_rad(r)=η L/4π r^2=3×10^-4η L_45/r_pc^2 ergs 
cm^-3, where r_pc is the radius in parsecs.   From Eq. (<ref>) one obtains α= 0.7 γ_uη  L_45 r_pc^-1. 
A more realistic estimate can be obtained using recent calculations of 
the ambient radiation intensity contributed by the various radiation sources surrounding a blazar jet <cit.>.  
It is found that for a prototypical blazar the intensity profile of radiation intercepted by the jet is flat up to a distance of about one parsec, 
with u_rad≃ 10^-3 erg cm^-3, followed by a gradual decline .   For an internal shock 
produced by colliding shells, Equation (<ref>) yields for the drag coefficient α≈ 3γ_u r_pc.   
Typically γ_u >10, implying significant drag on sub-parsec scales in blazars.




 §.§ Global energy conservation

The change over time in the total energy of the shocked ejecta (or shocked fast shell) must equal the difference between the energy injected through 
the reverse shock and the energy radiated away through the contact.   The energy accumulated behind the forward shock is neglected in our
approximate treatment. 
The time change of the total energy per steradian of the shocked ejecta is

    ∂_t E    =   ∂_t ∫_R_r(t)^R_c(t)T^00r^2 dr = ∫_R_r(t)^R_c(t)∂_t(T^00r^2) dr 
       +   R_c^2 T^00_s(R_c)V_c-R_r^2 T^00_s (R_r)V_s,

where R_r(t) and R_c(t) are, respectively, the radii of the reverse shock and the contact surface, and 
V_r =dR_r/dt, V_c=dR_c/dt are the corresponding velocities. 
Integrating Eq. (<ref>) from the shock, R_r(t), to the contact, R_c(t), one obtains:


    ∫_R_r(t)^R_c(t)∂_t(T^00r^2) dr +R_c^2 w_s(R_c)γ_s^2(R_c) v_s(R_c)
        -    R_r^2 w_s(R_r)γ_s^2(R_r) v_s(R_r)=∫_R_r(t)^R_c(t)S^0 r^2 dr.

Combining Eqs. (<ref>) and (<ref>), using the jump condition  (<ref>) with p_u=0 (cold ejecta) and the relation T^00=w γ^2 -p, and recalling that the fluid velocity is continuous across the contact, viz., v_s(R_c)=V_c, yields

    ∂_t E =R_r^2 ρ_uγ_u^2(v_u-V_r)-R_c^2p_s(R_c)V_c +∫_R_r(t)^R_c(t)S^0 r^2 dr.

The first term on the right hand side accounts for the power incident through the reverse shock, the second term for pdV
work at the contact and the last term for radiative losses. 




§ SELF-SIMILAR SOLUTIONS

Similarity solutions describing the interaction of a relativistic shell with an ambient medium, in the absence 
of radiative losses, 
were derived in Ref NS06, and their stability was subsequently analyzed <cit.>.    Such solutions can be obtained for a
freely expanding ejecta characterized by a velocity profile v_u=r/t at time t after the explosion, and  
a proper density of the form

    ρ_u=a_u/t^3γ_u^n,

where a_u is a normalization constant, and γ_u=1/√(1-v_u^2) is the corresponding Lorentz factor (it can be 
readily seen that the continuity equation is satisfied for this choice of ρ_u and v_u).   
Self-similarity requires that the Lorentz factors of the forward shock, reverse shock and the contact discontinuity have a similar 
time evolution, viz., Γ^2_f=At^-m, Γ^2_r=Bt^-m, Γ^2_c=Ct^-m, where 
A,B,C and m are constants determined upon matching 
the solutions in region 1 (shocked ejecta) and region 2 (shocked ambient medium) at the contact discontinuity<cit.>.  
For an ambient density profile of the form ρ_i∝ r^-k, the index m is given by <cit.>

    m=6-2k/n+2.


In cases where the shocked plasma is subject to strong radiative losses it is still possible to obtain self-similar solutions provided
the energy source term scales as S^0∝γ_s^2 w r^-1 (see Eq. (<ref>) and (<ref>)), that is, 
α in Eq. (<ref>) is constant.   We shall henceforth make this assumption even though it implies a somewhat artificial intensity
profile of the external radiation field.  We note, however, that as long as the deceleration length is smaller than the radius of the shock these
details are unimportant. This is confirmed by numerical simulations, presented in the next section.

As explained above, in the presence of a strong radiation drag the dynamics of the system is dictated by the interaction of the external radiation
field with the shocked plasma enclosed between the contact and the reverse shock.  To a good approximation one can then ignore the
contribution of the forward shock to the overall dynamics.  The evolution of the radiation-supported shock still satisfies  Γ^2_r=Bt^-m, 
Γ^2_c=Ct^-m, but the index m is now determined from global energy conservation, as will be shown below.   Now, 
to order O(Γ_r^-2) the trajectory 
of the reverse shock is given by

    R_r(t)=∫_0^t(1-1/2Γ_r^2) dt^'=t-t/2(m+1)Γ_r^2,

from which we obtain for the velocity of the ejecta crossing the shock: v_u(R_r)=R_r/t=1-1/[2(m+1)Γ_r^2].
The corresponding Lorentz factor is thus given, to the same order, by

    γ_u^2=(m+1)Γ_r^2,

and the density profile by

    ρ_u= a_u/(m+1)^n/2 B^3/mΓ_r^(6/m-n).

We adopt the similarity parameter introduced in Ref BMK76,

    χ=[1+2(m+1)Γ_r^2](1-r/t).

With this choice the reverse shock is located at χ=1, and the contact at χ_c=Γ^2_r/Γ_c^2=B/C<1.
Following Ref. NS06 we define the self-similar variables, G(χ), H(χ) and F(χ), such that at the reverse shock they satisfy the boundary conditions G(1)=H(1)=F(1)=1.  Upon solving the jump conditions (<ref>)-(<ref>), the shocked fluid quantities
can be expressed in terms of the self-similar variables as


    γ_s^2=qΓ_r^2G(χ),
          ρ_sγ_s=m q ρ_uγ_u/(m+1)(q-1)  H(χ),
          p_s=mρ_u/a(q-1)+2(1-√(q/(m+1))) F(χ),
           w_s=ρ_s+ap_s≡ K(χ)p_s,

where the parameter √(q) is the only positive root of the polynomial equation

    γ̂x^3+(2-γ̂)√(m+1) x^2-(2-γ̂)x-γ̂√(m+1) x=0.

Upon substituting these relations into Eqs. (<ref>)-(<ref>), and using Eq. (<ref>) with a constant drag coefficient, α= const [More generally, α can be taken to be any function of the similarity coordinate χ. In the presence of strong drag the width of the shocked layer is much smaller the its radius, hence we anticipate little variations in α across it.], we obtain the following set of equations for the self-similar variables:

    2(1   +   qGχ)∂_χln F-(1-qGχ)K∂_χln G
       =   mn-6-(m-2/3α)K/(m+1)qG, 
     
    2(1   -   qGχ)∂_χln F-γ̂(1+qGχ)∂_χln G
       =   6-mn+(m-4)γ̂-2/3(γ̂-1)
    α K/m+1qG,
    
    2(1   -   qGχ)∂_χln H - 2∂_χln G
       =   -(mn-m-2)/(m+1)qG,

subject to the boundary conditions G(1)=F(1)=H(1)=1.  These equations reduce to those derived in 
Ref. NS06 in the special case α=0.  
As can be seen, at the contact, χ=χ_c, the Lorentz factor must satisfy qG_cχ_c=1, where for short we denote G(χ_c)=G_c. 
This relation defines the limit of integration.  

The self-similar equations involve two eigenvalues; the index m and the location of the contact χ_c.   Thus, two conditions are needed
to find them. The first one is the relation 

    q G_cχ_c=1

mentioned above.   The second one is global energy conservation, Eq. (<ref>).
To order O(γ_s^-2) we have T_s^00=w_sγ_s^2-p_s≃ w_sγ_s^2,  r=t and dr=-t dχ/[2(m+1)Γ_r^2], 
from which we obtain

    E(t)   =   ∫_R_r(t)^R_c(t)T^00r^2dr = t^n m/2 q m a_u (1-√(q/(m+1)))/2(m+1)^n/2+1B^n/2[a(q-1)+2]
       ×   ∫_χ_c^1K(χ)F(χ)G(χ) dχ,

and ∂_t E=nm E/2t.  To the same order one has ∫_R_r^R_cS^0r^2dr = -α E/t.  Substituting these results into Eq. (<ref>) yields the constraint 

    ∫_χ_c^1K(χ)   F   (χ)G(χ) dχ=(m+1)/q(α+nm/2)
       ×   [ a(q-1)+2/1-√(q/(m+1))-2F_c  ],

where F_c=F(χ_c).



 §.§ Results

For a given choice of the drag coefficient α we guess the value of m and integrate Eqs. (<ref>)-(<ref>) from χ=1 to the point χ_c at which Eq. (<ref>) is satisfied.   We then check if Eq. (<ref>) is satisfied.  If not, we change the value of m and repeat the process until
a solution satisfying all constraints is found. 
Sample profiles are shown in Fig. <ref>  for n=1 and different values of α.  The radiation free case (α=0) is shown for a comparison.  It was computed using the full solution of the two-shock model described in Ref. NS06.  As seen, the width of the shocked layer, Δχ=1-χ_c, decreases
with increasing drag coefficient α, as naively expected.   Moreover, larger radiative losses lead to increased non-uniformity  of the
Lorentz factor and pressure in the shocked layer.  The divergence of the density at the contact 
is a basic feature of the similarity solutions <cit.> and occurs even in the absence of radiative losses (α=0), as seen in the upper panel of Fig. <ref>.





Figure <ref> depicts the dependence  of the index m on α.   For a comparison, the value of m obtained in the case α=0 for a blast wave propagating in a uniform density medium (k=0) is m=2 (see Eq. (<ref>) with k=0 and n=1).   

The 4-velocity of the unshocked fluid, as measured in the frame of the reverse shock,  is given by

    u_u^'=γ_uΓ_r(v_u-V_r)=m/2√(m+1),

and it is seen that the shock becomes substantially stronger as α increases.   
Note that the power dissipated behind the shock, ρ_uγ^'_uu^'_u, increases roughly linearly with the index m.
At sufficiently large drag the entire power incident through the shock is radiated away, and the solution becomes independent of α,
as seen in Fig. <ref>.  For our choice of parameters, specifically n=1, this occurs at α 50, for which m≃ 19 and u^'_u≃2.1.    We emphasize that this limit can be approached provided the Lorentz factor γ_u of the unshocked shell is sufficiently large.   
To be more concrete, Eq. (<ref>) implies that γ_u≃ 4.6Γ_r as α→∞, while our analysis is valid only for Γ_r>>1.   The Compton drag terms given in Eqs. (<ref>) and (<ref>) assume that the intensity of ambient radiation is highly 
beamed in the frame of the shocked fluid.   Once γ_s decelerates to modest values, γ_s1, the 
drag force is strongly reduced, ultimately becoming ineffective. 













§ NUMERICAL SIMULATIONS

Numerical simulations were performed using the PLUTO code <cit.>, that was modified to incorporate energy and momentum losses of the shocked plasma through Compton scattering.  To clarify the presentation we adopted an external intensity profile of the form u_s∝ r^-1 for which α in Eq. (<ref>) is constant.   Over the deceleration scale this profile is a reasonable approximation to the flat profile expected in blazars <cit.>.
We start with the basic spherical blast wave problem, releasing an ejecta with a uniform Lorentz factor γ_e = 10 into a stationary ambient medium
(henceforth, quantities in the unshocked ejecta are designated by subscript "e").
In our setup, both the ejecta and the ambient medium are taken to be uniform initially with a density ratio ρ_e/ρ_i = 100. 
The calculations are performed in "simulation units" in which a fluid element traveling with a unit velocity (the speed of light) passes a unit of length per 
unit of time. The initial impact occurs at a radius R_0 = 10^3 and time t = 0.  
Due to the collision a double-shock structure forms and a shocked layer is created. As the simulation progresses the shocked layer widens and 
Compton drag is then applied on the shocked fluid contained between the forward and the reverse shocks. 
In order to suppress artificial transients created by abrupt changes we apply the radiation drag gradually by increasing the drag coefficient over
time as α(t) = α (1- e^-t/t_0).  As our model assumes that the drag only affects the hot shocked plasma, the in-simulation drag is applied
only in the region where the shocked fluid velocity is in the range [0.1u_e, 0.9u_e], where u_e is the 4-velocity of the unshocked ejecta. 
Preliminary  runs have shown that for the radiation free case the shocked layer becomes sufficiently developed by t=30 (see Fig. <ref>), which is why the time constant for applying the radiation drag was chosen to be t_0=10 (by t=30 the drag reaches its maximum). For sufficiently high values of α, the shocked layer decelerates and the reverse shock quickly becomes radiation-supported.




 §.§ Test case
 
As a check, we ran a test case with α=0 for the same setup described above, and compared the results with an analytic solution
obtained under the thin shell approximation, whereby the shocked layers are assumed to be uniform.   Under this approximation
the shocked ejecta and shocked ambient medium have the same Lorentz factor, γ_s=Γ_c.  
In terms of the ratios  q_e=(γ_e/Γ_r)^2 and q=(γ_s/Γ_r)^2,
the jump conditions, Eqs. (<ref>) - (<ref>), yield

    2q_e(1-√(q/q_e))+(a(q-1)+2)(q-q_e/q+1)=0,

and

    q/q_e=3/4γ_e^2(ρ_e/ρ_i)(q_e-1)(√(q/q_e))/q(q-1)+2,

where ρ_e, q and q_e are functions of time.  
Solving these equations for the given initial conditions we obtain q_e=3.765, q=1.2 at t=0.  As the ejecta expands, its density just upstream 
of the reverse shock evolves as ρ_e(t)∝ [R_r(t)]^-2. Solving the above equations at any given time t using ρ_e(t) in Eq. (<ref>), 
one obtains the Lorentz factor of the contact discontinuity, Γ_c(t)=γ_e√(q(t)/q_e(t)).  The contact 4-velocity,
U_c(t)=√(Γ_c^2(t)-1), is shown as a dashed line in Fig. <ref>, along with the 4-velocity, pressure and density 
obtained from the simulations of the radiation free case (α=0), at different times. 
As seen from Fig. <ref>, the jump at the shock agrees well with the analytic result.  




 §.§ Simulation Results

For numerical reasons we were only able to run cases with modest values of α, however these suffice to illustrate the main trends.
Below we present results  for α =3, 6, 10, 15.   The initial conditions in all cases were set as described above.  




Fig <ref> displays snapshots of the plasma 4-velocity at different times for α=10, showing the evolution of the entire 
structure.   The x-axis gives the distance from the reverse shock, x=R-R_r(t), so that the reverse shock is located 
at x=0 at all times.   The vertical red and blue
lines indicate the location of the contact discontinuity x_c and the forward shock x_f, respectively.  
 The deceleration of the shocked fluid, that leads to strengthening of the reverse shock and weakening of the forward shock is clearly
seen.   The radiation force gives rise to compression of the fluid near the contact which is
communicated  to the reverse shock and decelerates it (Fig. <ref>).   
The shocked layer gradually expands as time progresses, but initially less than in
the case α=0.  The apparent sudden expansion of the shocked layer at late times (between t=600 and t=900) commences
when the shock becomes mildly relativistic, as can be observed from Fig. <ref>.   This is mainly due to the fact that the width of the
shocked layer evolves with time as t/Γ_r^2.  We note that at such small Lorentz 
factors (Γ_r<2.5 at t>600) our choice of Compton drag terms, that assume perfect beaming of the external radiation in the 
rest frame of the shocked fluid, overestimates
the actual drag.  In reality the drag is expected to be suppressed as the shocked fluid becomes mildly relativistic , so that the Lorentz factor
of the reverse shock will eventually saturate once Γ_r becomes sufficiently low. 




  
Figs <ref> shows the profiles of the 4-velocity, pressure and density around the shock at  the same simulation time, t=800, for four values of α 
as indicated in the figure label.  The radius of  the reverse shock at this time is also indicated for convenience, and it ranges from R_r≈1776 (or (R_r-R_0)/R_0 ≃ 0.77) for α=3 to R_r≈1716 for α=15.  The result of the run with no radiation drag (α=0), presented in Fig. <ref>, is plotted as a dotted-dashed line and is exhibited for a comparison. 
The contact location is found from the density profile, noticing that, as in the radiation-free case, the density peaks towards the contact and drops to 
a minimum right after it.
In the resulting structure we see the transition from "material - supported" to  "radiation-supported" shock; as α is increased the forward shock becomes progressively weaker and ultimately negligible, while the reverse shock strengthens.  This trend justifies the neglect of the shocked ambient pressure in the 
self-similar solution outlined in section <ref>. 
 Fig <ref> confirms that, as long as the shock is sufficiently relativistic, the width of the shocked layer shrinks and the shocked fluid velocity decreases with increasing radiation drag, in accord with the self-similar solution.  
 The widening of the shocked layer for the α=10 and α=15 cases stems from the transition to the mildly relativistic regime, as explained above.   
The increase of the pressure with increasing α further indicates strong compression by the radiation force. 
The formation of the cold dense shell near the contact, seen in the lower left panel of Fig. <ref>, results from the
rapid cooling of the compressed plasma. 

The left panel of figure <ref> exhibits the evolution of the 4-velocity of the shocked ejecta.  Gradual deceleration of the shocked flow over a timescale of the order of the cooling time is observed, as expected.   The evolution of the 4-velocity of the unshocked ejecta with respect to the reverse shock is 
exhibited in the right panel of figure <ref>, indicating a substantial increase in shock efficiency with increasing α. 











§ SUMMARY AND CONCLUSIONS


We studied the effect of radiation drag on the dynamics of shocks that form in relativistic outflows. Such situations are expected in cases where the outflow propagates through a quasi-isotropic, ambient radiation field, on scales at which the inverse Compton cooling time is significantly shorter than the outflow expansion time. The observations of the continuum emission in blazars suggest that these conditions may prevail on sub-parsec to parsec scales in those objects. 

For certain profiles of the external radiation intensity and the density of the unshocked ejecta we were able to find self-similar solutions 
of the radiation hydrodynamics equations, describing a radiation-supported shock.   
We also performed 1D numerical simulations of a uniform, spherical shell interacting with an ambient medium that contains cold gas and seed radiation.  For that purpose we used the PLUTO code, that we modified to incorporate energy and momentum losses of the shocked plasma through Compton scattering. 

In both, the analytical model and the simulation results, we find significant alteration of the shock dynamics when the ratio of dynamical time and Compton cooling time exceeds a factor of a few. Quite generally, substantial radiation drag leads to a faster deceleration, strengthening of the reverse shock and weakening of the forward shock. In the self-similar model with n=1 the 4-velocity of the upstream flow with respect to the shock increases from u^'_u≃ 0.5 for a uniform ambient density profile and no radiation (α=0), to about u^'_u≃ 2.1 in the asymptotic limit of maximum efficiency (obtained for α >50). The Lorentz factor in the asymptotic limit evolves roughly as Γ∝ t^-10.  

For numerical reasons the simulations were limited to modest values of α, but show similar trends. The numerical experiments enabled us to follow the evolution of the system from the onset of fluid collision, and demonstrate the gradual strengthening of the reverse shock (and weakening of the forward shock) over time. For convenience we adopted intensity profile (of the seed radiation) that scales as r^-1, which is close to the flat profile obtained for blazars from detailed calculations.  Since for a sufficiently large drag the deceleration scale is much smaller than the shock radius, the details of the intensity profile has little effect on the evolution.    For our choice of constant drag coefficient, the deceleration continues indefinitely.   In reality the drag will be strongly suppressed once the Lorentz factor of the emitting (shocked) fluid drops to values at which our "beaming approximation" breaks down.   In this regime strong cooling still ensues, but with little momentum losses.  Our late time results for the cases α=10 and α=15 are therefore not reliable.  More detailed calculations of the drag terms are needed to follow the evolution in the mildly relativistic regime.

The effect of radiation drag on the dynamics of the shock might have important implications for the resultant emission. On the one hand, the strengthening of the reverse shock leads to enhanced efficiency of internal shocks, in particular under conditions at which the reverse shock is sub-relativistic in the absence of external radiation.  On the other hand, the deceleration of the shocked fluid, that emits the observed radiation, leads to a dramatic change in the beaming factor. This temporal change in beaming factor can significantly alter the light curves, particularly for observers viewing the source off-axis.  Detailed calculations of variable emission from dragged shocks is beyond the scope of this paper, but our analysis suggests that detailed emission models, at least in blazars, should account for such effects.



This research was supported by  a grant from the Israel Science Foundation no. 1277/13.














14
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Joshi, Marscher, and Bottcher(2014)]JMB14
  author author M. Joshi, author A. P. Marscher,
   and author M. Bottcher,  title title Seed photon fields of
  blazars in the internal shock scenario, @noop  journal journal Astrophys. J. volume
  785, pages 132 (year 2014)NoStop
[Sikora et al.(2009)Sikora,
  Stawarz, Moderski, Nalewajko, and Madejski]SSMNM09
  author author M. Sikora, author L. Stawarz,
  author R. Moderski, author K. Nalewajko,  and author G. M. Madejski, title title Constraining emission models of luminous
  blazar sources, @noop  journal journal
  Astrophys. J. volume 704, pages
  38–50 (year 2009)NoStop
[Dermer and Schlickeiser(1993)]DS93
  author author C. D. Dermer and author R. Schlickeiser, title title Model for
  the high-energy emission from blazars, @noop  journal journal Astrophys. J. volume
  416, pages 458 (year 1993)NoStop
[Sikora, Begelman, and  Rees(1994)]SBR94
  author author M. Sikora, author M. C. Begelman,  and author M. J. Rees, title title Comptonization of
  diffuse ambient radiation by a relativistic jet: The source of gamma rays
  from blazars? @noop  journal journal
  Astrophys. J. volume 421, pages
  153–162 (year 1994)NoStop
[Blandford and Levinson(1995)]BL95
  author author R. D. Blandford and author A. Levinson, title title Pair cascades
  in extragalactic jets. 1: Gamma rays, @noop  journal
  journal Astrophys. J. volume 441,  pages 79–95 (year 1995)NoStop
[Ghisellini and Madau(1996)]GM96
  author author G. Ghisellini and author P. Madau, title title On the origin of
  the gamma-ray emission in blazars, @noop  journal
  journal Mon. Not. Roy. Ast. Soc. volume
  280, pages 67–76 (year 1996)NoStop
[Joshi and Bottcher(2011)]JB11
  author author M. Joshi and author M. Bottcher, title title Time-dependent
  radiation transfer in the internal shock model scenario for blazar jets,  @noop  journal journal Astrophys. J.  volume 727, pages 21 (year
  2011)NoStop
[Nakamura and Shigeyama(2006)]NS06
  author author K. Nakamura and author T. Shigeyama, title title Self-similar
  solutions for the interaction of relativistic ejecta with an ambient
  medium, @noop  journal journal
  Astrophys. J. volume 645, pages
  431–4325 (year 2006)NoStop
[Levinson(2010)]Le10
  author author A. Levinson, title title Relativistic
  rayleigh-taylor instability of a decelerating shell and its implications for
  gamma-ray bursts, @noop  journal journal GApFD volume 104, pages
  85–111 (year 2010)NoStop
[Golan and Levinson(2015)]GL15
  author author O. Golan and author A. Levinson, title title Blazar flares
  from compton dragged shells, @noop  journal
  journal Astrophys. J. volume 809,  pages 23 (year 2015)NoStop
[Levinson(2009)]Le09
  author author A. Levinson, title title Convective
  instability of a relativistic ejecta decelerated by a surrounding medium: An
  origin of magnetic fields in gamma-ray bursts? @noop  journal journal Astrophys. J. volume
  705, pages L213–L216 (year 2009)NoStop
[Blandford and McKee(1976)]BMK76
  author author R. D. Blandford and author C. F. McKee, title title Fluid dynamics of
  relativistic blast waves, @noop  journal journal Phys. Fluids volume 19, pages 1130–1138 (year 1976)NoStop
[Note1()]Note1
  note More generally, α can be taken to be any function of
  the similarity coordinate χ. In the presence of strong drag the width
  of the shocked layer is much smaller the its radius, hence we anticipate
  little variations in α across it.Stop
[Mignone et al.(2007)Mignone, Bodo, Massaglia, and  et al.]PLUTO07
  author author A. Mignone, author G. Bodo,
  author S. Massaglia,  and  author et al., title title Pluto: A numerical code for computational astrophysics,  @noop  journal journal Astrophys. J.
  Supp. volume 170, pages 228 (year 2007)NoStop



