Magnetic outflows from turbulent accretion disks: I. Vertical structure&secular evolution

Astrophysical disks are likely embedded in an ambient vertical magnetic field. This ambient field is known to drive magneto-rotational turbulence in the disk bulk but is also responsible for the launching of magnetized outflows at the origin of astrophysical jets. The vertical structure and long-term (secular) evolution of such a system lack quantitative predictions. It is nevertheless this secular evolution that is proposed to explain time variability in many accreting systems such as X-ray binaries. We compute and analyze global 3D ideal-MHD simulations of an accretion disk threaded by a large-scale magnetic field. We evaluate the role of the turbulent terms in the equilibrium of the system. We then compute the transport of mass, angular momentum, and magnetic fields in the disk to characterize its secular evolution. We perform a parameter survey to characterize the influence of disk properties on secular transport. We find that weakly magnetized disks drive jets that carry away a small fraction of the disk angular momentum. The mass-weighted accretion speed remains subsonic although, there is always an upper turbulent atmospheric region where transonic accretion takes place. We show that a strongly magnetized version of the magneto-rotational instability drives this turbulence. The disk structure is drastically different from the conventional hydrostatic picture. The magnetic field is always dragged inwards in the disk, at a velocity that increases with the disk magnetization. Beyond a threshold on the latter, the disk undergoes a profound radial readjustment. It leads to the formation of an inner accretion-ejection region with a supersonic mass-weighted accretion speed and where the magnetic field distribution becomes steady, near equipartition with the thermal pressure. This inner structure shares many properties with the Jet Emitting Disk model described by Ferreira (1997).


Introduction
The emission of jet-like outflows from accretion disks is ubiquitous in the Universe, as both jets and disks are observed around a variety of astrophysical objects and on a wide range of spatial scales. Jet-like outflows are observed being emitted from active galactic nuclei (AGN) and quasar (Merloni et al. 2003 and references therein;Blandford et al. 2019 ), X-ray binaries (Mirabel & Rodríguez 1999;Corbel et al. 2000;Gallo et al. 2003Gallo et al. , 2005, possibly also nova-like cataclysmic variables (Coppejans et al. 2015) and protoplanetary disks (Burrows et al. 1996;Hirth et al. 1997;Ray et al. 1996;Dougados et al. 2000;Bally et al. 2007). All of these jet-like outflows exhibit several common properties : they are highly collimated, supersonic and highly correlated with the measured accretion rate in their emitting disk (Cabrit et al. 1990;Hartigan et al. 1995;Serjeant et al. 1998;Markoff et al. 2003;. This shows that accretion and ejection are intrinsically linked. Recently, an additional component consisting of a more massive and slower outflow, called winds, have been observed around the same objects : AGN (Ci-cone et al. 2014), X-ray binaries (Tetarenko et al. 2018;Jiménez-Ibarra et al. 2019), proto-planetary disk (Tabone et al. 2017). Whether winds and jets share a common source is still an open debate.
The question of accretion in astrophysical disks is tightly linked to that of angular momentum transport: mass accretion is allowed only if angular momentum is removed from the accreted material. In this framework, it is customary to separate two mechanisms through which accretion can occur: turbulent torques, which transport angular momentum radially in the disk bulk (Shakura & Sunyaev 1973) and laminar torques, usually associated to outflows which carry angular momentum vertically away from the disk (Blandford & Payne 1982). In the simple hypothesis of a Keplerian disk embedded in a large scale vertical magnetic field, both mechanisms naturally emerge. First, the seminal work of Balbus & Hawley (1991) showed that the magneto-rotationnal instability (MRI) is triggered, leading to vigorous MHD turbulence and radial transport of angular momentum. Second, large scale magnetic fields spontaneously lead Article number, page 1 of 27 arXiv:2011.14782v1 [astro-ph.HE] 30 Nov 2020 A&A proofs: manuscript no. JoWinds to the launching of an outflow (Wardle & Königl 1993;Ferreira & Pelletier 1995;Konigl & Pudritz 2000;Ferreira 2002;Pudritz et al. 2007) which carries away angular momentum very efficiently. Historically, these two processes have been studied separately, but it is now becoming clear that these mechanisms are interdependent, and can even be considered as two faces of the same basic physical process .
The launching of outflows has been studied using self-similar models (Ferreira & Pelletier 1995;Li 1995;Ferreira 1997) or simplified 2.5D simulations (Murphy et al. 2010;Stepanovs & Fendt 2016). These effective models generically rely on azimuthally-averaged equations of motions, in which turbulence is modelled by an effective viscosity and resistivity. While they have been successful in elucidating the constraints needed for steady ejection (Ferreira & Pelletier 1995) or the relation between disk structure and jet properties (Ferreira 1997;Stepanovs & Fendt 2016;Jacquemin-Ide et al. 2019); their use of ad hoc turbulent profile limits their predictive power as they are dependent on the prescription used to model turbulence (see Jacquemin-Ide et al. 2019 for a discussion).
MRI-driven turbulence has been studied using both shearing box simulations (Hawley et al. 1995) and more recently global simulations (Sorathia et al. 2010). Most of these simulations focused on the configuration where there is no mean field threading the disk (i.e. no external source of magnetic field), thereby excluding the possibility of a magnetically-driven outflow. In these "zero net flux simulations", it is found that the MRI drives 3D MHD turbulence, which results in an efficient radial angular momentum transport and therefore accretion. It is only recently that, thanks to robust numerical methods and more abundant numerical resources, stratified shearing box models with a mean field have been developed. These models show both turbulent radial torques and magnetically-driven outflows (Suzuki & Inutsuka 2009;Fromang et al. 2013;Bai & Stone 2013). These models however suffer from several drawbacks due to the local nature of the shearing box model: the mass loss rate depends on the vertical domain size and the flow top/down symmetry is unphysical (Fromang et al. 2013). This implies that global curvature effects, absent from the shearing box model, are mandatory to obtain physically valid outflow solutions.
To this end, global numerical simulations combining MRIdriven turbulence and magnetically-driven outflows have been computed (Suzuki & Inutsuka 2014;Zhu & Stone 2018;Mishra et al. 2020). These simulations are known to exhibit an exotic global configuration. For instance, the stationary state found by Zhu & Stone (2018) possess 3 distinct features: a turbulent disk, a powerfully accreting atmosphere, mostly driven by a laminar Maxwell torque, and a tenuous wind that does not transport a considerable amount of angular momentum and mass.
These numerical results are in flagrant contradiction with usual two dimensional "effective" models of wind-emitting disk. In these models, accretion typically occurs in the disk bulk and angular momentum escapes through a magnetized outflow launched from the disk surface (Murphy et al. 2010;Stepanovs & Fendt 2016). Even though some effective models predict accretion preferentially located within the disk surface ( 3h: Guilet & Ogilvie 2013;Jacquemin-Ide et al. 2019), where h is the disk pressure scale-height, none of these models produce accretion in the disk atmosphere at z ∼ 10h, as is observed in the numerical simulations of Zhu & Stone (2018) and Mishra et al. (2020). Furthermore, Zhu & Stone (2018) (see their Fig. 10) and Mishra et al. (2020) (see their Fig. 7) show that this accreting atmosphere features turbulent torques at high magnetization, high above the disk. The dynamical role and the origin of this "second" turbulence it not known as of today.
While most of the literature has focused on mass and angular momentum transport, the question of magnetic field transport is also a key element if one is to predict the secular evolution of a disk subject to a magnetized outflow. It is this secular evolution which then leads to interesting observational evidence, like for example eruptions in X-ray binaries Marcel et al. 2018;Marcel et al. 2019) or dwarf novae (Scepi et al. 2019;Scepi et al. 2020). This question was tackled in the seminal work of Lubow et al. (1994), who showed how the magnetic field could be advected towards the inner regions of the disk by the accretion flow. However, they also show that in the case where accretion is driven by turbulent torques, the field advection is efficiently counter-balanced by field diffusion due to MRI turbulence. Rothstein & Lovelace (2008) later showed that the vertical structure plays a quintessential role on the secular transport of the magnetic field. Indeed, the upper layers of the disk, that are magnetically dominated and have a lower conductivity, could enhance the transport of magnetic field inwards. They showed that in the case of wind driven accretion it is possible to accumulate magnetic flux in the inner regions. Guilet & Ogilvie (2014) later showed that even in a purely viscous case, the vertical structure can help the transport of magnetic field towards the inner regions. However, all of these results rely on simplified models of turbulent transport, which are not entirely self-consistent. It is now mandatory to confirm the possibility of field transport from direct numerical simulations as it would then guide future theoretical and observational developments on this problem.
In this work we aim at understanding the dynamics of a disk subject to MRI-driven turbulence and magnetized winds, and link our understanding to the question of the secular evolution. The paper is organised as follows. Section 2 provides the fundamental equations, and the numerical methods for the simulations and the analysis. Section 3 describes in details our fiducial simulation: its vertical structure, the transport properties and secular evolution of the magnetic field, and the characterisation of the outflow. Section 4 focuses on the secular evolution of the density distribution and the formation of ring-like structures. In section 5 we vary the field strength and disk aspect ratio to study their impact on the processes discussed in section 3 and 4. We finally put in perspective our results and discuss their astrophysical implications in section 6.

Governing equations
The analysis done in this paper is performed using spherical coordinates (r,θ,ϕ). Even though this choice may lead to more cumbersome expressions, the spherical symmetry of the system justifies it. For convenience, we also define the cylindrical radius R = r sin θ and the cylindrical height z = r cos θ. We assume the disk is sufficiently ionized to be described in the ideal MHD approximation. In this limit the equations describing a magnetized fluid around a compact object of mass M can be written as where ρ, u, P, and B are respectively the plasma density, velocity, thermal pressure and magnetic field, Φ G = −GM/r is the gravitational potential, Γ is the heat capacity ratio which is fixed to 1.001 and J = c∇ × B/4π is the current. To close Eq. (1-4), we assume an ideal equation of state. We relax our system to a locally isothermal temperature profile T = T (R) ∝ 1/R using a prescribed cooling where T is the temperature of the plasma, T eff is the relaxation temperature and τ cool is the cooling time scale. To avoid the development of the vertical shear instability (VSI, Nelson et al. 2013) we use a thermal relaxation method on a fixed timescale. The VSI would otherwise appear in a strictly locally isothermal model. Hence we relax the temperature of our system on a time scale τ cool , equal to 0.1 times the local keplerian time scale. This cooling can be rewritten as which leads to an effective equation of state of the form where c s is the local isothermal sound speed defined from the Keplerian angular velocity Ω K = GM/R 3 and the disk geometrical scale height h. Finally, we define the poloidal velocity and magnetic field as u p = u r e r + u θ e θ B p = B r e r + B θ e θ

Numerical method
To solve Eq. (1-5) on a spherical grid, we use the conservative Godunov-type code PLUTO (Mignone et al. 2007). Constrained transport (Evans & Hawley 1988) ensures the solenoidal condition for the magnetic field is satisfied at machine precision. We use a second order Runge-Kutta method for the time step combined with a HLLD type solver with linear reconstruction to handle Riemann problems. The extent of our radial and toroidal domain varies from one simulation to the next, while our latitudinal extent is always [0, π]. We summarise the domain properties, as well as the main control parameters for each simulations in Tab.
(1). In the following, we use the innermost radius r in as the length unit, so that r in = R in = 1 in all simulations. We use periodic boundary conditions in the toroidal direction and zero-gradient in the radial direction (mass flow into the numerical domain is forbidden). As for the latitudinal direction we implement the same polar (z-axis) boundary condition as Zhu & Stone (2018).
We follow Nelson et al. (2013) for the hydrostatic initial condition. The initial density and temperature profiles are defined at the disk mid-plane as: where, for convenience, we use the cylindrical coordinates and we define R min = max (R, R in ) to avoid singularities at the pole. As explained above, the temperature is constant on cylinders and remains so during the entire the simulation thanks to the ad hoc cooling function. We can write the solution of Eq. (1-5) assuming a hydrostatic equilibrium (neglecting all velocities but the azimuthal one) and no magnetic fields using the initial density and temperature at the disk mid-plane defined above. We get where we have defined the Keplerian velocity V K = Ω K R. We then use the Keplerian angular velocity and the sound speed to define the disk scale height Article number, page 3 of 27  Table 1. Parameters, grid spacing and radial extension of the simulations discussed in the paper. Note that the latitudinal grid is not homogeneous, with grid spacing varing linearly (L) close to the disk and geometrically (S=stretched) close to the poles. We define where is the disk geometrical thickness. Following Zhu & Stone (2018), we assume a large scale vertical magnetic field that is constant with respect to z but follows a power law in R. In order to get a self-similar initial condition, we choose a radially constant initial plasma beta, in the disk midplane. The magnetic field radial profile is then related to that of temperature and density: To ensure that ∇ · B = 0 at t = 0, we initialize the magnetic field using the magnetic potential A. We can solve the equation Finally, we choose p = −3/2 and q = −1 for all our simulations, in the same way as Mishra et al. (2020).
To handle runaway magnetizations close to the polar boundary we implement a density floor in our simulation where ρ fl,0 = 10 −7 and a fl ∼ 1. Finally, to minimize integration time, we enforce a maximal Alfvénic velocity V A,max = V K (R in ), where V A = B/ 4πρ is computed with the instantaneous magnetic components. This maximal Alfvénic velocity is imposed by changing the local density and not the magnetic field, hence it acts as an additional density floor. In our simulations, this numerical artifice is only noticeable close to the axis as well as near to the inner boundary, r ∼ r in , which are ignored in our analysis.
We use a nonuniform grid ( Fig. 1): in the radial direction the grid points are logarithmically spaced, while in the latitudinal direction the grid is linearly spaced up to |z| = 5.5h and then geometrically stretched up to the axis. The grid is identical for all the simulation except for the simulation at = 0.05 where the grid is spread differently. For this simulation the grid is linearly spaced up to |z| = 3h and then stretched up until |π/2 − θ| > π/2 − 0.35 where the spacing becomes linear again for two grid points close to the polar boundary. This particular grid maximises the resolution within the disk and its surface as well as providing a buffer zone close to the polar boundary where the CFL condition is the smallest. Finally, in the toroidal direction the grid is linearly spaced. The resolutions for the different simulations per region and direction are displayed in Tab. (1). This choice of resolution allow us to resolve the MRI within the disk thanks to roughly 16 grid points per vertical scale height.
We integrate our simulations up to T end , see Tab.(1). We define our reference time scale T in = T K (R in ), corresponding to orbits of the innermost radii, where T K (R) = 2π/Ω K (R).

Averaging procedure
To quantify the turbulent state of the accreting system, we need to define a set of averages. This allows us to decompose our quantities in turbulent, δX, and laminar, X , terms. We use three kinds of average, on ϕ, on (ϕ, t), and on (ϕ, t) density weighted: In these definitions, ϕ 1 and ϕ 2 are taken as the boundaries of the whole computational domain (Tab. (1)). It is important to note that X X ρ , since we define our turbulent fluctuations from the · average With this definition, the quantities X (R, z) are equivalent to the ones used in 2d effective models such as self-similar solutions (see Jacquemin-Ide et al. (2019)). Hence, they allow for a direct comparison to said models and provide guidance for turbulence closures. Let us emphasize however that δX ρ 0, which leads to additional complications when dealing with high order correlation functions. Finally, note that since the boundary condition in the ϕ direction is periodic. The averages are taken from time t 1 to t 2 , with a temporal resolution of ∆t = 1/Ω K (R in ) 0.16T in and then azimuthally averaged over the whole domain (Tab. 1).

Turbulent decomposition
In order to study the impact of turbulence on the vertical structure of the accretion-ejection system, we use the averages defined above to rewrite the ideal MHD equations as a set of Reynolds-averaged MHD equations. In this approach, turbulence then appears in the averaged equation as an additional dynamical term influencing the mean flow. This is actually an extension of the famous Shakura & Sunyaev (1973) α-disk approach which can be derived as a Reynolds-averaged model for MHD turbulence (Balbus & Papaloizou 1999). While the α disk approach only focuses on angular momentum transport, our objective is to generalise this to the whole set of ideal MHD equations. To illustrate our procedure, we show its application to the equation of angular momentum transport below. A complete derivation for the whole set of ideal MHD equations can be found in Appendix A.
We project and rearrange (using Eq.1) the equation for toroidal momentum transport to find the equation of angular momentum transport We decompose the different quantities into fluctuating (turbulent) plus mean (averaged) terms. We will also call mean terms, laminar terms, in the same way as Béthune et al. (2017) (see also Mishra et al. 2020). We then combine it with the continuity equation (Eq. 1) and average the result with respect to time and ϕ. We eventually get For the sake of brevity, we have neglected the term with the least dynamical importance and assumed stationarity, the equation can be found in full form in Appendix A. We then define a fluctuating stress tensor as well as a laminar stress tensor: where we define and We can then define the usual turbulent alpha viscosity as where α la is related to the laminar torque and α tu to the turbulent one. The torques have different origins as one is a consequence of the turbulence while the latter is a consequence of the mean field geometry. It is important to note that MRI is expected to induce a radial turbulent torque while an MHD wind is expected to induce a latitudinal/vertical laminar torque. However, we see in Eq. (14) that the fluctuating turbulent and mean laminar terms have both a radial and a vertical component. Furthermore, in the present work we find 1 that the radial laminar torque can be important in the angular momentum budget. It is not clear whether such a torque is a consequence of the wind or the turbulence. This calls into question the a priori distinction of wind vs turbulent torques and may require a change of nomenclature within the community.

Global picture
We first present a global picture of our fiducial simulations. We show in figure 2 the mean poloidal stream lines and field lines, the Alfvén surface and the fast magneto-sonic surface of the fiducial β ini = 10 4 run. We define the Alfvén surface as the location where u p = V ap , where V ap = B p / 4π ρ and the fast magneto-sonic surface where In the top panel of Fig. (2) we show the whole domain of our numerical simulation. We find that the upper and lower hemispheres of the system are approximately symmetric. The fast magneto-soni (FM) surface is conical up to the border of our domain. In contrast, the Alfvénic surface remains conical up to a certain radii where its shape is modified, demonstrating that the outer radii of the domain has not yet reached a stationary state. This is corroborated by the lack of toroidal field within the outer regions of the disk.
In Fig. (2) (top right) we can clearly see some collimation happening for the inner field lines, close to the end of their trajectory. This may be the sign of some intrinsic MHD property, as it has been also observed in self-similar jet solutions (Ferreira 1997; Jacquemin-Ide et al. 2019). However, it is unclear whether this is a physical or a numerical effect, akin to a bias of our boundary conditions or a consequence of the density floor. Nevertheless, a study of collimation in our numerical models is outside the scope of this paper.
Since the FM surface corresponds to the point where the outflow becomes causally disconnected from the disk, we can use this surface to estimate the last radius that has achieved stationarity. We do this by computing the anchoring radius of the last field line that traverses the FM surface, this corresponds to R ∼ 35. Hence, we will focus the rest of our analysis to a domain within that radii, the stationary inner regions (grey rectangle in top of Fig. (2)).
We can then identify three zones shown in the bottom panels of  The region we identify as the accreting atmosphere exhibits turbulence which is not obvious from Fig. 2. To quantify the dynamical importance of this turbulence, we measure the ratio α tu /α la (Fig. 3, top). As expected, the disk is clearly dominated by the turbulent stress component. As we move through the accreting atmosphere from the disk surface, the stress becomes dominated by its laminar component before turning back to a turbulent stress region, as in the disk. Even higher still, in the wind region, the torque is laminar. The fact that the atmosphere is divided into two sub-regions, dominated respectively by a laminar and a turbulent stress, invites us to define a laminar and a turbulent atmosphere (LA and TA), which we need to delimit more precisely.
Under the assumptions of ideal MHD and steady state, the average poloidal velocity and magnetic field should be perfectly aligned (e.g. Mestel 1961). Hence, the value of the angle ψ between the averaged poloidal stream and field lines, identifies the regions where the turbulent terms are important. In Fig. 3 (bottom) we plot this angle as a function of the latitudinal coordinate. We also provide in Fig. 4 the behaviour of the mean velocity and mean magnetic fields as functions of the latitudinal coordinate.
In the disk, the mean poloidal velocities and field lines are approximately perpendicular (cos(ψ) = 0): the disk is turbulent. We identify the disk surface as the location where cos ψ = ±1. In the following, we define the disk as |z| < 3.5h (the central shaded area in Fig. 4 and Fig. 3 (bottom)), consistently with Fig. 3 (top). The disk surface coincides with the location where the average radial velocity becomes negative (Fig. 4), marking the beginning of the accreting atmosphere.
Above the disk, the average poloidal stream and field are aligned (cos(ψ) = ±1). This defines the laminar atmosphere previously identified with the stress, located at altitudes 3.5h < |z| 6.5h. As we move to even higher altitudes, cos(ψ) changes sign twice (Fig. 3). This suggests that ideal-MHD is broken in this region, which we identify as the turbulent atmosphere, located at altitudes 6.5h < |z| 12.5h. This region ends at the surface delimiting the base of the wind, where cos(ψ) = ±1, which coincides with the slow magneto-sonic (SM) surface defined as ). It is important to note that V S M is indeed constant in units of c s since the magnetic field is dominant inside the atmosphere (laminar and turbulent), β = 8π P / B 2 < 1. However, the mean poloidal magnetic field stays infra-thermal within the laminar atmosphere, β p > 1, where but becomes supra-thermal within the turbulent atmosphere, β p < 1.

Mass transport
Using the regions defined in section 3.1, we compute the fluxes through the boundaries of each region, delimited radially by r 1 = 8 and r 2 = 14. We then analyse the transport of mass and   : mean magnetic field profiles normalised to the vertical magnetic field in the disk mid-plane, B z (r, θ = π/2) as functions of the latitudinal coordinate. The grey zones correspond to the turbulent regions in figure 3. All profiles are radially averaged between r 1 = 8 and r 2 = 14 after being normalised. angular momentum within our system. The analysis is carried out in detail in appendices B.1 and B.2. In this section, we summarise the main results as they have already been described in details by Zhu & Stone (2018).
Within the disk the mass is latitudinally transported towards the atmosphere and radially transported outwards, as described by Zhu & Stone (2018). Moreover, the total mass within the disk section delimited by r 1 and r 2 decreases with time. This is the consequence of the secular self-organisation of the radial disk structure. Indeed, a closer look at the radial density profiles reveals the formation of ring-like structures. Our chosen disk section turns out to be located on a gap. The mechanism behind the formation of these structures is discussed in section 4.
In the accreting atmosphere the mass is transported radially inwards, i.e. gas is accreted. In the laminar region the material falls in along the mean poloidal field. Most of the actual accretion happens in the turbulent atmosphere, where material slips though the magnetic field thanks to a turbulent "effective" resistivity. In the latitudinal direction the atmosphere receives mass from the disk and looses a negligible amount to the outflow. In contrast to the disk, the transport of mass within the atmosphere is in steady state, as no mass excess or deficit is observed.

Angular momentum transport and accretion
Within the disk, angular momentum is transported radially outwards through turbulent stresses. Moreover, the disk receives angular momentum from the atmosphere through laminar latitudinal stresses, which explains the observed decretion in the disk (see also Zhu & Stone 2018). The accreting atmosphere losses angular momentum and transfer it mostly to the disk through a laminar stress. As previously shown by Zhu & Stone (2018), the wind contribution to the angular momentum budget is essentially negligible.
To understand the global behaviour of the accretion flow in our simulation, we compute the mass weighted accretion velocity, where θ SM,1 and θ SM,2 are respectively the angles where the flow becomes super SM in the upper and the lower hemisphere. These surfaces coincide with the end of the turbulent atmosphere. The mass weighted accretion velocity is related to the mass accretion rateṀ acc , indeed v acc is the velocity at which mass is transported. We then time average this quantity between t a = 318T in and valid for R ∈ [2, 14], which is clearly subsonic. The accretion velocity, v acc , follows a radial dependency very close to Keplerian. The negative sign shows that accretion in the turbulent atmosphere dominates over the disk decretion and selforganisation. Moreover, the accretion time scale we deduce from this: t acc ∼ 9.1×10 2 /Ω K (R) = 145T K is substantially longer than than the local dynamical timescale of MRI or ejection which are of the order of T K . Hence, there is a clear timescale separation between local dynamics and accretion.

Magnetic flux transport
The question of magnetic flux transport in disks is a corner stone of accretion theory (Lubow et al. 1994;Rothstein & Lovelace 2008;Guilet & Ogilvie 2014). To characterize this transport, we introduce the magnetic flux threading the disk as Article number, page 8 of 27 J. Jacquemin-Ide et al.: Magnetic outflows from turbulent accretion disks (23) which corresponds to an integration on a surface that goes from the axis (z = 1 and R = 0) down to the inner boundary of the disk (z = 0 and R = R in ) and then to radius R in the disk midplane. In our fiducial simulation, we first confirm that Ψ(R = R max ) is constant during the entire run, indicating that the total flux within the numerical domain is conserved. We show Ψ in Fig. (5) as a function of the radial coordinate. In this figure, each contour line corresponds to the midplane footpoint of a poloidal magnetic field line. We see that the field lines are advected inwards. Eventually, the magnetic flux accumulates in the central hole of the domain, but is not lost.
Given that magnetic flux accumulates towards the disk center, it is desirable to measure an effective velocity of magnetic field transport, v Ψ , from our simulation. We compute this effective advection speed assuming that the magnetic flux can be modelled following a simple advection equation, as proposed by Guilet & Ogilvie (2012) Note that this equation is only a phenomenological equation designed to evaluate the typical advection speed of the flow. Using the above definition, we compute v Ψ and then average it between t a = 318T in and valid for R ∈ [2, 11]. For R > 11 we cannot conclude on magnetic flux advection, since we did not integrate for long enough to get a measurable deviation from the initial condition. As expected from Fig. (5) the advection velocity is negative within the inner regions of the disk. Moreover, the advection velocity follows a radial dependency close to Keplerian. We cross-checked the numerical value of the magnetic field advection velocity by fitting the contours of Fig. (5) with where R Ψ is the position of the contour, R i is its initial position and t i is the time the system needs to stabilise the transient due to the initial condition and start advecting magnetic flux, a delay visible in Fig. (5). Equation (26) is the integrated form of an equation describing a field line advected at a velocity proportional to the local Keplerian velocity,Ṙ Ψ = v Ψ (R). These fits show a reasonable agreement with the data (Fig. 5, dashed line), which suggests that magnetic field lines are advected in a self-similar manner once the local equilibrium is reached (t i ∼ 30T K (R)).
Since v Ψ ∼ v acc , the time scale for the magnetic field transport, namely t Ψ = 10 3 /Ω K (R) = 160T K (R), is comparable to the accretion time scale and substantially longer than the dynamical time scale.

Dynamical equilibrium
Assuming the total poloidal flow acceleration is negligible we get the following radial and latitudinal equilibria from Eq. (A.11) where g = GM/r 2 and P B = B 2 + δB 2 /8π is the magnetic pressure, which contains a laminar, B 2 , and a turbulent, δB 2 , contribution. We neglect the tension due to the latitudinal magnetic field in Eq. (28) due to its small impact in the latitudinal equilibrium. Figure 6 shows the latitudinal profile of the radial and latitudinal equilibria, respectively normalized to g and gr.
We find that the disk region is in radial and latitudinal equilibrium thanks to the balance between the centrifugal force, u ϕ 2 /r or cot θ u ϕ 2 , the gravity and the thermal pressure gradient. Hence, the disk is in hydrostatic balance as is expected from a weakly magnetized disk. As we enter the accreting atmosphere the thermal pressure gradient becomes close to negligible, when compared with the magnetic forces. The radial and latitudinal equilibriums are then enforced by a balance between gravity, the centrifugal force, the magnetic pressure gradient and the hoop stress, B 2 ϕ /(4π ρ r) or cot θ B 2 ϕ /(4π ρ ). The radial equilibrium is maintained up to the base of the outflow (θ π/2 ± π/4). From this point on, MHD acceleration is no longer negligible and Eqs. (27,28) break down.
Within the laminar and turbulent atmosphere, the system is no longer in hydrostatic equilibrium, but in magnetostatic equilibrium. This magnetostatic equilibrium emerges from a combination of the laminar and turbulent magnetic pressure (Fig. 7). This situation is similar to what was proposed by Begelman et al.
. Left: radial equilibrium (eq. 27), right: latitudinal equilibrium (eq. 28) as a function of the latitudinal coordinate θ. The terms are averaged between r 1 = 8 and r 2 = 14 after being normalized to g (left) or gr (right). If the system is in equilibrium the sum of the terms, the black solid line , should be equal to 1 (left) or 0 (right). The shaded areas denotes the turbulent regions.
pressure emerges as a consequence of turbulence. Indeed, within this region, we see the emergence of a powerful toroidal field, 20 times the vertical field at the disk mid-plane ( Fig. 4 bottom), that changes sign within the turbulent atmosphere. The strong toroidal field within the laminar atmosphere is consistent with the fact that the laminar atmosphere is in ideal MHD. In this region, accretion drags the poloidal field lines inward. Since the fluid is plunging into a region with higher angular momentum, the toroidal speed increases, which leads to enhanced shear for the toroidal magnetic field, increasing its magnitude (Zhu & Stone 2018).
The emergent toroidal field is consumed by the dynamical equilibrium. Indeed, it is the gradient of the magnetic pressure that supports the turbulent atmosphere (Mishra et al. 2020) and achieves latitudinal equilibrium within the accreting atmosphere. Hence, as we move further up into the turbulent atmosphere the toroidal magnetic field decreases (Fig. 4). Within the turbulent atmosphere the toroidal field changes sign, which might be surprising at first sight. However, such a change of sign is mandatory for launching an MHD outflow, as first realised by Ferreira & Pelletier (1995). This allows to switch from the underlying braking torque to an accelerating MHD torque that drives the outflow. It corresponds also to a radial electric current which is circulating outwardly at the TA layer. The origin of this current is the electromotive force driven by the rotating disk across the large scale B z field (see also Ferreira 1997 and Jacquemin-Ide et al. 2019 for more details on the electric current vertical profile).
Overall, it is clear from Fig. (7) that the dynamical importance of the turbulent pressure cannot be understated. Indeed, Ferreira & Pelletier (1995) showed that the mass loaded onto the field lines depends crucially on the disk pressure which pushes material up to the "surface", and loads the field lines. It is clear from Fig. (7) that the laminar magnetic pressure is compressing the disk surface, as shown by Ferreira & Pelletier (1995). However, the effect of this compression seems nonexistent in Fig. (6). This is a consequence of the presence of a important turbulent magnetic pressure in the disk. Indeed, the compression of the disk surface due to the laminar magnetic pressure, B 2 , is compensated by the turbulent magnetic pressure, δB 2 within the disk (Fig. 7). This ensures that material can flow relatively freely from the disk to the atmosphere and enhances the mass loading into the laminar atmosphere region. The same process is at play for the transition between the turbulent atmosphere and the outflow. However, in this case it is the combined effect of thermal pressure and the turbulent magnetic pressure (Fig. 7) that allows the loading of the magnetic field lines at the base of the outflow. The reason for the reappearance of the turbulent pressure within the turbulent atmosphere will be described in section 3.5.

Super fast magneto-sonic wind
In the regions where the flow is approximately laminar and steady-state, it is possible to define a set of MHD invariants which consists of conserved physical quantities along the field lines. These invariants characterize MHD outflows and allow direct comparisons to be made with other outflow solutions.
Under these conditions, Eq. (1) reduces to: whereη measures how much mass is being loaded onto the field lines. It can be transformed into a dimensionless mass loading parameter κ ≡ηR S M Ω K (R S M )/B z,S M , where R S M is the anchoring radius of the magnetic field line and the SM stands for quantities evaluated at the slow magneto-sonic surface.
Following the same regime of approximation, we can also write the ϕ component of the induction equation as where Ω can be interpreted as the angular velocity of the magnetic field lines. This can be recast into a dimensionless number Similarly, the angular momentum invariant is deduced from Eq. (A.13) where R A is the Alfvén radius of the anchored field line, which is the radius where the flow following the field line becomes super Alfvénic. We can then define the dimensionless invariant λ = L/(RΩ K )| S M , known as the magnetic lever arm. It is a measure of the angular and poloidal acceleration of the wind by the disk mediated by the magnetic field. Finally, by projecting Eq.
(2) along the poloidal magnetic surface we can define the Bernoulli invariant where w is the thermal energy term which can be written as where l is the coordinate along the field line. It is important to note that, when heating is allowed, w is a function of the full stream line as a result of its integral form, therefore it cannot be computed from a single point within the field line as the other terms.
However, in our case, the heating contribution w is negligble so that the Bernoulli invariant is fully determined at the SM (launching) point. We normalize this invariant with respect to the gravitational energy at the mid-plane defining e = E E K,S M , where E K,S M = Ω 2 K (R S M )R 2 S M /2. We follow three different field lines originating at R 0 = [5, 6, 8] respectively (Fig. 8, top). We then compute the ideal MHD invariants using Eq. (29-32) across the 3 different field lines anchored at the radii R S M [3.5, 4, 7.5].
We show in figure 8 (bottom) the MHD invariants computed along the field lines. We find that the invariants are approximately constant once the outflow leaves the turbulent atmosphere. Note that the field lines closest to the central object are the ones that exhibit the least variability in their invariants. This is to be expected, as not all anchored field lines are created equal: the farther they are from the central object the longer their MHD invariants will take to converge.
Several features of the wind can be deduced from its MHD invariants (see also Tab.3).
-The rotation invariant ω 1, indicates that field lines are rotating close to the Keplerian speed of their anchoring radius, or slightly slower. This is consistent with the values expected in self-similar models (Jacquemin-Ide et al. 2019).
-The angular momentum invariant λ 4-6, so the wind is effectively free (λ > 3/2). This provides R 2 A /R 2 S M = λ/ω ∼ 8 which is similar to the value found in Zhu & Stone (2018). The magnetic lever arm can be used to estimate the terminal velocity of the outflow using u p∞ R S M Ω K (R S M ) √ 2λ − 3 (Blandford & Payne 1982). Using this expression we find u p∞ ∼ 3R S M Ω K (R S M ).
-The mass loading invariant κ ∼ 0.1 implies that the energetic content of the wind at the base of the outflow is magnetically dominated instead of being kinematically dominated, consistent with a jet-like outflow. Since the laminar atmosphere is in ideal MHD, we can also compute the value of the mass loading invariant within it, finding κ ∼ 8. The value of κ is an order of magnitude larger within the laminar atmosphere that within the outflow. This is consistent with a laminar atmosphere that is so heavily loaded with mass that it falls towards the central object. It is that fall that leads to the generation of a large scale toroidal magnetic field within this region.
The Bernoulli invariant allows us to characterise the wind energetics and its drivers. Fig.(9) shows the different contributions. The conservation of the Bernoulli invariant is striking, especially when compared to the other invariants. The positive sign of the Bernoulli invariant indicates that the flow is free from the potential well and continues its propagation up to infinity, consistently with λ > 3/2. When the flow reaches the end of our simulation box, its energy content is dominated by the kinetic component and the magnetic energy has been mostly consumed, so the outflow is close to its asymptotic state. This confirms what was deduced from the magnetic lever arm λ, mainly that u p∞ ∼ 3R S M Ω K (R S M ), since e ∼ 6.
It is also clear from figure 9 that the outflow is cold, i.e the thermal pressure term w is negligible in the outflow energetics. Indeed, at the wind launching point, the outflow is dominated by the magnetic energy, and even the gravitational energy is negligible at this location. Finally, we also compute the invariants on the south pole of our system but originating on the same radii, R 0 = [5, 6, 8], and find very similar values, confirming the symmetric nature of the system. Jacquemin-Ide et al. (2019) propose a criterion to distinguish between jets and winds, which essentially quantifies whether the outflow can self collimate thanks to magnetic forces ("jet") or not ("wind"). They use the definition of the jet magnetization defined as the ratio of the MHD poynting flux to the kinetic plus thermal energy flux at the base of the outflow. If σ > 1 the outflow can be categorized as a jet, in the contrary it is a wind.
Using figure 9, we derive σ ∼ 8 which is compatible with a jet like outflow and is consistent with the observed collimation.

Origin of turbulence in the atmosphere
As described in section 3.1 the accreting atmosphere is divided into two distinct regions: a turbulent and a laminar atmosphere. This configuration, where turbulence is localized in a layer above the disk is pretty surprising and has never been addressed so far. It suggests that within the laminar atmosphere, the source of turbulence must first be quenched and then reinvigorated further up, within the turbulent atmosphere. Because the turbulent atmosphere is strongly magnetized ( β p 1), several MHD instabilities could in principle be invoked to explain this. However, since the main source of free energy for turbulence is sill the shear, the MRI remains a very serious contender.
It is widely believed that the MRI is a weak field instability and that it is quenched once β p 1. This condition is however only geometrical and based on thin disks, requiring MRI modes to be confined within the vertical thickness of the disk (set by the hydrostatic equilibrium). If one considers a system with a size not determined by thermal pressure, this geometrical constraint vanishes and the MRI can in principle exist for low β p . This is the situation studied by Kim & Ostriker (2000), who confirmed that in the limit β p → 0 the MRI still exists. However, in contrast to the weak field regime, the MRI can be suppressed by a too strong toroidal field. The stability condition in this low beta limit reads . Maximum growth rate of the compressible MRI as a function of the poloidal plasma beta and the ratio between the toroidal and poloidal magnetic field, calculated using Eq. (C.2). The red points correspond to the values of the ratio of the magnetic field and the poloidal plasma beta for different regions within our simulation averaged between r 1 = 8 and r 2 = 12 (see Fig.4). We find that laminar regions are characterised by reduced MRI growth rates ( 0.25Ω K ).
. instability. For the sake of completeness we derive Eq. (35) as well as the dispersion relation for the compressible MRI, Eq. (C.2), in Appendix C. We then use this dispersion relation to compute the maximum growth rate of the MRI as a function of the poloidal plasma beta and the ratio between the toroidal and poloidal magnetic field (Fig. 10). To compare this to our numerical results, we use the profiles of the magnetic field components and pressure in our simulations to characterise the value of the MRI growth rate within each region. From Fig. (10) we find that the regions where MRI grows at a sufficient rate ( 0.25Ω K ) correspond to the zones described as turbulent in section 3.1, while the laminar regions show lower growth rates. We thus conclude that MRI could be the main driver of the turbulence within the turbulent atmosphere.
As explained in section 3.3 the emergence of a powerful toroidal field is consistent with the laminar atmosphere being in ideal MHD. Furthermore, the decrease of the toroidal magnetic field within the turbulent atmosphere is also consistent with the magnetostatic equilibrium of the atmosphere. Hence, we can deduce that the re-ignition of turbulence is required for the establishment of a stationary magnetostatic equilibrium.
The re-activation of turbulence leads to anomalous resistivities within the turbulent atmosphere that allows accretion through the poloidal field lines and a steady state transition between accretion and ejection (Ferreira & Pelletier 1995). It is therefore an essential ingredient for launching the wind.

Formation of ring like structures.
We discussed in section 3.2.1 that the disk region has not achieved a rigorous steady state. Indeed, the plasma within the disk is seen to concentrate into ring like structures, a special type of self-organisation. To understand the formation and evolution of the ring like structures, it is useful to vertically and azimuthally average the equation of mass conservation. Equation (1) leads to where Σ = ρ ϕ is the column density of the disk and the vertical and azimuthal average X is defined in Appendix B. We compute the vertical averages within the latitudinal wedge θ 2,1 = π/2 ± arctan(3.5 ), equivalent to the disk region (section 3.1). We show Σ in a space time diagram in the top panel of Fig. (11). We concentrate our analysis on the gaps, minimum of Σ, located at R [6.6, 12.3] as well as the rings, maximum of Σ located at R [9.5, 15]. Indeed, the ring located at R 2, and to a certain extent the one located at R 4, are within close range of the inner boundary which could bias their analysis. Nonetheless, we find that the innermost ring-like structure (R ∼ 4) starts to emerge at around 200 Keplerian orbits of the innermost radii (25 local orbits) while the farther rings (at R ∼ 10) appear at 800 Keplerian orbits of the innermost radii (25 local orbits). Hence, we can conclude that the local time scale of self-organisation t so (R) 25T K (R) = 160/Ω K (R).
Once formed, the ring structures survive for the entire duration of the simulation. They also exhibit a weak radial inwards migration, more evidently seen for the R ∼ 2 ring. We measure the migration speed of the R ∼ 2 ring to be v ring 5 × 10 −4 V K (R) by measuring the radial velocity of its density maximum. This velocity is twice smaller than v acc . It is unclear whether this drift is a bias of the boundary condition or a consequence of the intrinsic physics. It is important to note that while R ∼ 4 may seem to be migrating inwards we were not able to measure any radial drift.
The magnetic field is probably linked to the production of ring-like structures. To understand the role of the magnetic field in the self-organisation mechanism, we compute the vertically average vertical field b z = B Z ϕ /(7h) and show its space time diagram in Fig. 11, bottom. We find that the vertical magnetic field is accumulating within the gaps, and is anti-correlated with the rings. However, the magnetic structure is not as stable with time as the density structure.
To understand if the self-organisation of the disk into ringlike structures is a consequence of radial or latitudinal transport of mass, we need to turn back to Eq. (36). We average it with respect to time and divide it by the time average of the column density, Σ = ρ , and the Keplerian angular velocity to find the dimensionless equation where we defineσ w = sin θ ρu θ θ 2 θ 1 as the wind mass loss rate. To simplify the comparison of each term, we have normalized equation (37) to the Keplerian frequency. Figure (12) shows the wind mass loss rate as well as the derivative of the radial mass flux. We see that the wind mass loss rate,σ w , is essentially negligible, while the radial term, ∂ r r 2 ρu r shows a good correlation with the ring structures: the radial locations where the radial term is positive (divergent mass flow) correspond to minima of the column density, while the radial locations where the radial term is negative (convergent mass flow) correspond precisely to maximal column densities. We have checked that the mean radial mass flow forming these structures is driven by the radial stress, following the vertical average of Eq. (A.12). This implies that the origin of the selforganisation is due to α tu . From Fig. (12) we can probe the time scale, t so , on which this self-organisation takes places by looking at the amplitude of the radial divergence of the radial velocity in units of keplerian frequency. This leads to t so ∼ 200/Ω K (R) = 30T K (R) this corresponds quite well with the values measured in Fig. (11). Comparing the self organization time scale with the ones deduced from section 3.2.1 we find t so t acc t Ψ . Self-organisation into ring-like structures is not impeded by the global transport mass or magnetic field, it happens on shorter time scales. As mentioned above, the self-organisation is due to radial motions which are driven by the radial component of the turbulent stress tensor. It is tempting to interpret this as a viscous instability, as proposed by Lightman & Eardley (1974). Indeed, Lightman & Eardley (1974) showed that, for the system to be unstable to a viscous type instability, the torque needs to be a decreasing function of the surface density, ∂ log T tu, rϕ ∂ log Σ < 0 . In the case of a strongly diffusive field, the limit where the field is steady, this is equivalent to α tu ∝ β mid −1 , where we define and β mid = β p (r, θ = π/2).
However, in MRI turbulence, the measured scaling is closer to α tu ∝ β mid −1/2 (Salvesen et al. 2016). For the sake of com- pleteness we check the scaling for the turbulent transport of angular momentum in our simulation. In Fig.(13) we sample the values of the averaged turbulent parameter and the local mean plasma beta at the disk mid-plane for radii ranging from R = 4 to R = 16. We expect α tu to not only be a function of β mid , indeed in Fig.(13) we have the same value of α tu for different plasma beta. This is probably a consequence of the dependence of α tu on the local shear, d log Ω/ d log R. Indeed, self-organisation induces local deviation from the Keplerian shear, which in turn influences the efficiency of the MRI-driven angular momentum transport (Gressel et al. 2015;Pessah et al. 2007), weaker shear leading to lower values of α tu . It can be seen that the scaling of the turbulent angular momentum transport coefficient, α tu , verifies a scaling closer to −1/2 than −1.
The −1/2 scaling is inconsistent with a simple viscous type instability, hence the self-organisation of our simulation into a ring-like system is not purely viscous, and must probably be somehow connected to magnetic field transport in the disk.
In conclusion, the fact that the self-organisation is driven by radial motions of matter tends to disfavor a wind-driven instability (Riols & Lesur 2019) and tends to favor some kind of viscously-driven instability (Hawley 2001;Bai & Stone 2014). However, as explained above, a simple viscous instability could not explain the disk instability since the scaling of α tu on β mid is not steep enough. Finally, we note that the radial density structure of the outflow is unaffected by the ring-like structures found in the turbulent disk. The radial structures disappear in the upper layers (LA or TA), probably because the rings are washed out by the supersonic accretion in those regions.
Bai & Stone (2014) construct a phenomenological model for ring formation, that does not require a steep scaling for α tu , based on the evolution of the vertical magnetic field driven by an anisotropic anti-diffusivity. However, it is not evident that an anisotropic anti-diffusion emerges from MRI turbulence. Furthermore, recent measures of the mean field diffusivities of MRI turbulence do not agree with a negative anisotropic diffusivity (Gressel et al. 2015). Finally, Johansen et al. (2009) explore the possibility that stochastic perturbations driven by an inverse turbulent cascade of the magnetic components could be driving the self-organisation of the disk into ring-like structures as well as the formation of pressure bumps. However, their modeling concentrates on unstratified shearing box simulations with a zero-net flux. This leads to their ring-like structures being short lived, lasting at most 50 local orbits while our inner ring-like structure are still stable at 75 local orbits and show no sign of weakening. Nonetheless, a stochastic model taking into account the feedback of the vertical field evolution may be needed to explain the ring formation observed here.

Parameter exploration
To have a comprehensive understanding of how the properties described in section 3 and section 4 evolve with different fundamental parameters of the system we have produced 4 additional simulations, see Tab.
(1). In these simulations we explore the effects of the initial magnetic field strength (SB2, SB3, SB4), the geometrical thickness (SEp) and the limited extent of the toroidal domain (S2pi). In this section we compute the averages defined in section 2.2 and Appendix A for all simulations between t 1 and t 2 = T end (see tab 1) with a ∆t = 0.16T in .
We start our analysis by confirming that all of our simulations, except for simulation SB2, reach a final steady state quantitatively similar to the one described in section 3.1: a turbulent disk, an atmosphere with super-sonic accretion and a super-fast wind . Therefore we will address the weak-field simulations first, while the steady state of simulation SB2 as well as its properties will be described at the end of this section.

Weak-field simulations: S2pi, SB3, SEp
We start by validating the wedge approximation, ∆ϕ < 2π, for global simulations, S2pi harbours no significant differences when compared to SB4. Indeed, since S2pi and SB4 are identical with respect to the properties of their vertical structure, we will ignore S2pi while we discuss such properties.
We confirm that similarly to simulation SB4, all weak-field simulations reach a equilibrium within the disk and a magnetostatic equilibrium within laminar and turbulent atmospheres (section 3.3). The role of the turbulent magnetic pressure is as essential for the equilibrium and for the transitions (section 3.3) within those simulations as it is for SB4. The turbulence, in the same ways as for SB4, is also a consequence of strong field compressible MRI. This has been verified by checking the growth rate achieved within the turbulent regions (disk and atmosphere) using the procedure detailed in section 3.5. To quantify the evolution of the turbulent stratification as a function of the different parameters we compare the heights and extents of the turbulent disk, the laminar atmosphere and the turbulent atmosphere for the different simulations. We quantify the extent and height of the region by looking at cos ψ the angle between the poloidal velocity and the poloidal magnetic field (section 3.1). Then, by comparing the definition of the vertical domains (TA and LA) given by cos ψ with the one given by the ratio between the turbulent and laminar torques, |α tu /α la |, we cross check the validity of the domains.
We summarise our findings in Fig. (14), where we show the turbulent regions (disk and atmosphere) as the shaded regions for the different simulations. We see that when the magnetic field increases (from SB4 to SB3) the disk region decreases in size. This is because the quenching of the MRI is only possible when β p ∼ 1 (section 3.5). Hence, starting at a lower β mid limits the vertical extent of the disk and push the transition to the laminar atmosphere at lower altitudes.
The altitude of the jet launching surface is tightly related to the density profile. In Fig. (15) we show the density profiles of the different simulations normalised to their value at the disk mid-plane. We see that the simulations SB4 and S2pi converge to a similar density configuration. Moreover, we observe in Fig. (15) that as β ini decreases the density profile becomes shallower and as a consequences more mass is fed into the TA and the outflow. This is possibly related to the enhanced mass loading fulfilled by the turbulent magnetic pressure. Indeed, decreasing β ini leads to a decrease of β mid which increases the turbulence strength (Salvesen et al. 2016;Mishra et al. 2020). This increase in turbulent strength with the decrease of the plasma beta is also verified in our simulations.
From Fig. (14,15), we see that, as is expected when the the disk geometrical thickness decreases, the vertical extent of the turbulent disk decreases. However, it is striking that the turbulent atmosphere stays at roughly the same height compared to SB4. This shows that the turbulent atmosphere and its properties mostly depend on the magnetisation and not on the disk geometrical thickness. However, the density profile is affected, becoming less dense within the turbulent atmosphere. It is not clear if an even smaller could lead to the disappearance of the turbulent and laminar atmosphere.
The initial magnetic field as well as the disk geometrical thickness influence the transport of mass, angular momentum and magnetic field. To understand the impact of the different parameters on accretion, we compute the mass weighted accretion velocity, v acc , from Eq. (21) and then we average it between t a = 318T in and t b = 955T in with temporal resolution of ∆t = 0.8T in for all simulations.
We summarise the values in Tab.
(2). We see that, as the magnetic field increases, v acc also increases. Since the mean radial velocity in the accreting atmosphere does not vary a lot between simulations ( u r 0.3V K ), the dependency of v acc must be related to the evolution of the density profile. This is confirmed by Fig. (15): when the magnetic field strength increases, the density profile becomes shallower and the accreting atmosphere becomes denser. Similarly, the decrease of v acc as decreases is also a consequence of the density profile. Hence, the mass weighted accretion velocity is tightly related to the vertical density structure.
All simulations feature magnetic field transport towards the inner boundary (section 3.2.1 and Appendix D). To quantify the magnetic flux transport we compute v Ψ defined in Eq. (23) for all our simulations and then we time average it in the same manner as v acc . We show v Ψ , normalised to the Keplerian velocity, as a function of radius in Fig. (16) (left). We confirm that v Ψ follows a radial scaling close to the Keplerian velocity. It is clear, that as the β ini decreases v Ψ /V K increases. Indeed, v Ψ /V K seems to be proportional to β mid −1 . Moreover, v Ψ /V K strongly decreases when we decrease the disk geometrical thickness.
In Fig. (16) (right) we show v Ψ /v acc as a function of the radial coordinate. We observe that, except for the inner regions of SB2, all simulations follow approximately the same dependency, This is remarkable and must be related to the evolution of the vertical structure since, as explained above, the evolution of v acc is tightly related to it. In Tab.
(2) we show the accretion time scale as well as the magnetic field transport time scale computed from their respective velocities. In all simulation the the secular time scales t Ψ and t acc remain longer than the dynamical time scale T K . Nonetheless, a word of caution is appropriate since for simulation SB2 the accretion time scale t acc is of the order of the dynamical time scale, while t Ψ is comparable for the outer regions of SB2.
1.1 × 10 −3 144 160 10 4 SB3 1 × 10 −2 16 31 10 3 SEp 4 × 10 −4 400 1.1 × 10 3 10 4 SB2 (R > 8) 4 × 10 −2 4 5 4 × 10 2 SB2 (R < 4) 2 × 10 −1 1 265 1  To understand the impact of the different parameters on the outflows, we summarise in Tab. (3) the values of the MHD invariants, defined in section 3.4, computed from a field line originating at R 0 = 6 for the different simulations. We find that the invariants do not vary from one simulation to another. This may seem surprising. However, it should be kept in mind that the wind invariants depend on the flow dimensionless properties ( u r /V K , u ϕ /V K , B r / B θ , B r / B ϕ ) at the wind launching point, i.e. the top of the turbulent atmosphere. Since these dimensionless properties are similar between all weak-field simulations, the invariants should naturally be independent of β mid .
Ring-like structures are also found in all weak-field simulations. They are always a consequence of converging radial mass flows. In general the rings in the weak-field simulations feature a similar time-scale of self-organisation, t so ∼ 200/Ω K (R) = 32T K (R), consistent with the one measured in SB4. However, there are some quantitative differences: -The ring-like structures in run SB3 are accreted toward the inner region a lot faster than for the fiducial run SB4, v ring 3 × 10 −3 V K (R). This velocity is around a factor of 3 smaller than the local v acc (R) in this simulation. Moreover, the accreted rings are sufficiently far away (R ∼ 4 and R ∼ 6) from the inner boundary for us to believe that radial migration is not a bias of the inner boundary, in contrast with SB4. The ratios between the maxima and minima of density are very similar to the ones in the fiducial run, Σ max /Σ min 2. .
-The ring-like structures in run S2pi are less pronounced than in run SB4, Σ max /Σ min 1.5. This may be related to the limited toroidal extent of simulation SB4. Similarly to SB4, the radial migration of the ring-like structures for S2pi is inexistent or negligible.
-The ring-like structures in run SEp are very similar to the ones found in run SB4, the ratios between minima and maxima are similar, Σ max /Σ min 2. However, run SEp features Article number, page 17 of 27 A&A proofs: manuscript no. JoWinds many more ring-like structures, by about a factor 2. This is consistent with a viscous like instability bounded at small scale by the disk scale-height, in which the spacing between ring-like structures is expected to scale like h. Finally, radial migration of ring-like structures is also negligible in this simulation.
The Rossby Wave Instability (RWI) is known to produce nonaxisymmetric structures, like Rossby vortices or spiral shocks which smear out local density maxima. The RWI is triggered 2 when the inverse of the potential vorticity has a local maximum (Lovelace et al. 1999). To check whether our ring-like structures are unstable to the RWI, we show the inverse of potential vorticity as a function of radii in Fig. (17) for simulations S2pi and SB4. We find that for SB4, the local maxima of the inverse of the vorticity potential are well correlated with the local maxima of the column density profile. This correlation is less clear for run S2pi, as the ring-like structures are less pronounced. One would therefore conclude that SB4 rings should be RWI unstable and therefore develop nonaxisymmertic structures. However, no clear spiral structures are visible in SB4, while non-axisymmetric features similar to the ones shown by Mishra et al. (2020) are appreciable in S2pi. It is therefore temping to think that the RWI in SB4 is stabilised by some mechanism. One possibility is that the RWI is stabilised because of turbulent diffusion which damps RWI modes. Another one is that the limited azimuthal extension of SB4 (∆φ = π/2) filters out some of the most destructive low-m modes excited by the RWI. The presence of strong spiral waves in S2pi is in favour of the second hypothesis, but a proper study of the ring stability is needed to ascertain this conclusion, which is outside the scope of this work.

Strong field simulation: SB2
The main difference between SB2 and all weak-field simulations is that SB2 takes a much longer time to reach a reasonable steady state, the whole disk being subject to a drastic reorganization of the vertical magnetic field. At t end = 1910T in , only the region below R 7 has reached a steady state, while the outer regions continue to struggle readjusting the magnetic field distribution. The final stage of our simulation shows an internal quasi-stationary highly magnetized region with β mid 1 and an outer evolving region with β mid ∼ β ini . Hereafter, we define R J as the end of the equipartition field region.
The reorganisation as well as the lack of stationarity in the outer regions are a consequence of the enormous v Ψ in the outer regions. As shown in Fig.(16), v Ψ follows the same scaling (Eq.40) as for the weak field simulations, driving in this case a much faster field advection. This magnetic flux accumulation reaches a saturation in the inner region with R J 7, where turbulent field diffusion balances advection leading to v Ψ → 0. Because of the presence in our simulation of an important magnetic flux, this transition radius is increasing in time as more magnetic field is being brought in. This can be clearly seen in the lower left panel of Fig. D.1, where we measure a radial drift speed of R J to beṘ J ∼ 10 −3 V K (R in ). Note that this increase in R J is highly 2 Where for simplicity, we neglect the entropy profile. In practice, it will not change the location of the local maxima of F (r) due to our prescribed cooling.  Fig.(14). The quantities are averaged between t 1 = 1719T in and t 2 = 1910T in with a resolution of ∆t = 0.16T in .
. time variable, involving bursts and dramatic changes in the disk within an extended transition region beyond R J . Describing in details this transition region is beyond the scope of this work (see discussion next Section). We note however that the saturation value of the magnetic field reached in the inner steady zone is near equipartition. For a larger magnetic field strength, Ferreira & Pelletier (1995) showed that no vertical equilibrium could be found anymore because of the overwhelming magnetic compression. This is consistent with our own findings and advocates for a maximal mean magnetic field value set by the vertical equilibrium.
In both radially distinct regions the turbulent atmosphere collapses into the disk, following the same trend as the weak field simulations (Fig. 14). This leads to the disappearance of the laminar atmosphere and the formation of a thicker disk when compared to SB3. The merging of the turbulent disk with the turbulent atmosphere is directly related to the density profile (see Fig. 15). As described in the previous section, the turbulent magnetic pressure determines the steepness of the density profile. Thus, as the disk magnetization β mid −1 increases the amount of mass lifted by the turbulent magnetic pressure increases, leading to a massive turbulent atmosphere that merges with the disk. Figure (18) shows a zoomed-in outcome of our SB2 simulation, averaged between t 1 = 1719T in and t 2 = 1910T in . The shape of the poloidal magnetic field lines (black solid lines) anchored in the inner region is clearly different from that seen in weak field simulations (Fig.2). Its structure is remarkably similar to that originally invoked in Blandford & Payne (1982) and later computed by Ferreira (1997). The background color is the density and the dotted black lines delineate the turbulent disk (as in Fig. 14). The turbulence is consistent with the strong field compressible MRI (section 3.5) and leads to a turbulent disk region in magnetostatic equilibrium, much alike the turbulent atmosphere of weak-field simulations. Below R J ∼ 7, accretion occurs within the turbulent disk and the mass weighted accretion speed is supersonic. The disk rotation is near keplerian with u ϕ 0.8v K , namely a deviation of the order of as expected in near equipartition disks. A cold outflow is launched from the disk surface at z 3.5h, becoming soon a super Alfvénic (red dotted curve) and then super fast magnetosonic (red dashed curve) jet.
Overall, this inner accretion-ejection structure shares most properties with the Jet Emitting Disk model (Ferreira 1997). The main difference is the amount of ejected mass as the outflows obtained here are quite massive (Fig.15). As discussed above, this is a consequence of the turbulence itself since the turbulent magnetic pressure plays a major role in the magnetostatic vertical equilibrium and launching of the outflow. Finally, we report the lack of ring-like structures in our strong field simulation. This is presumably due to the fact that accretion is now too fast, t acc t so ,with the t so as obtained from the weak field simulations. Mass has no time to accumulate locally and lead to the formation of rings and gaps.

Discussion and conclusion
We performed 3D ideal MHD simulations of global accretion disks threaded by a large scale vertical magnetic field. Our simulations evolved for more than thousand orbits at the inner radius, allowing our accretion-ejection configuration to achieve a reasonable steady-state. We found that the disk structure is drastically dependent on its magnetization β mid −1 . (1) At disk magnetizations lower than 10 −3 , a non trivial vertical stratification sets in on a scale that can go up to z ∼ 10h, where an outflow that becomes super fast magnetosonic is launched. The disk is highly turbulent, mostly in hydrostatic equilibrium and the presence of an important magnetic turbulent pressure allows to lift up a strong amount of mass into a laminar ideal MHD atmosphere. This lifted mass is basically falling in towards the central object, dragging in the magnetic field which, in turn, provides a torque that transfers its angular momentum back to the underlying disk. When the magnetic field achieves equipartition in these upper layers, MRI is re-ignited and drives a turbulence. In this levitating turbulent atmosphere, transsonic accretion is achieved mostly through the turbulent torque and a small fraction of the mass is then loaded into the open field lines where a magnetically driven cold wind is launched. In contrast with shearing box simulations, the bipolar wind configuration is found to be fairly symmetric. We report also the formation in the disk of rings and gaps, on time scales shorter than the accretion time scale.
(2) The ambient magnetic field is found to be always dragged inwards in the disk, at a velocity which increases with the disk magnetization. The disk vertical structure is also seen to be deeply affected, with a progressive decrease in height of the disk turbulent atmosphere. However, the outflow properties, as described by the usual MHD invariants, remain remarkably constant. This shows that the outflow is determined by the physical conditions at the disk upper layers, which remain mostly constant, rather than by the disk mid-plane properties.
(3) However, beyond a threshold on the disk magnetization, located between 10 −3 and 10 −2 , the global accretion-ejection configuration undergoes a drastic readjustment. The magnetic field is being accumulated into the central regions until a global equilibrium is achieved involving both the disk and its magnetosphere. The inner disk reaches a steady-state (balance between field advection and diffusion) when its magnetization achieves unity and no more gaps and rings are found. Despite this strong field regime, the disk is turbulent and drives a super fast-magnetosonic outflow right from its surface z ∼ 3h. The size of this inner region keeps on increasing in time, as more magnetic flux is being added from the outer regions. The mass-weighted accretion speed is supersonic in the inner region. The transition region, where the disk material must become transsonic, is quite complex and its investigation is postponed for future work.
Some of our results have already been obtained or discussed in the literature. Previous numerical works (Zhu & Stone 2018;Mishra et al. 2020) found the same stationary states as in our simulations. However, these authors did not analyse the appearance of the second upper turbulent zone and how its properties are related to the vertical structure of the system. Another major difference is related to the magnetic field advection. Zhu & Stone (2018) report a small amount of magnetic flux accumulation through their inner boundary albeit without providing a clear analysis. We show that the field advection speed v Ψ follows a Keplerian scaling and increases with the disk magnetization. The strong field simulations of Mishra et al. (2020) are somewhat different as they do not detect considerable magnetic flux evolution. This is in strong contrast with our own results. Our guess is that it may be a consequence of their disk being two times thinner than ours ( = 0.05). Indeed, Fig. (16) clearly shows that the geometrical thickness has a tremendous impact on the magnetic field advection speed (see e.g. Lubow et al. (1994)). This is an interesting aspect with potentially strong astrophysical consequences and deserves further investigations.
One of the main results is the realisation of the utmost importance of the MRI in the strong field regime (Kim & Ostriker 2000) for accretion-ejection structures. While it is widely believed that the MRI can only exist in high β plasmas, Kim & Ostriker (2000) have shown that the MRI exist for arbitrary β provided that arbitrary long wavelength perturbations are allowed. In the β < 1 regime, the MRI can be quenched when the toroidal field gets too strong. We have shown that the turbulent and laminar regions of the atmosphere, which is a β < 1 region, is naturally explained by the MRI in this peculiar regime.
Magnetic field dragging seems to be a generic property of ideal MHD accretion disks threaded by a large scale magnetic field (provided that ∼ 0.1 or larger). Indeed, it is also reported in GRMHD simulations of accretion disks around a Kerr black hole (e.g. McKinney et al. (2012); White et al. (2019) and references therein, see also Liska et al. (2020)). In these simulations, the magnetic field is seen to be gradually advected and concentrated into the central region, leading to the build up of what has been termed a Magnetically Choked Accretion Flow (McKinney et al. 2012). As in our SB2 simulation, magnetic field dragging stops when the vertical field reaches an equipartition value, defining a transition radius with an outer disk at lower magnetization. As more magnetic flux is being advected, a violent relaxation occurs at this radius involving a magnetic Rayleigh-Taylor instability. It leads to the expulsion of the magnetic flux excess and a progressive increase in time of the transition radius, as long as some magnetic flux remains available in the simulation (McKinney et al. 2012). This situation is quite nicely consistent with the hybrid disk configuration proposed for transition disks (Combet & Ferreira 2008;Wang & Goodman 2017) and around black holes . Within this framework, an inner Jet Emitting Disk (JED, Ferreira (1997)) is established until a transition radius R J , beyond which an outer Standard Accretion Disk (SAD, Shakura & Sunyaev (1973)) is settled. Because of its supersonic accretion speed associated with the launching of powerful jets, a JED provides the physical conditions allowing to explain most observational properties of X-ray Binaries (see Marcel et al. (2019) and references therein).
We further note that our weakly magnetized disks (i.e. with a magnetization smaller than 10 −3 ) provide a highly subsonic mass-weighted accretion speed as expected in a SAD. However they differ from the SAD model in two important aspects: (i) the presence of outflows and (ii) the existence of an upper turbulent accreting zone.
As stated above, the properties of the super fastmagnetosonic jets found in our simulations are remarkably independent of the disk magnetization. Indeed, our outflows exhibit a magnetic lever arm λ ∼ 4 − 5 valid for 4 orders in magnitude in β mid −1 (see Tab. 3). The best observational constraints on jet kinematics are provided by proto-planetary systems, where the magnetic lever arm has been recently characterized in molecular observations of several outflows. In HH212 Tabone et al. (2017) measure λ 5.5 for a jet/wind launching zone extending from R 0,1 = 0.1 au to R 0,2 = 40 au. In HH30, Louvet et al. (2018) estimate λ 1.6 and a launching radius R 0 1.5 au, while de Valon et al. (2020) measure λ 1.6 and R 0 2 au for DG Tau B. This discrepancy in λ with our simulations could be explained by several aspects. If some additional heat is deposited at the base of the outflow, the disk ejection efficiency is enhanced which leads to the decrease of the magnetic lever arm (Casse & Ferreira 2000). Such heating could be a natural outcome of irradiation from the central object and/or local turbulent heating. Moreover, these molecular outflows are emitted from regions beyond 1 au. At those large distances, disk ionization becomes an issue and non ideal MHD effects come into play. The inclusion of ambipolar diffusion and Hall effect, which are believed to be present in the outer regions of proto-planetary disks, may change the outflow properties (Béthune et al. 2017).
The second important difference with the Shakura-Sunyaev SAD model is the existence of a turbulent accreting zone levitating above a weakly magnetized disk. Our simulations show that the disk has almost no accretion taking place within the densest body and an elevated supersonic accretion zone. This is somewhat reminiscent of the layered accretion model invoked in the outer regions of proto-planetary disks (Gammie 1996;Fleming & Stone 2003;Bai & Stone 2013). However, the physical reason is completely different. In our case, the disk is fully ionized and turbulent whereas in the dead zone picture non-ideal effects quench turbulence and lead to the formation of a MRI stable ("dead") zone in the disk. The interesting part is that, in both cases and as long as the vertical field magnetic is weak, accretion takes place at the disk "surface". This is a feature that is drastically different from the conventional view of Shakura-Sunyaev accretion disks. Moreover, this may lead to direct observational consequences on the emitted spectra of weakly magnetized ac-cretion disks. Indeed, mass accretes across the magnetic field lines and requires thereby some magnetic (turbulent) diffusivity (see however Zhu et al. 2020). As a consequence, a fraction of the accretion energy must also be dissipated locally in this lowdensity region. In the context of compact objects, we expect the production of an optically thin emission, usually referred to as a coronal emission. This "elevated accretion" disk configuration could thus provide an explanation for the "warm corona" component often seen in soft X rays in some active galactic nuclei (Done et al. 2012;Petrucci et al. 2018).
It is important to note that, to make the problem numerically tractable we have decided to ignore the effects of the energy equation. We did this by forcing a locally isothermal structure with a fixed disk scale, = 0.1 or = 0.05. Optically thick disks in AGN and X-ray binaries are expected to be geometrically thin with = 0.01 or less, which is incompatible with our simulation setup. Nonetheless, we believe that the accreting atmosphere should perdure at lower , as long as the disk magnetization is small. To properly study the impact of the energy equation, one would need to include the effects of radiative transfer into our simulations, which is clearly beyond the scope the present paper. Irradiation from the central object or disk regions (as in Begelman et al. 1983;Higginbottom & Proga 2015;Gressel et al. 2020, for instance) could lead to a quantitative modification of our picture. Indeed, including heat deposition by external irradiation will probably lead to more massive outflows (Casse & Ferreira 2000) and modify the vertical stratification leading to the turbulent accreting layer. The study of such thermo-magnetic winds is postponed for future work.
The parameter space explored in our simulations overlaps with the one explored by Jacquemin-Ide et al. (2019) using selfsimilar solutions. The outflow invariants computed in our simulations (Tab. 3) are found consistent with the borders of the parameter space they derive semi-analytically. However, while in our simulations λ ∼ 4 − 5 is roughly constant, they obtain values ranging from 2 to 100. As shown in their Figs.(4) and (5), the magnetic lever arm is independent of the disk magnetization and is strongly affected by the mass load κ. The reason for their wide range in λ stems therefore only from their capability to produce outflows with a mass load κ ranging from a few 10 −3 to almost unity, while our 3D simulations always provide κ ∼ 0.2 − 0.7. This discrepancy is related to our different vertical stratification. In our simulations, the levitating turbulent atmosphere acts as a buffer between the proper turbulent disk and the ideal MHD outflow. As discussed in Section 3.3, the existence of a dense, first laminar then turbulent, atmosphere is a consequence of the disk turbulence itself. Note that this turbulent accreting atmosphere is actually missing in all works where turbulence is not self-consistently computed and must therefore be prescribed: self-similar studies (Ferreira 1997; Jacquemin-Ide et al. 2019), 2.5D MHD simulations done with alpha prescriptions (Murphy et al. 2010;Stepanovs & Fendt 2016) and other semi-analytical models (Guilet & Ogilvie 2013). Our guess is that including the complex vertical stratification should allow to recover with these approaches an outflow mass load consistent with full 3D simulations. This requires however the use of numerically "educated" profiles for all relevant turbulent effects. This is postponed for future work.
Finally, weakly magnetized disks are subject to a secular instability which leads to the formation of long-lived ring-like structures in the disk surface density distribution, the magnetic field being concentrated in the gaps. These ring-like structures are also visible in the surface density profiles of Zhu & Stone (2018) and, like for us, they also seem to increase in number when the disk geometrical thickness is decreased. This is consistent with a viscous like instability bounded at small scale by the disk scale-height, for which the spacing between ring-like structures is expected to scale like h. On the other hand, we showed that this global self-organisation is unlikely due to the wind instability proposed by Riols & Lesur (2019). Yet, it remains unclear if it is purely a viscous instability, as the turbulent viscosity scaling α tu ( β mid ) is not steep enough to develop a viscous instability. When the disk magnetization increases above the threshold, these ring-like structures are smeared out and disappear as the accretion speed becomes supersonic. Our suspicion is that the ring-formation mechanism and the magnetic field transport could be related to each other. Finally, we show that the limited toroidal extent of our simulation SB4 and SB3 could have an impact on the long term evolution of the ring-like structures. Indeed, it is not clear if the rings are destroyed by non-axissymmetric instabilities, like the RWI.
In any case, our results advocate for the existence of a dichotomy in ideal MHD accretion disks. In disks where the magnetization remains below a threshold value, around 10 −3 − 10 −2 , ring-like structures should be present, and absent in disks with equipartition magnetic fields. It could thus be an observationally convenient way to probe the existence and strength of large scale magnetic fields in proto-planetary disks (see for instance de Valon et al. (2020)). access to the secular radial mass flux, indicating that mass is flowing on average outwards and that the disk is decreting. Interestingly, the flux budget indicates that ∆M = −6.19 < 0, implying that the disk portion is not in a stationary state and is loosing mass. This related to the appearance of ring-like structures within the disk region, a form of secular self-organisation. The domain, between r 1 = 8 and r 2 = 14, is located on a two distinct gaps, one at R = 9.5 and another at R = 15. The formation of ring-like structures is discussed in section 4. Turning now to the regions above the disk, we find that the laminar atmosphere is accreting, but less efficiently than the turbulent atmosphere, despite the fact that both regions receive approximately the same amount of mass from the outer radial boundary. This trend can be understood from the shape of the stream lines (Fig. 2), since the mass transported within the laminar atmosphere inevitably ends up being accreted in the turbulent atmosphere at a smaller radius. This is confirmed by Fig. B.2 which shows that the turbulent atmosphere is being fed mass by the laminar atmosphere. Eventually, most of the mass received by the turbulent atmosphere is radially accreted, the mass ejected being almost negligible. Note finally that the mass budget ∆M for the turbulent and laminar atmospheres is approximately null, corresponding to a steady system.
Overall, our analysis suggests that mass is being channeled from a decreting turbulent disk, through the laminar atmosphere along the field lines, and is accreted once in the turbulent atmosphere. To understand this exotic configuration, we need to look at the angular momentum transport in our domain.

Appendix B.2: Angular momentum flux
We integrate the equation of conservation of angular momentum with respect to time and space, as we did for the conservation of mass, to find where the fluxes definition can be found in appendix B. The lefthand side denotes angular momentum transported inwards by the accreting flow, while the right-hand side measure the angular momentum flux due to the torques (turbulent and laminar). The angular momentum fluxes are quasi-stationary when compared to the mass fluxes. Hence, we skip the examination of the time variability of these fluxes and average them across the whole time domain (t 1 = 5800/Ω K (R in ) to t 2 = 7000/Ω K (R in )). We summarize their values in figure B.3.
We can see that the angular momentum is transported radially outwards by the torques in all of the regions. Surprisingly, in the latitudinal direction, angular momentum is transported downwards: from the atmosphere to the disk. The disk is therefore azimuthally accelerated by the atmosphere, as already seen by Zhu & Stone (2018). Overall, this latitudinal flux of angular momentum implies that the disk is gaining angular momen-tum, hence decretion must occur 3 . Overall, one concludes that the latitudinal transport of angular momentum from the laminar atmosphere to the disk leads to the mass decretion within the disk found in the previous section. Higher up, the laminar atmosphere The different fluxes are defined in appendix B and they are integrated between t 1 and t 2 . The arrows represent the direction of the fluxes. In this figure we are only representing the fluxes that emerge from the Reynolds and magnetic stresses F L,r 2 , F L,r 1 , F L,θ 1 , F L,θ 2 not the flux of angular momentum though accretion. However, the sum of these fluxes (taking in to account their direction) leads to the value of the angular momentum being transported by the accretion. gets its angular momentum from the inner radius and the turbulent atmosphere. Indeed, the angular momentum is transported radially outwards. However, it is important to note that in this region the transport is dominated by laminar stresses (Fig. (3).
Finally, below the wind basis, the turbulent atmosphere looses angular momentum from all sides because of the torques. The wind itself however carries less than 10% of the angular shown in the lower left panel in Fig. (D.1), which is associated with the SB2 strong field simulation, corresponds in fact to a region where the disk magnetization β mid −1 is around unity.