Issue 
A&A
Volume 647, March 2021



Article Number  A192  
Number of page(s)  26  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202039322  
Published online  31 March 2021 
Magnetic outflows from turbulent accretion disks
I. Vertical structure and secular evolution
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: jonatan.jacquemin@univgrenoblealpes.fr
Received:
1
September
2020
Accepted:
30
November
2020
Context. Astrophysical disks are likely embedded in an ambient vertical magnetic field generated by its environment. This ambient field is known to drive magnetorotational turbulence in the disk bulk, but it is also responsible for launching magnetised outflows at the origin of astrophysical jets. Yet, the interplay between turbulence and outflows is not understood. In particular, the vertical structure and longterm (secular) evolution of such a system lack quantitative predictions. It is, nevertheless, this secular evolution which is proposed to explain time variability in many accreting systems such as FuOr, Xray binaries, and novae like systems.
Aims. We seek to constraint the structure and longterm evolution of turbulent astrophysical disks subject to magnetised outflows in the nonrelativistic regime. More specifically we aim to characterise the mechanism driving accretion, the dynamics of the disk atmosphere, the role played by the outflow, and the longterm evolution of mass and magnetic flux distributions.
Methods. We computed and analysed global 3D ideal magnetohydrynamic (MHD) simulations of an accretion disk threaded by a largescale magnetic field. We measured the turbulent state of the system by Reynolds averaging the ideal MHD equations and evaluate the role of the turbulent terms in the equilibrium of the system. We then computed the transport of mass, angular momentum, and magnetic fields in the disk to characterise its secular evolution. Finally, we performed a parameter exploration survey in order to characterise how the transport properties depend on the disk properties.
Results. We find that weakly magnetised disks drive jets that carry a small fraction of the disk angular momentum away. The massweighted accretion speed remains subsonic, although there is always an upper turbulent atmospheric region where transsonic accretion takes place. We show that this turbulence is driven by a strongly magnetised version of the magnetorotational instability. The internal disk structure therefore appears drastically different from the conventional hydrostatic picture. We expect that the turbulent atmosphere region will lead to nonthermal features in the emission spectra from compact objects. In addition, we show that the disk is subject to a secular viscoustype instability, which leads to the formation of longlived ringlike structures in the disk surface density distribution. This instability is likely connected to the magnetic field transport. Finally, we show that for all of the parameters explored, the ambient magnetic field is always dragged inward in the disk at a velocity which increases with the disk magnetisation. Beyond a threshold on the latter, the disk undergoes a profound radial readjustment. It leads to the formation of an inner accretionejection region with a supersonic massweighted accretion speed and where the magnetic field distribution becomes steady and reaches a magnitude near equipartition with the thermal pressure. This inner structure shares many properties with the jet emitting disk model. Overall, these results pave the way for quantitative selfconsistent secular disk models.
Key words: accretion, accretion disks / magnetohydrodynamics (MHD) / turbulence / protoplanetary disks / ISM: jets and outflows / Xrays: binaries
© J. JacqueminIde et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The emission of jetlike 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. Jetlike outflows are observed being emitted from active galactic nuclei (AGN) and quasar (Merloni et al. 2003 and references therein; Blandford et al. 2019), Xray binaries (Mirabel & Rodríguez 1999; Corbel et al. 2000; Gallo et al. 2003, 2005), and possibly also novalike 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 jetlike 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; Ferreira et al. 2006a). 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 (Cicone et al. 2014), Xray binaries (Tetarenko et al. 2018; JiménezIbarra et al. 2019), and protoplanetary disks (Tabone et al. 2017; de Valon et al. 2020). 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 with outflows which carry angular momentum vertically away from the disk (Blandford & Payne 1982). In the simple hypothesis of a Keplerian disk embedded in a largescale vertical magnetic field, both mechanisms naturally emerge. First, the seminal work of Balbus & Hawley (1991) showed that the magnetorotationnal instability (MRI) is triggered, leading to vigorous magnetohydrodynamic (MHD) turbulence and radial transport of angular momentum. Second, large scale magnetic fields spontaneously lead to the launching of an outflow (Wardle & Königl 1993; Ferreira & Pelletier 1995; Konigl & Pudritz 2000; Ferreira et al. 2002; Pudritz et al. 2007), which carries angular momentum away 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 (Lesur et al. 2013).
The launching of outflows has been studied using selfsimilar 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 azimuthallyaveraged 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; JacqueminIde 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 JacqueminIde et al. 2019 for a discussion).
MRIdriven 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 magneticallydriven 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 magneticallydriven 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 magneticallydriven 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 windemitting disk. In these models, accretion typically occurs in the disk bulk and angular momentum escapes through a magnetised 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; JacqueminIde et al. 2019), where h is the disk pressure scaleheight, none of these models produce accretion in the disk atmosphere at z ∼ 10 h, 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 magnetisation, 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 magnetised outflow. It is this secular evolution which then leads to interesting observational evidence, like for example eruptions in Xray binaries (Ferreira et al. 2006b; Marcel et al. 2018, 2019) or dwarf novae (Scepi et al. 2019, 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 counterbalanced 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 selfconsistent. 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 MRIdriven turbulence and magnetised 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 ringlike structures. In Sect. 5 we vary the field strength and disk aspect ratio to study their impact on the processes discussed in Sects. 3 and 4. We finally put in perspective our results and discuss their astrophysical implications in Sect. 6.
2. Physical and numerical framework
2.1. 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 ionised to be described in the ideal MHD approximation. In this limit the equations describing a magnetised 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 Eqs. (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 and the disk geometrical scale height h. Finally, we define the poloidal velocity and magnetic field as
2.2. Numerical method
To solve Eqs. (1)–(5) on a spherical grid, we use the conservative Godunovtype 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 RungeKutta 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 Table 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.
Parameters, grid spacing and radial extension of the simulations discussed in the paper.
We use periodic boundary conditions in the toroidal direction and zerogradient in the radial direction (mass flow into the numerical domain is forbidden). As for the latitudinal direction we implement the same polar (zaxis) 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 midplane 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 Eqs. (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 midplane defined above. We get
and
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
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 selfsimilar 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:
where
To ensure that ∇⋅B = 0 at t = 0, we initialise the magnetic field using the magnetic potential A. We can solve the equation B = ∇×A to find
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 magnetisations close to the polar boundary we implement a density floor in our simulation
where ρ_{fl, 0} = 10^{−7} and a_{fl} ∼ 1. Finally, to minimise integration time, we enforce a maximal Alfvénic velocity V_{A, max} = V_{K}(R_{in}), where 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 Table 1. This choice of resolution allow us to resolve the MRI within the disk thanks to roughly 16 grid points per vertical scale height.
Fig. 1. Computational grid. 
We integrate our simulations up to T_{end}, see Table 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).
2.3. Diagnostics and turbulent decomposition
2.3.1. 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 (Table 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 selfsimilar solutions (see JacqueminIde et al. 2019). Hence, they allow for a direct comparison to said models and provide guidance for turbulence closures. Let us emphasise however that ⟨δX⟩_{ρ} ≠ 0, which leads to additional complications when dealing with high order correlation functions. Finally, we 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 (Table 1).
2.3.2. Turbulent decomposition
In order to study the impact of turbulence on the vertical structure of the accretionejection system, we use the averages defined above to rewrite the ideal MHD equations as a set of Reynoldsaveraged 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 Reynoldsaveraged 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 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 and 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.
3. Simulation SB4, the fiducial case
3.1. Global picture
We first present a global picture of our fiducial simulations. We show in Fig. 2 the mean poloidal stream lines and field lines, the Alfvén surface and the fast magnetosonic surface of the fiducial β_{ini} = 10^{4} run. We define the Alfvén surface as the location where ⟨u_{p}⟩ = ⟨V_{ap}⟩, where and the fast magnetosonic surface where ⟨u_{p}⟩ = ⟨V_{fm}⟩, with , having defined the total Alfvénic velocity .
Fig. 2. Top, left: gas density and mean poloidal stream lines. The red dotted line corresponds to the Alfvénic surface, and the red dashed line corresponds to the fast magnetosonic surface. The colour of the poloidal stream lines correspond to the logarithm of their magnitude normalised to the sound speed. Top, right: RB_{φ} normalised to B_{i}R_{in}; and mean poloidal field lines. The grey square corresponds to the zoomed in region the bottom figure. Bottom: same as top but zoomed in the greyed region. The black dashed line indicates the surface where ⟨β⟩ = 8π⟨P⟩/⟨B⟩^{2} = 1. 
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 magnetosoni (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 selfsimilar jet solutions (Ferreira 1997; JacqueminIde 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 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 Fig. 2:
Firstly, A turbulent disk where the density is maximal and the poloidal flow is disorganised, indicating a turbulent flow. The toroidal component of the magnetic field exhibits similar disorganised structure, while the average poloidal field lines are remarkably straight.
Secondly, An accreting atmosphere where the poloidal flow is laminar and transonic. The field lines in this region seem to be curved by the poloidal flow. In addition, a strong and laminar toroidal field emerges.
Finally, a super fast magnetosonic outflow is ejected from the atmosphere, in the wind region. The poloidal field lines are approximately aligned with the outflow direction.
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 subregions, 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.
Fig. 3. Top: ratio α_{la}/α_{tu} denoting turbulent (red) and laminar (blue) regions. The dotted lines delimit these regions and represent z = 3.5h, z = 6.5h and z = 12.5h. The two different circles denote the radii used for the calculation of the fluxes (r_{1} = 8, r_{2} = 14), see Sects. 3.2.1 and Appendix B. TA and LA correspond for turbulent and laminar atmosphere, respectively. Bottom: cosine of the angle ψ between the mean poloidal velocity and the mean poloidal magnetic field as a function of the latitudinal coordinate. The grey zones correspond to the turbulent regions in the top panel, delimited by the dotted lines. 
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.
Fig. 4. Top: velocity profiles as functions of the latitudinal coordinate, θ, normalised to the local sound speed. Bottom: mean magnetic field profiles normalised to the vertical magnetic field in the disk midplane, ⟨B_{z}⟩(r, θ = π/2) as functions of the latitudinal coordinate. The grey zones correspond to the turbulent regions in Fig. 3. All profiles are radially averaged between r_{1} = 8 and r_{2} = 14 after being normalised. 
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 Figs. 4 and 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 idealMHD 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 magnetosonic (SM) surface defined as ⟨u_{p}⟩ = V_{SM}, where (Fig. 4 (top)).
It is important to note that V_{SM} 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 infrathermal within the laminar atmosphere, ⟨β_{p}⟩ > 1, where
but becomes suprathermal within the turbulent atmosphere, ⟨β_{p}⟩ < 1.
3.2. Transport properties
3.2.1. Mass transport
Using the regions defined in Sect. 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 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 selforganisation of the radial disk structure. Indeed, a closer look at the radial density profiles reveals the formation of ringlike structures. Our chosen disk section turns out to be located on a gap. The mechanism behind the formation of these structures is discussed in Sect. 4.
In the accreting atmosphere the mass is transported radially inwards, the 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.
3.2.2. 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 t_{b} = 955T_{in} with temporal resolution of Δt = 0.8T_{in}. We find
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.
3.2.3. 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 characterise this transport, we introduce the magnetic flux threading the disk as
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_{out}) 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.
Fig. 5. Magnetic flux threading the disk midplane Ψ, computed between the axis and R, as a function of time in innermost orbital units. Every contour line represents a field line. The dashed line shows an example of a fit using Eq. (26) with R_{i} = 4.4, t_{i} = 328T_{in} and . 
Given that magnetic flux accumulates towards the disk centre, 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)
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 t_{b} = 955T_{in}, which provides
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 crosschecked 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 selfsimilar 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.
3.3. 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 normalised to g and gr.
Fig. 6. 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 normalised 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. 
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 magnetised 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, or . 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. (2015), a disk in magnetostatic equilibrium where the magnetic 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 midplane (Fig. 4 bottom), that changes sign within the turbulent atmosphere.
Fig. 7. Comparison of the different pressure terms normalised to the thermal pressure within the disks midplane as functions of the latitudinal coordinate. The shaded areas denote the turbulent regions. 
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 inwards. 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 JacqueminIde 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 Sect. 3.5.
3.4. Super fast magnetosonic wind
In the regions where the flow is approximately laminar and steadystate, it is possible to define a set of MHD invariants which consists of conserved physical quantities along the field lines. These invariants characterise 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 , where R_{SM} is the anchoring radius of the magnetic field line and the SM stands for quantities evaluated at the slow magnetosonic 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 ω_{⋆} = Ω_{⋆}/Ω_{K}(R_{SM}).
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})_{SM}, 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 normalise this invariant with respect to the gravitational energy at the midplane defining , where .
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 Eqs. ((29)–(32)) across the 3 different field lines anchored at the radii R_{SM} ≃ [3.5, 4, 7.5].
Fig. 8. Top: mean poloidal magnetic field lines in the (R,z) plane. The blue lines define the different critical surfaces (see Sect. 3.1) and they are represented as vertical lines in the lower panel. The grey zones determine the turbulent zones defined in Sect. 3.1. The lower zone corresponds to the disk and the upper zone corresponds to the turbulent atmosphere. Bottom: MHD invariants calculated along the field lines of the upper panel as functions of the latitudinal coordinate. The line style of the MHD invariants has a one to one correspondence with the field line where the invariant was calculated. The grey zone corresponds to the end of the turbulent atmosphere. 
We show in Fig. 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. 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 Table 3).
Firstly, 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 selfsimilar models (JacqueminIde et al. 2019).
Secondly, the angular momentum invariant λ ≃ 4—6, so the wind is effectively free (λ > 3/2). This provides 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 (Blandford & Payne 1982). Using this expression we find u_{p∞} ∼ 3R_{SM}Ω_{K}(R_{SM}).
Finally, 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 jetlike 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. Figure 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_{SM}Ω_{K}(R_{SM}), since e ∼ 6.
Fig. 9. Bernoulli invariant and its components for the field line originating at R_{0} = 6 as functions of the latitudinal coordinate θ. All quantities are normalised to . The shaded grey area corresponds to the turbulent atmosphere, the end of the turbulent atmosphere coincides with the SM surfaces. The vertical lines correspond to the different critical surfaces. The vertical dashed line is the fast magnetosonic surface while the vertical dotted line is the Alfvénic surface, consistent with the notation Fig. 8. 
It is also clear from Fig. 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.
JacqueminIde 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 magnetisation
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 Fig. 9, we derive σ ∼ 8 which is compatible with a jet like outflow and is consistent with the observed collimation.
3.5. Origin of turbulence in the atmosphere
As described in Sect. 3.1 the accreting atmosphere is divided into two distinct regions: a turbulent and a laminar atmosphere. This configuration, where turbulence is localised 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 magnetised (⟨β_{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
This condition is verified in the laminar atmosphere (Fig. 4, bottom), indicating that it should be stable to the magnetorotational instability. For the sake of completeness we derive Eq. (35) as well as the dispersion relation for the compressible MRI, Eq. (C.2). 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 Sect. 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.
Fig. 10. 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}). 
As explained in Sect. 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 reignition of turbulence is required for the establishment of a stationary magnetostatic equilibrium.
The reactivation 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.
4. Formation of ring like structures
We discussed in Sect. 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 selforganisation. 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 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 (Sect. 3.1). We show Σ in a space time diagram in the top panel of Fig. 11.
Fig. 11. Space time diagram of the column density Σ (up) and the vertically averaged vertical magnetic field b_{z} (down) as functions of the time in units of orbits at the inner most radius and of the cylindrical radial coordinate. The solid black lines represent the gaps located at R ≃ [6.6, 12.3] while the dashed black lines represent the rings located at R ≃ [4, 9.5, 15]. 
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 ringlike 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 selforganisation 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 inward 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 ringlike structures. To understand the role of the magnetic field in the selforganisation mechanism, we compute the vertically average vertical field and show its space time diagram in Fig. 11, bottom. We find that the vertical magnetic field is accumulating within the gaps, and is anticorrelated with the rings. However, the magnetic structure is not as stable with time as the density structure.
To understand if the selforganisation 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 as the wind mass loss rate. To simplify the comparison of each term, we have normalised Eq. (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, , is essentially negligible, while the radial term, 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}.
Fig. 12. Radial (blue) and vertical (red) contributions to the mass evolution rate in local Ω_{K} units as a function of the cylindrical radii, R. A positive (negative) term leads to the formation of local gap (ring). The solid black lines represent the gaps located at R ≃ [6.6, 12.3] while the dashed black lines represent the rings located at R ≃ [9.5, 15]. We use symmetrical logscale that becomes linear within the interval [5 × 10^{−5}, −5 × 10^{−5}]. 
From Fig. 12 we can probe the time scale, t_{so}, on which this selforganisation 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 organisation time scale with the ones deduced from Sect. 3.2.1 we find t_{so} ≪ t_{acc} ≃ t_{Ψ}. Selforganisation into ringlike structures is not impeded by the global transport mass or magnetic field, it happens on shorter time scales.
As mentioned above, the selforganisation 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, . In the case of a strongly diffusive field, the limit where the field is steady, this is equivalent to , where we define
and
However, in MRI turbulence, the measured scaling is closer to (Salvesen et al. 2016). For the sake of completeness 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 midplane for radii ranging from R = 4 to R = 16. We expect to not only be a function of ⟨β_{mid}⟩, indeed in Fig. 13 we have the same value of for different plasma beta. This is probably a consequence of the dependence of on the local shear, d logΩ/d log R. Indeed, selforganisation induces local deviation from the Keplerian shear, which in turn influences the efficiency of the MRIdriven angular momentum transport (Gressel et al. 2015; Pessah et al. 2007), weaker shear leading to lower values of . It can be seen that the scaling of the turbulent angular momentum transport coefficient, , verifies a scaling closer to −1/2 than −1.
Fig. 13. Vertically (within the disk region) averaged turbulent angular momentum transport coefficient, as a function of the mean poloidal plasma, ⟨β_{mid}⟩, evaluated in the disk midplane for different cylindrical radii, R = [4, 16]. We also show different scaling laws for to compare them with our simulation. 
The −1/2 scaling is inconsistent with a simple viscous type instability, hence the selforganisation of our simulation into a ringlike system is not purely viscous, and must probably be somehow connected to magnetic field transport in the disk.
In conclusion, the fact that the selforganisation is driven by radial motions of matter tends to disfavour a winddriven instability (Riols & Lesur 2019) and tends to favour some kind of viscouslydriven instability (Hawley 2001; Bai & Stone 2014). However, as explained above, a simple viscous instability could not explain the disk instability since the scaling of on ⟨β_{mid}⟩ is not steep enough. Finally, we note that the radial density structure of the outflow is unaffected by the ringlike 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 , based on the evolution of the vertical magnetic field driven by an anisotropic antidiffusivity. However, it is not evident that an anisotropic antidiffusion 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 selforganisation of the disk into ringlike structures as well as the formation of pressure bumps. However, their modelling concentrates on unstratified shearing box simulations with a zeronet flux. This leads to their ringlike structures being short lived, lasting at most 50 local orbits while our inner ringlike 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.
5. Parameter exploration
To have a comprehensive understanding of how the properties described in Sects. 3 and 4 evolve with different fundamental parameters of the system we have produced 4 additional simulations, see Table 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 Sect. 2.2 and Appendix A for all simulations between t_{1} and t_{2} = T_{end} (see Table 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 Sect. 3.1: a turbulent disk, an atmosphere with supersonic accretion and a superfast wind. Therefore we address the weakfield simulations first, while the steady state of simulation SB2 as well as its properties will be described at the end of this section.
5.1. Weakfield 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 ignore S2pi while we discuss such properties.
We confirm that similarly to simulation SB4, all weakfield simulations reach a equilibrium within the disk and a magnetostatic equilibrium within laminar and turbulent atmospheres (Sect. 3.3). The role of the turbulent magnetic pressure is as essential for the equilibrium and for the transitions (Sect. 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 Sect. 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 (Sect. 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 (Sect. 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.
Fig. 14. Extent of the turbulent zones (disk and atmosphere) found in the different simulations, see Table 1. As the magnetic field increases, the turbulent layer goes down and eventually merges with the turbulent disk. 
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 midplane. 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.
Fig. 15. Mean density of the different simulations normalised to its value at the disk midplane, see Table 1. The symbols correspond to the height at which the flow reaches the SM speed. The mean profile are also averaged radially between r = 8 and r = 12. 
From Figs. 14 and 15, we see that, as is expected when 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 Table 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 (Sect. 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.
Fig. 16. Left: field advection velocity as a function of radius normalised to the Keplerian velocity. Right: ratio of field and mass advection velocities as a function of radius. We truncate the advection velocity where we can no longer measure the drift of the magnetic field lines. The truncation radii corresponds to R = 3.5, R = 10, R = 20, for SEp, SB4 and SB3 respectively. 
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 Table 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 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.
To understand the impact of the different parameters on the outflows, we summarise in Table 3 the values of the MHD invariants, defined in Sect. 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, that is the top of the turbulent atmosphere. Since these dimensionless properties are similar between all weakfield simulations, the invariants should naturally be independent of ⟨β_{mid}⟩.
Values of the MHD invariants for all simulations, measured in the upper ‘hemisphere’ for a field line originating at R_{0} = 6.
Ringlike structures are also found in all weakfield simulations. They are always a consequence of converging radial mass flows. In general the rings in the weakfield simulations feature a similar timescale of selforganisation, t_{so} ∼ 200/Ω_{K}(R) = 32T_{K}(R), consistent with the one measured in SB4. However, there are some quantitative differences:

The ringlike structures in run SB3 are accreted towards 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 ringlike 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 ringlike structures for S2pi is inexistent or negligible.

The ringlike 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 many more ringlike structures, by about a factor 2. This is consistent with a viscous like instability bounded at small scale by the disk scaleheight, in which the spacing between ringlike structures is expected to scale like h. Finally, radial migration of ringlike 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 ringlike 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 ringlike 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 nonaxisymmetric 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 lowm 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.
Fig. 17. Inverse of potential vorticity and surface density as functions of the radial coordinate, for simulations SB4 and S2pi, w_{θ} = ∇×⟨u⟩_{θ}.. 
5.2. Strong field simulation: SB2
The main difference between SB2 and all weakfield simulations is that SB2 takes a much longer time to reach a reasonable steady state, the whole disk being subject to a drastic reorganisation 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 quasistationary highly magnetised 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}). This increase in R_{J} is highly 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 magnetisation ⟨β_{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 zoomedin 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 colour 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 (Sect. 3.5) and leads to a turbulent disk region in magnetostatic equilibrium, much alike the turbulent atmosphere of weakfield 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.
Fig. 18. Appearance of the SB2 strong field simulation: logarithm of the density (background colour), mean poloidal magnetic field lines (black solid lines), Alfvén (red dotted line) and fast magnetosonic (red dashed line) critical surfaces. The black dotted lines represent the end of the disk region defined in Fig. 14. The quantities are averaged between t_{1} = 1719T_{in} and t_{2} = 1910T_{in} with a resolution of Δt = 0.16T_{in}. 
Overall, this inner accretionejection 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 ringlike 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.
6. 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 accretionejection configuration to achieve a reasonable steadystate. We found that the disk structure is drastically dependent on its magnetisation ⟨β_{mid}⟩^{−1}.
(1) At disk magnetisations 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 reignited 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 magnetisation. 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 midplane properties.
(3) However, beyond a threshold on the disk magnetisation, located between 10^{−3} and 10^{−2}, the global accretionejection 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 steadystate (balance between field advection and diffusion) when its magnetisation achieves unity and no more gaps and rings are found. Despite this strong field regime, the disk is turbulent and drives a super fastmagnetosonic 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 massweighted 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 magnetisation. 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 accretionejection 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 magnetisation. As more magnetic flux is being advected, a violent relaxation occurs at this radius involving a magnetic RayleighTaylor 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 (Ferreira et al. 2006b). 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 Xray Binaries (see Marcel et al. 2019 and references therein).
We further note that our weakly magnetised disks (i.e. with a magnetisation smaller than 10^{−3}) provide a highly subsonic massweighted 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 magnetisation. Indeed, our outflows exhibit a magnetic lever arm λ ∼ 4−5 valid for 4 orders in magnitude in ⟨β_{mid}⟩^{−1} (see Table 3). The best observational constraints on jet kinematics are provided by protoplanetary systems, where the magnetic lever arm has been recently characterised in molecular observations of several outflows. In HH212 Tabone et al. (2017) measure λ ≃ 5.5 for an outflow 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 local turbulent heating. Moreover, these molecular outflows are emitted from regions beyond 1 au. At those large distances, disk ionisation 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 protoplanetary disks, may change the outflow properties (Béthune et al. 2017).
The second important difference with the ShakuraSunyaev SAD model is the existence of a turbulent accreting zone levitating above a weakly magnetised 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 protoplanetary disks (Gammie 1996; Fleming & Stone 2003; Bai & Stone 2013). However, the physical reason is completely different. In our case, the disk is fully ionised and turbulent whereas in the dead zone picture nonideal 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 ShakuraSunyaev accretion disks. Moreover, this may lead to direct observational consequences on the emitted spectra of weakly magnetised accretion 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 Xray 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 magnetisation 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 thermomagnetic winds is postponed for future work.
The parameter space explored in our simulations overlaps with the one explored by JacqueminIde et al. (2019) using selfsimilar solutions. The outflow invariants computed in our simulations (Table 3) are found consistent with the borders of the parameter space they derive semianalytically. 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 magnetisation 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 Sect. 3.3, the existence of a dense, first laminar then turbulent, atmosphere is a consequence of the disk turbulence itself. This turbulent accreting atmosphere is actually missing in all works where turbulence is not selfconsistently computed and must therefore be prescribed: selfsimilar studies (Ferreira 1997; JacqueminIde et al. 2019), 2.5D MHD simulations done with alpha prescriptions (Murphy et al. 2010; Stepanovs & Fendt 2016) and other semianalytical 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 magnetised disks are subject to a secular instability which leads to the formation of longlived ringlike structures in the disk surface density distribution, the magnetic field being concentrated in the gaps. These ringlike 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 scaleheight, for which the spacing between ringlike structures is expected to scale like h. On the other hand, we showed that this global selforganisation 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 is not steep enough to develop a viscous instability. When the disk magnetisation increases above the threshold, these ringlike structures are smeared out and disappear as the accretion speed becomes supersonic. Our suspicion is that the ringformation 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 longterm evolution of the ringlike structures. Indeed, it is not clear if the rings are destroyed by nonaxissymmetric 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 magnetisation remains below a threshold value, around 10^{−3} − 10^{−2}, ringlike 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 protoplanetary disks (see for instance de Valon et al. 2020).
Similarly to Zhu & Stone (2018).
Acknowledgments
This work was granted access to the HPC resources of TGCC under the allocation A0080402231 attributed by GENCI (Grand Equipement National de Calcul Intensif). GL acknowledges support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme (Grant agreement No. 815559 (MHDiscs)). J.F. acknowledges support from the CNES french agency and CNRS PNHE. This research made use of the SciPy and NumPy libraries (Oliphant 2006; Van Der Walt et al. 2011; Virtanen et al. 2020). All the data reduction presented in this paper was performed using the GRICAD infrastructure (https://gricad.univgrenoblealpes.fr), which is supported by Grenoble research communities.
References
 Bai, X.N., & Stone, J. M. 2013, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2014, ApJ, 796, 31 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Bally, J., Reipurth, B., & Davis, C. J. 2007, Protostars and Planets V, 215 [Google Scholar]
 Begelman, M. C., McKee, C. F., & Shields, G. A. 1983, ApJ, 271, 70 [CrossRef] [Google Scholar]
 Begelman, M. C., Armitage, P. J., & Reynolds, C. S. 2015, ApJ, 809, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Cabrit, S., Edwards, S., Strom, S. E., & Strom, K. M. 1990, ApJ, 354, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Casse, F., & Ferreira, J. 2000, A&A, 361, 1178 [NASA ADS] [Google Scholar]
 Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Combet, C., & Ferreira, J. 2008, A&A, 479, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Coppejans, D. L., Körding, E. G., MillerJones, J. C. A., et al. 2015, MNRAS, 451, 3801 [NASA ADS] [CrossRef] [Google Scholar]
 Corbel, S., Fender, R. P., Tzioumis, A. K., et al. 2000, A&A, 359, 251 [NASA ADS] [Google Scholar]
 de Valon, A., Dougados, C., Cabrit, S., et al. 2020, A&A, 634, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Done, C., Davis, S. W., Jin, C., Blaes, O., & Ward, M. 2012, MNRAS, 420, 1848 [NASA ADS] [CrossRef] [Google Scholar]
 Dougados, C., Cabrit, S., Lavalley, C., & Ménard, F. 2000, A&A, 357, L61 [NASA ADS] [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Ferreira, J. 1997, A&A, 319, 340 [Google Scholar]
 Ferreira, J. 2002, EAS Publications Series, eds. J. Bouvier, & J.P. Zahn, 3, 229 [CrossRef] [EDP Sciences] [Google Scholar]
 Ferreira, J., & Pelletier, G. 1995, A&A, 295, 807 [NASA ADS] [Google Scholar]
 Ferreira, J., Dougados, C., & Cabrit, S. 2006a, A&A, 453, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferreira, J., Petrucci, P. O., Henri, G., Saugé, L., & Pelletier, G. 2006b, A&A, 447, 813 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 [NASA ADS] [CrossRef] [Google Scholar]
 Frieman, E., & Rotenberg, M. 1960, Rev. Mod. Phys., 32, 898 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallo, E., Fender, R. P., & Pooley, G. G. 2003, MNRAS, 344, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Gallo, E., Fender, R., & Kaiser, C. 2005, Interacting Binaries: Accretion Evolution, and Outcomes, 797, 189 [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [CrossRef] [Google Scholar]
 Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097 [NASA ADS] [CrossRef] [Google Scholar]
 Guilet, J., & Ogilvie, G. I. 2013, MNRAS, 430, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Guilet, J., & Ogilvie, G. I. 2014, MNRAS, 441, 852 [NASA ADS] [CrossRef] [Google Scholar]
 Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736 [Google Scholar]
 Hawley, J. F. 2001, ApJ, 554, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Heinemann, T., & Papaloizou, J. C. B. 2009a, MNRAS, 397, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Heinemann, T., & Papaloizou, J. C. B. 2009b, MNRAS, 397, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Higginbottom, N., & Proga, D. 2015, ApJ, 807, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Hirth, G. A., Mundt, R., & Solf, J. 1997, A&AS, 126, 437 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 JacqueminIde, J., Ferreira, J., & Lesur, G. 2019, MNRAS, 490, 3112 [CrossRef] [Google Scholar]
 JiménezIbarra, F., MuñozDarias, T., Casares, J., Armas Padilla, M., & CorralSantana, J. M. 2019, MNRAS, 489, 3420 [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, W.T., & Ostriker, E. C. 2000, ApJ, 540, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Konigl, A., & Pudritz, R. E. 2000, Protostars and Planets IV, 759 [Google Scholar]
 Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, Z.Y. 1995, ApJ, 444, 848 [NASA ADS] [CrossRef] [Google Scholar]
 Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656 [CrossRef] [Google Scholar]
 Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [Google Scholar]
 Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235 [NASA ADS] [Google Scholar]
 Marcel, G., Ferreira, J., Petrucci, P.O., et al. 2018, A&A, 617, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marcel, G., Ferreira, J., Clavel, M., et al. 2019, A&A, 626, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Markoff, S., Nowak, M., Corbel, S., Fender, R., & Falcke, H. 2003, A&A, 397, 645 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Bland ford, R. D. 2012, MNRAS, 423, 3083 [NASA ADS] [CrossRef] [Google Scholar]
 Merloni, A., Heinz, S., & di Matteo, T. 2003, MNRAS, 345, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L. 1961, MNRAS, 122, 473 [NASA ADS] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mirabel, I. F., & Rodríguez, L. F. 1999, ARA&A, 37, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Mishra, B., Begelman, M. C., Armitage, P. J., & Simon, J. B. 2020, MNRAS, 492, 1855 [NASA ADS] [CrossRef] [Google Scholar]
 Murphy, G. C., Ferreira, J., & Zanni, C. 2010, A&A, 512, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [NASA ADS] [CrossRef] [Google Scholar]
 Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (USA: Trelgol Publishing) [Google Scholar]
 Pessah, M. E., Chan, C.K., & Psaltis, D. 2007, ApJ, 668, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Petrucci, P.O., Ursini, F., De Rosa, A., et al. 2018, A&A, 611, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 277 [Google Scholar]
 Ray, T. P., Mundt, R., Dyson, J. E., Falle, S. A. E. G., & Raga, A. C. 1996, ApJ, 468, L103 [NASA ADS] [CrossRef] [Google Scholar]
 Riols, A., & Lesur, G. 2019, A&A, 625, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rothstein, D. M., & Lovelace, R. V. E. 2008, ApJ, 677, 1221 [NASA ADS] [CrossRef] [Google Scholar]
 Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857 [NASA ADS] [CrossRef] [Google Scholar]
 Scepi, N., Dubus, G., & Lesur, G. 2019, A&A, 626, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scepi, N., Lesur, G., Dubus, G., & JacqueminIde, J. 2020, A&A, 641, A133 [EDP Sciences] [Google Scholar]
 Serjeant, S., Rawlings, S., Lacy, M., et al. 1998, MNRAS, 294, 494 [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Sorathia, K. A., Reynolds, C. S., & Armitage, P. J. 2010, ApJ, 712, 1241 [NASA ADS] [CrossRef] [Google Scholar]
 Stepanovs, D., & Fendt, C. 2016, ApJ, 825, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Suzuki, T. K., & Inutsuka, S.I. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Suzuki, T. K., & Inutsuka, S.I. 2014, ApJ, 784, 121 [Google Scholar]
 Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tetarenko, B. E., Lasota, J.P., Heinke, C. O., Dubus, G., & Sivakoff, G. R. 2018, Nature, 554, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Wang, L., & Goodman, J. J. 2017, ApJ, 835, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Wardle, M., & Königl, A. 1993, ApJ, 410, 218 [NASA ADS] [CrossRef] [Google Scholar]
 White, C. J., Stone, J. M., & Quataert, E. 2019, ApJ, 874, 168 [CrossRef] [Google Scholar]
 Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Jiang, Y.F., & Stone, J. M. 2020, MNRAS, 495, 3494 [Google Scholar]
Appendix A: Reynolds averaged MHD equations
In this appendix we compute the Reynolds averages of the MHD equations. We proceed by decomposing the different terms into fluctuating (turbulent) terms plus the laminar (averaged) terms, using the averages defined in Sect. 2.3.1. Then we can average Eqs. (1)–(4)) with respect to time and φ. We start with the equation for the conservation of mass, Eq. (1), thanks to its simple form this is straightforward
It is important to note that the toroidal component of the divergence vanishes thanks to Eq. (12). Indeed, the averaged quantities are only function of the latitudinal, θ, and radial, r, coordinates. The term on the left hand side of the equation is associated with the nonsteadiness of the turbulence, if the simulation is quasi stationary this term should be negligible. Then we average the toroidal component of the induction equation, Eq. (4):
where the turbulent emf is defined as
We write the poloidal component of the induction equation as the toroidal component of the equation for the vector potential
where A is the vector potential of the magnetic field, ∇×A = B and G emerges to satisfy the gauge condition. Averaging this equation leads to
where the term ∇G is averaged out thanks to the condition of Eq. (12). We project, then we rearrange the equation of latitudinal momentum using u_{θ} = ⟨u_{θ}⟩ + δu_{θ} and the conservation of mass to find a not so usual expression:
where cotθ = cos θ/sin θ. We decompose the rest of the quantities into turbulent and laminar terms and average with respect to time and the toroidal coordinate to find:
where ⟨u⟩_{ρ} = ⟨u⟩ + ⟨δu⟩_{ρ} and
In the derivation of Eq. (A.7) we have used
In the same way as before the derivative with respect to the toroidal component vanishes thanks to Eq. (12). We repeat the same procedure for the radial projection of Eq. (2) to find
as before we then average the equation and find:
We project and rearrange the equation for toroidal momentum transport to find the usual equation of angular momentum transport:
we use the conservation of mass and then repeat the averaging procedure to find:
Appendix B: Computation of fluxes
In this section we show how we compute the different fluxes for the equation of conservation of angular momentum and the equation of conservation of mass. First we define the following useful average
We can then integrate Eq. (1) with respect to time and space to find
where we have introduced
we can redo the same exercice for Eq. (A.12) to find
There is an additional term F_{Mθ} when compared to the cylindrical case due to the spherical geometry. The angular momentum fluxes are defined with respect to an equilibrium velocity, Rω. We define
the deviations to the equilibrium velocity, we are free to choose the function ω(r, θ). Indeed, a proper choice of the equilibrium velocity, Rω, leads to a simplification of this term. we define the angular momentum fluxes as
We can simplify the term F_{Mθ} by choosing ω = Ω_{K}/sin^{2}θ. This leads to a mathematical form similar to the one obtained in the cylindrical coordinates. It is important to note that the function ω(r, θ) is different from the mean angular velocity of the flow ⟨Ω⟩,
B.1. Mass flux
To get a global understanding of mass accretion and ejection, let us define the fluxes through the different disk surfaces. Using the regions defined above (Sect. 3.1), we compute the fluxes through the boundaries of each regions, delimited radially by r_{1} = 8 and r_{2} = 14 (see Fig. 3, top). We then integrate Eq. (1) and get
where the mass fluxes are timeaveraged on a duration Δt which spans several local orbital periods. A strict steady state would translate into ΔM = 0. To illustrate timevariability, we first average the mass budget on a sliding window of length Δt = 200/Ω_{K}(R_{in}). We show in Fig. B.1 the fluxes calculated within the disk region. First, we find that there is a coherent latitudinal mass flux removing mass from the disk and transferring it to the atmosphere. In addition, the radial mass flux seems quasiperiodic, with a period of the order of the local orbital period. This is most probably due to the stochastic excitation of spiral density waves by MRIdriven turbulence (Heinemann & Papaloizou 2009a,b).
Fig. B.1. Mass fluxes inside the disk, defined for z< 3h as function of the local orbital time at r_{1} = 8. The different fluxes are defined in Appendix B, negative fluxes corresponding to mass being removed from the disk. 
The large variability of the radial fluxes on orbital timescales makes their secular interpretation obtuse. To characterise the mass budget inside the disk, we therefore average the mass fluxes over the whole time domain from t_{1} = 923T_{in} to t_{2} = 1114T_{in} and plot them as simple scalars (Fig. B.2). We recover the mass transfer from the disk to the atmosphere, but we now also have 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 ringlike structures within the disk region, a form of secular selforganisation. 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 ringlike structures is discussed in Sect. 4.
Fig. B.2. Mass fluxes inside of the whole domain, the regions shaded blue are turbulent regions. The different fluxes are defined in Appendix B. The arrows represent the direction of the fluxes. Since we are representing the different mass fluxes of a region of the disk, their sum (taking into account their direction) should be equal to the mass difference between t_{1} and t_{2}. 
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. 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.
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 righthand side measure the angular momentum flux due to the torques (turbulent and laminar). The angular momentum fluxes are quasistationary 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 summarise their values in Fig. B.3.
Fig. B.3. Angular momentum fluxes accross the whole domain, the regions shaded blue are turbulent regions. 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, r2}, F_{L, r1}, 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. 
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 momentum, 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 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 momentum extracted from this region, so most of it is actually transported radially outwards or transferred down to the laminar atmosphere. In contrast to the laminar atmosphere however, the transport is here due to turbulent torques (Fig. 3). This is be consistent with an MRIturbulent region, as shown in Sect. 3.5.
Appendix C: MRI in the low β regime
We rederive here in a simplified manner some of the results shown by Kim & Ostriker (2000) for the MRI in the strong field (or cold plasma) limit. Consider an ideal, isothermal plasma in differential rotation around a central object and set up a cylindrical coordinates system (R, φ, z). Assume that the plasma is initially at equilibrium, with angular velocity Ω(R) constant on cylinders, and constant density ρ_{0}, pressure (where c_{s} is the isothermal sound speed) and magnetic field B_{0}. Take a fiducial radius R_{0} and attach a cartesian frame (x, y, z) centred on R_{0} and rotating at Ω = Ω(R_{0}) so that x corresponds to the local radial direction R − R_{0} and y the local azimuthal direction φ. Let us then add a small perturbation to this plasma so that ρ = ρ_{0} + δρ and B = B_{0} + δb. In order to characterise the evolution of this perturbation, we use the displacement vector field ξ, for which the equation of motion in the rotating frame (x, y, z) reads
where we have introduced the rotation vector Ω = Ωe_{z}. One recognises in the last two terms the Coriolis force and the radial effective potential. The density and magnetic perturbations are easily deduced from (1) and (4)
which can be combined to the equation of motion on ξ to get
In this equation, we have introduced the Alfvén velocity vector and q ≡ −d logΩ/d log R. This equation is reminiscent of the traditional stability analysis of ideal plasmas (e.g., Frieman & Rotenberg 1960), with the addition of inertial and gravitational forces.
Although this equation is in principle general, we simplify it by considering a background field with only toroidal and vertical components B_{0} = B_{0, φ}e_{y} + B_{0, z}e_{z}, and perturbations with only a vertical spatial dependency. We then seek for growing (unstable) solutions of the form . One eventually gets the dispersion relation
where we have introduced the epicyclic frequency κ^{2} ≡ 2Ω^{2}(2 − q). Although general solutions can be derived for this dispersion relation, they are not really enlightening. It is more instructive to look at the weakly magnetised and strongly magnetised limits.
In the weakly magnetised case (V_{A}/c_{s} → 0) the first line of (C.2) dominates, and we recover^{4} the usual dispersion relation for the MRI in the incompressible limit (Balbus & Hawley 1991). The stability condition is then given by the constant term: . The usual thinking is that this stability condition implies that β < 1 disks are stable for the MRI. This statement is correct provided that 1/k < H ≃ c_{s}/Ω, that is the vertical perturbation wavelength fits in the hydrostatic disk scale height. However, in the disk atmosphere, the vertical scale height is not set by the hydrostatic equilibrium, and is actually much larger than this, so longer wavelength perturbations can in principle exist and be unstable in the atmosphere.
If we now take the strong field limit V_{A}/c_{s} → ∞, another picture starts to emerge. The first line of the dispersion relation then becomes negligible, and the stability condition now reads . Obviously, a weak enough V_{Aφ} is required for the MRI to survive independently of k, in sharp contrast to the weak field case. More specifically, a necessary condition for instability is . Let us finally point out that in this limit, Eq. (C.2) reduces to Eq. (46) in Kim & Ostriker (2000), which corresponds to their cold plasma limit.
Overall, we find that in the strongly magnetised regime, where V_{A} ≫ c_{s}, the condition of existence for the MRI is affected by the strength of the toroidal field, a strong enough field stabilising the plasma. This is confirmed by computing the growth rate γ in the general case from the 6th order polynomial (C.2). We illustrate this in Fig. 10, which shows the growth rate (maximised on k) as a function of the field strength V_{Aφ}/V_{Az} for the Keplerian case q = 3/2. We believe this is the main explanation to why turbulence disappears from the laminar atmosphere and then reappears once V_{Aφ}/V_{Az} is sufficiently low.
Appendix D: Transport of magnetic field and magnetic flux
We show in Fig. D.1 for various simulations the spacetime diagrams of the magnetic flux Ψ (defined in Sect. 3.2.1) and, for practical reasons, the ratio P_{b, mid}/P_{mid}, where
Fig. D.1. Magnetic field flux Ψ defined in Eq. (23) and ratio P_{b, mid}/P_{mid} (see text) as functions of time and the radial coordinate for the different simulations: SB4 (top left), SB3 (top right), SB2 (bottom left), and SEp (bottom right). We note that P_{b, mid}/P_{mid} > ⟨β_{mid}⟩^{−1} since the former contains also the turbulent magnetic pressure. 
and θ_{d2, d1} = π/2 ± arctan(h_{disk}/R). It is important however to note that the ratio P_{b, mid}/P_{mid} is larger than since it contains the turbulent field contribution ⟨δB^{2}⟩.
While the turbulent magnetic pressure measured at the mid plane in weakly magnetised disks remains always subthermal (see Fig. 7), it becomes suprathermal at large magnetisations. For instance, in the case of a disk with ⟨β_{mid}⟩^{−1} ≃ 1, we measure ⟨δB^{2}⟩(r, θ = π/2)/8π ∼ 10P_{0}. Thus, the yellow disk region 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 magnetisation ⟨β_{mid}⟩^{−1} is around unity.
All Tables
Parameters, grid spacing and radial extension of the simulations discussed in the paper.
Values of the MHD invariants for all simulations, measured in the upper ‘hemisphere’ for a field line originating at R_{0} = 6.
All Figures
Fig. 1. Computational grid. 

In the text 
Fig. 2. Top, left: gas density and mean poloidal stream lines. The red dotted line corresponds to the Alfvénic surface, and the red dashed line corresponds to the fast magnetosonic surface. The colour of the poloidal stream lines correspond to the logarithm of their magnitude normalised to the sound speed. Top, right: RB_{φ} normalised to B_{i}R_{in}; and mean poloidal field lines. The grey square corresponds to the zoomed in region the bottom figure. Bottom: same as top but zoomed in the greyed region. The black dashed line indicates the surface where ⟨β⟩ = 8π⟨P⟩/⟨B⟩^{2} = 1. 

In the text 
Fig. 3. Top: ratio α_{la}/α_{tu} denoting turbulent (red) and laminar (blue) regions. The dotted lines delimit these regions and represent z = 3.5h, z = 6.5h and z = 12.5h. The two different circles denote the radii used for the calculation of the fluxes (r_{1} = 8, r_{2} = 14), see Sects. 3.2.1 and Appendix B. TA and LA correspond for turbulent and laminar atmosphere, respectively. Bottom: cosine of the angle ψ between the mean poloidal velocity and the mean poloidal magnetic field as a function of the latitudinal coordinate. The grey zones correspond to the turbulent regions in the top panel, delimited by the dotted lines. 

In the text 
Fig. 4. Top: velocity profiles as functions of the latitudinal coordinate, θ, normalised to the local sound speed. Bottom: mean magnetic field profiles normalised to the vertical magnetic field in the disk midplane, ⟨B_{z}⟩(r, θ = π/2) as functions of the latitudinal coordinate. The grey zones correspond to the turbulent regions in Fig. 3. All profiles are radially averaged between r_{1} = 8 and r_{2} = 14 after being normalised. 

In the text 
Fig. 5. Magnetic flux threading the disk midplane Ψ, computed between the axis and R, as a function of time in innermost orbital units. Every contour line represents a field line. The dashed line shows an example of a fit using Eq. (26) with R_{i} = 4.4, t_{i} = 328T_{in} and . 

In the text 
Fig. 6. 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 normalised 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. 

In the text 
Fig. 7. Comparison of the different pressure terms normalised to the thermal pressure within the disks midplane as functions of the latitudinal coordinate. The shaded areas denote the turbulent regions. 

In the text 
Fig. 8. Top: mean poloidal magnetic field lines in the (R,z) plane. The blue lines define the different critical surfaces (see Sect. 3.1) and they are represented as vertical lines in the lower panel. The grey zones determine the turbulent zones defined in Sect. 3.1. The lower zone corresponds to the disk and the upper zone corresponds to the turbulent atmosphere. Bottom: MHD invariants calculated along the field lines of the upper panel as functions of the latitudinal coordinate. The line style of the MHD invariants has a one to one correspondence with the field line where the invariant was calculated. The grey zone corresponds to the end of the turbulent atmosphere. 

In the text 
Fig. 9. Bernoulli invariant and its components for the field line originating at R_{0} = 6 as functions of the latitudinal coordinate θ. All quantities are normalised to . The shaded grey area corresponds to the turbulent atmosphere, the end of the turbulent atmosphere coincides with the SM surfaces. The vertical lines correspond to the different critical surfaces. The vertical dashed line is the fast magnetosonic surface while the vertical dotted line is the Alfvénic surface, consistent with the notation Fig. 8. 

In the text 
Fig. 10. 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}). 

In the text 
Fig. 11. Space time diagram of the column density Σ (up) and the vertically averaged vertical magnetic field b_{z} (down) as functions of the time in units of orbits at the inner most radius and of the cylindrical radial coordinate. The solid black lines represent the gaps located at R ≃ [6.6, 12.3] while the dashed black lines represent the rings located at R ≃ [4, 9.5, 15]. 

In the text 
Fig. 12. Radial (blue) and vertical (red) contributions to the mass evolution rate in local Ω_{K} units as a function of the cylindrical radii, R. A positive (negative) term leads to the formation of local gap (ring). The solid black lines represent the gaps located at R ≃ [6.6, 12.3] while the dashed black lines represent the rings located at R ≃ [9.5, 15]. We use symmetrical logscale that becomes linear within the interval [5 × 10^{−5}, −5 × 10^{−5}]. 

In the text 
Fig. 13. Vertically (within the disk region) averaged turbulent angular momentum transport coefficient, as a function of the mean poloidal plasma, ⟨β_{mid}⟩, evaluated in the disk midplane for different cylindrical radii, R = [4, 16]. We also show different scaling laws for to compare them with our simulation. 

In the text 
Fig. 14. Extent of the turbulent zones (disk and atmosphere) found in the different simulations, see Table 1. As the magnetic field increases, the turbulent layer goes down and eventually merges with the turbulent disk. 

In the text 
Fig. 15. Mean density of the different simulations normalised to its value at the disk midplane, see Table 1. The symbols correspond to the height at which the flow reaches the SM speed. The mean profile are also averaged radially between r = 8 and r = 12. 

In the text 
Fig. 16. Left: field advection velocity as a function of radius normalised to the Keplerian velocity. Right: ratio of field and mass advection velocities as a function of radius. We truncate the advection velocity where we can no longer measure the drift of the magnetic field lines. The truncation radii corresponds to R = 3.5, R = 10, R = 20, for SEp, SB4 and SB3 respectively. 

In the text 
Fig. 17. Inverse of potential vorticity and surface density as functions of the radial coordinate, for simulations SB4 and S2pi, w_{θ} = ∇×⟨u⟩_{θ}.. 

In the text 
Fig. 18. Appearance of the SB2 strong field simulation: logarithm of the density (background colour), mean poloidal magnetic field lines (black solid lines), Alfvén (red dotted line) and fast magnetosonic (red dashed line) critical surfaces. The black dotted lines represent the end of the disk region defined in Fig. 14. The quantities are averaged between t_{1} = 1719T_{in} and t_{2} = 1910T_{in} with a resolution of Δt = 0.16T_{in}. 

In the text 
Fig. B.1. Mass fluxes inside the disk, defined for z< 3h as function of the local orbital time at r_{1} = 8. The different fluxes are defined in Appendix B, negative fluxes corresponding to mass being removed from the disk. 

In the text 
Fig. B.2. Mass fluxes inside of the whole domain, the regions shaded blue are turbulent regions. The different fluxes are defined in Appendix B. The arrows represent the direction of the fluxes. Since we are representing the different mass fluxes of a region of the disk, their sum (taking into account their direction) should be equal to the mass difference between t_{1} and t_{2}. 

In the text 
Fig. B.3. Angular momentum fluxes accross the whole domain, the regions shaded blue are turbulent regions. 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, r2}, F_{L, r1}, 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. 

In the text 
Fig. D.1. Magnetic field flux Ψ defined in Eq. (23) and ratio P_{b, mid}/P_{mid} (see text) as functions of time and the radial coordinate for the different simulations: SB4 (top left), SB3 (top right), SB2 (bottom left), and SEp (bottom right). We note that P_{b, mid}/P_{mid} > ⟨β_{mid}⟩^{−1} since the former contains also the turbulent magnetic pressure. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.