© J. Jacquemin-Ide et al. 2021

1. Introduction

The emission of jet-like outflows from accretion disks is ubiquitous in the Universe, as both jets and disks are observed around a variety of astrophysical objects and on a wide range of spatial scales. Jet-like outflows are observed being emitted from active galactic nuclei (AGN) and quasar (Merloni et al. 2003 and references therein; Blandford et al. 2019), X-ray binaries (Mirabel & Rodríguez 1999; Corbel et al. 2000; Gallo et al. 2003, 2005), and possibly also nova-like cataclysmic variables (Coppejans et al. 2015) and protoplanetary disks (Burrows et al. 1996; Hirth et al. 1997; Ray et al. 1996; Dougados et al. 2000; Bally et al. 2007). All of these jet-like outflows exhibit several common properties: They are highly collimated, supersonic and highly correlated with the measured accretion rate in their emitting disk (Cabrit et al. 1990; Hartigan et al. 1995; Serjeant et al. 1998; Markoff et al. 2003; 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), X-ray binaries (Tetarenko et al. 2018; Jiménez-Ibarra et al. 2019), and proto-planetary 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 large-scale vertical magnetic field, both mechanisms naturally emerge. First, the seminal work of Balbus & Hawley (1991) showed that the magneto-rotationnal instability (MRI) is triggered, leading to vigorous 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 self-similar models (Ferreira & Pelletier 1995; Li 1995; Ferreira 1997) or simplified 2.5D simulations (Murphy et al. 2010; Stepanovs & Fendt 2016). These effective models generically rely on azimuthally-averaged equations of motions, in which turbulence is modelled by an effective viscosity and resistivity. While they have been successful in elucidating the constraints needed for steady ejection (Ferreira & Pelletier 1995) or the relation between disk structure and jet properties (Ferreira 1997; Stepanovs & Fendt 2016; Jacquemin-Ide et al. 2019); their use of ad hoc turbulent profile limits their predictive power as they are dependent on the prescription used to model turbulence (see Jacquemin-Ide et al. 2019 for a discussion).

MRI-driven turbulence has been studied using both shearing box simulations (Hawley et al. 1995) and more recently global simulations (Sorathia et al. 2010). Most of these simulations focused on the configuration where there is no mean field threading the disk (i.e. no external source of magnetic field), thereby excluding the possibility of a magnetically-driven outflow. In these ‘zero net flux simulations’, it is found that the MRI drives 3D MHD turbulence, which results in an efficient radial angular momentum transport and therefore accretion. It is only recently that, thanks to robust numerical methods and more abundant numerical resources, stratified shearing box models with a mean field have been developed. These models show both turbulent radial torques and magnetically-driven outflows (Suzuki & Inutsuka 2009; Fromang et al. 2013; Bai & Stone 2013). These models however suffer from several drawbacks due to the local nature of the shearing box model: the mass loss rate depends on the vertical domain size and the flow top/down symmetry is unphysical (Fromang et al. 2013). This implies that global curvature effects, absent from the shearing box model, are mandatory to obtain physically valid outflow solutions.

To this end, global numerical simulations combining MRI-driven turbulence and magnetically-driven outflows have been computed (Suzuki & Inutsuka 2014; Zhu & Stone 2018; Mishra et al. 2020). These simulations are known to exhibit an exotic global configuration. For instance, the stationary state found by Zhu & Stone (2018) possess 3 distinct features: a turbulent disk, a powerfully accreting atmosphere, mostly driven by a laminar Maxwell torque, and a tenuous wind that does not transport a considerable amount of angular momentum and mass.

These numerical results are in flagrant contradiction with usual two dimensional ‘effective’ models of wind-emitting disk. In these models, accretion typically occurs in the disk bulk and angular momentum escapes through a 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; Jacquemin-Ide et al. 2019), where h is the disk pressure scale-height, none of these models produce accretion in the disk atmosphere at z ∼ 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 X-ray 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 counter-balanced by field diffusion due to MRI turbulence. Rothstein & Lovelace (2008) later showed that the vertical structure plays a quintessential role on the secular transport of the magnetic field. Indeed, the upper layers of the disk, that are magnetically dominated and have a lower conductivity, could enhance the transport of magnetic field inwards. They showed that in the case of wind driven accretion it is possible to accumulate magnetic flux in the inner regions. Guilet & Ogilvie (2014) later showed that even in a purely viscous case, the vertical structure can help the transport of magnetic field towards the inner regions. However, all of these results rely on simplified models of turbulent transport, which are not entirely self-consistent. It is now mandatory to confirm the possibility of field transport from direct numerical simulations as it would then guide future theoretical and observational developments on this problem.

In this work we aim at understanding the dynamics of a disk subject to MRI-driven turbulence and 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 ring-like 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

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·\left[\rho \mathit{u}\right]=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \rho (\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}·\mathbf{\nabla })\mathit{u})=-\mathbf{\nabla }P+\rho \mathbf{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}+\frac{1}{c}\mathit{J}\times \mathit{B},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial P}{\partial t}+\mathit{u}·\mathbf{\nabla }P+\mathrm{\Gamma }P\mathbf{\nabla }·\mathit{u}=-\mathrm{\Lambda }\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}=\mathbf{\nabla }\times [\mathit{u}\times \mathit{B}],\hfill \end{array}$$(4)

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

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{P}{T}\frac{T-T_{\mathrm{eff}}}{{\tau }_{\mathrm{cool}}}\end{array}$$

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

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{P-\rho c_{\mathrm{s}}^2(R)}{{\tau }_{\mathrm{cool}}}\end{array}$$

which leads to an effective equation of state of the form

$$\begin{array}{c}\hfill P=c_{\mathrm{s}}^2(R)\rho ,\end{array}$$(5)

where c_s is the local isothermal sound speed

$$\begin{array}{c}\hfill c_{\mathrm{s}}=h{\mathrm{\Omega }}_{\mathrm{K}}\end{array}$$(6)

defined from the Keplerian angular velocity ${\mathrm{\Omega }}_{\mathrm{K}}=\sqrt{GM/R^3}$ and the disk geometrical scale height h. Finally, we define the poloidal velocity and magnetic field as

$$\begin{array}{cc}\hfill {\mathit{u}}_{\mathbf{p}}& =u_r{\mathit{e}}_{\mathit{r}}+u_{\theta }{\mathit{e}}_{\mathit{\theta }}\hfill \\ \hfill {\mathit{B}}_{\mathbf{p}}& =B_r{\mathit{e}}_{\mathit{r}}+B_{\theta }{\mathit{e}}_{\mathit{\theta }}.\hfill \end{array}$$

2.2. Numerical method

To solve Eqs. (1)–(5) on a spherical grid, we use the conservative Godunov-type code PLUTO (Mignone et al. 2007). Constrained transport (Evans & Hawley 1988) ensures the solenoidal condition for the magnetic field is satisfied at machine precision. We use a second order Runge-Kutta method for the time step combined with a HLLD type solver with linear reconstruction to handle Riemann problems.

The extent of our radial and toroidal domain varies from one simulation to the next, while our latitudinal extent is always [0, π]. We summarise the domain properties, as well as the main control parameters for each simulations in 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.

We use periodic boundary conditions in the toroidal direction and zero-gradient in the radial direction (mass flow into the numerical domain is forbidden). As for the latitudinal direction we implement the same polar (z-axis) boundary condition as Zhu & Stone (2018).

We follow Nelson et al. (2013) for the hydrostatic initial condition. The initial density and temperature profiles are defined at the disk mid-plane as:

$$\begin{array}{cc}\hfill {\rho }_i(R,z=0)& ={\rho }_i(R=R_{\mathrm{in}},z=0){\left(\frac{R_{\min }}{R_{\mathrm{in}}}\right)}^p\hfill \\ \hfill T_i(R)& =T_i(R=R_{\mathrm{in}}){\left(\frac{R_{\min }}{R_{\mathrm{in}}}\right)}^q\hfill \end{array}$$

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 mid-plane defined above. We get

$$\begin{array}{c}\hfill {\rho }_i(R,z)={\rho }_i(R,z=0)exp\left[\frac{{\mathrm{\Omega }}_{\mathrm{K}}^2(R_{\min })R_{\min }^3}{c_{\mathrm{s}}^2(R_{\min })}(\frac{1}{\sqrt{R_{\min }^2+z^2}}-\frac{1}{R_{\min }})\right]\end{array}$$

and

$$\begin{array}{c}\hfill u_{\phi ,i}(R,z)=V_{\mathrm{K}}{[(q+1)-q\frac{R_{\min }}{\sqrt{R_{\min }^2+z^2}}+\frac{c_{\mathrm{s}}^2(R_{\min })}{V_{\mathrm{K}}^2(R_{\min })}(p+q)]}^{1/2}\end{array}$$

where we have defined the Keplerian velocity V_K = Ω_KR. We then use the Keplerian angular velocity and the sound speed to define the disk scale height

$$\begin{array}{c}\hfill h=c_{\mathrm{s}}/{\mathrm{\Omega }}_{\mathrm{K}}=ϵR,\end{array}$$(7)

where ϵ is the disk geometrical thickness.

Following Zhu & Stone (2018), we assume a large scale vertical magnetic field that is constant with respect to z but follows a power law in R. In order to get a self-similar initial condition, we choose a radially constant initial plasma beta,

$$\begin{array}{c}\hfill {\beta }_{\mathrm{ini}}=\frac{8\pi P_i}{B_{z,i}^2},\end{array}$$(8)

in the disk midplane. The magnetic field radial profile is then related to that of temperature and density:

$$\begin{array}{c}\hfill B_{z,i}(R)=B_i{\left(\frac{R_{\min }}{R_{\mathrm{in}}}\right)}^{\frac{p+q}{2}}\end{array}$$(9)

where

$$\begin{array}{c}\hfill B_i=\sqrt{8\pi P_i(R=R_{\mathrm{in}},z=0)/{\beta }_{\mathrm{ini}}}.\end{array}$$(10)

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

$$\begin{array}{c}\hfill A_{\phi }=\{\begin{array}{c}\frac{R{B}_i}{2}\mathrm{for}R\le R_{\mathrm{in}},\hfill \\ B_i{\left(\frac{R_{\min }}{R_{\mathrm{in}}}\right)}^{\frac{p+q}{2}}\frac{2R}{(p+q)+4}+B_i\frac{R_{\mathrm{in}}^2}{R}(\frac{1}{2}-\frac{2}{(p+q)+4})\mathrm{for}R>R_{\mathrm{in}}.\hfill \end{array}\end{array}$$

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

$$\begin{array}{c}\hfill {\rho }_{\mathrm{fl}}=\max ({\rho }_{\mathrm{fl},0}{\left(\frac{R_{\min }}{R_{\mathrm{in}}}\right)}^p\frac{1}{z^2+a_{\mathrm{fl}}{ϵ}^2{R}_{\min }^2},{10}^{-9})\end{array}$$

where ρ_fl, 0 = 10⁻⁷ and a_fl ∼ 1. Finally, to minimise integration time, we enforce a maximal Alfvénic velocity V_A, max = V_K(R_in), where $V_{\mathrm{A}}=B/\sqrt{4\pi \rho }$ 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.

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:

$$\begin{array}{cc}\hfill {\langle X\rangle }_{\phi }& =\frac{1}{\mathrm{\Delta }\phi }{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}X\mathrm{d}\phi ,}\hfill \\ \hfill \langle X\rangle & =\frac{1}{\mathrm{\Delta }\phi }{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}\frac{1}{t_2-t_1}\underset{t_1}{\overset{t_2}{\int }}X\mathrm{d}\phi \mathrm{d}t,}\hfill \\ \hfill {\langle X\rangle }_{\rho }& =\frac{1}{\mathrm{\Delta }\phi \langle \rho \rangle }{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}\frac{1}{t_2-t_1}\underset{t_1}{\overset{t_2}{\int }}\rho X\mathrm{d}\phi \mathrm{d}t=\frac{1}{\langle \rho \rangle }\langle X\rho \rangle .}\hfill \end{array}$$

In these definitions, φ₁ and φ₂ 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

$$\begin{array}{c}\hfill X=\langle X\rangle +\delta X.\end{array}$$(11)

With this definition, the quantities ⟨X⟩(R, z) are equivalent to the ones used in 2d effective models such as self-similar solutions (see Jacquemin-Ide et al. 2019). Hence, they allow for a direct comparison to said models and provide guidance for turbulence closures. Let us emphasise however that ⟨δX⟩_ρ ≠ 0, which leads to additional complications when dealing with high order correlation functions. Finally, we note that

$$\begin{array}{c}\hfill \langle \frac{\partial X}{\partial \phi }\rangle =0\end{array}$$(12)

since the boundary condition in the φ direction is periodic. The averages are taken from time t₁ to t₂, 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 accretion-ejection system, we use the averages defined above to rewrite the ideal MHD equations as a set of Reynolds-averaged MHD equations. In this approach, turbulence then appears in the averaged equation as an additional dynamical term influencing the mean flow. This is actually an extension of the famous Shakura & Sunyaev (1973)α-disk approach which can be derived as a Reynolds-averaged model for MHD turbulence (Balbus & Papaloizou 1999). While the α disk approach only focuses on angular momentum transport, our objective is to generalise this to the whole set of ideal MHD equations. To illustrate our procedure, we show its application to the equation of angular momentum transport below. A complete derivation for the whole set of ideal MHD equations can be found in Appendix A.

We project and rearrange (using Eq. (1)) the equation for toroidal momentum transport to find the equation of angular momentum transport

$$\begin{array}{c}\hfill \frac{\partial R{u}_{\phi }\rho }{\partial t}+\mathbf{\nabla }·R[\rho u_{\phi }\mathit{u}-\frac{B_{\phi }\mathit{B}}{4\pi }]=-\frac{\partial }{\partial \phi }(\frac{B^2}{8\pi }+P).\end{array}$$(13)

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

$$\begin{array}{cc}& \langle \rho \rangle {\langle {\mathit{u}}_{\mathit{p}}\rangle }_{\rho }·\mathbf{\nabla }R\langle u_{\phi }\rangle +\mathbf{\nabla }·R[\langle \rho \rangle {\langle \delta u_{\phi }\mathit{\delta }{\mathit{u}}_{\mathit{p}}\rangle }_{\rho }\hfill \\ \hfill & -\frac{1}{4\pi }\langle \delta B_{\phi }\mathit{\delta }{\mathit{B}}_{\mathit{p}}\rangle -\frac{1}{4\pi }\langle B_{\phi }\rangle \langle {\mathit{B}}_{\mathit{p}}\rangle ]=0.\hfill \end{array}$$(14)

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:

$$\begin{array}{c}\hfill \mathcal{T}={\mathcal{T}}_{\mathbf{tu}}+{\mathcal{T}}_{\mathbf{la}},\end{array}$$(15)

where we define

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathbf{la}}=-\frac{1}{4\pi }\langle B_{\phi }\rangle \langle {\mathit{B}}_{\mathit{p}}\rangle ,\end{array}$$(16)

and

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathbf{tu}}=\langle \rho \rangle {\langle \delta u_{\phi }\mathit{\delta }{\mathit{u}}_{\mathit{p}}\rangle }_{\rho }-\frac{1}{4\pi }\langle \delta B_{\phi }\mathit{\delta }{\mathit{B}}_{\mathit{p}}\rangle .\end{array}$$(17)

We can then define the usual turbulent alpha viscosity as

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{ef}}& =\frac{{\mathcal{T}}_{r\phi }}{\langle P\rangle },\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{ef}}& ={\alpha }_{\mathrm{la}}+{\alpha }_{\mathrm{tu}}\hfill \end{array}$$(19)

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¹ 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 magneto-sonic surface of the fiducial β_ini = 10⁴ run. We define the Alfvén surface as the location where ⟨u_p⟩ = ⟨V_ap⟩, where $V_{\mathit{ap}}=\langle B_p\rangle /\sqrt{4\pi \langle \rho \rangle }$ and the fast magneto-sonic surface where ⟨u_p⟩ = ⟨V_fm⟩, with $V_{\mathrm{fm}}^2=\frac{1}{2}(c_{\mathrm{s}}^2+V_a^2+|V_a^2-c_{\mathrm{s}}^2|)$, having defined the total Alfvénic velocity $V_{\mathrm{a}}=\langle B\rangle /\sqrt{4\pi \langle \rho \rangle }$.

In the top panel of Fig. 2 we show the whole domain of our numerical simulation. We find that the upper and lower hemispheres of the system are approximately symmetric. The fast magneto-soni (FM) surface is conical up to the border of our domain. In contrast, the Alfvénic surface remains conical up to a certain radii where its shape is modified, demonstrating that the outer radii of the domain has not yet reached a stationary state. This is corroborated by the lack of toroidal field within the outer regions of the disk.

In Fig. 2 (top right) we can clearly see some collimation happening for the inner field lines, close to the end of their trajectory. This may be the sign of some intrinsic MHD property, as it has been also observed in self-similar jet solutions (Ferreira 1997; Jacquemin-Ide et al. 2019). However, it is unclear whether this is a physical or a numerical effect, akin to a bias of our boundary conditions or a consequence of the density floor. Nevertheless, a study of collimation in our numerical models is outside the scope of this paper.

Since the FM surface corresponds to the point where the outflow becomes causally disconnected from the disk, we can use this surface to estimate the last radius that has achieved stationarity. We do this by computing the anchoring radius of the last field line that traverses the FM surface, this corresponds to R ∼ 35. Hence, we 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 magneto-sonic 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 sub-regions, dominated respectively by a laminar and a turbulent stress, invites us to define a laminar and a turbulent atmosphere (LA and TA), which we need to delimit more precisely.

Under the assumptions of ideal MHD and steady state, the average poloidal velocity and magnetic field should be perfectly aligned (e.g., Mestel 1961). Hence, the value of the angle ψ between the averaged poloidal stream and field lines,

$$\begin{array}{c}\hfill cos\psi =\frac{\langle {\mathit{u}}_{\mathit{p}}\rangle ·\langle {\mathit{B}}_{\mathit{p}}\rangle }{\langle u_p\rangle \langle B_p\rangle },\end{array}$$

identifies the regions where the turbulent terms are important. In Fig. 3 (bottom) we plot this angle as a function of the latitudinal coordinate. We also provide in Fig. 4 the behaviour of the mean velocity and mean magnetic fields as functions of the latitudinal coordinate.

In the disk, the mean poloidal velocities and field lines are approximately perpendicular (cos(ψ) = 0): the disk is turbulent. We identify the disk surface as the location where cos ψ = ±1. In the following, we define the disk as |z|< 3.5h (the central shaded area in 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 ideal-MHD is broken in this region, which we identify as the turbulent atmosphere, located at altitudes 6.5h < |z|≲12.5h. This region ends at the surface delimiting the base of the wind, where cos(ψ) = ± 1, which coincides with the slow magneto-sonic (SM) surface defined as ⟨u_p⟩ = V_SM, where $V_{\mathrm{SM}}^2=\frac{1}{2}(c_{\mathrm{s}}^2+V_a^2-|V_a^2-c_{\mathrm{s}}^2|)$ (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⟩² < 1. However, the mean poloidal magnetic field stays infra-thermal within the laminar atmosphere, ⟨β_p⟩ > 1, where

$$\begin{array}{c}\hfill \langle {\beta }_p\rangle =\frac{8\pi \langle P\rangle }{{\langle B_p\rangle }^2},\end{array}$$(20)

but becomes supra-thermal 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₁ = 8 and r₂ = 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₁ and r₂ decreases with time. This is the consequence of the secular self-organisation of the radial disk structure. Indeed, a closer look at the radial density profiles reveals the formation of ring-like structures. Our chosen disk section turns out to be located on a gap. The mechanism behind the formation of these structures is discussed in 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,

$$\begin{array}{c}\hfill v_{\mathrm{acc}}=\frac{\underset{{\theta }_{\mathrm{SM},1}}{\overset{{\theta }_{\mathrm{SM},2}}{\int }}rsin\theta \langle u_r\rho \rangle \mathrm{d}\theta .}{\underset{{\theta }_{\mathrm{SM},1}}{\overset{{\theta }_{\mathrm{SM},2}}{\int }}rsin\theta \langle \rho \rangle \mathrm{d}\theta .},\end{array}$$(21)

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

$$\begin{array}{c}\hfill v_{\mathrm{acc}}(R)\simeq -1.1\times {10}^{-3}{V}_{\mathrm{K}}(R),\end{array}$$(22)

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 self-organisation. Moreover, the accretion time scale we deduce from this: t_acc ∼ 9.1 × 10²/Ω_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

$$\begin{array}{c}\hfill \mathrm{\Psi }(R,t)={\displaystyle \underset{0}{\overset{\pi /2}{\int }}{r}^{2}sin\theta {\langle B_r\rangle }_{\phi }(r=R_{\text{in}},\theta ,t)\mathrm{d}\theta +\underset{R_{\mathrm{in}}}{\overset{R}{\int }}R{\langle B_z\rangle }_{\phi }(r,\theta =\pi /2,t)\mathrm{d}R}\end{array}$$(23)

which corresponds to an integration on a surface that goes from the axis (z = 1 and R = 0) down to the inner boundary of the disk (z = 0 and R = R_in) and then to radius R in the disk midplane.

In our fiducial simulation, we first confirm that Ψ(R = R_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.

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)

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Psi }}{\partial t}+v_{\mathrm{\Psi }}\frac{\partial \mathrm{\Psi }}{\partial R}=0.\end{array}$$(24)

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

$$\begin{array}{c}\hfill v_{\mathrm{\Psi }}\sim -1\times {10}^{-3}{V}_{\mathrm{K}},\end{array}$$(25)

valid for R ∈ [2, 11]. For R > 11 we cannot conclude on magnetic flux advection, since we did not integrate for long enough to get a measurable deviation from the initial condition. As expected from Fig. 5 the advection velocity is negative within the inner regions of the disk. Moreover, the advection velocity follows a radial dependency close to Keplerian. We cross-checked the numerical value of the magnetic field advection velocity by fitting the contours of Fig. 5 with

$$\begin{array}{c}\hfill R_{\mathrm{\Psi }}(t)=R_i{[1-\frac{3{\mathrm{\Omega }}_{\mathrm{K}}(R_i)}{2}\frac{v_{\mathrm{\Psi }}}{V_{\mathrm{K}}}(t-t_i)]}^{\frac{2}{3}}\end{array}$$(26)

where R_Ψ is the position of the contour, R_i is its initial position and t_i is the time the system needs to stabilise the transient due to the initial condition and start advecting magnetic flux, a delay visible in Fig. 5. Equation (26) is the integrated form of an equation describing a field line advected at a velocity proportional to the local Keplerian velocity, Ṙ_Ψ = v_Ψ(R). These fits show a reasonable agreement with the data (Fig. 5, dashed line), which suggests that magnetic field lines are advected in a self-similar manner once the local equilibrium is reached (t_i ∼ 30T_K(R)).

Since v_Ψ ∼ v_acc, the time scale for the magnetic field transport, namely t_Ψ = 10³/Ω_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)

$$\begin{array}{cc}& \frac{{\langle u_{\phi }\rangle }^2}{r}-\frac{{\langle B_{\phi }\rangle }^2+\langle \delta B_{\phi }^2\rangle }{4\pi r\langle \rho \rangle }-\frac{1}{\langle \rho \rangle }\frac{\partial \langle P\rangle }{\partial r}\hfill \\ \hfill & -\frac{1}{\langle \rho \rangle }\frac{\partial \langle P_B\rangle }{\partial r}+\frac{1}{4\pi \langle \rho \rangle }\langle {\mathit{B}}_{\mathit{p}}\rangle ·\mathbf{\nabla }\langle B_r\rangle =g,\hfill \end{array}$$(27)

$$\begin{array}{cc}& cot\theta [{\langle u_{\phi }\rangle }^2-\frac{{\langle B_{\phi }\rangle }^2+\langle \delta B_{\phi }^2\rangle }{4\pi \langle \rho \rangle }]-\frac{1}{\langle \rho \rangle }\frac{\partial \langle P\rangle }{\partial \theta }-\frac{1}{\langle \rho \rangle }\frac{\partial \langle P_B\rangle }{\partial \theta }=0,\hfill \end{array}$$(28)

where g = GM/r² and ⟨P_B⟩ = (⟨B⟩²+⟨δB²⟩)/8π is the magnetic pressure, which contains a laminar, ⟨B⟩², and a turbulent, ⟨δB²⟩, 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.

We find that the disk region is in radial and latitudinal equilibrium thanks to the balance between the centrifugal force, ⟨u_φ⟩²/r or cotθ⟨u_φ⟩², 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, $\langle B_{\phi }^2/(4\pi \rangle \rho r)$ or $\cot \theta \langle B_{\phi }^2/(4\pi \rangle \rho )$. 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 mid-plane (Fig. 4 bottom), that changes sign within the turbulent atmosphere.

The strong toroidal field within the laminar atmosphere is consistent with the fact that the laminar atmosphere is in ideal MHD. In this region, accretion drags the poloidal field lines 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 Jacquemin-Ide et al. 2019 for more details on the electric current vertical profile).

Overall, it is clear from Fig. 7 that the dynamical importance of the turbulent pressure cannot be understated. Indeed, Ferreira & Pelletier (1995) showed that the mass loaded onto the field lines depends crucially on the disk pressure which pushes material up to the ‘surface’, and loads the field lines. It is clear from Fig. 7 that the laminar magnetic pressure is compressing the disk surface, as shown by Ferreira & Pelletier (1995). However, the effect of this compression seems nonexistent in Fig. 6. This is a consequence of the presence of a important turbulent magnetic pressure in the disk. Indeed, the compression of the disk surface due to the laminar magnetic pressure, ⟨B⟩², is compensated by the turbulent magnetic pressure, ⟨δB²⟩ 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 magneto-sonic wind

In the regions where the flow is approximately laminar and steady-state, it is possible to define a set of MHD invariants which consists of conserved physical quantities along the field lines. These invariants characterise MHD outflows and allow direct comparisons to be made with other outflow solutions.

Under these conditions, Eq. (1) reduces to:

$$\begin{array}{c}\hfill \langle {\mathit{u}}_{\mathit{p}}\rangle =\frac{\stackrel{\sim }{\eta }}{4\pi \langle \rho \rangle }\langle {\mathit{B}}_{\mathit{p}}\rangle ,\end{array}$$(29)

where $\stackrel{\sim }{\eta }$ measures how much mass is being loaded onto the field lines. It can be transformed into a dimensionless mass loading parameter $\kappa \equiv \stackrel{\sim }{\eta }{R}_{\mathrm{SM}}{\mathrm{\Omega }}_{\mathrm{K}}(R_{\mathrm{SM}})/B_{z,\mathrm{SM}}$, where R_SM is the anchoring radius of the magnetic field line and the SM stands for quantities evaluated at the slow magneto-sonic surface.

Following the same regime of approximation, we can also write the φ component of the induction equation as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\star }=\langle \mathrm{\Omega }\rangle -\stackrel{\sim }{\eta }\frac{\langle B_{\phi }\rangle }{4\pi \langle \rho \rangle R},\end{array}$$(30)

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)

$$\begin{array}{c}\hfill L={\mathrm{\Omega }}_{\star }{R}_{\mathrm{A}}^2=\langle \mathrm{\Omega }\rangle R^2-\frac{R}{\stackrel{\sim }{\eta }}\langle B_{\phi }\rangle ,\end{array}$$(31)

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

$$\begin{array}{c}\hfill E=\frac{{\langle u\rangle }^2}{2}+{\mathrm{\Phi }}_G-{\mathrm{\Omega }}_{\star }\frac{R\langle B_{\phi }\rangle }{\stackrel{\sim }{\eta }}+w,\end{array}$$(32)

where w is the thermal energy term which can be written as

$$\begin{array}{c}\hfill w={\displaystyle \int }\frac{\mathbf{\nabla }\langle P\rangle }{\langle \rho \rangle }·\mathrm{d}\mathit{l},\end{array}$$(33)

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 mid-plane defining $e=\frac{E}{E_{\mathrm{K},\mathrm{SM}}}$, where $E_{\text{K},\text{SM}}={\mathrm{\Omega }}_{\text{K}}^2(R_{\text{SM}})R_{\text{SM}}^2/2$.

We follow three different field lines originating at R₀ = [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].

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 self-similar models (Jacquemin-Ide et al. 2019).

Secondly, the angular momentum invariant λ ≃ 4—6, so the wind is effectively free (λ > 3/2). This provides $R_{\text{A}}^2/R_{\text{SM}}^2=\lambda /{\omega }_{\star }~8$ which is similar to the value found in Zhu & Stone (2018). The magnetic lever arm can be used to estimate the terminal velocity of the outflow using $u_{p\infty }\simeq R_{\mathrm{SM}}{\mathrm{\Omega }}_{\mathrm{K}}(R_{\mathrm{SM}})\sqrt{2\lambda -3}$ (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 jet-like outflow. Since the laminar atmosphere is in ideal MHD, we can also compute the value of the mass loading invariant within it, finding κ ∼ 8. The value of κ is an order of magnitude larger within the laminar atmosphere that within the outflow. This is consistent with a laminar atmosphere that is so heavily loaded with mass that it falls towards the central object. It is that fall that leads to the generation of a large scale toroidal magnetic field within this region.

The Bernoulli invariant allows us to characterise the wind energetics and its drivers. 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.

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₀ = [5, 6, 8], and find very similar values, confirming the symmetric nature of the system.

Jacquemin-Ide et al. (2019) propose a criterion to distinguish between jets and winds, which essentially quantifies whether the outflow can self collimate thanks to magnetic forces (‘jet’) or not (‘wind’). They use the definition of the jet magnetisation

$$\begin{array}{c}\hfill \sigma ={\frac{-{\mathrm{\Omega }}_{\ast }R\langle B_{\phi }\rangle \langle B_p\rangle }{(\frac{{\langle u\rangle }^2}{2}+H)\langle \rho \rangle \langle u_p\rangle 4\pi }|}_{\mathrm{SM}},\end{array}$$(34)

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

$$\begin{array}{c}\hfill \frac{{\langle B_{\phi }\rangle }^2}{{\langle B_z\rangle }^2}>\frac{3}{4}\mathrm{with}\langle {\beta }_p\rangle \to 0\Rightarrow \mathrm{stability}.\end{array}$$(35)

This condition is verified in the laminar atmosphere (Fig. 4, bottom), indicating that it should be stable to the magneto-rotational 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.

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 re-ignition of turbulence is required for the establishment of a stationary magnetostatic equilibrium.

The re-activation of turbulence leads to anomalous resistivities within the turbulent atmosphere that allows accretion through the poloidal field lines and a steady state transition between accretion and ejection (Ferreira & Pelletier 1995). It is therefore an essential ingredient for launching the wind.

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 self-organisation. To understand the formation and evolution of the ring like structures, it is useful to vertically and azimuthally average the equation of mass conservation. Equation (1) leads to

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Sigma }}{\partial t}+\frac{1}{r^2}\frac{\partial r^2\overline{{\langle \rho u_r\rangle }_{\phi }}}{\partial r}+{[sin\theta {\langle \rho u_{\theta }\rangle }_{\phi }]}_{{\theta }_1}^{{\theta }_2}=0,\end{array}$$(36)

where $\mathrm{\Sigma }=\overline{{\langle \rho \rangle }_{\phi }}$ is the column density of the disk and the vertical and azimuthal average $\overline{X}$ is defined in Appendix B. We compute the vertical averages within the latitudinal wedge θ_2, 1 = π/2 ± arctan(3.5ϵ), equivalent to the disk region (Sect. 3.1). We show Σ in a space time diagram in the top panel of Fig. 11.

We concentrate our analysis on the gaps, minimum of Σ, located at R ≃ [6.6, 12.3] as well as the rings, maximum of Σ located at R ≃ [9.5, 15]. Indeed, the ring located at R ≃ 2, and to a certain extent the one located at R ≃ 4, are within close range of the inner boundary which could bias their analysis. Nonetheless, we find that the innermost ring-like structure (R ∼ 4) starts to emerge at around 200 Keplerian orbits of the innermost radii (25 local orbits) while the farther rings (at R ∼ 10) appear at 800 Keplerian orbits of the innermost radii (25 local orbits). Hence, we can conclude that the local time scale of self-organisation t_so(R) ≃ 25T_K(R) = 160/Ω_K(R).

Once formed, the ring structures survive for the entire duration of the simulation. They also exhibit a weak radial 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⁻⁴V_K(R) by measuring the radial velocity of its density maximum. This velocity is twice smaller than v_acc. It is unclear whether this drift is a bias of the boundary condition or a consequence of the intrinsic physics. It is important to note that while R ∼ 4 may seem to be migrating inwards we were not able to measure any radial drift.

The magnetic field is probably linked to the production of ring-like structures. To understand the role of the magnetic field in the self-organisation mechanism, we compute the vertically average vertical field $b_z=\overline{{\langle B_Z\rangle }_{\phi }}/(7h)$ and show its space time diagram in Fig. 11, bottom. We find that the vertical magnetic field is accumulating within the gaps, and is anti-correlated with the rings. However, the magnetic structure is not as stable with time as the density structure.

To understand if the self-organisation of the disk into ring-like 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, $\langle \mathrm{\Sigma }\rangle =\overline{\langle \rho \rangle }$, and the Keplerian angular velocity to find the dimensionless equation

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }\mathrm{\Sigma }}{\langle \mathrm{\Sigma }\rangle \mathrm{\Delta }t{\mathrm{\Omega }}_{\mathrm{K}}}+\frac{1}{\langle \mathrm{\Sigma }\rangle r^2{\mathrm{\Omega }}_{\mathrm{K}}}\frac{\partial r^2\overline{\langle \rho u_r\rangle }}{\partial r}+\frac{{\dot{\sigma }}_w}{\langle \mathrm{\Sigma }\rangle {\mathrm{\Omega }}_{\mathrm{K}}}=0\end{array}$$(37)

where we define ${\dot{\sigma }}_{\mathit{w}}={[sin\theta \langle \rho u_{\theta }\rangle ]}_{{\theta }_1}^{{\theta }_2}$ 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, ${\dot{\sigma }}_{\mathit{w}}$, is essentially negligible, while the radial term, ${\partial }_r{r}^2\overline{\langle \rho u_r\rangle }$ 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 self-organisation is due to α_tu.

From Fig. 12 we can probe the time scale, t_so, on which this self-organisation takes places by looking at the amplitude of the radial divergence of the radial velocity in units of keplerian frequency. This leads to t_so ∼ 200/Ω_K(R) = 30T_K(R) this corresponds quite well with the values measured in Fig. 11. Comparing the self organisation time scale with the ones deduced from Sect. 3.2.1 we find t_so ≪ t_acc ≃ t_Ψ. Self-organisation into ring-like structures is not impeded by the global transport mass or magnetic field, it happens on shorter time scales.

As mentioned above, the self-organisation is due to radial motions which are driven by the radial component of the turbulent stress tensor. It is tempting to interpret this as a viscous instability, as proposed by Lightman & Eardley (1974). Indeed, Lightman & Eardley (1974) showed that, for the system to be unstable to a viscous type instability, the torque needs to be a decreasing function of the surface density, $\frac{\partial log\overline{T_{\mathrm{tu},r\phi }}}{\partial log\mathrm{\Sigma }}<0$. In the case of a strongly diffusive field, the limit where the field is steady, this is equivalent to $\overline{{\alpha }_{\mathrm{tu}}}\propto {\langle {\beta }_{\mathrm{mid}}\rangle }^{-1}$, where we define

$$\begin{array}{c}\hfill \overline{{\alpha }_{\mathrm{tu}}}=\frac{\overline{T_{\mathrm{tu},r\phi }}}{\overline{P}},\end{array}$$(38)

and

$$\begin{array}{c}\hfill \langle {\beta }_{\mathrm{mid}}\rangle =\langle {\beta }_p\rangle (r,\theta =\pi /2).\end{array}$$(39)

However, in MRI turbulence, the measured scaling is closer to $\overline{{\alpha }_{\mathrm{tu}}}\propto {\langle {\beta }_{\mathrm{mid}}\rangle }^{-1/2}$ (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 mid-plane for radii ranging from R = 4 to R = 16. We expect $\overline{{\alpha }_{\mathrm{tu}}}$ to not only be a function of ⟨β_mid⟩, indeed in Fig. 13 we have the same value of $\overline{{\alpha }_{\mathrm{tu}}}$ for different plasma beta. This is probably a consequence of the dependence of $\overline{{\alpha }_{\mathrm{tu}}}$ on the local shear, d logΩ/d log R. Indeed, self-organisation induces local deviation from the Keplerian shear, which in turn influences the efficiency of the MRI-driven angular momentum transport (Gressel et al. 2015; Pessah et al. 2007), weaker shear leading to lower values of $\overline{{\alpha }_{\mathrm{tu}}}$. It can be seen that the scaling of the turbulent angular momentum transport coefficient, $\overline{{\alpha }_{\mathrm{tu}}}$, verifies a scaling closer to −1/2 than −1.

The −1/2 scaling is inconsistent with a simple viscous type instability, hence the self-organisation of our simulation into a ring-like system is not purely viscous, and must probably be somehow connected to magnetic field transport in the disk.

In conclusion, the fact that the self-organisation is driven by radial motions of matter tends to disfavour a wind-driven instability (Riols & Lesur 2019) and tends to favour some kind of viscously-driven instability (Hawley 2001; Bai & Stone 2014). However, as explained above, a simple viscous instability could not explain the disk instability since the scaling of $\overline{{\alpha }_{\mathrm{tu}}}$ on ⟨β_mid⟩ is not steep enough. Finally, we note that the radial density structure of the outflow is unaffected by the ring-like structures found in the turbulent disk. The radial structures disappear in the upper layers (LA or TA), probably because the rings are washed out by the supersonic accretion in those regions.

Bai & Stone (2014) construct a phenomenological model for ring formation, that does not require a steep scaling for $\overline{{\alpha }_{\mathrm{tu}}}$, based on the evolution of the vertical magnetic field driven by an anisotropic anti-diffusivity. However, it is not evident that an anisotropic anti-diffusion emerges from MRI turbulence. Furthermore, recent measures of the mean field diffusivities of MRI turbulence do not agree with a negative anisotropic diffusivity (Gressel et al. 2015).

Finally, Johansen et al. (2009) explore the possibility that stochastic perturbations driven by an inverse turbulent cascade of the magnetic components could be driving the self-organisation of the disk into ring-like structures as well as the formation of pressure bumps. However, their modelling concentrates on unstratified shearing box simulations with a zero-net flux. This leads to their ring-like structures being short lived, lasting at most 50 local orbits while our inner ring-like structure are still stable at 75 local orbits and show no sign of weakening. Nonetheless, a stochastic model taking into account the feedback of the vertical field evolution may be needed to explain the ring formation observed here.

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₁ and t₂ = 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 super-sonic accretion and a super-fast wind. Therefore we address the weak-field simulations first, while the steady state of simulation SB2 as well as its properties will be described at the end of this section.

5.1. Weak-field simulations: S2pi, SB3, SEp

We start by validating the wedge approximation, Δφ < 2π, for global simulations, S2pi harbours no significant differences when compared to SB4. Indeed, since S2pi and SB4 are identical with respect to the properties of their vertical structure, we ignore S2pi while we discuss such properties.

We confirm that similarly to simulation SB4, all weak-field simulations reach a equilibrium within the disk and a magnetostatic equilibrium within laminar and turbulent atmospheres (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.

The altitude of the jet launching surface is tightly related to the density profile. In Fig. 15 we show the density profiles of the different simulations normalised to their value at the disk mid-plane. We see that the simulations SB4 and S2pi converge to a similar density configuration. Moreover, we observe in Fig. 15 that as β_ini decreases the density profile becomes shallower and as a consequences more mass is fed into the TA and the outflow. This is possibly related to the enhanced mass loading fulfilled by the turbulent magnetic pressure. Indeed, decreasing β_ini leads to a decrease of ⟨β_mid⟩ which increases the turbulence strength (Salvesen et al. 2016; Mishra et al. 2020). This increase in turbulent strength with the decrease of the plasma beta is also verified in our simulations.

From 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⟩⁻¹. Moreover, v_Ψ/V_K strongly decreases when we decrease the disk geometrical thickness.

In Fig. 16 (right) we show v_Ψ/v_acc as a function of the radial coordinate. We observe that, except for the inner regions of SB2, all simulations follow approximately the same dependency,

$$\begin{array}{c}\hfill v_{\mathrm{\Psi }}\simeq v_{\mathrm{acc}}.\end{array}$$(40)

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₀  =  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 weak-field simulations, the invariants should naturally be independent of ⟨β_mid⟩.

Ring-like structures are also found in all weak-field simulations. They are always a consequence of converging radial mass flows. In general the rings in the weak-field simulations feature a similar time-scale of self-organisation, t_so ∼ 200/Ω_K(R) = 32T_K(R), consistent with the one measured in SB4. However, there are some quantitative differences:

• The ring-like structures in run SB3 are accreted towards the inner region a lot faster than for the fiducial run SB4, v_ring ≃ 3 × 10⁻³V_K(R). This velocity is around a factor of 3 smaller than the local v_acc(R) in this simulation. Moreover, the accreted rings are sufficiently far away (R ∼ 4 and R ∼ 6) from the inner boundary for us to believe that radial migration is not a bias of the inner boundary, in contrast with SB4. The ratios between the maxima and minima of density are very similar to the ones in the fiducial run, Σ_max/Σ_min ≃ 2.

• The ring-like structures in run S2pi are less pronounced than in run SB4, Σ_max/Σ_min ≃ 1.5. This may be related to the limited toroidal extent of simulation SB4. Similarly to SB4, the radial migration of the ring-like structures for S2pi is inexistent or negligible.

• The ring-like structures in run SEp are very similar to the ones found in run SB4, the ratios between minima and maxima are similar, Σ_max/Σ_min ≃ 2. However, run SEp features many more ring-like structures, by about a factor 2. This is consistent with a viscous like instability bounded at small scale by the disk scale-height, in which the spacing between ring-like structures is expected to scale like h. Finally, radial migration of ring-like structures is also negligible in this simulation.

The Rossby wave instability (RWI) is known to produce non-axisymmetric structures, like Rossby vortices or spiral shocks which smear out local density maxima. The RWI is triggered² when the inverse of the potential vorticity

$$\begin{array}{c}\hfill \mathcal{F}(r)=\frac{\langle \mathrm{\Sigma }\rangle }{\mathbf{\nabla }\times \langle \mathit{u}\rangle {|}_{\theta }},\end{array}$$(41)

has a local maximum (Lovelace et al. 1999). To check whether our ring-like structures are unstable to the RWI, we show the inverse of potential vorticity as a function of radii in Fig. 17 for simulations S2pi and SB4. We find that for SB4, the local maxima of the inverse of the vorticity potential are well correlated with the local maxima of the column density profile. This correlation is less clear for run S2pi, as the ring-like structures are less pronounced. One would therefore conclude that SB4 rings should be RWI unstable and therefore develop non-axisymmertic structures. However, no clear spiral structures are visible in SB4, while non-axisymmetric features similar to the ones shown by Mishra et al. (2020) are appreciable in S2pi. It is therefore temping to think that the RWI in SB4 is stabilised by some mechanism. One possibility is that the RWI is stabilised because of turbulent diffusion which damps RWI modes. Another one is that the limited azimuthal extension of SB4 (Δφ = π/2) filters out some of the most destructive low-m modes excited by the RWI. The presence of strong spiral waves in S2pi is in favour of the second hypothesis, but a proper study of the ring stability is needed to ascertain this conclusion, which is outside the scope of this work.

5.2. Strong field simulation: SB2

The main difference between SB2 and all weak-field simulations is that SB2 takes a much longer time to reach a reasonable steady state, the whole disk being subject to a drastic 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 quasi-stationary 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⁻³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⟩⁻¹ increases the amount of mass lifted by the turbulent magnetic pressure increases, leading to a massive turbulent atmosphere that merges with the disk.

Figure 18 shows a zoomed-in outcome of our SB2 simulation, averaged between t₁ = 1719T_in and t₂ = 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 weak-field simulations. Below R_J ∼ 7, accretion occurs within the turbulent disk and the mass weighted accretion speed is supersonic. The disk rotation is near keplerian with ⟨u_φ⟩ ≃ 0.8v_K, namely a deviation of the order of ϵ as expected in near equipartition disks. A cold outflow is launched from the disk surface at z ≃ 3.5h, becoming soon a super Alfvénic (red dotted curve) and then super fast magnetosonic (red dashed curve) jet.

Overall, this inner accretion-ejection structure shares most properties with the Jet Emitting Disk model (Ferreira 1997). The main difference is the amount of ejected mass as the outflows obtained here are quite massive (Fig. 15). As discussed above, this is a consequence of the turbulence itself since the turbulent magnetic pressure plays a major role in the magnetostatic vertical equilibrium and launching of the outflow. Finally, we report the lack of ring-like structures in our strong field simulation. This is presumably due to the fact that accretion is now too fast, t_acc ≪ t_so, with the t_so as obtained from the weak field simulations. Mass has no time to accumulate locally and lead to the formation of rings and gaps.

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 accretion-ejection configuration to achieve a reasonable steady-state. We found that the disk structure is drastically dependent on its magnetisation ⟨β_mid⟩⁻¹.

(1) At disk magnetisations lower than 10⁻³, a non trivial vertical stratification sets in on a scale that can go up to z ∼ 10h, where an outflow that becomes super fast magnetosonic is launched. The disk is highly turbulent, mostly in hydrostatic equilibrium and the presence of an important magnetic turbulent pressure allows to lift up a strong amount of mass into a laminar ideal MHD atmosphere. This lifted mass is basically falling in towards the central object, dragging in the magnetic field which, in turn, provides a torque that transfers its angular momentum back to the underlying disk. When the magnetic field achieves equipartition in these upper layers, MRI is re-ignited and drives a turbulence. In this levitating turbulent atmosphere, transsonic accretion is achieved mostly through the turbulent torque and a small fraction of the mass is then loaded into the open field lines where a magnetically driven cold wind is launched. In contrast with shearing box simulations, the bipolar wind configuration is found to be fairly symmetric. We report also the formation in the disk of rings and gaps, on time scales shorter than the accretion time scale.

(2) The ambient magnetic field is found to be always dragged inwards in the disk, at a velocity which increases with the disk 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 mid-plane properties.

(3) However, beyond a threshold on the disk magnetisation, located between 10⁻³ and 10⁻², the global accretion-ejection configuration undergoes a drastic readjustment. The magnetic field is being accumulated into the central regions until a global equilibrium is achieved involving both the disk and its magnetosphere. The inner disk reaches a steady-state (balance between field advection and diffusion) when its magnetisation achieves unity and no more gaps and rings are found. Despite this strong field regime, the disk is turbulent and drives a super fast-magnetosonic outflow right from its surface z ∼ 3h. The size of this inner region keeps on increasing in time, as more magnetic flux is being added from the outer regions. The mass-weighted accretion speed is supersonic in the inner region. The transition region, where the disk material must become transsonic, is quite complex and its investigation is postponed for future work.

Some of our results have already been obtained or discussed in the literature. Previous numerical works (Zhu & Stone 2018; Mishra et al. 2020) found the same stationary states as in our simulations. However, these authors did not analyse the appearance of the second upper turbulent zone and how its properties are related to the vertical structure of the system. Another major difference is related to the magnetic field advection. Zhu & Stone (2018) report a small amount of magnetic flux accumulation through their inner boundary albeit without providing a clear analysis. We show that the field advection speed v_Ψ follows a Keplerian scaling and increases with the disk 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 accretion-ejection structures. While it is widely believed that the MRI can only exist in high β plasmas, Kim & Ostriker (2000) have shown that the MRI exist for arbitrary β provided that arbitrary long wavelength perturbations are allowed. In the β < 1 regime, the MRI can be quenched when the toroidal field gets too strong. We have shown that the turbulent and laminar regions of the atmosphere, which is a β < 1 region, is naturally explained by the MRI in this peculiar regime.

Magnetic field dragging seems to be a generic property of ideal MHD accretion disks threaded by a large scale magnetic field (provided that ϵ ∼ 0.1 or larger). Indeed, it is also reported in GRMHD simulations of accretion disks around a Kerr black hole (e.g., McKinney et al. 2012; White et al. 2019 and references therein, see also Liska et al. 2020). In these simulations, the magnetic field is seen to be gradually advected and concentrated into the central region, leading to the build up of what has been termed a Magnetically Choked Accretion Flow (McKinney et al. 2012). As in our SB2 simulation, magnetic field dragging stops when the vertical field reaches an equipartition value, defining a transition radius with an outer disk at lower magnetisation. As more magnetic flux is being advected, a violent relaxation occurs at this radius involving a magnetic Rayleigh-Taylor instability. It leads to the expulsion of the magnetic flux excess and a progressive increase in time of the transition radius, as long as some magnetic flux remains available in the simulation (McKinney et al. 2012). This situation is quite nicely consistent with the hybrid disk configuration proposed for transition disks (Combet & Ferreira 2008; Wang & Goodman 2017) and around black holes (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 X-ray 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⁻³) provide a highly subsonic mass-weighted accretion speed as expected in a SAD. However they differ from the SAD model in two important aspects: (i) the presence of outflows and (ii) the existence of an upper turbulent accreting zone.

As stated above, the properties of the super fast-magnetosonic 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⟩⁻¹ (see Table 3). The best observational constraints on jet kinematics are provided by proto-planetary 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₀ ≃ 1.5 au, while de Valon et al. (2020) measure λ ≃ 1.6 and R₀ ≃ 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 proto-planetary disks, may change the outflow properties (Béthune et al. 2017).

The second important difference with the Shakura-Sunyaev SAD model is the existence of a turbulent accreting zone levitating above a weakly 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 proto-planetary disks (Gammie 1996; Fleming & Stone 2003; Bai & Stone 2013). However, the physical reason is completely different. In our case, the disk is fully ionised and turbulent whereas in the dead zone picture non-ideal effects quench turbulence and lead to the formation of a MRI stable (‘dead’) zone in the disk. The interesting part is that, in both cases and as long as the vertical field magnetic is weak, accretion takes place at the disk ‘surface’. This is a feature that is drastically different from the conventional view of Shakura-Sunyaev accretion disks. Moreover, this may lead to direct observational consequences on the emitted spectra of weakly 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 low-density region. In the context of compact objects, we expect the production of an optically thin emission, usually referred to as a coronal emission. This ‘elevated accretion’ disk configuration could thus provide an explanation for the ‘warm corona’ component often seen in soft X rays in some active galactic nuclei (Done et al. 2012; Petrucci et al. 2018).

It is important to note that, to make the problem numerically tractable we have decided to ignore the effects of the energy equation. We did this by forcing a locally isothermal structure with a fixed disk scale, ϵ = 0.1 or ϵ = 0.05. Optically thick disks in AGN and X-ray binaries are expected to be geometrically thin with ϵ = 0.01 or less, which is incompatible with our simulation setup. Nonetheless, we believe that the accreting atmosphere should perdure at lower ϵ, as long as the disk 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 thermo-magnetic winds is postponed for future work.

The parameter space explored in our simulations overlaps with the one explored by Jacquemin-Ide et al. (2019) using self-similar solutions. The outflow invariants computed in our simulations (Table 3) are found consistent with the borders of the parameter space they derive semi-analytically. However, while in our simulations λ ∼ 4−5 is roughly constant, they obtain values ranging from 2 to 100. As shown in their Figs. 4 and 5, the magnetic lever arm is independent of the disk 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⁻³ 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 self-consistently computed and must therefore be prescribed: self-similar studies (Ferreira 1997; Jacquemin-Ide et al. 2019), 2.5D MHD simulations done with alpha prescriptions (Murphy et al. 2010; Stepanovs & Fendt 2016) and other semi-analytical models (Guilet & Ogilvie 2013). Our guess is that including the complex vertical stratification should allow to recover with these approaches an outflow mass load consistent with full 3D simulations. This requires however the use of numerically ‘educated’ profiles for all relevant turbulent effects. This is postponed for future work.

Finally, weakly magnetised disks are subject to a secular instability which leads to the formation of long-lived ring-like structures in the disk surface density distribution, the magnetic field being concentrated in the gaps. These ring-like structures are also visible in the surface density profiles of Zhu & Stone (2018) and, like for us, they also seem to increase in number when the disk geometrical thickness is decreased. This is consistent with a viscous like instability bounded at small scale by the disk scale-height, for which the spacing between ring-like structures is expected to scale like h. On the other hand, we showed that this global self-organisation is unlikely due to the wind instability proposed by Riols & Lesur (2019). Yet, it remains unclear if it is purely a viscous instability, as the turbulent viscosity scaling ${\overline{\alpha }}_{\mathrm{tu}}(\langle {\beta }_{\mathrm{mid}}\rangle )$ is not steep enough to develop a viscous instability. When the disk magnetisation increases above the threshold, these ring-like structures are smeared out and disappear as the accretion speed becomes supersonic. Our suspicion is that the ring-formation mechanism and the magnetic field transport could be related to each other. Finally, we show that the limited toroidal extent of our simulation SB4 and SB3 could have an impact on the long-term evolution of the ring-like structures. Indeed, it is not clear if the rings are destroyed by non-axis-symmetric 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⁻³ − 10⁻², ring-like structures should be present, and absent in disks with equipartition magnetic fields. It could thus be an observationally convenient way to probe the existence and strength of large scale magnetic fields in proto-planetary disks (see for instance de Valon et al. 2020).

---

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.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.

References

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

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }t}{\left[{\langle \delta \rho \rangle }_{\phi }\right]}_{t_1}^{t_2}+\mathbf{\nabla }·\langle \rho \rangle [\langle \mathit{u}\rangle +{\langle \delta \mathit{u}\rangle }_{\rho }]=0.\end{array}$$(A.1)

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 non-steadiness 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):

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }t}{\left[{\langle \delta B_{\phi }\rangle }_{\phi }\right]}_{t_1}^{t_2}={\mathbf{\nabla }\times [\langle \mathit{u}\rangle \times \langle \mathit{B}\rangle +\mathcal{E}]|}_{\phi }\end{array}$$(A.2)

where the turbulent emf is defined as

$$\begin{array}{c}\hfill \mathcal{E}=\langle \mathit{\delta }\mathit{u}\times \mathit{\delta }\mathit{B}\rangle .\end{array}$$(A.3)

We write the poloidal component of the induction equation as the toroidal component of the equation for the vector potential

$$\begin{array}{c}\hfill \frac{\partial A_{\phi }}{\partial t}={\mathit{u}}_{\mathit{p}}\times {\mathit{B}}_{\mathit{p}}+\frac{1}{R}\frac{\partial G}{\partial \phi },\end{array}$$(A.4)

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

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }t}{\left[{\langle \delta A_{\phi }\rangle }_{\phi }\right]}_{t_1}^{t_2}=\langle {\mathit{u}}_{\mathit{p}}\rangle \times \langle {\mathit{B}}_{\mathit{p}}\rangle +{\mathcal{E}}_{\phi },\end{array}$$(A.5)

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:

$$\begin{array}{cc}& \frac{\partial r\delta u_{\theta }\rho }{\partial t}+\rho {\mathit{u}}_{\mathbf{p}}·\mathbf{\nabla }r\langle u_{\theta }\rangle +\mathbf{\nabla }·r[\rho \delta u_{\theta }\mathit{u}-\frac{1}{4\pi }{B}_{\theta }\mathit{B}]\hfill \\ \hfill & =cot\theta [\rho u_{\phi }^2-\frac{1}{4\pi }{B}_{\phi }^2]-\frac{\partial }{\partial \theta }(\frac{B^2}{8\pi }+P),\hfill \end{array}$$(A.6)

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:

$$\begin{array}{cc}& \frac{1}{\mathrm{\Delta }t}{\left[r{\langle \rho \delta u_{\theta }\rangle }_{\phi }\right]}_{t_1}^{t_2}+\langle \rho \rangle {\langle {\mathit{u}}_{\mathit{p}}\rangle }_{\rho }·\mathbf{\nabla }r\langle u_{\theta }\rangle -\frac{1}{4\pi }\langle {\mathit{B}}_{\mathit{p}}\rangle ·\mathbf{\nabla }r\langle B_{\theta }\rangle \hfill \\ \hfill & +\mathbf{\nabla }·r[\langle \rho \rangle {\langle \delta u_{\theta }{\mathit{u}}_{\mathit{p}}\rangle }_{\rho }-\frac{1}{4\pi }\langle \delta B_{\theta }\mathit{\delta }{\mathit{B}}_{\mathit{p}}\rangle ]=\hfill \\ \hfill & -\frac{\partial }{\partial \theta }[\langle P\rangle +\frac{1}{8\pi }({\langle B\rangle }^2+\langle \delta B^2\rangle )]\hfill \\ \hfill & +cot\theta [\langle \rho \rangle ({\langle u_{\phi }\rangle }^2+\langle \delta u_{\phi }^2\rangle )-\frac{1}{4\pi }({\langle B_{\phi }\rangle }^2+\langle \delta B_{\phi }^2\rangle )]\hfill \end{array}$$(A.7)

where ⟨u⟩_ρ = ⟨u⟩ + ⟨δu⟩_ρ and

$$\begin{array}{c}\hfill {\langle \delta u_{\theta }{\mathit{u}}_{\mathit{p}}\rangle }_{\rho }={\langle \delta u_{\theta }\delta {\mathit{u}}_{\mathit{p}}\rangle }_{\rho }+\langle {\mathit{u}}_{\mathit{p}}\rangle {\langle \delta u_{\theta }\rangle }_{\rho }.\end{array}$$(A.8)

In the derivation of Eq. (A.7) we have used

$$\begin{array}{c}\hfill \langle \mathbf{\nabla }·\mathit{B}{B}_i\rangle =\langle {\mathit{B}}_{\mathit{p}}\rangle ·\mathbf{\nabla }\langle B_i\rangle +\mathbf{\nabla }·\langle \delta {\mathit{B}}_{\mathit{p}}\delta B_i\rangle .\end{array}$$(A.9)

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

$$\begin{array}{cc}& \frac{\partial \delta u_r\rho }{\partial t}+\rho {\mathit{u}}_{\mathit{p}}·\mathbf{\nabla }\langle u_r\rangle +\mathbf{\nabla }·[\rho \delta u_r\mathit{u}-\frac{1}{4\pi }{B}_r\mathit{B}]\hfill \\ \hfill & =\frac{1}{r}[\rho (u_{\phi }^2+u_{\theta }^2)-\frac{1}{4\pi }(B_{\phi }^2+B_{\theta }^2)]-\frac{\partial }{\partial r}(\frac{B^2}{8\pi }+P)-\rho g,\hfill \end{array}$$(A.10)

as before we then average the equation and find:

$$\begin{array}{cc}& \frac{1}{\mathrm{\Delta }t}{\left[{\langle \rho \delta u_r\rangle }_{\phi }\right]}_{t_1}^{t_2}+\langle \rho \rangle {\langle {\mathit{u}}_{\mathit{p}}\rangle }_{\rho }·\mathbf{\nabla }\langle u_r\rangle -\frac{1}{4\pi }\langle {\mathit{B}}_{\mathit{p}}\rangle ·\mathbf{\nabla }\langle B_r\rangle \hfill \\ \hfill & +\mathbf{\nabla }·[\langle \rho \rangle {\langle \delta u_r{\mathit{u}}_{\mathit{p}}\rangle }_{\rho }-\frac{1}{4\pi }\langle \delta B_r\mathit{\delta }{\mathit{B}}_{\mathit{p}}\rangle ]=-\frac{\partial }{\partial r}[\langle P\rangle +\frac{1}{8\pi }({\langle B\rangle }^2+\langle \delta B^2\rangle )]\hfill \\ \hfill & +\frac{\langle \rho \rangle }{r}[{\langle u_{\phi }\rangle }^2+\langle \delta u_{\phi }^2\rangle +{\langle u_{\theta }\rangle }^2+\langle \delta u_{\theta }^2\rangle ]-\langle \rho \rangle g\hfill \\ \hfill & -\frac{1}{4\pi r}[{\langle B_{\phi }\rangle }^2+\langle \delta B_{\phi }^2\rangle +{\langle B_{\theta }\rangle }^2+\langle \delta B_{\theta }^2\rangle ].\hfill \end{array}$$(A.11)

We project and rearrange the equation for toroidal momentum transport to find the usual equation of angular momentum transport:

$$\begin{array}{c}\hfill \frac{\partial R{u}_{\phi }\rho }{\partial t}+\mathbf{\nabla }·R[\rho u_{\phi }\mathit{u}-\frac{1}{4\pi }{B}_{\phi }\mathit{B}]=-\frac{\partial }{\partial \phi }(\frac{B^2}{8\pi }+P)\end{array}$$(A.12)

we use the conservation of mass and then repeat the averaging procedure to find:

$$\begin{array}{cc}& \frac{1}{\mathrm{\Delta }t}{\left[R{\langle \rho \delta u_{\phi }\rangle }_{\phi }\right]}_{t_1}^{t_2}+\langle \rho \rangle {\langle {\mathit{u}}_{\mathit{p}}\rangle }_{\rho }·\mathbf{\nabla }R\langle u_{\phi }\rangle +\hfill \\ \hfill & \mathbf{\nabla }·R[\langle \rho \rangle {\langle \delta u_{\phi }{\mathit{u}}_{\mathit{p}}\rangle }_{\rho }-\frac{1}{4\pi }\langle \delta B_{\phi }\mathit{\delta }{\mathit{B}}_{\mathit{p}}\rangle -\frac{1}{4\pi }\langle B_{\phi }\rangle \langle {\mathit{B}}_{\mathit{p}}\rangle ]=0.\hfill \end{array}$$(A.13)

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

$$\begin{array}{c}\hfill \overline{X}={\displaystyle \underset{{\theta }_1}{\overset{{\theta }_2}{\int }}rsin\theta X\mathrm{d}\theta .}\end{array}$$(B.1)

We can then integrate Eq. (1) with respect to time and space to find

$$\begin{array}{c}\hfill \mathrm{\Delta }M=F_{r_2}-F_{r_1}+F_{{\theta }_1}-F_{{\theta }_2}\end{array}$$(B.2)

where we have introduced

$$\begin{array}{cc}\hfill F_{r_i}=& -{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{r}_i{\left[\overline{{\langle \rho u_r\rangle }_{\phi }}\right]|}_{r_i}\mathrm{d}t,}\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill F_{{\theta }_i}=& {\displaystyle \underset{t_1}{\overset{t_2}{\int }}\underset{r_1}{\overset{r_2}{\int }}rsin\theta {{\langle \rho u_{\theta }\rangle }_{\phi }|}_{{\theta }_i}\mathrm{d}r\mathrm{d}t}\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill \mathrm{\Delta }M=& {\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\left[r\overline{{\langle \rho \rangle }_{\phi }}\right]}_{t_1}^{t_2}\mathrm{d}r,}\hfill \end{array}$$(B.5)

we can redo the same exercice for Eq. (A.12) to find

$$\begin{array}{c}\hfill F_{M_r}+F_{M_{\theta }}=F_{r_2}-F_{r_1}+F_{{\theta }_1}-F_{{\theta }_2}.\end{array}$$(B.6)

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

$$\begin{array}{c}\hfill {}_{\phi }=u_{\phi }-R\omega (r,\theta ),\end{array}$$(B.7)

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

$$\begin{array}{cc}\hfill F_{r_i}=& {\displaystyle \underset{t_1}{\overset{t_2}{\int }}{r}_i^2{[\overline{{\langle \rho v_{\phi }{u}_r\rangle }_{\phi }sin\theta }-\overline{\frac{{\langle B_{\phi }{B}_r\rangle }_{\phi }sin\theta }{4\pi }}]|}_{r_i}\mathrm{d}t,}\hfill \end{array}$$(B.8)

$$\begin{array}{cc}\hfill F_{{\theta }_i}=& {\displaystyle \underset{t_1}{\overset{t_2}{\int }}\underset{r_1}{\overset{r_2}{\int }}{r}^2{sin}^2\theta {[{\langle \rho u_{\theta }{v}_{\phi }\rangle }_{\phi }-\frac{1}{4\pi }{\langle B_{\phi }{B}_{\theta }\rangle }_{\phi }]|}_{{\theta }_i}\mathrm{d}r\mathrm{d}t,}\hfill \end{array}$$(B.9)

$$\begin{array}{cc}\hfill F_{M_r}=& {\displaystyle \underset{t_1}{\overset{t_2}{\int }}\underset{r_1}{\overset{r_2}{\int }}r\overline{{\langle \rho u_r\rangle }_{\phi }}\frac{\partial \omega r^2{sin}^2\theta }{\partial r}\mathrm{d}r\mathrm{d}t,}\hfill \end{array}$$(B.10)

$$\begin{array}{cc}\hfill F_{M_{\theta }}=& {\displaystyle \underset{t_1}{\overset{t_2}{\int }}\underset{r_1}{\overset{r_2}{\int }}\overline{{\langle \rho u_{\theta }\rangle }_{\phi }}\frac{\partial \omega r^2{sin}^2\theta }{\partial \theta }\mathrm{d}r\mathrm{d}t.}\hfill \end{array}$$(B.11)

We can simplify the term F_Mθ by choosing ω = Ω_K/sin²θ. 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₁ = 8 and r₂ = 14 (see Fig. 3, top). We then integrate Eq. (1) and get

$$\begin{array}{c}\hfill \mathrm{\Delta }M=F_{m,r_2}-F_{m,r_1}+F_{m,{\theta }_1}-F_{m,{\theta }_2},\end{array}$$(B.12)

where the mass fluxes are time-averaged on a duration Δt which spans several local orbital periods. A strict steady state would translate into ΔM = 0. To illustrate time-variability, 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 quasi-periodic, 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 MRI-driven turbulence (Heinemann & Papaloizou 2009a,b).

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₁ = 923T_in to t₂ = 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 ring-like structures within the disk region, a form of secular self-organisation. The domain, between r₁ = 8 and r₂ = 14, is located on a two distinct gaps, one at R = 9.5 and another at R = 15. The formation of ring-like structures is discussed in Sect. 4.

Turning now to the regions above the disk, we find that the laminar atmosphere is accreting, but less efficiently than the turbulent atmosphere, despite the fact that both regions receive approximately the same amount of mass from the outer radial boundary. This trend can be understood from the shape of the stream lines (Fig. 2), since the mass transported within the laminar atmosphere inevitably ends up being accreted in the turbulent atmosphere at a smaller radius. This is confirmed by Fig. B.2 which shows that the turbulent atmosphere is being fed mass by the laminar atmosphere. Eventually, most of the mass received by the turbulent atmosphere is radially accreted, the mass ejected being almost negligible. 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

$$\begin{array}{c}\hfill F_{L,M_r}=F_{L,r_2}-F_{L,r_1}+F_{L,{\theta }_1}-F_{L,{\theta }_2},\end{array}$$(B.13)

where the fluxes definition can be found in Appendix B. The left-hand side denotes angular momentum transported inwards by the accreting flow, while the right-hand side measure the angular momentum flux due to the torques (turbulent and laminar). The angular momentum fluxes are quasi-stationary when compared to the mass fluxes. Hence, we skip the examination of the time variability of these fluxes and average them across the whole time domain (t₁ = 5800/Ω_K(R_in) to t₂ = 7000/Ω_K(R_in)). We summarise their values in Fig. B.3.

We can see that the angular momentum is transported radially outwards by the torques in all of the regions. Surprisingly, in the latitudinal direction, angular momentum is transported downwards: from the atmosphere to the disk. The disk is therefore azimuthally accelerated by the atmosphere, as already seen by Zhu & Stone (2018). Overall, this latitudinal flux of angular momentum implies that the disk is gaining angular momentum, hence decretion must occur³. 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 MRI-turbulent 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 ρ₀, pressure $c_{\text{s}}^2{\rho }_0$ (where c_s is the isothermal sound speed) and magnetic field B₀. Take a fiducial radius R₀ and attach a cartesian frame (x, y, z) centred on R₀ and rotating at Ω = Ω(R₀) so that x corresponds to the local radial direction R − R₀ and y the local azimuthal direction φ. Let us then add a small perturbation to this plasma so that ρ = ρ₀ + δρ and B = B₀ + δ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

$$\begin{array}{cc}\hfill {\rho }_0\frac{D^2\mathit{\xi }}{D{t}^2}& =-c_{\mathrm{s}}^2\mathbf{\nabla }\delta \rho -\frac{1}{4\pi }\mathbf{\nabla }({\mathit{B}}_0·\delta \mathit{b})+\frac{1}{4\pi }{\mathit{B}}_0·\mathbf{\nabla }\delta \mathit{b}\hfill \\ \hfill & -2{\rho }_0\mathbf{\Omega }\mathbf{\times }\frac{D\mathit{\xi }}{\mathit{Dt}}-2{\rho }_0\mathrm{\Omega }\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}logR}{\xi }_x{\mathit{e}}_x,\hfill \end{array}$$(C.1)

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)

$$\begin{array}{cc}\hfill \delta \rho & =-{\rho }_0\mathbf{\nabla }·\mathit{\xi },\hfill \\ \hfill \delta \mathit{b}& ={\mathit{B}}_0·\mathbf{\nabla }\mathit{\xi }-{\mathit{B}}_0(\mathbf{\nabla }·\mathrm{\nabla }\xi ),\hfill \end{array}$$

which can be combined to the equation of motion on ξ to get

$$\begin{array}{cc}\hfill \frac{D^2\mathit{\xi }}{D{t}^2}& =(c_{\mathrm{s}}^2+V_{\mathrm{A}}^2)\mathbf{\nabla }(\mathbf{\nabla }·\mathit{\xi })-\mathbf{\nabla }[{\mathit{V}}_{\mathrm{A}}·\mathbf{\nabla }({\mathit{V}}_{\mathrm{A}}·\mathit{\xi })]\hfill \\ \hfill & +{\mathit{V}}_{\mathrm{A}}·\mathbf{\nabla }[{\mathit{V}}_{\mathrm{A}}·\mathbf{\nabla }\mathit{\xi }]-{\mathit{V}}_{\mathrm{A}}[{\mathit{V}}_{\mathrm{A}}·\mathbf{\nabla }(\mathbf{\nabla }·\mathit{\xi })]\hfill \\ \hfill & -2\mathbf{\Omega }\mathbf{\times }\frac{D\mathit{\xi }}{\mathit{Dt}}+2q{\mathrm{\Omega }}^2{\xi }_x{\mathit{e}}_x.\hfill \end{array}$$

In this equation, we have introduced the Alfvén velocity vector ${\mathit{V}}_{\text{A}}\equiv {\mathit{B}}_0/\sqrt{4\pi {\rho }_0}$ 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₀ = B_0, φe_y + B_0, ze_z, and perturbations with only a vertical spatial dependency. We then seek for growing (unstable) solutions of the form $\mathit{\xi }(z,t)=\widehat{\mathit{\xi }}exp(\gamma t+ikz)$. One eventually gets the dispersion relation

$$\begin{array}{cc}\hfill c_{\mathrm{s}}^2{k}^2[{\gamma }^4+(2V_{\mathrm{A}z}^2{k}^2+{\kappa }^2){\gamma }^2+V_{\mathrm{A}z}^2& k^2(V_{\mathrm{A}z}^2{k}^2-2q{\mathrm{\Omega }}^2)]\hfill \\ \hfill +{\gamma }^6+{\gamma }^4[(V_{\mathrm{A}}^2& +V_{\mathrm{A}z}^2)k^2+{\kappa }^2]\hfill \\ \hfill +{\gamma }^2{V}_{\mathrm{A}}^2{k}^2& (V_{\mathrm{A}z}^2{k}^2+4{\mathrm{\Omega }}^2\frac{V_{\mathrm{A}\phi }^2}{V_{\mathrm{A}}^2}-2q{\mathrm{\Omega }}^2)=0\hfill \end{array}$$(C.2)

where we have introduced the epicyclic frequency κ² ≡ 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⁴ the usual dispersion relation for the MRI in the incompressible limit (Balbus & Hawley 1991). The stability condition is then given by the constant term: $V_{\text{A}z}{}^2{k}^2-2q{\mathrm{\Omega }}^2<0\to \text{instability}$. 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 $V_{\text{A}z}{}^2{k}^2+4{\mathrm{\Omega }}^2{V}_{\text{A}\phi }{}^2/V_{\text{A}}{}^2-2q{\Omega }^2<0\to \text{instability}$. 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 $V_{\text{A}{V}_{\text{A}}rphi}{}^2/V_{\text{A}}{}^2<q/2$. 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 space-time diagrams of the magnetic flux Ψ (defined in Sect. 3.2.1) and, for practical reasons, the ratio P_b, mid/P_mid, where

$$\begin{array}{cc}\hfill P_{b,\mathrm{mid}}& =\frac{1}{8\pi }{[\frac{1}{2h_{\mathrm{disk}}}{\displaystyle \underset{{\theta }_{d1}}{\overset{{\theta }_{d2}}{\int }}}rsin\theta {\langle B_z\rangle }_{\phi }\mathrm{d}\theta ]}^2,\hfill \end{array}$$(D.1)

$$\begin{array}{cc}\hfill P_{\mathrm{mid}}& =\frac{1}{2h_{\mathrm{disk}}}{\displaystyle \underset{{\theta }_{d1}}{\overset{{\theta }_{d2}}{\int }}rsin\theta {\langle P\rangle }_{\phi }\mathrm{d}\theta ,}\hfill \end{array}$$(D.2)

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 ${\beta }_{\text{mid}}^{-1}$ since it contains the turbulent field contribution ⟨δB²⟩.

While the turbulent magnetic pressure measured at the mid plane in weakly magnetised disks remains always sub-thermal (see Fig. 7), it becomes supra-thermal at large magnetisations. For instance, in the case of a disk with ⟨β_mid⟩⁻¹ ≃ 1, we measure ⟨δB²⟩(r, θ = π/2)/8π ∼ 10P₀. 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⟩⁻¹ is around unity.

All Tables

All Figures