© The Authors 2023

1. Introduction

Astrophysical jets are commonly observed in most, if not all, types of accreting sources. They are emitted from young stellar objects (YSOs; Bally et al. 2007; Ray et al. 2007; Ray & Ferreira 2021), active galactic nuclei (AGNs) and quasars (Boccardi et al. 2017), close interacting binary systems (Fender & Gallo 2014; Tudor et al. 2017), and even post-AGB stars (Bollen et al. 2017). Despite the different central objects (be it a black hole, a protostar, a white dwarf, or a neutron star), these jets share several properties: (i) they are supersonic collimated outflows with small opening angles, (ii) the asymptotic speeds scale with the escape speed from the potential well of the central object, and (iii) they carry away a sizeable fraction of the power released in the accretion disk. As the only common feature shared by all these different astrophysical objects is the existence of an accretion disk, it is natural to seek a jet model that is related to the disk and not to the central engine. This universal approach is further consistent with the accretion–ejection correlations observed in these objects (see e.g., Merloni & Fabian 2003; Corbel et al. 2003; Gallo et al. 2004; Coriat et al. 2011; Ferreira et al. 2006; Cabrit 2007 and references therein).

Despite these general common trends, astrophysical jets do show some differences in their collimation properties. For instance, the core-brightened extragalactic jets, classified as FRI jets after Fanaroff & Riley (1974), appear conical and show large-scale wiggles (see e.g., Laing & Bridle 2013, 2014 and references therein). On the contrary, the edge-brightened FRII jets appear nearly cylindrical, with a terminal hotspot (Laing et al. 1994; Boccardi et al. 2017). Most of the jets imaged with very long baseline interferometry do not appear as continuous flows, but can be modeled as a sum of discrete features, known as blobs or knots, usually associated with shocks (Zensus 1997). Those shocks are assumed to originate either from pressure mismatches at the jet boundary with the external medium or from major changes at the base of the flow (e.g., new plasma ejections or directional changes), with some of these knots being stationary features (e.g., Lister et al. 2009, 2013; Walker et al. 2018; Doi et al. 2018; Park et al. 2019).

On the other hand, jets from young forming stars do not seem to have such a clear FRI/FRII dichotomy and often display evidence of a conflictual interaction (shocks) with the ambient cloud medium (Reipurth & Bally 2001). This might be consistent with the suspicion that the FR dichotomy would only be a consequence of the jet interaction with its environment, with low-power jets remaining undisrupted and forming hotspots in lower mass hosts (Mingo et al. 2019). However, protostellar jets might also have intrinsic collimation properties different from those of extragalactic jets, possibly because they are nonrelativistic outflows. This is an open question.

Since the seminal model of Blandford & Payne (1982; hereafter BP82), it is known that a large-scale vertical magnetic field threading an accretion disk is capable of accelerating the loaded disk material up to super-fast magnetosonic speeds. This acceleration, usually termed magneto-centrifugal, goes along with an asymptotic collimation of the ejected plasma thanks to the magnetic tension associated with the toroidal magnetic field (hoop stress). In this semi-analytical model, a self-similar ansatz has been used allowing the full set of stationary ideal magnetohydrodynamic (MHD) equations to be solved. Later, this self-similar jet model was generalized in different ways by altering the magnetic field distribution (Contopoulos & Lovelace 1994; Ostriker 1997), thermal effects (Vlahakis et al. 2000; Ceccobello et al. 2018) and was even extended to the relativistic regime (Li et al. 1992; Vlahakis & Königl 2003; Polko et al. 2010, 2014). However, it is unclear whether or not self-similarity affects the overall jet collimation properties. Not only are both the axis and the jet-ambient medium region not taken into account, but the final outcome of the jets (i.e., acceleration efficiency, jet kinematics and opening angle, presence of radial oscillations, or even shocks) may well also be impacted by the imposed geometry.

Using the only class of self-similar jet models smoothly connected to a quasi-Keplerian accretion disk, Ferreira (1997; hereafter F97) showed that these super-fast magnetosonic jets systematically undergo a refocusing toward the axis (see also Polko et al. 2010). Such a recollimation is due to the dominant effect of the internal hoop stress and has nothing to do with a pressure mismatch at the jet–ambient medium interface proposed to explain knotty features in extragalactic jets (Komissarov & Falle 1998; Perucho & Martí 2007; Perucho 2020). According to F97, recollimation would be generic to MHD jets anchored over a large range of Keplerian accretion disks. This is indeed verified for warm outflows (Casse & Ferreira 2000a) and weak magnetic fields (Jacquemin-Ide et al. 2019).

While MHD recollimation is also seen in nonself-similar works (e.g., Pelletier & Pudritz 1992), other strong assumptions are usually made, leaving the question of the jet asymptotics open. Heyvaerts & Norman (1989) used another approach based on the electric poloidal current (or Poynting flux) still present at infinity. These authors showed that any stationary axisymmetric magnetized jet will collimate at large distances from the source to paraboloids or cylinders, depending on whether or not the asymptotic electric current vanishes. This important theorem was later generalized (Heyvaerts & Norman 2003a) by taking into account the issue of current closure and its effect on the geometry of the solution (Okamoto 2001, 2003). However, the theorem only addresses the asymptotic electric current, and it is unclear how much of this current is actually left as no simple connection with the source can be made.

Connecting the asymptotic electric current to the source is naturally done with time-dependent MHD simulations. Those reaching the largest spatial scales treat the accretion disk as a boundary condition, allowing the jet dynamics to be studied independently of the disk (Ustyugova et al. 1995, 1999; Ouyed & Pudritz 1997a,b, 1999; Krasnopolsky et al. 1999, 2003; Ouyed et al. 2003; Anderson et al. 2005, 2006; Fendt 2006; Pudritz et al. 2006; Porth & Fendt 2010; Porth & Komissarov 2015; Staff et al. 2010, 2015; Stute et al. 2014; Barniol Duran et al. 2017; Tesileanu et al. 2014; Tchekhovskoy & Bromberg 2016; Ramsey & Clarke 2019). The drawback of these platform jet simulations is their huge degree of freedom, which is attributable to the fact that several distributions must be specified at the lower injection boundary. It has therefore been very difficult to determine the exact generic results on jet collimation that can be derived from them.

In summary, despite many theoretical and numerical studies, no connection has been firmly established between the jet-launching conditions and the jet-collimation properties at observable scales. This work is the first of a series designed to bridge this gap. Our approach here is to assess whether the general results obtained within the self-similar framework still hold in full 2D time-dependent simulations. We will address in particular whether the existence of recollimation shocks is indeed unavoidable for the physical conditions expected in Keplerian accretion disks, as proposed by F97.

As a consequence of this approach, we focus only on steady-state jets, allowing us to directly confront our simulations with MHD jet theory. It is clear that most if not all astrophysical jets exhibit time-dependant features; see for example Cheung et al. (2007) for M87 or Bally et al. (2007) for young stars. However, our goal is not to reproduce a specific astrophysical jet, but instead to deduce the generic behaviors of MHD jets emitted from Keplerian accretion disks.

The paper is organized as follows. Section 2 describes our numerical setup and boundary conditions, which mimic an axial spine (related to the central object) surrounded by a self-similar cold jet. As analytical studies require huge spatial and temporal scales, a special temporal numerical scheme has been designed. Our reference simulation, which corresponds to a typical BP82 jet, is described in length in Sect. 3. We show that recollimation shocks are indeed obtained in agreement with the analytical theory. This is the first time that such shocks are obtained self-consistently, showing that these are not artificial biases due to the mathematical ansatz used, but consequences of the jet-launching conditions. A parametric study is presented in Sect. 4, where we vary the magnetic flux exponent α and the jet mass load κ, confirming the striking qualitative correspondance between our numerical simulations and analytical solutions. In particular, we show that the asymptotic jet collimation depends mostly on the exponent α. However, the existence of an axial spine introduces quantitative differences hinting at a possible role of the central object in affecting the collimation properties of the jets emitted by the surrounding disk. Our results are finally confronted to the wealth of previous 2D numerical simulations in Sect. 5 and we conclude in Sect. 6.

2. MHD simulations of jets from Keplerian disks

2.1. Physical framework and governing equations

We intend to study the collimation properties of magnetically driven jets emitted from Keplerian accretion disks, as depicted in Fig. 1. The disk is settled from an inner radius R_d to an outer radius R_ext = 5650.4R_d and is assumed to be orbiting around a central object of mass M located at the center of our coordinate system. The disk itself is not computed and we assume that it behaves like a JED, with consistent prescribed boundary conditions. As we use a spherical grid, the central object as well as its interaction with the disk are assumed to occur inside a sphere of radius R_d (the green zone in Fig. 1). This sphere defines the inner boundary discussed below.

Fig. 1. Sketch of our computational domain. The central object and its interaction with the innermost disk are located below the inner boundary at Rd (green region), the near-Keplerian jet-emitting disk (JED) being established from Rd to the end of the domain Rext. An axial outflow (the spine) is emitted from the central regions (in red) and the jet is emitted from the JED (in blue). The solid purple line represents a recollimation shock surface starting on the axis at a height Zc. For each point N lying on this surface, we use local poloidal unit vectors (e⊥, e∥), respectively perpendicular and parallel to the shock surface. Also, at any point M inside the domain, we either use spherical (eR,eθ,eϕ) or cylindrical (er,eϕ,ez) coordinates.

We further assume that a large-scale magnetic field is threading both the disk and the central object. The existence of this field allows the production of two outflows, one from the disk (blue region in Fig. 1) and one from the central spherical region (red region in Fig. 1). Hereafter, we always refer to the disk-emitted outflow as the “jet”, and to the outflow emitted from the spherical region (green zone in Fig. 1) as the “spine”. As our goal is to focus on the dynamics of the jet itself, we try to limit the influence of the spine as much as possible.

Two systems of coordinates centered on the mass M are used, spherical (R, θ, ϕ) and cylindrical (r, ϕ, z). Both the spine and the jet are assumed to be in ideal MHD and we numerically solve the usual set of MHD equations. This includes mass conservation

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{u})=0,\end{array}$$(1)

where ρ is the density and u the flow velocity, and the momentum equation is

$$\begin{array}{c}\hfill \frac{\partial \rho \mathit{u}}{\partial t}+\mathrm{\nabla }·[\rho \mathit{u}\mathit{u}+(P+\frac{\mathit{B}·\mathit{B}}{2{\mu }_o})-\frac{\mathit{B}\mathit{B}}{{\mu }_o}]=-\rho \mathrm{\nabla }{\mathrm{\Phi }}_G,\end{array}$$(2)

where P is the thermal pressure, B the magnetic field, and Φ_G = −GM/R the gravitational potential due to the central mass.

The evolution of the magnetic field is determined by the induction equation,

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}+\mathrm{\nabla }\times (\mathit{B}\times \mathit{u})=0.\end{array}$$(3)

As we focus on highly supersonic flows, we decided to derive the pressure P and internal energy by solving the entropy equation,

$$\begin{array}{c}\hfill \frac{\partial \rho S}{\partial t}+\mathrm{\nabla }·\left(\rho S\mathit{u}\right)=0,\end{array}$$(4)

where S = P/ρ^Γ is the specific entropy and Γ = 1.25 is the polytropic index (the same for all our simulations). This simple advection equation guarantees that the pressure does not assume nonphysical (e.g., negative) values. But on the other hand, it does fail to provide the correct entropy jump in the shocks. However, as long as the thermal energy of the flow remains negligible compared to the kinetic and magnetic energy, this should not present a problem.

Thanks to axisymmetry, the poloidal magnetic field can be computed using the magnetic flux function Ψ (which corresponds to R sin θ A_ϕ, where A is the vector potential),

$$\begin{array}{c}\hfill B_R=\frac{1}{R^{2}sin\theta }\frac{\partial \mathrm{\Psi }}{\partial \theta }{B}_{\theta }=-\frac{1}{Rsin\theta }\frac{\partial \mathrm{\Psi }}{\partial R},\end{array}$$(5)

which already verifies ∇ ⋅ B = 0. An axisymmetric magnetized jet can therefore be seen as a bunch of poloidal magnetic surfaces defined by Ψ(R, θ) = constant, nested around each other and anchored on the disk for the jet and in the central object for the spine.

In steady state, Eqs. (1) to (4) lead to the existence of the following five MHD invariants, namely quantities that remain constant along each magnetic surface (Weber & Davis 1967):

• the mass flux to magnetic flux ratio η(Ψ) = μ_oρu_p/B_p,

• the rotation rate of the magnetic surface Ω_*(Ψ) = Ω − ηB_ϕ/(μ₀ρr),

• the total specific angular momentum carried away by that surface L(Ψ) = Ωr² − rB_ϕ/η,

• the Bernoulli invariant $E(\mathrm{\Psi })=\frac{u^2}{2}+H+{\mathrm{\Phi }}_G-{\mathrm{\Omega }}_{\ast }r{B}_{\varphi }/\eta$,

• the specific entropy S(Ψ) = P/ρ^Γ,

where Ω = u_ϕ/r and $H=\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }-1}\frac{P}{\rho }$ is the specific enthalpy. We make use of these relations when designing boundary conditions.

2.2. Numerical setup

We solve the above set of equations using the MHD code PLUTO¹ (Mignone et al. 2007). We configured PLUTO to use a second-order linear spatial reconstruction with a monotonized-centered limiter on all the variables. This method provides the steeper linear reconstruction compatible with the stability requirements of the scheme. A flatter and more diffusive linear reconstruction is employed in a few cells around the rotation axis to dampen numerical spurious effects that typically appear in these zones due to the discretization of the equations around the geometrical singularity of the axis. The HLLD Riemann solver of Miyoshi & Kusano (2005) is employed to compute the intercell fluxes. This solver is one of the best suited to properly capture Alfvén waves, a crucial element in properly modeling trans-Alfvénic flows. So as to match the order of the spatial reconstruction, we chose a second-order Runge-Kutta scheme to advance the equations in time. The ∇ ⋅ B = 0 condition is ensured by employing a constrained transport (CT) scheme, enforcing that constraint at machine accuracy.

The two-dimensional computational domain is discretized using spherical coordinates (R, θ) assuming axisymmetry around the rotation axis of the disk. The domain encompasses a spherical sector going from the polar axis (θ = 0) to the surface of the disk that is assumed to be θ = π/2 for simplicity, and is resolved with N_θ = 266 points in the θ direction. The cell size in the θ direction is mostly uniform, but decreases on a few cells near the axis. This is essential to our setup, as the expected collimation shocks are formed near the axis: an overly low resolution in this zone would prevent their formation. In the radial direction, the grid goes from R_d to R_ext = 5650.4R_d with N_R = 1408 points in a logarithmic spacing (ΔR ∝ R) so as to make sure the cells remain approximately square (ΔR ≈ RΔθ) far from the axis.

We choose such a huge numerical domain because our goal is to capture the recollimation shocks predicted in the self-similar solutions of F97. According to their Fig. 6, those shocks may occur at altitudes spanning from several hundred to a thousand times the jet launching radius. Our spherical grid with a ratio of 5650.4 therefore provides a suitable range to observe such shocks. The drawback of these huge spatial scales is of course the terrible contrast in timescales. The Keplerian time scaling in r^3/2 means that in order to compute a full orbit at the outer disk edge, the inner one should have completed over 4 × 10⁵ orbits. This would be barely affordable if we were to use a standard evolutionary scheme. In order to achieve such long timescales, we designed a specific method that accelerates the numerical integration using larger and larger time steps to evolve the equations as the solution starts to converge towards a steady-state. This method is very successful and allowed us to significantly boost the evolution of our jets (see Appendix A).

2.3. Initial conditions

Our initial magnetic field is assumed to be potential, which leads to a second-order partial differential equation on Ψ(R, θ). In order to represent suitable self-similar solutions, we solve this equation by assuming

$$\begin{array}{c}\hfill \mathrm{\Psi }={\mathrm{\Psi }}_d{(R/R_d)}^{\alpha }\mathrm{\Phi }(\theta ),\end{array}$$(6)

where the function Φ(θ) has been determined assuming that the initial field is potential, that is, current-free and force-free (J_ϕ = 0).

The exponent α is a free parameter of the model leading to B_R ∝ B_θ ∝ R^α − 2. For α = 0, field lines are conical, for α = 1 they are parabolic, and α = 2 describes a constant (straight) vertical field. The seminal BP82 solution is for α = 3/4.

As the magnetic field is potential, no magnetic force is initially imposed on the plasma. It is therefore assumed to be in spherically symmetric hydrostatic equilibrium (u = 0) with dP/dR = −ρGM/R². We choose the following trivial solution:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \rho & ={\rho }_a{\left(\frac{R}{R_d}\right)}^{2\alpha -3}\hfill \\ \hfill P& =\frac{1}{4-2\alpha }\frac{{\rho }_{a}GM}{R_d}{\left(\frac{R}{R_d}\right)}^{2\alpha -4}.\hfill \end{array}\end{array}$$(7)

The sound speed C_s is defined as $C_s^2=\partial P/\partial \rho =\mathrm{\Gamma }P/\rho$. In the following, ρ_a refers to the density at the axis immediately above the central sphere.

2.4. Boundary conditions

Boundary conditions must be imposed at the polar axis (θ = 0, R from R_d to R_ext), at the outer frontier (R = R_ext, θ from 0 to π/2), at the JED surface (θ = π/2, R from R_d to R_ext), and at the spine boundary (R = R_d, θ from 0 to π/2). On the polar axis, usual proper reflecting boundary conditions are imposed on all quantities. The special treatment done for the other three boundaries is described, especially for the JED and spine boundaries where mass is being injected.

2.4.1. Outer boundary (R = R_ext)

“Outflow” conditions are imposed at the outer frontier: for ρ, P, B_R, B_θ, RB_ϕ, u_R, u_θ and u_ϕ, the gradient along the radial direction is conserved, and we use the Van Leer slope limiter to avoid spurious oscillations. Additionally, we enforce a positive toroidal Lorentz force on the subalfvénic part of this boundary.

2.4.2. Jet generation: the jet-emitting disk (θ = π/2)

We need to specify eight quantities (ρ, u, B, P) that must be representative of the fields expected at the surface of a JED. As the lifted material gets accelerated along a field line, its poloidal velocity will become larger than the slow magnetosonic V_sm, poloidal Alfvén V_Ap, and fast magnetosonic V_fm phase speeds. Crossing each of these critical speeds defines a regularity condition that determines one quantity at the jet basis, therefore leaving five free functions to be specified. However, we wish to control the mass loss from the JED, which requires that the injected outflow be already super-slow magnetosonic (hereafter super-SM). We therefore have to impose six functions at the JED boundary, leaving two free to adjust over time, B_ϕ and B_R, the latter controlling the magnetic field bending.

Our choice of boundary conditions at θ = π/2 (so that R = r) is therefore as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \rho & ={\rho }_d{\left(\frac{R}{R_d}\right)}^{2\alpha -3}\hfill \\ \hfill P& ={\rho }_d\frac{C_{s_d}^2}{\mathrm{\Gamma }}{\left(\frac{R}{R_d}\right)}^{2\alpha -4}\hfill \\ \hfill B_{\theta }& =-B_d{\left(\frac{R}{R_d}\right)}^{\alpha -2}\hfill \\ \hfill u_{\theta }& =-u_d{\left(\frac{R}{R_d}\right)}^{-1/2}\hfill \\ \hfill u_R& =u_{\theta }\frac{B_R}{B_{\theta }}\hfill \\ \hfill u_{\varphi }& ={\mathrm{\Omega }}_{\ast }r+u_{\theta }\frac{B_{\varphi }}{B_{\theta }},\hfill \end{array}\end{array}$$(8)

where Ω_* is the angular velocity of the magnetic surfaces (an MHD invariant in steady-state). We assume ${\mathrm{\Omega }}_{\ast }={\mathrm{\Omega }}_K=\sqrt{GM/r^3}$, in agreement with a near Keplerian accretion disk, leaving four normalizing quantities, ρ_d, C_sd, B_d, and u_d, to be specified at R_d. These distributions are consistent with a self-similar JED and describe an ideal steady MHD flow with u_p ∥ B_p², anchored on a disk that imposes magnetic field lines rotating at the Keplerian angular velocity Ω_K. We note that the fixed component of the magnetic field threading the disk (B_θ) is actually the initial condition to conserve the magnetic flux injected into the computational domain, and only B_R is allowed to vary in response to the jet dynamics.

In order to pick up values at R_d that are consistent with the jet calculations performed by BP82 or F97, we express the JED boundary conditions as a function of four dimensionless parameters: (1) the jet density ρ_d is fixed with respect to the density at the polar axis using ρ_d = δρ_a; (2) the disk sound speed (temperature) is defined relative to the Keplerian speed with ϵ = C_sd/V_Kd; (3) the magnetic field strength B_d is controlled by measuring the θ component of the poloidal Aflvén speed with respect to the Keplerian speed, namely $\mu =V_{\mathit{Ad}}/V_{\mathit{Kd}}=B_d\sqrt{R_d/({\mu }_o{\rho }_{d}GM)}$; and (4) the (vertical) injection velocity u_d can be determined with the well-known mass-loading parameter κ introduced by BP82 using

$$\begin{array}{c}\hfill \kappa =\frac{{\mu }_o{\rho }_d{u}_d{V}_{\mathit{Kd}}}{B_d^2}=\frac{u_d{V}_{\mathit{Kd}}}{V_{\mathit{Ad}}^2}=\frac{u_d}{V_{\mathit{Kd}}}\frac{1}{{\mu }^2}.\end{array}$$(9)

By fixing u_d/V_Kd = 0.1 for all the simulations, we obtain κ = 0.1/μ². In order to be able to fix the value of the injection speed u_d and therefore the JED mass flux, we must require that u_p > V_sm. As we are mostly interested in producing cold MHD outflows, we assume ϵ = 0.01 so that the θ component of the sonic Mach number is M_sθ = u_d/C_s = u_d/V_Kdϵ = 10. As the total poloidal speed at the jet boundary is larger than u_d, the sonic Mach number M_s = u_p/C_s > 10. As the poloidal Alfvén speed at the disk surface is much larger than the sound speed, C_s > V_sm and M_s > 1 is enough to warrant a super-SM condition.

We decided to vary the mass load and the disk Alfvén speed by only changing the disk density ρ_d (and keeping the injection speed u_d and the disk magnetic field B_d constant for all the simulations). As a consequence, the density contrast δ can be expressed as a function of μ (or κ). We assume the relation δ = 100/μ² = 1000κ. We highlight the fact that with our parametrization the JED boundary conditions are determined by only one dimensionless parameter, typically κ, while the other two free parameters μ and δ are determined as a function of κ, and ϵ is fixed for all the simulations.

2.4.3. Spine generation: the central object (R = R_d)

In the spine, we follow a similar methodology to that in the JED and specify six quantities along the inner spherical boundary at R_d. This again leaves two quantities that are free to evolve, B_ϕ and B_θ. In order to conserve the magnetic flux injected into the computational domain, we fix B_R(θ) to its initial value. We note that, as the B_θ(R) profile is fixed along the JED boundary (θ = π/2) and B_R(θ) is kept constant in time along the spine boundary (R = R_d), the total poloidal field and its inclination B_R/B_θ do not change with time at the inner radius of the disk (R = R_d, θ = π/2). The strength of the magnetic field is already determined by the value of μ chosen in the JED. As the outflowing material leaving the central region is in ideal MHD and we are looking for a steady jet, one has u_θ = u_RB_θ/B_R. This leaves us with the four distributions ρ, C_s, u_R, and u_ϕ to be specified along θ.

If the central object possesses its own magnetosphere, then R_d might be considered as the disk truncation radius. What would be encapsulated within R_d could then be a complex combination of a stellar wind plus any type of magnetospheric wind (steady or not; see for instance Zanni & Ferreira 2013 and references therein). If the central object is instead a black hole, then R_d might be considered as the innermost stable circular orbit and what is hidden inside R_d would highly depend on the black hole spin. While a nonrotating black hole would provide no outflow, a rather strong magnetic flux concentration is seen to occur in GRMHD simulations of spinning black holes, leading to the generation of powerful outflows through the Blandford-Znajek process (see e.g., Blandford & Znajek 1977; Tchekhovskoy et al. 2010; Liska et al. 2018 and references therein).

However, our goal is to study the outcome of the jet emitted from the disk. We therefore decided to minimize the influence of the spine as much as possible. This was found to be an almost impossible task; details are given below. As pointed out in early works on magnetized rotating objects (e.g., Ferreira & Pelletier 1995), the jet power depends on the available electromotive force (emf) e = ∫E_m ⋅ dl = ∫(u × B_p)⋅dl. While the disk provides an emf e_disk ≃ ∫Ω_KrB_zdr, the central region provides e_obj ≃ ∫ΩrB_RR_ddθ. An obvious way to decrease e_obj is therefore to allow Ω to decrease as one goes from the disk to the pole. We therefore use (in agreement with steady-state ideal MHD) u_ϕ = Ω_*r  +  u_RB_ϕ/B_R, with magnetic surfaces rotating as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\ast }={\mathrm{\Omega }}_a(1-f(\theta ))+{\mathrm{\Omega }}_{\mathit{Kd}}f(\theta ),\end{array}$$(10)

where f(θ) is a spline function varying smoothly from zero at θ = 0 to unity θ = π/2 (see Appendix B). Most of the simulations presented in this paper were done with Ω_a = 0 (but not all, see Sect. 4.3). This choice is consistent with a nonrotating black hole but also with an innermost disk radius (our R_d) well below the co-rotation radius in the case of a star.

The fixed radial speed is defined through the sonic Mach number M_sR, by u_R = M_sRC_s. For M_sR we assume a constant value along θ that can be derived from the JED boundary conditions by assuming its continuity at the inner disk radius R_d, M_sR = M_sθ|B_R/B_θ|_d = 10|B_R/B_θ|_d > 1. As the field inclination at the inner disk radius |B_R/B_θ|_d is constant, also M_sR does not change with time. The sound speed at the base of the spine is computed as

$$\begin{array}{c}\hfill C_s=C_{\mathit{sa}}(1-f(\theta ))+C_{\mathit{sd}}f(\theta ),\end{array}$$(11)

where the sound speed on the axis C_sa is computed so as to verify the Bernoulli integral E_a = E(θ = 0) at the axis. As the MHD contribution vanishes on the axis, one directly obtains

$$\begin{array}{c}\hfill C_{\mathit{sa}}^2=\frac{\mathit{GM}}{R_d}\frac{1+e_a}{\frac{1}{2}{M}_{\mathit{sR}}^2+\frac{1}{\mathrm{\Gamma }-1}}\end{array}$$(12)

and

$$\begin{array}{c}\hfill u_{\mathit{Ra}}^2=\frac{\mathit{GM}}{R_d}\frac{1+e_a}{\frac{1}{2}+\frac{1}{M_{\mathit{sR}}^2(\mathrm{\Gamma }-1)}},\end{array}$$(13)

where u_Ra is the injection radial speed on the axis and e_a = E_aR_d/GM is the Bernoulli integral normalized to the gravitational energy at R_d and will be used as a parameter to fix the axial spine temperature. We note that the normalized Bernoulli integral for the jet at R_d writes e_d = λ_d − 3/2 + ϵ²/(Γ − 1), where $\lambda =L(\mathrm{\Psi })/\sqrt{GM{r}_o}$ is the magnetic lever arm parameter, measured here at the anchoring radius r_o = R_d. As our jets are cold, enthalpy plays no role and e_d is mostly determined by λ (which is known only once the simulation has converged to a steady state). For our simulations, we expect a λ of around 10 (see our parameter space Fig. 15). We therefore fix e_a = 2 in order to obtain a spine with a smaller energetic content than the surrounding jet. We note that, with our choice of parameters, the injection speed along the axis, set in Eq. (13), is higher than the escape speed. As the flow is cold and there is no magneto-centrifugal acceleration along the symmetry axis, the flow will gradually slow down along R in the spine from this very high speed in its core. The spine flow can cross the Alfvén and fast-magnetosonic critical points due to a decrease of the magnetic field intensity, not thanks to a flow acceleration.

Finally, for the density, we need to smoothly connect its axial value ρ_a to the much larger value injected at the disk surface ρ_d. We choose to do this by computing ρ(θ) = ηB_R/(μ_ou_R), with the MHD invariant η following

$$\begin{array}{c}\hfill \eta ={\eta }_a(1-f(\theta ))+{\eta }_{d}f(\theta ),\end{array}$$(14)

with η_a and η_d being fully determined (see Appendix B). This method ensures that the mass flux to magnetic flux ratio has a smooth variation from the disk to the axis. For numerical stability reasons, as the strongest magnetic field is on the axis, the density in the code is normalized to ρ_a, providing a dimensionless density at the axis of 1.

2.5. Summary of parameters and normalization

Each simulation is entirely determined by the following dimensionless parameters and quantities: the radial exponent of the magnetic field α; the Bernoulli parameter e_a (equal to 2 for most cases) and the spine angular frequency on the axis Ω_a (equal to 0 in most simulations); and the cold jet parameters: κ, $\mu =\sqrt{0.1/\kappa }$, δ = 1000κ, ϵ = 0.01.

With our choices, we ensure that the injected flow is everywhere super-SM and that the main emf is due the JED, which is magnetically launching a cold jet. The profiles of several quantities along the magnetic flux near the lower boundary (inner spherical boundary and disk) reached by our reference simulation K2 at the final time can be seen in Fig. B.1.

This leaves us with only two free parameters, α and κ. We do not explore their whole range here but keep them within the parameter space of jets from JEDs as obtained by F97 but also by the solutions of Contopoulos & Lovelace (1994) and the simulations of Ouyed & Pudritz (1997a). The radial exponent α is varied from 10/16 to 15/16. In a strict self-similarity, this exponent must be consistent with the underlying disk, namely α = (12 + 8ξ)/16, where ξ is the disk ejection efficiency defined with the disk accretion rate as Ṁ_a(r) ∝ r^ξ (Ferreira & Pelletier 1995). However, our simulations are not strictly self-similar because of the presence of the axis and its spine, and so we also explore a slightly smaller α than the fiducial BP82 value α = 12/16. As discussed further below, values of α ≥ 1 are numerically problematic. The mass load parameter κ is varied between 0.05 and 1, which is a range globally consistent with BP82 and F97 jets, both solutions leading to a flow recollimation toward the axis.

The MHD equations have been solved with PLUTO and the results will be presented in dimensionless units. Unless otherwise specified, lengths are given in units of R_d, velocities in units of $V_{\mathit{Kd}}=\sqrt{GM/R_d}$, time in units of T_d = R_d/V_Kd, densities in units of ρ_a, magnetic fields in units of $B_d=V_{\mathit{Kd}}\sqrt{{\mu }_o{\rho }_a}$, mass fluxes in units of ${\dot{M}}_d={\rho }_a{R}_d^2{V}_{\mathit{Kd}}$ and powers in units of $P_d={\rho }_a{R}_d^2{V}_{Kd}^3$. In order to be more specific, we translate these quantities for the case of a young star, assuming a star of one solar mass with an innermost disk radius R_d = 0.1 au, namely

$$\begin{array}{cc}\hfill V_{\mathit{Kd}}& =94.3{\left(\frac{M}{M_{\odot }}\right)}^{1/2}{\left(\frac{R_d}{0.1\times \mathrm{au}}\right)}^{-1/2}{\mathrm{km}\mathrm{s}}^{-1}\hfill \\ \hfill {\dot{M}}_d& =3.3\times {10}^{-10}\left(\frac{{\rho }_a}{{10}^{-15}{\mathrm{g}\mathrm{cm}}^{-3}}\right){\left(\frac{M}{M_{\odot }}\right)}^{1/2}{\left(\frac{R_d}{0.1\mathrm{au}}\right)}^{3/2}{M}_{\odot }{\mathrm{yr}}^{-1}\hfill \\ \hfill P_{\mathrm{jet}}& =6.7\times {10}^{41}\left(\frac{{\rho }_a}{{10}^{-15}{\mathrm{g}\mathrm{cm}}^{-3}}\right){\left(\frac{M}{M_{\odot }}\right)}^{3/2}{\left(\frac{R_d}{0.1\mathrm{au}}\right)}^{1/2}W\hfill \\ \hfill B_d& =10.6{\left(\frac{{\rho }_a}{{10}^{-15}{\mathrm{g}\mathrm{cm}}^{-3}}\right)}^{1/2}{\left(\frac{M}{M_{\odot }}\right)}^{1/2}{\left(\frac{R_d}{0.1\mathrm{au}}\right)}^{-1/2}\mathrm{G}\hfill \\ \hfill T_d& =1.8{\left(\frac{M}{M_{\odot }}\right)}^{-1/2}{\left(\frac{R_d}{0.1\mathrm{au}}\right)}^{3/2}\mathrm{days}.\hfill \end{array}$$(15)

The list of all the simulations performed in this paper is provided in Table 1, with their input parameters α and κ and several quantities that are measured at the final stage t_end of the simulation. As explained in Sect. 2.4.1, the values of μ and δ are dictated by the values of κ. As discussed below, all our simulations display several recollimation shocks. In the table, we provide only the altitude (measured at the axis) of the first main recollimation shock Z_shock. As stationary jets require them to become super-Afvénic and super-fast magnetosonic (hereafter super-A and super-FM, respectively), we also display the colatitudes ${\theta }_A^{\text{ext}}$ and ${\theta }_{\mathrm{FM}}^{\mathrm{ext}}$ of the intersection of the outer boundary R_ext and the Alfvén and FM surfaces, respectively. The last super-FM magnetic surface (defining the jet) can then be followed down to the disk, allowing us to identify the largest anchoring radius r_o,FM that we consider in the JED. This allows us to measure the mass flux emitted from the JED as ${\dot{M}}_{\mathrm{jet}}={\int }_{R_d}^{r_{o,\mathrm{FM}}}\rho \mathit{u}·\mathbf{d}\mathit{S}$ and compare it with the mass loss emitted from the spine only ${\dot{M}}_{\mathrm{spine}}={\int }_0^{\pi /2}\rho \mathit{u}·\mathbf{d}\mathit{S}$. We also compute the power emanating from the jet $P_{\text{jet}}={\int }_{R_d}^{r_{o,\text{FM}}}\rho E\mathit{u}\cdot \mathbf{d}\mathit{S}$ and compare it to the power emanating from the spine $P_{\text{spine}}={\int }_0^{\pi /2}\rho E\mathit{u}\cdot \mathbf{d}\mathit{S}$.

Simulations K1 to A5 were performed with a nonrotating spine, namely Ω_a = 0 and e_a = 2. Our reference simulation K2 is extensively analyzed in the following section. This reference was repeated with a lower resolution in K2l –all other things being equal– to verify numerical convergence. Section 4 addresses the influence of κ (simulations K1 to K5) and α (simulations A1 to A5). In Sect. 4.3, the effect of a rotating spine (simulation SP) is briefly addressed.

3. The Blandford & Payne case

3.1. Overview

In this section, we discuss our reference simulation K2 performed for the BP82 α = 3/4 magnetic field distribution and a mass-loading parameter κ = 0.1. It was run up to t_end = 6.5 10⁵ T_d and has reached a steady-state in a sizable fraction of our computational domain (a quarter of an orbit has been done at r_o = R_ext).

Figure 2 displays the final stage reached by K2 at t_end. The black solid lines are the poloidal field lines, the dotted red line is the Alfvén surface (where the Alfvénic Mach number m = u_p/V_Ap is equal to unity) and the dashed red line is the FM surface (where the FM Mach number n = u_p/V_fm = 1). The left panel shows our simulation on the full computational domain, with a close-up view on the scale used by (Fendt 2006) in the right panel. The background color is the logarithm of n on the left and the logarithm of the density on the right. The last magnetic surface characterizing the super-FM jet is anchored at r_o, FM = 323 in the JED, and the critical surfaces (A and FM) both achieve a conical shape over a sizable fraction of the domain, which is characteristic of a self-similar steady-state situation. Our spine also achieves super-FM speed at an altitude of z ∼ 260.

Fig. 2. Snapshot at tend of our Blandford & Payne simulation K2. Left: global view with field lines anchored on the disk at ro = 3; 15; 40; 80; 160; 320; 600; 1000; 1500, where the background is the logarithm of the FM mach number n. Right: close-up view of the innermost regions, where the background is the logarithm of the density. In both panels, black solid lines are the poloidal magnetic surfaces, the yellow solid lines are isocontours of the poloidal electrical current, and the red dashed (resp. dotted) line is the FM (resp. Alfvén) critical surface.

The poloidal velocity vectors can be seen in Fig. 3. The velocity decreases radially very rapidly, mirroring the injection conditions, going from 3.5 at the spine to 0.2 at the edge of the super-FM zone (in V_Kd units). The white lines show streamlines inside which 50%, 75%, and 100% (from left to right) of the total super-FM mass outflow (spine + jet) rate is being carried in. These lines are anchored in the disk at r₀ = 10, r_o = 66, and r_o, FM = 323, respectively. As dṀ/dR = 2πRV_z falls off very rapidly, this plot shows that even ejection from a very large radial domain may be observationally dominated by the innermost, highly collimated regions up to r_o ∼ 10, with the outer “wide angle wind” probably remaining barely detectable.

Fig. 3. Snapshot of our reference simulation K2 at tend. We use the same color coding as in Fig. 2, left. The black arrows show the poloidal velocity. The white lines are streamlines inside which (from left to right) 50%, 75%, and 100% of the super-FM (spine+jet) mass outflow rate is carried in.

The yellow solid lines are isocontours of the poloidal electric current I = 2πrB_ϕ/μ_o. These contours are very useful as they allow us to grasp several important features of the simulation: (1) The typical butterfly shape of the initial accelerating closed electric circuit can be seen up to a spherical radius R ∼ 3000. (2) For disk radii r_o ≳ 2000, the electric current flowing out of the disk reaches the outer boundary (most of it in the sub-A regime at high colatitudes) and re-enters in the jet at smaller colatitudes, in the super-FM regime. (3) More importantly, several current sheets can be clearly seen (as an accumulation of current lines), highlighting the existence of several standing (stationary) recollimation shocks. To our knowledge, this is the first time that simulations of super-FM jets exhibit the patterns predicted in analytical jet studies. Justifications for this assessment are provided in Sect. 5.2.

These shocks are best seen in Fig. 4, which presents a zoom onto the region of interest. Five shocks (highlighted in colors) can be seen starting near the polar axis, following approximately the expected shape of the MHD characteristics in self-similar jets (see Figs. 3 in Vlahakis et al. 2000; Ferreira & Casse 2004). They are located at Z₁ = 1850, Z₂ = 2000, Z₃ = 2160, Z₄ = 2372, and Z₅ = 2634. Only two of these shocks span a significant lateral portion of the jet (those best seen also in the left panel in Fig. 2). The first one (in red) leaves the axis at an altitude Z₁ = 1850 (labeled Z_shock in Table 1) and stays within our domain, finally merging with the FM surface (red dashed curve) around (r = 2500, z = 3800). The second shock starts at Z₅ and leaves the simulation domain at (r = 1800, z = 5200), and is therefore not fully captured by our simulation. For this reason, only the first shock is extensively described here.

Fig. 4. Close-up view of our reference simulation K2 at tend showing the shock forming region, with field lines anchored at ro = 1.2; 2; 3; 4; 5; 7; 9. The five shocks are highlighted in red, orange, cyan, blue, and purple. We use the same color coding as in the left panel of Fig. 2.

It can be seen that all shocks do occur only after the magnetic surface has started to bend toward the axis (with a decreasing cylindrical component B_r), and give rise to a sudden outward refraction of the surface with its ouflowing material. This, plus the fact that their positions remain steady in time (see below), justifies our use of the name “standing recollimation shocks” to describe them. While these standing recollimation shocks appear quite generic for our set of simulations, we stress that fantastically large spatial and temporal scales are required to see them.

This simulation was also performed with a two-times-smaller resolution (see K2l in Table 1). We also observed standing recollimation shocks that are similar, although with less complexity.

3.2. Quasi steady-state jet and spine

Jet production is a very rapid process that scales with the local Keplerian timescale (especially As μ is constant with the radius). It is therefore an inside-out build-up of the jet with its associated electric circuit, until the innermost jet regions (including the spine) reach the outer boundary. This will take a time of typically t_ext(r_o = 1)∼R_ext/V_z, with the maximal jet speed $V_z\simeq \sqrt{2\lambda -3}{V}_K$. According to BP82, cold outflows with moderate inclinations and κ = 0.1 reach λ ∼ 10 (see their Fig. 2), which is exactly what we also obtain (see Fig. 15, with a minimum value λ = 11). This leads to a dynamical time t_ext ∼ 1.3 × 10³ (in T_d units) and so the innermost jet regions achieve an asymptotic state very early. However, the transverse (radial) equilibrium of the outflowing plasma is still slowly adjusting, because, as time increases, more and more of the outer magnetic surfaces achieve their asymptotic state, providing an outer pressure and modifying the global jet transverse equilibrium. This, in turn, necessarily modifies the shape of the magnetic surfaces as well as the associated poloidal electric circuit.

This means that the MHD invariants along each magnetic surface can always be defined (each surface is quasi-steady), but that they are also slowly evolving in time as the global magnetic structure evolves. For our simulation K2, the last super-FM magnetic surface is anchored at r_o, FM = 323, defining a local Keplerian time T_K = 5.8 × 10³. At that distance, the boundary is located at Z_ext ∼ R_ext cos θ_FM ∼ 4430, leading to a dynamical time t_ext(r_o, FM)∼Z_ext/V_z(r_o, FM)∼1.6 × 10⁴ (with λ = 14) as the speed distribution on the disk is Keplerian. We can therefore expect MHD invariants within our jet to evolve much less only after a time ∼10⁴.

This slight evolution of the jet quantities over time is illustrated in Fig. 5. We chose to look at global quantities, such as the radius r_o, FM of the last super-FM surface, the jet mass-loss rate Ṁ_jet, and the two colatitudes ${\theta }_A^{\text{ext}}$ and ${\theta }_{\mathrm{FM}}^{\mathrm{ext}}$ that define the position of the two critical surfaces. We pick up their values at t = 5.1 × 10⁵ and plot their evolution by normalizing them to this “initial” value. It can be seen from Fig. 5 that their evolution is quite obvious: r_o, FM keeps on increasing, leading to an increase in Ṁ_jet and a decrease in both ${\theta }_A^{\text{ext}}$ and ${\theta }_{\mathrm{FM}}^{\mathrm{ext}}$. But the relative variations are less than 3% for r_o, FM and 1% for the other quantities.

Fig. 5. Late evolution of several global jet quantities for the simulation K2: the radius ro, FM of the last super-FM surface, the jet mass-loss rate Ṁjet, and the two colatitudes ${\theta }_A^{\text{ext}}$ and ${\theta }_{\mathrm{FM}}^{\mathrm{ext}}$ that define the position of the two critical surfaces. A slight drift from their initial value is indeed observed. The values provided in Table 1 are those achieved at the final time.

We therefore consider that our simulation K2 has achieved a relatively global steady-state. In physical units (and for our choice of axial density ρ_a), the jet mass loss is about 2 × 10⁻⁷ M_⊙ yr⁻¹ with a magnetic field around 10 G at 0.1 au. The spine mass loss is only ∼10% of the jet mass loss, and so we can safely presume that the dynamics are mostly controlled by the JED, as expected. However, as the spine power is comparable to the jet power (P_spine/P_jet = 0.81), the impact of the spine on the collimation and topology of the electric field cannot be neglected. The influence of the spine is detailed in Sect. 4.3.

Figure 6 shows the various contributions to the Bernoulli integral E(ψ) along a magnetic surface anchored at r_o = 100 at the final stage t_end. It can be seen that E is indeed conserved and that jet acceleration follows the classical pattern (Casse & Ferreira 2000a): the kinetic energy (green) increases thanks to the magnetic acceleration, leading to a decrease in the magnetic contribution (magenta). Enthalpy (red) is negligible in this cold outflow. The presence of the shock is clearly seen around Z = 3800: the flow is suddenly slowed down and the energy is transferred back to the magnetic field, in agreement with the Rankine-Hugoniot jump conditions (see Appendix C). Beyond the shock, MHD acceleration is resumed but, at the edge of our domain, the magnetic field still maintains around 45% of the initial available energy.

Fig. 6. Evolution of the various energy contributions along a magnetic surface of anchoring radius r0 = 100 for the simulation K2 at tend: the Bernoulli invariant E, the gravitational potential ΦG, the total specific kinetic energy u2/2, the enthalpy H, and the magnetic energy −Ω*rBΦ/η. The absicssa is the altitude Z(Ψ).

The evolution along a magnetic surface of the five MHD invariants is shown in Fig. 7 for two surfaces, one anchored at r_o = 100 (left) and the other at r_o = 1000 (right). In order to plot η, Ω_*, L, E, and S on the same figure, we normalize each quantity by its initial (at the disk surface) value at t_end. On the left, variations of the invariants can indeed be seen but only at the shock located at Z = 3800; they remain very small, much less than 1% for all but the entropy (which is conserved to machine accuracy). On the right plot, the field line is anchored beyond r_o, FM and the flow remains sub-FM while crossing no shock. Variations of the invariants are again observed, but always less than 0.3%. This shows that the PLUTO code is quite efficient and the MHD solution is indeed steady.

Fig. 7. Evolution of the MHD invariants along field lines of two different anchoring radii ro at tend for simulation K2. All invariants have been normalized to their values at ro. The absicssa is the altitude Z(Ψ).

3.3. Jet collimation

Before analyzing the shocks in the following section, let us briefly discuss jet collimation. Figure 8 shows the initial magnetic field configuration (in dotted lines) along with the final configuration (solid lines) obtained at t_end. Each color is associated with a different anchoring radius r_o, allowing us to see the evolution from the initial potential field and the final full MHD solution. This plot clearly illustrates how magnetic collimation works. As the poloidal electric circuit responsible for the collimation must be closed, its sense of circulation must change within the whole outflow (defined as both the jet and its spine). The poloidal current density J_p is therefore downward in the inner regions and outward in the outer jet regions. As a consequence, field lines anchored up to r_o ∼ 25 are focused to the polar axis (Z-pinch due a pole-ward J_p × B_ϕ force as J_z < 0), while field lines anchored beyond r_o ∼ 30 are de-collimated and pulled out (because of the outward action of the same J_p × B_ϕ force as J_z > 0). See the right panel of Fig. 2 for the topology of J.

Fig. 8. Evolution of several magnetic field lines during the simulation computation, for different anchoring radii R0 and for the reference simulation K2. The field lines at the first output of the simulation (initial conditions) are shown in dotted lines. The field lines at the last output of the simulation (final state) are shown as full lines.

The inner self-collimated jet region can only exist thanks to the existence of these outer jet regions that are pulled back and out. The final state of the jet collimation, namely the asymptotic jet radius achieved by these inner regions (the densest ones, possibly responsible for the observed astrophysical jets), is therefore also a consequence of the transverse equilibrium achieved by the outflow outskirts with the ambient medium. This balance is mathematically described by the Grad-Shafranov equation and expresses the action of the poloidal electric currents – how they are flowing and how electric circuits are closed within the jet – on the shape of the magnetic surfaces (Heyvaerts & Norman 1989, 2003a,b; F97; Okamoto 2001). We return to this point below.

3.4. Standing recollimation shocks

Our reference simulation K2 ends with five standing shocks, of which only the first (red in Fig. 4), starting on the axis around Z₁ = 1850, is studied hereafter. This is for two reasons: (1) Contrary to the orange and cyan shocks, the red shock surface covers a large extension in the jet itself. It is therefore most probably related to the dynamics of our quasi self-similar jet, whereas these smaller shocks may be related to the spine–jet interaction. (2) It remains far away from the outer boundary (it ends up at the FM surface at a point r = 2500, z = 3800), which is clearly not the case for the purple shock and also probably not for the blue one.

3.4.1. Rankine-Hugoniot conditions

All shocks are due to the flow heading toward the axis (with decreasing, usually negative B_r and u_r components) at a super-FM velocity, resulting in a sudden jump in all flow quantities with an outward refraction of the magnetic surface (this can be seen in Fig. 8). The Rankine-Hugoniot jump conditions (see Appendix C for more details) are of course satisfied with the shock-capturing scheme of PLUTO and MHD invariants are conserved (up to some accuracy, as discussed above).

The left panel of Fig. 9 displays the evolution of several quantities along the red shock surface. It can be seen that the compression rate χ (green curve), defined as the ratio of the post-shock to the pre-shock densities, is larger than unity while remaining small (≤1.3, see red curve in the right panel), despite a very large Alfvénic Mach number m ∼ 100 (orange). This is probably because the shock is oblique and the incoming jet material does not reach very large super-FM Mach numbers n_⊥ (defined as the ratio of the flow velocity normal to the shock surface to the FM phase speed in that same direction; blue curve). In fact, n_⊥ ≤ n = u_p/V_fm, and at these distances the jet has reached its asymptotic state with a maximum velocity of $u_p\sim \sqrt{2\lambda -3}{V}_{\mathit{Ko}}$ (V_Ko = Ω_*r_o the Keplerian speed at the anchoring radius). Assuming B_p ≪ |B_ϕ|, m ≫ 1 and a jet widening such that r ≫ r_A where r_A is the Alfvén cylindrical radius along a flow line where m = 1, leads to

$$\begin{array}{c}\hfill n\simeq \frac{u_p}{V_{A\varphi }}=\frac{u_p}{{\mathrm{\Omega }}_{\ast }r}\frac{m}{1-r_A^2/r^2}\simeq m\frac{r_o}{r}\sqrt{2\lambda -3}\sim m\frac{r_A}{r},\end{array}$$(16)

Fig. 9. Distributions along the shock of several quantities for K2 at tend. Left: normal incident FM Mach number n⊥, Alfvénic Mach number m, compression rate χ, relative variations of the toroidal magnetic field δBϕ and plasma angular velocity −δΩ, and total deviation δi (in rad) of the poloidal magnetic field line for the main recollimation shock. Right: compression rates χ of all shocks appearing in Fig. 4, using the same color code. The main shock corresponds to the red curve. See Appendix C for more details on this figure.

which shows that the asymptotic FM Mach number critically depends on how much the jet widens (see Pelletier & Pudritz 1992 and Sect. 5 in F97). In our case, n ∼ 4 at the outer edge of the spine–jet interface, in agreement with self-similar studies. Following the main red shock along growing r, the incident angle³ of the magnetic field lines on the shock front decreases until turning into a normal shock on its external edge (r ∼ 2000). Hence, on this edge, the shock front coincides with the FM critical surface n = 1. As the shock becomes normal, n_⊥ → n = 1 and the shock vanishes, with a compression rate χ going to 1.

The three other curves in the left panel of Fig. 9 describe other modifications in jet dynamics. The brown curve is the magnetic field line deviation at the shock front, δi = i₂ − i₁, where i is the flow incidence angle to the normal to the shock surface (subscripts 1 and 2 refer to the pre- and post-shock zones, respectively). The maximum deviation of 0.07 rad = 4° is very small, in agreement with the small compression rate.

The purple curve describes the relative variation of the flow rotation δΩ = (Ω₂ − Ω₁)/Ω₁, which is always negative. The shock introduces a sudden brake in the azimuthal speed, meaning that the compressed shocked material is always rotating less. As t.e detection of rotation signatures in YSO jets is an important tool for retrieving fundamental jet properties (see e.g., Anderson et al. 2003; Ferreira et al. 2006; Louvet et al. 2018; Tabone et al. 2020), recollimation shocks appear to be a very interesting means to lower the jet apparent rotation. However, the rather weak shock found here only introduces a decrease of ∼20% at the outer edge of the shock.

The plasma loss of its angular momentum at the shock is of course compensated for by a gain of magnetic field (the angular moment is a MHD invariant). This means that the magnetic field lines are more twisted after the shock than before, as illustrated in the red curve showing δB_ϕ = (B_ϕ2 − B_ϕ1)/B_ϕ1 > 0. The shock surface acts therefore as a current sheet with an electric current density flowing outwardly (in the spherical R direction).

3.4.2. Two families of shocks

The right panel of Fig. 9 displays the compression rate χ for the five shocks seen in Fig. 2, using the same color code. All shocks but the red one have larger compression factors near the axis. The orange and cyan shocks merge with the main red one (respectively at r ∼ 500 and r ∼ 900), leading to an increase in its compression rate χ. The large blue shock has the same behavior as the orange and cyan but remains alone (i.e., not merging with the red) with χ converging to 1, while the purple shock seems to have a similar behavior to the red one, maintaining a larger value for χ. Although these last two shocks are probably affected by their proximity with the outer boundary, it seems that two classes of recollimation shocks are at stake.

The first class (represented by the red and purple shocks) corresponds to the recollimation shock predicted in self-similar studies (FP97, Polko et al. 2010). The reason for their existence is the hoop-stress that becomes dominant as the jet widens, leading to a magnetic focusing toward the axis. As shown in FP97, such a situation always arises in the super-FM regime, meaning that a shock is the only possibility for the converging flow to bounce away. However, as long as no dissipation is introduced, such a situation will repeat. Indeed, after the flow refraction due to the shock, the magnetic field starts to accelerate the plasma again, the magnetic surface widens and the same situation repeats. One would therefore expect periodic oscillations and shocks on a vertical scale H_R (measured on the axis). Figure 8 provides some evidence of this pattern for the field lines anchored at r_o = 20, 30 or 40. The first recollimation shock (red) is quite far away from the disk (Z_shock = 1850), but the second shock (purple) occurs at Z₅ ≃ 2634. A much larger computational domain would be necessary in order to clearly assess a periodic pattern.

The second class of recollimation shocks (represented by the orange, cyan, and blue shocks) are limited to the vicinity of the spine–jet interface and are thereby a consequence of a radial equilibrium mismatch between these two super-FM outflows. The transverse equilibrium of a magnetic surface is provided by projecting the stationary momentum equation in the direction perpendicular to that surface, leading to the equation (FP97)

$$\begin{array}{c}\hfill (1-m^2)\frac{B_p^2}{{\mu }_0\mathcal{R}}-{\mathrm{\nabla }}_{\perp }(P+\frac{B^2}{2{\mu }_0})-\rho {\mathrm{\nabla }}_{\perp }{\mathrm{\Phi }}_G+(\rho {\mathrm{\Omega }}^{2}r-\frac{B_{\varphi }^2}{{\mu }_{0}r}){\mathrm{\nabla }}_{\perp }r=0.\end{array}$$(17)

Here, ∇_⊥ = e_⊥ ⋅ ∇ provides the gradient perpendicular to a magnetic surface with e_⊥ = ∇Ψ/|∇Ψ| and $B_p^2/\mathcal{R}={\mathit{e}}_{\perp }\cdot ({\mathit{B}}_p\cdot \nabla {\mathit{B}}_p)$, measures the local curvature radius ℛ of the magnetic surface. In the asymptotic region where these small recollimation shocks are observed, the field lines are almost vertical and gravity is negligible. The above equation therefore reduces to

$$\begin{array}{c}\hfill -\frac{\partial }{\partial r}(P+\frac{B^2}{2{\mu }_0})+\rho {\mathrm{\Omega }}^{2}r-\frac{B_{\varphi }^2}{{\mu }_{0}r}=0.\end{array}$$(18)

Looking at Fig. 10, where the various forces are plotted as a function of the cylindrical radius at a constant height Z = 2400, it can be clearly seen that the dominant force is the hoop stress $-B_{\varphi }^2/({\mu }_{o}r)$ (purple curve) near the spine-jet interface (dashed vertical line), defined as the magnetic field line anchored at r_o = R_d. That pinching force overcomes the others (the above sum is actually nonzero and negative), indicating that the field lines do have a curvature and are actually converging toward the axis (therefore Eq. (18) is too simplistic). Nevertheless, this behavior of the forces is consistent with the stationary shape of the magnetic surfaces seen at Z = 2400 at those radii. Further out, downstream, the fifth (purple) shock will make the field lines bounce back again. This tells us that we are witnessing radial oscillations of the radius of the spine driven by a mismatch between the dominant forces (hoop stress, magnetic pressure, and centrifugal term).

Fig. 10. Radial distribution of the radial accelerations and their sum at the altitude Z = 2400 for the simulation K2 at tend. The vertical dashed line corresponds to the spine–jet interface, namely the field line anchored at ro = Rd.

This oscillatory behavior may be a generic feature of MHD outflows from a central rotator as shown by Vlahakis & Tsinganos (1997), because of the different scaling with the radius of the pinching force and the centrifugal force (as in the self-similar jet; F97). But it may also be triggered by the pinching due to the outer jet recollimation, namely the spine–jet interface response to the global jet recollimation. Indeed, no spine–jet shock is seen before the onset of the main recollimation shock. We further note that the five shocks are located at a slightly increasing distance from each other. Indeed, ΔZ₁₂ = Z₂ − Z₁ = 150, ΔZ₂₃ = 160, ΔZ₃₄ = 212, and ΔZ₄₅ = 262, which is the sign of some damping of the spatial oscillations at the spine–jet interface (see e.g., Vlahakis & Tsinganos 1997). This corresponds to three spine–jet shocks (orange, cyan, and blue) located between the two large jet recollimation shocks (red and purple), as can be seen in Fig. 4.

Despite the fact that the magnetic surfaces are in a steady state, it is useful to look at this spatial oscillatory pattern as the nonlinear outcome of transverse waves. Immediately after a shock, the flow is again outwardly accelerated leading to a widening of the magnetic surface and its refocusing toward the axis with the unavoidable shock. It will therefore take a time Δt_z = ΔZ/u_z to reach the next shock. On the other hand, any radial unbalance triggered immediately after the shock gives rise to a transverse (radial) FM wave that bounces back on a timescale of Δt_r = 2r/V_fm measured at the spine–jet interface. In steady-state, these two times must be the same, which requires ΔZ ∼ 2nr, where n is the FM Mach number measured at the spine–jet interface. At Z = 2400, the width of the spine is r ≃ 40, and n ≃ 3, providing the correct order of magnitude for ΔZ ∼ 240. This is also verified for all other shocks. This correspondence strengthens the idea that these small shocks are actually triggered by the first large recollimation shock.

3.5. Electric circuits

The existence of these shocks drastically affects the poloidal electric circuits that go along with MHD acceleration (and of course collimation). This can be seen in Fig. 11, where several interesting circuits are shown in color. Each poloidal circuit corresponds to an isocontour of rB_ϕ, the arrows indicating their flowing direction.

Fig. 11. Plot of the poloidal electric circuits at tend for simulation K2. The two red curves are the critical surfaces, Alfvén (dotted) and FM (dashed). The yellow curves are the poloidal electric circuits, defined as isocontours of rBϕ, where the arrow indicate the direction of the poloidal current density Jp. Four circuits are highlighted in particular: (1) the envelope of the inner accelerating current in white (rBϕ = −2.06), (2) the outermost circuit still fully enclosed within the domain in blue (Bϕ = −2.005), (3) a circuit closed outside the domain in orange (rBϕ = −1.80), and (4) a post-shock accelerating circuit in purple (also with rBϕ = −2.06).

The white contour marks the last electric circuit that flows below the first recollimation shock and defines the envelope of the initial accelerating current. It links the disk emf with the accelerated jet plasma and flows back to the disk along the spine. This is due to the fact that the largest electric potential difference is with the axis.

The blue circuit is the last electric circuit fully enclosed within the computational domain. The current flows out of the disk (further away than the previous circuit) and makes a large loop that goes beyond (downstream) the main recollimation shock, flowing back on the axis below Z₅ until it encounters the smaller shock at Z₄ (blue curve in Fig. 4). As a shock behaves as an emf, with an outwardly (positive J_r) electric current flowing along its surface, the blue electric current gets around it and goes back to the axis where it meets the next shock surface at Z₃ (cyan in Fig. 4). As that shock merges with the main recollimation shock, the blue electric current flows along these two surfaces, gets around the main shock (near the point r ∼ z ∼ 3000 in Fig. 11) and returns back to the disk via the spine, where it joins the white circuit below Z₁.

As mentioned above, the outflowing plasma gets re-accelerated after each shock. This requires a local accelerating electric circuit, which is naturally enclosed within two recollimation shocks. One such circuit is exemplified by the purple contour in Fig. 11. It actually has the same rB_ϕ value as the white one, but is enclosed between the two large recollimation shocks. As the small shocks (orange and cyan in Fig. 4) merge with the main one, the purple accelerating circuit is the envelope of the current used to go from the first main shock to the second one flowing back to the spine just before Z₅ (and getting around the shock starting at Z₄ near the point r ∼ 2000, z ∼ 4000, like the previous blue electric circuit).

These three examples of electric circuits (white, purple, and even blue) are fully closed within the computational domain and therefore do not contribute to any further asymptotic collimation. However, it can be seen that the electric current outflowing from the disk beyond r_o ∼ 2000 leaves the computational domain and is supposedly closed by the inflowing current that enters the computational domain at small colatitudes. One example of such an electric circuit is shown by the orange curve in Fig. 11. This inward electric current is responsible for the inner jet collimation at large distances, say at Z ≳ 3000, and is seen to flow back to the disk along the spine, which acts as a conductor. This implies that the asymptotic jet collimation is here somewhat controlled by an electric circuit that is not fully self-consistent. Indeed, this electric circuit is actually determined by the boundary conditions at R_ext and there is no guarantee that its evolution is consistent with the disk emf. This is of course unavoidable but it may have an impact on the collimation properties of numerical jets (see discussion in Sect. 5). Moreover, as this current embraces very large spatial scales, very long timescales are consistently implied and may lead to evolution of the jet transverse equilibrium on those scales.

3.6. Time evolution

Figure 12 shows the evolution in time of the vertical height Z (measured at the axis) of all shocks found in the simulation. As discussed previously, it takes a time t_ext(1)∼10³ for the innermost jet (anchored at r_o = R_d) to reach the boundary of our computational box. During this early evolution with t < t_ext(1), the detected shocks correspond to the first bow shock where the jet front meets the initial unperturbed ambient medium. Once the jet has reached the boundary, the spine is in a steady-state while the jet transverse equilibrium continues to evolve due to the self-similar increase in its width. Indeed, the time to reach the outer boundary for a magnetic surface anchored at r_o in the disk grows as $t_{\text{ext}}(r_o)\propto r_o^{1/2}$. It therefore takes a time ∼10⁴ to achieve a steady ejection from r_o = r_o, FM = 323, where r_o, FM is the maximum radius of the field lines that achieve a super-FM flow speed. The global jet structure is only expected to have a steady state beyond that time.

Fig. 12. Altitude Z of the different shocks (measured at the axis) as a function of time (in Td units) for simulation K2. The vertical lines correspond to the six times ti used in Fig. 13: t1 = 551, t2 = 2.08 × 103, t3 = 8.51 × 103, t4 = 1.99 × 104, t5 = 1.05 × 105, and t6 = 1.58 × 105. The last vertical line is tend = 6.51 × 105.

The vertical lines in Fig. 12 trace six times t₁ = 551, t₂ = 2.08  ×  10³, t₃  =  8.51  ×  10³, t₄  =  1.99  ×  10⁴, t₅  =  1.05  ×  10⁵, and t₆  =  1.58  ×  10⁵. The snapshots corresponding to each of these times are shown in Fig. 13, allowing us to see the global jet evolution. The times t₁ and t₂ have been chosen to enclose t_ext(1). The bow shock with the ambient medium is clearly seen, as is the outward (radial) evolution of the jet width. At t₂, several shocks near Z ∼ 2000 can be seen in both figures. The jet radial equilibrium is clearly not yet steady. However, Fig. 12 shows that while the positions of the shocks (as measured on the axis) are already close to their final values, their final number is not yet settled.

Fig. 13. Snapshots of our reference simulation K2 at different times (given in Td units). From top to bottom, left to right: t1 = 551, t2 = 2.08 × 103, t3 = 8.51 × 103, t4 = 1.99 × 104, t5 = 1.05 × 105, and t6 = 1.58 × 105. The background color is the logarithm of the density, black lines are the magnetic surfaces, red lines the Alfvén (dotted) and FM (dashed) surfaces, and yellow curves are isocontours of the poloidal electric current.

Four standing recollimation shocks seem to settle somewhere between t₃ = 8.51 × 10³ and t₄ = 1.99 × 10⁴, in agreement with our previous estimate; they can be clearly seen in Fig. 13, where some transient shocks located further up at Z > 4000 at t₃ have disappeared at t₄. Also, given the huge spatial scales involved, most of our JED is still evolving. For instance, while there have already been 3183 orbital periods at R_d at t₄, the disk has done only half of an orbit at r_o = 323. This can be seen in the shape of the FM surface, which has not yet reached its steady-state configuration (conical).

Beyond t₄, the global flow is slowly evolving in time in some adiabatic way, with four standing recollimation shocks. The evolution of the jet outer regions and the progressive evolution of the A and FM surfaces to their conical shapes seems to produce no obvious evolution in the shocks until t₅. At that time, a dramatic evolution is triggered with the appearance of shocks beyond Z₄. Figure 12 clearly shows this pattern with the appearance of a fifth shock; its altitude Z₅ evolves in time, consistently with Z₄ until a steady-state is finally reached approximately near t₆. Our final state t_end shows no relevant difference in the positions of the five shocks. This evolution of the distribution of the shocks has only slightly affected the position of the farthest shock Z₄, leading to the final regular distance ΔZ discussed above.

The appearance of a fifth shock at t₅ leading – after a transient phase ending at t₆ – to a new steady state jet configuration is illustrated in Fig. 14. The time evolution of the cylindrical radius of the field line anchored at r_o = 3 (in blue) is measured at a constant height Z = 3500. It can be seen that this radius is steadily slowly decreasing in time, going from r ∼ 150 near t₄ down to r ∼ 135 at t₅ where some fluctuations are suddenly triggered. These oscillations describe a time-dependent behavior which ends at t₆, with a new radial balance found at a smaller radius r ≃ 125. Globally, this evolution describes a magnetic surface that is first slowly getting more and more confined, and then enters an unstable situation until the formation of another (tighter) equilibrium. This behavior is consistent with the evolution of the electric current (red curve) that flows within that magnetic surface, which is seen to first steadily (although very slowly) increase until achieving a final value.

Fig. 14. Time evolution of the cylindrical radius r measured at Z = 3500 of the magnetic surface anchored at ro = 3 (blue curve) and the electric current I = rBϕ (red) flowing within that surface for simulation K2. The two vertical dashed lines correspond to t5 = 1.05 × 105 and t6 = 1.58 × 105.

In these inner regions, transverse equilibrium of the cold jet is mostly achieved by the poloidal magnetic pressure balancing the toroidal pressure and hoop stress (see Eq. (18) and Fig. 10). In this Bennett relation, one gets

$$\begin{array}{c}\hfill {(r{B}_{\varphi })}^2={(r{B}_z)}^2+{\displaystyle {\int }_0^{r}2B_z^{2}r\mathrm{d}r.}\end{array}$$(19)

This relation shows that the toroidal field B_ϕ is compelled to follow the same scaling as B_z in order to maintain the jet transverse equilibrium. If we assume that the self-similar radial scaling for the vertical field B_z is recovered at large distances, we get (rB_ϕ)² ∝ r^{2(α − 1)}. As a consequence, for α < 1 (which is the case here), whenever the electric current is increasing, the radius of the magnetic surface decreases (as in a Z-pinch). This scaling provides ΔI/I = (α − 1)Δr/r = −0.25Δr/r, which is consistent with the evolution seen in Fig. 14.

This slow increase in time of the electric current flowing near the axis is a natural consequence of the increase of the disk emf e_disk as the outer disk regions achieve a steady state. Indeed, $e_{\text{disk}}\simeq {\int }_{R_d}^{r_{\text{max}}}{\mathrm{\Omega }}_{K}r{B}_z\text{d}r$ increases with r_max, which is the maximal radius in the disk that achieved a steady state. This increase in the available current is expected to stop when no further relevant emf is added. The available current is determined at the disk surface by the crossing of the Alfvén critical point, because it is that point that fixes the available total specific angular momentum carried away. One can therefore estimate the time where the current should level off as the time when the outermost magnetic field line reached the Alfvén point. Figure 14 shows that this is achieved approximatively around t₆ with r_max ∼ 10³, corresponding to a full orbital period at r_max equal to $t=2\pi r_{\max }^{3/2}=2\times {10}^5$ (we note that the time for the flow ejected at r_max to reach the outer boundary is comparable). After a time of a few 10⁵ T_d, our simulation has finally achieved a global steady state.

4. Parameter dependence

In the previous section, we showed that our reference simulation K2 behaves qualitatively like the self-similar analytical calculations of Blandford & Payne (1982) and F97, with a refocusing toward the axis and the formation of a recollimation shock. However, the presence of the axial spine breaks down the self-similarity and introduces additional shocks localized at the spine–jet interface.

To further understand the behavior of these shocks, we conducted a parameter study in κ and α, the jet mass load and the radial exponent of the disk magnetic flux, respectively. We finally ran one simulation with the same JED parameters as in our reference simulation K2, but with a rotating spine. All our simulations are presented in Table 1.

4.1. Influence of the mass loading parameter κ

In this section, we present our simulations K1 to K5, obtained with the same parameters as K2 except for the mass load κ which is varied from κ = 0.05 to κ = 1. Our parameter range in κ is slightly smaller than the one achieved by Blandford & Payne (1982) which goes down to κ = 0.01. It can be seen in Fig. 15, which represents the magnetic lever arm λ as a function of the mass-loading parameter κ. Each simulation is obtained for a unique value of κ but as the simulation is not strictly self-similar, we obtain a range in λ: the larger the anchoring radius, the larger the field line inclination at the disk surface and the larger the magnetic lever arm λ. To ease the comparison with Fig. 2 of Blandford & Payne (1982), for each simulation, we computed the anchoring radius r_o at which the field line inclination ${\xi }_0^{\prime }=B_r/B_z$ at the disk surface is equal to 1.4, 1.5, and 1.65. This allowed us to draw in our Fig. 15 iso-contours of ${\xi }_0^{\prime }$, which are in agreement with the above expectations and analytical self-similar calculations (see also Fig. 3 in F97).

Fig. 15. Jet parameter space λ(κ) at the final stage of our simulations K1–K5 with α = 3/4. Each simulation is obtained for a unique mass loading κ and gives rise to a distribution of the magnetic lever arm λ with the radius: green, red, and blue dots correspond to anchoring radii ro = 5, 50, and 500, respectively. The solid curves are obtained for constant values (indicated in the panel) of the initial magnetic field inclination ${\xi }_0^{\prime }=B_r/B_z$ at the disk surface.

All simulations achieve a steady-state and exhibit roughly the same behavior as K2, as can be seen in Fig. 16. From top to bottom, κ increases from 0.05 to 1, the left panels showing the whole computational domain with the two critical surfaces in red (A, dotted and FM dashed) and the right panels providing a close up view around the shock-forming regions near the axis. Table 1 provides the value of several jet quantities: the spine mass-loss rate stays around 10% of the jet mass-loss rate, despite the net increase in Ṁ_jet (∝κ, factor 20 increase). Similarly, while the jet power P_jet scales in κ, the spine power stays around 80% of the jet power.

Fig. 16. Influence of the mass-loading parameter κ on the final stage of jets obtained with α = 3/4. The color background is the logarithm of the FM Mach number n, black solid lines are field lines, yellow lines are isocontours of the electric current rBϕ and the red dashed (resp. dotted) curve is the FM (resp. Alfvén) critical surface. The left panels show the whole domain and the right panels a close-up view around the shock-formation regions. In the left panel, the field lines anchoring radii are ro = 3; 15; 40; 80; 160; 320; 600; 1000; 1500. In the right panel, the field lines anchoring radii are ro = 1.2; 2; 3; 4; 5; 7; 9; 11; 13; 15.

The first observation is that the altitude of the main recollimation shock (the one merging with the FM surface) decreases globally with κ. This is quantitatively shown in Fig. 17, where Z_shock moves from 2150 down to 700. The same evolution occurs for the altitude Z_tip where the main shock merges with the FM surface. Globally, as κ increases, the whole jet structure decreases toward the disk.

Fig. 17. Influence of κ (orange lines) and α (blue lines) on the altitude of the main recollimation shock. This is done by measuring two altitudes for each shock: its height at the axis (Zshock, solid lines) and the altitude of its outer edge (Ztip, dashed lines). The scale for κ is indicated below, while the scale for α is above.

Such behavior is consistent with the self-similar calculations obtained by F97. Indeed, as evidenced in his Fig. 6, the denser the jet (larger κ), the sooner (smaller altitudes) recollimation takes place. This can be understood qualitatively by the fact that λ = 1 + q/κ, where q = |B_ϕ/B_z| is the magnetic shear at the disk surface (F97). Now, as the mass load κ increases, the magnetic lever arm λ must decrease (see also Fig. 15). This translates into magnetic surfaces that open less, a less efficient magneto-centrifugal acceleration, and recollimation shocks that are not only closer to the disk but also show a smaller compression rate χ due to a smaller FM Mach number n.

However, the physical scales implied are very different. Here, a factor 20 difference in κ leads to a decrease in Z_shock by a mere factor 3, with a minimum value of 700. In F97, the mass-loading parameter κ ∼ ξ varies from 0.01 to 0.05 only (for a constant disk aspect ratio, see his Fig. 3) but leads to variations in recollimation altitudes that span six decades. Our lowest Z_shock obtained for κ = 1 is still much farther away than the minimum height of ∼10 found for κ ∼ 0.05 in F97.

This discrepancy can of course be attributed to the very different injection properties. Indeed, our numerical simulations assume a supersonic flow while the self-similar calculations compute the disk structure and outflows are found to only be super-SM (and still subsonic) at the disk surface. However, our guess is that the huge difference in the shock position is probably due to the existence of the spine, which breaks down the self-similarity. Indeed, recollimation is due to the dominant hoop stress, and while in our case all quantities are leveling off on the axis, strictly self-similar solutions have an axial electric current that grows without limits. For instance, at a cylindrical distance r = 0.1 from the axis at the spine basis, our B_z remains comparable to the disk field, B_ϕ goes to zero, and the normalized Bernoulli integral e has decreased by a factor 5 (see Fig. B.1). Self-similar solutions, on the contrary, have fields and a Bernoulli integral increasing respectively by a factor 10^5/4 = 17.8 and 10. This suggests that the conditions assumed on the axis most certainly affect the overall jet collimation properties. We return to this aspect later on by changing the spine properties.

A second interesting aspect is the appearance of a second ensemble of shocks arriving at higher altitudes, namely Z > 3000 for κ = 0.5 and Z > 2000 for κ = 1. This second group of shocks is also composed of two large recollimation shocks separated by smaller ones. The distance (measured at the axis) between the two large shocks is comparable to the width of the first group of shocks. As discussed in Sect. 3.3, this hints to the fact that each group is caused by the global jet recollimation dynamics, which should be periodic in the Z-direction on a scale H_R. Looking at the simulation K4 with κ = 0.5, this would give H_R ∼ 1650, while the width of each group is around W ∼ 300 (the first group of shocks being located between 1150 and 1450, and the second between 3100 and 3400). As long as no dissipation is introduced, such a periodic behavior should continue in a box of infinite size.

4.2. Influence of the magnetic field distribution α

In this section, we vary the magnetic flux distribution exponent α in our simulations A1 to A5, keeping all other parameters (see Table 1) constant. A strict mathematical self-similarity links the magnetic field distribution α with the disk density in such a way that α = (12 + 8ξ)/16, where ξ is the disk ejection efficiency and is related to the disk accretion rate Ṁ_a(r) ∝ r^ξ (Ferreira & Pelletier 1995). As long as material is only outflowing from the disk (namely ξ > 0) and jet power is only released from accretion (namely ξ < 1, F97), this leads to the unavoidable constraint 12/16 < α < 20/16.

Only α = 12/16 can be compared to the cold jet solutions of F97, assuming an ejection efficiency ξ < 0.1. Larger values of α > 12/16 would require a disk ejection efficiency of ξ = 0.125 up to 0.5. These values are only achievable in analytical studies by introducing an additional heat deposition at the disk surface (magneto-thermal flows, Casse & Ferreira 2000b) and/or a much smaller magnetic field strength (Jacquemin-Ide et al. 2019). However, the physical processes required to get these solutions are missing in our simulation setup. We nevertheless vary α in order to allow a comparison with the self-similar jet solutions found by Contopoulos & Lovelace (1994), who tested the effects of α ranging from 0.5 to 1.02. We did not succeed to numerically obtain steady-state solutions for α ≥ 1 for reasons that are discussed below, and show only simulations with α ranging from 10/16 = 0.625 to 15/16 = 0.937.

All our A1 to A5 simulations reach a steady state with the same overall behavior as K2, namely the existence of two main recollimation shocks within the computational domain, separated by several smaller standing shocks located at the spine–jet interface. This can be clearly seen in Fig. 18 where, on the other hand, a trend with α can be observed. Indeed, the radial extension of the shocks decreases with α. This can also be seen in Fig. 17 which shows that the altitude Z_tip of the main recollimation shock decreases with α, while its altitude Z_shock at the axis barely changes. For simulations A1 (α = 10/16) and A2 (α = 11/16), the main recollimation shock ends out of the box, and therefore Z_tip cannot be defined. However, the value of Z_shock remains similar (see Table 1). This is a geometrical effect due to the fact that, as α increases, the magnetic field configuration goes from a highly inclined magnetic configuration (α = 0 corresponds to a monopole) to one that is much less inclined (α = 2 is a purely vertical field). This can be seen in the shape of the magnetic field lines (black solid lines) in Fig. 18. This geometrical effect translates into a smaller incidence angle near the axis and therefore to weaker shocks (the incidence becomes normal and n_⊥ decreases to unity).

Fig. 18. Influence of the magnetic field distribution α on the final stage of jets obtained with κ = 0.1. We use the same notations, colors, and field lines anchoring radii as in Fig. 13.

For α = 10/16 and 11/16, the shocks still exist but the MHD characteristics are much more vertical than in K2. A larger box would probably be necessary to recover the K2 behavior. The opposite trend can be seen for α = 14/16, with less vertical MHD characteristics allowing now the second main recollimation shock to merge with the FM surface within the domain.

Table 1 also shows that as α increases, the last radius on the disk giving rise to a super-FM flow increases and the colatitudes (measured at the outer boundary) ${\theta }_{\mathrm{FM}}^{\mathrm{ext}}$ and ${\theta }_A^{\mathrm{ext}}$ of the critical FM and A surfaces decrease. These results are a natural consequence of the magnetic field distribution becoming more vertical as α increases. As the jet mass loss Ṁ_jet and jet power P_jet are computed up to r_0,FM which increases with α, Ṁ_jet and P_jet increase with α. Still, the mass loss increases even when computing up to a fixed radius. Indeed, the density decreases less with an increasing α (ρ ∝ r^{2α − 3}) and the outer disk regions contribute more to the mass flux and jet power. Moreover, as the distribution in density is flatter with an increasing α, it is only natural that Ṁ_spine/Ṁ_jet and P_spine/P_jet decrease when α increases.

In summary, we find that the altitude of the shocks barely changes with α, which is in strong contrast with Contopoulos & Lovelace (1994). Indeed, their Table 1 shows that as α increases from 0.5 to a critical value 0.856, their self-similar jet becomes super-FM and undergoes a recollimation at a distance that increases by several decades (as in F97). As discussed above, we believe that this discrepancy is due to our nonstrict self-similar scaling (which forbids the unlimited growth of the inner electric current and the subsequential Z-pinch in self-similar solutions) and the presence of the spine. Contopoulos & Lovelace (1994) also report that their solutions with α > 0.856 remain sub-FM, while we clearly achieve super-FM flows up to α = 0.937 = 15/16. This is again probably a difference in our jet radial balance, leading to a slightly different jet acceleration efficiency. Finally, their solutions with α = 1, 1.01, and 1.02 remain sub-FM but evolve through a series of radial oscillations at logarithmically equal distances in Z.

Our simulation A5 is the simulation with the largest value α = 15/16, of close to unity. Although the final integration time t_end is quite comparable with the other simulations, its appearance clearly shows that the global configuration is still far from achieving a steady state like K2. This can be seen in Fig. 18, in the shape of the critical surfaces but also on the isocontours of the poloidal electric circuit (yellow curves) that are still struggling to find their final state. This is normal and points to a numerical difficulty in computing MHD codes when the magnetic configuration is such that α ≥ 1.

Indeed, it has been argued above that for self-similar boundary conditions, the jet transverse balance imposes a toroidal magnetic field scaling with the vertical magnetic field. This leads to an electric current at the disk surface behaving as I = rB_ϕ ∝ r^α − 1. Magnetic configurations with α < 1 correspond to a poloidal current density leaving the disk surface and closing along or near the spine (where it flows back to the disk), whereas configurations with α > 1 correspond to an inward poloidal current density, with current closure being done only at the outskirts of the outflow (F97). As div J = 0 in ideal MHD, all electric circuits must be closed. Let us define a radius in the disk r_I such that for r < r_I, the electric current flows down into the disk whereas it flows out of it for r > r_I. For α < 1, r_I is always larger but close to R_d implying very short timescales. As discussed in Sect. 3.6, as time is evolved, the outer disk regions provide more current that struggles to reach the innermost disk radius. But a global radial balance can be achieved consistently with the electric current closure condition because the local time near r_I is small. On the contrary, configurations with α > 1 have r_I that is constantly increasing in time (as t^2/3), leading to an electric circuit that freezes in time and therefore to a transverse MHD balance that takes a much longer time to achieve steady state. We observe this behavior for all values of α > 1 and none of these simulations has achieved a steady state.

The limiting value α = 1 (close to our A5 simulation) would correspond to r_I = R_d and absolutely no electric current flowing out of the disk until some outer radius. Current closure could only be done through the spine and the outer jet interface with the ambient medium. But then, no magnetic acceleration would be possible as no electric current could be used along the magnetic surfaces. This would correspond to an exact force-free field configuration fully determined by the chosen boundary conditions. To compute such a situation, boundary conditions must be designed where the jet launching region has a finite extent. However, we doubt that the outcome would be a force-free field unless explicitly enforced. This is beyond the scope of this paper.

4.3. Influence of a rotating central object

In this section, we do not intend to fully explore the physical parameters of the spine but instead wish to probe whether the spine, despite its small spatial extent, has indeed a profound influence on the overall jet dynamics. All simulations K1-K5 and A1-A5 were done with the same nonrotating central object in order to minimize its emf and numerically follow the outcome of a jet emitted from an outer self-similar disk (see Sect. 2.4.2). Our choice of parameters gives rise to a spine carrying typically 10% of the mass flux, and therefore providing only a small contribution to the overall outflow. Nevertheless, this spine carries a large fraction of the emitted power, even superior to that of the disk for the simulation A1. The spine plays an important role in introducing extra standing shocks at its interface with the jet, but is probably also determining the altitude where the first large recollimation shock occurs. Indeed, as discussed in the previous sections, the amount of electric current that is flowing along the innermost axial regions (along the spine and the inner jet) is what determines the strength of the Z-pinch acting upon the jet and thereby the altitude Z_shock.

In order to probe this idea, we ran another simulation with a rotating object (simulation SP in Table 1). We chose an object rotating at the same angular velocity as the innermost disk radius R_d, namely Ω_a = Ω_Kd. This is for instance representative of a star–disk interaction where the disk truncation radius is located at the co-rotation radius. By doing so, the emf due to the central object becomes non-negligible and we expect a stronger poloidal electric current. However, care must be taken as enhancing the hoop stress may also lead to an overwhelming radial pinch. To prevent this and get somewhat closer to the self-similar conditions, we also increase the value of the Bernoulli integral on the axis and use e_a = 10 (e_a = 2 for the other simulations). We note that the Bernoulli invariant from the innermost disk region is e_d ≃ λ_d − 3/2 ∼ 10. This translates mostly into a thermal pressure that is five times larger than previously. Thus, our new conditions for the spine provide a rotation and a specific energy that are only comparable to those at the inner jet, not much larger as in a self-similar situation.

Figure 19 shows the final outcome of this new simulation. It achieves a global steady state with the same features as in our reference simulation K2. However, the shocks are, as expected, localized at lower altitudes, allowing a second set of large recollimation shocks to appear near Z ∼ 3700. The first large recollimation shock appears at Z_shock = 1300, which is significantly smaller than Z_shock = 1850 obtained for a nonrotating central object. This result confirms the role of the central object in shaping, through its spine, the collimation properties of the jets emitted by the surrounding disk. This is very promising and deserves further investigation.

Fig. 19. Snapshot at tend of our SP simulation with a rotating spine, α = 3/4 and κ = 0.1. We use the same color coding as in Fig. 2. The magnetic field lines (black solid lines) are anchored at the same disk radii.

5. Discussion

5.1. Caveats

This paper provides some novel information on the collimation of jets emitted from self-similar magnetized disks. However, we would like the reader to pay attention to several caveats:

As we tried to replicate a Blandford & Payne process, we mostly explored the ejection conditions on the disk and not on the source, with the exception of the rotation (see Sect. 4.3). For instance, we kept a value of sonic mach number M_S = 10 on the whole ejection zone, thus only creating cold jets. It would be interesting to provide a real exploration of the spine parameters to better understand the role of the central object in the collimation.

The shocks displayed in our simulations are rather weak (χ ∼ 2 near the axis, χ ∼ 1.3 in the jet). This is probably due to our isentropic scheme, but as the simulations only reach small values of n, the compression rates are intrinsically low (see Appendix C). Moreover, as shown in Eq. (C.4) the impact of those shocks on jet angular velocity in YSOs is probably weak, at least for jets achieving m² ≫ 1.

Our adiabatic solutions do not allow energy dissipation, and so the shocks should go on, one above the other, with each shock producing its own local accelerating circuit. This is showcased by the presence of a second subset of shocks in simulations K4, K5 (see Fig. 16), and SP (see Fig. 19). Thus, this setup does not allow the presence of a “true” asymptotic circuit that would extend up to infinity on the jet axis.

These simulations are highly dependent on the numerical setup. In order to capture the expected shocks, we used an HLLD solver, switching to a more diffusive HLL solver and MINMOD linear spatial reconstruction in regions of extremely low density and very high Alfvén speed. For the same reason, we used a higher resolution in θ around the axis to resolve the shocks. Still, we were only able to reach very long timescales thanks to a novel method that boosts the numerical integration (see Appendix A).

5.2. Comparison with other numerical works

In this section, we only compare our findings with previous 2D platform simulations of nonrelativistic jets. Indeed, relativistic jets develop an electric force that deeply affects the asymptotic collimation, forbidding a direct comparison with our nonrelativistic setup. We also disregard 3D jet simulations as they usually introduce a whole new phenomenology related to jet instabilities that are not present in our work. Making 3D simulations of our jets is planned for future work.

Platform simulations of jets have a great many degrees of freedom and it is therefore very difficult to determine the exact generic results on jet collimation that can be derived from them. Even if the whole injection domain is chosen to be sub-SM, three free distributions must be chosen at the boundary (assumed to be the disk surface), usually B_z(r),ρ(r), and u_z(r). The disk being assumed to be Keplerian, most previous works used field lines rotating at Keplerian speeds and u_z = V_injV_K, where V_inj is a small dimensionless number. We use the following notations B_z ∝ r^α − 2 and ρ ∝ r^−α_ρ leading to

$$\begin{array}{c}\hfill \kappa =\frac{V_{\mathrm{inj}}}{{\mu }^2}={\kappa }_d{r}^{{\alpha }_{\kappa }},\mu =\frac{V_{\mathit{Az}}}{V_K}={\mu }_d{r}^{{\alpha }_{\mu }},\end{array}$$(20)

with α_κ = −2α_μ = −2α + α_ρ + 3 and the cylindrical radius is normalized to the inner radius R_d. For a given magnetic field distribution, the way the mass is injected in the outflow (or how the magnetic energy must vary within the disk) is an important quantity allowing to compare the various jet models. Our injection conditions (Eq. (8)) have α_ρ = 2α − 3, and therefore α_μ = α_κ = 0, which is in agreement with self-similar studies (BP82; F97). For our explored range in α < 1, our jets are always dominated by the mass flux emitted from the innermost disk regions. As discussed previously, the spatial and temporal scales are also very important in order to obtain recollimation shocks. We recall that our spherical domain goes up to R_ext = 5650R_d covered with 266 × 1408 zones and lasts T_f ≳ 10⁵ T_d.

The box size and timescales achieved in the pioneering works of Ouyed & Pudritz (1997a,b, 1999) were of course quite small, with a cylindrical domain (z, r) = (80, 20) in R_d units and a resolution of (500, 200) cells, with the simulations lasting up to T_f ∼ 500. These authors studied mostly α = 1 and α = 2 magnetic configurations, assuming α_ρ = 3/2, μ_d = 0.01 and with no injected spine. Their jets have therefore a steeply decreasing κ (or increasing μ) with α_κ = 3/2 − 2α, providing situations very different from ours. The authors argued that the nature of the outflow (steady or not) is mostly determined by the mass load κ, with unsteady jets containing shocks and associated knots arising at small values of κ_d ∼ 10⁻². While these shocks are indeed due to jet material being focused toward the axis, Ouyed & Pudritz (1997a,b, 1999) did not report any steady-state situation. Our own simulations show that their timescales were still too short to warrant a transverse jet balance, especially for α ≥ 1. Moreover, it remains unclear as to whether these knots were indeed a consequence of a small mass load κ_d or due to the boundary conditions used at the jet basis, which were too numerous and therefore over-determined the outflow dynamics (see discussions in e.g., Bogovalov 1997; Krasnopolsky et al. 1999; Ramsey & Clarke 2019).

Ustyugova et al. (1999) showed that if the simulation region is elongated in the z-direction, then Mach cones may be partially directed inside the domain, leading to an artificial influence (usually collimation) on the flow. Using a domain (z, r) = (200, 170) with 100 × 100 cells, these authors showed that this effect can be reduced with a square or spherical grid.

Pudritz et al. (2006) extended their work by exploring a larger range in α = 1, 3/4, 1/2, 1/4, using κ = 5r^{3/2 − 2α}, μ_d = 0.01 and no spine. These latter authors argued that the collimation of a jet depends on its radial current distribution, which in turn is prescribed by the mass load. Simulations with α = 1, 3/4 would collimate to cylinders due to a decreasing κ leading to a large B_ϕ, whereas simulations with α = 1/2, 1/4 with an increasing κ would produce a smaller B_ϕ and jets closer to wide-range outflows with parabolic collimation. However, our simulations show that the physical scales needed to observe the correct asymptotic state are much larger than those achieved in these early simulations. Moreover, it is indeed correct that self-collimation depends on B_ϕ, which, in a magnetic jet that carries away the disk angular momentum, namely −rB_ϕ/η ∝ Ωr², varies as B_ϕ ∝ κ(r)B_z at the disk surface. In their case, this expression leads to B_ϕ ∝ r^{−1/2 − α}, which is indeed more steeply decreasing with α. (However, our guess is that the collimation observed within the box of Pudritz et al. 2006 is mostly a consequence of the potential magnetic field configuration used as the initial condition, as also illustrated in Fig. 8.) Therefore, the smaller α, the wider (less collimated at a fixed distance) the jet.

The influence of the magnetic field profile α on the asymptotic jet collimation has also been investigated: Fendt (2006) performed 40 simulations in a larger cylindrical grid (z, r) = (300, 150) with 256 × 256 cells, with the simulations lasting up to T_f ∼ 300 to 5 10³ (for those achieving a steady state over at least 50% of the grid). He explored a wide range in α_ρ from 0.3 to 2 and in α from 0.5 to 1.8, using the same boundary conditions as Ouyed & Pudritz (1997b), with V_inj = 10⁻³, κ_d = 5, no spine, and μ_d varying between 0.1 and 2.67. Fendt (2006) confirmed that the degree of collimation is decreasing for a decreasing α regardless of α_ρ, in agreement with our suspicion that the overall MHD collimation trivially follows the potential field configuration (see also Sect. 3.2). For α > 1.6, no steady-state jet is actually found, with a wavy radial pattern evolving along the outflow. This is consistent with our finding that for α ≥ 1 the timescales for reaching stationarity become overwhelmingly long and also with the existence of radially oscillating, sub-FM analytical solutions for α ≥ 1 (Contopoulos & Lovelace 1994). Fendt (2006) also reports a degree of jet collimation increasing with the jet magnetization exponent, namely with μ_σ = −α_κ − 3/2 (see his Eq. 10). Now, of the 40 simulations, only 6 have α_κ = α_μ = 0 and 8 have α_κ > 0, meaning that most simulations describe a mass loading decreasing with the radius. None of the simulations show standing recollimation shocks, even in the BP82 case obtained with μ_d = 0.177 (we use μ_d = 1). Putting aside this difference, our Fig. 12 shows that around their final time T_f = 5000, we observed shocks only around Z ∼ 1900, far beyond their computational domain.

Krasnopolsky et al. (1999) used a cylindrical grid (z, r) = (80, 40) with 256 × 128 cells, with the simulations lasting up to T_f ∼ 170, introducing a ballistic axial flow below R_d (the spine), injected close to the escape speed and surrounded by a disk wind. These authors used the correct number of boundary conditions and, by testing the effects of adjusting the size of the box, they showed the drastic importance of the amount of magnetic flux becoming super-A within the box on the overall flow collimation. They studied mostly α = 1/2 and 3/4 with μ_d = 4 and rather flat density distributions leading to α_κ > 0, from 2 to 3/2. The authors do not report any time-dependent behavior seen in previous studies, which they attribute to both the existence of their sub-FM inner spine (where magneto-centrifugal acceleration is inefficient) and the correct treatment of boundary conditions. This latter work was extended by Krasnopolsky et al. (2003) on a much larger box (z, r) = (10³, 10³) with 190 × 210 zones, with the simulations lasting an unspecified time T_f. They only studied the case α = 1/2, with ejection from a finite zone r_o = R_d and r_o = 10R_d, yielding α_ρ = 1 (α_κ = 1) or α_ρ = 3 (α_κ = −1). The authors found that the collimation degree of this finite jet is improved for a steeper density profile, namely with a decreasing mass load with the radius, as discussed above. Krasnopolsky et al. (2003) report neither recollimation toward the jet axis nor radial oscillations, and attributed this behavior to their nonself-similar scaling. Our own results show instead that recollimation should be seen farther out (beyond their box) and that radial oscillations are expected only for α > 1.

Using the same grid and numerical setup as Krasnopolsky et al. (2003), Anderson et al. (2005) studied the effect of κ_d on the collimation of a cold BP82 jet model with α = 3/4 and α_ρ = 3/2 (thus α_κ = 0). These authors varied κ_d from 6.3 10⁻³ to 19 assuming that ejection takes place only from r_o = R_d and r_o = 10 R_d (but enforcing B_z to zero at the edge of the launching region), while we assumed ejection from the whole disk and varied κ_d only from 5 × 10⁻² to 1. Despite the truncation due to the limited ejection range and the (almost) purely radial magnetic field at the edge of the launching region, Anderson et al. (2005) recover the same results as in steady-state jet theory (FP97): jets become increasingly open as κ_d decreases (see discussion in Sect. 4.1). Anderson et al. (2005) do not report any recollimation shock (although wiggles can be seen in their Fig. 4) but again, our shocks fall below Z = 1000 (within their box) only for κ_d ∼ 1 (see Fig. 17). We conclude that their box was too small to observe any standing recollimation shock. The authors report the inability to reach steady state (the timescale T_f is unspecified) for κ_d larger than unity, when field lines start to oscillate and produce ripples that propagate outward. This behavior is consistent with analytical studies and is related to the capability to produce super-A flows when they are heavily loaded (or have a weak magnetic field). Indeed, magnetically driven cold flows are possible only up to κ ∼ 1, leading to a magnetic lever arm λ ∼ 2. For larger mass loads (and smaller λ), gravity plays an important role, with the Alfvén surface getting closer to the disk, requiring the field lines to be bent by much more than the fiducial 30° at the disk surface⁴ (see Fig. 4 and discussion around the Grad-Shafranov equation in Jacquemin-Ide et al. 2019).

The largest axisymmetric simulations have been provided by Ramsey & Clarke (2011, 2019), using nine levels of AMR in a cylindrical grid (z, r) = (8 × 10⁴, 5 × 10³)R_d with simulations lasting up to T_f ∼ 6 × 10⁴. These authors computed the propagation and evolution of eight jets up to observable scales, defined with varying mass loads κ_d from 5 × 10⁻² to 32 and α_ρ = 3/2, α = 1 (thus a decreasing mass load with α_κ = −1/2). In the simulations of these latter authors, mass is injected with V_inj = 10⁻³ and there is no injected spine as in Ouyed & Pudritz (1997a), although a spine naturally emerges. In all simulations, Ramsey & Clarke (2011, 2019) observe that regions beyond r_o ∼ 10 R_d fail to displace the hot atmosphere and that the outflow is stifled, despite the decrease in κ. This is actually consistent with our previous discussion for simulations with α ≥ 1, which take a much longer timescale to reach steady state. Nevertheless, as the inner parts of the outflow evolve on much shorter timescales, some quasi-stationary situation can settle (see their Sect. 5.3). With no surprise, this is the case for small mass loads, while knots appear for κ_d = 0.5 (simulation E) and are recurrent (quasi-periodic) for κ_d = 2. These knots are not to be compared with our standing recollimation shocks, as none of the MHD invariants are constant along field lines passing through them. They are made of plasmoids launched from R_d ≲ r_o ≲ 2R_d, where gas is both dense and hot. The knot formation mechanism is here directly related to the jet-launching process from this innermost disk region. Indeed, in this region, the field line bending is insufficient to drive the massive injected material, until a sufficiently strong toroidal field builds up and lifts the matter, in agreement with steady-state theory of massive outflows (F97; Jacquemin-Ide et al. 2019). The regularity of knot spacing is indicative of a simple oscillator related to the necessary build up of a strong toroidal field. These plasmoids are magnetically confined by the surrounding poloidal magnetic field, follow the path of the jet, and eventually merge together. For larger mass loads (κ_d = 8 and 32, simulations G and H), the outflows are fully unsteady while keeping their global structure (probably because of their 2D nature, as destroying instabilities such as kink or Kelvin-Helmholtz require 3D, as argued by the authors).

To our knowledge, no previous jet simulation has shown the existence of standing recollimation shocks, either because the computational domain was too small and/or the simulation timescales were too short. These limitations are even worse of course for simulations that do take into account the disk physics, as they must also struggle to follow the disk and the mass-loading process.

The first of these simulations computed an accretion–ejection configuration with α = 3/4 and α_ρ = 3/2 (the BP82 case) within a cylindrical grid (z, r) = (80, 40) on a time T_f = 251 only (Casse & Keppens 2002, 2004). On these timescales, the mass-loading process is computed, leading to the inside-out establishment of self-similar conditions with α_κ = 0. Further simulations, carried out with the same initial configuration but exploring various disk parameters, were computed on slightly extended scales, a grid (z, r) = (120, 40) on a timescale of T_f = 400 (Zanni et al. 2007; Tzeferacos et al. 2009, 2013) and a grid (z, r) = (180, 50) on a timescale of T_f = 5.6 × 10³ (Sheikhnezami et al. 2012). As most of these works were focused on the disk physics and less on the jet dynamics, they provided little information about the latter. The simulations of Stepanovs & Fendt (2016) were done on a spherical grid up to R_ext = 1500 with (N_R × N_θ) = (600 × 128) zones and up to T_f = 10⁴, for the same BP82 initial configuration. Such scales would be relevant for the appearance of recollimation shocks but they only show close-up views below R = 30 and focus instead on the accrection–ejection correlations. However, the long timescales allow us to see a radial redistribution of both the vertical magnetic field and the disk density (as both evolve on accretion timescales; Jacquemin-Ide et al. 2019), thereby modifying the initial strict self-similar conditions.

The time evolution of the disk magnetic field distribution has been reported previously (Murphy et al. 2009). These simulations were done in a cylindrical grid (z, r) = (120, 40) up to a time T_f ≃ 6 × 10³, and using α = 1/4 with α_ρ = 3/2. Such an initial magnetic field distribution leads to a magnetic energy density on the disk midplane that decreases very rapidly (∝r⁻¹), meaning that a super-FM ejection (with proper MHD invariants) only takes place up to a certain radius r_o ∼ 5 (Murphy et al. 2010). This latter study focused on this ejection from a limited zone within the disk and little was mentioned about the jets. However, we report that on the long timescale of the simulation, the magnetic field is seen to slowly evolve within the disk, leading to some readjustments also in the jet. How such a modification affects the jet transverse balance and possible standing recollimation shocks is an open issue that deserves further investigation.

We note that standing recollimation shocks have already been discussed in steady-state 2D jet simulations built upon analytical self-similar solutions. In these works, a cylindrical box is used, which starts at a z_o well above the disk (say z from z_o = 10 to 210 and r from 0 to 100 in units of R_d). This allows the whole domain to be filled with either only a self-similar BP82 jet model (Gracia et al. 2006; Stute et al. 2008) or a combination of an axial (meridionally self-similar) stellar wind surrounded by a BP82 jet model (Matsakos et al. 2008, 2009). The numerical procedure, which evolves the MHD equations over time for a set of boundary conditions, allows a stationary solution to be rapidly obtained on timescales of T_f ∼ 40 to 10³. A weak recollimation shock is always found between the axial flow and the BP82 jet, which fulfills most properties discussed in our paper. However, in strong contrast with our own work, the existence of this shock is unavoidable in these works and is directly imposed by the boundary conditions. Indeed, the outflow is already super-FM at the injection altitude z_o for all radii below r_o ∼ 6 (see for instance Fig. 1 in Matsakos et al. 2008), while field lines are already being focused toward the axis.

5.3. Astrophysical consequences

In this paper, we showcase one mechanism enabling the creation of a recollimating jet and its subsequent shocks. There are other models explaining the creation of such shocks. They could be triggered for instance by a sudden mismatch between the jet and the ambient medium pressure. Studying FRII jets such as those from the radio galaxy Cygnus A, Komissarov & Falle (1998) proposed that the jet confinement and its consequential shocks are caused by the thermal pressure of an external cocoon. For the case of FRI jets, in Perucho & Martí (2007) the jet expands until it becomes under-pressured with respect to the ambient medium, and then recollimates and generates shocks, unless a turbulent mixing layer at its interface with the ambient medium forbids its formation (Perucho 2020). In any case, such shocks happen much farther away than in our case and depend critically on the ambient pressure distribution.

On the contrary, the jets in our simulations are intrinsically collimated by the self-induced hoop stress (see Fig. 10). As shown in FP97 for self-similar cold models and proven here in full 2D time-dependent simulations, this force will lead the cold jets toward the axis, leading to the formation of standing recollimation shocks. Such a mechanism should therefore apply regardless of the external medium and around various astrophysical objects.

Extragalactic jets imaged by VLBI display knots of enhanced emission that could be associated with shocks (as they play an important role for the production of nonthermal emission). While most of these features are moving, some of them appear stationary (Lister et al. 2009, 2013; Doi et al. 2018 and Boccardi et al. 2017 for a review). The closely studied M87 jet is a particularly interesting case. It contains several moving and stationary bright features near the HST-1 complex (Asada & Nakamura 2012; Walker et al. 2018; Park et al. 2019), whose origin may be due to pressure imbalance when the jet reaches the Bondi radius. This distance is actually larger than the scales reached by our simulations. However, these are Newtonian and it is unclear whether or not relativistic effects (in particular the decollimating force due to the electric field) would push the recollimation scale farther out. In any case, we note that our nonrelativistic simulations provide shocks that are located on the same scale as the closest features in the M87 jet (see Fig. 2 of Asada & Nakamura 2012).

Protostellar jets also present some interesting features along the flow usually interpreted as being bow shocks, as in HH212 (Lee et al. 2017) or HH30 (Louvet et al. 2018). Their origin remains highly debated, either instabilities triggered during jet propagation or variability induced by a time-dependent jet production mechanism (as advocated for instance in HH212 by the remarkable jet–counter-jet symmetry; see Tabone et al. 2018). However, we suspect that whenever a jet undergoes an MHD recollimation shock that refracts the jet away from the axis, more shocks are to be expected downstream (and probably affected by the external pressure distribution). MHD recollimation may therefore provide an intrinsic means to trigger jet variability on observable scales. Stationary emission features are sometimes indeed detected, as in HH154 (Bonito et al. 2011). These features are located from a few tens to a few hundreds of astronomical units from the source, a distance comparable to the altitude of the first standing recollimation MHD shock. This is worthy of further investigation.

6. Conclusion

We present axisymmetric simulations of nonrelativistic MHD jets launched from a Keplerian platform. These are the first to show the formation of standing recollimation shocks, at large distances from the source. These recollimation shocks are intrinsic to the MHD collimation process and have been proposed as a natural outcome of self-similar jet-launching conditions (F97; Polko et al. 2010). Because they were never seen in previous MHD simulations of jets, the suspicion grew that recollimation would be a bias due to the self-similar ansatz. It turns out that the physical scales required to capture these shocks are much larger than those used in previous works. Using unprecedentedly large space and temporal scales allowed us to firmly demonstrate the existence of such internal standing shocks and thereby bridge the gap between analytical and numerical approaches.

We analyzed the conditions of formation of these recollimation shocks and show that they qualitatively follow the behavior demonstrated in analytical studies, namely that they get closer to the source as the mass load increases. We also confirm that the magnetic field distribution in the disk (B_z ∝ r^α − 2) is the key quantity shaping the asymptotic jet collimation. For our self-similar ejection setup, this MHD collimation closely follows the trend satisfied by the potential field: the larger the α the stronger the collimation. However, no steady-state solution is obtained for α ≥ 1, because of the difficulty in establishing a stationary self-consistent poloidal electric circuit along the outer jet regions. As the magnetic field distribution is very likely to evolve on the accretion time scale, we expect jet signatures to vary as well (see e.g., discussion in Barnier et al. 2022).

Despite their qualitative agreement with analytical studies, our results reveal an undeniable impact of the central axial flow on the jet asymptotics. This inner spine is not related to the Keplerian disk but instead to the central object and its interaction with the surrounding disk. Indeed, the spine carries a poloidal electric current responsible for the innermost jet collimation. However, it may also introduce extra localized spine–jet interactions, leading potentially to disruptive instabilities (like kink and/or Kelvin-Helmholtz) or, on the contrary, to global jet stabilization in 3D. Going to 3D is therefore necessary in order to assess the role of the inner spine and the possible persistence of recollimation shocks. In any case, our results confirm the role of the central object in shaping, through its spine, the collimation properties of the jets emitted by the surrounding disk. This is a very interesting topic that merits further investigation.

These internal recollimation shocks introduce several interesting features: (i) an enhanced emission likely seen as stationary knots in astrophysical jets; (ii) a sudden decrease in the rotation rate of the ejected material, and (iii) a possible electric decoupling between the pre-shock and the post-shock regions. This is of especially great interest as these shocks occur at observable distances, typically ∼150 − 200 au in the case of a YSO. However, our setup also assumes ejection up to several hundreds of astronomical units, which is clearly inconsistent with derived jet kinematics (see e.g., Ferreira et al. 2006; Tabone et al. 2020 and references therein). Simulations with ejection from only a finite zone within the disk (the JED) must therefore be carried out in order to verify whether MHD recollimation shocks are indeed maintained. This is a work in progress.

Acknowledgments

We thank the referee for providing thoughtful comments on the manuscript. The authors acknowledge financial support from the CNES French space agency and PNHE program of French CNRS. All the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.

References

Appendix A: Evolution on long timescales

In numerical simulations, the time increment is fixed by the–Friedrichs–Lewy condition Δt < Δx/C_max, where C_max is the maximal wave speed in the cell (in our case u + v_FM) and Δx is the cell size. In a standard simulation, the time increment of all cells is chosen by taking the absolute minimum of all the time increments in the full computational domain. In our simulations, this time step is set by the smallest cells around the inner spherical boundary at R_d, which also happen to have the strongest field and the highest Alfvén speed. However, these cells near the source are also the ones that converge the fastest to a stationary solution. Thus, the cells that we consider are converging to a steady state (meaning the relative variation of the density in one integration step is smaller than an arbitrarily small parameter) and are not used to determine the time increment of the cells that are still evolving in time. The time increment used to evolve the evolving cells is computed by taking a minimum over only the cells that have not converged yet. The cells that have converged to a steady state are still integrated in time using their own local time increment so as to ensure the stability of the integration and to be able to capture any perturbation that could possibly alter their steady condition. As the cells that converge the fastest are those characterized by the shortest time increment, the time increments used to evolve the cells that are still evolving in time and have not converged yet to a steady state become larger and larger.

It is important to point out that the stationary solutions obtained with this time boost are also a solution of the standard nonaccelerated algorithm.

Figure A.1 shows the gain in computing time obtained thanks to the time boost. The acceleration factor, defined as the ratio between the physical time reached using the time boost and the physical time that would have been achieved using the standard CFL condition, is plotted versus the progressive numbering of the outputted files.

Fig. A.1. Evolution of the acceleration for the simulations K2, K5, and A5.

Without the acceleration due to our handling of the CFL condition, the time interval between two outputs would have been constant, and the physical time of the solution would have been proportional to the output number. Any increase in the acceleration factor means that another batch of cells has converged. This means the time increment of the cells that are still evolving in time becomes larger, thus increasing the timescales reached.

This increase is clearly visible after the 300th output. Using the output number as a proxy for the computational cost of a simulation, this figure clearly shows that, at the end of the integration, the time boost enables us to reach timescales at least two or three orders of magnitude larger than using a standard CFL condition, without increasing the computational cost of the simulation. Analogously, without employing the time boost, we would have required two or three orders of magnitude more CPU hours to reach the same timescales. Our approach enabled us to produce simulations that would have consumed much more computing time otherwise. The reference simulation K2 consumed 725 CPU hours, but without the time boost it would have required almost two million CPU hours. This enabled us to work on simply 64 processors kindly provided by GRICAD (Grenoble Alpes Recherche - Infrastructure de Calcul Intensif et de Données).

For the evolution of the acceleration factor with the mass load, we can see in Table 1 that the simulations with κ closest to that of Blandford & Payne (1982) converge the fastest, reaching larger timescales at the end of the computation. Nevertheless, all seem to reach comparable convergence speeds: in Figure A.1 we see that the simulation K5 reaches an acceleration factor similar to that of K2 at the final output. For the evolution with the magnetic field, we can clearly see that the higher the α and thus the flatter the profile of the vertical magnetic field, the slower the simulation converges. For higher values of α, the jet is initially more collimated as B_r/B_z is higher on the disk. Therefore, the field lines further on the disk have a higher impact, retarding the global convergence. That is why the simulation A5 with α = 15/16 has not yet converged. As an instability develops, some cells that were previously stable become unstable, hence the decrease in acceleration.

Appendix B: Boundary conditions

The plots represented in Figure B.1 illustrate the injection boundary conditions we have chosen for the reference simulation K2. The Bernouillli invariant E, the poloidal magnetic field B_Φ, the vertical magnetic field B_z, the mass to magnetic flux ratio η = μ₀ρv_p/B_p, the rotation speed of the magnetic surfaces Ω_⋆ = Ω − ηB_Φ/(μ₀ρr), the speed of sound $C_s=\sqrt{\gamma P/\rho ,}$ and vertical Alfvénic speed $V_{A_z}=B_z/\sqrt{{\mu }_0\rho }$ are plotted on the first cell over the injection boundary. The toroidal magnetic field goes to zero on the axis for symmetry reasons, and $|B_{\mathrm{\Phi }}/B_z|\gtrsim 1$ on the disk : The JED magnetic field is weakly toroidal. The launching conditions are very cold as V_Az/C_s ∼ 10² on both the source and the disk.

Fig. B.1. Conditions on the lowermost cells of the reference simulation K2. Here all parameters are traced along the first cell above the lower boundary : R = 1 and θ ∈ [0; π/2] on the source and then θ = π/2 and R ∈ [0; Rext] on the disk. Shown are the Bernouilli invariant E, the toroidal and vertical magnetic field BΦ and Bz, the mass to magnetic flux ratio η, the rotation speed of magnetic surfaces Ω*, the speed of sound Cs and the vertical Alfvénic speed VAz, over the magnetic flux Ψ. The black vertical line corresponds to the flux anchored at (R = 1,θ = π/2), at the source/disk interface.

The reader may observe the power-law dependency with the magnetic flux on the whole disk (Ψ > 10, after the black vertical line) for all parameters but the electric current. This is directly induced by the self-similar ansatz. However, the torroidal magnetic field B_Φ breaks the power-law dependency and shows a swift decrease for Ψ > 2.10⁶. Of all eight variables, only B_R and B_Φ are free at the injection boundary (see Equation (8)), and need to cross a characteristic surface to be fixed. For the toroidal current, it is the Alfvénic surface. As all magnetic surfaces over Ψ ≳ 10⁹ never cross the Alfvénic surface, the current can never be fixed. Thus, the simulation cannot ever be stationary in this region.

As the disk is a self-similar jet-emitting disk, all dimensionless parameters are assumed to be independent of the radius (Blandford & Payne 1982). All these parameters are regrouped in section 2.4.

As explained in section 2.3.2.2, we defined a spline function f(θ) equal to zero on axis (θ = 0) and one at the inner disk radius (θ = π/2) to smoothly connect the axis values with the inner disk ones: f(θ)≡(3sin²θ − 2sin³θ)^3/2.

The injection speed (V_R on the source, −V_θ on the disk) is fixed at κμ² = 0.1 even when varying κ. However, even though the injection Mach number M_s = u_p/C_s is assumed to be the same all along the boundary, in order to account for the varying inclination of the magnetic surfaces, its value is modified with the parameter α from simulation to simulation. Its variation with α is shown in Table 1.

Appendix C: Rankine-Hugoniot jump conditions

In this section, we write the Rankine-Hugoniot jump conditions valid for standing, adiabatic recollimation shocks. Contrary to Ouyed & Pudritz (1993), we take into account the toroidal magnetic field as the shocks arise when that component is dominant. The local jump [A]=A₂ − A₁ between a pre-shock quantity A₁ and its post-shock value A₂ are expressed in the rest frame as

$$\begin{array}{cc}\hfill [\rho u_{\perp }]=0& \hfill \\ \hfill [\rho u_{\perp }(\frac{u^2}{2}+H)+\frac{B^2}{{\mu }_o}{u}_{\perp }-\frac{\mathit{u}·\mathit{B}}{{\mu }_o}{B}_{\perp }]=0& \hfill \\ \hfill [P+\rho u_{\perp }^2+\frac{B_{\Vert }^2-B_{\perp }^2}{2{\mu }_o}]=0& \hfill \\ \hfill [\rho u_{\perp }{\mathit{u}}_{\Vert }-\frac{B_{\perp }}{{\mu }_o}{\mathit{B}}_{\Vert }]=0& \\ \hfill [B_{\perp }]=0& \hfill \\ \hfill [B_{\perp }{\mathit{u}}_{\Vert }-u_{\perp }{\mathit{B}}_{\Vert }]=0,& \hfill \end{array}$$(C.1)

where $H=C_s^2/(\mathrm{\Gamma }-1)$ is the enthalpy and u_⊥, B_⊥ (respectively u_∥, B_∥ ) are the normal (respectively tangential) components to the shock surface. The shock is axisymmetric, and so the tangential component of the magnetic field B_‖ = B_te_t + B_ϕe_ϕ, whereas the poloidal component is B_p = B_te_t + B_⊥e_⊥, with the unit vectors (e_⊥,e_t,e_ϕ) defining a local orthonormal coordinate system. As these jump conditions express the conservation of mass, angular momentum, and energy in ideal MHD, the five MHD invariants along a given magnetic surface (η, Ω_*, L, E, S) are therefore also conserved (see Fig.7).

In the case of a shock, the mass flux through the surface is nonzero, which requires B_⊥2 = B_⊥1 ≠ 0 and leads to

$$\begin{array}{c}\hfill {\mathit{B}}_{\Vert 2}=\frac{m^2-1}{m^2/\chi -1}{\mathit{B}}_{\Vert 1},\end{array}$$(C.2)

where m = u_⊥/V_Ap, ⊥ = u_p/V_Ap is the Alfvénic Mach number of the incoming (pre-shock) flow and χ = ρ₂/ρ₁ = u_⊥1/u_⊥2 is the shock compression rate. This equation shows that there are three nontrivial discontinuities with χ ≥ 1: (1) an oblique shock with m² > χ > 1, (2) a normal shock with m² = χ > 1 (requiring B_‖1 = 0), and (3) an Alfvén shear discontinuity with m² = χ = 1 (allowing an arbitrary jump between the two tangential field components). The oblique shock is the only case studied here.

After some algebra all post-shock quantities can be expressed as function of the pre-shock ones, in particular

$$\begin{array}{cc}& \frac{B_{\varphi 2}}{B_{\varphi 1}}=\frac{B_{t2}}{B_{t1}}=\chi \frac{m^2-1}{m^2-\chi }\hfill \\ \hfill & \frac{u_{\varphi 2}}{u_{\varphi 1}}=\frac{m^2-1}{m^2-\chi }\frac{m^2{r}_A^2-\chi r^2}{m^2{r}_A^2-r^2}\hfill \\ \hfill & \frac{P_2}{P_1}=1+\mathrm{\Gamma }{m}_s^2(\chi -1)(\frac{1}{\chi }+\frac{b^2}{2}\frac{2\chi -m^2(1+\chi )}{{(\chi -m^2)}^2})\hfill \\ \hfill & \frac{T_2}{T_1}=\frac{1}{\chi }\frac{P_2}{P_1},\hfill \end{array}$$(C.3)

where the sonic Mach number m_s = u_⊥/C_s and magnetic shear b² = (B_∥/B_⊥)² are computed in the pre-shock region. Of particular interest are the relative variations of the toroidal magnetic field component δB_ϕ = B_ϕ2/B_ϕ1 − 1 and the plasma angular velocity δΩ = Ω₂/Ω₁ − 1, as well as the total deflection angle of the poloidal magnetic surface δi = i₂ − i₁ where tani = B_t/B_⊥, which read

$$\begin{array}{cc}\hfill \delta B_{\varphi }& =(\chi -1)\frac{m^2}{m^2-\chi }\hfill \\ \hfill -\delta \mathrm{\Omega }& =\frac{\chi -1}{m^2-\chi }\frac{m^2(r^2-r_A^2)}{m^2{r}_A^2-r^2}\le \frac{\chi -1}{m^2-\chi }\hfill \\ \hfill tan\delta i& =\frac{m^2(\chi -1)}{m^2-\chi }\frac{tan{i}_1}{1+\chi {tan}^2{i}_1\frac{m^2-1}{m^2-\chi }}.\hfill \end{array}$$(C.4)

These quantities are plotted in Fig. 9. The compression rate χ is the solution of the cubic polynomial equation

$$\begin{array}{c}\hfill -A{\chi }^3+B{\chi }^2-C\chi +D=0\end{array}$$(C.5)

with

$$\begin{array}{cc}\hfill A& =1+b^2+\frac{1+{\chi }_o}{\mathrm{\Gamma }{m}_s^2}\hfill \\ \hfill B& ={\chi }_o(1+b^2)+2m^2(1+\frac{1+{\chi }_o}{\mathrm{\Gamma }{m}_s^2}-b^2\frac{{\chi }_o-3}{4})\hfill \\ \hfill C& =m^2(2{\chi }_o+b^2\frac{1+{\chi }_o}{2}+m^2(1+\frac{1+{\chi }_o}{\mathrm{\Gamma }{m}_s^2}))\hfill \\ \hfill D& ={\chi }_o{m}^4,\hfill \end{array}$$

where χ_o = (Γ + 1)/(Γ − 1) is the maximal compression ratio for a hydrodynamic shock. Equation C.5 has one positive root only for an incoming super-FM flow, namely for n_⊥ = u_⊥/V_fm, ⊥ larger than unity.

We are dealing here with supersonic (m_s > > 1) and super-A (m > > 1) cold jets, where the dominant magnetic field is the toroidal one (b² ≃ (B_ϕ/B_⊥)² > > 1). The FM Mach number in the normal direction therefore writes n_⊥ ≃ mV_Ap, ⊥/V_Aϕ = m/b, which leads to the simplified equation for χ

$$\begin{array}{c}\hfill \frac{{\chi }_o-3}{2}{\chi }^2+(\frac{1+{\chi }_o}{2}+n_{\perp }^2)\chi -{\chi }_o{n}_{\perp }^2=0.\end{array}$$(C.6)

This shows that whenever jets reach a very large FM Mach number n_⊥, a large compression rate χ ≃ χ_o is possible. But this is never achieved in our case. Indeed, the poloidal FM Mach number n = u_p/V_FM, p (> n_⊥) writes

$$\begin{array}{c}\hfill n^2={\omega }_A\frac{B_{\mathit{pA}}{r}_A^2}{B_p{r}^2}\frac{1-1/m^2}{1-r_A^2/r^2}{\left(\frac{u_p}{{\mathrm{\Omega }}_{\ast }{r}_A}\right)}^3\sim {\omega }_A{\left(\frac{u_p}{{\mathrm{\Omega }}_{\ast }{r}_A}\right)}^3,\end{array}$$(C.7)

where ω_A = Ω_*r_A/V_Ap, A is the fastness parameter introduced in F97 (ratio at the Alfvén point of the speed of the magnetic rotator to the poloidal Alfvén speed). For magneto-centrifugal jets like ours, with m² > > 1, r > > r_A and achieving their maximal velocity $u_p\sim \sqrt{2}{\mathrm{\Omega }}_{\ast }{r}_A$, the FM Mach number is n² ∼ ω_A, which is larger than but of the order of unity (see also Krasnopolsky et al. 2003). As a consequence, we expect rather weak shocks as illustrated by the small values of χ achieved along the various shocks (see Fig. 9).

Figure C.1 plots the theoretical solution χ_th of Eq. C.5 (in yellow) computed along the main recollimation shock of our simulation K2 and compares it with the ratio χ = ρ₂/ρ₁ (in blue) directly measured (see Fig. 4). The correspondence is very good, with discrepancies remaining below a few percent. The two regions where larger differences are obtained correspond to the positions where the two smaller shocks (triggered at the spine-jet interface) merge with the main shock: the orange one near r ∼ 500 and the cyan one near r ∼ 900.

Fig. C.1. Distribution at tend along the main recollimation shock of the compression ratio for simulation K2. The yellow curve is the theoretical solution χth of Eq. C.5, computed using the pre-shock quantities, while the blue curve is the ratio χ = ρ2/ρ1.

The shocks were detected by following all magnetic field lines anchored on the disk and looking for discontinuities. This is not obvious in a discrete grid. To do so, we computed the derivative of the toroidal magnetic field (δB_ϕ) over the curvilinear abscissa along the field line, as shocks are best seen with the electric poloidal current and explored its local extrema. We checked that a different approach, based on the calculation of the refraction angle δi of the poloidal magnetic surface, produces very similar results. This gave us the shock locations used to produce the plots in Figures 2 and 17. As PLUTO has a shock capturing scheme, each shock is resolved and has a finite width. To determine the shock width, we checked that the density was growing within the shock as expected. Then, still following the field line, we looked for the closest local minimum and maximum in density. The positions of these extrema allowed us to compute the values of the pre-shock and the post-shock quantities, respectively. These were finally used to compute the parameters leading to the Figures 9 and C.1.

All Tables

All Figures

Fig. 1. Sketch of our computational domain. The central object and its interaction with the innermost disk are located below the inner boundary at Rd (green region), the near-Keplerian jet-emitting disk (JED) being established from Rd to the end of the domain Rext. An axial outflow (the spine) is emitted from the central regions (in red) and the jet is emitted from the JED (in blue). The solid purple line represents a recollimation shock surface starting on the axis at a height Zc. For each point N lying on this surface, we use local poloidal unit vectors (e⊥, e∥), respectively perpendicular and parallel to the shock surface. Also, at any point M inside the domain, we either use spherical (eR,eθ,eϕ) or cylindrical (er,eϕ,ez) coordinates.

In the text

Fig. 2. Snapshot at tend of our Blandford & Payne simulation K2. Left: global view with field lines anchored on the disk at ro = 3; 15; 40; 80; 160; 320; 600; 1000; 1500, where the background is the logarithm of the FM mach number n. Right: close-up view of the innermost regions, where the background is the logarithm of the density. In both panels, black solid lines are the poloidal magnetic surfaces, the yellow solid lines are isocontours of the poloidal electrical current, and the red dashed (resp. dotted) line is the FM (resp. Alfvén) critical surface.

In the text

	Fig. 3. Snapshot of our reference simulation K2 at tend. We use the same color coding as in Fig. 2, left. The black arrows show the poloidal velocity. The white lines are streamlines inside which (from left to right) 50%, 75%, and 100% of the super-FM (spine+jet) mass outflow rate is carried in.
In the text

	Fig. 4. Close-up view of our reference simulation K2 at tend showing the shock forming region, with field lines anchored at ro = 1.2; 2; 3; 4; 5; 7; 9. The five shocks are highlighted in red, orange, cyan, blue, and purple. We use the same color coding as in the left panel of Fig. 2.
In the text

Fig. 5. Late evolution of several global jet quantities for the simulation K2: the radius ro, FM of the last super-FM surface, the jet mass-loss rate Ṁjet, and the two colatitudes ${\theta }_A^{\text{ext}}$ and ${\theta }_{\mathrm{FM}}^{\mathrm{ext}}$ that define the position of the two critical surfaces. A slight drift from their initial value is indeed observed. The values provided in Table 1 are those achieved at the final time.

In the text

	Fig. 6. Evolution of the various energy contributions along a magnetic surface of anchoring radius r0 = 100 for the simulation K2 at tend: the Bernoulli invariant E, the gravitational potential ΦG, the total specific kinetic energy u2/2, the enthalpy H, and the magnetic energy −Ω*rBΦ/η. The absicssa is the altitude Z(Ψ).
In the text

	Fig. 7. Evolution of the MHD invariants along field lines of two different anchoring radii ro at tend for simulation K2. All invariants have been normalized to their values at ro. The absicssa is the altitude Z(Ψ).
In the text

	Fig. 8. Evolution of several magnetic field lines during the simulation computation, for different anchoring radii R0 and for the reference simulation K2. The field lines at the first output of the simulation (initial conditions) are shown in dotted lines. The field lines at the last output of the simulation (final state) are shown as full lines.
In the text

Fig. 9. Distributions along the shock of several quantities for K2 at tend. Left: normal incident FM Mach number n⊥, Alfvénic Mach number m, compression rate χ, relative variations of the toroidal magnetic field δBϕ and plasma angular velocity −δΩ, and total deviation δi (in rad) of the poloidal magnetic field line for the main recollimation shock. Right: compression rates χ of all shocks appearing in Fig. 4, using the same color code. The main shock corresponds to the red curve. See Appendix C for more details on this figure.

In the text

	Fig. 10. Radial distribution of the radial accelerations and their sum at the altitude Z = 2400 for the simulation K2 at tend. The vertical dashed line corresponds to the spine–jet interface, namely the field line anchored at ro = Rd.
In the text

Fig. 11. Plot of the poloidal electric circuits at tend for simulation K2. The two red curves are the critical surfaces, Alfvén (dotted) and FM (dashed). The yellow curves are the poloidal electric circuits, defined as isocontours of rBϕ, where the arrow indicate the direction of the poloidal current density Jp. Four circuits are highlighted in particular: (1) the envelope of the inner accelerating current in white (rBϕ = −2.06), (2) the outermost circuit still fully enclosed within the domain in blue (Bϕ = −2.005), (3) a circuit closed outside the domain in orange (rBϕ = −1.80), and (4) a post-shock accelerating circuit in purple (also with rBϕ = −2.06).

In the text

	Fig. 12. Altitude Z of the different shocks (measured at the axis) as a function of time (in Td units) for simulation K2. The vertical lines correspond to the six times ti used in Fig. 13: t1 = 551, t2 = 2.08 × 103, t3 = 8.51 × 103, t4 = 1.99 × 104, t5 = 1.05 × 105, and t6 = 1.58 × 105. The last vertical line is tend = 6.51 × 105.
In the text

Fig. 13. Snapshots of our reference simulation K2 at different times (given in Td units). From top to bottom, left to right: t1 = 551, t2 = 2.08 × 103, t3 = 8.51 × 103, t4 = 1.99 × 104, t5 = 1.05 × 105, and t6 = 1.58 × 105. The background color is the logarithm of the density, black lines are the magnetic surfaces, red lines the Alfvén (dotted) and FM (dashed) surfaces, and yellow curves are isocontours of the poloidal electric current.

In the text

	Fig. 14. Time evolution of the cylindrical radius r measured at Z = 3500 of the magnetic surface anchored at ro = 3 (blue curve) and the electric current I = rBϕ (red) flowing within that surface for simulation K2. The two vertical dashed lines correspond to t5 = 1.05 × 105 and t6 = 1.58 × 105.
In the text

Fig. 15. Jet parameter space λ(κ) at the final stage of our simulations K1–K5 with α = 3/4. Each simulation is obtained for a unique mass loading κ and gives rise to a distribution of the magnetic lever arm λ with the radius: green, red, and blue dots correspond to anchoring radii ro = 5, 50, and 500, respectively. The solid curves are obtained for constant values (indicated in the panel) of the initial magnetic field inclination ${\xi }_0^{\prime }=B_r/B_z$ at the disk surface.

In the text

Fig. 16. Influence of the mass-loading parameter κ on the final stage of jets obtained with α = 3/4. The color background is the logarithm of the FM Mach number n, black solid lines are field lines, yellow lines are isocontours of the electric current rBϕ and the red dashed (resp. dotted) curve is the FM (resp. Alfvén) critical surface. The left panels show the whole domain and the right panels a close-up view around the shock-formation regions. In the left panel, the field lines anchoring radii are ro = 3; 15; 40; 80; 160; 320; 600; 1000; 1500. In the right panel, the field lines anchoring radii are ro = 1.2; 2; 3; 4; 5; 7; 9; 11; 13; 15.

In the text

	Fig. 17. Influence of κ (orange lines) and α (blue lines) on the altitude of the main recollimation shock. This is done by measuring two altitudes for each shock: its height at the axis (Zshock, solid lines) and the altitude of its outer edge (Ztip, dashed lines). The scale for κ is indicated below, while the scale for α is above.
In the text

	Fig. 18. Influence of the magnetic field distribution α on the final stage of jets obtained with κ = 0.1. We use the same notations, colors, and field lines anchoring radii as in Fig. 13.
In the text

	Fig. 19. Snapshot at tend of our SP simulation with a rotating spine, α = 3/4 and κ = 0.1. We use the same color coding as in Fig. 2. The magnetic field lines (black solid lines) are anchored at the same disk radii.
In the text