Magnetically-driven jets and winds from weakly magnetized accretion disks
J. Jacquemin-Ide J. Ferreira G. Lesur


firstpage–lastpage
Persistence of Excitation in  Reproducing Kernel Hilbert Spaces, Positive Limit Sets, and Smooth Manifolds
    [
    December 30, 2023
==========================================================================================================





Semi-analytical models of disk outflows have successfully described magnetically-driven, self-confined super-Alfvénic jets from near Keplerian accretion disks. These Jet Emitting Disks are possible for high levels of disk magnetization μ defined as μ=2/β where beta is the usual plasma parameter. In near-equipartition JEDs, accretion is supersonic and jets carry away most of the disk angular momentum. However, these solutions prove difficult to compare with cutting edge numerical simulations, for the reason that numerical simulations show wind-like outflows but in the domain of small magnetization. In this work, we present for the first time self-similar isothermal solutions for accretion-ejection structures at small magnetization levels. We elucidate the role of MRI-like structures in the acceleration processes that drive this new class of solutions. The disk magnetization μ is the main control parameter: massive outflows driven by the pressure of the toroidal magnetic field are obtained up to μ∼ 10^-2, while more tenuous centrifugally-driven outflows are obtained at larger μ values. The generalized parameter space and the astrophysical consequences are discussed. We believe that these new solutions could be a stepping stone in understanding the way astrophysical disks drive either winds or jets. Defining jets as self-confined outflows and winds as uncollimated outflows, we propose a simple analytical criterion based on the initial energy content of the outflow, to discriminate jets from winds. We show that jet solution are achieved at all magnetization level, while winds could be obtained only in weakly magnetized disks that feature heating. 






 





black hole physics –
accretion, accretion disks –
magnetohydrodynamics (MHD) – 
ISM: jets and outflows –
X-rays: binaries







§ INTRODUCTION




Jets are observed from a wide variety of astrophysical objects. They are emitted from the central regions of young stellar objects where a protostar is being born <cit.>, from the central core of active galactic nuclei and quasars (<cit.> and references therein ), and also from the compact object of a binary system, be it a black hole, a neutron star or even a white dwarf <cit.>. These jets are detected by different means in each environment, mostly in radio around compact objects (interpreted as self-absorbed synchroton emission) and in emission lines (from radio to optical) in young stellar objects. While jets from compact objets are relativistic, those from young forming stars have speeds ranging from 100 to 600 km/s. They do however share some properties: they are both supersonic, have small opening angles already close to the source and they exhibit a tight correlation with the underlying accretion disk <cit.>. These are evidences for two important aspects: (i) the acceleration process must also be related to the jet confinement issue, calling therefore for an initial self-collimation instead of an external confinement; (ii) the mass and/or also power that are feeding the jets must be related to the accretion activity.  


It is now accepted that large scale magnetic fields anchored on a rotating object are a necessary ingredient for launching self-confined outflows <cit.>. This rotating object could be either the central object or the surrounding accretion disk. However, given the universality of the process, it sounds reasonable to rely on their common denominator, namely the accretion disk as proposed by <cit.>. The physical ingredients that need to be included in a self-consistent steady-state accretion-ejection model are then the following: (1) a near-keplerian disk of plasma surrounding a central mass; (2) a large scale vertical magnetic field threading the disk; (3) the possibility of mass diffusion through the field, so that the disk material can accrete onto the central object while leaving the magnetic field behind. Accretion is then driven by two possible mechanisms: the turbulent torque as proposed by <cit.> and the laminar torque due to the jets themselves <cit.>. While the latter is a natural consequence of the presence of magnetic jets, the former requires the existence of some self-sustained turbulence within the disk. 


Since the seminal work of <cit.>, it is now well known that magnetized disks are unstable to the magnetorotational instability (hereafter MRI). The saturated state of the MRI is a 3D magnetohydrodynamic (MHD) turbulence giving rise to an anomalous radial transport of angular momentum, that can indeed be described by a turbulent viscosity (see <cit.> and references therein). 
While it is now quite commonly argued that magnetic fields are present in accretion disks, the question of their topology remains open. 
Indeed, 3D global numerical simulations show that as long as there is no large scale vertical magnetic field, MRI is present (thus accretion proceeds) but no jet is launched (<cit.>, see however <cit.>). 
Thus, the process of launching jets from accretion disks, that will carry away mass, energy and angular momentum, requires the presence of a large scale vertical field. 
In that case then, whatever the relative importance of the jet torque to the turbulent torque, steady accretion is achieved only if the disk material is allowed to diffuse through the magnetic field. 

The origin of this diffusion remains one of the less studied aspects of accretion-ejection theory. One possibility is the influence of non-ideal MHD processes, such as Ohmic resistivity, ambipolar diffusion and Hall effect <cit.>. However, while relevant in outer regions of protostellar accretion disks, innermost disk  regions and disks around active galactic nuclei and X-ray Binary are ionized enough and these effects vanish. One needs therefore to rely on another source for diffusion, the MHD turbulence itself. 


Building upon this idea, <cit.> analyzed the general conditions for designing self-similar steady-state models of accretion-ejection structures. Within this model, most of the disk plasma accretes in a resistive (turbulent) MHD region around the disk midplane, while a fraction is deviated vertically and is smoothly connected to an ideal MHD zone, where it crosses the usual MHD critical points. <cit.> obtained the first outflow solutions becoming super-Slow Magnetosonic and showed, for the first time, that the required magnetic field needs to be smaller than but close to equipartition with the total (gas plus radiation) pressure, see also <cit.>. This result was then generalized to super-Alfvénic <cit.> and super-Fast Magnetosonic <cit.> jets. In these highly magnetized solutions, termed Jet Emitting Disks (hereafter JED), the inclusion of a turbulent viscous torque appears to play no significant role <cit.>, most of the disk angular momentum being extracted by the jets. On the other hand, the mass loaded in the jets happens to be highly dependent on the thermodynamic conditions at the disk upper layers: allowing for some heat deposition (coronal heating) is shown to lead to a significant enhancement of the ejected mass <cit.>. 


Most results shown in these early works have been confirmed by other groups, in particular using 2D numerical simulations of "alpha"  disks, where viscosity and magnetic diffusivity are prescribed using an alpha law (see eg. <cit.>). The main caveat of these alpha-disk simulations is their possible inconsistency with MHD turbulence, since all anomalous transport coefficients are parametrized. In order to probe the analytical results on accretion-ejection structures, 3D global MHD simulations of turbulent accretion disks with large scale magnetic fields must be done. 

But achieving reliable simulations of this kind is a fantastic task. Indeed, MHD turbulence in a vertically stratified disk needs to be properly followed to make sure that the simulations have converged. Moreover, when a vertical magnetic field is included, mass loss is systematically observed from the disk surface, loaded field lines become bent and ejection is obtained. As a consequence, the size of the computational domain needs to be large enough so that boundary conditions do not affect (or not too much) the outcome of the simulation. These are the main reasons why shearing box simulations could hardly address MRI with a non zero net magnetic flux (see for instance <cit.> and references therein). Therefore, for quite a long time, the main focus of MRI studies was the measurement of the Shakura-Sunayev alpha parameter describing the turbulent viscosity and the influence of non ideal MHD effects. It became clear only recently that these MRI-driven outflows would carry away some angular momentum as well, possibly affecting the structure of weakly ionized accretion disks <cit.>. 


Converged global simulations of accretion disks threaded by a weak vertical magnetic field have been obtained for a plasma beta around 10^4 <cit.>. Super-Alfvénic flows are systematically obtained and, in the case of <cit.>, there are even indications of some collimation occurring within the domain. Despite the presence of the vertical laminar torque due to these outflows, most of the disk angular momentum is transported outwardly in the radial direction. As a consequence, the power carried away by these "winds"  remains a small fraction of the released accretion power. But the existence of super-Alfvénic outflows from weakly magnetized accretion disks is in contradiction with the analytical (JED) model. Although self-similarity introduces some biases in the flow solutions, it does allow to take into account all dynamical terms. Thus, the results of those global simulations motivated us to revisit the analytical theory of accretion-ejection structures and to seek for new outflow solutions at low disk magnetization levels.  
 
 
The paper is organized as follows. Section 2 provides the governing equations and assumptions allowing to describe steady-state accretion disks driving jets. The JED parameters are introduced and those describing the MHD turbulence are discussed in the framework of MRI simulations. It will be shown that a condition, used to obtain the previously published solutions, must be disregarded in order to be consistent with global MRI simulations. New solutions, obtained at low magnetization levels, are then indeed naturally obtained. Section 3 describes the new parameter space of super-Slow Magnetosonic (SM) flows. Although there is no MRI in our steady-state calculations, it is shown that these winds are a natural outcome of MRI-like modes or MRI channel flows in stratified unbounded flows. Super-Alfvénic flows are then a subset of these super-SM solutions and their properties are shown in Section 4. Section 5 analyses the effects of the turbulence parameters on the solutions and in particular on the disk mass loss, with a possible clear distinction between jets and winds. Some caveats of our study and comparison with other works are then presented in Section 6. We conclude in Section 7 by discussing some astrophysical implications.  





§ DESCRIBING ACCRETION-EJECTION STRUCTURES







 §.§ Governing equations



Accretion-Ejection structures are described in the framework of axisymmetric MHD. The plasma velocity and magnetic field can be decomposed into poloidal and toroidal components, 𝐮=+Ω r and 𝐁=+B_ϕ respectively. The poloidal magnetic field can then be written 

    = 1/r a ×

where a(r,z) is related to the vector potential 𝐀 by a=rA_ϕ and is the magnetic flux function. A poloidal magnetic surface is defined by a constant vertical flux, namely a(r,z)=a(r_o,0) where r_o is the cylindrical anchoring radius of the magnetic surface. The magnetic field topology is then assumed bipolar with an even symmetry with respect to the disk equatorial plane. This translates into an even function a(r,z) in z and an odd function B_ϕ.  The set of MHD equations are then the following <cit.>

    (ρ𝐮)=0 
       ρ(𝐮_𝐩·)𝐮_𝐩=-P+ρΦ_G+ J_ϕ a/r -  (r B_ϕ)^2/2 μ_o r^2
       [ρΩ r^2  -  rB_ϕ/μ_0 -rT]=0
       ν_mJ_ϕ=1/μ_0×
       (ν'_m/r^2rB_ϕ)=1/r(B_ϕ-Ω r)

with ρ the total mass density, P the thermal pressure of the plasma, 𝐉=1/μ_0𝐁 the plasma electric current density, ν_m and ν'_m are the anomalous resistivities (poloidal and toroidal respectively), Φ_G= - GM/√(r^2+z^2) the gravitational potential of the central object of mass M (ignoring the self gravitation of the disk) and T= 𝒯_rϕ𝐞_r where 𝒯_rϕ is the radial stress of turbulent origin, associated to an anomalous viscosity ν_v <cit.>.

This set of equations is closed using the equation of state for a perfect gas. In this paper, the temperature will be assumed to be remain constant along each magnetic surface (isothermal solutions). The specific form of this equation of state is detailed in Appendix <ref>.





 §.§ MHD turbulence and transport coefficients


As discussed earlier, the disk is assumed to be fully turbulent and that such a turbulence can be described using a mean field approach with anomalous transport coefficients. This is the alpha-disk description introduced by <cit.>. The disk is then defined as the densest region around the equatorial plane, of scale-height h(r), inside which accretion and turbulence are taking place. 

The real disk scale height is provided by the vertical balance between gravity, magnetic forces and plasma pressure support. It is thus only known once a full solution (including the energy equation) is computed. It is however practical to define the hydrostatic scale height such that C_s = Ω_K h, where C_s is the midplane isothermal sound speed and  Ω_K = √(GM/r^3) the Keplerian angular velocity. While accurately providing the scale height of standard accretion disks, it is a slight overestimate in the case of strongly magnetized disks, usually by a factor 2 or so <cit.>. In this work, the disk aspect ratio ϵ= h/r = C_s/Ω_Kr will be used as a free parameter for prescribing the temperature at the disk equatorial plane.   

Whatever the instability that triggers and sustains the MHD turbulence, we assume that it translates into a viscosity as well as a magnetic diffusivity, namely that the turbulent electromotive force is proportional to the mean electric current density. All transport coefficients are then assumed to vanish outside the disk, the jet region being described within the ideal MHD regime. Our description must then allow for a smooth transition from a resistive viscous MHD regime (the disk) to an ideal MHD regime (the jet) on a few disk height scales. For simplicity, will use the same gaussian profile for the vertical behavior of all transport coefficients (see Appendix <ref>).   

The amplitude of the turbulent transport coefficients is then prescribed as follows: 
 

  * Viscosity ν_v: following  <cit.> it is chosen as ν_v = α_v C_s h, where C_s is the midplane sound speed and α_v is the usual turbulence parameter. 
   

  * Poloidal diffusivity  ν_m: it is the magnetic diffusivity acting on the poloidal magnetic field (or J_ϕ), responsible for the bending of the magnetic field (as measured for instance by the ratio B_r^+/B_z at the disk surface) and allowing steady accretion through the magnetic field. Following the initial prescription made in <cit.>, we use ν_m= α_m V_A h, where V_A is the Alfvén velocity at the disk midplane and α_m a constant.  


  * Toroidal diffusivity  ν'_m: it is the magnetic diffusivity acting on the toroidal magnetic field (or J_p). It is responsible for the magnetic shear (as measured by the ratio -B_ϕ^+/B_z at the disk surface) and allowing for a steady rotation. Given our lack of knowledge on these aspects in turbulent MHD disks, we follow  <cit.> and use ν'_m= ν_m/ χ_m, where χ_m is a measure of a possible anisotropy. 



In the first studies of JEDs <cit.>, only the jet torque was taken into account and the relevant turbulent parameters were α_m and χ_m. When the viscous torque was included in the equations <cit.>, the Shakura-Sunyaev α_v parameter was computed using the effective magnetic Prandtl number  𝒫_m = ν_v/ν_m, namely 

    α_v= α_m 𝒫_m μ^1/2

where μ= V_A^2/C_s^2= B^2/μ_o P is the disk magnetization measured at the disk midplane (note that μ=2/β where β is the usual plasma beta only in accretion disks dominated by the gas pressure). Hence, assuming a constant 𝒫_m leads to a Shakura-Sunyaev α_v parameter scaling with the disk magnetization, whereas the magnetic diffusivity α_m remains a constant. All results published so far on JEDs verify this property. It is striking to realize that MRI studies actually provide 

    α_v = α_o μ^1/2

where α_o ≃ 7 according to <cit.> and the scaling law α_v ∝μ^1/2 previously identified <cit.>. Identifying Eq.(<ref>) to Eq.(<ref>) leads to α_m = α_o/𝒫_m which is indeed a constant. Now, measuring the magnetic diffusivity in turbulent accretion disks is a tricky task, especially in global simulations. To our knowledge, only few works attempted it and found, using different approaches, 𝒫_m and χ_m both slightly larger than, but of the order, unity <cit.>, consistent with estimates done in global simulations <cit.>. Note however that the anisotropy parameter χ_m has been measured only in one configuration (radial diffusion of a vertical field compared to the vertical diffusion of a toroidal field). Because of the absence of data in the more general case (vertical diffusion of a radial field for instance) the value of the anisotropy parameter χ_m will be used here as a free parameter.

To conclude, using 𝒫_m, χ_m and α_m ≥ 1 as free constants and α_v= α_m 𝒫_m μ^1/2 is actually consistent with our current knowledge of MRI-driven turbulence.





 §.§ Parameters and method of resolution


The full set of MHD equations (<ref>)-(<ref>) is solved using a self-similar Ansatz. Since gravity is expected to be the leading energy source, all other quantities will have to follow the same mathematical dependency. Such an approximation brings of course a lot of caveats but it allows to take into account all dynamical terms in the equations. For a newtonian potential, this translates into seeking solutions of the form
 
    A(r,z) = A_o (r/r_o)^ζ_Af_A(x)

where f_A(x) is the profile of any quantity A, expressed with the self-similar variable x= z/h(r)= z/ϵ r. In this radial self-similarity, a constant x corresponds to a cone and the values of the exponents ζ_A are obtained by solving algebraic equations. Note that, as with many other disk models, all quantities are then a power-law of the radius. In a JED, disk mass loss must be allowed and quantified. This is done by allowing a radial dependence of the disk accretion rate 
 
    Ṁ_a (r)= - 2 r ∫_-h^+h dz ρ u_r ∝ r^ξ

where 1≥ξ >0 is the disk ejection efficiency. The bigger ξ the larger the amount of ejected matter, while ξ=0 describes a standard accretion disk. The other exponents ζ_A can then be expressed as function of ξ, including the magnetic field distribution <cit.>. See Appendix <ref> for more details. 

The global energy budget of geometrically thin accretion-ejection structures, established between an inner disk r_i and an outer disk r_e 
writes P_acc= 2 P_jet + P_diss. In this expression, P_jet is the power leaving the disk and carried away by each jet while P_diss is the power that is released within the disk through turbulent dissipation and giving rise to the disk luminosity. The accretion power is

    P_acc =  [ G M Ṁ_a(r)/2r ]^r_i_r_e = G M Ṁ_a(r_i)/2r_i -  G M Ṁ_a(r_e)/2r_e

and its amplitude depends thereby on how much mass is leaving the disk.  In this global budget, advection of energy into the central object scales as P_adv∝ (h/r)^2 P_acc and has thus been neglected, as well as any external source of energy (such as irradiation from a central source). As a consequence, requiring that jet launching and disk luminosity are both powered by the release of mechanical energy leads to the constraint ξ <1. We will come back to this constraint later. 

The self-similar antsatz allows for a full description in the poloidal plane and is therefore required when dealing with the deviation of the flow from accretion to ejection. One consequence of self-similarity is that all local dimensionless quantities must be real constants, defining thereby the parameters of the solution. The list of the 7 JED parameters, evaluated at the disk midplane, is then
 
 
    [             ϵ = h/r                          α_m = ν_m/V_Ah                        μ = B_0^2/μ_0P_0;      ξ =dln_a/dln r                          𝒫_m  = ν_v/ν_m                     p = J_ϕ 0/B_0/μ_0 h;                                             χ_m  = ν_m/ν'_m                                         ]


Here, p controls the toroidal electric current density at the disk mid plane. It is a measure of the bending of the magnetic surface at the disk mid plane, resulting from the interplay between advection and turbulent diffusion. It provides also a rough estimate of the bending of the field lines at the disk surface, namely B^+_r/B_z ∼ h μ_o J_ϕ 0/B_z ∼ p.  

In the above list, three parameters (α_m, 𝒫_m, χ_m) are unavoidable as they describe the MHD turbulence. We will use 𝒫_m=1 in this work and explore the other two, for the values χ_m=[0.01,0.1,1,2] and α_m=[0.8,1,2,8]. Our reference set of parameters will be (α_m=1, χ_m=1). The disk aspect ratio ϵ should be computed using the energy equation. But, as said before, this is not done here and we will instead fix it to the common value ϵ=0.1.

Since two parameters will be constrained by the crossing of two critical points (see below), this leaves one free parameter and we choose the disk ejection efficiency ξ. Thus, for a given set (ϵ, α_m, χ_m, 𝒫_m), we compute the values of the toroidal current p and disk magnetization μ that are necessary to allow for a JED with the desired value ξ. The ejection index ξ will thus be varied from the smallest value allowing for a solution to ξ=1. For illustrative purposes, we will display the resulting parameter space showing ξ as function of the disk magnetization μ. Table <ref> contains a list of the disk parameters evaluated at the disk mid-plane as well as their type (constrained or free). We have also included other useful quantities that might be needed for the comprehension of the dynamical properties. 



Thanks to the method of variable separation, the set of PDEs is transformed into a set of ODEs on the functions f_A (see Appendix <ref> for their expressions). These equations can then be numerically solved from the disk mid plane (x=0) to infinity using a Burlish-Stoer method for stiff equations. The resolution is done in the same way as in <cit.>. The integration starts at x=0 with a guess for the parameters (μ, p). This guess of parameters allows us to define the vertical boundary conditions at the disk mid-plane, for example the toroidal and radial current as well as the accretion speed (u_r(x=0)<0). To be consistent with the assumption of bipolar magnetic topology we need to choose B_r(x=0)=B_ϕ(x=0)=0. The initial conditions for all fields are explicitly defined in Appendix <ref>. 

After properly setting up the boundary conditions the integration is propagated upwards using the resistive viscous MHD equations. As we move upward, the accretion flow is deviated and becomes parallel to the poloidal magnetic field. When this is achieved with enough accuracy, we switch to ideal MHD equations. In this regime, magnetic forces are more effective and try to accelerate the flow up to a super-slow magnetosonic (SM) speed. 

Fulfilling the regularity condition is not necessarily achieved for our initial choice of parameters. This condition is going to constrain the magnetization μ for a given value p. If μ is too large the flow will be accelerated too efficiently, which results in a shock. If μ is too small the acceleration will not be efficient enough and the flow falls back to the disk. By fine-tuning the value of μ, one can approach the critical point enough to safely make a leapfrog. This is done by conserving the various MHD invariants (see their definition in section <ref>). Once super-SM, the flow is still accelerated by the magnetic force and needs to become super-Alfvénic (A). 

This condition is going to constrain the parameter p, in the same way as before. If p is too small, the magnetic tension overcomes the centrifugal push and the magnetic surface closes back to the axis (B_r → 0). If p is too large, centrifugal acceleration is now too efficient leading to a vanishing toroidal field. By fine-tuning the parameter p, one can approach the Alfvén critical point close enough to jump beyond it and propagate the solution farther out (again by conserving the MHD invariants). It is important to note that a full super-SM solution needs to be computed from the origin (hence a new critical μ found) each time p is changed. This can be computationally demanding when the size of the explored parameter space is considerable.


 


While previously published JED solutions were found for a magnetization μ∈[0.1;0.8], we now wish to reproduce the results of global simulations and achieve super-A jets with magnetization values as low as 10^-4. In this regime, MRI is active and should be the source of the required MHD turbulence. However, around μ∼ 10^-4 and bellow, MRI dynamo becomes significant  <cit.>. Since such an effect is not included in our calculations, we restrict ourselves to solutions with a magnetization no smaller than μ∼ 5×10^-4 (see however <cit.> and <cit.> for the study of the influence of a mean field dynamo).

When seeking for the critical value of μ for obtaining super-SM flows, we realized that our previous numerical procedure was explicitly forbidding spatial oscillations in the magnetic field within the disk. The physical justification for this choice was that oscillating magnetic fields would give rise to an oscillating velocity field as well, leading most probably to an unstable situation (through e.g. Kelvin-Helmholtz instability). Since we were interested only in steady-state configurations, such solutions have been simply disregarded. But a careful look at global simulations (ie Fig.12 in <cit.> or Fig.6 in <cit.>) shows that this situation is actually realized, with B_r first becoming negative in the disk upper layers before becoming positive at higher altitude. To be consistent with these simulations, we thus relaxed our previous constraint and allowed now for negative radial fields within the disk. The fact that all previous JED solutions have been obtained only for μ>0.1 is a direct consequence of the explicit requirement (within our code) that the poloidal magnetic field has a monotonous vertical behavior within the resistive MHD disk zone. As will be shown in the next sections, relaxing this constraint (ie, removing any condition on B_r), allows for new solutions at much smaller magnetization levels.   Although we still recover the previous ones at near equipartition fields, we will mainly focus our attention on the new ones. It turns out that the asymptotic behavior of the associated jets is not different than that described in <cit.> and <cit.>, namely with a recollimation towards the axis. Our main interest will therefore be on the disk physics unveiled by these new solutions.   





§ SUPER-SM FLOWS
 




 §.§ The super-SM parameter space




Figure <ref> shows the parameter space for our fiducial case, obtained without any restriction imposed neither on p nor on μ. Each point corresponds to a flow that goes smoothly from the resistive MHD disk to the ideal MHD flow regime and becomes super-SM.  We have been able to extend the parameter space in μ by 4 orders of magnitude. The new enlarged parameter space entails several features:


  * The existence of distinct and well separated islands, namely zones in the plane p-μ where solutions can be found. The old parameter space found in <cit.> corresponds to the top island (with their ℛ_m=p/ϵ).



  * A monotonous behavior p(ξ), as discussed in <cit.>. This is mostly due to the vertical disk balance which leads to smaller ejected mass (ξ) when p is increased. Although a deeper examination shows that the rate at which p(ξ) varies with ξ depends also on μ and the turbulent parameters, we will not discuss further this already known trend.  


  * For a given ξ, the function μ(p) is bi-valued in some islands for μ<0.1. This is a signature of two distinct vertical equilibria. One branch is associated with a dominant toroidal field at small μ, whereas the other has a dominant radial field at larger μ. This will be further discussed in section <ref>.


  * The range [ξ_min,ξ_max] of possible super-SM solutions varies with the magnetization: both ξ_max and ξ_min increase at small magnetization. This will also be further detailed in section <ref>. Note that we restrained ourselves to values ξ≤ 1 because isothermal flows with ξ >1 would have a negative energy and could not describe unbounded flows <cit.>. Putting aside this issue, super-SM flows could nevertheless be achieved with larger ξ (up to 2 or more, <cit.>).   


Figure <ref> illustrates also the reason why solutions at small magnetization were difficult to find. The existence of forbidden zones in μ between islands (in particular at high magnetization levels) introduces a discontinuity that requires to jump to much smaller values in μ for a given ξ. In this forbidden zone, magnetic fields exhibit spatial oscillations without allowing for super-SM flows. This discontinuity made it hard to believe that the parameter space continued beyond what was already explored. 






In order to better understand these new super-SM solutions, we plot the vertical profiles of several quantities as function of the variable s=z/h along a magnetic surface for two solutions (Fig <ref>). They have the same parameter set (α_m=1, χ_m=1, 𝒫_m=1, ϵ=0.1), the same ejection index ξ=0.1 but were obtained with μ=6.7×10^-2 (top row) and μ=5.7×10^-3 (bottom). While the former exhibits one spatial oscillation in B_r, the second has 3 spatial oscillations. The different islands correspond therefore to different spatial oscillation modes in the radial magnetic field, starting with n=0 (no oscillation) for μ>0.1. For example the island located at 3×10^-2<μ<0.1 corresponds to n=1 spatial oscillation while 3×10^-3<μ<8×10^-3 corresponds to n=3 spatial oscillations. Furthermore, we can see that for decreasing μ both the extent of the islands, δμ, and the distance between the islands, Δμ, get smaller and smaller.  

This oscillatory behavior is also visible in all the other disk quantities namely, the other components of the magnetic field, the density and the velocities. These spatial oscillations start above the disk mid-plane, they seem to exhibit a constant wavelength and decay very rapidly before the SM critical point (shown as a vertical blue line). These spatial oscillations are therefore localized at the disk surface, overriding the resistive and ideal MHD regimes.  Putting aside the spatial oscillations, these solutions behave like the previous ones. A fraction of the disk mass (controlled by ξ) is deflected vertically by the combined effect of the thermal and magnetic (toroidal) pressure gradients and is ejected vertically. Within our isothermal situation with ϵ=0.1, the flow is energetically "cold" and the <cit.> criterion for cold ejection applies. As a consequence, the poloidal field lines are indeed bent by more than 30 degrees with respect to the vertical axis. This bending is actually more pronounced for these oscillating solutions, especially since magnetic compression is less of a danger for the disk vertical equilibrium at low magnetization. As analyzed in Appendix <ref>, this larger initial bending allows solutions to meet the Alfvén point at smaller altitudes. Once they become super-A, nothing seems to distinguish solutions at large n from the previously published n=0 solutions.    




 §.§ MRI-driven magnetic winds





The oscillatory solutions obtained at low magnetization are actually a generalization of the "exotic solutions" of <cit.>. 
The existence of spatial oscillations in all quantities, occurring at low magnetization levels, is a manifestation of saturated MRI-like modes or channel flows, as described for instance in <cit.>. However, a word of caution is appropriate. 
The solutions found in our work are not linear unstable MRI modes as the equations we solve are stationary. They are exact non-linear solutions of the MHD equations, which exhibit physics similar to that of MRI modes in the disc. This is not surprising since MRI modes are known to spontaneously saturate into wind-like solutions <cit.>. 
Nevertheless, as will be shown below, a linear approach allows to grasp the complex non-linear physics.
Since MRI is an ideal MHD instability, these modes will tend to develop only when the Alfvén time scale L/V_Az becomes smaller than the diffusion time scale L^2/ν_m  over a length L. A second necessary condition is that there is enough room to allow for spatial oscillations on that scale. This requires that the fastest growing MRI mode has a wavelength λMRI∼ L/n, where n is the number of spatial oscillations. The fact that the spatial oscillations need to stop when the radial magnetic field is positive imposes an integer number of spatial oscillations. Taking L ∼ h as an estimate of the relevant vertical dynamical scale provides the following crude conditions for the appearance of spatial oscillations 
    Rmag≡hV_Az/ν_m > 1
    λMRI/h≡ 2B_z/B_0√(μρ_0/ρ)∼1/n

We checked that these two conditions are indeed always verified in our solutions. The disk mid plane is always too diffusive for α_m > 1 and no spatial oscillations are present regardless of the magnetization μ. But since both the density and diffusivity decrease vertically, Rmag becomes large enough and spatial oscillations can then develop on a length scale of order h. Indeed, spatial oscillations are possible at the disk surface (z≃ h) since the flow is already in ideal MHD. Hence, at the disk surface the channel mode kicks in and produces the oscillatory behavior, as can be observed in figure <ref>.
The expression of λMRI harbors several features: 
 


  * The number of possible spatial oscillations depends mostly on μ, the smaller μ the larger n. However, this is only an order of magnitude estimate and there is some interval δμ around an average value μ_n allowing for super-SM solutions. This can be done by playing with the vertical profile of the density, namely the toroidal current parameter p and the ejection index ξ.     


  * As the density decreases, the wavelength increases and becomes eventually larger than the local dynamical scale. At the SM point, V_A >> C_s and the plasma beta writes β_SM≃ V_Az^2/2 V_SM^2 ≪ 1, showing that the magnetic tension becomes too large and quenches the instability. Spatial oscillations are therefore limited between the disk upper layers and the SM point, located a few disk scale heights (x_SM∼ 2 - 3). 



The above properties explain the existence of the islands seen in Fig. <ref>, as well as the fact that their spacing in μ decreases with μ. Indeed, the MRI wavelength and the SM point can be approximately related by n λMRI∼ h x_SM which leads to μ_n=f(x_SM)/n^2α_m^2, where f(x_SM) is a complicated function of the altitude x_SM  and μ_n is the value of μ for a given n. Since x_SM is weakly dependent of the magnetization μ (it is mostly related to the diffusivity scale height), this condition writes

    μ_n/μ_n+1∼n+1n^2

In order to test the generality of this expression, we took three different MHD solutions obtained with 𝒫_m=1, ϵ=0.1 but with different values for (α_m, χ_m). For a fixed value of the parameter p, we varied the disk ejection efficiency ξ and computed super-SM solutions, spanning thereby the various islands in μ. We could therefore associate the number of spatial oscillations n to a precise value μ_n, which is the critical value required to get a super-SM outflow. We could thus obtain the ratio μ_n/μ_n+1 as function of n, as shown in Figure <ref>. This plot demonstrates that the above simple analytical estimate is actually accurate and provides further evidence that the origin of the spatial oscillations it is indeed a saturated vertically stratified MRI-like mode, or channel mode, in an unbounded flow.  It explains also why the islands become closer and thinner as μ decreases (i.e. Δμ= μ_n - μ_n+1 decreases). The size of the islands becomes very small when μ approaches 10^-4. 
 
 

 
Since <cit.> it is well known that cold (negligible jet enthalpy) ejection requires at the disk surface a radial magnetic field component comparable to or larger than the vertical field. For near equipartition solutions with n=0, the generation of the magnetic field geometry at the disk surface is a natural outcome of the transfer of the disk angular momentum to the jet base, as described in <cit.>. This process can also be seen within the framework of MRI: as μ increases so does the MRI wavelength, resulting in the transport of the disk angular momentum along the field lines and fueling the jets <cit.>. There is therefore a possible continuous transition when the disk is thread by a large scale vertical field: from a turbulent (mostly radial) angular momentum transport at low μ to a laminar (mostly vertical) transport at large μ. Near equipartition field solutions with μ≥ 0.1 give rise to mostly centrifugally-driven super-A outflows <cit.>. This does not hold anymore at low magnetization levels and the existence of these channel modes in the disk upper layers is of utmost importance for cold ejection. Indeed, it is the existence of these spatial oscillations that actually leads to the generation of the radial and toroidal magnetic field components required for ejection. In other words, MRI-like spatial oscillations provide the bending necessary for ejection. This will be further discussed in Section <ref>.

The spatial oscillations exhibited by our solutions may seem in contradiction with the resistive profile used. Indeed, parasitic instabilities (such as Kelvin-Helmholtz) may be triggered and lead to a disruption of these channel modes (see for instance <cit.>). It is therefore possible that our solutions will ultimately lead to some internal rearrangement, namely a modification of the turbulent profiles. However, only numerical simulations could elucidate this issue. This will be further discussed in sec <ref>.


To summarize, channel modes triggered above the disk mid plane build up the magnetic field components allowing the ejection of cold disk material in the form of a super-SM ideal MHD flow. These oscillating modes are confined on a few disk scale heights, between Rmag>1 and λMRI/h≤ 1. But these modes are also constrained by the imposed boundary conditions, both at the SM point and at the disk mid plane. While the conditions at the SM point are quite general  (B_r and B_ϕ must be respectively positive and negative), those imposed at the disk mid plane are questionable. Our self-similar solutions have been computed assuming a symmetric magnetic structure such that at z=0 (i) B_r= B_ϕ=0 and (ii) u_r <0 (inward accretion motion). This clearly forbids other modes like for instance those leading to an outward decretion motion at the disk midplane. Allowing for such a boundary condition would lead to a supplementary half wavelength for instance (n+ 1/2)λMRI∼ h x_SM). More importantly, breaking the z-symmetry could also allow other modes with B_ϕ=0 located above or below z=0, while not significantly modifying the physics involved.  

Computing such solutions is beyond the scope of the present paper. We nevertheless argue that the parameter space shown in Fig. <ref> is actually a subset of the real parameter space of super-SM accretion-ejection structures. Indeed, since these new solutions would make use of the same type of modes, we believe that they would simply fill-in the forbidden zones between the islands. This property will be used in Section <ref>, where the effects of the turbulence parameters on the parameter space will be explored. 






§ FROM SUPER-SM TO SUPER-A FLOWS
 





 §.§ MHD jet invariants




Not all of the super-SM solutions shown in Fig. <ref> lead to steady-state outflows. In order to achieve that, they need to become super-Alfvénic (super-A) as well. The theory of steady-state MHD jets makes use of the existence of MHD invariants defined in ideal MHD along each magnetic surface of flux a. An axisymmetric, isothermal magnetic surface requires 6 boundary conditions at the base and features 4 integrals of motion (not counting the temperature): the two remaining quantities are thus determined by the SM and A regularity conditions.   

In ideal MHD mass conservation Eq.(<ref>) becomes

    = η(a)/μ_0ρ

where η(a) is the first invariant and describes the mass to magnetic flux ratio (η= μ_o ∂Ṁ_j/∂Φ, where Ṁ_j is the mass flux in one jet and Φ its magnetic flux). 

The induction equation (<ref>) writes

    Ω_*(a) = Ω -ηB_ϕ/μ_0ρ r

where Ω_* is the rotation of the magnetic surface. Since the field lines are anchored on the accretion disk, they rotate 
at roughly the same rate as the disk material at the jet base, namely Ω_*≃Ω_SM. 



The disk angular momentum conservation Eq.(<ref>) becomes

    L=Ω_*r_A^2 = Ω r^2-r B_ϕ/η

where L is the total specific angular momentum carried away by both matter and the magnetic field and r_A is the cylindrical radius where the flow becomes super-A.  

Finally, the projection along the magnetic surface of Eq.(<ref>) leads to 

    E(a) = u^2/2+H+Φ_G-Ω_*rB_ϕ/η

where E(a) is the Bernoulli integral and describes the total specific energy carried away along the magnetic surface. The enthalpy, defined as H = p/ρ namely H= C_s^2 lnρ for isothermal flows, can be safely neglected in analytical estimates (since -H/Φ_G ∝ϵ^2). Note however that our numerical resolution solves the full set of MHD equations including all terms (see Appendix <ref>).

Since our super-SM solutions are in ideal MHD regime, the MHD invariants are already determined. It is therefore convenient to express them as function of the underlying disk parameters. To do so, we normalize these 4 invariants by quantities defined at the anchoring radius r_o of the magnetic surface, at the disk equatorial plane. This leads to the following dimensionless parameters for cold jets launched from thin accretion disks[To derive the expression of κ, mass conservation Eq.(<ref>) is written as d Ṁ_a/dr = 2 d Ṁ_j/dr, leading to the useful relation ρ^+ u_z^+ ≃ξϵρ_o u_o valid at the disk surface, while Eq.(<ref>) gives u_o= α_m p V_A.]     


    ω   ≡Ω_*/Ω_Ko
    κ   ≡ηΩ_Kor_o/B_o≃ξα_m p/μ^1/2
    λ   ≡Ω_* r_A^2/Ω_Ko r_o^2≃ 1 + 1/2ξΛ̅/ 1+ Λ̅
    
        e   ≡E/Ω_Ko^2 r_o^2/2≃ 2λ - 3 + Θ

where Ω_Ko= √(GM/r_o^3) is the Keplerian angular velocity and Θ=2H_SM/Ω_Ko^2 r_o^2 is the normalized enthalpy at the SM point. In geometrically thin disks, ω is always close to unity, but a significant deviation may occur for thicker disks <cit.>. For simplicity, we assumed ω=1 in the expression of the dimensionless specific energy e (see its exact expression in Appendix <ref>).

The term Λ̅= M_zϕ/M_rϕ is the ratio of the vertical (jet) torque 

    M_zϕ =  2 r B_z B_ϕ/μ_0|_SM

exerted at the disk surfaces to the total radial torque 

    M_rϕ = 2_0^z_SM1/rr[r^2(T_rϕ+B_rB_ϕ/μ_0)] z

acting within the disk. The latter torque includes thereby both the jet (laminar) and turbulent (viscous) contributions to the radial transport of angular momentum. The disk angular momentum conservation equation (<ref>) writes M_zϕ + M_rϕ = 2_0^z_SMρ𝐮_𝐩·Ω r^2 z≃ - Ṁ_a Ω_K/4. Using the definitions of the magnetic lever arm and ξ leads then to Eq.(<ref>). Note that it is a generalization of the relation found in <cit.>, where the radial transport of angular momentum by the laminar torque was negligible (parameter Λ in their Eq. 33).    
 
 
 


The mass load κ and magnetic lever arm parameter λ have been first introduced by  <cit.> and extensively used since then. They are related by λ≃ 1 + | B_ϕ/B_o |_SM/κ <cit.>. The value of the toroidal magnetic field depends on the underlying disk physics and is discussed next section.  In any case, the outflow must have a positive energy e >0, which requires λ > 3/2 when ω≃ 1 and  for cold flows with Θ≪1, namely when thermal effects (such as a warm corona) can be neglected.   





 §.§ The magnetic shear as a function of the disk magnetization μ





Using Eq.(<ref>) and (<ref>), we can easily compute the jet invariants for all our super-SM solutions and put them in the classical κ-λ plane of super-A solutions. This is shown in Fig. <ref>. The islands seen in Fig. <ref> appear as well here, with a clear trend in κ: the smaller μ the larger κ. The super-SM constraint, that determines μ(p) for a given ξ, allows to obtain almost all possible values of ξ up to 1, with κ∝ξ. But reaching large values of κ can only be done by switching to another island. The jet mass load is thus a function κ (ξ,μ) with an approximate linear dependence on ξ. For a given mass loss ξ, increasing κ can be done by decreasing the magnetic field strength μ. Matter dominated super-SM flows with κ >1 become thus achievable.   


On the other hand, the magnetic lever arm λ does not appear to be strongly dependent on μ but mostly on ξ. This can be clearly seen in Fig. (<ref>), where λ is plotted as function of ξ for all our super-SM solutions. Indeed, λ = 1 + R/2ξ, where R= Λ̅/(1+ Λ̅) is a rather weak function of μ and ξ (the small dispersion in λ doesn't seem to depend on μ). This is remarkable as ξ and μ span respectively 2.5 and 4 decades. Such a behavior must therefore be the outcome of some intrinsic physics. Using this result and Eq.(<ref>) leads to the necessary constraint on the  magnetic shear at the disk surface 

    |B_ϕ/B_0|_SM = p α_m/2 μ^1/2 R  ∝μ^-1/2

since R is a weakly varying function and p has a small range. 



This scaling of the toroidal magnetic field can be understood the following way. For a rather wide range in disk conditions, the flow must become super-SM near the surface, namely u_z^+ ∼ V_SM∼ C_s V_Az/V_A where V_A is the total Alfvén speed. The vertical velocity u_z^+ is provided by the unbalance in the vertical forces around the dis surface, which is quite difficult to estimate. Another way to grasp it is to look at the Ohm's law (Eq. <ref>) at the turning point where the radial velocity vanishes, right before the SM point. At this particular point u_z^+ B_r^+= ν_m^+ ∂ B_r/∂ z, which provides the scaling u_z^+ ∼ν_m^+/h. This simple relation tells us that mass loading in jets is a diffusion process and that the initial jet velocity is directly related to the strength of the poloidal magnetic diffusion. Using now u_z^+ ∼ V_SM leads to 

    α_m^2μ≃ F^2_SM1/1+B_r,SMB_0^2+B_ϕ,SMB_0^2

where F^2_SM is a function depending on the vertical profiles of the temperature and the magnetic diffusivity. It stems from this expression that, in order for this condition to remain valid at all μ, the magnetic shear |B_ϕ,SM|/B_0 must indeed scale as μ^-1/2. It is therefore the SM constraint itself that imposes such a scaling: it guarantees that, whatever μ, cold super-SM solutions can be found.

We can use the SM constraint on the magnetic shear, Eq.(<ref>), to derive an approximate expression for the vertical torque, Eq.(<ref>):

    M_zϕ/r P_o≃μ|B_ϕ/B_0|_SM≃κ (λ-1)μ∝μ^1/2

this scaling is consistent with Figure <ref>. Therefore, the wind stress can be easily modeled as a function of the magnetization and the plasma pressure at the disk mid-plane. This prescription could be useful for including the effects of wind driven accretion in hydrodynamic models.





 §.§ The super-A parameter space






As discussed previously, steady-state solutions are only those that have the capability to produce super-Alfvénic flows. Magnetic acceleration can be seen as some centrifugal effect, the frozen-in jet plasma being accelerated because magnetic field lines are rotating faster than the jet material. This can be illustrated using Eq.(<ref>,<ref>), leading to Ω = Ω_*(1-g) where

    g=m^2/m^2-1(1-r_A^2/r^2)

with r_A the Alfvén radius and m= u_p/V_Ap the poloidal Alfvén Mach number. The function g measures the discrepancy between the two angular velocities and is related to the poloidal current flowing in the jet <cit.>. Starting from a tiny value at the disk surface, this function increases as the flow gets accelerated. It can then be seen that when the flow becomes super-A, namely m=1, a regularity condition r=r_A must be fulfilled. 




The Bernoulli equation can be interpreted as providing the velocity that matter reaches for a given "magnetic funnel". At infinity, if acceleration is so efficient that the magnetic energy becomes negligible, Eq. (<ref>) gives the maximum jet velocity u_p∞≃√(2 E)≃Ω_Kor_o √(2 λ-3) for cold flows. But the shape r(z) of this magnetic funnel, or more precisely the jet transverse equilibrium, is provided by the Grad-Shafranov equation

    ·[( m ^  2  - 1 )  a /μ _  o  r ^  2 ] = ρ{ E  - ΩΩ _  *  r _  A  ^  2  + ( Ω r ^  2  - Ω _  *  r _  A  ^  2 ) Ω _  * }
    +  B _ ϕ ^  2  + m ^  2  B _  p  ^  2 /μ _  o lnη

where d/da = a·/a·a <cit.>. This equation is obtained by projecting Eq. (<ref>) across the magnetic surfaces and, in practice, is not used to solve the jet physics. It does however provide a useful tool to derive the analytical constraint to get trans-Alfvénic flows (see Appendix <ref>). 
 

A lengthy calculation allows to express g_A explicitly as function of the MHD invariants and the Alfvén position angle Ψ_A. There are two branches that can be both connected to the accretion disk, one of them being chosen according to the interplay between λ and  κ.  <cit.> pointed out that, for a given mass load κ, there must be a minimum λmin that they estimated such that κλmin (2λmin- 3)^1/2 = 1. The generalization of this estimate is that, for a given magnetic lever arm λ, there must be a minimum mass load κmin for a cold flow such that  

    κmin^2λ^3g_B^2= 1   with   g_B^2 = 1 - 3/λ+   2 sinΨ_A/λ^3/2

(Eq. (<ref>) in Appendix <ref>, with ω=1). This expression provides the two solid blue curves shown in Fig. <ref>, computed in two extreme cases for the location of the Alfvén surface: near the SM surface with x_SM=2, namely  Ψ_A=/2-arctan(2h/r) (lower curve) and much further out with Ψ_A=/3 for the upper one (as in typical near-equipartition cold flows). For a given magnetic lever arm λ (mostly determined by ξ), κ must be large enough. This rules out all super-SM solutions located at the left-hand side of these curves.   

It can be seen right away that the parameter space of super-A solutions for near-equipartition fields (n=0) will be very small, with mass loads κ smaller than 0.1, disk ejection efficiencies ξ smaller than 0.08 and magnetic lever arms larger than ∼ 7, in agreement with <cit.>. On the other hand, MRI-like driven flows at small μ seem to allow mass loads larger than unity with small λ. These isothermal flows could be of great interest for reproducing dense outflows (ξ >0.1) with low asymptotic speeds (u_p∞∼Ω_Ko r_o). In other words a weakly magnetized disk is better suited for a a massive ejection than a near-equipartition disk, not only by providing the necessary bending at the disk surface through the MRI-like mechanism, but also by bringing the Alfvén surface closer to the disk (see discussion on the role of Ψ_A in Appendix <ref>).  
    
Nevertheless, not all super-SM solutions located at the right hand side of the solid curves in Fig. <ref> can become super-A. For a given ξ, λ is roughly determined and so is the cylindrical radius r_A=r_o (λ/ω)^1/2 of the Alfvén point. But Eq. (<ref>) shows that κ depends also on p and μ. For a given toroidal current density p, the SM constraint provides μ so that κ is fixed. On the other hand, p determines also  the radial magnetic field component at the disk surface and thereby the initial jet bending (see fig. <ref>). As a consequence, playing with p not only affects the disk vertical equilibrium (SM point) but also this initial jet angle. Not all couples (κ, λ) fulfill the Grad-Shafranov equation (<ref>). If a solution is not possible, this means that there is no altitude z_A of the Alfvén point that can be found starting from the conditions provided at the base of the jet (SM point). Changing the value of p leads to a slight modification of κ, λ as well as the jet angle at the SM point, allowing thereby to (possibly) meet the A condition. This translates into an adaptation of the altitude z_A of the Alfvén point (thus the angle Ψ_A) according to the disk conditions. 

The necessary condition g_B^2>0 highlights this aspect. When λ is large, r_A/r_o is large and gravity plays a negligible role so that g_B^2>0 is satisfied whenever λ > 3 (dashed line in Fig. <ref>). But gravity cannot be neglected anymore at lower values of λ, as can be seen in Eq. (<ref>). In that case, providing g_B^2>0 requires to increase sinΨ_A, namely to bring the Alfvén surface closer to the disk surface. The closer it is to the disk and the less energy is been consumed to reach it. We have been able to find super-A solutions with λ as small as 1.6 from low magnetized accretion disks.  


Figure <ref> shows the parameter space ξ(μ) for  isothermal super-A solutions in our fiducial case. It is a subset of 
the SM parameter space shown Fig. <ref>. The islands are now seen as almost vertical stripes in μ with a range in ejection index ξ. We recover the same results as <cit.> for near equipartition fields (n=0, right) but with a significant enlargement in μ by almost 4 orders of magnitude (n=8, left). 

The color scale indicates the ratio of the radial to the toroidal magnetic field components at the SM point. While n=0 solutions are clearly dominated by the radial component, the toroidal field becomes gradually dominant as n increases (μ decreases). This is of course consistent with the scaling |B_ϕ/B_o|∝μ^-1/2 imposed by the SM regularity condition. However, it highlights a possible dichotomy between "magnetic tower" jets <cit.>, where ejection is due to a dominant B_ϕ field, and "centrifugally-driven" jets <cit.>, where a dominant radial field is of utmost importance. As already pointed out in <cit.>, these are two expressions of the same magnetic process. However, the dependence ξ(μ) is quite different for the two cases within each island. It can be seen for instance that for n=0 the ejection index ξ increases when μ increases (although in a very limited range and for ξ < 0.08). On the other hand, above n=3 (μ < 10^-2 for our fiducial case), it is the other way around: ξ decreases for increasing μ (although for 0.4 > ξ > 0.08). The functional dependence ξ(μ) can thus be seen as a fingerprint of the dominant ejection mode. 

As discussed earlier, we expect to find other solutions by changing the boundary conditions at the disk mid plane. These solutions should be  located between the islands appearing in Fig. <ref>, possibly filling-in the actual forbidden area. However, and for the same reasons, the general contour of the parameter space (shown as the solid line) should not be modified. The upper and lower contours describe the curves ξmax(μ) and ξmin(μ) respectively.      


The minimum value of the ejection index ξmin(μ) increases when the magnetization decreases. This a fossil feature of the SM constraint (Fig. <ref>). Indeed, when μ decreases the increasing toroidal magnetic field (Eq. <ref>) leads to a stronger vertical push and to a larger quantity of ejected plasma. Furthermore, as μ decreases, the number n of spatial oscillations before the SM point increases as well, enforcing thereby the disk to be subjected to them at deeper and more massive layers.
The terminal velocity of the outflow is going to be linked to the value of λ, which is a function of mostly ξ. Thus, the maximal terminal velocity will be determined by ξmin, which depends mostly on the disk magnetization μ. 


As shown in <cit.>, the maximum ejection index ξmax is determined by the Alfvénic constraint. It is interesting to see that it has a non-monotonous behavior, first increasing with μ until the maximum value ξmax=0.35  for n=2, and then decreasing down to ξmax=0.08 for n=0. As illustrated in  Fig. <ref>, increasing ξ leads to a decrease in λ. Low values of λ are possible only for large values of κ, which are accessible only by decreasing μ. As a consequence, the Alfvén surface comes closer to the disk (z_A/r_A decreases as ξ increases, see Fig.<ref>). However, if we keep decreasing μ, the magnetic energy available in the jet becomes also smaller and jet acceleration less efficient. The Alfvén surface moves away from the disk (z_A/r_A increases as ξ increases), requiring thereby a larger magnetic lever arm λ to get super-A flows so that ξmax decreases. 






§ GENERAL ACCRETION-EJECTION PROPERTIES




We showed that, for our fiducial parameter set, there is a maximum value ξ_max= 0.35, obtained with n=2 for μ∼ 10^-2, while the minimum value ξ_min=5× 10^-3 is obtained with n=0 for μ= 0.5. In this section we explore the effect of the turbulence parameter α_m on the existence of cold super-A flows. More specifically, we investigate how turbulence affects the contours of the parameter space, namely the curves ξmin(μ) and ξmax(μ). To avoid confusion, we will not plot the points corresponding to each solutions found, but mark only the contours of the parameter spaces.   

The exploration of the anisotropy parameter χ_m is done in Appendix <ref>. While the turbulence level parameter α_m affects all magnetic field components, χ_m affects only the toroidal field. Besides, as will be shown, the condition for jet launching from a thin disk introduces the extra link  χ_m ∼α_m^2, that can also be written as ν'_m ∼α_m^-1 V_A h. We thus focus here only on the effect of α_m.   
   




 §.§ Effect of the turbulence level α_m




Figure <ref> shows that α_m has a huge impact on the parameter space of super-A outflows, obtained here with ϵ=0.1, χ_m=1,  P_m=1. Two important trends arise with α_m:


    
  * When α_m increases above unity, ξmax is barely modified while the curve ξmin increases, leading to a shrinking of the parameter space. Furthermore, as α_m increases, solutions of same ξ are displaced to smaller μ. Note that we did not explore values larger than α_m=8, as it corresponds to the scaling deduced from shearing box simulations <cit.>.

    
  * When α_m decreases below unity most of our super-A solutions disappear. For α_m=0.8 only two solutions are found,  one with n=0 and the other with n=3. The fact that no solution can be found for α_m smaller than unity has been already reported in <cit.>. 


The displacement to smaller magnetizations as α_m increases arises naturally from the SM constraint (Eq. <ref>). Since the velocity at the disk surface is related to the poloidal diffusion, increasing α_m requires to decrease μ (so that μα_m^2 remains approximately constant). The drastic diminution of the parameter space as α_m varies is related to the jet launching condition and requires a deeper examination.  

Magnetic ejection occurs only if the jet torque F_ϕ= J_z B_r - J_r B_z ∼ - J_rB_z switches sign and becomes positive at the disk surface <cit.>. This requires therefore that J_r decreases on a disk scale height, which can be guaranteed only if

    Γ = 3/2χ_m/α_m^2p/p-𝒫_mϵ≲1

where Γ controls the vertical scale of the emf in the induction equation (see Eq. (<ref>) and Appendix <ref> for more details). This condition implies that for solutions with p of order unity, the toroidal current density must adapt to the turbulent properties of the disk with p ∼χ_m /α_m^2 (in turn, this also implies χ_m ∼α_m^2). Thus, when α_m increases p needs to decrease. Since p controls the radial magnetic field at the disk surface (see fig. <ref>), a decrease in p leads to a decrease of the magnetic vertical compression, thus to a larger mass loss rate from the disk. This feedback on the disk vertical balance explains why the curve ξmin increases when α_m increases.  

In order to keep Γ near unity, as α_m increases one gets p → pmin= P_m ϵ which could be very small. This might be an indication that, for larger values of α_m, the MHD solution would eventually try to reverse the sign of the accretion speed, with an outward motion at the disk mid plane (u_r >0 and J_ϕ <0). Such a situation, seen in numerical simulations, is actually forbidden by our assumed boundary condition. This is a general symptom that was discussed in section <ref> and will also be touched upon in section <ref>. On the contrary when α_m decreases Γ quickly becomes larger than unity since p/p-𝒫_mϵ is bounded by one. Hence, toroidal field induction becomes highly inefficient and the torque F_ϕ remains negative, providing no magnetic acceleration. As a consequence, solutions are mostly wiped out when α_m becomes smaller than unity. 







 §.§ From jets to winds






Despite the existence of biases introduced by our prescription of MHD turbulence, we would like to summarize here how the disk magnetization μ affects some properties such as (1) the disk ejection efficiency, (2) the accretion Mach number, (3) the angular momentum transport and global energy budget and (4) some jet/wind properties.   


The disk ejection efficiency ξ is one of the most important quantities in Jet Emitting Disks (JED) as it provides the link between disk and jet properties. It is defined as Ṁ_a ∝ r^ξ, so it can be measured for instance in numerical simulations that converge to power-law accretion disks (assuming they reach a steady state). Since the local disk magnetization μ is also easily measurable, simulations can be directly compared to our parameter space ξ(μ). However, as shown above, one needs to make sure that the turbulence parameters α_m, χ_m,  P_m are comparable, as they do influence considerably the final outcome. Note also that our parameter space has a maximal size for α_m ∼χ_m ∼ P_m ∼ 1. For these values, ν_v ∼ν_m ∼ν'_m ∼ V_A h and both modes of ejection are at play, allowing thereby to maximally explore the parameter space.

Also, our calculations have been done for cold outflows (isothermal magnetic surfaces) only and it is known that thermal effects may drastically enhance ξ <cit.>. While the largest value of ξ is imposed by the jet physics (Alfvén constraint), the smallest value depends on the disk physics (vertical equilibrium and toroidal field induction). For cold jets, we found a minimum ξmin≃ 5  10^-3 and a maximum ξmax≃ 0.47, with a clear tendency of reaching larger ξ with low-μ solutions. Contrary to previous near-equipartition solutions, cold massive outflows with ξ∼ 0.1 are possible as long as spatial oscillations are allowed. But cold super-A solutions with ξ > 0.5 remain out of reach.   


The accretion Mach number is defined at the disk equatorial plane as m_s= - u_r/C_s, where u_r is the radial accretion speed due to the turbulent and jet torques. However, m_s is also related to turbulence by

    m_s = p α_m μ^1/2

through Ohm's law (Eq. <ref>). Given the small range in p, it becomes obvious that m_s ∝μ^1/2 whatever the dominant torque. No wonder then that supersonic accretion becomes possible only for near-equipartition μ>0.1 solutions. Such a high accretion speed has profound consequences: not only accretion time scales are much shorter than in usual standard accretion disks, but it may lead to optically thin accretion disks with observable features in young stellar objects <cit.> or X-ray Binaries <cit.>. One might however question if defining m_s at the disk equatorial plane remains relevant in the case of MRI-like active disks. Indeed, not only the radial velocity is prone to spatial oscillations, but we expect a larger accretion speed at higher altitude. We checked this by computing the density weighted Mach number within the disk,  

    m̃_s =  Ṁ_a/2 r C_s Σ   =  -1/C_s∫_0^x_SMu_rρ x/∫_0^x_SMρ x,

from the disk mid plane up to the SM point. For all solutions found, m̃_s is never larger than m_s by more than a factor 3. Looking at fig. (<ref>) this may seem surprising, but the spatial oscillations tend to compensate each other leading to very little difference with m_s. This important result confirms that only n=0 solutions provide supersonic accretion. A second interesting aspect is the influence of the turbulence strength α_m on m_s. While Eq. (<ref>) seems to imply that m_s increases with α_m, it behaves in the opposite way. This is because, as the diffusion of the poloidal field (α_m) increases, the toroidal current density J_ϕ decreases with p ∝ 1/α_m^2 and so m_s decreases.  


Accretion is due to the vertical and radial torques acting in conjunction. To better understand the disk angular momentum transport  it is interesting to look at the ratio Λ̅= M_zϕ/M_rϕ of the vertical (jet) torque to the radial (laminar and turbulent) torques acting on the disk (see Sect. <ref> for their expression). It can be analytically estimated by neglecting the radial laminar  contribution, leading to 

    Λ̅∼ . - 2 B_ϕ B_z/μ_o α_v ϵ P_o |_SM∼2 μ^1/2/α_m  P_m ϵ |  B_ϕ/B_z  |_SM∝2/ϵ P_m α_m

where the last expression comes from the SM constraint | B_ϕ/B_z| ∝μ^-1/2. It can be seen that Λ̅ is only weakly dependent on μ and decreases when α_m increases. This is illustrated in Fig. <ref>. Previous cold solutions at near-equipartition fields with α_m=χ_m=  P_m=1 had Λ̅∼ϵ^-1 = 10.  Here, as μ decreases Λ̅ decreases also because of the increasing effect of the radial torque due to the laminar magnetic field. Large spatial oscillations of the magnetic field within the disk lead to radial transport of the disk angular momentum within the resistive layers that will not be carried away by the jets (channel modes). However, because of the MRI scaling of α_v with μ (Eq. <ref>), Λ̅ remains larger than unity even at μ∼ 10^-3 for α_m=1. This is no more the case for α_m=8, where  solutions with Λ̅ as low as 0.3 can be obtained (Fig. <ref>, right). Further increasing α_m would thus allow to produce accretion-ejection structures with massive winds (ξ >0.1) that carry away a negligible fraction of the released accretion power. We should however remain cautious as our results depend on our assumed vertical profiles for the turbulent coefficients (see next section).  


Finally, there is a need to address the asymptotic properties of the new low-μ solutions, namely jet speed and collimation. In any case, the fraction of the initial energy that remains stored within the magnetic structure depends on the jet transverse equilibrium (Eq. <ref>). For cold jets, the maximum asymptotic speed is reached only if the magnetic structure converts almost all its energy into kinetic energy, so that u_p∞= √(2E). Even in that case, the value of the Bernoulli invariant E(a) is given by Eq. (<ref>) and depends on the magnetic lever arm parameter λ. For λ >3, the position of the Alfvén point has no impact and the asymptotic jet speed can be estimated with the usual expression √(2λ-3) (in units of the Keplerian speed at the anchoring radius r_o). But for smaller λ, the Alfvén surface is much closer to the disk which introduces a deviation from this expression, as shown in Fig. <ref>. The difficulty with  low-μ solutions is to relate λ to the disk ejection efficiency ξ, since it requires the knowledge of the vertical to radial torques ratio Λ̅ (see Eq. <ref>). If we take for instance our fiducial parameter set, Fig. (<ref>) shows that λ(ξ) deviates the most from the usual expression λ = 1 + 1/2ξ for ξ∼ 0.1, corresponding neatly to where Λ̅ is the smallest. It is thus more problematic to derive ξ directly from the jet asymptotic speed in the case of low-μ JEDs. 

A clear distinction between MHD winds and jets, which are both super-A MHD flows, can nevertheless be made. Indeed, λ determines the importance of the initial magnetic reservoir feeding the jets. Defining the initial jet magnetization σ as the ratio of the MHD poloidal Poynting flux  to the kinetic plus thermal (enthapy) energy flux (measured at the jet base taken as the SM point) leads to the general useful relation for cold flows

    σ =  . - Ω_* r B_ϕ B_p/(u^2/2+H) ρ u_p μ_o |_SM≃ 2ω (λ -1) ≃ω/ξΛ̅/1+Λ̅

Jets are characterized by high speeds (large λ) that can become farther out self-confined thanks to the dominant hoop-stress wrt to both plasma pressure gradient and centrifugal terms. Jets are therefore Poynting-flux dominated flows with σ >1. On the other hand, winds are low speed MHD flows with small λ, with almost no collimation besides that introduced by the external pressure. Winds are therefore matter-dominated flows with σ <1. 

We plotted in fig. (<ref>) the contour σ=1 computed using ω=1. All solutions obtained with α_m=1 (left) lie well above this contour and display σ >1. They are therefore representative of self-confined jets, whatever the disk magnetization μ. Figure (<ref>) shows the magnetic surfaces of the two low-μ solutions displayed in Fig. (<ref>). It can be seen that their asymptotic behavior is very similar to the high-μ solutions: the magnetic surfaces first widen before undergoing a recollimation towards the axis (see for instance Fig. 6 and discussion Sect. 5 in <cit.>). 

But fig. (<ref>) reveals also that some low-μ solutions obtained with α_m=8 do cross the σ=1 contour. These solutions have λ very close to the limiting value 3/2 and reach σ=1 because both ξ is quite large and Λ̅ is small, in agreement with Eq. (<ref>). Although these solutions do cross the Alfvén point, they meet soon after the modified fast magnetosonic (FM) surface (and stop). However, according to the analysis done in <cit.>, getting super-FM solutions requires to play with the jet energy equation, which is forbidden with isothermal flows. Moreover, isothermal flows from thin accretion disks are cold and there is no way to provide energy to the outflow when λ becomes too small. The only possibility to get winds is then to include relevant thermal effects acting already at the disk surface. Taking into account such a warm corona can be done following the method used in <cit.>, leading to the build up of a relevant Θ term in the Bernoulli equation (<ref>)). Our guess is therefore that warm low μ solutions obtained with α_m=8 will provide proper MHD wind solutions. This is postponed to future work.      






§ DISCUSSION





 §.§ Comparison with numerical simulations



The first paper showing the existence of super-A flows from low magnetized accretion disks was <cit.> and was then extended by a large numerical survey in μ done by <cit.>. Both works used ϵ=0.1 but while the former included viscosity with α_v=0.9 (all stress components), the latter neglected it with P_m=0, so that the disk angular momentum removal is only done by the jet torque. In terms of magnetic diffusivity, the former used χ_m=1 and α_m starting from 20 and increasing with the radius as μ decreases, whereas the latter used χ_m=2, α_m ≃ 3 and a constant μ across the accretion disk. Nevertheless, super-A jets were found in both works with μ as low as 10^-4, the physics of ejection following the description done here for MRI-like driven outflows. 

The dependences in μ of the MHD invariants κ, λ as well as the accretion Mach number m_s shown in <cit.> are also followed by our solutions, although viscosity has been neglected in their work. Indeed, m_s ∝μ^1/2 is an analytical result and for μ∼ 10^-3, they found jets with κ≃ 15, λ≃1.8, for μ∼ 10^-2 κ≃ 4, λ≃ 2 and for μ∼ 10^-1 κ≃ 1, λ≃ 3 (using the proper normalization κ=k̃/μ^1/2ϵ). These values are consistent with our Fig. <ref> (derived with α_m=1, χ_m=1,  P_m=1). Moreover, they find a similar scaling for |B_ϕ/B_r| which has lead them to deduce a dichotomy between the magnetic tower and magneto-centrifugal solutions, with a critical magnetization μ≃ 0.01 separating them. This is also consistent with our work, even though we find that a proper differentiation between solutions should also include the disk ejection efficiency ξ (our Fig. <ref>). 

The main discrepancy between our semi-analytical work and these numerical "alpha" simulations lies in the disk ejection efficiency. It can be derived knowing λ and the jet torque ratio Λ̅  (Eq.<ref>). Since the simulations of <cit.> have Λ̅/(1+Λ̅)=1, we obtain ξ∼ 0.62, 0.5, 0.25 for μ∼ 10^-3, 10^-2, 0.1 respectively (these values are consistent with ξ derived using the independent expression κ≃ξ m_s/μ). These values of ξ are larger than those obtained in our case. This cannot be an effect of the self-similarity since the simulations themselves do converge to such a situation (see for instance the conical shape of the critical surfaces in their Fig. 1). On the contrary, we believe that it may be due to the fact that these numerical simulations are not cold outflows. Using the same normalization as these authors, the Bernoulli invariant (<ref>) writes e= λ- 3/2 + Θ/2. Their Fig. 6 clearly shows that the initial thermal content Θ is non negligible and becomes even of the order unity for several simulations. Since it has been demonstrated that heat deposition at the disk upper layers leads to enhanced mass loss <cit.>, we argue that this is the main cause of the observed discrepancy in ξ (see also discussion on a possible numerical bias p9-10 in <cit.>). 

For that same reason, we cannot compare our cold solutions with the 3D numerical simulations done by <cit.>, as they assumed the existence of a heating term acting at the disk surface (leading to a huge ejection efficiency ξ∼ 1). Furthermore, non-ideal terms (Hall and ambipolar diffusion) have been assumed in the induction equation and it is not clear yet how this affects the generation of the toroidal magnetic field and thereby the whole ejection process.  

The 3D simulations of vertically isothermal, ideal MHD disks with a magnetization μ∼ 10^-3 done by <cit.> should in principle be comparable to our own work. Indeed, they obtain a MRI active accretion disk giving rise to a super-A flow launched from r_i=0.5 to r_e=5 with λ∼ 10-14 and carrying a mass fraction 2 Ṁ_j/Ṁ_a ∼ 0.4%≃ξln (r_e/r_i). These values correspond to a very small ejection efficiency ξ∼ 1.7-2  × 10^-3 which, since they report Λ̅∼ 5% only of the disk angular momentum transport due to the wind, is indeed consistent with Eq.(<ref>). Within our framework, such a small value of Λ̅ would require to increase α_m up to ∼ 15 (according to Eq. (<ref>), since the authors report α_v= 0.5 and P_m ∼ 1). However, Fig.(<ref>) shows that for α_m=8 we obtain even more massive solutions with  ξ>0.08. The reason of this discrepancy is probably due to the MHD turbulence itself. Indeed, our semi-analytical solutions are mathematically exact but depend on the vertical profiles that are assumed for all the transport coefficients. 

The work of <cit.> but also <cit.> show that global 3D simulations that include the wind torque have a turbulent diffusivity and radial angular momentum transport up to almost 10 pressure scale heights. This is quite surprising as the original consensus was that the disk (defined by its scale height) would be responsible for turbulence. It seems on the contrary that MHD turbulence, under the presence of a large scale B_z field at low magnetization levels, is able to persist and self-sustain high above the disk. As a consequence, the resistive disk survives at larger altitudes and ideal MHD ejection occurs only further up, decreasing thereby the disk mass loss ξ. This new feature can be easily incorporated within the self-similar framework through the use of different vertical profiles for ν_v, ν_m and ν'_m. This deserves definitely further investigation.  
    
<cit.> used local shearing box simulations to infer the angular momentum transfer rate due to a magnetized wind in the context of dwarf novae. The resulting prescription has then been used to predict the secular evolution of such a system <cit.>. <cit.> found that |B_ϕ/B_o| ≃ 30 (a constant) at the disc surface for all μ≳ 10^-4 while we find | B_ϕ/B_o|∝μ^-1/2. This difference could have several origins: while <cit.> defines the surface as the altitude where B_ϕ is maximum, we define the surface at the SM point. More importantly, the shearing box does not satisfy Eq. (<ref>) since the poloidal field is allowed to be advected radially <cit.>. This implies that shearing box solutions are not strictly speaking secular stationary solutions, but also that the scaling (<ref>) is not valid in a shearing box. This is yet an other illustration of the limitations of shearing box solutions to describe winds. This difference will likely have a significant impact on the dynamics of these systems as our scaling predicts a stronger wind at low magnetization compared to shearing box solutions. 

Finally, we show that it is possible to derive a scaling of the wind stress as a function of the disk pressure from our solutions. We find that M_zϕ∝μ^1/2 (Fig. <ref>), which can be of use in secular disc models to include the effect of a magnetic wind on the disk evolution.
 



 §.§ Caveats




This work is subject to biases arising from (1) stationarity  (2) our imposed geometry and (3) the prescriptions used for the MHD turbulence. We briefly discuss these important points below.  

The validity of the steady-state assumption depends on the time scales considered. In thin accretion disks, the local dynamical time scale is of the order of the Keplerian period. Another important time scale is the accretion time scale, which is 1/(ϵ m_s) times longer than the dynamical time. Such a huge difference in these two scales allows to make a simple ordering. Our steady-state solutions can thus be considered valid on scales that are longer than the local dynamical time but smaller than the accretion time. On this longer time scale, both mass and magnetic fields can be advected inward/diffused outward, leading to modifications of the radial profiles. This is illustrated for example in <cit.>, where the disk magnetization, μ, is seen to evolve on these long time scales. Note also that one might incorporate magnetic flux advection within a steady-state approach by modifying Eq. (<ref>) (see for example <cit.>). It is however unclear why such an effect should also follow a self-similar scaling.


What we call geometry covers actually different assumptions. The first obvious one is self-similarity and has been already discussed extensively in <cit.>. Clearly, jet asymptotics are influenced since neither the inner (jet axis) nor the outer (jet boundary) regions can be described within this mathematical formulation. However, jets should be nevertheless well approximated by these solutions when the critical surfaces are close to cones, a situation which arises whenever the jet emitting disk is established over a large radial extent (r_e>> r_i and Ṁ_a ∝ r^ξ). Note however that most of the mathematical relations between jet and disk parameters are general and can be used to interpret and understand steady-state 3D simulations. For that same reason, we expect our parameter space to be only weakly affected by self-similarity.  
 
The second aspect is the z-symmetry imposed on the bipolar magnetic field structure, even for the flux function a(r,z) and odd for B_ϕ. While such symmetric fields seem to be realized in <cit.> simulations, the work of <cit.> shows that accretion disks, at least in the non ideal case, can accommodate a plethora of different symmetries. 

A final geometrical aspect is the boundary conditions imposed at the disk equatorial plane. As discussed  in section <ref>, the number of allowed MRI-like wavelengths depends critically on them. In our case, the assumed accretion at z=0 drives a positive electromotive force (p or J_ϕ >0) leading to the generation of a positive radial magnetic component B_r within the disk. So, by assumption, solutions with an outward radial motion at the mid plane (hence B_r <0) have been discarded. Such solutions could accommodate as well with the required outward bending at the disk surface, but probably with an extra half MRI-like wavelength. This can be seen in <cit.> for instance, where u_r(z=0)>0 due to a dominant and positive radial magnetic torque at the disk mid plane. We do seem to recover this behavior. Indeed, as α_m increases (hence α_m → 10 or so, as measured in MRI simulations), both the sonic Mach number m_s and B_r (i.e. the toroidal current p) decrease. Our lack of accretion-ejection solutions for large α_m could thus be an indication that the boundary condition for the accretion speed (sign of p) must be changed. This slight modification can be easily implemented within a self-similar approach.     


As discussed in Sect. 2.2, our work assumes a magnetic diffusion of turbulent origin. Indeed, the existence of self-confined jets in a wide range of objects advocates for a universal mechanism that would be independent of the physical conditions within disks, and in particular of their ionization degree. The natural source for magnetic diffusion in low ionized plasmas is ambipolar diffusion ν_AD. Would it be dominant, the induction equation of the toroidal magnetic field would be deeply affected, modifying significantly the generation of the toroidal field and possibly our results. Ambipolar diffusion can be estimated in the disk mid-plane as ν_AD≃ V_A^2/ν_ni, where ν_ni= 1/τ_ni is the neutral-ion collision frequency. Comparing this expression with our turbulent prescription ν_m= α_m V_A h shows that ν_AD is negligible whenever α_m ≫μ^1/2 A_m^-1, where A_m = V_A^2/ν_ADΩ_K= 1/Ω_K τ_ni is the Elsasser number. Thus, for Elsasser numbers around unity <cit.>, this estimate shows that ambipolar diffusion can be safely neglected with respect to turbulent diffusion (as long as turbulence is going on). This would be in agreement with the universality of accretion-ejection. However, because of the lack of ionisation in certain regions, circumstellar accretion disks are known to harbor dead zones with accretion occurring only at the disk surface <cit.>. We believe that the interdependent accretion-ejection structure will be mostly the same (ie same link between parameters), with an offset from the disk equatorial plane. Note that such a layered accretion structure could be actually described within a self-similar approach, by designing proper vertical turbulent profiles and boundary conditions at the disk equatorial plane. This is left for future work.  

The second major caveat of our work is related to the prescriptions used for the MHD turbulence, namely the viscosity and magnetic diffusivities. As argued before, our choice of α_v and α_m are consistent with our current knowledge of MRI and the way the stress (viscosity) scales with the initial magnetization μ (<cit.> and references therein). We would like to stress however that knowledge on the turbulent diffusion of magnetic fields is scarse. Global simulations <cit.>) and shearing box studies <cit.> report an effective magnetic Prandtl number 𝒫_m of order unity, but this is far from being fully assessed. Besides, the anisotropy χ_m of MRI turbulence has been measured only in non stratified shearing box setups, i.e. a very idealized configuration.  

Even though our prescriptions and scalings agree qualitatively with MRI turbulence, the vertical profiles of the turbulent coefficients used in our work (a Gaussian exp(-x^2)), do not seem to be in agreement with recent numerical studies of accretion disk turbulence. Not only the scale height for turbulence is much larger than the disk scale height <cit.>, but also the profile of the turbulent "viscous" stress does not seem to simply scale as ν_v ρ (see for instance Fig. 5 in <cit.>). In addition, several potentially important processes related to turbulence have been largely ignored in our model and in particular the pressure due to turbulent magnetic fluctuations. This term is known to strongly affect the disc vertical equilibrium <cit.> for μ≳ 10^-3, which could increase dramatically the disc thickness and therefore the quantitative predictions of our model. This advocates therefore for the use of more elaborate closure prescriptions, possibly educated from 3D simulations of MRI turbulence. Note that enhanced diffusion at higher altitudes (due for instance to parasitic instabilities such as Kelvin-Helmholtz <cit.>) could smooth out the spatial oscillations, building up a magnetic configuration closer to that obtained in 3D simulations. We leave these modifications of the self-similar prescriptions for the future.  

	




§ CONCLUSION



Motivated by recent global 3D simulations of accretion disks threaded by a weak vertical magnetic field and showing the launching of jets, we revisited the self-similar accretion-ejection solutions for cold (isothermal) magnetic surfaces. By allowing spatial oscillations of all quantities within the disk, we have been able to extend the previous parameter space by 4 orders in magnitude in the disk magnetization μ, namely from μ= 10^-4 to almost unity.   

We recovered the previous solutions and found a new class of MRI-like driven outflows from weakly magnetized disks, in agreements with some simulations. The role of MRI-like spatial oscillations is shown to be essential in order to provide the required bending of the poloidal field lines at the disk surface. Cold outflows from weakly magnetized accretion disks have the tendency to be more massive than their strong field (near equipartition) counterpart, leading to a critical Alfvén surface closer to the disk surface.  

There is a continuity in behavior as μ increases. Low μ isothermal solutions are quite massive with a typical ejection index ξ∼ 0.1 (increasing with μ) and are mostly driven by the pressure of the toroidal field. The previously published high μ solutions are much less massive, with a typical ejection index ξ∼ 0.01 (decreasing with μ), thus faster and mostly centrifugally driven. These are however two manifestations of the same magnetic acceleration process, linking accretion to ejection in an interdependent way.      

It is striking to realize that the confusion between the Shakura-Sunyaev viscosity parameter α_v <1 and the turbulence parameter α_m has led to restrict all past self-similar papers to values α_m ≤ 1. However, modern 3D simulations seem to imply α_m larger than unity instead. We explored, for the first time, the accretion-ejection behavior under such circumstances, allowing to reach the turning point situation where more angular momentum is being transported radially within the disk than vertically into the jets.    

We propose a simple criterion, based on the initial jet magnetization σ, allowing to discriminate between winds and jets. Solutions with σ>1 are Poynting flux-dominated jets, representative of fast, tenuous and self-collimated outflows, whereas solutions with σ<1 are matter-dominated winds, namely massive, slow and weakly collimated outflows. While the first kind (jets) is achievable at all magnetization levels (and are shown here), only weakly magnetized disks (μ < 10^-2 or less) could provide winds. However, some energy input must be added in order to provide a positive Bernoulli integral. A further development would thus be to include heating at the disk upper layers, 
as in <cit.>, mimicking the existence of irradiation from a central source. This is known to dramatically enhance the mass loss ξ as well, further decreasing σ and allowing for magneto-thermal winds. 

The disk magnetization μ appears to be the main control parameter for determining the intrinsic accretion-ejection properties, such as MHD turbulence and the fraction of the disk angular momentum that is transported by the jets. The existence of some external illumination would then be an extra factor allowing to change jet/wind properties through mass loss enhancement (ξ). Combining these two properties allows to draw an interesting framework. As argued in <cit.>, it is reasonable to expect μ to be a decreasing function of the radius, the magnetic field being dragged in by the  accretion flow. If the innermost regions reach near equipartition (μ between 0.1 and 0.8), then a proper JED solution can be established, with supersonic accretion and fast self-confined jets (low ξ). Such inner regions would have clear astrophysical signatures, in young stellar objects <cit.> and around compact objects <cit.>. The outer disk regions could have a much lower magnetization μ <<1 and accrete at a subsonic pace, while launching massive winds whenever an efficient irradiation is present. Whether or not such Wind Emitting Disks or WEDs are generic in astrophysics requires further investigation.  

 



§ ACKNOWLEDGEMENTS

We thank the referee for providing thoughtful comments on the manuscript. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 815559 (MHDiscs))







mnras
 














§ SELF-SIMILAR EQUATIONS



For the sake of completeness, the full set of MHD equations solved are reported in this section. We define the self-similar functions f_i(x) with x= z/h and h=ϵ r
	
    ρ   =ρ_orr_o^ζ_4f_4 ,        
    		P=P_orr_o^ζ_10f_10 ,
    
    		T   =T_orr_o^ζ_7f_7 ,        
    		u_z=ϵ u_orr_o^ζ_3f_3 ,
    
    		u_r   =-u_orr_o^ζ_2f_2 ,        Ω=Ω_orr_o^ζ_5f_5 ,
    
    		B_ϕ   =qB_orr_o^ζ_1-1f_1 ,    
    		a(r,z)=a_o rr_o^βψ

where the subscript "o" stand for a quantity evaluated at the disk equatorial plane (x=0). Here, q= μ_o J_r0h/B_0 is the normalized radial current density, a_o= B_o r_o^2/β the magnetic flux with B_o the vertical field component, Ω_o= δ_o Ω_Ko the angular velocity, P_o= ρ_o Ω^2_Ko h^2 and u_o= m_s C_s, with C_s = Ω_Ko h defining thereby the accretion Mach number m_s.  

The shape of a magnetic surface anchored at r_o is defined by a(r,z)=a_o and is provided by r= r_oψ^-1/β.  The three transport coefficients ν_v, ν_m, ν'_m (see Sect. <ref>) are written ν_A= ν_Aorr_o^ζ_8f_8 where the profile is a simple Gaussian[Note that <cit.> used f_8(x)=exp(-2x^2), while <cit.> used  f_8(x)=exp(-0.5x^2), namely ideal MHD starting sooner. This may explain why the latter found more massive jets than the former.] f_8(x)=exp(-x^2). Inserting these self-similar functions into the set of PDE (<ref>-<ref>) allows to separate them into an algebraic set of equations on the exponents ζ_i and a set of ODEs on the functions f_i. This leads to the unique solution for a near-Keplerian, gas supported, accretion disk

    β   = 3/4+ξ/2,        ζ_1=ξ/2-1/4,         ζ_2=ζ_3=-1/2
    ζ_4    = ξ-3/2,         ζ_5=-3/2,           ζ_7=-1
    ζ_8    = 1/2,              ζ_10= ξ-5/2

where ξ is the exponent of the disk accretion rate Ṁ_a ∝ r^ξ. 	

Defining  f̃_4 = ln f_4 and f_i'= df_i/dx, allows to express mass conservation and the equation of state as the following ODEs

    f̃'_4 (f_3 + x f_2)     =  (ξ - 1)f_2 - f'_3 - x f'_2  
    
    f_10   =   f_4f_7

Similarly, the radial, vertical and toroidal momentum transport equations become respectively

    m_s^2 ϵ^2 f_4(-ζ_2 f_2^2 + f'_2 (f_3 + xf_2))    =  -f_4 δ_0^2 f_5^2 + f_4/(1+x^2ϵ^2)^3/2
          + ϵ^2 ( ζ_10 f_10 - x f'_10) 
          + μ q^2 ϵ^2f_1 (ζ_1 f_1 - xf'_1) 
          + μΔ'ψ/β(ψ - xψ'/β) 
    
    m_s^2 ϵ^2 f_4(-ζ_3 f_2 f_3 + f'_3(f_3 + x f_2))     =  -x f_4/(1+ x^2ϵ^2)^3/2 - f'_10
         - μ q^2 f_1 f'_1 - μψ'Δ'ψ/β^2 ϵ^2
    
    2f_4(f_3 + x f_2)f'_5 - f_2 f_4 f_5     =   Λ/1+Λ(ψ f'_1 - ζ_1/β f_1 ψ' )
         - 1/1+Λfturb

where the modified laplacian (toroidal current density) is

    Δ'ψ =  ψ”(1+ϵ^2 x^2) + ϵ^2[(2β -3)x ψ' - β(2-β)ψ]

and Λ=p/𝒫_mϵ-1 is the ratio of the magnetic to the viscous torque at the disk mid plane. The function 
fturb= f_4f_8 is the prescription used for the turbulent stress. While the above ODEs are valid both in the disk and in the ideal MHD jet regime, the induction equation requires to deal with each regime in a separate way.  

Within the resistive disk, Ohm's law (<ref>) and the induction equation (<ref>) become respectively

    f_8 Δ'ψ   =   -ℛ_m ϵ^2 ( β f_2 ψ - ψ'(f_3 + xf_2)) 
    
    (f_8 f'_1)'     = ϵ^2 x(f_8(ζ_1 f_1 - x f'_1))' -ϵ^2 f_8(ζ_1 f_1 - xf'_1)(β - 5/2) 
          - χ_m ℛ_mδ_o/qm_s(3/2βψ' f_5 + ψ f'_5)
      + χ_m ℛ_m ϵ^2 β f_1 f_2  
          + χ_m ℛ_m ϵ^2(f_3 + x f_2)(f'_1 - f_1 f̃'_4)

where ℛ_m=r u_o/ν_mo= p/ϵ is the magnetic Reynolds number. When the ideal MHD regime becomes relevant, these equations write respectively

    (f_3 + xf_2) ψ' = βψ f_2 
       (f_3 + xf_2)(f'_1 - f_1 f̃'_4) = δ_o/qm_sϵ^2(3/2βψ' f_5 + ψ f'_5) - β f_1 f_2

We need to complement this set of ODEs with an energy equation providing f_7. Isothermal magnetic surfaces are represented by T/T_o= 1 along each surface anchored at a radius r_o, which translates into f_7= ψ^-1/β. The system of ODEs requires the following boundary values 

    f _  1  ( 0 )     = f _  3  ( 0 ) = 0 
     
    f _  2  ( 0 )     = f _  4  ( 0 ) = f _  5  ( 0 ) = ψ ( 0 ) = f _  7  ( 0 ) = f _  10  ( 0 ) = 1
     
    f _  2  ^ ' ( 0 )     = f _  4  ^ ' ( 0 ) = f _  5  ^ ' ( 0 ) = ψ ^ ' ( 0 ) = f _  7  ^ ' ( 0 ) = f _  10  ^ ' ( 0 ) = 0
     
    f _  1  ^ ' ( 0 )     = - 1 
     
    f _  3  ^ ' ( 0 )     = ξ - 1


We thus get a complete set of equations that can be formally written as M.X=P, where M is a matrix and P a vector depending only on the variable x and the functions f_i, while X is a vector of their derivatives f'_i. Propagating the equations requires to get X= M^-1P, where M^-1 can only be computed as long as the determinant of the matrix M does not vanish. This occurs at the disk equatorial plane x=0 (which is a fixed point of nodal type) and at each critical point of the outflow (see <cit.> for more details). The integration cannot therefore start at x=0 and a Taylor expansion must be made.






§ THE GRAD-SHAFRANOV CONSTRAINT




The Grad-Shafranov equation (<ref>) or GSE requires a regularity condition at the Alfvén point, where the Alfvénic Mach number m= u_p/V_Ap is equal to unity. Its position (r_A, z_A) is labelled by the angle Ψ_A such Ψ_A= z_A/r_A. This regularity condition provides the value g_A= g (Ψ_A,κ,λ, ω,e) of the amount of the poloidal electric current that remains to be used beyond this point, normalized to the current available at the jet base <cit.>. It is therefore known for a given position Ψ_A of the Alfvén point and the MHD invariants κ,λ, ω,e. The angle Ψ_A is however not known a priori. This is due to the fact that the GSE (<ref>) is a PDE of mixed type, where the shape and position of the Alfvén surface need to be imposed. In a time-dependent problem, it would naturally emerge by taking into account all causal connections. Within a self-similar approach of jets, this translates into a conical Alfvén surface with a possible choice of the angle Ψ_A (as done for instance in <cit.>). In our case however, we do have an extra geometrical constraint since the ideal MHD jet is connected to the accretion disk. Our integration starts from the disk equatorial plane and is propagated upwards (increasing variable x=z/ϵ r) so that a trans-A solution can be found only if, at some point  x_A, g(x_A) =g_A is verified, which then fixes also Ψ_A. So the position Ψ_A of the Alfvén surface emerges also as function of the disk parameters. 

Once g_A is known, all jet quantities can be computed at the Alfvén surface as function of their mid plane value and jet invariants: Ω_A= ωΩ_Ko (1-g_A), ρ_A= μϵ^2 κρ_o, u_pA=   _AΩ_Ko r_o, B_pA= κ_A B_o, |B_ϕ/B_p|_A= g_A(ωλ)^1/2/_A.  In these expressions, the value _A gives the poloidal velocity at the Alfvén point and is obtained from the Bernoulli invariant (<ref>), namely _A^2= ωλ (g_B^2 - g_A^2) where, for a cold flow (negligible enthalpy), one has
 
    [              g_B^2=1+ 2 e'/ωλ                               e' =- 2+ω^2/2 + sinΨ_A √(ω/λ) ]

Here, g_B represents the maximum value for the acceleration efficiency g_A as imposed by energy conservation. Since g_B^2>0, this gives a constraint on the minimum energy reaching the Alfvén point, namely a minimum magnetic lever arm.  The closer the Alfvén surface and the smaller λ can be. We found super-A solutions close to the limiting value λ=3/2. For λ>3, the position of the Alfvén surface plays no role anymore.   

The velocity at the Alfvén point can be seen as a lower limit of the jet terminal velocity u_p∞, while the upper limit is  √(2λ - 3) for cold flows (in units of the Keplerian speed at the footpoint). These two limits are shown in Fig. <ref> for our super-A solutions obtained in the fiducial case. For the lower limit, we used for simplicity ω=1, sinΨ_A=1 and g_A=0. Although rather crude, these two limits do bracket quite well our numerical solutions.  




Writing the GSE at the Alfvén point leads to 

    ℰ -  .  a · m^ 2 /μ _  o  r_A^2ρ_A  |_A + g_AΩ_*Ω _  * r_A^2 + (1-g_A)Ω_*r_A^2Ω_*
    +  .  B _ ϕ ^  2  +B _  p ^  2 /μ _  o ρ_A |_A lnη  =0

where  ℰ(a)=E(a)-Ω_*^2r_A^2. Computing the derivatives of the MHD invariants is quite simple within the self-similar ansatz, leading to  

    A=ζ_A/βA/a_o=ζ_A/βA/B_or_o^2

for an invariant A of radial exponent ζ_A. 

Self-similarity introduces a geometrical constraint by imposing that the Alfvén surface is a cone. Along a magnetic surface, one has necessarily B_z/B_o - z/r B_r/B_o= (r/r_o)^-2. Defining the local jet opening angle as tanθ= B_r/B_z, allows to write 

    cosθ -  z/rsinθ = B_o r_o^2 /B_p r^2

which is verified everywhere along a magnetic surface, and in particular at the Alfvén point.  

Making use of this and remembering that m^2 is only a function of the self-similar variable x, leads after some algebra to 

    .  a · m^ 2  |_A = 2/g_A(r_A/r_0)^2B_pA^2/B_0( cosθ_A ω/κλ_A-1)

Inserting this expression into the GS constraint provides

    g_A(gGS-g_A)=(g_B^2-g_A^2) (1-cosθ_A ω/κλ_A )

where 

    gGS=3/4 - 2+ω^2/4ωλ - ζ_4/4 g_B^2

is another maximal value for g_A, imposed by the jet transverse equilibrium. The constant ζ_4= ξ -3/2 is the radial exponent of the density (it comes from the lnη term). Noting that the jet opening angle writes  

    cosθ_A = ω/κλ_A+z_A/r_A√(1+(z_A/r_A)^2-(ω/κλ_A)^2)/1+(z_A/r_A)^2

and inserting it into Eq. (<ref>) allows finally to express the GS constraint into a quadratic equation on X= g_A/g_B

    (k^2c^2 + cos^2Ψ_A)X^2 - 2cX(k^2 - sin^2 Ψ_A) + k^2-1/k^2(k^2- sin^2 Ψ_A) = 0

where c= gGS/g_B and k^2= κ^2/κmin^2 with the minimum mass load κmin defined with

    κmin^2λ^3 g_B^2=ω
 
Which corresponds to the absolute lower limit for getting a super-A flow. Indeed, for g_A=0 the GS constraint can only be satisfied for a minimum value k^2=1 (see Eq. <ref>). Equation (<ref>) is none other than a generalization of equation (3.1) in <cit.>. 

The GS constraint (<ref>) shows that there are always two positive roots 

    . g_A/g_B |_±=   c(k^2 - sin^2 Ψ_A) ±cosΨ_A √((k^2 - sin^2 Ψ_A)(c^2- k^2-1/k^2))/k^2c^2 + cos^2 Ψ_A

In the limit k^2=1 (super-A solutions found along the solid curve in Fig. <ref>), the GS constraint provides g_A=0 and 

    g_A = g_GS2 cos^2Ψ_A/c^2 + cos^2 Ψ_A


When λ increases (ξ decreases), both g_B and the contrast factor c tend to unity so that solutions with large g_A become possible. When λ decreases (ξ increases) but remains bigger than 3, g_B goes to zero very slowly. Hence,the Alfvén surface can still get closer to the pole (Ψ_A decreases), to allow for a more efficient acceleration (larger g_A). But when λ≲3 (ξ≳0.25) and decreases, g_B goes to zero unless the Alfvén surface moves closer to the disk (Ψ_A increases, see Eq. <ref>). This non monotonous behavior, seen only for super-A solutions at low magnetization levels, is illustrated in Figure <ref>. It can also be seen that Ψ_A increases as the disk magnetization μ decreases. For a given ξ (hence λ), decreasing μ leads to an increase in k and g_A is then maximized by increasing Ψ_A. An intuitive picture would be that the smaller μ, the larger the MRI-like induced magnetic bending at the disk surface (see fig. <ref>). A larger inclination enhances jet acceleration allowing to meet the Alfvén point at a smaller altitude z_A.    
  
    
  







§ THE TOROIDAL FIELD AT THE DISK SURFACE




One of the most important quantities in JEDs is the importance of the toroidal magnetic field component at the disk surface (taken here as the SM point). Using the self-similar expression, the magnetic shear writes

    . B_ϕ/B_z |_SM = q f_1 (x_SM) = μ_o J_ro h/B_o f_1 (x_SM)

showing that it depends on both (1) the amount of the radial current density J_ro flowing within the disk and (2) the vertical profile f_1 (x_SM).  

The conducting disk behaves like an unipolar inductor (Faraday disk), where rotation through a magnetic vertical field gives rise to an electromotive force driving a radial electric current J_ro. The disk drives therefore two electric circuits, corresponding each to one of the jets. Therefore, the value of J_ro (as measured by the shear parameter q) is related to the global electric circuit designed by the existence of these two jets. On the other hand, the disk angular momentum equation (<ref>) requires that q= α_m δ_o (p- P_m ϵ)/2μ^1/2, which shows already the "natural" tendency to have | B_ϕ/B_z |_SM∝μ^-1/2.  

The second important element is the vertical profile of the radial current density J_r, which determines the value f_1 (x_SM).  The induction equation (<ref>) writes

    ν'_mJ_r(z) = ν'_moJ_r0  + μ_or∫_0^z z ·Ω

where J_z and B_ϕ advection have been neglected, in agreement with the thin disk approximation. In the absence of any shear (rigid rotation), one would have ν'_mJ_r(z) = ν'_moJ_r0 and the vertical profile of J_r would only depend on the profile of the turbulent diffusivity ν'_m. In thin accretion disks, the radial shear is dominant and a Taylor expansion of the rhs of the above equation leads to

    ν'_mJ_r =  ν'_moJ_r0 ( 1 - Γ x^2  )     Γ=3/2χ_m/α_m^2p/p-𝒫_mϵ

The differential rotation is therefore counteracting the J_ro electric current, which is crucial to deviate the current towards the disk surface and allow current closure in the jets. Now, the amplitude of this effect, as measured by Γ, is challenged  by the vertical decrease of the turbulent diffusivity ν'_m. Indeed, if Γ>1, J_r is going to quickly tend to zero, despite the decrease of ν'_m. On the contrary, if Γ<1,  J_r is going first to reach a plateau (or may even increase) before decreasing to 0 (at a higher altitude)[Note that we used a profile exp(-x^2) for ν'_m. Using exp(-τ x^2) would have lead instead to a comparison between Γ and τ, Γ> τ giving rise to a small value of f_1 (x_SM) (see Appendix B in <cit.>).]. In the former case, f_1 (x_SM) is small (remember that J_r = - ∂ B_ϕ/∂ z), whereas in the latter case f_1 (x_SM) can be quite large. This is illustrated in Fig. <ref>, where different values of Γ have been obtained by playing with different χ_m. Note that in our approach, Γ is not free but depends on the turbulence parameters α_m,  P_m and χ_m. 

As argued in <cit.>, the magnetic (jet) torque F_ϕ= J_z B_r - J_r B_z ≃ - J_r B_z must change sign around the disk surface so that magnetic acceleration can take place. This requires therefore Γ to be of order unity, allowing J_r to conveniently decrease to zero neither too close to the equatorial plane, nor too far away.  Assuming Γ∼ 1 then leads to χ_m ∼α_m^2 providing ν'_m ∼α_m^-1 V_A h. But this is some optimal estimate as solutions can be found with Γ≠ 1.  

For a given value of the ratio χ_m/α_m^2, the toroidal current parameter p can be adapted in order to get super-SM solutions. When χ_m decreases (and/or α_m increases), solutions with Γ of order unity require p → pmin, with pmin= P_m ϵ. Rather small values of χ_m are thus allowed since the ratio p/(p-pmin) can be very large. On the contrary, when χ_m increases (and/or α_m decreases), p needs to increase but the ratio  p/(p-pmin) is bounded by 1. This explains why the parameter space is disappearing so abruptly in this case.  







§ EFFECT OF THE TURBULENCE ANISOTROPY Χ_M




Figure <ref> shows the effect of χ_m on the existence of super-A flows obtained with ϵ=0.1, α_m=1,  P_m=1. The parameter space shrinks when χ_m decreases below unity, with again a displacement to smaller values of μ, but with a more complex behavior on ξmin and ξmax. For χ_m=2, we find the same behavior as for α_m=0.8, namely only two sets of separated solutions, one at n=0 and the other at n=3. Above this value, we found no super-A solution.

Increasing χ_m leads to the same kind of behavior as reducing α_m. Hence, the parameter space of χ_m=2 and α_m=0.8 are similar (see fig. <ref>). Using the same argument, increasing χ_m forbids the outflow acceleration and  the majority of solutions are wiped out. 

Decreasing χ_m has however a much less pronounced effect as increasing α_m. This is due to the fact that the latter is controlling all magnetic field components while the former only affects the toroidal field. The parameter space for χ_m=0.1 appears quite similar to the fiducial case. However, it can be seen that decreasing χ_m to 0.01 leads to a much reduced parameter space, shifted to lower μ and smaller ξmax. A magnetic shear |B_ϕ/B_z| too large produces a strong vertical pinch on the disk so that solutions tend to have both smaller ξ and μ (large ξ become forbidden for a given μ). However, solutions with a magnetization μ > 10^-2 become now impossible. Indeed, the only possibility to lower the magnetic compression at high μ would be to reduce also the magnetic bending, namely B_r/B_z. 



In practice, our super-A solutions require a toroidal current parameter p → pmin= P_m ϵ (see Appendix <ref>). This points to a situation where the field lines would have a different curvature at the disk equatorial plane, namely J_ϕ <0 (thus u_r >0) at z=0. This situation is forbidden by our choice of boundary conditions and hints to a clear bias of our solutions for χ_m <<1. Whether or not such anisotropy is physically relevant remains however to be assessed.      




