Exact solution of one dimensional relativistic jet with relativistic equation of state
Raj Kishor Joshi, Indranil Chattopadhyay, Dongsu Ryu, Lallan Yadav
aa r X i v : . [ a s t r o - ph . H E ] F e b MNRAS , 1–19 (2020) Preprint 8 February 2021 Compiled using MNRAS L A TEX style file v3.0
Exact solution of one dimensional relativistic jet with relativisticequation of state
Raj Kishor Joshi , , Indranil Chattopadhyay ⋆ , Dongsu Ryu , Lallan Yadav Aryabhatta Research Institute of Observational Sciences (ARIES), Manora Peak, Nainital-263002, India. Department of Physics, School of Natural Science, UNIST, Ulsan 44919, Korea. Deen Dayal Upadhyaya Gorakhpur University, Gorakhpur, Uttar Pradesh-273009, India
ABSTRACT
We study the evolution of one-dimensional relativistic jets, using the exact solution of the Riemann problem forrelativistic flows. For this purpose, we solve equations for the ideal special relativistic fluid composed of dissimilarparticles in flat space-time and the thermodynamics of fluid is governed by a relativistic equation of state. We obtainthe exact solution of jets impinging on denser ambient media. The time variation of the cross-section of the jet-head is modeled and incorporated. We present the initial condition that gives rise to a reverse shock. If the jet-headcross-section increases in time, the jet propagation speed slows down significantly and the reverse-shock may recedeopposite to the propagation direction of the jet. We show that the composition of jet and ambient medium can affectthe jet solution significantly. For instance, the propagation speed depends on the composition and is maximum fora pair-dominated jet, rather than a pure electron-positron or electron-proton jet. The propagation direction of thereverse-shock may also strongly depend on the composition of the jet.
Key words: hydrodynamics; relativistic processes; shock waves; stars: jets; galaxies: jets
CONTENTS ⋆ Email: [email protected] (IC)
Astrophysical jets were first discovered by Curtis (Curtis1918) while studying the M87 galaxy in optical wavelength.With the advent of radio telescopes, the subject of astro-physical jets started to flourish. Now astrophysical jets arestudied at least in three wavebands, radio, optical, and X-rays. All these detailed observations revealed that astrophys-ical jets are quite common and are associated with a varietyof objects like AGNs (active galactic nuclei), microquasars,GRBs (gamma-ray bursts), and even YSOs (young stellarobjects), albeit in a vastly varied length scales and ener-getics. Out of these, the jets associated with AGNs, micro-quasars, and GRBs are relativistic in terms of bulk speedas well as, temperature. After analysing the observations, ageneral picture emerged that the jets originate from a re-gion very close to the compact objects like black holes re-siding at the heart of AGNs and microquasars (for exam-ple in M87 Junor et. al. 1999; Doelemann et. al. 2012). Fromthe detailed observations of microquasars, the jet states werefound to correlate with the accretion states, suggesting thatthe jets originated from the accretion discs (Galo et. al. 2003;Fender et. al. 2010; Rushton et. al. 2010). However, even af-ter the vast improvement in observational facilities, there isno consensus about the composition, formation, and collima-tion of the jets.Investigations of astrophysical jets, on one hand, involvethe formation problem and on the other, the propagationproblem. In this paper, we are interested to study the propa- © Joshi, Chattopadhyay, Ryu & Yadav gation properties of jets. One of the earliest models for the jetflow was proposed to be a version of the de Laval nozzle byBlandford & Rees (1974). This paper outlined the basic jetstructure of a supersonic jet beam far away from the centralobject. The jet after interaction with the ambient medium,develops a reverse shock in the jet beam, a contact disconti-nuity, and a forward shock. The separation surface betweenthe jet and the ambient medium is called the contact disconti-nuity (CD). The forward shock (FS) is the shock transition inthe ambient medium as the jet head rams through the ambi-ent medium. And the reverse shock (RS) is in the jet beam.RS and FS are on either side of a CD. Marscher & Grear(1985) associated the flares in the astrophysical jets with trav-eling shocks along an adiabatic conical jet.Relativistic jets are studied by solving the equations of rel-ativistic hydrodynamics (RHD). Early attempts involved us-ing explicit finite difference code with monotonic transport(Wilson 1972). However these kind of codes depended on ar-tificial viscosity techniques to capture shocks. Studies of theinteraction of the jet with the ambient medium and the evo-lution of jet morphology got a tremendous boost with theadvent of high resolution shock capturing (or, HRSC) nu-merical simulation methods (van Putten 1993; Marti et. al.1994; Duncan & Hughes 1994; Marti et. al. 1995; Aloy et al.2003; Walg et al. 2013). HRSC methods solve the fluid equa-tions in the conserved form. The flow variables in the observerframe, also called ‘state vectors’, are updated in time by eval-uating the fluxes of these state vectors at the cell surfaces.These fluxes are computed by means of exact or approxi-mate Riemann solvers. RHD codes based on exact Riemannsolvers were developed by Marti & M¨uller (1996); Wen et. al.(1997). There are great variety of numerical schemes availableusing approximate Riemann solvers. Balsara (1994) devel-oped an approximate Riemann solver using the jump condi-tions in the shock frame while Dai & Woodward (1997) fol-lowed a similar path but by using conditions for obliquesshocks. Eulderink (1993); Eulderink & Mellema (1995) ex-tended the Roe-type approximate Riemann solver(Roe 1981)to RHD. Higher order reconstructions were also proposed(Marquina et. al. 1992; Dolezal & Wong 1995). More detailedaccount on the development of HRSC methods can be foundin Marti & M¨uller (2003, 2015).A global picture of the flow with multiple shocks within ajet has emerged from these simulations.Most of the simulations were based on the ideal gasequation of state (with a notable exception of Scheck et. al.2002) which is a reasonable approximation when the flowremains sub-relativistic or extremely relativistic, but thejets travel over a long distance and the jet material cango through a transition from the relativistic to the non-relativistic regime and vice versa. It may be noted that aflow can be called relativistic on the account when its bulkvelocity v ∼ c (where c is the speed of light in vacuum) orif the thermal energy of the flow kT ( k is the Boltzmannconstant and T is the temperature) is of the order, or greaterthan the rest mass energy of the gas particles. Assumingrelativistic Maxwell-Boltzmann distribution of particles, theenergy density of the flow or the equation of state (EoS) ofthe fluid, was computed independently by various authors(Chandrasekhar 1938; Synge 1957; Cox & Giuli 1968), whichwere a combination of modified Bessel’s function of variouskinds. This EoS is relativistically perfect and is abbreviated as RP. By using the recurrence relation, Vyas et al. (2015)showed the equivalence between the different forms of theRP EoS, obtained by various authors mentioned above.One of the features of such relativistic EoS is that oneneed not specify any adiabatic index to describe the ther-modynamics of the fluid. The adiabatic index is a functionof temperature and can be automatically obtained if thetemperature is known. Although it is possible to implementthe RP in numerical simulation codes, but the presenceof the modified Bessel function makes it computationallyexpensive (Falle & Komissarov 1996). Additionally, Taub(1948) in his seminal paper of relativistic shock adiabat,also obtained a relation between thermodynamic variablesof the fluid in the form of a fundamental inequality whichany EoS of relativistic fluid should obey, and is called the‘Taub’s inequality’ or TI. Any approximate EoS, apart frombeing a very good fit of RP, should simultaneously obey TItoo. To circumvent the problem of using RP at an extracomputational cost, a number of approximate relativisticEoS were proposed by various authors (Mignone et. al.2005; Ryu et. al. 2006). The EoS used by Mignone et. al. isactually the lower limit of TI. On the other hand, Ryu et. al.(2006) proposed a more accurate approximate EoS, andyet the energy density or the enthalpy of the fluid is analgebraic function of pressure and mass density. One cancompute the adiabatic index (Γ) from this EoS and it gavecorrect asymptotic values at non-relativistic and relativistictemperatures. Chattopadhyay & Ryu (2009), extended theEoS of Ryu et. al. for the fluids composed of dissimilarparticles and the EoS is abbreviated as CR. In this work,we use CR EoS which has been applied to a variety ofastrophysical scenarios (Chattopadhyay & Chakrabarti2011; Cielo et. al. 2014; Chattopadhyay & Kumar 2016;Vyas & Chattopadhyay 2018; Sarkar & Chattopadhyay2019; Sarkar, Chattopadhyay, Laurent 2020;Singh & Chattopadhyay 2019; Dihingia et al. 2018).An astrophysical jet is a relativistic beam plying throughan ambient medium. The beam of the jet, because it is rel-ativistic should be less dense than the surrounding medium.Therefore, the propagation of jets through a medium is essen-tially the time evolution of an initial discontinuity betweentwo states of a fluid which on one side is lighter and fast, anddenser and static on the other side. There might be a pres-sure/composition jump across the initial discontinuity or thepressure/composition may also be uniform. Such an initialvalue problem is generally known as the Riemann problem.Depending on the physical condition of the initial disconti-nuity, a Riemann problem might evolve into a shock-tubeproblem, oppositely moving shock waves or oppositely mov-ing rarefaction waves, a wall shock, etc. In the Newtonianregime, solutions of the Riemann problem played an impor-tant role in testing several hydrodynamic codes (Sod 1978).Most of the modern hydrodynamic codes in the Newtonianregime are based on exact or approximate Riemann problemsolutions (LeVeque 1992), such that building a better Rie-mann solver (exact or approximate), has become a very im-portant field of research (Toro 1997). Codes developed basedon the above techniques and in the Newtonian regime, havebeen used extensively to study diverse fields like the forma-tion of stars in the galactic plane (Kim & Ostriker 2001), su-pernovae ejecta propagation (Arnett et. al. 1989), cosmolog- MNRAS , 1–19 (2020) ical simulations (Ryu et al. 1993), accretion discs (Ryu et al.1995; Ryu et. al. 1997), etc.In 1994, Marti & M¨uller derived analytical solutions of theRiemann problem for special relativistic fluid, but only forthe flow with velocity component normal to the initial dis-continuity. This paper helped to test and even develop vari-ous numerical simulation codes for relativistic fluid (also seeLora-Clavijo et. al. 2013). However, because of the existenceof the upper limit of the fluid velocity, the velocity compo-nents are not entirely independent of each other. As a re-sult, the form of the eigenvalues of a flow with normal andtangential velocity components with respect to the discon-tinuity, are different from that of a flow with only a nor-mal velocity component. Pons et. al. (2000) generalized theanalytical Riemann problem of relativistic fluid obtained byMarti & M¨uller for arbitrary tangential velocity components.The Riemann problem associated with relativistic jets is char-acterised by an FS, a CD or the jet head, and the RS some-where in the jet beam. The solution of Riemann problem witha fixed adiabatic index EoS is relatively easier, but we will dis-cuss later that such a solution is non-trivial with a relativisticEoS like CR. Although, exact solutions of the Riemann prob-lem are generally used to test simulation codes, but theseproblems resemble astrophysical scenarios, so if judiciouslyused, exact solutions of the Riemann problem can be usedto study astrophysical problems as well (Harpole & Hawke2019).In this paper, we solve the Riemann problem associatedwith a relativistic jet. The thermodynamics of the flow is de-scribed by CR EoS. All the basic features of the relativisticjet are discussed considering an electron-proton flow. In thispaper, we would like to find the necessary condition for whichthe initial discontinuity develops into two shocks, the reverseshock in the jet beam and a forward shock in the ambientmedium ahead of the jet head. We would like to investigatehow the expanding jet-head cross-section affects the evolutionof the jet. We also study the effect of jet fluid compositionon the overall jet evolution. We investigate the effect of thecomposition of the jet beam on jet evolution, while keep-ing the ambient composition same. We compared the effectof the ambient medium composition on jet evolution, whilekeeping the jet composition same. For the sake of complete-ness, we have given a short account of the solution of a rel-ativistic shock tube problem and wall-shock problem in theAppendix. We also compared the exact solution of the wall-shock with the relativistic TVD simulation code (Ryu et. al.2006; Chattopadhyay et. al. 2013) in the Appendix.In section 2, we present the governing equations. In section2.1, we describe the CR EoS for fluid composed of dissimilarparticles. In section 2.2, the structure of relativistic shocks isdescribed. We discuss the methodology to obtain the solutionin section 3. In section 4.1, we derive the condition for theformation of two shocks of a jet. In section 4.2, we discuss theeffect of an expanding cross-section of the jet head, on thestructure of the jet. In sections 4.3, 4.3.1, 4.3.2, we discussthe effect of fluid composition on the evolution of jets. Andfinally, in section 5 we discuss the highlights and summarizethe results. In addition, we also present exact solutions oftwo types of relativistic Riemann problem, one is the shock-contact-rarefaction fan problem or a shock-tube problem andtwo, shock-contact-shock with CR EoS in Appendix A.
We study ideal, relativistic fluid in flat space-time, and theenergy-momentum tensor of such a fluid is given by; T µν = ρhu µ u ν + pη µν , (1)where ρ , p , and h are the rest-mass density, local pressure,and the specific enthalpy of the fluid, respectively. We followthe convention where the Greek indices represent space-timecomponents of the vectors and tensors. Contravariant compo-nents of the four-velocity are represented by u µ and η µν arethe components of Minkowski metric tensor in the Cartesiancoordinates. η µν = diag(-1,1,1,1) (2)Four-velocity of the fluid satisfies the normalization conditioni.e. u µ u µ = − ρu ν ) , ν = 0 (4) T µν , ν = 0 (5)Using the normalization condition (3) we can write four-velocity as u µ = γ (1 , v x , v y , v z ) [ v i = v i components of three velocity](6)where γ is the Lorentz factor γ = (1 − v ) − (7)and v = ( v x ) + ( v y ) + ( v z ) (8)In Minkowski space-time, the equations of relativistic hy-drodynamics can be written in conservative form ∂ t U + ∂ i F ( i ) = 0 (9)Where U and F ( i ) ( i ≡ x, y, z) are the vectors and fluxes ofthe conserved variables, respectively. U = ( D, M x , M y , M z , E ) T (10) F ( i ) = ( Dv i , M x v i + pδ xi , M y v i + pδ yi , M z v i + pδ zi , ( E + p ) v i ) T (11)where conserved variables D , M i , and E denote the massdensity, momentum density, and energy density, respectively. Throughout this paper we will use unit system where speed oflight c is set to unity, unless mentioned otherwiseMNRAS , 1–19 (2020) Joshi, Chattopadhyay, Ryu & Yadav
These conserved variables can also be written in terms ofprimitive variables, D = ργM i = ρhγ v i E = ρhγ − p (12)The equation of state (EoS) is used to close the set of equa-tions (9). The EoS can be written in form e = e ( p, ρ ) (13)where e is the energy density in the local frame.The set of equations of motion (4, 5 or 9) are hyperbolic innature and admits five real eigenvalues. Three of which aredegenerate and are the entropy mode, while the first and thelast one are non-degenerate and are the acoustic modes. Theeigenvalues are: β = v x (1 − c s ) − c s √ (1 − v )[1 − v c s − ( v x ) (1 − c s )]1 − v c s β = v x β = v x β = v x β = v x (1 − c s )+ c s √ (1 − v )[1 − v c s − ( v x ) (1 − c s )]1 − v c s (14)Here c s = p ( ∂p/∂e ) s is the adiabatic, relativistic soundspeed. The full eigenstructure with right and left eigenvec-tors for relativistic fluid has also been obtained previouslyfor a general EoS (see Ryu et. al. 2006). We use CR EoS (Chattopadhyay & Ryu 2009) for the fluidscomposed of electrons, positrons, and protons. The CR EoSis of the following form e = Σ i (cid:18) n i m i c + p i p i + 3 n i m i c p i + 2 n i m i c (cid:19) (15)In Eq. 15, c is used explicitly. The index i represents thevarious species that constitute the fluid and c is the speedof light in vacuum. In this paper, we consider the fluid tobe composed of electrons, protons, and positrons of variousproportion. The above equation can be represented in theunit system where c = 1, as, e = ρf, (16)where, f = 1 + (2 − ξ )Θ (cid:20)
9Θ + 6 /τ
6Θ + 8 /τ (cid:21) + ξ Θ (cid:20)
9Θ + 6 /ητ
6Θ + 8 /ητ (cid:21) (17)In the above equations ρ = Σ i n i m i = n e − m e − (2 − ξ + ξ/η ),where ξ = n p /n e − , η = m e − /m p and n e − , n p , m e − and m p are the electron number density, the proton number density,the electron rest mass, and proton rest mass. Moreover, theratio of the pressure and the local rest energy density of fluid, is a measure of temperature Θ = p/ρ and τ = 2 − ξ + ξ/η .The specific enthalpy is given as h = ( e + p ) /ρ = f + Θ; (18)The expression for polytropic index N is given as N = ρ ∂h∂p − ∂f∂ Θ = 6 h (2 − ξ ) +24Θ /τ +8 /τ (6Θ+8 /τ ) i (19)+6 ξ h +24Θ / ( ητ )+8 / ( ητ ) { / ( τη ) } i The polytropic index is not a constant but a function ofΘ and ξ . It may be noted that, N → ≫
1; while N → / ≪
1, therefore N approaches asymptoticvalues at very high and low temperatures. It may be notedfurther that, for ξ = 0 (i. e., single species gas), the expressionof the polytropic index is exactly same as that was presentedin Ryu et. al. (2006). And the adiabatic index Γ isΓ = 1 + 1 N (20)The variation of adiabatic index Γ with respect to Θ is shownin Fig. (1). The composition of the fluid is marked in thelegend, electron-positron or ξ = 0 (solid, blue), equal propor-tion of positrons and protons or ξ = 0 . ξ = 1 . >
10 is the ultra-relativistic tempera-ture regime i. e. Γ → /
3, for a fluid of any composition ξ .On the other hand, Θ < − is the non-relativistic temper-ature regime, where Γ ∼ / ξ >
0. But for ξ = 0, Γ → / < ∼ few × − . Clearly, the thermo-dynamics of the flow depends on both the temperature andthe composition ξ of the flow. It may be remembered thatthe conversion between the absolute temperature T and Θis given by T = τ m e − Θ / k . The expression of sound speedstarting from the first principle can be written as, c s = 1 h ∂p∂ρ = − ρNh ∂h∂ρ = ΓΘ h (21)For non-relativistic temperatures or Θ ≪ h → → /
3, sound speed approaches non-relativistic value c s → p /
3. In contrast, for Θ ≫ h ∼ ∼ /
3, so c s → / √
3, once again the sound speed achieves asymptoticvalues for ultra-relativistic and non-relativistic temperaturelimit.
The information about the states on both sides of a shockis obtained by the jump conditions based on the continu-ity of mass and energy-momentum fluxes. These conditionsare known as the Rankine-Hugoniot (RH) conditions. Therelativistic version of RH conditions was obtained by Taub(1948). The relativistic RH conditions are given as (also seePons et. al. 2000)[ ρu µ ] n µ = 0 (22)[ T µν ] n ν = 0 (23) Although the EoS in equation 16 is exactly same as that inChattopadhyay & Ryu (2009), but in this paper, τ is included inthe definition of Θ and f MNRAS , 1–19 (2020) −6 −5 −4 −3 −2 −1 Θ Γ Γ = 5/3 Γ = 4/3 ξ = 0.0ξ = 0.5ξ = 1.0
Figure 1.
Variation of adiabatic index with Θ for fluids with com-position parameter ξ = 0 . , ξ = 0 . , ξ =1 . where n µ is a unitary normal vector to shock surface Σ andwe have used the notation. Here,[ F ] = F a − F b (24) F a and F b are the values of the function F on the either sidesof the shock surface Σ. Considering Σ to be normal to the xaxis and using the unitarity of n µ , we can write it as n µ = γ s ( V s , , , , (25)where V s is the speed of shock (the speed of surface Σ) and γ s is the Lorentz factor of the shock. γ s = 1 p − V s2 (26)We can introduce the invariant relativistic mass flux j acrossthe shock as j = γ s D a ( V s − v xa ) = γ s D b ( V s − v xb ) (27)In the above the suffix ‘a’ and ‘b’ denoted the quantites oneither side fo the shock. The positive (negative) value of j represents the shock propagating to right (left). Multiplyingequation (23) with n µ and using the definition of j (equa-tion 27), we obtain the expression for relativistic mass flux interms of pressure, enthalpy and density j = − [ p ][ h/ρ ] (28)We can write the RH conditions equations (22) and (23) interms of conserved quantities D, M j , and E as[ v x ] = − jγ s (cid:20) D (cid:21) (29)[ p ] = jγ s (cid:20) M x D (cid:21) (30) (cid:20) M y D (cid:21) = 0 (31) RS CD FS42 5 ( (cid:1) ,p ,v ) Figure 2.
Cartoon diagram of the one dimensional jet (1DJ) struc-ture. The jet is divided into six regions, where region 1 is the orig-inal beam of the jet, and region 6 is the original ambient medium.Region 2 is the location of reverse shock (RS), and region 5 is thelocation of the forward shock (FS). The CD is the location of thejet head (JH). The region 3 is the shocked jet region and region 4is the shocked ambient medium. Generally the cross-section mayvary across every region. (cid:20) M z D (cid:21) = 0 (32)[ v x p ] = jγ s (cid:20) ED (cid:21) (33)Equations (31) and (32) imply hγv y,z = constant (34) v y,zb = h a γ a v y,za (cid:20) − ( v xb ) h b + ( h a γ a v ta ) (cid:21) / (35)Where v t is the absolute value of tangential velocity of theflow. v t = p ( v y ) + ( v z ) (36)The expression for normal flow velocity v x is obtained byusing equations (29), (30) and (33) v xb = (cid:16) h a γ a v xa + γ s ( p b − p a ) j (cid:17)h h a γ a + ( p b − p a ) (cid:16) γ s v xa j + ρ a γ a (cid:17)i (37)The expression for the shock velocity is obtained using thedefinition of mass flux V s ± = ρ a γ a v xa ± | j | p j + ρ a γ a (1 − ( v xa ) ) ρ a γ a + j (38) V s+ ( V s − ) corresponds to shock propagating towards right(left).Multiplying equation (23) first by ( hu µ ) a and then by( hu µ ) b and adding both the expressions results in[ h ] = (cid:18) h b ρ b + h a ρ a (cid:19) [ p ] (39)Equation (39) is known as Taub adiabat. For the generalcase Taub adiabat is solved along with the EoS to obtain h b as a function of p b . A jet is characterized by low density, high velocity (high v x ),‘geometrically narrow’ material travelling with relativistic ve-locity through a denser, colder, static medium. We consider MNRAS , 1–19 (2020)
Joshi, Chattopadhyay, Ryu & Yadav xt ρ , p , v , p , v *3 ,p *3 ,v * ρ *4 ,p *4 ,v * (a) xt ρ , p , v , p , v *3 ,p *3 ,v * ρ *4 ,p *4 ,v * (b) Figure 3.
Schematic diagram for: (a) relativistic jet (1DJ) — At t = 0, region 1 is the beam of the jet, and region 6 is the ambi-ent medium. The initial discontinuity evolves into a forward shock(FS), and contact discontinuity (CD) and a reverse shock (RS). (b)Shock tube problem — Region 1 and 6 are the initial left and rightstate and are marked as suffix of the flow variables in the respec-tive regions. The initial discontinuity evolves into a CD flanked byan FS and an RF. the relativistic jet as the time evolution of a one dimensional,initial discontinuity. The defining distinction being that theinitial left state of the fluid is considered as the jet beam,which is lighter than the ambient medium, moving with rela-tivistic speed to the right. The surging jet beam drives a shockon the ambient medium which also moves in the direction ofpropagation and is called the forward shock (FS). The initialdiscontinuity evolves with local jet speed and is the jet-head(JH) or contact discontinuity (CD). The CD moves slowerthan the jet beam and drives a shock in the jet beam and iscalled the reverse shock (RS).Therefore, one dimensional jet (1DJ) consists of a forwardand a reverse shock separated by a CD. The solution is rep-resented as six regions (see Fig. 2), which are: • Region 1: It is the initial unperturbed state which ischaracterised by injection parameters of jet beam, ρ , p and v (i. e., v x ). • Region 2: is a thin surface of RS, which may move in thesame direction as FS, but depending on the physical conditionmay also move in the opposite direction. • Region 3: It is the region between the RS and the CD.The state is represented as ρ ∗ , p ∗ and v ∗ . • Region 4: It is the region between CD and FS. There is adensity jump across CD and the flow variables in this regionare represented as ρ ∗ , p ∗ and v ∗ . • Region 5 is the surface of FS travelling to the right. • Region 6: This region corresponds to the ambientmedium which has not been influenced by the forward shockwave and the flow variables in this region are represented as ρ , p , & v .It may be noted that regions 2 and 5 are thin shock surfaces.However, depending on the initial conditions both thesesurfaces may become rarefaction fan with non-negligiblespatial extent (see, Figs. 3a,b for qualitative argument andFigs. 4a-d for a quantitative understanding). Therefore as apart of general approach, these surfaces are called regions.In Fig. (3a), we present the schematic diagram of the time evolution of these surfaces in 1DJ. Figure (3a) shows thatthe location of CD, FS, and RS at any time will be on thethree straight lines and the values of the flow variables in thesix regions will be as mentioned above. It is clear that thetime evolution of 1DJ is similar to an initial value problemwhich is also known as the Riemann problem. Riemannproblems are of various kinds, one of the most popularexample of one, is the so-called shock-tube (ST) problem.In Fig. (3b), we present the schematic diagram of an STproblem for comparison. It may be remembered that ST toois an initial value problem, however, the initial conditionsare such (high pressure/density) that FS and CD forms butinstead of an RS, a rarefaction fan (RF) is formed. The RFmoves in a direction opposite to FS and CD. The detailedsolution procedure for a couple of other Riemann problemsis presented in the Appendix. In this paper, we present thesolution for 1DJ problem.The time evolution and solution of the 1DJ problem, de-pends crucially on the measurement of the post shock veloc-ity. Recalling equation (37) we obtain an expression for ( v x ) ∗ in terms of v x for a RS and an expression for ( v x ) ∗ in termsof v x for a FS.( v x ) ∗ = (cid:18) h γ v x + γ s ( p ∗ − p ) j (cid:19)(cid:18) h γ + ( p ∗ − p ) (cid:18) γ s v x j + 1 ρ γ (cid:19)(cid:19) − (40)( v x ) ∗ = (cid:18) h γ v x + γ s ( p ∗ − p ) j (cid:19)(cid:18) h γ + ( p ∗ − p ) (cid:18) γ s v x j + 1 ρ γ (cid:19)(cid:19) − (41)For equation (40) j is the negative root of equation (28) andfor equation (41) it is the positive root of equation (28). Thecontinuity of velocity across JH is the key to obtain the solu-tion of this problem( v x ) ∗ − ( v x ) ∗ = 0 (42)Equation (42) is solved for p ∗ using the iterative root finderand rest of the quantities can be calculated once p ∗ isobtained. Densities in Region 3 and 4 are obtained fromthe Taub adiabat for RS and right FS respectively. The( v x ) ∗ = ( v x ) ∗ = ( v x ) ∗ is the speed V j with which the jet-heador CD is propagating.xThe continuity of the normal component of velocity fol-lows from the fact that the mass flux across the CD is zero. j = γ j D a ( V j − v xa ) = γ j D b ( V j − v xb ) = 0 (43)Where V j is the velocity of jet head or the CD and γ j is theLorentz factor of jet head. As the jet expands it can result in discontinuous cross sec-tional area across the jet head. To obtain the equations thatgovern the dynamics of the flow when the area across CD/JH
MNRAS , 1–19 (2020) is not same, we use the momentum flux balance across theCD/JH (Mizuta et. al. 2004) A b (cid:2) ρ b h b γ j γ b ( v b − V j ) + p b (cid:3) = A a (cid:2) ρ a h a γ j γ a ( v a − V j ) + p a (cid:3) (44)Using equation (43) with (44) we obtain A b p b = A a p a (45)And the velocity balance condition across the jet head is givenas v ∗ ( p ∗ ) − v ∗ ( p ∗ ) = 0 (46)Assuming p ∗ = p ∗ , equation (46) is a transcendental equationfor variable p ∗ , as p ∗ and p ∗ are related by equation (45) as A p ∗ = A p ∗ .The jet kinetic luminosity is related with the jet cross-section L j = γ ( e + p ) v πy , (47)Where quantities with suffix ‘1’ are the jet beam variableswhere the jet beam velocity v = v x , and y is the cross-sectional dimension of the jet beam, so A = πy . The essential condition for the formation of two shock frontsis that the pressure in the intermediate state p ∗ (refer to Fig3) should be greater than the pressure in the initial states.Rezzolla & Zanotti (2001) have shown that the intermediatepressure is a function of relative velocities between the initialstates so there will be a limiting value of relative velocity forthe occurrence of two shock fronts. In this case, we assume thecross-sectional area of the flow to be invariant across all theregions. To obtain the value of limiting velocity we start byassuming the scenario with a left and right shock separated bycontact discontinuity. We assume the pressure of jet beam p is greater than the pressure of ambient medium p . Assumingthere is an RS, the relative velocity of the pre-shock region (1)and the post-shock region (3) is given by (Rezzolla & Zanotti2001) v = s ( p ∗ − p )( e − e )( e + p ∗ )( e + p ) (48)Similarly, the relative velocity (between region 4 and 6) aheadof FS is given as v = − s ( p ∗ − p )( e − e )( e + p ∗ )( e + p ) (49)It may be remembered, v = v = v ∗ across the CD. Hencethe relative velocity of between two initial states is given as( v ) S = v − v − v v (50)For the limiting case where there exist an RS, the minimumvalue of pressure in the intermediate region i.e. p ∗ , is givenas p ∗ = max(p , p ) = p (51) Hence the limiting value of relative velocity is given as v lim = s ( p − p )( e − e )( e + p )( e + p ) (52)As the ambient medium is at rest, therefore the rel-ative velocity between the initial states is ( v ) S =( v − v ) / (1 − v v ) = v . For a particular density and pres-sure configuration, if the injection velocity v is greater thanthe limiting velocity v lim (see, equation 52), then the injectedjet beam evolves with an FS and an RS, otherwise, it evolvesas a standard shock-tube problem with an FS and an RFseparated by CD. In Fig. (4a-d) we plot flows with two injec-tion velocities, one with v < v lim (Fig. 4 a, c) and two with v > v lim (Fig. 4b, d). The initial parameters are given by p = 100 , ρ = 0 . , p = 1 . ρ = 100 , and v = 0 (53)The composition parameter of beam and ambient mediumis same ( ξ = 1 . v lim = 0 .
66. The thresholdlevel of injection velocity is represented by the dashed line inFigs. (4c & d). Figure (4a & c) show the variation of densityand velocity, respectively, as functions of position with beaminjection velocity v = 0 .
60. In Figs. (4 b & d) the variationof ρ and v x are for the case v = 0 .
70. It is clear that theinitial discontinuity with the injection v < v lim , evolves intoa shock tube problem with an CD flanked by an FS aheadand an RF behind it (Fig. 4a & c). The RF is zoomed in Fig.(4a) for clarity. The CD and FS going to the right while RFto the left. On the other hand, in Fig. (4b & d) v > v lim , theinitial discontinuity evolves into a solution which has a CDflanked by an FS and RS. The RS is zoomed in the inset ofFig. (4b). For a pressure matched jet ( p = p ) the threshold value ofrelative velocity from equation (52) is obtained as v lim = 0 (54)Hence for any non zero injection velocity v , the pressurematched jet will always evolve as two shock fronts separatedby a contact discontinuity. Figure (5) shows the flow vari-ables for a jet beam with initial parameters ρ = 1 . , p =1 . , v = 0 . ρ = 100 . , p = 1 . , v = 0 .
0. We have assumed the uni-form cross sectional flow for this case.The density jumps in Fig. (5a), show the positions of RS,CD, and FS (from left to right, respectively). The shocked re-gion of the jet beam is in between CD and RS and the kineticenergy of the jet beam is converted to thermal energy in thepost-shock region of RS. This shock heated region containsthe most relativistic gas (Γ ∼ / As the jet ploughs through the ambient medium, the jet headcross-section may expand. To employ the expansion effect we
MNRAS , 1–19 (2020)
Joshi, Chattopadhyay, Ryu & Yadav −1 ρ (a) (b) x v x (c) x (d) Figure 4.
Comparison of a shock-tube (a, c) and 1DJ solution (b, d). The injection velocities are v = 0 . v = 0 . p = 100 , ρ = 0 . , p = 1 . , ρ = 100 and the location of the initial discontinuity is at x = 0 .
25. The plot is at time t = 0 .
4. The composition parameter of injected beam and ambient medium is taken to be ξ = 1 . ρ (a) p (b) x v x (c) x Γ (d) Figure 5.
Flow variables (a) ρ ,(b) p , (c) v x and (d) Γ as functions of x at t = 1 . ξ = 1 .
0. Initial conditions: ρ = 1 . , p = 1 . , v = 0 . ρ = 100 . , p = 1 . , v = 0 . x = 0 . , 1–19 (2020) t A (a) t −0.2−0.10.00.10.20.30.4 v (b) V FS V j V RS t p (c) t ρ (d) Figure 6. (a) Variation of area ratio across the jet head with respect to time. (b) Variation of the RS, JH and FS velocities as functionsof time for an expanding electron-proton jet with injected Lorentz factor 10 and p = p = 0 . , ρ = 1 . , ρ = 500 . x = 0 .
25. Thevariation of pressure and density in region 4 which is the post shock region of FS is shown in panels (c) and (d) respectively. ρ (a) t=1.0t=4.0t=7.0t=10.4t=10.6 −1 p (b) x v x (c) x Γ (d) Figure 7.
Flow variables ρ ( a ) , p ( b ) , v x ( c ) and Γ( d ) as functions of x at different time for an expanding electron-proton jet. The jet solutionare for the same parameters as in Figs. (6). assume that the area across JH varies in time given by afunction A = 1 + C exp ( t − t )1 + exp ( t − t ) (55)Where C is a constant which sets up and upper limit on A and t limits the onset time for expansion. Figure (6a) isplotted for C = 20 . t = 15 .
0. From equation 45 onecan conclude that across the JH p A = p A (56)Hence the pressure for an expanding jet also exhibits a jumpacross the CD unlike the Riemann problem where the areaacross the CD was same. It may however be noted that some of the additional effects like back flow of jet material andgeneration of transverse flow components from the region be-tween RS and CD has not been considered for simplicity.Figure (6a) shows the variation of area across JH ( A ) withrespect to time t . One may refer to Fig. (2), where we haveshown expansion across the jet head, we have assumed thebeam area A to be unity and A = A , also the area acrossthe reverse shock surface is same on both sides ( A = A ). Inall other sections, we have assumed invariant jet cross-section.In Fig. (6b) we have plotted the propagation velocities V RS (dashed-dotted, blue), V j (solid, red) and V FS (dashed, black)for a pressure matched electron-proton jet propagating in ahomogeneous electron-proton ambient medium. The varia- MNRAS , 1–19 (2020) Joshi, Chattopadhyay, Ryu & Yadav
100 200 300 400 500 600 Θ /Θ t rev (a) p = 0.1, ρ = 1.0p = 0.01, ρ = 0.1 0 2 4 6 8 10 x −1 ρ (b) p = 0.1, ρ = 1.0, ρ = 300.0p = 0.01, ρ = 0.1, ρ = 30.0 Figure 8. (a) We plot the time t rev taken by RS to move in the opposite direction of JH or CD as a function of Θ / Θ . (b) The densityprofile ρ w.r.t x of the two jets marked by ⋆ ( ≡ t = 12 .
08) in panel (a). The injection speed of both the electron-proton ( ξ = 1) jets v = 0 . p = p = p = 0 . ρ = 1 for one jet model (solid, blue) and p = p = p = 0 .
01 & ρ = 0 . ρ for the jets, scale according to the ratio Θ / Θ . Θ /Θ t rev ξ = 0.0ξ = 0.5ξ = 1.0 Figure 9.
Variation of t rev as a function of Θ / Θ for expandingjets with v = 0 . , p = p = 0 . , ρ = 1 .
0, but for composition ξ = 0 (dashed, red), ξ = 0 . ξ = 1 . tion of pressure and density in post-shock region of FS (re-gion 4) with respect to time also shown plotted in panel (c)and (d) respectively. Initial discontinuity is at x = 0 .
25 andthe initial flow parameters are given as ρ = 1 . , ρ = 500 . , p = p = 0 . , v = 0 . , v = 0 . V RS may become nega-tive or in other words, the RS may start to move backward, asis shown in Fig. (6b), where the RS starts to move backwardat t ∼ .
4. The corresponding solutions ( ρ, p, v x , & Γ) atdifferent epochs are presented in Figs. (7a-d). The overall jetstructure advances forward (to right), however at t > . V RS <
0. Therefore Fig. (7c) shows thatthe RS moves forward from t = 0 → .
4, but it movesbackward t = 10 . → . t = 10 . t = 10 . F S and CD continue to move forward. For the initial stages ofjet propagation ( t . p and initial density contrast ρ /ρ would give different values of V j , V FS and V RS . It isinteresting to note that certain physical condition of the en-vironment of the jet may cause the RS to revert back! In Fig.(8a), we plot the time t rev at which the RS starts to move inthe reverse direction as a function of initial ratio of Θ / Θ ,for two electron-proton jets with the same injection Lorentzfactor γ = 10 (i. e., v = 0 . p = 0 . , & ρ = 1 . p = 0 .
01 & ρ = 0 . / Θ , initial ambient density ρ scales similarlyfor both the jets but have different values. For example, thejet with p = 0 . / Θ = 300 implies ρ = 300while for the jet p = 0 .
01 (open circle, red), the same Θ ratioimplies ρ = 30. In this figure, we have considered Θ / Θ ratio from 50 to 600. The two curves exactly match for theentire range. So even if the pressure and the density of thejets are not similar, the RS reverts back at the same timefor a given Θ / Θ ratio. It may be noted that, there existsa limiting value of Θ / Θ ( ∼ t rev decreases drastically. For any value of Θ / Θ greater than this limiting value, the RS will move opposite tothe direction of propagation from the start. In Fig. (8b), wecompare ρ as a function of x of the jet cases, for parameterscorresponding to the black star on the curve in Fig. (8a). Al-though the density distribution differs between the two jets MNRAS , 1–19 (2020) being considered here, but the location of RS, CD and FScoincide. Therefore even if the ρ and p are different for thesame ratio of the Θ, the jet evolution is exactly similar. Itmay however be noted the with different p and ρ , the jet ki-netic luminosity L j differs along the curve. In contrast, if thecomposition parameter is varied for the jet and the beam, t rev is different even if v and Θ / Θ are the same. In Fig. (9), weplot the time of reversal of the RS t rev with the ratio Θ / Θ for jets with v = 0 . , p = p = 0 . , ρ = 1 .
0, howeverfor jets for three composition parameter ξ = 0 (dashed, red), ξ = 0 . ξ = 1 . ξ . The thermal state of the fluid depends upon the ratio
T /m ( m is the mass of the constituent particles of gas) hence thechange in composition parameter will affect the thermal state.To study the effect of the composition parameter on jet mor-phology we study two scenarios. In the first case, we keepthe composition parameter of ambient medium and jet beamsame and for the second scenario, we study the case in whichthe composition of ambient matter differs from the jet com-position. To study the effect of composition we have assumedthe uniform cross-sectional flow in all the cases. Figure (10) compares the solution snapshots of a relativisticjet at t = 1 . ξ = 0 . ξ = 0 . ξ = 1 . ρ = 1 . , ρ = 100 . , p = p = 1 . , v = 0 . ξ = 0 . ξ = 0 . ξ = 0 . ξ = 0 . Figure (11) compares the snapshot of flow variables at t = 1 . ξ = 1 .
0) is injectedinto ambient medium with different composition parameters like ξ = 0 (dashed, red), ξ = 0 . ξ = 1 . ρ = 1 . , ρ = 100 . , p = p = 1 . , v = 0 . ξ = 1 . ξ = 0, has thewidest and tallest high pressure region bounded by FS andRS. It also has the widest high density shell (the region be-tween CD and FS), although the height of the shell is small-est compared to the other two cases. The FS travels withthe highest velocity for ξ = 0 because of the low resistanceoffered by the medium. The shock heating is maximum forelectron, positron & proton plasma resulting in the lowestadiabatic index in comparison to others. Since the jet beaminitial conditions are exactly same in all the three cases, theambient medium with ξ = 0 is the least hot (Γ ∼ /
3) andtherefore offers the least resistance to the jet, which causesthe difference in the jet structure as well.We would now study the effect of composition of a jet beamflowing through the same ambient medium. We consider thecold ambient medium of a static electron-proton fluid. Figure(12) compares the flow variables for jets of different composi-tion traveling in the ambient medium composed of electronsand protons. The composition of the jet beam are ξ = 0 . . . ρ = 1 . , p = 0 . , ρ = 100 . , p = 0 . , v = 0 . ξ = 1 . ρ and p values. The difference inRS and FS velocities determines the size of post shock re-gion. The post shock region in blazar jets is considered as theflaring region. This region is broadest for electron-positronjet. The composition of the jet beam can also affect the directionof propagation of the reverse shock. In Fig. (13) we com-pare the flow variables of a pressure matched jet beam withLorentz factor of 10 injected into an electron-proton ambientmedium. The initial flow parameters are given as ρ = 1 . , p = 0 . , ρ = 100 . , p = 0 . , v = 0 . x = 0 .
25. Figure (13) showsthat for pure electron-positron jet beam i.e. ξ = 0 . V RS (Fig. 14a), V j (Fig. 14b) and V FS (Fig. 14c)with the composition parameter of the jet beam ξ . For theseflow parameters, the reverse shock moves in the opposite di-rection for the pair plasma jet. However, V RS increases rapidlywith the addition of protons and maximizes for ξ ∼ . V RS < ξ > ∼ .
24. The effect of composition on
MNRAS , 1–19 (2020) Joshi, Chattopadhyay, Ryu & Yadav ρ (a) ξ=0.0ξ=0.5ξ=1.0 p (b) x v x (c) x Γ (d) Figure 10.
Flow variables ρ ( a ) , p ( b ) , v x ( c ) and Γ( d ) as functions of x at t = 1 .
0. The composition parameter of ambient medium is sameas that of jet beam and mentioned in the legends in panel (a) . The initial condition is given by ρ = 1 . , ρ = 100 . , p = p = 1 . , v =0 . x = 0 . various propagation velocities is quite clear, although the ef-fect is more pronounced on V RS . Figure (14b) shows that forthese parameters the electron-positron jet is not the fastestjet, instead, a jet with composition ξ ∼ .
18 is the fastest.Finally, in Fig. (15), we plot the injection speed v as afunction of Θ / Θ for jets of various jet composition ξ = 0(solid, red), ξ = 0 . ξ = 1 . v for different Θ / Θ , below which rar-efaction fan (RF) form, and above which reverse shock RSforms but with negative velocity. The upper curve representsthe locus of the injection speed above which the RS moveswith positive velocity. Assuming the width of the beam y as1 pc, and the ρ as 10 − g cm − , we used equation (47) tomap the jet luminosity L j (right label) from the values of v .The limiting velocity v lim does not depend on ξ which is thelower curve and therefore coincides for all ξ . However, theupper curve, which is the limiting injection speed for V RS > In this paper, we obtained the exact solution of time-dependent, one-dimensional relativistic jets plying throughexternal media (or ambient media), by solving equations forthe ideal special relativistic fluid composed of dissimilar par-ticles and the thermodynamics of the fluid governed by arelativistic equation of state. The relativistic EoS we used iscalled the CR EoS, is an approximation of the ChandrasekharEoS. It is difficult to obtain a general Riemann solution witha general EoS and has been discussed in the introduction. Forthe sake of completeness, we presented the solution method- ology and the solutions for three Riemann problems — rel-ativistic shock tube, shock-contact-shock problem and thewall-shock problem in the appendix. The intended impactof this work is firstly to provide exact solutions for numericalsimulation codes to compare their results with. But more thanthat, we studied various properties of the time-dependent rel-ativistic jets, which are in general glossed over. We addressedthe issue of the properties of reverse shock. When a reverse-shock will form and when will it not. When the reverse shockmove in the opposite direction to the jet head and when alongwith it. We wanted to study the effect of expanding jet-headcross-section on the jet solution. We also wanted to study theeffect of the composition of the jet or the ambient medium onthe jet solutions and on all these counts we did get interest-ing results. In the process, we modified the method to solvefor jet-like Riemann problem to incorporate the cross-sectioneffect across the jet-head. However, it should be noted thatmany features of a relativistic jet, like generation of trans-verse velocity near the jet head, back flow etc are not in-corporated in this model for simplicity. Such omission mayproduce a higher estimate of jet propagation speed, althougha suitable modeling of the expanding jet geometry may pro-duce reasonable estimates. The change in cross-section acrossthe jet-head is an approximate attempt to mimic the tran-verse expansion of the jet and has been discussed in onepart of this paper (section 4.2). Such modifications acrossCD have been attempted earlier (Begelman & Cioffi 1989;Mizuta et. al. 2004).In this paper, we obtained the condition for the formationof reverse shock and showed that it crucially depends on theinjection speed and the initial density contrast between thejet and the ambient medium, although for pressure matchedjet that critical injection speed turns out to be zero, in other
MNRAS , 1–19 (2020) ρ (a) ξ =0.0ξ =0.5ξ =1.0 p (b) x v x (c) x Γ (d) Figure 11.
Flow variables (a) ρ , (b) p , (c) v x and (a) Γ as functions of x at t = 1 . ξ = 0 . , ξ = 0 . , ξ = 1 . x = 0 . ξ = 1 . words, in a pressure matched jet a reverse shock will alwaysform. However, if the jet-head cross-section is larger thanthe beam cross-section, then the reverse shock may travel inthe reverse direction. For the particular case of a pressurematched jet which is launched with a Lorentz factor 10 withan initial density contrast of 500, if the cross-section acrossthe jet head remains invariant for some initial time and thenexpands for some time and settles into a constant value, thenwe showed that initially RS, CD, and FS moved in the pos-itive direction, after a code time t > ∼ .
4, RS reversed back.And then as the cross-section settles to a certain uniformvalue, then the RS also settles to a uniform negative veloc-ity. We also showed that the time of reversal of RS t rev isunique for a given ratio of the initial Θ of the beam and theambient medium, if the injection speed and composition issame. However, for different values of ξ , t rev will be differenteven if v and initial Θ ratio is same. Therefore, ξ affects thejet evolution. In general, the propagation speed depends onthe beam injection velocities, jet-head cross-section as wellas initial density contrast across the surface of separation.We also showed that the composition of the jet also affectsthe jet structure as it evolves in time, even if all the initialmacroscopic fluid variables are same. The composition of thejet plasma affects the jet flow variables like density and pres-sure and also the propagation velocities like V RS , V j and V FS too. If jets of different composition flow through an ambi-ent medium of a certain composition, then the jet with lessbaryons are pushed back. Depending on initial conditions ifthe jet beam and jet head cross-sections are different then,we showed that for both baryon poor and baryon loaded jets,the reverse shock can go back, even though the jet-head ismoving in the positive direction. The composition of the jetnot only quantitatively affect the jet flow like different den- sity, pressure and velocity distribution, but can also producequalitatively different effects like shocks moving in oppositedirection for one composition and shocks moving in the samedirection for another composition. We also estimated the in-jection speed for an initial pressure mismatch jet, for whichrarefaction fan (RF) will form instead of the reverse-shockRS and also the upper limit of injection speed for which theRS has negative velocity. Assuming the jet cross section ofthe order of 1 pc, these injection speeds correspond to a jetkinetic luminosity of around 10 erg s − . Interestingly thisupper limit of injection speed depends on the compositionparameter even for uniform cross-section jet. Our results in-dicate that in FR II jets typically with jet kinetic luminositiesabove 10 erg s − , both FS and RS would move forward. Onthe other hand, in FR I jets with lower luminosities, eitherbackwardly moving RS or RF could form. Hence, the numer-ical simulations of FR I jets need special care, partly becauseof it, but also because other processes including the entrain-ment of ambient material seems to become important (e. g.,Perucho & Mart´ı 2007). DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
ACKNOWLEDGMENTS
The work of DR was supported by the NationalResearch Foundation (NRF) of Korea through grants
MNRAS , 1–19 (2020) Joshi, Chattopadhyay, Ryu & Yadav ρ (a) ξ =0.0ξ =0.5ξ =1.0 p (b) x v x (c) x Γ (d) Figure 12.
Flow variables ρ ( a ) , p ( b ) , v x ( c ) and Γ( d ) as functions of x at t = 1 . ξ =0 . , ξ = 0 . , ξ = 1 . ξ = 1 .
0. The initial discontinuity is at x = 0 . ρ (a) ξ =0.0ξ =0.005ξ =0.2 p (b) v x (c) Γ (d) Figure 13.
Flow variables ρ ( a ) , p ( b ) , v x ( c ) and Γ( d ) as functions of x at t = 1 . REFERENCES
Aloy M.- ´A., Mart´ı J.-M., G´omez J.-L., Agudo I., M¨uller E., Ib´a˜nezJ.-M., 2003, ApJL, 585, L109Arnett D., Fryxell B., Muller E., 1989, ApJ, 341, L63.Balsara D.S., 1994, Journ. Comput. Phys., 114, 284Begelman M. C., Cioffi D. F., 1989, ApJL, 345, L21MNRAS , 1–19 (2020) ξ -0.030-0.020-0.0100.0000.010 V RS (a) ξ V j (b) ξ V FS (c) Figure 14.
The propagation velocities of (a) the RS i. e., V RS , (b) the CD i. e., V j and (c) the FS i. e., V FS are plotted as a function of ξ at t = 1 . Θ /Θ ξ = 0.0ξ = 0.5ξ = 1.0 L j ( × e r g ss − ) Figure 15.
Injection speed v ( L j ) as a function of Θ / Θ , for ρ = 1 . p = 1 . p = 0 .
1. The curves corresponds beamcomposition ξ = 0 (solid), ξ = 0 . ξ = 1 . ξ = 1 .
0. The top curves (red) are upper value of v (or L j )for which the RS goes in the opposite direction to the JH and thebottom curves (blue and merged for all ξ ) represent the upperlimit of v (or L j ) for which RF is formed instead of RS.Blandford R. D., Rees M. J., 1974, MNRAS, 169, 395Chandrasekhar S., 1938, An Introduction to the Study of StellarStructure . Dover, New York.Cox J. P., Giuli R. T., 1968,
Principles of Stellar Structure, Vol.2.
Gordon and Breach Science Publishers, New York.Chattopadhyay I., Ryu D., 2009, ApJ, 694, 492Chattopadhyay I., Chakrabarti S. K., 2011, IJMPD, 20, 159Chattopadhyay I., Ryu D., Jang H., 2013, ASInc, 9, 13.Chattopadhyay I., Kumar R., 2016, MNRAS, 459, 3792Cielo S., et. al., 2014, MNRAS, 439, 2903.Curtis H. D., 1918, Lick Obs. Publ., 13, 31Dai W., Woodward, P. R., 1997, SIAM J. Si. Stat. Comput., 18,982Dihingia I. K., Das S., Maity D., Chakrabarti S., 2018, PhRvD,98, 083004 Dolezal A., Wong S. S. M., 1995, Journ. Comput. Phys., 120, 266Doeleman S. S., et al., 2012, Science, 338, 355Duncan G. C., Hughes P. A., 1994, ApJ, 436, L119Eulderink F., 1993, Numerical Relativistic Hydrodynamics, PhDThesis, (Rijksuniverteit te Leiden, Leiden, Netherlands)Eulderink F., Mellema G., 1995, A&AS, 110, 587Fender R. P., Gallo E., Russell D., 2010, MNRAS, 406, 1425Falle S. A. E. G., Komissarov S. S., 1996, MNRAS, 278, 586Gallo E., Fender R. P., Pooley G., 2003, MNRAS, 344, 60Harpole A., Hawke I., 2019, ApJ, 884, 110Junor W., Biretta J. A., Livio M., 1999, Nature, 401, 891Kamm, J. R. 2015
An exact, compressible one-dimensional Rie-mann solver for General, convex equations of state , LossAlamos Report.Kim W. T., Ostriker E. C., 2001, ApJ, 559, 70.LeVeque R. J., 1992
Numerical Methods for Conservation Laws ,2nd edn., Birkh¨auser.Lora-Clavijo F. D., Cruz-Perez J. P., Guzman F. S., Gonzalez J. A.,2013, Revista Mexicana de Fisica E, 59, 28.Marquina A., et. al., 1992 A& A, 258, 566,Marti J. M., Muller E., 1994, JFM, 258, 317Marti J. M., Muller R., 1994, Ibanez J. M., 1994, A&A, 281, L9Marti J. M., M¨uller E., Font J. A., Ibanez J. M. 1995, ApJ, 448,L105Marti J. M., M¨uller E., 1996, Journ. Comp. Phys., 123, 1Marti J. M., M¨uller E., 2003, LRR, 6, 7Marti J. M., Mu¨ller E., 2015, LRCA, 1, 3Matzner C. D., 2003, MNRAS, 345, 575Marscher A. P., Gear W. K., 1985, ApJ, 298, 114Mart´ı J. M., M¨uller E., Font J. A., Ib´a˜nez J. M. Z., Marquina A.,1997, ApJ, 479, 151Mignone A., Plewa T., Bodo G., 2005, ApJS, 160, 199Mizuta A., Yamada S., Takabe H., 2004, ApJ, 606, 804Pons J. A.,Mart´ı J. M., M¨uller E., 2000, JFM, 422, 125Rezzolla L., Zanotti O., 2001, JFM, 449, 395Roe P. L., 1981, Journ. Comput. Phys., 43, 357Rushton A., Spencer R., Fender R., Pooley G., 2010, A&A, 524,A29Ryu D., Ostriker J. P., Kang H., Cen R., 1993, ApJ, 414, 1Ryu D., Brown G. L., Ostriker J. P., Loeb A., 1995, ApJ, 452, 364Ryu D., Chakrabarti S. K., Molteni D., 1997, ApJ, 474, 378Ryu D., Chattopadhyay I., Choi E., 2006, ApJS, 166, 410Sarkar S., Chattopadhyay I., 2019, IJMPD, 28, 1950037Sarkar S., Chattopadhyay I., Laurent P., 2020, A&A, 642, A209Scheck L., et. al., 2002, MNRAS, 331, 615Sod G. A., 1978, JCoPh, 27, 1Singh K., Chattopadhyay I., 2019, MNRAS, 488, 5713.MNRAS , 1–19 (2020) Joshi, Chattopadhyay, Ryu & Yadav
Synge J. L., 1957,
The Relativistic Gas . North Holland Publ. Co.,Amsterdam.Taub A. H., 1948, PhRv, 74, 328Toro, E. 1997
Riemann Solvers and Numerical Methods for FluidDynamics . Springer.van Putten M. H. P. M., 1993, ApJ, 408, L21Vyas M. K., Kumar R., Mandal S., Chattopadhyay I., 2015, MN-RAS, 453, 2992Vyas M. K., Chattopadhyay I., 2018, A& A, 614, A51Perucho M., Mart´ı J. M., 2007, MNRAS, 382, 526Walg S., Achterberg A., Markoff S., Keppens R., Meliani Z., 2013,MNRAS, 433, 1453Wen L., Panaitescu A., Laguna P., 1997, ApJ, 486, 919Wilson J. R., 1972, ApJ, 173, 431
APPENDIX A: EXACT SOLUTION OF RIEMANNPROBLEM WITH CR EOS
Riemann problem is the time evolution of an initial disconti-nuity between two non-uniform states. Depending upon theinitial conditions it can evolve into different combinations ofrarefaction and shock waves and the possible scenarios are(i) Rarefaction-Contact-Shock (RCS)(ii) Shock-Contact-Shock (SCS)(iii) Rarefaction-Contact-Rarefaction (RCR)Exact solutions of Riemann problem have been extensivelyused to test the accuracy of numerical simulation codes. Here,we present solutions for three type of Riemann problems,(i) RCS (as shown in Fig. 3b in the main text) which isknown as the shock-tube test, (ii) SCS which is mathemat-ical model of a situation when two fluids collide head onsupersonically, and (iii) wall-shock (WS) problem which isphysically equivalent to a supersonic flow which hits a wall.We also compare the analytical WS solution with that of arelativistic one-dimensional relativistic TVD code obtainedby Ryu et. al. (2006); Chattopadhyay et. al. (2013). The so-lution to the Riemann problem is obtained by solving thejump condition across shock, and solving an ordinary differ-ential equation arising out of self similarity condition of therarefaction waves.
A1 RCS
The RCS problem has a rarefaction fan, a contact disconti-nuity and a forward shock.
A1.1 Shock front evaluation
The evolution of shock wave is governed by Rankine-Hugoniotjump conditions as described in section (2.2). The post shockvelocities are obatined as v y,zb = h a γ a v y,za (cid:20) − ( v xb ) h b + ( h a γ a v ta ) (cid:21) / (A1)Where subscript a ( b ) denotes the state ahead (behind) theshock and v t is the absolute value of tangential velocity ofthe flow. v t = p ( v y ) + ( v z ) (A2) The normal component of “star state” velocity v x ∗ s prior tothe shock front is by equation (37). v x ∗ s = [ h a γ a v xa + γ s ( p ∗ − p a ) /j ][ h a γ a + ( p ∗ − p a ) { γ s v xa /j + 1 / ( ρ a γ a ) } ] (A3)Here ∗ and a are the states behind and ahead to the shockfront. For the right (left) shock, a will be the initial right(left) state. As from equation (28) j = j ( ρ a , p a , ρ ∗ , p ∗ ) (A4)we need the value of density in star state region ρ ∗ in orderto calculate v x ∗ s . This value of density can be calculated bysolving the Taub adiabat (eq. 39) for ρ ∗ using any iterativeroot finder method (we have used bisection method) for agiven value of p ∗ . Now we have to match this p ∗ and v ∗ s withthe starred pressure and velocity obtained from RF tail. A1.2 Relativistic Rarefaction Waves
Rarefaction waves are represented by the self-similar solu-tions of the flow equations. All the quantities describing thefluid depend on the variable α = ( x − x ) /t , where x is theposition of the interface separating the initial left and rightstates (see, Pons et. al. 2000, for details). Substitution of thederivatives of x and t in terms of derivative of α (equation 9)results in( v x − α ) dρdα + (cid:8) ργ v x ( v x − α ) + ρ (cid:9) dv x dα + ργ ( v x − α ) (cid:18) v y dv y dα + v z dv z dα (cid:19) = 0 (A5) ρhγ ( v x − α ) dv x dα + (1 − v x α ) dpdα = 0 (A6) ρhγ ( v x − α ) dv y dα − v y α dpdα = 0 (A7) ρhγ ( v x − α ) dv z dα − v z α dpdα = 0 (A8)From equation (A7) and (A8) we conclude that if there is notangential velocity in the initial state, no tangential flow willdevelop inside the rarefaction. Since the process along α isisentropic dpdα = hc s dρdα = ρ dhdα (A9)The determinant of the system (A5)-(A8) vanishes for thenon-trivial solution, and one obtains either α = β or α = β (see, equation 14). We are following the convention where theplus (minus) sign corresponds to the rarefaction wave prop-agating to right (left). We can reduce the system (A5)-(A8)to an ordinary differential equation (Marti & M¨uller 1994;Pons et. al. 2000) ρhγ ( v x − α ) dv x + (1 − αv x ) dp = 0 (A10)and two algebraic conditions hγv y = constant (A11) hγv z = constant (A12) MNRAS , 1–19 (2020) Equations (A11) and (A12) are similar to the equation(34) obtained for the shock. Hence the expression for thetransverse velocity given in equation (A1) can also be usedto calculate the transverse velocity components prior to therarefaction waves.Using the value of β , from equation (14), equation (A10)can be written in the form dv x dp = ± ρhγ c s p g ( α ± , v x , v t ) (A13) dv x dρ = ± c s ργ p g ( α ± , v x , v t ) (A14)Where g ( α ± , v x , v t ) = ( v t ) ( α ± − − α ± v x ) (A15)We have considered α − = β and α + = β . The sign ± cor-responds to the sign taken in equation (14). Normal velocityprior to the rarefaction can be calculated by integrating theequation (A14). For ID or fixed Γ EoS, the method to find outthe solution of the Riemann problem is relatively easy. Theequation of motion describing the RF (equation A14) admitsan analytical solution for ID EoS, and is known as the Rie-mann Invariant. And therefore, starting from the left state(region 1), we use the Riemann Invariant to obtain the flowvariables in region 3. While we use the shock conditions to ob-tain the flow variables in region 4 in terms of those in region6. Then equating velocity and pressure in region 3 and 4, weobtain a polynomial of the pressure. Solving which we recon-struct the full solution. However, for a general EoS, equation(A14) on integration does not admit an analytical solution.Therefore, the solution of Riemann problem is not trivial forgeneral EoS like CR. In the following, we present the generalmethod to solve the Riemann problem (see, Kamm 2015, forfluids with Newtonian equations of motion). A1.3 Rarefaction Fan evaluation
For right going FS and left going RF (refer to Fig. 3b), it isclear that the tail of RF (region 2) is adjacent to the starredstate (region 3) and the head adjoins the initial state (region1). From Eq. (A9) we obtain the relation between ρ, p, and c s , dpdρ = h ( ρ, p ) c s (A16)To obtain the p ∗ starting from initial state of region 1, equa-tion (A16) is numerically integrated from initial state withknown pressure and density to the final star state with given p ∗ and unknown density ρ ∗ . The problem of determining theunknown star state density is addressed computationally byfollowing the steps mentioned below,(i) We solve the differential equation by RK-4 method fora constant density step size δρ for M number of integrationsteps such that¯ ρ ∗ = ρ M = ρ M − + δρ (A17)where ¯ ρ ∗ is the provisional value of star state density. As thecorrect value of δρ is unknown, consequently we do not knowthe correct value of ρ ∗ . (ii) If I ( δρ, ρ a , p a ) be the value of pressure obtained bynumerical integration then this value of pressure should beequal to p ∗ , as the pressure remains continuous across thecontact. p ∗ − I ( δρ, ρ a , p a ) = 0 (A18)We can calculate the actual step size by solving the equa-tion (A18). Once the correct value of step size is known thenormal component of the flow velocity adjacent to the tail iscalculated by integrating the equation (A14)( v xr ) ∗ = v a ± Z ρ ∗ ,p ∗ ρ a ,p a c s ργ p g ( α ± , v x , v t ) dρ (A19)And the corresponding spatial location for the m th integra-tion step is given by x ( m ) = x + α ± ( ρ ( m ) , p ( m ) , v ( m ) ) t (A20)The final step of integration ( M th step) for equation (A20)corresponds to the position of tail. Here α + = β and α − = β (from equation 14).As the normal component of the flow velocity across thecontact discontinuity is continuous( v x ∗ r ) − ( v x ∗ s ) = 0 (A21)The value of “star state” pressure is obtained by solvingequation (A21) by any iterative root finder method. Once p ∗ is computed the tangential velocity components can becalculated using equation (35) for shock wave and equations(A11, A12) for the rarefaction wave.The solution of Riemann problem is given in table 1 andsome of these solutions are shown in fig (A1). The existenceof the upper limit of speed in relativity, causes various ve-locity components to be related to each other through theLorentz factor, solutions depend strongly on the differentcombinations of initial v t . For high values of v t , the val-ues of v x are low. Therefore, solutions with v t = 0 is fasterand hotter than solutions with v t = 0. The solutions fromrelativistic TVD codes proposed before (Ryu et. al. 2006;Chattopadhyay et. al. 2013) have been compared with theseRiemann solutions and they agree well. A2 Shock-Contact-Shock (SCS)
The collision of two streams is the physical scenario which canbe associated with the shock-contact-shock (SCS) Riemannproblem. Referring back to of Fig. 3a, SCS can be representedif regions 2 and 5 both are surfaces of the shock waves, insteadof 2 being RF. Equation (A3) provides an expression for v x in terms of v x for a left moving shock and an expression for v x in terms of v x for a right moving shock. v x = (cid:18) h γ v x + γ s ( p ∗ − p ) j (cid:19)(cid:18) h γ + ( p ∗ − p ) (cid:18) γ s v x j + 1 ρ γ (cid:19)(cid:19) − (A22) v x = (cid:18) h γ v x + γ s ( p ∗ − p ) j (cid:19)(cid:18) h γ + ( p ∗ − p ) (cid:18) γ s v x j + 1 ρ γ (cid:19)(cid:19) − (A23) MNRAS , 1–19 (2020) Joshi, Chattopadhyay, Ryu & Yadav v t v t p ∗ v x ∗ ρ ρ v t v t V s x H x t Table A1.
Solution of Shock tube problem problem at t = 0 .
25 in spatial domain x ∈ [0 ,
1] with initial data p = 10 . p = 0 . ρ = 5 . ρ = 1 . v x = 0 . v x = 0 . x = 0 . x H and x T are the positions of rarefaction head andtail respectively. In all cases we have taken composition parameter ξ = 1 . ρ v t1 =0.00v t1 =0.90v t1 =0.99 p v t1 =0.00v t1 =0.90v t1 =0.99 v x v t1 =0.00v t1 =0.90v t1 =0.99 x v t v t1 =0.00v t1 =0.90v t1 =0.99 x Γ v t1 =0.00v t1 =0.90v t1 =0.99 x h v t1 =0.00v t1 =0.90v t1 =0.99 Figure A1.
Comparison of ρ , p , v x , v t , Γ, and h , all plotted at the time t = 0 .
25. Initial conditions ρ = 5 , ρ = 1, p = 10 , p = 0 . v x = v x = 0 . v t = 0 (solid, red), v t = 0 . v t = 0 .
99 (dash-dotted, black) and for all cases initial v t = 0. The fluidcomposition is ξ = 1 . For equation (A22) j is the negative root of equation (28)and for equation (A23) it is the positive root of equation(28). Across the contact discontinuity v x − v x = 0 (A24)Equation (A24) is solved for p ∗ using the iterative root finderand rest of the quantities can be calculated once p ∗ is ob-tained. Densities in region 3 and 4 are obtained from theTaub’s adiabat for left and right shock respectively. The so-lution ( ρ, p, v x , v t , Γ , & h ) is presented in Fig. (A2a-f) fortwo time snap t = 0 .
25 (solid, dashed) and t = 0 . v x . The initial conditionused for this particular model is ρ = 10 . , p = 10 , v x = 0 . , v t = 0 .
2; and ρ = 10 . , p = 10 . , v x = − . , v t = 0 . (A25)The two shock surfaces are at the either side of x = x , andthe CD is marked by the jump in v t . A3 Solution of relativistic wall shock problem
The relativistic wall shock problem is a mathematical modelof the collision between a fluid of density ρ moving with ex-tremely high velocity v (presently along right) and a reflect-ing wall. The fluid after being reflected back, gets compressed MNRAS , 1–19 (2020) ρ (a) p (b) -0.50-0.250.000.250.50 v x (c) t = 0.25t = 0.50 x v t (d) x Γ (e) x h (f) Figure A2.
SCS:Flow variables (a) ρ , (b) p , (c) v x , (d) v t , (e) Γ and (f) h as a function of x , for Shock-Contact-Shock case at t = 0 .
25 andfor fluid composition ξ = 1 . and eventually a reverse shock is generated. The shock movesin the opposite direction of the fluid leaving behind a hot andcompressed fluid with zero velocity in post shock region. Thefluid density, pressure and velocity of fluid in post shock re-gion are represented by ρ , p and v x respectively. As velocityin post shock region v x = 0, we need only two equations tofind out ρ and p . We can use equations (37) and (39) tofind out the post shock density and pressure. From equation(37) we have h γ v x + γ s ( p − p ) j = 0 (A26)As shock is moving towards the left direction j should betaken as negative root of equation (28) j = − (cid:18) − ( p − p )( h /ρ − h /ρ ) (cid:19) / (A27)From equation (39) h − h = (cid:18) h ρ + h ρ (cid:19) ( p − p ) (A28)We can solve the equations (A26) and (A28) using theNewton-Raphson method to obtain ρ and p . The shockvelocity is calculated using (38) V s = ρ γ v − | j | p j + ρ γ (1 − v ) ρ γ + j (A29)If x is the location of wall then the shock location after time t is given by x s = x + V s t (A30)In Fig. A3a-c, we plot the flow variables ρ , p , and v x asa function of x . The solid red lines show the exact solutionsand the solutions marked by black open circles are obtainedusing relativistic TVD code of Chattopadhyay et. al. (2013).The initial conditions for this problem are ρ = 1 . , p = 1 . , and v x = 0 . ρ (a) ExactSimulation510 p (b) x v x (c) Figure A3.
WS: Solution of the wall shock problem with initialcondition ρ = 1 . p = 1 .
0, and v = 0 .
8. The solution is obtainedfor wall located at x = 1 .
0, composition parameter ξ = 0 . t = 0 . The solution contains only one discontinuity in the form ofa shock jump, across which ρ , p and v x are discontinuous. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000