Dynamics of Relativistic Shock Waves Subject to a Strong Radiation Drag: Similarity Solutions and Numerical Simulations
DDynamics of Relativistic Shock Waves Subject to a StrongRadiation Drag: Similarity Solutions and NumericalSimulations
Ilia Leitus and Amir Levinson a) Raymond and Beverly Sackler School of Physics & Astronomy, Tel Aviv University,Tel Aviv 69978, Israel (Dated: 8 November 2018)
We examine the effect of Compton drag on the dynamics of a relativisticshock wave. Similarity solutions describing a radiation-supported shock areobtained for certain profiles of the external radiation intensity and the densityof the unshocked ejecta, and are compared with 1D numerical simulations ofa blast wave expanding into an ambient medium containing isotropic seedradiation. Both the analytic model and the simulations indicate that underrealistic conditions the radiation drag should strongly affect the dynamics andstructure of the shocked layer. In particular, our calculations show significantstrengthening of the reverse shock and weakening of the forward shock overtime of the order of the inverse Compton cooling time. We conclude that theeffect of radiation drag on the evolution of the emitting plasma should affectthe resultant light curves and, conceivably, spectra of the observed emissionduring strong blazar flares.
I. 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 byan accreting black hole or a magnetar. The interaction of these outflows with their environ-ments results in formation of expanding cocoons and shocks, that decelerate the flow andlead to lower energy emissions over a large range of scales, e.g., afterglow emission in GRBsand radio lobes in AGNs and micro-quasars. In certain circumstances the relativistic jetsinteract also with ambient radiation fields emitted by the surrounding gas. This interactioncan considerably alter the dynamics of dissipative fluid shells at relatively small scales, andis likely to be the origin of the variable gamma-energy emission observed in certain classesof high-energy sources.Direct evidence for such interactions is found in blazars, where the properties of theambient radiation field can be measured. This radiation is contributed primarily by theaccretion disk around the black hole, by gaseous clouds in the broad line region (BLR) and,at larger scales, by a dusty molecular torus . The conventional wisdom has been thatthe high-energy emission seen in powerful blazars results from inverse Compton scatteringof ambient seed photons by non-thermal electrons accelerated in dissipative regions in thejet . 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 lightcurves of flared emission . However, the effect of the radiation drag on the structure anddynamics of the shock has been neglected.In this paper we construct a model for the dynamics of a relativistic shock wave in thepresence of intense radiation field. In section II we derive the governing equations of aradiation-supported shock. In section III we present self-similar solutions, obtained underseveral simplifying assumptions regarding the properties of the ejecta and the intensity pro-file of the external radiation field, and study the dependency of the evolved structure on the a) Electronic mail: [email protected] a r X i v : . [ a s t r o - ph . H E ] J u l magnitude of the radiation drag. In section IV we present results of numerical simulationsof a uniform, spherical ejecta propagating through an ambient medium containing isotropicseed radiation. The simulations enable analysis of the temporal evolution from the initialstage under more realistic conditions. While the analysis is motivated by the application toblazars, mainly because there we have sufficient information to asses the conditions in thesource, it might also be relevant to other high-energy astrophysical systems. II. 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 ambientmedium and collision of fast and slow shells (internal shocks). In the absence of externalforces, the dynamics of the double-shock system is dictated solely by the properties of theunshocked media. In the presence of external radiation field, the shocked layers are subjectto radiation drag that can alter the dynamics significantly. When the drag is strong enoughit 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 shellin case of internal shocks) is supported by the radiation pressure imposed on it. Thus, asa simple approximation one can ignore the entire region beyond the contact discontinuity,and obtain solutions describing a radiation-supported shock. This is the treatment adoptedin section III below. The validity of this approximation is checked in section IV usingnumerical simulations.Below we adopt the following notation: The velocities of the reverse shock and the contactsurface 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 referto the unshocked and shocked ejecta, respectively. A. Flow equations
The dynamics of a relativistic fluid is governed by the continuity equation, ∂ µ ( ρu µ ) = 0 , (1)and energy-momentum equations, ∂∂x α T αµ = S µ , (2)here expressed in terms of the energy-momentum tensor T µν = w u µ u ν + g µν p, (3)where ρ , p and w are, respectively, the proper density, pressure and specific enthalpy, u µ = γ (1 , v , v , v ) is the fluid 4-velocity, and g µν is the metric tensor, here taken to be theMinikowsky metric g µν = diag( − , , , w = ρ + ˆ γ ˆ γ − p ≡ ρ + ap, (4)with a = 4 for a relativistically hot plasma and a = 5 / r, θ, φ ). Equations (1) and (2) for the shocked fluid then reduce to: ∂ t ( ργ ) + 1 r ∂ r ( r ργv ) = 0 , (5) ∂ t ( wγ − p ) + 1 r ∂ r ( r wγ v ) = S , (6) w d ln γdt + dpdt − γ − ∂ t p = v ( S r − vS ) , (7)here d/dt = ∂ r + v∂ r denotes the convective derivative.A common approach is to seek solutions of Eqs. (5)-(7) inside the shocked layer, given theproperties of the unshocked flows upstream of the forward and reverse shocks as input. Thejump conditions at the forward and reverse shocks then serve as boundary conditions for theshocked flow equations. In addition, the fluid velocity and total momentum flux must becontinuous across the contact discontinuity surface. In the approximate treatment adoptedhere, whereby the reverse shock is fully supported by radiation and the layer enclosedbetween the contact and the forward shock can be ignored, the location of the contact isfixed by the requirement that the total energy is conserved, given formally by Eq. (19)below. B. Jump conditions at the reverse shock
For the range of conditions envisaged here the shocked plasma is optically thin. Theshock, in this case, is collisionless, i.e., mediated by collective plasma processes (as opposedto radiation mediated shocks), and its width, roughly the skin depth, is much smaller thanthe Thomson length. Thus, the shock transition layer can be treated as a discontinuity inthe flow parameters. The local jump conditions at the reverse shock can be expressed interms of the local shock velocity V r = dR r /dt , where R r ( t ) is the shock radius at time t ,as ρ u γ u ( v u − V r ) = ρ s γ s ( v s − V r ) , (8) w u γ u v u ( v u − V r ) + p u = w s γ s v s ( v s − V r ) + p s , (9) w u γ u ( v u − V r ) + p u V r = w s γ s ( v s − V r ) + p s V r . (10)We shall henceforth assume that the unshocked medium is cold and set p u = 0. The jumpconditions at the forward shock can, in principle, be derived in a similar manner, but areredundant under our radiation-supported shock approximation. C. Compton drag terms
Behind the shock electrons are heated and accelerated to relativistic energies. We denotethe energy distribution of electrons there by F e ( γ e ), where γ e = (cid:15)/m e c is the dimensionlesselectron energy, as measured in the fluid frame. The total electron proper density is thengiven by n e = (cid:82) F e ( γ e ) dγ e , and the average energy and second energy moment by < γ e > = 1 n e (cid:90) γ e F e ( γ e ) dγ e , < γ e > = 1 n e (cid:90) γ e F e ( γ e ) dγ e . (11)If only the thermal electrons are taken into account one expects m e < γ e > (cid:39) m p γ (cid:48) u / γ (cid:48) u = γ u Γ r (1 − v u V r ) is the Lorentz factor of theupstream flow, as measured in the shock frame, and it is assumed that γ (cid:48) u >
1. If thecontribution of non-thermal electrons is substantial, then < γ e > may be larger. Since theshocked electron plasma is relativistic we have m e c n e < γ e > = w < γ e > .The derivation of the source terms is given in Ref 10. In terms of the energy density ofthe radiation intercepted by the flow, u rad ( r ), and the quantities defined above they can beexpressed as S = − γ s < γ e > u rad σ T n e , (12)and S r = v s S + S / γ s . (13)From Eqs. (7) and (13) it is readily seen that the deceleration time of the shocked fluidis given by t dec = − wγ s /v s S = 2 γ s t (cid:48) IC , where t (cid:48) IC (cid:39) m e c/ (4 σ T u (cid:48) rad < γ e > ) is theinverse Compton cooling time of an electron having an energy m e c < γ e > , as measured inthe fluid rest frame, and u (cid:48) rad = 4 γ s u rad /
3. Substantial deceleration occurs if t dec is muchshorter than the outflow time, t f = r/c . For convenience, we define the drag coefficient as α = 3 t f /t dec = 83 m e c γ s < γ e > σ T u rad r. (14)As noted above, for thermal electrons γ s < γ e > (cid:39) ( m p / m e ) γ (cid:48) u γ s = ( m p / m e ) γ u . If onlythe contribution of thermal electrons is accounted for, then α = 8 m p σ T m e c γ u u rad r. (15)The energy loss term can now be expressed as: S = − αγ s wr − . (16)To get an estimate of the value of α anticipated in blazars, suppose that a fraction η ofthe nuclear luminosity, L = 10 L erg/s, is scattered across the jet. The correspondingenergy density of the radiation intercepted by the jet is roughly u rad ( r ) = ηL/ πr =3 × − ηL /r pc ergs cm − , where r pc is the radius in parsecs. From Eq. (15) one obtains α = 0 . γ u ηL r − pc . A more realistic estimate can be obtained using recent calculations ofthe ambient radiation intensity contributed by the various radiation sources surroundinga blazar jet . It is found that for a prototypical blazar the intensity profile of radiationintercepted by the jet is flat up to a distance of about one parsec, with u rad (cid:39) − ergcm − , followed by a gradual decline . For an internal shock produced by colliding shells,Equation (15) yields for the drag coefficient α ≈ γ u r pc . Typically γ u >
10, implyingsignificant drag on sub-parsec scales in blazars.
D. 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 theenergy radiated away through the contact. The energy accumulated behind the forwardshock is neglected in our approximate treatment. The time change of the total energy persteradian of the shocked ejecta is ∂ t E = ∂ t (cid:90) R c ( t ) R r ( t ) T r dr = (cid:90) R c ( t ) R r ( t ) ∂ t ( T r ) dr + R c T s ( R c ) V c − R r T s ( R r ) V s , (17)where R r ( t ) and R c ( t ) are, respectively, the radii of the reverse shock and the contactsurface, and V r = dR r /dt , V c = dR c /dt are the corresponding velocities. Integrating Eq.(6) from the shock, R r ( t ), to the contact, R c ( t ), one obtains: (cid:90) R c ( t ) R r ( t ) ∂ t ( T r ) dr + R c w s ( R c ) γ s ( R c ) v s ( R c ) − R r w s ( R r ) γ s ( R r ) v s ( R r ) = (cid:90) R c ( t ) R r ( t ) S r dr. (18)Combining Eqs. (17) and (18), using the jump condition (10) with p u = 0 (cold ejecta) andthe relation T = wγ − p , and recalling that the fluid velocity is continuous across thecontact, viz., v s ( R c ) = V c , yields ∂ t E = R r ρ u γ u ( v u − V r ) − R c p s ( R c ) V c + (cid:90) R c ( t ) R r ( t ) S r dr. (19)The first term on the right hand side accounts for the power incident through the reverseshock, the second term for pdV work at the contact and the last term for radiative losses. III. SELF-SIMILAR SOLUTIONS
Similarity solutions describing the interaction of a relativistic shell with an ambientmedium, in the absence of radiative losses, were derived in Ref 8, and their stability wassubsequently analyzed . Such solutions can be obtained for a freely expanding ejectacharacterized by a velocity profile v u = r/t at time t after the explosion, and a properdensity of the form ρ u = a u t γ nu , (20)where a u is a normalization constant, and γ u = 1 / (cid:112) − v u is the corresponding Lorentzfactor (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 andthe contact discontinuity have a similar time evolution, viz., Γ f = At − m , Γ r = Bt − m , Γ c = Ct − m , where A, B, C and m are constants determined upon matching the solutions in region1 (shocked ejecta) and region 2 (shocked ambient medium) at the contact discontinuity .For an ambient density profile of the form ρ i ∝ r − k , the index m is given by m = 6 − kn + 2 . (21)In cases where the shocked plasma is subject to strong radiative losses it is still possible toobtain self-similar solutions provided the energy source term scales as S ∝ γ s wr − (see Eq.(6) and (7)), that is, α in Eq. (16) is constant. We shall henceforth make this assumptioneven 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 theshock these details are unimportant. This is confirmed by numerical simulations, presentedin the next section.As explained above, in the presence of a strong radiation drag the dynamics of thesystem is dictated by the interaction of the external radiation field with the shocked plasmaenclosed between the contact and the reverse shock. To a good approximation one can thenignore the contribution of the forward shock to the overall dynamics. The evolution of theradiation-supported shock still satisfies Γ r = Bt − m , Γ c = Ct − m , but the index m is nowdetermined from global energy conservation, as will be shown below. Now, to order O (Γ − r )the trajectory of the reverse shock is given by R r ( t ) = (cid:90) t (cid:18) − r (cid:19) d t (cid:48) = t − t m + 1)Γ r , (22)from which we obtain for the velocity of the ejecta crossing the shock: v u ( R r ) = R r /t =1 − / [2( m + 1)Γ r ]. The corresponding Lorentz factor is thus given, to the same order, by γ u = ( m + 1)Γ r , (23)and the density profile by ρ u = a u ( m + 1) n/ B /m Γ (6 /m − n ) r . (24)We adopt the similarity parameter introduced in Ref 12, χ = [1 + 2( m + 1)Γ r ](1 − r/t ) . (25)With this choice the reverse shock is located at χ = 1, and the contact at χ c = Γ r / Γ c = B/C <
1. Following Ref. 8 we define the self-similar variables, G ( χ ), H ( χ ) and F ( χ ), suchthat at the reverse shock they satisfy the boundary conditions G (1) = H (1) = F (1) = 1.Upon solving the jump conditions (8)-(10), the shocked fluid quantities can be expressed interms of the self-similar variables as γ s = q Γ r G ( χ ) , (26) ρ s γ s = mqρ u γ u ( m + 1)( q − H ( χ ) , (27) p s = mρ u a ( q −
1) + 2 (1 − (cid:112) q/ ( m + 1)) F ( χ ) , (28) w s = ρ s + ap s ≡ K ( χ ) p s , (29)where the parameter √ q is the only positive root of the polynomial equationˆ γx + (2 − ˆ γ ) √ m + 1 x − (2 − ˆ γ ) x − ˆ γ √ m + 1 x = 0 . (30)Upon substituting these relations into Eqs. (5)-(7), and using Eq. (16) with a constantdrag coefficient, α = const , we obtain the following set of equations for the self-similarvariables: 2(1 + qGχ ) ∂ χ ln F − (1 − qGχ ) K∂ χ ln G = mn − − ( m − α ) K ( m + 1) qG, (31)2(1 − qGχ ) ∂ χ ln F − ˆ γ (1 + qGχ ) ∂ χ ln G = 6 − mn + ( m − γ − (ˆ γ − αKm + 1 qG, (32)2(1 − qGχ ) ∂ χ ln H − ∂ χ ln G = − ( mn − m − m + 1) qG, (33)subject to the boundary conditions G (1) = F (1) = H (1) = 1. These equations reduce tothose derived in Ref. 8 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 . Thisrelation defines the limit of integration.The self-similar equations involve two eigenvalues; the index m and the location of thecontact χ c . Thus, two conditions are needed to find them. The first one is the relation qG c χ c = 1 (34)mentioned above. The second one is global energy conservation, Eq. (19). To order O ( γ − s )we have T s = w s γ s − p s (cid:39) w s γ s , r = t and dr = − tdχ/ [2( m + 1)Γ r ], from which we obtain E ( t ) = (cid:90) R c ( t ) R r ( t ) T r dr = t nm/ qma u (1 − (cid:112) q/ ( m + 1))2( m + 1) n/ B n/ [ a ( q −
1) + 2] × (cid:90) χ c K ( χ ) F ( χ ) G ( χ ) dχ, (35)and ∂ t E = nmE/ t . To the same order one has (cid:82) R c R r S r dr = − αE/t . Substituting theseresults into Eq. (19) yields the constraint (cid:90) χ c K ( χ ) F ( χ ) G ( χ ) dχ = ( m + 1) q ( α + nm/ × (cid:34) a ( q −
1) + 21 − (cid:112) q/ ( m + 1) − F c (cid:35) , (36)where F c = F ( χ c ). A. Results
For a given choice of the drag coefficient α we guess the value of m and integrate Eqs.(31)-(33) from χ = 1 to the point χ c at which Eq. (34) is satisfied. We then check if Eq.(36) is satisfied. If not, we change the value of m and repeat the process until a solutionsatisfying all constraints is found. Sample profiles are shown in Fig. 1 for n = 1 and differentvalues of α . The radiation free case ( α = 0) is shown for a comparison. It was computedusing the full solution of the two-shock model described in Ref. 8. As seen, the width ofthe shocked layer, ∆ χ = 1 − χ c , decreases with increasing drag coefficient α , as naivelyexpected. Moreover, larger radiative losses lead to increased non-uniformity of the Lorentzfactor and pressure in the shocked layer. The divergence of the density at the contact is abasic feature of the similarity solutions and occurs even in the absence of radiative losses( α = 0), as seen in the upper panel of Fig. 1.Figure 2 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. (21) with k = 0 and n = 1). The 4-velocity of the unshockedfluid, as measured in the frame of the reverse shock, is given by u (cid:48) u = γ u Γ r ( v u − V r ) = m √ m + 1 , (37)and it is seen that the shock becomes substantially stronger as α increases. Note that thepower dissipated behind the shock, ρ u γ (cid:48) u u (cid:48) u , increases roughly linearly with the index m .At sufficiently large drag the entire power incident through the shock is radiated away, andthe solution becomes independent of α , as seen in Fig. 2. For our choice of parameters,specifically n = 1, this occurs at α > ∼
50, for which m (cid:39)
19 and u (cid:48) u (cid:39) .
1. We emphasizethat this limit can be approached provided the Lorentz factor γ u of the unshocked shellis sufficiently large. To be more concrete, Eq. (23) implies that γ u (cid:39) . r as α → ∞ ,while our analysis is valid only for Γ r >>
1. The Compton drag terms given in Eqs. (12)and (13) assume that the intensity of ambient radiation is highly beamed in the frame ofthe shocked fluid. Once γ s decelerates to modest values, γ s > ∼
1, the drag force is stronglyreduced, ultimately becoming ineffective.
IV. NUMERICAL SIMULATIONS
Numerical simulations were performed using the PLUTO code , that was modified toincorporate energy and momentum losses of the shocked plasma through Compton scat-tering. To clarify the presentation we adopted an external intensity profile of the form u s ∝ r − for which α in Eq. (16) is constant. Over the deceleration scale this profile is areasonable approximation to the flat profile expected in blazars . We start with the basicspherical blast wave problem, releasing an ejecta with a uniform Lorentz factor γ e = 10into a stationary ambient medium (henceforth, quantities in the unshocked ejecta are des-ignated by subscript ”e”). In our setup, both the ejecta and the ambient medium are takento be uniform initially with a density ratio ρ e /ρ i = 100. The calculations are performed in FIG. 1. Profiles of the normalized density (top panel), pressure (middle panel) and Lorentz factor(bottom panle), plotted as functions of the self-similar parameter χ , for different values of the dragcoefficient α , as indicated in the upper panel.FIG. 2. Dependence of the index m on α . ”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 = 10 and time t = 0. Due to the collision a double-shock structure forms and a shocked layeris created. As the simulation progresses the shocked layer widens and Compton drag isthen 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 radiationdrag gradually by increasing the drag coefficient over time as α ( t ) = α (1 − e − t/t ). As ourmodel assumes that the drag only affects the hot shocked plasma, the in-simulation dragis applied only in the region where the shocked fluid velocity is in the range [0 . u e , . u e ],where u e is the 4-velocity of the unshocked ejecta. Preliminary runs have shown that forthe radiation free case the shocked layer becomes sufficiently developed by t = 30 (see Fig.3), which is why the time constant for applying the radiation drag was chosen to be t = 10(by t = 30 the drag reaches its maximum). For sufficiently high values of α , the shockedlayer decelerates and the reverse shock quickly becomes radiation-supported. A. Test case
As a check, we ran a test case with α = 0 for the same setup described above, andcompared the results with an analytic solution obtained under the thin shell approximation,whereby the shocked layers are assumed to be uniform. Under this approximation theshocked ejecta and shocked ambient medium have the same Lorentz factor, γ s = Γ c . Interms of the ratios q e = ( γ e / Γ r ) and q = ( γ s / Γ r ) , the jump conditions, Eqs. (8) - (10),yield 2 q e (cid:18) − (cid:114) qq e (cid:19) + ( a ( q −
1) + 2) (cid:18) q − q e q + 1 (cid:19) = 0 , (38)and qq e = 34 γ e (cid:18) ρ e ρ i (cid:19) ( q e − (cid:16)(cid:113) qq e (cid:17) q ( q −
1) + 2 , (39)where ρ e , q and q e are functions of time. Solving these equations for the given initialconditions we obtain q e = 3 . q = 1 . t = 0. As the ejecta expands, its density justupstream of the reverse shock evolves as ρ e ( t ) ∝ [ R r ( t )] − . Solving the above equationsat any given time t using ρ e ( t ) in Eq. (39), one obtains the Lorentz factor of the contactdiscontinuity, Γ c ( t ) = γ e (cid:113) q ( t ) q e ( t ) . The contact 4-velocity, U c ( t ) = (cid:112) Γ c ( t ) −
1, is shown asa dashed line in Fig. 3, along with the 4-velocity, pressure and density obtained from thesimulations of the radiation free case ( α = 0), at different times. As seen from Fig. 3, thejump at the shock agrees well with the analytic result. B. Simulation Results
For numerical reasons we were only able to run cases with modest values of α , howeverthese suffice to illustrate the main trends. Below we present results for α = 3 , , ,
15. Theinitial conditions in all cases were set as described above.Fig 4 displays snapshots of the plasma 4-velocity at different times for α = 10, showingthe 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 redand 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 thereverse shock and weakening of the forward shock is clearly seen. The radiation force givesrise to compression of the fluid near the contact which is communicated to the reverse shockand decelerates it (Fig. 5). The shocked layer gradually expands as time progresses, butinitially less than in the case α = 0. The apparent sudden expansion of the shocked layerat late times (between t = 600 and t = 900) commences when the shock becomes mildlyrelativistic, as can be observed from Fig. 5. This is mainly due to the fact that the widthof the shocked layer evolves with time as t/ Γ r . We note that at such small Lorentz factors(Γ r < . t > FIG. 3. Profiles of the 4-velocity (blue solid line), pressure (black solid line), and density (red solidline) obtained from the test run of a radiation free system ( α = 0), plotted as a functions of thedistance from the reverse shock, x = R − R r ( t ). The value of the contact 4-velocity U c , computedfrom the analytic model, is marked by the dashed line. As seen, the agreement is excellent.FIG. 4. Profiles of the plasma 4-velocity (solid line) and pressure (dashed line) at different times,for the case α = 10. For clarity of presentation, the pressure is scaled by a factor of 5. Timeprogresses from top to bottom panels, as indicated. drag. In reality the drag is expected to be suppressed as the shocked fluid becomes mildlyrelativistic , so that the Lorentz factor of the reverse shock will eventually saturate once Γ r becomes sufficiently low.Figs 6 shows the profiles of the 4-velocity, pressure and density around the shock at thesame simulation time, t = 800, for four values of α as indicated in the figure label. Theradius of the reverse shock at this time is also indicated for convenience, and it ranges from R r ≈ R r − R ) /R (cid:39) .
77) for α = 3 to R r ≈ α = 15. The result ofthe run with no radiation drag ( α = 0), presented in Fig. 3, is plotted as a dotted-dashedline and is exhibited for a comparison. The contact location is found from the densityprofile, noticing that, as in the radiation-free case, the density peaks towards the contact1 FIG. 5. Time evolution of the reverse shock Lorentz factor. Up to t ≈
30 the shock structure isnot sufficiently developed and the determination of Γ r is uncertain. 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 shockbecomes progressively weaker and ultimately negligible, while the reverse shock strengthens.This trend justifies the neglect of the shocked ambient pressure in the self-similar solutionoutlined in section III. Fig 6 confirms that, as long as the shock is sufficiently relativistic, thewidth of the shocked layer shrinks and the shocked fluid velocity decreases with increasingradiation drag, in accord with the self-similar solution. The widening of the shocked layerfor 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 strongcompression by the radiation force. The formation of the cold dense shell near the contact,seen in the lower left panel of Fig. 6, results from the rapid cooling of the compressedplasma.The left panel of figure 7 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 isobserved, as expected. The evolution of the 4-velocity of the unshocked ejecta with respectto the reverse shock is exhibited in the right panel of figure 7, indicating a substantialincrease in shock efficiency with increasing α . V. SUMMARY AND CONCLUSIONS
We studied the effect of radiation drag on the dynamics of shocks that form in relativisticoutflows. Such situations are expected in cases where the outflow propagates through aquasi-isotropic, ambient radiation field, on scales at which the inverse Compton coolingtime is significantly shorter than the outflow expansion time. The observations of thecontinuum emission in blazars suggest that these conditions may prevail on sub-parsec toparsec scales in those objects.For certain profiles of the external radiation intensity and the density of the unshockedejecta we were able to find self-similar solutions of the radiation hydrodynamics equations,describing a radiation-supported shock. We also performed 1D numerical simulations ofa uniform, spherical shell interacting with an ambient medium that contains cold gas andseed radiation. For that purpose we used the PLUTO code, that we modified to incorporateenergy and momentum losses of the shocked plasma through Compton scattering.In both, the analytical model and the simulation results, we find significant alteration ofthe shock dynamics when the ratio of dynamical time and Compton cooling time exceedsa 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-similarmodel with n = 1 the 4-velocity of the upstream flow with respect to the shock increasesfrom u (cid:48) u (cid:39) . α = 0), to about u (cid:48) u (cid:39) . α > ∝ t − .For numerical reasons the simulations were limited to modest values of α , but show similartrends. The numerical experiments enabled us to follow the evolution of the system from the2onset of fluid collision, and demonstrate the gradual strengthening of the reverse shock (andweakening of the forward shock) over time. For convenience we adopted intensity profile(of the seed radiation) that scales as r − , which is close to the flat profile obtained forblazars from detailed calculations. Since for a sufficiently large drag the deceleration scaleis much smaller than the shock radius, the details of the intensity profile has little effecton the evolution. For our choice of constant drag coefficient, the deceleration continues FIG. 6. Profiles of the 4-velocity (upper panel), density (lower left panel) and pressure (lower rightpanel) at simulation time t = 800, plotted as a functions of the distance from the reverse shock, fordifferent values of the drag coefficient α . The location of the reverse shock is at x = 0. The locationsof the contact and the forward shock are marked by red and blue vertical lines respectively. Thethin dotted-dashed lines correspond to the case α = 0, and are shown for comparison. As seen, thestrengthening of the reverse shock and weakening of the forward shock become more substantialas α increases. Strong compression by the radiation force is also seen in the pressure and densityplots. Note the logarithmic scale of the y-axis in the density plots. FIG. 7. Time evolution of the 4-velocity of the shocked ejecta (left panel), and of the unshockedejecta, as measured in the shock frame (right panel). indefinitely. In reality the drag will be strongly suppressed once the Lorentz factor of theemitting (shocked) fluid drops to values at which our ”beaming approximation” breaksdown. In this regime strong cooling still ensues, but with little momentum losses. Ourlate time results for the cases α = 10 and α = 15 are therefore not reliable. More detailedcalculations of the drag terms are needed to follow the evolution in the mildly relativisticregime.The effect of radiation drag on the dynamics of the shock might have important implica-tions for the resultant emission. On the one hand, the strengthening of the reverse shockleads to enhanced efficiency of internal shocks, in particular under conditions at which thereverse 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 dramaticchange in the beaming factor. This temporal change in beaming factor can significantlyalter the light curves, particularly for observers viewing the source off-axis. Detailed calcu-lations of variable emission from dragged shocks is beyond the scope of this paper, but ouranalysis suggests that detailed emission models, at least in blazars, should account for sucheffects. ACKNOWLEDGMENTS
This research was supported by a grant from the Israel Science Foundation no. 1277/13. M. Joshi, A. P. Marscher, and M. Bottcher, “Seed photon fields of blazars in the internal shock scenario,”Astrophys. J. , 132 (2014). M. Sikora, L. Stawarz, R. Moderski, K. Nalewajko, and G. M. Madejski, “Constraining emission modelsof luminous blazar sources,” Astrophys. J. , 38–50 (2009). C. D. Dermer and R. Schlickeiser, “Model for the high-energy emission from blazars,” Astrophys. J. ,458 (1993). M. Sikora, M. C. Begelman, and M. J. Rees, “Comptonization of diffuse ambient radiation by a relativisticjet: The source of gamma rays from blazars?” Astrophys. J. , 153–162 (1994). R. D. Blandford and A. Levinson, “Pair cascades in extragalactic jets. 1: Gamma rays,” Astrophys. J. , 79–95 (1995). G. Ghisellini and P. Madau, “On the origin of the gamma-ray emission in blazars,” Mon. Not. Roy. Ast.Soc. , 67–76 (1996). M. Joshi and M. Bottcher, “Time-dependent radiation transfer in the internal shock model scenario forblazar jets,” Astrophys. J. , 21 (2011). K. Nakamura and T. Shigeyama, “Self-similar solutions for the interaction of relativistic ejecta with anambient medium,” Astrophys. J. , 431–4325 (2006). A. Levinson, “Relativistic rayleigh-taylor instability of a decelerating shell and its implications for gamma-ray bursts,” GApFD , 85–111 (2010). O. Golan and A. Levinson, “Blazar flares from compton dragged shells,” Astrophys. J. , 23 (2015). A. Levinson, “Convective instability of a relativistic ejecta decelerated by a surrounding medium: Anorigin of magnetic fields in gamma-ray bursts?” Astrophys. J. , L213–L216 (2009). R. D. Blandford and C. F. McKee, “Fluid dynamics of relativistic blast waves,” Phys. Fluids , 1130–1138 (1976). More generally, α can be taken to be any function of the similarity coordinate χ . In the presence of strongdrag the width of the shocked layer is much smaller the its radius, hence we anticipate little variations in α across it. A. Mignone, G. Bodo, S. Massaglia, and et al., “Pluto: A numerical code for computational astrophysics,”Astrophys. J. Supp.170