Stationary relativistic jets

In this paper we describe a simple numerical approach which allows to study the structure of steady-state axisymmetric relativistic jets using one-dimensional time-dependent simulations. It is based on the fact that for narrow jets with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$v_{z}\approx c$\end{document}vz≈c the steady-state equations of relativistic magnetohydrodynamics can be accurately approximated by the one-dimensional time-dependent equations after the substitution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$z=ct$\end{document}z=ct. Since only the time-dependent codes are now publicly available this is a valuable and efficient alternative to the development of a high-specialised code for the time-independent equations. The approach is also much cheaper and more robust compared to the relaxation method. We tested this technique against numerical and analytical solutions found in literature as well as solutions we obtained using the relaxation method and found it sufficiently accurate. In the process, we discovered the reason for the failure of the self-similar analytical model of the jet reconfinement in relatively flat atmospheres and elucidated the nature of radial oscillations of steady-state jets.


Introduction
Highly collimated flows of plasma from compact objects of stellar mass, like young stars, neutron stars and black holes, as well as supermassive black holes residing in the centers of active galaxies is a wide-spread phenomenon which has been and will remain the focal point of many research programs, both observational and theoretical.Some features of these cosmic jets, like moving knots, are best described using time-dependent fluid models.However, most of these jets have sufficiently regular global structure, which is indicative of steady production and propagation and promotes development of stationary models.Such models are also easier to analyze, and they are very helpful in our attempts to figure out the key factors of the jet physics.
The simplest approach to steady-state flows is to completely ignore the variation of flow parameters across the jet.This allows to reduce the complicated system of non-linear partial differential equations (PDEs) describing the jet dynamics to a set of ordinary differential equations (ODEs) which can be integrated more easily [e.g.1,2].A similar reduction in the dimensionality is achieved in self-similar models, where unknown functions depend only on a combination of independent variables known as a self-similar variable.This also allows to reduce the original PDEs to a set of ODEs [e.g.3,4].While providing important test cases and useful insights, this approach is not sufficiently robust -boundary and other conditions that select such exceptional solutions are not always present in nature.
As it is well known to engineers working on aircraft jet engines, supersonic jets naturally develop quasi-periodic stationary chains of internal shocks, similar to what arXiv:1504.07534v1[astro-ph.HE] 28 Apr 2015 is shown in Figure 1.These shocks emerge as apart the adjustment of the jet pressure to that of the surrounding air.Interestingly, bright knots are often seen in cosmic jets and they are often interpreted as shocks [e.g.[5][6][7][8][9].Some of these knots are known to be travelling and they must be part of the jet's non-stationary dynamics.Others appear to be static and hence connected to the underlying quasi-steady-state structure of these cosmic jets.Quite often, the knots form quasi-periodic chains, reminiscent of those seen in aerodynamic jets.If the similarity is not accidental, then these knots are also related to the process of pressure adjustment.In particular, we expect the powerful cosmic jets to be expanding freely soon after leaving their central engines and to become confined by external pressure again only much later [e.g.6,10].The first shock driven into the jet by the external pressure is called the reconfinement shock.The growing observational evidence of stationary knots in cosmic jets, there has been a increase of interest to the reconfinement process among theorists in recent years [e.g 11 -17].One of the key aims of these studies was to come up with approximate analytical or semi-analytical solutions for the structure of steady-state jets.
Obviously, such shocked flows cannot be described by one-dimensional (1D) and self-similar models, which we mentioned earlier, and more complex, at least twodimensional (2D), models have to be applied instead.The system of steady-state equations of compressible fluid dynamics, not to mention magnetohydrodynamics, is already very complicated and generally requires numerical treatment.One of the ways of finding its solutions involves integration of the original time-dependent equations in anticipation that if the boundary conditions are time-independent then the time-dependent numerical solution will naturally evolve towards a steady-state [e.g.[18][19][20].One clear advantage of this approach is that it allows to use standard codes for time-dependent fluid dynamics.Such codes are now well advanced and widely available.However, this type of the relaxation approach is characterized by slow convergence and hence rather expensive.Moreover, only stable steady-state solutions can be found this way but these may not exist.
In order to speed up the convergence, one can use other relaxation methods, which are developed specifically for integrating steady-state equations [e.g.21].They often involve a relaxation variable which is called "pseudo time".However, this time evolution is not realistic but designed to drive solutions towards a steadystate in the fastest way possible.The only disadvantage of this approach is that it involves development of a specialized computer code dedicated to solving only steady-state problems.The authors are not aware of such codes for relativistic hydro-and magnetohydrodynamics.
For supersonic flows, the system of steady-state equations turns out to be hyperbolic, with one of spatial coordinates playing the role of time [22].(In the case of magnetic jets, the speed of sound is replaced with the fast magneto-sonic speed and we classify flows as sub-, trans-, or super-sonic based on its value compared to the flow speed.)In this case, one can find steady-state solutions utilizing numerical methods which were designed specifically for hyperbolic systems, like the method of characteristics or "marching" schemes.These methods have been used in the past in applications to relativistic jets [e.g. 6, [23][24][25][26] but publicly available codes do not exist yet.Their development is as time-consuming as that of time-dependent codes whereas the range of applications is much more limited.This explains their current unavailability.Moreover, when flow becomes subsonic, even very locally, this approach fails.
In this paper, we propose a new approach, which allows to find approximate numerical steady-state jet solutions rather cheaply and using widely available computer codes.To be more precise, we focus on highly relativistic narrow axisymmetric jets and show that in this regime the 2D steady-state equations of Special Relativistic MHD (SRMHD) are well approximated by 1D time-dependent equations of SRMHD.Like in the standard marching schemes, the spatial coordinate along the jet plays the role of time.This allows us to find steady-state structure of axisymmetric jets by carrying out basic 1D SRMHD simulations, which can be done with very high resolution even on a very basic personal computer.In such simulations, no special effort is needed to preserve the magnetic field divergence-free and the computational errors associated with multi-dimensionality are eliminated.As the result, more extreme conditions can be tackled.Here we focus only on relativistic jets, because of our interest to AGN and GRB jets, but we see no reason why this approach cannot be applied to non-relativistic hypersonic jets as well.Our approach is closely related to the so-called "frozen pulse" approximation, which also utilizes the similarity between the steady-state and time-dependent equations describing ultra-relativistic flows [27][28][29].In this approximation, the steady-state equations are used to analyze the dynamics of time-dependent flows.
In order to study the potential of this new approach we have carried out a number of test simulations and compared the results obtained in this way with both analytical models and numerical solutions obtained with more traditional methods.The results are very encouraging and allow us to conclude that this method is viable and can be used in a wide range of astrophysical applications.

Approximation
We start by writing down the time-dependent equations of Special Relativistic Magnetohydrodynamics (SRMHD).In this section we use units where the speed of light c = 1 and the factor 1/4π does not appear in the expression for the electromagnetic energy density.The components of vectors and tensors are given in normalized bases.The evolution equations of SRMHD include the continuity equation the Faraday equation and the energy-momentum equation where is the total stress-energy-momentum tensor, is stress-energy-momentum tensor of matter and the components of the electromagnetic stress-energy-momentum tensor are In these equations, B and E are the vectors of magnetic and electric fields respectively, p, ρ and w are the thermodynamic pressure, rest-mass density of matter and relativistic enthalpy of matter respectively, v is the velocity vector, Γ is the Lorentz factor and g is the metric tensor of space.These equations are to be supplemented with Equation of State w = w(ρ, p) and the Ohm's law of ideal MHD Finally, the magnetic field is divergence-free In this analysis, we focus on axisymmetric jets and adopt cylindrical coordinate system with the z axis coincident with the jet symmetry axis.We consider only narrow jets, so that r z 1 .
We also constrain ourselves with a relatively simple magnetic configurations where the divergence-free condition leads to In axisymmetry, the steady-state Faraday equation implies E φ = 0.When combined with Eq.9, this result yields Thus, the radial components of both the magnetic field and the velocity vectors are small compared to their axial components.We also assume that v φ 1.In fact, in the case of magnetically accelerated jets, v φ (r lc /r) when r r lc , the radius of light cylinder (See eq.66 in KVKB09).Thus, this is a good approximation for astrophysical jets.For a highly relativistic flow, the condition v z v r , v φ means Following the standard flux freezing argument, along the jet B φ /B z (r j /r lc ) −1 , where and r j is the jet radius.(This argument does not apply to turbulent jets, which are non-axisymmetric and allow non-trivial conversion of components.)Hence one may argue that far away from the central engine In order to introduce the key idea of our approach we consider first the steadystate continuity equation: Using the condition (12) we may replace v z with unity.This makes Eq.14 identical to the 1D time-dependent version of the continuity equation.In order to stress this point we replace z with t and write: Similarly, all 2D steady-state equations can be approximated by the corresponding 1D time-dependent equations.Let us show this for the equations of magnetic field.The 1D version of the divergence free condition reads ∂ r (rB r ) = 0 or rB r = const.
Thus if B r vanishes outside of the jet, which is expected when it is in direct contact with ISM, then one has to put B r = 0 everywhere.As we shell see, the terms involving B r are sub-dominant in all other equations and hence this is a reasonable simplification.Moreover, once the 1D solution is found, one can substitute the determined B z (r, z) into the 2D divergence free condition and solve it for B r (r, z).The result can then be used to verify that B r (r, z) B z (r, z).The φ component of the Faraday equation can be written as where i = r, z.In steady-state, the first term vanishes, the next two terms are of the order B z v φ /z and small compared to the last two terms, which are of the order B φ v z /z.Removing these small terms we obtain the approximate steady-state equation Finally, we replace v z with unity, z with t, and obtain This is indeed the 1D version of the φ component of eq.16.Now consider the z component of the Faraday equation, Since v z is very much constant, ∂ i v z v z /x i .Due to this, the second and third terms are small compared to the last two terms.Removing them we obtain the approximate steady-state equation Now once again we replace v z with unity and z with t to obtain which is the 1D version of the z component of eq.16.
Finally, we analyze the energy-momentum equations.These can be written as so the steady-state versions are These already have the same form as the 1D time-dependent equations, so we only need to show that Let us start with the hydrodynamic contribution.First, we notice that T tz hd = wΓ 2 v z wΓ 2 as v z 1 .
Thus, T zt hd T tt hd .Then we notice that hd .Now we inspect the electromagnetic contributions.First, we find good estimates for the components of electric field.From eq.9 it follows that and In fact, it is easy to show that Indeed, for magnetically accelerated jets B φ ΩrB z (e.g.KVKB09) for r r lc .Hence Using these estimates we find that and hence T zt em T tt em .Moreover, and hence T zz em T tz em as well.Next we show that T zφ em T tφ em .Indeed, Finally, we show that T zr em T tr em .First, we find straight away that Since E z E r E z B φ , we only need to show that B z B r is significantly smaller compared to these terms.This is indeed the case as B z B r v r B 2 z whereas using eqs.25 and 26 we obtain z .Thus, within our approximation the steady-state 2D equation of energymomentum reduces to which is the 1D time-dependent energy-momentum equation.
Given that in relativistic fluid dynamics small differences between the magnitudes of energy and momentum may result in huge variations of Lorentz factor and even lead to inconsistency, one could feel uneasy about the approximations we make.However, the final result is exactly the system of 1D time-dependent SRMHD and this means that self-consistency is not compromised.For example, the flow speed will not exceed the speed of light because of the errors of our approximation.

Numerical Implementation
The analysis of Section 2 shows that as long as they are applied to narrow jets with high Lorentz factor, the axisymmetric steady-state equations of SRMHD are very close to 1D time-dependent equations of SRMHD in cylindrical geometry.This suggests that it may be possible to use time-dependent simulations with 1D SRMHD codes to study the 2D structure of steady-state jet solutions.However in order to be able to do this, we also need to find a way of accommodating the 2D boundary conditions of steady-state problems in such simulations.
For 2D supersonic flows we need to fix all flow parameters at the jet inlet and impose some conditions at the jet boundary, consistent with it being a stationary contact wave.No boundary conditions are needed for the outlet boundary -its flow parameters are part of the solution.In the corresponding 1D problem, the 2D boundary conditions at the inlet boundary simply become the initial conditions of the 1D Cauchy problem.The final 1D solution corresponds to the slice of the 2D solution at the outlet boundary.As to the contact discontinuity at the 2D jet boundary, the situation is not that trivial.
Suppose that the total pressure at this boundary is a function of z, p = p b (z).When we replace z with t this becomes p = p b (t).Thus we need somehow to impose time-dependent boundary conditions.In the simulations presented below, the following approach was utilized: 1) we extend the computational domain so that it includes the external gas, 2) we track the point separating the jet from the external gas and 3) we reset the external gas parameters according to the prescribed functions of time every computational time step.
In practice, the jet radius can be located using the passive scalar τ , which satisfies the advection equation The jet material is initialized with τ = 1 and we set τ = 0 in the external gas.However, even after one time-step few cells at the separation between the two regions will have 0 < τ < 1 due to numerical diffusion.Thus, in order to separate the jet from the external gas one better use a value of the scalar in this range.We find that τ t = 10 −2 serves well as such a threshold.During the reset, we restore the initial values for the scalar with the separation at the new jet boundary in order to limit the effect of numerical diffusion.After the reset, the 1D jet boundary is no longer a contact but a more general discontinuity.In particular, the jet plasma will generally have radial velocity component.If it is positive, but in the external gas it is set to zero, then a shock wave will launched into the jet when this discontinuity is resolved.If it is negative, then this will be a rarefaction wave.On the one hand, this reflects how the information about changing environment is communicated to the interior of a steady-state jet.On the other hand, in 1D simulations the strength of the emitted wave depends on the external density -higher density, and hence lower temperature, will result in stronger waves moving into the jet.This is obviously not so for 2D steady-state jets, which react only to the external pressure.Thus additional measures need to be undertaken.First, in order to negate the effect of the radial velocity jump at the jet boundary, the radial velocity of the external gas is reset not to zero but to its value at the last jet cell.Second, in order reduce the role of the external gas inertia, it helps to set its density to a low value, so that its sound speed becomes relativistic.Although we have not tried this, one could set the polytropic index of the external gas to Γ = 2, which would make the sound speed of ultra-relativistically hot gas equal to the speed of light.

Bowman's jet
To test the validity of our approach, we first use our method to reproduce the numerical steady-state solutions for supersonic unmagnetized jets obtained by Bowman [25, B94] using the marching scheme described in [24].In this study pressurematched uniform jets with zero opening angle are injected into an atmosphere with the pressure distribution with z s = 10, z c = 50.According to this equation, the external pressure initially decreases almost as fast as ∝ z −2 but at z > z c becomes uniform.The initial jet radius r 0 = 1 and the injection nozzle is located at z = z s .The equation of state is that of Synge [30] for an electron-proton plasma.The initial jet pressure p j = p 0 .
For the comparison we selected the model with the Mach number M j = 15 and the initial temperature T j = √ 10 × 10 13 K.At such a high temperature the EOS of electro-proton plasma is almost the same as that of the pure proton gas.The latter was used in our simulations.
Bowman's solution is shown in the top part of figure 1.As the external pressure decreases rapidly, the jet quickly becomes under-expanded and enters the phase of almost free expansion.When it enters the outer region of constant pressure it becomes over-expanded and a reconfinement shock is pushed towards its axis, where it gets reflected.Gas passed though these two shocks becomes hot and its pressure rises.As a result, the jet becomes somewhat under-expanded again and begins to expand for the second time.Then it becomes over-expanded again and another shock is pushed into the jet and so on.
In the bottom part of this figure, we show the results of our 1D simulations for this jet using exactly the same visualization technique as in the original paper.The agreement between the two solutions is quite remarkable.A very good match for the maximal radial extension and the oscillation-length of the jet is obtained.The successive reconfinement shocks are somewhat sharper than in B94, most likely due to the application of a shock-capturing scheme.We checked our approach against other numerical models of B94 as well.In all models, the results for profile of jet radius and Mach number are in good agreement.Noticeable but still minor differences arise only for the colder models, most likely due to the different equation of state used in our simulations.

Self-similar models of jet reconfinement
The problem of reconfinement of initially free-expanding steady-state jets is quite important and a number of authors have tried to find simple analytic of semianalytic solutions.Falle [31] and Komissarov & Falle [10] used the Kompaneets approximation, which assumes that the gas pressure immediately downstream of The Kompaneets approximation is accurate only for very narrow jets.To improve on it, one also has to take into account the variation of the gas pressure across the shocked layer [11].In our second test, we compare our results with the semianalytical model by [15, thereafter KBB12], who assumed self-similarity of the flow in this layer.This assumption is more suitable for the case where the reconfinement shock never reaches the jet axis, because otherwise the distance where this occurs sets a characteristic length scale.
KBB12 studied jets with ultra-relativistic equation of state (w = 4p, γ = 4/3 ), propagating in a power-law atmosphere, These jets emerge from a nozzle at z = z 0 with the Lorentz factor Γ 0 , opening angle θ 0 = 1/Γ 0 and pressure p 0 .The initial velocity distribution correspond to a conical flow originating from z = 0 and hence the initial jet radius r 0 = z 0 tan(θ 0 ).They could only find self-consistent solutions for 8/3 ≤ κ < 4 and later argued that for κ < 8/3 the entropy of the shocked layer must increase with the distance along the jet in order for the solution to be consistent with the energy conservation [17].They proposed that this additional heating is caused by multiple shocks driven into the flow as it gradually collimates.
We selected the KBB12 model with κ = 8/3 and Γ 0 = 50 and made simulations on a uniform grid with only 300 cells (each run took only several CPU minutes on a laptop using only one core of its processor).Our results are shown in the first panel of figure 2, which should be compared with figure 7 in KBB12.Again we find a very good agreement between the models -at z = 9 we have got the jet radius r j ≈ 0.114 and the shock radius r s ≈ 0.07, whereas in KBB12 r j = 0.110 and r s = 0.064.
In order to understand the difference between the cases with κ > 8/3 and κ < 8/3, we also computed models with κ = 7/3 and 2 -the evolution of the Lorentz factor in these models is shown in the second and third panels of figure 2 respectively.In these plots we see no evidence of the additional shocks proposed in [17].Neither could we find them in plots of other parameters.However, figure 2 suggests that in the models with κ = 7/3 and 2 the reconfinement shock is much stronger than in the model with κ = 8/3.Moreover, the shock strength is increasing with the distance along the jet.As the result, the entropy of the shocked layer in the models with κ = 7/3 and 2 is higher and its mean value across the layer is growing with the distance.This is why the self-similar model of KBB12, which assumed isentropy of the flow in the layer, failed for κ < 8/3.In contrast, in the model with κ = 8/3 the mean entropy of the layer remains fairly constant (see fig. 3).Based on these results, we conclude that the value of κ = 8/3 is not special, but the accuracy of the constant-entropy approximation used in KBB12 greatly reduces as κ decreases.

Magnetized jets. 1D versus 2D solutions
The steady-state structure of magnetized jets is more complex, mainly due to the non-trivial contribution of the magnetic tension to the force balance.A number of authors have tackled this problem analytically using various approximations [e.g.16,[32][33][34].However, none of these studies deliver a model suitable for detailed testing of our numerical approach.Dubal & Pantano [35] studied the steady-state structure of relativistic jets with azimuthal magnetic field using the method of characteristics.This would be a good test case but the setup of their simulations is ambiguous.We have tried several variants of the setup but each time failed to reproduce the results.The mechanisms of magnetic collimation and acceleration of relativistic jets were studied numerically by [36,37] and [38] using a "rigid wall" outer boundary.While this allows for a well-controlled experiment, Komissarov et al. [19] have shown that the connection between the shape of the boundary and the external pressure gradient is not straightforward, with significant degeneracy.For this reason, we concluded that in the magnetic case the best way of testing the performance of our 1D approach would be via new 2D axisymmetric time-dependent simulations using the relativistic AMRVAC code [39].
The problem we selected for this test is similar in its setup to the one described in Sec.4.2 as it also involves a jet propagating through the atmosphere with the powerlaw pressure distribution (30), and the nozzle is still located at z = z 0 .However, this time the jet is magnetized and the rest mass density of its particles is not negligibly small.The jet structure at the inlet is that of a cylindrical jet in magnetostatic equilibrium, which satisfies the following force balance equation where b φ = B φ /Γ is the azimuthal component of the magnetic field as measured in the fluid frame using normalized basis and p t is the sum of the gas pressure and the magnetic pressure due to the axial magnetic field B z [40].Equation (31) has infinitely many solutions -given a particular distribution for b φ (r) one can solve this equation for the corresponding distribution of the pressure p t (r).We adopted the "core-envelope" solution of Komissarov [40]: where r j is the jet radius and r m is the radius of its core (Note a typo in the expression for α in [40].).As one can see, the core is pinched and in the envelope the magnetic field is force-free.This may be combined with any distribution of density and axial velocity.We imposed ρ = ρ 0 and with ν = 8; this gives an almost "top-hat" velocity profile.The velocity vector is set to be aligned with the jet axis, so v r = v φ = 0.This solution is illustrated in Figure 4. We considered two models, A and B. In the models A, the magnetic field is purely azimuthal and the other parameters are r j = 1, r m = 0.37, b m = 1, ρ 0 = 1, z 0 = 1, β m = 0.34, Γ 0 = 10.The local magnetization parameter in this model does not exceed σ max = 0.7 and thus the jet is only moderately magnetized.The jet core is relativistically hot, with the gas pressure reaching p max = ρ at the axis, which opens the possibility of efficient hydrodynamic acceleration once the jet is allowed to expand.In the simulations we used the adiabatic equation of state w = ρ+(γ/γ−1)p with γ = 4/3.
In model B, this configuration is modified to include nonvanishing longitudinal magnetic field B z .In particular, we considered the case where the gas pressure p = αp 0 everywhere within the jet and which keeps p t unchanged.In this model, the magnetic field is force-free not only in the envelope but also in the core.The other parameters of this model that differ from those of model A are ρ 0 = 0.05 and β m = 0.14.The corresponding magnetization is much higher, with σ max = 17.Model B turned out too stiff for our 2D code, but presented no problems in 1D simulations.For this reason we compare here the 1D and 2D results for model A only.In these simulations we used the atmosphere with κ = 1.The computational domain is 20 r j in the radial direction and 800 r j in the axial direction.
First, let us describe the overall jet structure found in these simulations.Initially, as the jet enters the region of rapidly declining external pressure, it expands rapidly and a rarefaction wave moves towards its axis.Eventually, the jet becomes overexpanded, its expansion slows down, and a reconfinement shock sets in.It reaches the axis at z ≈ 400, gets reflected and then returns to the jet boundary at z ≈ 700 (see Figure 5).
To quantify the convergence of the 1D simulations we carried out simulations with three different resolution and used this data to determine the grid-convergence index where f 1 , f 2 are solutions with doubled and quadrupled resolution compared to f 0 and |f a − f b | 1 is the difference between two solutions in the L 1 norm.We found that η ≈ 1, as this is expected for a TVD scheme in the presence of discontinuities.At 6400 grid cells in the radial direction, the density contours become visually unchanged on the scale of figure 5.The 1D solution with this resolution was used for comparison with the results of our 2D simulations.In what follows we refer to it as the "converged" 1D solution.
The initial solution in our 2D simulations was constructed via interpolation of the converged 1D solution onto the 2D cylindrical grid.Since we did not include gravity to balance the pressure gradient in the external atmosphere, in order to preserve the atmosphere in its initial state the atmospheric parameters were reset to their initial values every time step, just like this was done in the 1D case.In order to test the convergence of 2D solutions, we made three runs with doubled resolution, N r = 400, 800, and 1600 cells in the radial direction.The number of cell in the axial direction was always twice the number of cells in the radial one.
Typically, the 2D solutions exhibited some evolution at first but then quickly settled into a stationary state.For example in the case of N r = 400, the timestep-totimestep relative variation of the conserved flow variables dropped below 6 × 10 −6 at t = 1000 and remained approximately constant thereafter.Furthermore, the relative L 1 error of density between times t = 1000 and t = 3000 was 2.8 × 10 −4 , indicating that a stationary state had been reached.The 2D solutions converge with the grid-convergence index η > 1.25 over the entire simulated time.
The difference between the converged 1D solution and the relaxed 2D solutions with with N r = 800 (dotted lines) and N r = 1600 (dashed lines) is illustrated in Figure 5 which shows the mass density distribution.One can see that the 2D solutions are very close to the 1D solution and that the difference decreases with the resolution of 2D runs.To quantify the difference between the relaxed 2D solutions and the converged 1D approximate solution we introduce the parameter For the 2D solution with N r = 400 cells in the radial direction we obtain δρ 6%, for N r = 800 δρ 4.3% and for N r = 1600 the relative error decreased to δρ 3.2%.This shows that the approximation error of our 1D approach is at the level of no more than 3%.Converged 1D solution for a stationary magnetized jet (solid lines) and two corresponding 2D solutions found via the relaxation method, one with the 1600 × 3200-resolution (dashed lines) and one with the 800 × 1600-resolution (dotted lines).The lines are 10 rest-frame density contours consecutively spaced by the factor of one half from the starting value of ρmax = 1.The solution involves a reconfinement shock which reaches the jet axis at z ≈ 400.

Magnetized jets in power-law atmospheres
Komissarov et al. [19] derived an approximate equation for the radius of highly magnetized jets, in the limit where it strongly exceeds that of the light cylinder.Using this equation they concluded that in the case of power-law atmosphere with 0 < κ < 2 the jet radius increases as Lyubarsky [33] developed the theory of Poynting-dominated jets further and using more accurate analysis found that the expansion is modulated by oscillations with the wavelength growing as These oscillations can be understood as a standing magneto-sonic wave bouncing across the jet.Indeed, denote the wave speed as a m .Then the jet crossing time is τ c = r j /a m in the co-moving jet frame and t c = Γτ c in the rest frame of the atmosphere.As the wave is advected along the jet almost at the speed of light the wavelength of the associated structure is Since for the jets considered in [19,33] a m c and Γ ∝ r j we obtain λ ∝ r 2 j and using Eq.39 recover Eq.40.The results (39,40) are well suited for testing of our approach.To this aim, we carried out additional 1D simulations with models A and B described in Sec.4.3.Since in model A the jet is not Poynting-dominated, it allows us to explore the regime not covered in [33].To see how sensitive these results may be to the assumptions made in [19,33] let us consider unmagnetized relativistic jets.From the mass conservation law we obtain r j ∝ (Γρ) −2 .For relativistically cold jets with p ρc 2 we have Γ const and thus whereas for the hot jets Γ ∝ r j and thus where we put γ = 4/3.The last result is the same as for the Poynting-dominated jets.Even for the cold jets the difference is rather minor, e.g. for γ = 5/3 the index in Eq.42 differs from κ/4 only by κ/20 and for γ = 4/3 by κ/8.
In order find λ we note that for cold jets a 2 m ∝ (p/ρ) ∝ z −κ(1−1/γ) and hence Eq.41 yields Eq.40 independently of the value of γ.For hot jets, a m const and Eq.41 still leads to Eq.40 if we use γ = 4/3.Thus, the law (40) for the wavelength of oscillations is very robust.
Figure 6 illustrates the overall jet structure in model A and its response to changes in the parameter κ of the external atmosphere.One can see that this weakly magnetized jet also shows a combination of secular expansion and oscillations.These oscillations appear to be a generic feature of the adjustment process of supersonic jets to variations of external pressure, which occurs by means of magneto-sonic waves travelling across the jet.In the very beginning, the decrease of external pressure makes the jet under-expanded and a rarefaction wave is launched from the jet boundary towards the jet axis.Behind this wave the radial velocity is positive and the flow expands.The rarefaction reduces the jet pressure and at some point it becomes over-expanded.Now a compression wave is driven inside the jet.Behind this wave the jet expansion slows down and eventually turns into a contraction.The contraction increases the jet pressure and at some point it becomes under-expanded again and then the whole cycle repeats.
The deviation from the force-balance corresponding to the secular jet expansion is due to the finite propagation speed of the waves -as they move across the jet they are also advected downstream by the supersonic flow.As the result, the jet interior reacts to the changes in the external pressure with a delay.It keeps expanding when the internal pressure is already too low and keeps contracting when it is already too As κ increases, the wavelength of the oscillation increases as well.This is expected as the more rapid overall expansion of the jet in an atmosphere with larger κ means that it takes longer for a magneto-sonic wave to traverse the jet, not only as the result of the larger jet radius but also as the result of its higher Mach number (and hence smaller Mach angle).Overall, this is very similar to the well-known evolution of under-expanded supersonic jets studied in laboratories.Normally, their compressive transverse waves steepen into shocks.In our model A with κ = 1 we also detect shocks, but they become progressively weaker, suggesting that they may disappear further out along the jet.For κ = 0.5, shocks do not form at all.The exact reason for this in not yet clear.
Figure 7 shows the evolution of other flow parameters in model A with κ = 1.Both the secular and oscillatory behavior of the jet radius are mirrored in the variation of the Lorentz factor.The secular expansion leads to secular increase of the Lorentz factor as both the thermal and the magnetic energy are being converted into the kinetic energy of the flow.The thermal acceleration is most pronounced in the jet core, which is relativistically hot at the inlet.The conversion of magnetic energy is supported by the overall decrease of the magnetization parameter σ along the jet (the right panel of figure 7).The oscillations of the jet radius lead to additional increase of the Lorentz factor during the expansion phase and its decrease upon contraction.
Figure 7 shows the same parameters for the highly magnetized jet of model B with with κ = 1.In this model, the reconfinement shock is no longer present.This may be related to the fact that in this model the fast magneto-sonic speed is higher and the corresponding jet Mach number is lower, at the inlet M lower wavelength of the jet oscillations as it takes less time for the waves to traverse the jet.The theoretical predictions for the secular evolution of the jet radius and the wavelength of its oscillations are put to a quantitative test in figure 9, which shows the jet radius rescaled according to its expected secular evolution against z 1−κ/2 .In such plots, the mean jet radius and the wavelength of oscillations should remain constant.For the highly magnetized jet of model B the scaling factor is z κ/4 and for the low magnetized jet of model A it is z 3κ/8 , as appropriate for a cold hydrodynamic jet with γ = 4/3.In general, we obtain a very good agreement with the theoretical scalings for the mean jet radius, both in the low-and high-magnetization limit.A small departure from the z κ/4 -scaling is observed for case B with κ = 1 -it expands slightly faster.This could be because the jet magnetization is not sufficiently high and decreases more rapidly with distance than in the atmosphere with κ = 0.5.The evolution of the wavelength scaling is also in a very good agreement with the theory -the residual error is between 0.7% and 3.4%.

Conclusions
In this paper we presented a novel numerical approach, which can be used to determine the structure of steady-state relativistic jets.It is based on the similar- ity between the two-dimensional steady-state equations and the one-dimensional time-dependent equations of SRMHD with the cylindrical symmetry in problems involving narrow highly-relativistic (v z ≈ c) flows.Such similarity has already been utilized in the so-called "frozen pulse" approximation where dynamics of timedependent relativistic flows is analyzed using the steady-state equations [27][28][29].
Here we do the opposite construct approximate steady-state solutions via numerical integration of the time-dependent equations.The main advantage of this approach is utilitarian.First, it allows us to use computer codes for relativistic MHD (or hydrodynamics in the case of unmagnetized flows), which are now widely available, in place of highly-specialized codes for integrating steady-state equations, which are not openly available at the moment.Moreover, the reduced dimensionality means that the computational facilities can be very modest -a basic laptop will suffice.In contrast, the relaxation method based on integration of two-dimensional time-dependent equations can be computationally quite expensive.
We compared numerical solutions obtained with this approach with analytical models and numerical solutions obtained with other techniques.The considered problems involved a variety of flows both magnetized and unmagnetized, with different equations of state and external conditions.The results show that the method is sufficiently accurate and robust.In both models the expected average expansion is captured quite well.To show that the oscillation wavelength scales as λ ∝ z κ/2 , the x-axis has been rescaled accordingly.In order to visually separate the curves corresponding to different values of κ, they have been shifted up by a factor of κ.
Although we focused only on relativistic flows, we see no reason why this approach cannot be applied to non-relativistic hypersonic flows.For such flows, the axial velocity of bulk motion plays the role of the speed of light in the substitution z = ct used in our derivations.
As a byproduct of our test simulations, we obtained two results of astrophysical interest.We demonstrated that the failure of the self-similar model of the jet reconfinement in power-law atmospheres with the index κ < 8/3 [15] is rooted in the assumption of isentropy of the shocked layer, which is made in this model.In reality, the reconfinement shock becomes stronger with the distance along the jet, resulting in a strong spatial variation of the entropy.We also found that the radial oscillations of steady-state jets, discovered in the analytical models of Poynting-dominated jets [33] is a generic part of the jet adjustment to the space-variable external pressure and not specific to the high-magnetization regime only.The oscillations are standing waves induced by the variation.
The steady-state solutions are useful for elucidating some key factors in flow dynamics and may closely describe some of the observed phenomena in astrophysical jets.However, they are often subject to various instabilities which may dramatically modify the flow properties.Most instability studies, both analytical and numerical, deal with very simple problems where the steady-state solution is readily available.In more realistic setup, the issue of finding the steady-state solution, which can then be subjected to perturbations, becomes more involved and this is where our method can be applied in the instability studies.

Figure 1
Figure 1 Reconfinement of the M j = 15, T j = √ 10 × 10 13 K jet.The top panel is a reproduction of Figure 3 from B94.The bottom panel shows the solution obtained with our method.In each panel, the top halves show 50 pressure contours (spaced by the factor of 1.18) and the bottom halves show the temperature parameter τ ≡ ρh/(ρh − p) in 50 contours (spaced by the factor of 1.003).The light grey lines are streamlines.

Figure 4
Figure 4 Initial radial structure of the magnetized jets in the test simulations described in Sec.4.3.

Figure 5
Figure5Converged 1D solution for a stationary magnetized jet (solid lines) and two corresponding 2D solutions found via the relaxation method, one with the 1600 × 3200-resolution (dashed lines) and one with the 800 × 1600-resolution (dotted lines).The lines are 10 rest-frame density contours consecutively spaced by the factor of one half from the starting value of ρmax = 1.The solution involves a reconfinement shock which reaches the jet axis at z ≈ 400.

Figure 6
Figure 6 Structure of steady-state magnetized jets obtained via time-dependent 1D simulations.The plots show the co-moving density distribution for model A with κ = 0.5 and κ = 1.The distance along the vertical axis is defined as z = ct/r j , where r j is the initial jet radius.The white contour shows the jet boundary, using the passive scalar.

Figure 7
Figure7shows the same parameters for the highly magnetized jet of model B with with κ = 1.In this model, the reconfinement shock is no longer present.This may be related to the fact that in this model the fast magneto-sonic speed is higher and the corresponding jet Mach number is lower, at the inlet M 3 compared to M 10 in model A. The lower Mach number is also responsible for the observed

Figure 8
Figure8As figure7for model B with κ = 1.As the Poynting-flux vanishes on the axis (and the thermal energy is negligible), we obtain a hollow jet with fastest region away from the axis.Due to the increased fast-magneto-sonic speed (thus lower Mach-number) compared to the case of figure7, no reconfinement-shock forms and the jet-oscillation frequency is increased.

4 Figure 9
Figure9Compensated jet-expansion laws for models (top) and B (bottom).In both models the expected average expansion is captured quite well.To show that the oscillation wavelength scales as λ ∝ z κ/2 , the x-axis has been rescaled accordingly.In order to visually separate the curves corresponding to different values of κ, they have been shifted up by a factor of κ.