© The Authors 2024

1. Introduction

Accretion disks are found in a wide variety of astrophysical objects: around supermassive black holes (BHs) at the center of active galaxies (AGNs) and quasars (Blandford et al. 2019; Lu et al. 2023), around every low-mass young stellar object (YSO) where planet formation occurs (Agra-Amboage et al. 2014; de Valon et al. 2020), around interacting binary systems hosting either a white dwarf (cataclysmic variables), a neutron star or a BH (X-ray binaries) (Ponti et al. 2016; Tudor et al. 2017), or even some post-AGB stars (Gorlova et al. 2015; Bollen et al. 2017). For many years, accretion was believed to be mainly driven by turbulence, which would allow an outward radial transport of the disk angular momentum, while mass would accrete onto the central object (Shakura & Sunyaev 1973). The quest for the instability that would allow nearly Keplerian accretion flows to become turbulent led to the conclusion that accretion disks must be magnetized (Balbus & Hawley 1991). The magneto-rotational instability (MRI) is indeed a very robust and fast (growing on dynamical timescales) instability that leads to fully developed 3D turbulence (Balbus 2003). It has thus become clear that accretion disks must be described as magnetized plasmas, namely within the framework of magnetohydrodynamics (MHD).

Simultaneously, highly collimated bipolar supersonic jets have also been discovered emerging from the same astrophysical objects. Mostly detected in the radio band around compact objects or in radio to optical emission lines in YSOs, they share the property of being emitted from the innermost regions close to the central object. Moreover, both compact sources and YSOs show a tight correlation between accretion and ejection signatures (Corbel et al. 2013; Nisini et al. 2018). These properties appear to be universal, as they do not depend on the nature of the central object (BH or stellar object), thus favoring the scenario where jets are emitted from the accretion disks. Since, in all resolved sources, jet collimation appears to be effective already very close to the source, jet acceleration and collimation must work hand in hand. This can occur naturally if a large-scale vertical magnetic field defines a magnetic channel along which the ejected plasma flows. Such a solution was first proposed by Lovelace (1976) and Blandford (1976), and the coupled acceleration and collimation dynamics of jets emitted from the surface of a Keplerian accretion disk was later demonstrated by Blandford & Payne (1982).

There have been many attempts to analytically link accretion disks to magnetized outflows, and it is beyond the scope of this introduction to cite them all. Most models either rely on strong approximations (e.g., neglected forces, negligible wind mass loss, and simplified wind equations, Wardle & Koenigl 1993; Narayan & Yi 1995; Li 1995; Ogilvie & Livio 1998; Bai et al. 2016 to cite only a few) or use a self-similar variable separation approach that allows all dynamical terms to be included in the equations (Ferreira & Pelletier 1993, 1995).

Only the anomalous turbulent coefficients need to be specified to close the system of equations describing the accretion-ejection theory. The turbulent transport coefficients resulting from MRI turbulence can be derived from measurements conducted in 3D shearing box simulations (Hawley et al. 1995; Salvesen et al. 2016) or 3D global simulations (Zhu & Stone 2018; Jacquemin-Ide et al. 2021). The Shakura & Sunyaev (1973) angular momentum transport coefficient has indeed been computed numerous times, and it is a well-constrained function of the local magnetic field strength (Salvesen et al. 2016). In contrast, turbulence resistivity has been measured sparingly, and no clear dependency on the field strength has been identified (Lesur & Longaretti 2009; Guan & Gammie 2009; Fromang & Stone 2009). The situation becomes even more obtuse when pertaining to the complete resistivity tensor (Gressel & Pessah 2015).

Magnetized accretion-ejection structures describe both the (turbulent) viscous and resistive MHD disk and its two (laminar) ideal MHD jets as a single, mathematically connected, system. The smooth crossing of the Alfvén point determines the position of this critical point as well as the strength of the toroidal field at the disk surface. This, in turn, fixes not only the torque allowing accretion but also the vertical magnetic compression acting on the disk, and hence the disk magnetization μ that is consistent with the disk ejection efficiency ξ and the smooth crossing of the slow-magnetosonic critical point (see Ferreira 1997 for details).

The goal of the theory is then to provide the relation ξ(μ), where the disk ejection efficiency is defined as Ṁ_a(r) ∝ r^ξ and Ṁ_a(r) is the disk accretion rate measured at a given radius. The theory aims to directly link the properties of the disk to the outflow. This is crucial for understanding the long-term evolution of an accretion disk and connecting observed outflow properties to the disk itself. Models of disk evolution, considering mass loss through outflows, must directly prescribe the ejection efficiency, which would otherwise remain unconstrained without an accretion-ejection theory (Tabone et al. 2022). Ejection efficiency is also vital for comparing theories of magnetized outflows to observations of disk winds in low-mass X-ray binaries (LMXRBs), leading to constraints on disk properties (Chakravorty et al. 2016, 2023; Ranjan Datta et al. 2024). Moreover, the mass ejection index significantly influences the history of X-ray binaries (XRBs) by altering mass and angular momentum loss, thus affecting the secular evolution of orbital separation (Gallegos-Garcia et al. 2024).

The first self-consistent accretion-ejection solutions were found for isothermal (Ferreira 1997, hereafter F97) or adiabatic (Casse & Ferreira 2000a) cold outflows. The disk magnetization μ was found to lie in a very narrow range, between 0.1 and 0.8, with a typical disk ejection efficiency ξ ∼ 0.01, leading to very fast and tenuous collimated jets. The accretion disk associated with these solutions has thus been termed jet emitting disk (JED, Ferreira et al. 2006). In a JED, accretion is supersonic and the disk is much less radiative than any other usual accretion disk solution (due to energy transfer to jets), with important implications for instance in X-ray binary cycles (Ferreira et al. 2006; Marcel et al. 2018, 2022 and references therein). We note that the presence of supersonic accretion was subsequently found in 3D global simulations of highly magnetized disks (Jacquemin-Ide et al. 2021; Scepi et al. 2024).

Until recently, the minimum value of the disk magnetization was constrained to order equipartition within the JED theory. At μ <  10⁻², MRI sets in despite the presence of anomalous diffusivities, and solutions display spatial (vertical) oscillations. Such oscillations are actually nonlinear channel modes due to the MRI (Jacquemin-Ide et al. 2019). Despite a smaller magnetic field, these oscillations help to drive the outflow by building up a stronger toroidal field at the disk surface. Since the vertical pinching of the disk is smaller, these solutions allow for denser isothermal (cold) outflows with a typical ejection efficiency ξ ∼ 0.1. In this case, the denser outflow is also much slower, so the accretion disk associated with these weakly magnetized solutions has been termed wind emitting disk (WED). In a WED, accretion is always subsonic and the magnetic configuration is reminiscent of the so-called magnetic towers (Lynden-Bell 2003). As discussed in Jacquemin-Ide et al. (2019), it is doubtful that oscillations would survive within the disk, since secondary instabilities (such as Kelvin-Helmholtz or Rayleigh-Taylor) would most likely be triggered, leading to enhanced anomalous transport within the turbulent region. This is therefore a good incentive to look at global 3D numerical solutions where MHD turbulence is fully resolved. Another reason to look at global simulations is that local shearing box simulations are subject to several biases when simulating stratified disks with vertical magnetic fields, which could skew the vertical profiles of the anomalous transport coefficients (Lesur et al. 2013).

Nowadays, general relativistic magnetohydrodynamic (GRMHD) simulations are commonly categorized using the magnetically arrested disk (MAD) and the standard and normal evolution (SANE) terminology (e.g., Event Horizon Telescope Collaboration 2021). A MAD is defined by the magnetic flux at the horizon of the BH¹. Specifically, it is considered that the system enters the MAD state when the flux normalized by the accretion rate reaches its maximal value (see e.g., Tchekhovskoy et al. 2011; Narayan et al. 2022 and references therein). On the contrary, whenever the magnetic flux on the BH is less than that specific value, the numerical outcome has been termed SANE. However, the numerical convergence to one state or another depends mostly on two things: (1) the initial conditions (how much initial magnetic flux is available in the computational domain) and (2) how long the simulation runs (since the flux appears to be advected inward, leading to a growth of the magnetization of the inner regions). The reason why the definitions of the MAD and SANE state rely on the dynamics of the inner BH and not on the disk structure is twofold: (1) for most BH spin, the Blandford & Znajek (1977) jet energetically outshines any contribution from the disk; and (2) due to numerical limitations, only the inner regions, r <  20 r_g, would have achieved a complete steady state, making a study of the large-scale turbulent disk dynamics out of reach. There is no explanation yet of the physical state reached by the accretion disk in each case. It is nevertheless tempting to associate numerical MADs and SANEs with their semi-analytical counterparts, JEDs and WEDs, due to them being associated with different field strengths.

It is only recently that the turbulent structure of a disk threaded by a large-scale vertical field has started to be thoroughly analyzed (Scepi et al. 2024; Manikantan et al. 2024), following similar efforts in nonrelativistic simulations (Zhu & Stone 2018; Jacquemin-Ide et al. 2021). Those global numerical investigations on the turbulent structure have shown that the role played by the turbulent magnetic pressure is of paramount importance in shaping the disk vertical balance (Jacquemin-Ide et al. 2021; Scepi et al. 2024). This is a major effect that has not yet been studied in semi-analytical approaches of magnetized accretion-ejection structures.

The aim of this paper is precisely to investigate the influence of turbulent magnetic pressure in highly magnetized accretion disks, namely in JEDs. Section 2 recalls the assumptions, equations, and parameters needed to fully describe JEDs. In particular, we discuss how MHD turbulence is taken into account in our self-similar steady-state mean-field approach. Section 3 focuses on analyzing the impact of this additional pressure term on the disk structure, by considering only the slow magnetosonic constraint. The full parameter space, including the additional Alfvénic constraint, is analyzed in Section 4. While Sections 3 and 4 focus on a fiducial set of parameters, in Section 5 we vary the disk thickness and the level of the magnetic diffusivities in order to explore their effects on the JED response to a turbulent magnetic pressure. Section 6 summarizes our results and discusses their implications, in particular providing comparisons with some published results of 3D numerical simulations. Finally, we conclude in Section 7.

2. Accretion-ejection theory

2.1. Magnetohydrodynamic equations

The accretion-ejection theory assumes axisymmetry and stationarity. In this framework, a disk of plasma orbits in near-Keplerian motion around a central object of mass M. The disk self-gravity is neglected and the gravitational potential of the central object is Newtonian. The disk is assumed to be threaded by a large-scale vertical magnetic field. The whole system is assumed to have a plane symmetry so that z = 0 describes the disk midplane and the bipolar jets are symmetric. Using the cylindrical coordinates (r, ϕ, z) the plasma velocity and magnetic field are decomposed into poloidal and toroidal components, namely u = u_p + Ωre_ϕ and B = B_p + B_ϕe_ϕ, where Ω is the plasma angular velocity. The poloidal magnetic field writes

$$\begin{array}{c}\hfill {\mathit{B}}_p=\frac{\mathrm{\nabla }a}{r}\times {\mathit{e}}_{\varphi },\end{array}$$(1)

where a(r, z) is the magnetic flux function (a(r, z)=a_o describes a magnetic surface of constant poloidal magnetic flux). Hence, a jet can be seen as magnetic surfaces anchored over a radial extent of the accretion disk and nested around each other.

The set of partial differential equations describing a steady-state magnetized accretion-ejection structure writes

$$\begin{array}{cc}& \mathbf{\nabla }·\left(\rho {\mathit{u}}_p\right)=0\hfill \end{array}$$(2)

$$\begin{array}{cc}& \rho (\mathit{u}·\mathbf{\nabla })\mathit{u}=-\mathbf{\nabla }(P+P_{\mathrm{turb}})+\rho \mathbf{\nabla }{\mathrm{\Phi }}_G+\mathit{J}\times \mathit{B}+\mathbf{\nabla }·\mathcal{T}\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\eta }_m{J}_{\varphi }{\mathit{e}}_{\varphi }={\mathit{u}}_p\times {\mathit{B}}_p\hfill \end{array}$$(4)

$$\begin{array}{cc}& \mathbf{\nabla }·\left(\frac{{\nu }_m^{\prime }}{r^2}\mathbf{\nabla }\left[r{B}_{\varphi }\right]\right)=\mathbf{\nabla }·\left(\frac{1}{r}[B_{\varphi }{\mathit{u}}_p-\mathrm{\Omega }r{\mathit{B}}_p]\right)\hfill \end{array}$$(5)

$$\begin{array}{cc}& P=\rho C_s^2,\hfill \end{array}$$(6)

and accounts for, respectively, mass and momentum conservation, Ohm’s law, toroidal magnetic field induction and the equation of state (Ferreira & Pelletier 1995; Casse & Ferreira 2000a). In these equations, ρ is the plasma density, P the total (gas + radiation) pressure, ${\mathrm{\Phi }}_G=-GM/\sqrt{r^2+z^2}$ the gravitational potential of the central object, J = ∇ × B/μ₀ the electric current density and C_s the sound speed. In this paper, radiation pressure is neglected and the plasma is assumed to be thermalized (identical electronic and ionic temperatures) so that C_s² = 2ρkT/m_p, where m_p is the proton mass. An energy equation should then be provided to compute the plasma temperature T. However, we will only focus on isothermal structures where T remains constant along magnetic field lines but is allowed to vary radially.

The disk is expected to be fully turbulent (Balbus 2003), so the above equations must be understood as a mean-field approach only. Such turbulence is further assumed to give rise to enhanced (anomalous) transport coefficients such as a viscosity ν_v, magnetic diffusivity in the poloidal plane ν_m (and resistivity η_m = μ_oν_m, where μ_o is the vacuum permeability) and in the toroidal direction ν′_m. The presence of an anomalous viscosity implies the existence of a viscous stress tensor 𝒯 (see next section for its expression). The novelty of this study is the inclusion of a turbulent magnetic pressure P_turb in Eq. (3). As we will show, this term indeed plays a major role and brings about significant changes in the overall accretion-ejection picture.

Magnetohydrodynamic turbulence is expected to decline vertically as the ejected plasma moves away from the disk. By definition, the flow transits from resistive MHD to ideal MHD at the disk surface, allowing for a smooth connexion to a super-magnetosonic (hereafter super-SM) and, further out, super-Alfvénic (hereafter super-A) outflow. Once in ideal MHD, the above equations lead to (F97)

$$\begin{array}{cc}& {\mathit{u}}_p=\frac{\eta (a)}{{\mu }_o\rho }{\mathit{B}}_p\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\mathrm{\Omega }}_{\ast }(a)=\mathrm{\Omega }-\eta \frac{B_{\varphi }}{{\mu }_0\rho r}=\frac{\mathrm{\Omega }}{(1-g)},\mathrm{with}g=\frac{m^2}{m^2-1}(1-\frac{r_A^2}{r^2})\hfill \end{array}$$(8)

$$\begin{array}{cc}& L(a)={\mathrm{\Omega }}_{\ast }{r}_A^2=\mathrm{\Omega }{r}^2-\frac{r{B}_{\varphi }}{\eta }\hfill \end{array}$$(9)

$$\begin{array}{cc}& E=\frac{u^2}{2}+H+{\mathrm{\Phi }}_G-{\mathrm{\Omega }}_{\ast }\frac{r{B}_{\varphi }}{\eta }+E_{\mathrm{turb}},\hfill \end{array}$$(10)

where η(a) is the mass flux to magnetic flux ratio, Ω_*(a) the rotation rate of the magnetic surface, L(a) the total specific angular momentum carried along that surface, m = u_p/V_Ap is the poloidal Alfvénic Mach number, r_A the cylindrical radius where the poloidal flow speed meets the poloidal Alfvén speed (A critical point, m = 1). These 3 quantities are MHD invariants along any magnetic surface anchored on a radius r_o at the disk midplane of constant magnetic flux a (along with the plasma temperature in our isothermal case).

The last quantity is the total specific energy E and is also an MHD invariant (Bernoulli integral) along the flow whenever energy is neither gained nor lost (see e.g., Casse & Ferreira 2000b when additional heating is added). It states that jet acceleration is a conversion of the (initially dominant) magnetic energy and enthalpy reservoir H into jet kinetic energy. In this paper, since turbulent pressure has been added, there is an extra energy reservoir

$$\begin{array}{c}\hfill E_{\mathrm{turb}}=-{\int }_{s^{+}}^s\frac{\partial P_{\mathrm{turb}}}{\partial s}\frac{\mathrm{d}s}{\rho },\end{array}$$(11)

where the integration starts from the disk surface (taken as the transition s⁺ between the resistive and the ideal MHD regime) and is done with the curvilinear coordinate s along a given magnetic surface. Since the turbulence is assumed to decrease vertically, E_turb >  0. However, it turns out to be always negligible in all the cases studied here so that it can be safely discarded, and using E(a) as an invariant is not too bad an approximation.

2.2. MHD turbulence prescriptions

In the standard accretion disk theory, turbulence leads to an anomalous viscosity ν_v = α_vC_sh, where α_v is a free parameter of the model (Shakura & Sunyaev 1973). Building upon this idea and anticipating that turbulence should be of MHD origin, Ferreira & Pelletier (1993) proposed that the anomalous poloidal magnetic diffusivity in the disk midplane should be written as

$$\begin{array}{c}\hfill {\nu }_m={\alpha }_m{V}_{A}h.\end{array}$$(12)

Such a scaling allows to link the two transport coefficients by introducing the (effective) magnetic Prandtl number

$$\begin{array}{c}\hfill {\mathcal{P}}_m=\frac{{\nu }_v}{{\nu }_m}=\frac{{\alpha }_v}{{\alpha }_m}{\mu }^{-1/2}.\end{array}$$(13)

It turns out that MRI-driven MHD turbulence (1) can indeed be safely described by a viscous stress tensor 𝒯 with a dominant component ${\mathcal{T}}_{r\varphi }=\rho {\nu }_{v}r\frac{\partial \mathrm{\Omega }}{\partial r}$ (Balbus 2003 and references therein) and (2) that the Shakura-Sunyaev parameter follows the scaling

$$\begin{array}{c}\hfill {\alpha }_v\simeq 8{\mu }^{1/2},\end{array}$$(14)

at least in shearing box simulations (Hawley et al. 1995; Salvesen et al. 2016). The prescriptions used in the theory are therefore consistent with MRI-driven turbulence and lead to α_m𝒫_m ∼ 8, namely a true constant independent of the field strength. Note however that its value seems to depend on simulation details such as disk aspect ratio or local (shearing box) versus global simulations. We will therefore assume that a product α_m𝒫_m around unity is acceptable.

Measuring the magnetic diffusivity in 3D simulations is a very difficult task and has been only rarely done. But both local (Lesur & Longaretti 2009; Guan & Gammie 2009; Fromang & Stone 2009) and global (Zhu & Stone 2018; Jacquemin Ide 2021; Jacquemin-Ide et al. 2021) simulations tend to show that 𝒫_m ∼ 1 − 3. For the rest of the paper, and for the sake of simplicity, we will use 𝒫_m = 1 and allow α_m to vary (F97; Jacquemin-Ide et al. 2019).

However, it is still unclear whether MHD turbulence in accretion disks leads to a fully isotropic resistivity tensor, that is, if the turbulent diffusion is the same in all directions (see however Lesur & Longaretti 2009; Gressel & Pessah 2015). Ferreira & Pelletier (1995) showed that some slight anisotropy is needed in order to provide the best conditions for cold jets to be launched. This translates into a diffusivity ν′_m in the toroidal direction possibly larger than ν_m, with the introduction of an anisotropy parameter

$$\begin{array}{c}\hfill {\chi }_m=\frac{{\nu }_m}{{\nu }_m^{\prime }}\end{array}$$(15)

(see e.g., the parametric studies in Casse & Ferreira 2000a; Jacquemin-Ide et al. 2019). Nevertheless, and again for the sake of simplicity, we will keep χ_m = 1 in our study.

The existence of a turbulent magnetic pressure in disks is a natural outcome of MRI. Its magnitude (as measured at the disk midplane) depends of course on the initial (vertical) magnetic field. Following Salvesen et al. (2016), who provided analytical fits from local simulations, we will use the scaling

$$\begin{array}{c}\hfill P_{\mathrm{turb}}=\frac{⟨\delta B^2⟩}{2{\mu }_o}=\rho \delta V_A^2={\alpha }_P\sqrt{\mu }P,\end{array}$$(16)

at the disk midplane, where ⟨δB²⟩ is some local average. The fits made in Salvesen et al. (2016) provided a value α_P ∼ 26 which is quite high and shows how important this component is in local (shearing box) simulations. This scaling has been validated in both global MHD and GRMHD simulations at low magnetization values (Jacquemin-Ide et al. 2021; Mishra et al. 2020), and has been observed in a global MHD simulation at higher magnetization values (Jacquemin Ide 2021). The values of α_P derived from these simulations usually tend to be smaller and may also depend on h/r (e.g., α_P ∼ 4; Jacquemin Ide 2021). We will therefore use α_P as a free parameter. Note that in our axisymmetric approach, this turbulent pressure only plays a role in the radial and the vertical momentum equations. In agreement with numerical simulations, we expect the dominant contribution to this pressure to arise from the toroidal magnetic field. This allows therefore to use the same term in both equations.

Since this turbulent pressure adds up with the kinetic pressure, we expect an easier vertical balance against gravity and laminar magnetic compression or, alternatively, some enhanced mass-loss rate from the disk. In the radial direction, the turbulent pressure should also provide some support against gravity. Both situations actually correspond to a modification of the local sound speed. Indeed, since the prescription amounts to P_turb ∝ P ∝ ρ, the modified sound speed writes

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_s^2=(1+{\alpha }_P{\mu }^{1/2})C_s^2=C_s^2+\delta V_A^2.\end{array}$$(17)

This modified sound speed can actually be interpreted as some fast magnetosonic speed. While in weakly magnetized WEDs compressible turbulence does not cause a significant deviation of the sound speed, this is no longer the case in strongly magnetized JEDs. As a consequence, taking this turbulence into account shifts to higher altitudes the (modified) sonic critical point that is encountered in resistive MHD (Ferreira & Pelletier 1995).

As the matter leaves the disk and enters the ideal MHD regime, the turbulence is assumed to fade away, leading to negligible anomalous transport coefficients. In practice, the resistive-ideal MHD transition occurs only when the LHS resistive terms in Eqs. (4) and (5) as well as the viscous torque in Eq. (3) become smaller than 1% of the others. But there is no physical reason to neglect the turbulent magnetic pressure, at least at the base of the ideal region. We therefore keep the turbulent term in the ideal MHD equations and this is the reason why the Bernoulli integral (Eq. 10) has a turbulent contribution. This also introduces a modification of the SM and FM critical speeds defined in ideal MHD regime, which can be written as

$$\begin{array}{c}\hfill \begin{array}{cc}& V_{\mathrm{SM},\mathrm{FM}}^2=\frac{1}{2}({\stackrel{\sim }{C}}_s^2+V_A^2\pm \sqrt{{({\stackrel{\sim }{C}}_s^2+V_A^2)}^2-4{\stackrel{\sim }{C}}_s^2{V}_{\mathit{An}}^2})\hfill \\ \hfill & V_A^2=V_{\mathit{Ap}}^2+V_{A\varphi }^2\hfill \\ \hfill & V_{\mathit{An}}^2=\frac{{({\mathit{B}}_p·\mathit{n})}^2}{{\mu }_0\rho },\hfill \end{array}\end{array}$$(18)

where the SM (resp. FM) speed is defined with the minus (resp. plus) sign, V_A is the total Alfvén speed and V_An is the projection of the Alfvén speed in the self-similar direction n (see Ferreira & Pelletier 1995 and references therein for more details). Thus, the turbulent magnetic pressure could potentially play an important role at the base of the jet where the flow must become super-SM if it introduces a significant modification of the critical phase speed V_SM. In any case, due to a prescribed decaying vertical profile of the turbulence (see below), such an effect remains quite negligible in the study done in this paper.

2.3. Self-similar ansatz

The set of MHD equations (2)–(6) forms a system of nonlinear partial differential equations. Solving such a system while keeping all dynamical terms in the equations can be done by a variable separation method (Ferreira & Pelletier 1993, 1995). Since gravity is indeed the dominant contribution, a steady-state solution over a large radial extension can only be obtained if each term follows the same scaling with the cylindrical variables (r, z). Therefore, each physical quantity A(r, z) is assumed to follow the self-similar form

$$\begin{array}{c}\hfill A(r,z)=A_o{\left(\frac{r}{r_o}\right)}^{{\zeta }_A}{f}_A(x),\end{array}$$(19)

namely a power law of the radius with exponent ζ_A times a function f_A(x) of the self-similar variable x = z/h(r), with a disk aspect ratio ϵ = h/r a constant of the radius. Actually, keeping all the dynamical terms in the equations leads to the mathematical constraint that all dimensionless quantities must remain constants of the radius. But this is the price to pay for obtaining exact MHD solutions. The subscript “o” stands for a quantity evaluated at the disk mid-plane x = 0. For instance, the magnetic flux function writes a(r, z)=a_o(r/r_o)^βΨ(x) so the shape of a magnetic surface with constant magnetic flux a(r, z)=a_o is given by r = r_oΨ^−1/β, where r_o is the anchoring radius at the disk midplane.

The self-similar ansatz replaces the system (2)–(6) with a system of algebraic equations on the radial exponents ζ_A and a system of nonlinear ordinary differential equations (ODEs) on the functions f_A(x). All the exponents ζ_A are then expressed as a function of the disk ejection index ξ, namely

$$\begin{array}{c}\hfill \begin{array}{cc}& {\dot{M}}_a\propto r^{\xi }a\propto r^{\frac{3}{4}+\frac{\xi }{2}}{u}_i\propto r^{-1/2}\hfill \\ \hfill & \rho \propto r^{\xi -3/2}{B}_i\propto r^{\frac{-5}{4}+\frac{\xi }{2}}T\propto r^{-1},\hfill \end{array}\end{array}$$(20)

where the subscript “i” stands for the three components (r, ϕ, z). These exponents are important because they tell us that the magnetic field radial distribution must be consistent with the density distribution, which of course must be closely related to the disk mass-loss rate. Any discrepancy between these exponents unavoidably leads to a time-dependent evolution of the system. This has to be kept in mind when comparing them to the outcome of numerical simulations for instance.

According to the self-similar approach, vertical profiles must also be assumed for all turbulent quantities. This is probably the least controlled prescription, as they have been very seldom provided in the literature. Moreover, the local MHD turbulence and its associated vertical profiles are most certainly self-adapting to the global accretion-ejection interrelations, while in our approach they must be prescribed ab initio without any feedback from the laminar solution. Deeper investigations using vertical profiles derived from global 3D simulations will be done in future work. Here, we limit ourselves to a simple Gaussian profile m(x)=e^−x2/2x_t2 of scale x_t as done previously (Ferreira & Pelletier 1993, 1995; Casse & Ferreira 2000a; Jacquemin-Ide et al. 2019). This leads to a decrease of all turbulence-related terms on a scale x_t ∼ 1. More precisely, the self-similar profiles of the magnetic diffusivities, viscous torque and turbulent magnetic pressure follow

$$\begin{array}{c}\hfill \begin{array}{cc}& {\nu }_m\propto {\nu }_m^{\prime }\propto m(x)\hfill \\ \hfill & \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\mathcal{T}}_{r\varphi }\right)\propto m(x)f_{\rho }(x)\hfill \\ \hfill & P_{\mathrm{turb}}\propto m(x)f_P(x).\hfill \end{array}\end{array}$$(21)

However, such simplified profiles are roughly consistent with simulations of strongly magnetized accretion disks, namely JEDs with μ ∼ 1 (Jacquemin-Ide et al. 2021; Scepi et al. 2024). On the contrary, the vertical profiles seen in weakly magnetized accretion disks, WEDs with μ ∼ 10⁻³ or less, are clearly different from simple Gaussians (Zhu & Stone 2018; Jacquemin Ide 2021). This is consistent with the finding that decaying Gaussian profiles lead to vertical oscillations of all quantities in WEDs, as discussed in Jacquemin-Ide et al. (2019). For this reason, this paper will focus only on JEDs, that is, solutions without vertical oscillations.

2.4. Summary of parameters and methodology

We numerically solve the system of ODEs on the functions f_A using a predictor-corrector integrator for stiff equations. Starting from the disk midplane, the integration propagates upward using first the resistive MHD equations. When the conditions for ideal MHD are satisfied (resistive and viscous torque negligible wrt the other terms), we switch to the ideal MHD equations. No significant jumps are seen in this transition. Once in ideal MHD, the plasma poloidal velocity is parallel to the poloidal magnetic field and the flow must cross the SM critical point. The regularity condition that allows a smooth transition to the super-SM regime fixes a quantity defined at the disk midplane. We choose to adjust the disk magnetization μ since it finely controls the disk vertical balance for a given ejection efficiency ξ (see Ferreira & Pelletier 1995 for more details).

As the super-SM flow is further accelerated (mostly by magnetic means, since our isothermal outflows have a negligible enthalpy), it must become super-A. Again, this is done by fine-tuning another physical quantity defined at the disk midplane. We choose to adjust the curvature of the poloidal field lines (parameter p associated with the toroidal electric current density J_ϕo), which is loosely related to the magnetic field bending at the disk surface. Once a new value of p has been chosen, the whole procedure of finding the new μ allowing for trans-SM solutions must be repeated (see F97 for more details). An example of a typical JED super-Alfvénic solution is illustrated in Fig. 1.

Fig. 1. 3D view of a typical super-A JED solution. The blue lines show the magnetic field lines and the white lines show the streamlines along the same magnetic surfaces, anchored at r0 = 1.0 and r0 = 5.0.

We are not looking for super-FM outflows. Not only is this last critical point deeply biased by the self-similar ansatz (unlike the first two, see discussions in Blandford & Payne 1982; Ferreira & Pelletier 1995; Bogovalov 1997; Vlahakis et al. 2000), but the regularity condition can be satisfied by adjusting the jet adiabatic index beyond the Alfvén point and thus does not provide a strong constraint on the disk physics (see Ferreira & Casse 2004). In any case, all super-A solutions become super-FM in the conventional sense, that is, the poloidal kinetic energy always becomes larger than the total magnetic energy (F97).

In summary, a magnetized accretion-ejection solution leading to super-A isothermal (γ = 1) outflows is described by the set of six parameters

$$\begin{array}{c}\hfill \begin{array}{cc}& \xi =\frac{\mathrm{dln}({\dot{M}}_a)}{\mathrm{dln}(r)}ϵ=\frac{h}{r}{\alpha }_m=\frac{{\nu }_m}{V_{A}h}\hfill \\ \hfill & {\chi }_m=\frac{{\nu }_m}{{\nu }_m^{\prime }}{\mathcal{P}}_m=\frac{{\nu }_v}{{\nu }_m}{\alpha }_P=\frac{P_{\mathrm{turb}}}{\sqrt{\mu }P}\hfill \end{array}\end{array}$$(22)

defined at the disk midplane. Four of them (α_m, 𝒫_m, χ_m, α_P) describe the MHD turbulence and are therefore unavoidable, unless future work on MHD turbulence provides their scaling as a function of, for instance, the disk aspect ratio ϵ and the disk magnetization μ. The disk aspect ratio ϵ is a real free parameter since no proper energy equation is solved. This is a drawback of the theory, but it guarantees the use of self-similarity and allows us to gain valuable insight into the effect of disk thickness on jet launching.

The physical ingredient that allows both accretion and super-A ejection is the vertical laminar magnetic field, whose strength is measured at the disk midplane by the disk magnetization parameter

$$\begin{array}{c}\hfill \mu =\frac{V_A^2}{C_s^2}=\frac{B_{\mathit{zo}}^2}{{\mu }_o{P}_o}·\end{array}$$(23)

The theory thus provides the disk ejection efficiency ξ(μ) for a set of prescribed values (ϵ, α_m, 𝒫_m, χ_m, α_P). In practice, all possible values for ξ between 10⁻⁴ and 1 are scanned (F97), looking for the value μ (and curvature p) that leads to super-A outflows.

The aim of this work is to study the effect of turbulent magnetic pressure on isothermal JED solutions. This is done in two steps. First (Sect. 3), we look at how both the JED vertical and radial balance respond to an increasingly large α_P. Since strict stationarity requires satisfying two constraints (SM and A), it is important to assess how each of them responds to the presence of turbulent pressure in the disk. Moreover, the crossing of the A point may be strongly affected by heat deposition or pressure effects above the disk (Ferreira & Casse 2004), both effects being neglected in our isothermal approach. In a second step (Sect. 4) we consider both SM and A constraints, thus providing the final parameter space and its modification with α_P.

The analyzes carried out in Sections 3 and 4 are performed for our fiducial set (ϵ = 0.1, α_m = 𝒫_m = χ_m = 1), while α_P is freely varied. Section 5 explores the effect of α_P for different disk thicknesses ϵ = 0.01, 0.1, 0.2 and for different turbulence levels α_m = 0.5, 1, 3.

3. Impact on the disk structure: super-SM parameter space

In this section, we investigate the effect of the turbulent magnetic pressure on the JED internal structure for our fiducial parameter set defined by (ϵ = 0.1, α_m = 𝒫_m = χ_m = 1). As a first approach, we impose a value of the ejection index ξ and vary the turbulence parameter α_P to understand the difference introduced by the turbulent pressure. As a second step, we explore the whole parameter set for ξ and analyze how the turbulent pressure modifies the super-SM JED parameter space.

3.1. Comparative analysis

Figure 2 shows two super-SM solutions with the same ejection efficiency ξ = 3 × 10⁻³, one obtained without turbulent pressure (α_P = 0, top panel) and corresponding to previously published solutions, and the other with turbulent pressure (α_P = 2, bottom panel). At first glance, the main visible effect is to make the disk puffier, which makes perfect sense since turbulence provides additional support against vertical compression (due to both gravity and laminar magnetic field compression). However, α_P = 0 requires μ = 0.6 to get a trans-SM outflow, whereas α_P = 2 requires μ = 3.9. This is easily understood from the disk vertical balance: the stronger the laminar field (larger μ), the stronger the vertical compression. Thus, imposing the same ξ while introducing an extra, outwardly directed vertical force, leads to the natural increase of μ. However, the impact of the turbulent pressure is more complex and affects all quantities.

Fig. 2. Poloidal cross-section of two super-SM JED solutions with the same ξ = 3 × 10−3, obtained with ϵ = 0.1, αm = 𝒫m = χm = 1, without (αP = 0, top panel) and with (αP = 2, bottom panel) a turbulent magnetic pressure. The background color represents the logarithm of the density normalized to its midplane value and the white lines show the magnetic surfaces. See text for more details.

Figure 3 shows the vertical profiles along x = z/h(r) of various disk quantities. The left panels display the classical JED solution with α_P = 0, while the right panels are for α_P = 2. In both cases, the light gray area corresponds to the resistive MHD region, the white one to the ideal MHD regime where the flow reaches the super-SM speed. Clearly, introducing a turbulent pressure leads to a thicker disk, although the thermal scale height (ϵ = h/r) remains the same: all disk quantities follow the turbulent pressure. We can understand this new disk structure from simple analytical considerations.

Fig. 3. Vertical profiles along x = z/h of the two super-SM solutions shown in Fig. 2, without (left) and with (right) turbulent magnetic pressure. From top to bottom: gas P, laminar horizontal Plam = (Br2 + Bϕ2)/2μo and turbulent Pturb magnetic pressures; velocity components ui (normalized to Keplerian speed ΩK0r0); electric current densities Ji (normalized to Bzo/μoh); laminar magnetic field components Bi (normalized to Bzo). The light gray area corresponds to the resistive MHD region and the vertical line corresponds to the position of the SM point. Note that αP = 0 requires μ = 0.6, while αP = 2 requires μ = 3.9. See text for more details.

The fact that the disk becomes puffier can be seen directly on the vertical profile of the gas pressure P (which is the same as the density profile, Fig. 3a). There is almost an inversion of the profile, which can be obtained directly by making a Taylor expansion of the vertical equilibrium Eq. (3) near the disk midplane:

$$\begin{array}{c}\hfill \frac{\partial lnP}{\partial x}\simeq -\frac{1-{\alpha }_P\sqrt{\mu }{x}_t^{-2}+\mu (p^2+q^2)}{1+{\alpha }_P\sqrt{\mu }}x=-Ax,\end{array}$$(24)

where J_ro = qB_zo/(μ_oh) and J_ϕo = pB_zo/(μ_oh) are respectively the radial and toroidal electric current densities flowing at the disk midplane. For α_P = 0, only the thermal and laminar magnetic horizontal pressure gradients (the terms 1 + μ(p² + q²)) play a role in the vertical equilibrium, so that the thermal pressure P always decreases (A >  0). In that case, the isothermal disk scale height h(r), defined as P_o = ρ_oΩ_K²h², is indeed a good proxy for the actual disk vertical magnetostatic scale height (Ferreira & Pelletier 1995). This is no longer the case with turbulent pressure. For α_P = 2, the turbulent magnetic pressure gradient $-{\alpha }_p\sqrt{\mu }{x}_t^{-2}$ balances the other two pressure gradients, so that the thermal pressure P increases near the disk midplane (A <  0). The disk material is thus pushed toward the disk surface, making it thicker and more massive, as illustrated in Fig. 2.

Another direct consequence of this extra pressure can be seen in the disk rotation rate at the midplane (Fig. 4). The deviation from the Keplerian rotation rate at the disk midplane is provided by the radial force balance

Fig. 4. Influence of αP on the super-SM JED parameter space: disk magnetization μ, deviation to Keplerian rotation δ, disk ejection efficiency ξ, thickness expansion hd = zid/h, accretion sonic Mach number ms, and laminar (jets) to turbulent (viscous) torque ratio Λ. Colored areas correspond to regions accessible by super-SM solutions. The curves represent the average solution for a given αP. See text for more details.

$$\begin{array}{c}\hfill \delta =\frac{{\mathrm{\Omega }}_o}{{\mathrm{\Omega }}_{\mathit{Ko}}}\simeq \sqrt{1-(1+{\alpha }_P\sqrt{\mu })(\frac{5}{2}-\xi ){ϵ}^2-p\mu ϵ-\frac{m_s^2}{2}{ϵ}^2},\end{array}$$(25)

where m_s = −u_ro/C_so = −u_ro/(Ω_Koh) is the accretion sonic Mach number. In the RHS of this equation, the forces leading to a deviation are (going from left to right), the radial (thermal and turbulent) pressure gradient, the laminar magnetic tension and the radial acceleration. At near equipartition field, the turbulent radial pressure gradient becomes relevant and the disk is significantly sub-Keplerian (δ ∼ 0.6 here).

Cold jets are only possible if magnetic field lines at the disk surface are bent enough, namely if B_r/B_z is greater than or on the order of unity (Blandford & Payne 1982, see Figs. 2 and 3d). Since B_r ≃ μ_o∫₀^zJ_ϕdz, the field line inclination at the disk surface is roughly B_r/B_z ≃ −∫(u_r/ν_m)dz ∼ p∫(f_ur/m)dx. As the disk inflates due to turbulence, the integral becomes larger, so p must decrease. In other words, field lines become straighter at the midplane (J_ϕo decreases), and J_ϕ is enhanced at higher altitudes. Although only moderate for the solution shown in Fig. 3c, this trend can be significantly more pronounced for other solutions.

Correspondingly, the midplane accretion speed (m_s) decreases, which implies that the total torque acting on the disk must also decrease despite an increase in the disk magnetization μ to maintain the same ejection rate ξ. Accretion occurs due to the presence of a laminar (jets) torque and a turbulent (viscous) torque, namely

$$\begin{array}{c}\hfill m_s=\frac{2q\mu }{\delta }+{\alpha }_vϵ.\end{array}$$(26)

The ratio of the jet torques to the viscous torque scales as Λ ∼ 2/(ϵα_m𝒫_m) (Eq. 35 in Jacquemin-Ide et al. 2019) and the viscous torque is therefore always negligible within our fiducial parameter set. In order for the midplane accretion speed to decrease, the jet torque −J_roB_zo must decrease despite an increase in B_zo. This dramatic decrease in J_ro (straighter field lines) is clearly seen in Fig. 3c for α_P = 2.

The evolution of the radial electric current density J_r (Eq. (5)) writes

$$\begin{array}{c}\hfill {\eta }_m^{\prime }{J}_r\simeq {\eta }_{\mathit{mo}}^{\prime }{J}_{\mathit{ro}}+r{\int }_0^z{B}_r\frac{\partial \mathrm{\Omega }}{\partial r}\mathrm{d}z.\end{array}$$(27)

The first term on the RHS of this equation results from the e.m.f. induced by the rotation of a conductive disk (of finite resistivity) amidst a large-scale vertical field (Barlow’s wheel effect). The second term is the disk feedback, always present whenever there is a velocity shear. Magneto-centrifugally driven jets are only possible if the laminar (jet) torque accelerates the ejected material, namely if F_ϕ ≃ −J_rB_z becomes positive at the disk surface (Ferreira & Pelletier 1995). This occurs only if J_r decreases vertically on the same scale, a situation achieved whenever the shear-induced term in Eq. (27) balances the Barlow term within the disk.

In a classical JED with α_P = 0, the vertical decrease of J_r on a scale height (which goes along a vertical decrease of Ω), leads to an almost linear growth of the toroidal field that achieves a value B_ϕ⁺ = −μ_o∫J_rdz ∼ −μ_oJ_roh at the disk surface, namely q ∼ 1. Therefore all three magnetic field components are comparable at the disk surface (Fig. 3d, left). But when turbulent pressure is present, rotation becomes smaller, the poloidal field lines are straighter and B_r remains smaller inside the disk, two effects that heavily reduce the shear-induced term. This pushes further up the point in the disk where acceleration takes place (F_ϕ >  0) and may even lead also to some profile inversion (barely visible in the right panel of Fig. 3d, but representative of these turbulent solutions). Indeed, since $\eta {\prime }_m{J}_r\sim \eta {\prime }_{\mathit{mo}}{J}_{\mathit{ro}}$, a decreasing magnetic diffusivity leads to an increase of the radial current density at higher altitudes. The general solution, which is self-consistent with the disk vertical equilibrium, is therefore to reduce the midplane current density J_ro, with a tendency to shift the path of the radial electric current to both disk surfaces. This reduction has far-reaching consequences because it inhibits the generation of the toroidal magnetic field B_ϕ, which is essential since it directly contributes to the accretion (torque) and ejection power (MHD Poynting flux).

3.2. SM-parameter space

We now study how the turbulent magnetic pressure modifies the super-SM parameter space of JEDs with our fiducial parameter set (ϵ = 0.1, α_m = 𝒫_m = χ_m = 1). For each value of α_P varying from 0 to 10, we explore all possible values of ξ and derive the value μ that allows the outflowing plasma to become super-SM.

It is important to note that no Alfvénic constraint has been used yet so that the field line bending (parameter p) is freely imposed. In practice, we also explore all possible values of p for each ξ. As a consequence, the number of solutions found for each ξ is variable for each α_P. This allows us to compute an “average” Super-SM solution for each α_P, which simply shows the average value for that quantity (based on the number of solutions found). The influence of α_P on the resulting parameter space is shown in Fig. 4. The shaded zones for each quantity correspond to regions that are accessible to super-SM solutions, regardless of their number. The mean solution is then represented by a point which, as α_P increases, gives an indication of how the JED structure adapts to α_P.

The actual size of the accretion disk can be very precisely defined as the altitude where the radial velocity is zero, namely the transition from inward (accretion) to outward (ejection) motion. Such a transition occurs in the resistive MHD region, and only further out does the flow transit to ideal MHD and becomes super-SM. The central right panel therefore shows the dependence of the actual disk size h_d = z_id/h as a function of α_P: clearly, the disk becomes puffier with increasing turbulent magnetic pressure, by almost a factor 2 when α_P = 10. The reason for this has been discussed previously and is related to the fact that this turbulent pressure provides additional vertical support against gravity and magnetic (laminar) compression. Note that the range where super-SM solutions are found becomes very narrow beyond α_P ∼ 7. This is an indication that finding super-A solutions at large turbulent pressures will hardly be feasible (see next section).

An unsuspected consequence of the disk becoming thicker with α_P is the fact that super-SM outflows also require increasingly large disk magnetizations μ (Fig. 4). This is actually an implication of our methodology which focuses only on non-oscillating solutions within the disk. To forbid spatial oscillations, one must require that the MRI wavelength be comparable to the disk vertical height, which translates into having comparable Alfvén and sound speeds. But our prescription of compressible MHD turbulence leads to a modification of the sound speed or, in other words, to an increase of the actual size (h_d >  1) of the disk. Thus, to recover a non-oscillating solution one needs to balance $V_A\sim {\stackrel{\sim }{C}}_s$ and to redefine an MRI-related magnetization parameter

$$\begin{array}{c}\hfill {\mu }_{\mathrm{MRI}}=\frac{V_A^2}{{\stackrel{\sim }{C}}_s^2}=\frac{\mu }{1+{\alpha }_p\sqrt{\mu }}·\end{array}$$(28)

Thus, requiring to deal only with non-oscillating solutions, namely keeping μ_MRI constant while increasing α_P, will necessarily lead to also increase μ as well.

However, this trend reaches a limit around α_P ∼ 4, where μ levels off. This is because when μ becomes too large, not only the vertical balance is compromised (overwhelming laminar magnetic compression), but also the radial one (see Eq. (25)). The top right panel in Fig. 4 shows indeed that the deviation δ = Ω_o/Ω_K from the Keplerian rotation decreases dramatically reaching almost δ ∼ 0.5. Since gravity does not change, lowering the plasma rotation rate is clearly not convenient for launching cold outflows. This value α_P ∼ 4 marks therefore the point where the super-SM region starts to shrink.

The effect of turbulent magnetic pressure on the disk ejection efficiency ξ is not obvious. At first sight, one would hope to enhance ξ thanks to this extra outward push. This effect can indeed be seen in the central left panel in Fig. 4: the minimum value for ξ does indeed increase with α_P. Close to 10⁻⁴ for α_P = 0, it cannot be less than a few 10⁻² for α_P = 10. However, the maximum value of ξ that allows for super-SM outflows also decreases: this is because the disk is rotating slower and the magnetic laminar compression becomes too strong (it scales as μ, while turbulence scales as μ^1/2). Overall, the average super-SM JED solution (solid line) does not seem to see its ejection efficiency deviating significantly from the classical (α_P = 0) value, namely ξ ∼ 10⁻² for isothermal flows.

Another counterintuitive effect of the turbulent magnetic pressure is to reduce the midplane accretion velocity, despite the increase in disk magnetization μ. The bottom left panel in Fig. 4 indeed shows a decreasing sonic Mach number m_s with increasing α_P (although it remains larger than, or on the order of, unity). This is due to the decreasing influence of the jet torque wrt the turbulent (viscous) torque. The bottom right panel clearly shows a steady decrease in the ratio Λ of these two torques with α_P. This is a direct consequence of the turbulent JED becoming unable to generate a strong toroidal magnetic field at the disk surface. As discussed previously, not only is the radial electric current density J_ro (parameter q in Eq. (26)) flowing in the midplane to be lowered in the presence of turbulent pressure, but the vertically decreasing magnetic diffusivity does not allow for a large toroidal field at the disk surface. As a consequence, jets play a less and less dominant role in accretion.

To conclude, our analysis of the inclusion of a turbulent magnetic pressure in JEDs shows that the disk becomes puffier, more sub-Keplerian, and less conductive, with field lines more vertical inside the disk. The generation of magnetic bending (B_r) and shear (B_ϕ) occurs mostly at the (more elevated) disk surface, with no significant change in the disk ejection efficiency. Overall, a magnetic turbulent pressure reduces the impact of the jets on the accretion.

4. Impact on the disk+jets: Super-A solutions

4.1. Comparative analysis

Figure 5 shows a 3D view of the two JED solutions discussed previously, obtained with our fiducial parameter set (ϵ = 0.1, α_m = 𝒫_m = χ_m = 1). Our goal is to assess how the presence of a turbulent magnetic pressure impacts the dynamics of the outflowing super-A jet and, in turn, affects the underlying disk. Since both solutions have the same ejection efficiency (ξ = 3 × 10⁻³), they should have roughly the same Alfvén radius r_A, since it depends mostly on ξ (F97, see also below). But a simple look at the figure reveals that the turbulent JED launches much less open (i.e., more collimated) jets with an Alfvén point located at a higher altitude z_A. Also, in agreement with the super-SM analysis, the magnetic field lines appear much less twisted, a sign that the jets may eventually become exhausted of any poloidal electric current, leading to a vanishing B_ϕ.

Fig. 5. 3D view of the JED solutions shown in Fig. 2 with the same ejection index (ξ = 3 × 10−3), obtained with our fiducial set (ϵ = 0.1, αm = 𝒫m = χm = 1), without (αP = 0, μ = 0.6) and with (αP = 2, μ = 3.9) turbulent magnetic pressure. The color represents the logarithm of the density normalized to its midplane value, the white lines show the magnetic field lines and the blue lines show the streamlines along the same magnetic surfaces, anchored at ro = 2 and ro = 3. The white dotted lines show the position of the Alfvén point on these surfaces.

Figure 6 shows the profiles of several quantities along a magnetic surface (constant a) for the two previous super-A solutions. The top row displays the various characteristic velocities normalized to the Keplerian speed at the anchoring radius r_o: the flow speed in the self-similar direction V = u_p ⋅ n, the modified sound speed and the modified A, SM, and FM phase speeds defined by Eq. (18) (see also Ferreira & Pelletier 1995). Due to the Gaussian decrease of the turbulence, the influence of the turbulent pressure on the sound speed rapidly decreases as the plasma is lifted out of the disk. Thus, once in ideal MHD, only the real thermal energy remains (${\stackrel{\sim }{C}}_s=C_s$) and the sound speed remains constant along the field lines, in agreement with our assumption. This can also be seen in the Bernoulli integral (bottom panels, Fig. 6): the turbulent energy content E_turb remains completely negligible. In the isothermal case, the enthalpy is also always negligible and the jets are therefore cold, the dominant energy being initially stored in the magnetic field E_mag and progressively transferred to the kinetic energy E_kin of the plasma.

Fig. 6. Vertical profiles along x = z/h(a) of the two super-A solutions shown in Fig. 5, without (left) and with (right) turbulent magnetic pressure. Top: characteristic velocities normalized to the Keplerian velocity at the disk midplane (flow speed V, modified sound speed ${\stackrel{\sim }{C}}_s$, slow and fast magnetosonic speeds VSM, VFM and Alfvén speed VAn). Bottom: contributions to the Bernoulli invariant E normalized to ΩKo2ro2/2. The light gray area corresponds to the resistive MHD region and the vertical lines corresponds to the position of the SM and A point respectively. See text for details.

4.2. A-parameter space

The top left panel in Fig. 7 shows how the usual JED parameter space ξ(μ) evolves with α_P. For each value of α_P, we recover the characteristic almost vertical behavior, namely the ability to achieve a wide range of ejection efficiencies ξ in a very narrow interval of μ. Three important features emerge immediately as α_P increases: (1) μ must be increasingly large, (2) ξ decreases drastically, (3) until no super-A solution can be found beyond α_P = 4. The fact that μ must increase with α_P is a consequence of keeping the MRI wavelength of nearly the same size as the puffed disk (as described by Eq. (28)). It is therefore a consequence of the disk vertical equilibrium, dealt with the SM critical point. However, it is interesting to note the existence of a “sweet spot” for α_P = 1, associated with the largest interval in ξ. But, globally, the super-A parameter space, which is a subset of the previous SM parameter space, shrinks with α_P, with a global decrease of the maximum ejection efficiency ξ_max. Such a decrease in ξ_max was already present in the SM parameter space (see Fig. 4), but it is greatly amplified by the Alfvénic constraint.

Fig. 7. Influence of the magnetic turbulent pressure αP on super-A solutions obtained with our fiducial set. Top left: JED parameter space ξ(μ). Top right: magnetic lever arm λ as function of ξ. The dotted lines correspond to λ = a[1 + 1/(2ξ)] with a = 0.9 and a = 0.5. Bottom left: fraction b of the released accretion power transferred to the two jets as function of the magnetic surface rotation ω anchored at a radius ro. Bottom right: maximum cylindrical radius rt/ro reached before recollimation, as function of λ. The dotted line is a fit rt/ro = 1.08λ0.526.

This can be understood by the fact that the disk rotation rate δ decreases with α_P, resulting in a global decrease in the rotation rate

$$\begin{array}{c}\hfill \omega =\frac{{\mathrm{\Omega }}_{\ast }(a)}{{\mathrm{\Omega }}_{\mathit{Ko}}}\end{array}$$(29)

of the magnetic surfaces (bottom left panel in Fig. 7). Since the magneto-centrifugal acceleration is the only means for launching our isothermal (cold) outflows, less and less disk material can be accelerated to super-A speeds as α_P increases. This interpretation is validated by the neat correlation b(ω) for all α_P, where b = 2P_jet/P_acc is the ratio of the total power (kinetic, thermal, potential, and magnetic) carried away by the two jets to the released accretion power. The reason is simply that the main energy content is magnetic and that it scales directly with Ω_*, namely P_jet = ∫ρu_pE ⋅ dS ≃ ∫S_MHD ⋅ dS, where the MHD Poynting flux is S_MHD = −Ω_*rB_ϕB_p/μ₀. As guessed from the study of the SM parameter space, magnetic turbulence quite significantly reduces the amount of laminar magnetic energy feeding the two jets.

The total specific angular momentum L(a) carried along a magnetic surface (Eq. (9)), normalized by the Keplerian angular momentum at the anchoring radius r_o, provides the so-called magnetic lever arm parameter λ introduced by Blandford & Payne (1982). This important jet quantity is actually related to the disk parameters

$$\begin{array}{c}\hfill \lambda =\frac{L(a)}{{\mathrm{\Omega }}_{\mathit{Ko}}{r}_o^2}\simeq \omega +\frac{1}{2\xi }\frac{\mathrm{\Lambda }}{1+\mathrm{\Lambda }},\end{array}$$(30)

which is a generalization of the classical relation λ ≃ 1 + 1/2ξ (F97). Note that most JED solutions have Λ ≫ 1. The top right panel in Fig. 7 shows how the dependency λ(ξ) evolves with α_P. The two bracketing dotted lines are the curves λ = a(1 + 1/2ξ), with a = 0.5 and a = 0.9. They show that the analytical link is not broken and that any measure of λ can give a fair estimate of ξ. The main difficulty, however, is the decrease of the disk rotation rate (hence ω) as α_P increases, leading to a decrease in λ at constant ξ.

Nevertheless, turbulent magnetic pressure in the disk leads to a dramatic difference in the jet collimation properties, as illustrated in the bottom right panel in Fig. 7. This plot displays the jet widening, namely the ratio r_t/r_o where r_t is the maximum cylindrical radius reached by a magnetic surface anchored at r_o before recollimating toward the axis (F97), as a function of the magnetic lever arm λ and for different values of the turbulent magnetic pressure parameter α_P. The gray dotted line is a numerical fit r_t/r_o = 1.08λ^0.526, namely r_t ≃ r_A. Here we recover the existence of the sweet spot at α_P = 1, which allows us to obtain the widest types of behavior, namely from solutions no larger than r_A to solutions opening up to 10⁵r_o. But it can be readily seen that as α_P increases, there are fewer and fewer solutions that propagate much further than the Alfven point. At α_P = 4, all solutions found merely become super-A before undergoing recollimation. Such solutions, termed “current-free” in F97, become super-A at the expense of almost all the available electric current flowing out of the disk.

Figure 8 shows the evolution of several jet parameters as a function of the turbulent magnetic pressure intensity α_P. The colored areas correspond to regions accessible by super-A solutions, while the curves represent the evolution of the mean solution obtained by averaging the values of all the solutions found for a given α_P. These figures allow us to better understand the influence of the disk turbulent pressure on the jet parameters and dynamics.

Fig. 8. Influence of αP on the super-A JED parameter space: jet mass load κ, rotation rate δ of the magnetic surface, magnetic lever arm λ, total specific energy e carried along the magnetic surface, altitude zA/ro of the Alfvén point, and the electric current gA still available at that point. Colored areas correspond to regions accessible by super-A solutions. The curves represent the average solution for a given αP. See text for more details.

The jet mass load, defined as (Blandford & Payne 1982)

$$\begin{array}{c}\hfill \kappa =\eta (a)\frac{{\mathrm{\Omega }}_{\mathit{Ko}}{r}_o}{B_{\mathit{zo}}}\simeq \frac{\xi m_s}{\mu }=\xi {\alpha }_{m}p{\mu }^{-1/2},\end{array}$$(31)

is an MHD invariant that can be related to the disk quantities using ξ. The value of p, which is a measure of the toroidal electric current density flowing at the disk midplane and loosely related to the field lines bending at the disk surface, is determined by the smooth crossing of the Alfvén critical point. The larger κ the more kinetically dominated the outflow is. It turns out that for JED (namely strongly magnetized) solutions κ ∼ ξ (F97; Casse & Ferreira 2000a), so we show below only ξ. The top left panel in Fig. 8 proves that only very tenuous, magnetically dominated cold outflows can become super-A. This is amplified by both the increase in μ and the corresponding decrease of ξ as α_P increases.

As shown above, this is mostly due to the dramatic decrease in the rotation of the magnetic surface ω with α_P. Indeed, as μ increases, the midplane plasma rotation δ decreases and, despite the MHD acceleration starting at higher altitudes, the plasma hardly becomes super-Keplerian (Ferreira & Pelletier 1995). As a consequence, the centrifugal term appearing in the radial momentum equation does not overcome gravity and only the radial component of the laminar magnetic force F_r = J_ϕB_z − J_zB_ϕ must become capable of driving an outwardly directed flow. While super-SM outflows are possible at α_P = 10 with δ ≃ 0.55, no super-A solution is found with α_P >  4 (for our fiducial parameter set), the lowest value being ω ≃ 0.63 (corresponding to δ ≃ 0.63 as well).

For cold outflows, and neglecting any additional energy E_turb associated with the turbulent magnetic pressure, Eq. (10) leads to a dimensionless Bernoulli invariant

$$\begin{array}{c}\hfill e=\frac{2E(a)}{{\mathrm{\Omega }}_{\mathit{Ko}}^2{r}_o^2}\simeq \omega \lambda -2+\frac{{\omega }^2}{2}·\end{array}$$(32)

The exact value e(a) is shown at the central right panel in Fig. 8. Not surprisingly for cold outflows, it fulfills the above equation and closely follows the evolution of λ, which is itself directly related to the ejection efficiency ξ. Our outflows from turbulent disks are not only gravitationally unbound (e >  0), but they are also highly magnetized (e ≫ 1). The ratio of the (laminar) MHD Poynting flux to the matter-energy flux at the SM point writes

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{SM}}={\frac{-{\mathrm{\Omega }}_{\ast }r{B}_{\varphi }{B}_p}{(\frac{u^2}{2}+H){\mu }_0\rho u_p}|}_{\mathrm{SM}}\simeq 2\omega (\lambda -1)\end{array}$$(33)

for cold outflows (Jacquemin-Ide et al. 2019), and is therefore much larger than unity. One might naively expect such highly magnetized jets to propagate far from the disk. However, this is clearly not the case, as they do not extend much beyond the Alfvén point, as shown in the bottom left panel of Fig. 8.

Fig. 9. Influence of αP on the laminar magnetic field components at the SM critical point for super-A solutions obtained with our fiducial parameter set. Left: ratio Br/Bz showing the bending. Right: ratio Bϕ/Bz showing the magnetic shear. Colored areas correspond to regions accessible by super-A solutions. The curves represent the mean solution for a given αP.

This can be understood by looking at how the altitude z_A/r_o of the Alfvén point evolves with α_P (Fig. 8). As the turbulent magnetic pressure increases, the disk becomes puffier, and its magnetization μ increases. This implies a larger Alfvén velocity and thus would require a longer acceleration length to allow the jet poloidal speed to reach the Alfvén velocity. But a look at the SM parameter space (hence the disk vertical balance, Fig. 4) shows that the ejection efficiency tends to remain nearly the same so that the cylindrical Alfven radius r_A/r_o does not change much. As a consequence, increasing the acceleration distance can only be done by increasing z_A/r_o. Clearly, the interplay between the disk vertical balance (SM point) and the jet acceleration (A point) requires playing on both sides, namely increasing both r_A (thus decreasing ξ) and z_A. For a given ξ, this ends up in a more collimated outflow emitted from a turbulent JED than from a nonturbulent JED, as illustrated in Fig. 5.

As the Alfvén point is pushed further away from the disk with increasing α_P, the field lines at the disk surface become more and more straight (Fig. 9, left panel). This appears to be a geometrical consequence, but it is due to the fact that in the sub-A region, the plasma is dominated by the magnetic field and in particular by the poloidal component. This magnetic tension always tends to close the magnetic surface, thus preventing plasma acceleration. For our fiducial parameter set, at α_P >  4 the magnetic bending is such that the energetic criterion for cold outflows (requiring B_r/B_z >  1 at the disk surface, Blandford & Payne 1982) is no longer satisfied and no super-A solution can be found.

This decrease in the magnetic field bending at the disk surface goes along with a decrease in the toroidal magnetic field component (Fig. 9). This is consistent with the decrease in the total energy feeding the jets (fraction b, Fig. 7). Note that this behavior is due to the fact that we only consider JED structures with μ larger than unity. Indeed, for weakly magnetized accretion-ejection structures (WEDs) with μ ≪ 1, a large toroidal (laminar) magnetic field could be generated (see Jacquemin-Ide et al. 2019), but these solutions are beyond the scope of the present paper.

This tremendous decrease in magnetic shear at the disk surface goes along with a much lower acceleration efficiency, consistent with the increasing distance z_A from the Alfvén point. As a proxy for this efficiency, one can use the function $g=1-\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\ast }}$, which measures the angular velocity drift between the plasma and the magnetic surface (Pelletier & Pudritz 1992). It increases monotonically from g ≪ 1 at the disk surface, while its value at the Alfvén point,

$$\begin{array}{c}\hfill g_A\simeq \frac{r_A{B}_{\varphi ,A}}{r_{\mathrm{SM}}{B}_{\varphi ,\mathrm{SM}}}\end{array}$$(34)

determines the fate of the jet (F97). Indeed, g_A is the ratio of the electric poloidal current I = rB_ϕ (enclosed in the magnetic surface) still flowing at the Alfvén point to the one provided at the disk surface (taken here as the SM point). Solutions with g_A >  1/2 are current-carrying outflows where the acceleration to super-A speeds has been done with not much expense of the initial magnetic energy: at the Alfvén point, the magnetic field still has a large reservoir of angular momentum (and energy) to transfer back to the plasma. On the contrary, solutions with g_A <  1/2 have exhausted all the available magnetic energy just to reach super-A speeds.

While current-carrying jets open wide and propagate much farther out than the Alfvén point, current-free jets with low g_A tend to recollimate immediately upon reaching this point. This is mostly due to the dominant effect of the poloidal magnetic tension (F97). Figure 8 shows that as α_P increases, g_A decreases down to ∼0.2, leading therefore mostly to current-free solutions at larger α_P for our fiducial parameter set.

5. Parameter study

5.1. Effect of the disk geometrical thickness

The results described in the previous sections were obtained for a fiducial set (ϵ = 0.1, α_m = 𝒫_m = χ_m = 1). In this section, we analyze the effect of the magnetic turbulent pressure on isothermal JED solutions by exploring two other values of the disk aspect ratio: ϵ = 0.01, characteristic of an optically thick, geometrically thin disk and ϵ = 0.2, characteristic of a slim/optically thick or geometrically thick, optically thin disk (see e.g., Marcel et al. 2018). The other parameters are kept constant, namely (α_m = 𝒫_m = χ_m = 1). As before, we discuss first the impact of the turbulent magnetic pressure on the disk structure (SM constraint only) and then the constraint imposed to obtain outflows with super-A speeds (both SM and A constraints).

Figure 10 shows the evolution of several disk quantities as a function of the turbulent magnetic pressure intensity α_P for different disk geometrical thicknesses ϵ. Colored areas correspond to regions accessible by super-SM solutions. The curves represent the evolution of the average solution obtained by averaging the values of all solutions found for a given α_P.

Fig. 10. Influence of αP on the super-SM JED parameter space for various disk aspect ratios ϵ = h/r obtained with (αm = 𝒫m = χm = 1): disk magnetization μ, deviation to Keplerian rotation δ, disk ejection efficiency ξ, thickness expansion hd = zid/h, accretion sonic Mach number ms, and laminar (jets) to turbulent (viscous) torque ratio Λ. Colored areas correspond to regions accessible by super-SM solutions with different aspect ratios: ϵ = 0.01 (blue), ϵ = 0.1 (yellow), ϵ = 0.2 (red). The corresponding curves represent the average solution for a given αP.

As expected the magnetization μ first increases and then saturates in the same way regardless of the geometrical thickness ϵ. This is because, as the disk inflates (i.e., h_d = z_id/h) due to the increasing importance of the magnetic turbulent pressure, the condition for the absence of MRI channel modes (spatial oscillations) remains fulfilled (Eq. (28)).

This increase in the laminar magnetic field has tremendous consequences on the midplane disk rotation rate. In JEDs (strongly magnetized accretion disks), the dominant contribution to the deviation from the Keplerian rotation law is the laminar magnetic radial tension, which is proportional to μϵ (Eq. (25)). But the radial gradient of the turbulent magnetic pressure provides a contribution that scales as α_Pμ^1/2ϵ², so it also becomes non-negligible at large μ and high α_P. It is therefore not surprising that the deviation remains tiny in thin accretion disks, namely δ = Ω_o/Ω_Ko ∼ 1 at ϵ = 0.01 (blue, Fig. 10). But as the disk thickens, the disk becomes more and more sub-Keplerian as α_P increases. For a thick disk with ϵ = 0.2 (red, Fig. 10), δ ∼ 0.2 at α_P = 4 and no super-SM solution can be found beyond this turbulence level. For even thicker disks (ϵ >  0.2) no super-SM (isothermal) solution can be found, even with α_P = 0 (F97; Casse & Ferreira 2000a). The reason is quite simple: for cold jets, only magnetic driving allows acceleration so that whenever the plasma rotation becomes too low, the gravitational attraction always wins.

For the values of ϵ explored here, there is no significant effect on the values of the average disk ejection efficiency ξ around ∼10⁻² (central left panel of Fig. 10). However, the range of allowed ξ for super-SM outflows changes significantly with ϵ and with α_P. It is noteworthy to see how much larger the SM parameter space is for the thin disk, even with values ξ >  1 possible. We will come back to this point in Sect. 6.2.

The midplane accretion sonic Mach number m_s = −u_ro/C_s = α_mpμ^1/2 remains always larger than unity (supersonic accretion) for all values of the disk thickness. However, the dependency with the turbulent magnetic pressure is different for thin (ϵ = 0.01) and for slim or thick disks (ϵ = 0.1 and 0.2): while in the thin case m_s increases with α_P, it decreases in thicker disks (Fig. 10). For thick disks, which become even thicker as α_P increases, it has been argued that both radial and azimuthal electric currents are pushed toward the disk surface, leading to a correspondingly decrease of their value at the disk midplane (see Sect. 3.1). As a consequence, the equatorial field line curvature p and jet torque become smaller leading to a decrease of m_s with α_P (see also Eq. (26)). This also translates into structures where less and less disk angular momentum is carried away by the outflows (the torques ratio Λ decreases). This is no longer true in thin disks, where the turbulent magnetic pressure has a smaller impact on the disk final size (since Λ scales with 1/ϵ). As a consequence, the current densities can remain mostly at the disk equatorial plane, which translate into magnetic configurations much more bent and sheared (both p and q of order unity). The increase of m_s with α_P is therefore mostly a consequence of its dependence as μ^1/2. We recover here the classical result obtained with α_P = 0 that both the laminar (jet) magnetic torque and the turbulent (viscous) torque scale as μ^1/2, so that their ratio Λ does not evolve with α_P and is proportional to ϵ⁻¹ (Casse & Ferreira 2000a; Jacquemin-Ide et al. 2019).

Figure 11 shows the evolution with α_P of the parameter space for super-A solutions for different geometrical thicknesses ϵ. As expected, it is much smaller than the SM parameter space studied previously. There is no solution beyond α_P = 1 for ϵ = 0.2 and beyond α_P = 4 for ϵ = 0.1. The reason lies in the fact that the magnetic surface rotates much too slowly to magnetically drive cold (isothermal) outflows (see the ω panel), in agreement with our analysis of the disk rotation rate at the equatorial plane.

Fig. 11. Influence of αP on the super-A JED parameter space for various disk aspect ratio ϵ = h/r obtained with (αm = 𝒫m = χm = 1): disk ejection efficiency ξ, rotation rate ω of the magnetic surface, magnetic lever arm λ, disk magnetization μ, altitude zA/ro of the Alfvén point, and the electric current gA still available at that point. Colored areas correspond to regions accessible by super-A solutions with different aspect ratios: ϵ = 0.01 (blue), ϵ = 0.1 (yellow), ϵ = 0.2 (red). The corresponding curves represent the average solution for a given αP.

Despite the supplementary Alfvénic constraint, there is no better contrast in ξ (or in the magnetic lever arm λ) achieved for different disk thicknesses. Indeed, the allowed ranges in ξ stay nearly the same with a typical ξ ∼ 10⁻³. Note that the jet mass load κ still always follows closely ξ.

The most striking feature is the very different behavior of the thin disk solution (ϵ = 0.01, blue), as shown in the plots of the disk magnetization μ, the altitude z_A/r_o of the Alfvén point and the energy reservoir g_A. The typical thin disk solution (curve in the plot) maintains almost the same magnetization μ regardless of α_P, in strong contrast to the other two thicker solutions. This can be understood by the fact that the slight increase in disk thickness due to α_P has not significantly changed the total pressure². As a consequence, the accretion-ejection can adapt its magnetic geometry while keeping roughly the same μ. This has several consequences. Instead of allowing more mass (larger ξ) to be ejected, the system favors solutions where the field lines become more and more bent at the disk surface as α_P increases. While a typical solution with α_P = 0 has B_r/B_z ∼ 2 at the SM point, it increases until B_r/B_z ∼ 5 for α_P = 10. Magneto-centrifugal acceleration is thus enhanced, and typical outflows reach the Alfvén point at about the same location (see bottom left panel in Fig. 11), but without expending much of the initial magnetic energy. In contrast to thicker disks with ϵ = 0.1 and 0.2, whose jets’ remaining energy g_A at the Alfvén point decreases with α_P, jets from a thin disk with ϵ = 0.01 see their energy increase and remain larger than 1/2.

5.2. Effect of the magnetic diffusivity

We now analyze how the turbulence level parameter α_m, which determines both anomalous magnetic diffusivities and viscosity, affects the JED response to the turbulent magnetic pressure α_P. Since classical JEDs are always dominated by the laminar jets torque, we will mostly interpret the role played by α_m as a modification of the disk magnetic diffusivities. In this context, it has already been shown that α_m must be large enough (of order unity) to build up a significant toroidal field at the disk surface (F97; Jacquemin-Ide et al. 2019). For this reason, we restrict our study to two other values α_m = 0.5 and α_m = 3, while keeping the same set (ϵ = 0.1, 𝒫_m = χ_m = 1).

The parameter space obtained only with the SM constraint is displayed in Fig. 12 for the three values of α_m. There is no significant modification introduced by α_m on the necessary increase of the disk magnetization μ as a response to an increase of the magnetic turbulent pressure α_P (same reason based on μ_MRI). The same lack of strong influence can be seen in the allowed ranges of the disk ejection efficiency ξ (always < 0.1), which are nearly the same for all α_m.

Fig. 12. Influence of αP on the super-SM JED parameter space for various magnetic diffusivity levels αm and (ϵ = 0.1, 𝒫m = χm = 1): disk magnetization μ, deviation to Keplerian rotation δ, disk ejection efficiency ξ, thickness expansion hd = zid/h, accretion sonic Mach number ms, and laminar (jets) to turbulent (viscous) torque ratio Λ. Colored areas correspond to regions accessible by super-SM solutions with different magnetic diffusivity levels: αm = 0.5 (blue), αm = 1 (yellow), αm = 3 (red). The corresponding curves represent the average solution for a given αP.

The turbulence level parameter α_m has a rather monotonous effect on several disk quantities, regardless of the amount of magnetic turbulent pressure. Figure 12 shows that as α_m increases, the deviation from the Keplerian rotation rate decreases (i.e., δ increases). This is because increasing the magnetic diffusivity leads to straighter poloidal and toroidal field lines. As a consequence, the radial magnetic tension becomes less pronounced, and δ increases (see Eq. (25)). For the same reason, the disk is less vertically pinched by the laminar magnetic pressure gradient (Eq. (24)), leading to more diffusive disks already becoming thicker. This effect is naturally maintained as α_P increases and the disk becomes puffier (h_d = z_id/h increases, central right panel in Fig. 12). Finally, since both laminar (jets) and turbulent (viscous) torques scale as μ^1/2, one expects Λ ∝ 1/α_m, so that increasing α_m increases the effect of the turbulent torque, as observed indeed. The decrease in Λ is then amplified as α_P increases. Thus, the stronger the turbulence (large α_m and α_P), the less energy and angular momentum the jets extract from the JED.

How the accretion sonic Mach number m_s behaves with turbulence is less straightforward. For α_m = 0.5 and α_m = 1, the behavior is similar and monotonous: a larger α_m leads to a more diffusive disk and a larger m_s. This can be interpreted as the fact that it is easier to diffuse through the vertical field lines so that the accretion speed can be higher. As α_P increases, leading to a decrease in the toroidal magnetic field at the disk surface and thus a decrease of the jet torque, m_s decreases. But for a rather large α_m = 3, it is quite surprising to get a smallerm_s at α_P = 0 (red curve, Fig. 12). This shows that the field lines inside the disk become too straight so the toroidal field is mostly generated only at the disk surface. As a consequence, accretion is lowered at the disk midplane (lack of torque) despite a very large diffusion. The portion where m_s increases with α_P is due to the increase of the disk magnetization μ, which boosts the torques. Eventually, this is overcome by the reducing effect of the magnetic turbulent pressure on the toroidal field, leading to a decrease of m_s with α_P.

For high magnetic diffusivity and turbulent magnetic pressure, no solution can be found because the ejection of matter is so delayed that the resistive MHD region expands until the SM point penetrates the disk. This configuration cannot provide any solution because the magnetization of the disk is no longer constrained.

Finally, Figure 13 presents the impact of α_m on the JED response to an increase of α_P, taking into account both SM and Alfvénic constraints. As seen in the SM-parameter space, the influence of α_m is quite similar for α_m = 0.5 and α_m = 1.

Fig. 13. Influence of αP on the super-A JED parameter space for various turbulence level αm obtained with (ϵ = 0.1, 𝒫m = χm = 1): disk ejection efficiency ξ, rotation rate ω of the magnetic surface, magnetic lever arm λ, disk magnetization μ, altitude zA/ro of the Alfvén point, and the electric current gA still available at that point. Colored areas correspond to regions accessible by super-A solutions for different magnetic diffusivity levels: αm = 0.5 (blue), αm = 1 (yellow), αm = 3 (red). The corresponding curves represent the average solution for a given αP.

Although super-A solutions can be found up to α_P = 4 for α_m = 1 and only α_P = 1 for α_m = 0.5, the typical ejection index remains qualitatively similar (a few 10⁻³ or mass load κ ∼ 10⁻³, leading to a typical magnetic lever arm λ ∼ 10²). The reason why super-A solutions cannot be found with increasing α_P is that because the magnetic field curvature at the disk surface becomes too small to allow cold ejection: the ratio B_r/B_z at the SM point decreases below unity. As α_P increases, the disk magnetization μ increases accordingly, the rotation rate ω of the magnetic surface decreases (following the decrease of the midplane rotation δ), and the altitude of the Alfvén point z_A/r_o increases (best seen with α_m = 1). All these behaviors are consistent with the increasing disk magnetization and the increasing dominant effect of the poloidal magnetic tension on the jet sub-A region dynamics. This leads to jets that do not propagate much further than the Alfvén point, as illustrated by the decrease seen in g_A (bottom right, Fig. 13).

As already seen with the SM constraint, the behavior is profoundly different for large α_m = 3. The first surprising result is that no super-A cold solution has been found for α_P = 0. The reason is that at such high magnetic diffusivity, the field lines in the disk are too straight to allow cold ejection. However, as α_P increases, the generation of electric currents at the disk surface layers is enhanced, leading to an increase in the ratio B_r/B_z and thus to the existence of super-A outflows. Thus, without turbulent magnetic pressure, no jet would be produced here. Quite interestingly, and in agreement with the SM constraint, the disk magnetization μ increases less with α_P, allowing to achieve super-A outflows closer to the disk (smaller altitude z_A/r_o, of course increasing as μ increases) and potentially propagating much further out from the Aflvén point (g_A increasing and greater than 1/2). However, for α_P >  4, ω becomes too small, and producing super-A jets requires to use of more and more magnetic energy, leading to a decrease in g_A until there is no super-A solution.

6. Discussion

6.1. Summary and caveats

We have studied the effect of including a turbulent magnetic pressure in a strongly magnetized accretion-ejection structure (a JED). It turns out that this extra term, which appears only in the vertical and radial momentum equations, has a tremendous effect on the JED parameter space and the cold (negligible enthalpy) jet properties.

Qualitatively, JEDs with a turbulent magnetic pressure become puffier and appear much less electrically conductive. They tend to force both radial (J_r) and toroidal (J_ϕ) electric current densities to flow at the disk surface, leading to much more straight magnetic (laminar) field configurations within the disk, the bending required for energetically allowing cold jets to be launched occurring mostly at the disk surface (see Fig. 14). As a consequence, the laminar toroidal magnetic field is also smaller and developed mostly in the upper layers of the disk. This has several consequences. While the midplane accretion remains supersonic, the jet torque is lowered within the disk and the accretion speed is faster in the disk upper layers. Overall, jets from these new solutions carry away less accretion energy and disk angular momentum and also become more collimated due to a stronger poloidal magnetic tension acting in the sub-A region.

Fig. 14. Schematic representation of the effect of the turbulent magnetic pressure. The turbulent magnetic pressure makes the JED puffier and less electrically conducting, pushing both radial and toroidal electric currents to flow at the disk surface. Field lines become more straight within the disk and their bending as well as their shear occur at the surface. Asymmetry between the two bipolar jets is much easier to achieve because the two poloidal currents can be decoupled within the disk.

There are however two important caveats that need to be discussed, both intimately related to MHD turbulence. The first one is our choice of the vertical profiles of the anomalous (turbulent) quantities. For the sake of simplicity, but also due to the lack of thorough analyzes in the literature, we choose to use simple Gaussian profiles decreasing at a fixed scale. Such Gaussian profiles are not realized in 3D global MHD simulations when the disk is weakly magnetized (Zhu & Stone 2018; Jacquemin Ide 2021; Jacquemin-Ide et al. 2021). This numerical situation corresponds to WED solutions found at magnetizations below μ ∼ 0.1, where spatial oscillations (MRI channel modes) have been found (Jacquemin-Ide et al. 2019). A Gaussian profile, or a seemingly decreasing profile on a disk scale height, has been however observed in simulations at large magnetization with μ ∼ 1 (Jacquemin Ide 2021). This is the main reason why we restricted our analysis to strongly magnetized JEDs. Our finding that α_P leads to a systematic shrinking of the parameter space may be a consequence of the fact that we do not adjust the vertical scale of the anomalous transport coefficients, that is, take into account the feedback of the mean flow on the turbulence. Indeed, the prescribed thermal disk scale height h, used as the vertical scale of the turbulence, is no longer a valid proxy for the real disk scale height h_d >  h. This can lead to solutions that did not become super-SM or super-A because they ejected matter too early inside the disk. Taking this new disk scale height into account could be achieved through an iterative method. However, this approach would introduce additional numerical complexities, while the actual profiles of turbulent diffusivities are close to, but not exactly, Gaussian. We have decided to postpone this analysis to future work, where we will use profiles educated by 3D numerical simulations.

A second important aspect that severely limits our parameter space is the simplifying assumption of isothermal (cold) outflows. As α_P increases, the tension associated with the vertical magnetic field becomes increasingly strong, leading to sub-Keplerian accretion disks. This dramatic decrease of the centrifugal drive (and the larger ϵ = h/r, the smaller ω), due to the field lines becoming more and more straight at the disk surface, forbids steady-state ejection to take place. Since Blandford & Payne (1982) it is well known that enough poloidal field curvature (bending) is a requirement only for purely magnetically driven outflows, that is, outflows with a negligible enthalpy. Casse & Ferreira (2000b) have shown that any additional heat input at the disk surface layers deeply modifies the disk vertical structure, allowing for instance to enhance the disk mass loss (they provided a warm super-A solution with ξ ∼ 0.5). Such surface heating may naturally arise from irradiation coming from a hot central source, which could be either the innermost disk regions or for instance the accretion shock onto the hard surface of the central object. But, as these authors point out, it may also be due to turbulence itself. Now, it is interesting to see that recent 3D numerical GRMHD simulations do display signs of an anomalous dissipation above the disk (see Fig. 12 in Scepi et al. 2024). It is still unclear how such a local dissipation, which is withdrawn from the numerical domain at each time step, actually affects the outflow dynamics.

Clearly, MHD turbulence in the disk does more than just provide anomalous transport coefficients and a turbulent magnetic pressure. Another example of a physical process neglected here is turbulent mass diffusion, which could greatly enhance ξ. Note that a numerical (due to the grid) mass diffusion can nevertheless occur in nonturbulent simulations, possibly leading to discrepancies between theoretical expectations and numerical experiments (see Sect. 3.5 in Murphy et al. 2010). However, addressing all these points is clearly beyond the scope of the present paper.

6.2. Implications

Putting aside the previous limitations, let us now discuss our results more quantitatively along with their implications. For our fiducial parameter set (ϵ = 0.1, α_m = 𝒫_m = χ_m = 1), classical (i.e., α_P = 0) isothermal JED solutions typically have ξ ∼ 10⁻², μ ∼ 0.5 and propagate much beyond the Alfvén point. When turbulent pressure is included up to its maximum value α_P = 4, the ejection efficiency ξ decreases around 10⁻³, μ increases up to 10, the disk is significantly sub-Keplerian (δ ∼ 0.6), jets carry only ∼25% of the released accretion power and recollimate immediately after crossing the Alfvén point. It is therefore highly doubtful that these solutions could be representative of powerful magnetically dominated astrophysical jets propagating much further away from the source. But dealing with the fate of recollimation jets cannot be done within the framework of self-similarity (see for instance Jannaud et al. 2023).

However, regardless of the limitations discussed previously, the turbulent pressure is highly dependent on the disk aspect ratio ϵ. Specifically, thin disk solutions with ϵ = 0.01 (and α_m = 𝒫_m = χ_m = 1) are achieved even up to α_P = 10, the disk remains nearly Keplerian with μ ∼ 1 and jets propagate much further out than thicker disks. They do however have a very low ejection efficiency ξ ∼ 10⁻³. Additionally, the turbulent pressure allows for new solutions with high diffusivity level α_m that were not achievable before. For instance, solutions with α_m = 3 (and ϵ = 0.1, 𝒫_m = χ_m = 1) are sub-Keplerian with μ ∼ 2 and still propagate far beyond the Alfvén point with ξ ∼ 10⁻².

It therefore appears that obtaining jets that propagate far beyond the Alfvén point can be achieved by playing with ϵ and/or α_m for JEDs at μ near equipartition with an important turbulent pressure. It turns out that increasing α_m is the best way to enhance the jet mass load, that is, increasing ξ or decreasing the magnetic lever arm λ. More importantly, classical JED solutions (α_P = 0) with α_m = 3 could not be found (in agreement with Jacquemin-Ide et al. 2019). In other words, the presence of a magnetic turbulent pressure (α_P ≠ 0) appears necessary to obtain more massive, steady-state super-A cold jets at large α_m. The physical reason behind this efficient collaboration between α_m and α_P lies in the accretion-ejection geometry depicted in Fig. 14.

Global 3D MHD simulations of highly magnetized jets, referred to as MAD in the literature, do provide some hints that such a geometry is indeed observed in the numerical experiments. The question now is how to determine in which theoretical parameter space those simulations ended up. In fact, the vast majority of numerical papers do not provide the value of the disk magnetization μ. Instead, the beta plasma is usually given but it barely allows to recover our control parameter μ. Indeed the beta plasma is most of the times the ratio of disk averages (not only averages done at the midplane) and, more importantly, it always includes the turbulent magnetic pressure.

To our knowledge, the only GRMHD paper allowing to infer μ, α_P and ϵ is the work of Scepi et al. (2024). The values provided are measured at 7 r_g so quite close to the black hole, where general relativistic effects might be already at work. Nevertheless we indicate them as they provide some insight. For the three simulations, the parameter sets (ϵ, μ, α_P) are: (0.03, 0.070, 15.1), (0.1, 0.045, 4.7) and (0.3, 0.015, 5.7). It can be seen from these values that α_P seems to vary with the disk aspect ratio ϵ, being larger in the thinner disk. It is also worth noting that the JED solutions considered in this paper are actually much more magnetized. This is an indication that we should seek for solutions with large α_m values (or relax the non-oscillatory condition). Unfortunately, there is no measure done in the numerical simulations of the α_m parameter. If we assume that scalings derived from shearing box simulations remain representative of global simulations, then one should seek for α_m𝒫_m ∼ 8 (see Eq. (14)). Looking for new solutions of that kind that could best reproduce MAD simulations (including the vertical profiles and outflow behavior) is postponed for future work.

Another long-standing discrepancy between theoretical models and numerical experiments is the value of the disk ejection efficiency. Commonly measured values for ξ varied from 0.5 to 1 (see e.g., McKinney et al. 2012), while in more recent simulations Manikantan et al. (2024) found ξ = 0.4 for example. Scepi et al. (2024) measure values that vary greatly whether the wind is cooled down or allowed to stay warm (due to turbulent heat dissipation beyond the disk surface). For warm winds, ξ goes from 0.12 for ϵ = 0.03 to ξ = 1.07 for ϵ = 0.3, while for cooled winds, ξ ∼ 0.3 − 07 for this range in ϵ. Larger disk ejection efficiencies for warm outflows are consistent with theoretical expectations (Casse & Ferreira 2000a) and our results here show that the allowed values for ξ may indeed depend on the disk aspect ratio.

In agreement with the disk global energy budget established between an inner radius r_in and an outer radius r_out, the released accretion power, defined as

$$\begin{array}{c}\hfill P_{\mathrm{acc}}={\left[\frac{GM{\dot{M}}_a(r)}{2r}\right]}_{r_{\mathrm{out}}}^{r_{\mathrm{in}}}=\frac{GM{\dot{M}}_a(r_{\mathrm{in}})}{2r_{\mathrm{in}}}(1-{\left(\frac{r_{\mathrm{out}}}{r_{\mathrm{in}}}\right)}^{1-\xi })\end{array}$$(35)

is positive only if ξ <  1. Any value of ξ larger than unity shows that either the disk is not yet in steady-state, or that some extra energy is being transferred at r_in from the central object. For instance, super-SM solutions with ξ up to 10 are perfectly possible (see Fig. 10), but their Bernoulli invariant e is negative and they cannot become super-A. In numerical simulations, that would correspond to unsteady, very massive bursts of matter that fall back to the disk at larger radii.

To ensure that simulations have indeed reached a steady-state, one should seek for the consistency between the disk ejection efficiency ξ derived from the radial profile of Ṁ_a(r) and other disk quantities, such as those in Eq. (20). It is for instance somewhat troublesome that simulations in Scepi et al. (2024) display a magnetic flux a ∝ r^3/4 (see their Fig. 4) regardless of the disk aspect ratio ϵ, which in theory would necessarily require ξ ≪ 1. Note however that these measurements are done quite close to the black hole (up to 10 r_g), where deviations from the Newtonian radial profiles are obviously expected.

Another possibility, related to turbulence, that might explain the very large disk ejection efficiencies reported in global 3D simulations is the presence of a turbulent mass diffusion. Addressing this last point may however be particularly tricky in numerical simulations. In any case, deeper comparisons between 3D numerical experiments and theory must be carried out.

Note that the discrepancy in ξ also holds for 2D MHD simulations of “alpha” disks, namely simulations where anomalous transport coefficients have been used to mimic turbulence (Casse & Keppens 2002; Zanni et al. 2007; Tzeferacos et al. 2009; Murphy et al. 2010; Stepanovs & Fendt 2016). No clear explanation has been proposed yet, except for a possible numerical effect (see discussion in Murphy et al. 2010).

7. Conclusion

In this paper, we have presented a major upgrade on the theory of highly magnetized (near equipartition) accretion-ejection structures, also called JEDs. We analyzed the effect of turbulent magnetic pressure on the disk structure and its impact on the overall interrelationship between the disk and its jets. Although this pressure appears only in the radial and vertical momentum equations, it plays a major role and deeply affects the JED parameter space.

The turbulent magnetic pressure makes the disk puffier and less electrically conductive, forcing both radial and toroidal electric currents to flow at the disk surface. Field lines become more straight inside the disk, and their bending and shearing occur mostly at the surface. Accretion remains supersonic but the speed increases at the disk upper layers. For the usual values of the turbulence parameters explored so far in classical (i.e., no turbulent magnetic pressure α_P = 0) JED models, the inclusion of turbulent pressure leads to a dramatic decrease in the cold jet power and angular momentum extraction from the disk.

However, recent 3D global simulations tend to show that MHD turbulence provides simultaneously large anomalous magnetic diffusivities and viscosity as well as a large magnetic turbulent pressure. When combined in the semi-analytical model, these large values allow for the existence of new super-Alfvénic jet solutions that were previously impossible to obtain.

The next obvious step is to better characterize the vertical profiles involved in MHD turbulence. They may indeed differ from the Gaussian profiles used in this work, which focuses on strongly magnetized disks (μ ∼ 1). Also associated with these profiles, some heat dissipation may be always present, providing an additional enthalpy reservoir for the outflows and allowing, for instance, more mass to be ejected. Taking into account non-Gaussian profiles is an absolute necessity in the case of WEDs (with μ ≪ 1). This is postponed for future work.

Our mathematical description of accretion-ejection structures allows us to incorporate both laminar and turbulent effects. However, it is now necessary to use more accurate values and more realistic vertical profiles for the turbulence, based on the results of converged 3D numerical simulations.

The JED and WED concepts have emerged on the theoretical side, while on the computational side MAD and SANE configurations have received much attention. A direct comparison between these costly 3D numerical experiments and steady-state theory has finally become feasible.

Our results demonstrate that JEDs provide a state-of-the-art mathematical description of their MAD numerical counterparts. However, further efforts from both sides are needed to firmly establish this point.

Acknowledgments

The authors thank Nicolas Scepi for sharing his results and interesting discussions and the anonymous referee for the very insightful comments that have improved the quality of the paper. NZ and JF acknowledge financial support from the CNES French space agency and PNHE, PNPS programs of French CNRS. JJ acknowledges support by the NSF AST-2009884, NASA80NSSC21K1746 and NASA XMM-Newton 80NSSC22K0799 grants.

References

All Figures

	Fig. 1. 3D view of a typical super-A JED solution. The blue lines show the magnetic field lines and the white lines show the streamlines along the same magnetic surfaces, anchored at r0 = 1.0 and r0 = 5.0.
In the text

	Fig. 2. Poloidal cross-section of two super-SM JED solutions with the same ξ = 3 × 10−3, obtained with ϵ = 0.1, αm = 𝒫m = χm = 1, without (αP = 0, top panel) and with (αP = 2, bottom panel) a turbulent magnetic pressure. The background color represents the logarithm of the density normalized to its midplane value and the white lines show the magnetic surfaces. See text for more details.
In the text

Fig. 3. Vertical profiles along x = z/h of the two super-SM solutions shown in Fig. 2, without (left) and with (right) turbulent magnetic pressure. From top to bottom: gas P, laminar horizontal Plam = (Br2 + Bϕ2)/2μo and turbulent Pturb magnetic pressures; velocity components ui (normalized to Keplerian speed ΩK0r0); electric current densities Ji (normalized to Bzo/μoh); laminar magnetic field components Bi (normalized to Bzo). The light gray area corresponds to the resistive MHD region and the vertical line corresponds to the position of the SM point. Note that αP = 0 requires μ = 0.6, while αP = 2 requires μ = 3.9. See text for more details.

In the text

Fig. 4. Influence of αP on the super-SM JED parameter space: disk magnetization μ, deviation to Keplerian rotation δ, disk ejection efficiency ξ, thickness expansion hd = zid/h, accretion sonic Mach number ms, and laminar (jets) to turbulent (viscous) torque ratio Λ. Colored areas correspond to regions accessible by super-SM solutions. The curves represent the average solution for a given αP. See text for more details.

In the text

Fig. 5. 3D view of the JED solutions shown in Fig. 2 with the same ejection index (ξ = 3 × 10−3), obtained with our fiducial set (ϵ = 0.1, αm = 𝒫m = χm = 1), without (αP = 0, μ = 0.6) and with (αP = 2, μ = 3.9) turbulent magnetic pressure. The color represents the logarithm of the density normalized to its midplane value, the white lines show the magnetic field lines and the blue lines show the streamlines along the same magnetic surfaces, anchored at ro = 2 and ro = 3. The white dotted lines show the position of the Alfvén point on these surfaces.

In the text

Fig. 6. Vertical profiles along x = z/h(a) of the two super-A solutions shown in Fig. 5, without (left) and with (right) turbulent magnetic pressure. Top: characteristic velocities normalized to the Keplerian velocity at the disk midplane (flow speed V, modified sound speed ${\stackrel{\sim }{C}}_s$, slow and fast magnetosonic speeds VSM, VFM and Alfvén speed VAn). Bottom: contributions to the Bernoulli invariant E normalized to ΩKo2ro2/2. The light gray area corresponds to the resistive MHD region and the vertical lines corresponds to the position of the SM and A point respectively. See text for details.

In the text

Fig. 7. Influence of the magnetic turbulent pressure αP on super-A solutions obtained with our fiducial set. Top left: JED parameter space ξ(μ). Top right: magnetic lever arm λ as function of ξ. The dotted lines correspond to λ = a[1 + 1/(2ξ)] with a = 0.9 and a = 0.5. Bottom left: fraction b of the released accretion power transferred to the two jets as function of the magnetic surface rotation ω anchored at a radius ro. Bottom right: maximum cylindrical radius rt/ro reached before recollimation, as function of λ. The dotted line is a fit rt/ro = 1.08λ0.526.

In the text

Fig. 8. Influence of αP on the super-A JED parameter space: jet mass load κ, rotation rate δ of the magnetic surface, magnetic lever arm λ, total specific energy e carried along the magnetic surface, altitude zA/ro of the Alfvén point, and the electric current gA still available at that point. Colored areas correspond to regions accessible by super-A solutions. The curves represent the average solution for a given αP. See text for more details.

In the text

	Fig. 9. Influence of αP on the laminar magnetic field components at the SM critical point for super-A solutions obtained with our fiducial parameter set. Left: ratio Br/Bz showing the bending. Right: ratio Bϕ/Bz showing the magnetic shear. Colored areas correspond to regions accessible by super-A solutions. The curves represent the mean solution for a given αP.
In the text

Fig. 10. Influence of αP on the super-SM JED parameter space for various disk aspect ratios ϵ = h/r obtained with (αm = 𝒫m = χm = 1): disk magnetization μ, deviation to Keplerian rotation δ, disk ejection efficiency ξ, thickness expansion hd = zid/h, accretion sonic Mach number ms, and laminar (jets) to turbulent (viscous) torque ratio Λ. Colored areas correspond to regions accessible by super-SM solutions with different aspect ratios: ϵ = 0.01 (blue), ϵ = 0.1 (yellow), ϵ = 0.2 (red). The corresponding curves represent the average solution for a given αP.

In the text

Fig. 11. Influence of αP on the super-A JED parameter space for various disk aspect ratio ϵ = h/r obtained with (αm = 𝒫m = χm = 1): disk ejection efficiency ξ, rotation rate ω of the magnetic surface, magnetic lever arm λ, disk magnetization μ, altitude zA/ro of the Alfvén point, and the electric current gA still available at that point. Colored areas correspond to regions accessible by super-A solutions with different aspect ratios: ϵ = 0.01 (blue), ϵ = 0.1 (yellow), ϵ = 0.2 (red). The corresponding curves represent the average solution for a given αP.

In the text

Fig. 12. Influence of αP on the super-SM JED parameter space for various magnetic diffusivity levels αm and (ϵ = 0.1, 𝒫m = χm = 1): disk magnetization μ, deviation to Keplerian rotation δ, disk ejection efficiency ξ, thickness expansion hd = zid/h, accretion sonic Mach number ms, and laminar (jets) to turbulent (viscous) torque ratio Λ. Colored areas correspond to regions accessible by super-SM solutions with different magnetic diffusivity levels: αm = 0.5 (blue), αm = 1 (yellow), αm = 3 (red). The corresponding curves represent the average solution for a given αP.

In the text

Fig. 13. Influence of αP on the super-A JED parameter space for various turbulence level αm obtained with (ϵ = 0.1, 𝒫m = χm = 1): disk ejection efficiency ξ, rotation rate ω of the magnetic surface, magnetic lever arm λ, disk magnetization μ, altitude zA/ro of the Alfvén point, and the electric current gA still available at that point. Colored areas correspond to regions accessible by super-A solutions for different magnetic diffusivity levels: αm = 0.5 (blue), αm = 1 (yellow), αm = 3 (red). The corresponding curves represent the average solution for a given αP.

In the text

Fig. 14. Schematic representation of the effect of the turbulent magnetic pressure. The turbulent magnetic pressure makes the JED puffier and less electrically conducting, pushing both radial and toroidal electric currents to flow at the disk surface. Field lines become more straight within the disk and their bending as well as their shear occur at the surface. Asymmetry between the two bipolar jets is much easier to achieve because the two poloidal currents can be decoupled within the disk.

In the text