© G. R. J. Lesur 2021

1 Introduction

Protoplanetary discs are relatively cold and dense objects, which typically last for a few million years around young stellar objects (YSOs). However, it is well known that these discs are primarily ‘accretion discs’, in which matter is slowly falling onto the central star. This accretion rate has now been measured in dozens of objects, and it typically lies in the range of 10⁻¹⁰–10⁻⁷ M_⊙ yr⁻¹. Accretion in astrophysical discs is usually thought to be due to a magnetohydrodynamic (MHD) instability, the magneto-rotational instability (MRI, Balbus & Hawley 1991), which transports angular momentum outwards and mass inwards. It was quickly realised that applying the MRI in the context of protoplanetary discs is a difficult task because of the dramatically low ionisation fraction expected in these objects. The inclusion of Ohmic (Gammie 1996) and ambipolar diffusions (Perez-Becker & Chiang 2011) led to the conclusion that MRI cannot account for the observed accretion rates in these discs by two to three orders of magnitude.

The lack of a proper mechanism to trigger turbulence and angular momentum transport in protoplanetary discs revived an old idea that was mostly abandoned since the finding of the MRI: magnetised winds. While the disc bulk might be too diffusive to sustain MHD turbulence, it is still weakly coupled to the ambient large-scale field originating from the parent cloud of the YSO. This ambient field can in principle be enough to trigger a magnetic breaking of the disc, leading to mass accretion and the formation of a wind. This idea was originally proposed in the context of YSOs by Wardle & Koenigl (1993). It was later revived through numerical simulations, first in local shearing box setups (Bai & Stone 2013; Simon et al. 2013; Lesur et al. 2014) and then in global simulations (Gressel et al. 2015, 2020; Béthune et al. 2017; Bai 2017; Wang et al. 2019). All of these simulations include Ohmic and ambipolar diffusion, while some of them also include the Hall effect, which is thought to be important in the densest regions of the disc (around 1 AU). The strength of these diffusions is in turn computed from ionisation and chemical models of various complexity. In addition, almost all of these models include some sort of heating of the disc atmosphere, designed to mimic non-thermal radiation heating that is expected in these systems (e.g. Thi et al. 2019). Each published model is of very high complexity, and only a handful of simulations can be performed to explore a (very) limited subspace of the parameter space.

In this context, the purpose of this work is to take a step back and simplify the physics in order to explore the coupling between protoplanetary discs and magnetised winds in a more or less systematic manner. This work essentially follows the approach initiated by Ferreira (1997) on globally self-similar models that was later extended to viscous (Casse & Ferreira 2000a) and weakly magnetised (Jacquemin-Ide et al. 2019) discs. I use a relatively different technic to derive the wind solutions, however, which I present in the next section, along with the physical model. I then focus on a fiducial series of wind solutions and explore their physical properties. Finally, I vary the disc diffusivities to offer an overview of the more exotic configurations before concluding.

2 Model

In the following, I consider a magnetised disc orbiting a central object of mass M. The disc is only subject to the gravitational pull of the central object (self-gravity being neglected) and to the Lorentz force due toelectrical currents. The disc is assumed to be weakly ionised, hence I consider a generalised Ohm law that includes Ohmic and ambipolar diffusivities. These diffusivities are prescribed in Sect. 2.3.

2.1 Equations

In the following, I use either spherical (r, θ, φ) or cylindrical (R, φ, z) coordinates, depending on the context. The cylindrical radius is written in upper case, and the spherical radius is lower case. I solve the non-ideal MHD equations, here in spherical coordinates, $${\partial }_t\rho +\nabla \cdot \rho v=0,$$(1) $${\partial }_t\rho v+\nabla \cdot \rho vv=-\nabla \rho c_{\text{s}}^2+\frac{J\times B}{c}-\frac{GM\rho }{r^2}{e}_r,$$(2) $${\partial }_{t}B=-\nabla \times \mathcal{E}$$(3) $$\mathcal{E}=-v\times B+\frac{4\pi }{c}({\eta }_{\text{O}}J+{\eta }_{\text{H}}J\times \widehat{B}-{\eta }_{\text{A}}J\times \widehat{B}\times \widehat{B}),$$(4)

where I havedefined the isothermal sound speed c_s(r), assuming the flow was locally isothermal, the plasma current J = c∇ × B∕4π, the gravitational constant G the magnetic field direction $\widehat{B}$, the electromotive field $\mathcal{E},$ and the Ohmic, Hall, and ambipolar diffusivities η_O, η_H, and η_A. The diffusivities are also functions of space and magnetic field strength. I stress that ‘no turbulence is assumed’ in this model, only molecular magnetic diffusivities that result from the low ionisation fraction of the plasma. If accretion occurs, it is only the result of torques that are self-consistently computed in the model.

2.2 Self-similar ansatz

In the context of protoplanetary discs, the strong diffusivities are known to almost suppress all non-axisymmetric structures (Béthune et al. 2017; Bai 2017) in the regions outside 1 AU. For this reason, I assume that the flow is 2.5D: I conserve three components for v and B but neglect φ derivatives.

In order to simplify the problem even more, I assume that the flow is globally self-similar. This implies that with increasing distance from the central object, the flow “looks” the same. This approach has several advantages: first it avoids issues with the inner boundary conditions, which are usually problematic in global numerical models. Second, it allows us to systematically explore the parameter space at exquisite resolution with limited numerical cost because the problem essentially becomes a 1D problem. I follow Ferreira & Pelletier (1993) and define the self-similar scaling for any field Q as $$Q(r,\theta )={(\frac{r}{r_0})}^{{\gamma }_Q}\tilde{Q}(\theta ),$$(5)

where γ_Q is the self-similar exponent, and $\tilde{Q}$ is a 1D function that completely determines the flow. The scaling of the gravitational force as 1∕r² imposes the self-similar scaling to all of the other components, that is, $${\gamma }_v=-\frac{1}{2};{\gamma }_B=-\frac{5}{4};{\gamma }_{\rho }=-\frac{3}{2};{\gamma }_{\eta }=\frac{1}{2}.$$

Similarly to the velocity field, the sound speed c_s is proportional to r^−1∕2, and $\widetilde{c_{\text{s}}}$ is a prescribed function of θ. The self-similar scaling is here written in spherical coordinates, but can be transformed into cylindrical (R, z, φ) coordinates by noting that r = Rsin(θ) and z∕R = tan⁻¹(θ), $$Q(R,z)={(\frac{R}{R_0})}^{{\gamma }_Q}\widehat{Q}(z/R),$$

where $\widehat{Q}(z/R)={\sin }^{{\gamma }_Q}(\theta )\tilde{Q}(\theta )$.

In the following, it will be useful to define the Keplerian velocity $v_{\text{K}}(R)\equiv \sqrt{GM/R}$ and the Keplerian angular frequency Ω_K ≡ v_K∕R. I also use the disc geometrical thickness h ≡ c_s(R, z = 0)∕Ω_K(R). By construction from the self-similar scaling, h∕R ≡ ε is a constant of the problem.

2.3 Diffusivities

The diffusivities η are usually computed from complex thermo-chemical networks (e.g. Thi et al. 2019). It is customary to use dimensionless numbers to quantify the diffusivities: the Ohmic Reynolds number Rm and the ambipolar Elsasser number Λ_A, defined as $$\text{Rm}\equiv \frac{{\Omega }_K{h}^2}{{\eta }_{\text{O}}}$$(6) $${\Lambda }_{\text{A}}\equiv \frac{V_{\text{A}}^2}{{\Omega }_K{\eta }_{\text{A}}}$$(7)

these dimensionless numbers are more useful than the more traditionally used Elsasser numbers because in most cases, they do not depend on the magnetic field strength (Wardle & Ng 1999). Hence, they only depend on the gas properties (density, temperature, and composition).

The use of the self-similar ansatz implies that Rm and Λ_A must be functions of z∕R (or θ) only. This is not true in full chemical models, except for ambipolar diffusion, which is known to be of the order of unity across a wide range of scales (see e.g. Thi et al. 2019, Fig. 8, middle panel). In the following, I therefore mostly focus on models dominated by ambipolar diffusion, and introduce Ohmic diffusion subsequently.

Following the Thi et al. (2019) models, I prescribe the ambipolar diffusion profile to be $${\Lambda }_{\text{A}}(z/R)={\Lambda }_{\text{A}}{}_0\exp [\frac{z^4}{{(\lambda h)}^4}],$$(8)

where λ is a free parameter quantifying the thickness of the non-ideal layer in units of h.

In models including Ohmic diffusion, I use the following diffusivity profile: $$\text{Rm}(z/R)={\text{Rm}}_0\exp (\frac{z^4}{{(\lambda h)}^4}){(\frac{\widehat{\rho }(z/h)}{\widehat{\rho }(0)})}^{-1}.$$

This profile,and in particular the dependence on the density ratio, is chosen to be consistent with the ambipolar diffusion profile, assuming the plasma is made of two types of charged species (such as electrons and molecular ions).

2.4 Numerical method

In contrast to commonly employed self-similar approaches (Casse & Ferreira 2000a; Jacquemin-Ide et al. 2019) where stationary equations are solved with shooting methods through the critical points of the flow, here I solved the time-dependent Eqs. (1)–(3) using the code PLUTO, a finite-volume, shock-capturing scheme (Mignone et al. 2007). I used a spherical domain, with the shape of a shell with only one grid point in the radial direction and 2048 points distributed homogeneously in the θ direction. I chose the grid to be centred on r = r₀ = 1, and the shell extended from θ = 0.15 to θ = π − 0.15. The radial boundary conditions were set enforcing the self-similar relations described in Sect. 2.2. For the boundary in the θ direction, I used standard outflow boundary conditions.

The initial condition was a disc in hydrostatic equilibrium threaded by a large-scale vertical (z) magnetic field whose initial value was set to have β = 10⁵ in the disc mid-plane. This initial condition is strongly unstable and the disc very quickly (in less than 10 Ω⁻¹) launched an outflow before reaching a quasi steady-state.

Given that the code is time-dependent and my choice of boundary conditions, the disc mass is not necessarily conserved. A fraction of the accreted material can be lost in the wind, leading to a slow decrease of the disc mass. In previous models, this has been taken into account by slightly adjusting the power exponent γ_ρ to ensure that mass is constant in the disc. This is not possible in our case because the solution is dynamically evolving: γ_ρ would need to be changed as a function of time. Instead, I therefore chose to renormalise the mass at each time step to keep the total mass inside the domain constant.

For similar reasons, the total magnetic flux threading the disc can evolve on secular timescales in this approach. This is not allowed in commonly used self-similar approaches, which usually assume that the total toroidal electromotive field is null. Here, I chose not toenforce such a constraint. This allowed me to measure the transport of magnetic flux, but it also implies that the magnetic flux evolves with time. As for the density, I therefore multiplied the each component of the field by a fixed factor at each time step, adjusted so that the magnetic flux threading the disc was equal to the desired value.

2.5 Accretion theory and diagnostics

Several diagnostics can be derived from self-similar solution. The most useful diagnostics are related to accretion theory, that is, how the disc surface density and mass accretion rate evolve with space and time. I therefore developed the mass, angular momentum, and magnetic flux conservation equations as $$\frac{\partial \Sigma }{\partial t}+\frac{1}{R}\frac{\partial }{\partial R}\frac{{\dot{M}}_{\text{acc}}}{2\pi }={[\overline{\rho v_z}]}_{-z_0}^{z_0},$$(9) $$\frac{{\dot{M}}_{\text{acc}}{v}_{\text{K}}}{4\pi }=\frac{\partial }{\partial R}\left(R^2{\displaystyle {\int }_{-z_0}^{z_0}\text{d}}z\overline{T_{R\phi }}\right)+R^2{[\overline{T_{z\phi }}]}_{-z_0}^{z_0},$$(10) $$\frac{\partial \overline{B_{z0}}}{\partial t}=-\frac{1}{R}\frac{\partial }{\partial R}R\overline{{\mathcal{E}}_{\phi 0}},$$(11)

where the overline denotes an azimuthal and ensemble average (in the case of a self-similar stationary solution, this average is not strictly needed), $\Sigma ={\displaystyle {\int }_{-z_0}^{z_0}\text{d}}z\overline{\rho }$ is the disc surface density, z₀ is the disc surface, ${\dot{M}}_{\text{acc}}\equiv -2\pi R{\displaystyle {\int }_{-z_0}^{z_0}\text{d}}z\overline{\rho v_R}$ is the mass accretion rate, T_xϕ = ρv_x(v_φ − v_K) − B_xB_φ∕4π is the stress φ component, and B_z0 = B_z(z = 0) and ${\mathcal{E}}_{\phi 0}$ are the vertical magnetic field strength and azimuthal EMF in the disc plane, respectively. It should be noted that in principle, accretion theory only require the first two relations. However, it is now well established that in a wind-driven disc, stresses and mass-loss rates are also a function of the mean poloidal field strength. In this context, accretion theory therefore has to be supplemented by Eq. (11).

Although Eqs. (9)–(11) fully describe the secular evolution of any disc at hand, it should be realised that the right hand-side terms are a priori unknown. This closure problem is well known in the disc community, and is usually solved by using the α disc paradigm when only the T_rφ term is present. In the case of a wind-driven disc, four terms are actually present. Following the α disc idea, I therefore define four parameters, $${\zeta }_{\pm }\equiv \pm \frac{\overline{\rho v_z}(\pm z_0)}{\Sigma {\Omega }_{\text{K}}},$$(12) $$\alpha \equiv \frac{{\displaystyle {\int }_{-z_0}^{z_0}\text{d}}z\overline{T_{R\phi }}}{{\displaystyle {\int }_{-z_0}^{z_0}\text{d}}z\overline{P}}=\frac{{\displaystyle {\int }_{z_0}^{z_0}\text{d}}z\overline{T_{R\phi }}}{\Sigma {\Omega }_{\text{K}}^2{h}^2},$$(13) $${\upsilon }_{\pm }\equiv \pm \frac{\overline{T_{z\phi }}(\pm z_0)}{\Sigma {\Omega }_{\text{K}}^{2}h},$$(14) $$v_{\text{B}}\equiv \frac{\overline{{\mathcal{E}}_{\phi 0}}}{{\Omega }_{\text{K}}h\overline{B_{z0}}},$$(15)

where I have assumed that the disc was isothermal to define α. With these definitions, a positive v_B implies that the field is ‘transported outwards’. The mass-loss rate ζ is comparable to the definition of Scepi et al. (2018), up to a factor of the order of unity. In numerical applications, I use z₀ = 6h.

It is also customary to measure the magnetic field strength as a function of the plasma β parameter. In this manuscript, I define the plasma β on the strength of the mean vertical field threading the disc mid-plane, $$\beta \equiv \frac{8\pi \overline{\rho }(z=0)c_{\text{s}}^2}{{\overline{B_{z0}}}^2}.$$(16)

The transport coefficients ζ, α, υ, and v_B being dimensionless coefficients, they are expected to depend only on dimensionless numbers. We therefore expect them to depend on the magnetic field strength β, but also on the diffusivity of the disc Rm and Λ_A. In principle, when all of these dependencies are known, accretion theory is said to be complete and the evolution of a wind-driven protoplanetary disc can be computed in a way similar to the historical α disc. Therefore, determining these dependencies is the main objective of this article.

Fig. 1 Space-time diagram showing the evolution for the first 1000 Ω−1 of Bθ (top), Bφ (middle), and vθ (bottom) for the fiducial run. The system quickly reaches steady state. In order to continue the solution as a function of the field strength, the field is slowly increased every 300 Ω−1.

3 Fiducial simulation

3.1 Time evolution

In the following, I define the fiducial run as a simulation with ambipolar diffusion only, Λ_A = 1 and λ = 3, so that the non-ideal MHD zone extends up to 3h above the disc mid-plane. The disc aspect ratio was fixed to ε = 0.1 and the disc and atmosphere were assumed to be locally isothermal, that is, c_s (R, Z) = c_s(R, z = 0). As for all the simulations, the fiducial run was started with a mid-plane β =10⁵. I have also run this simulation at half-resolution (i.e. 1024 points in θ) to confirm the convergence of the solution. The two solutions differ by less than 1%, so I am confident that the full resolution simulations presented here are perfectly resolved.

While most of the solutions presented in this article are quasi-steady state, I briefly discuss the space-time evolution of the simulations used to compute the solutions. This space-time evolution is shown in Fig. 1 for the fiducial run. The system clearly rapidly relaxes into steady state. During the first 100 Ω⁻¹, I applied a linear damping to the equations of motion to damp epicyclic oscillations more rapidly than the system would naturally do. These oscillations are due to the sudden launching of the wind, and are a spurious result of my wind-free initial conditions (we used a similar procedure in the shearing box Lesur et al. 2013).

After the first 100 Ω⁻¹, the system reaches a quasi-steady state, from which I can compute the required information: angular momentum transport, mass flux, etc. In order to continue the solution as a function of β, I slowly increased the total field strength between 300 Ω⁻¹ and 400 Ω⁻¹, and let the system reach its new equilibrium for the next 200 Ω⁻¹. I then repeated this procedure until the code stopped because of singular points in the solution. In doing so, I reached β =35 in the fiducial run in steps of 10% decrease in β. In the fiducial setup, I also continued the solutions from β = 10⁵ up to β = 10⁸ in steps of 10% increase to study the domain in which the solutions obtained by this procedure are valid. Hence, the domain explored for the fiducial setup is wider than for the other solutions (Sect. 4). In the following, I average the quasi-steady flow obtained for the last 100 Ω⁻¹ of “relaxation” for each magnetisation and use the result to compute the flow properties.

3.2 Flow topology

The 2D topology of the flow can be deduced from the 1D “shell” solutions of the simulation by reconstructing the 2D fields using the self-similar scalings (5). Using this procedure, I computed the flow topology for the two extreme cases (β = 10⁵ and β = 35) and for an intermediate case (β = 10³) shown in Fig. 2. In this fiducial simulation, the flow was symmetric with respect to the disc mid-plane, as seen here. The strongly magnetised solution produces a super-sonic accretion flow, while the weaker magnetised case tends to produce subsonic accretion (v_p∕c_s < 0.1). In the weak-field case, the field lines are almost straight, indicating that the decoupling between the field and the flow is efficient, thanks to ambipolar diffusion. As the field becomes stronger, however, the field lines are clearly pinched around the mid-plane. This indicates an efficient field dragging by the accretion flow.

The wind was super fast-magnetosonic in all cases. The fast surface was closer to the disc surface when the magnetisation became stronger, while the Alfvén surface tended to stay at the same altitude (around 4.5h). Generally speaking, stronger magnetisation leads to faster and more massive winds (for an identical disc). Interestingly, while the shape of the stream and field lines differ significantly in the disc, their aspect is very similar in the wind region starting above the disc surface. In all of these cases, the flow showed signs of collimation towards the axis.

3.3 Transport properties

The transport properties defined in Sect. 2.5 can be evaluated as a function of the magnetisation for the fiducial run. This gives the dependencies shown in Fig. 3. The resulting transport coefficients are remarkably close to power laws of β (except for v_B, which shows a shallower dependence on β than a power law). Several important trends can be deduced from this result.

First, it should be pointed out that because the solutions are all steady state and laminar, the α measured here is a ‘pure laminar stress’. It is not due to any form of turbulence because it is just the result of the large-scale wind that transports angular momentum both radially and vertically. Nevertheless, it can be compared to what is found in turbulent discs. The dependence of α on β⁻¹ is steeper than the one usually found in ideal-MHD shearing-box simulations, which usually give α ∝β^−1∕2 (e.g. Salvesen et al. 2016; Scepi et al. 2018). However, the ideal and non-ideal values for α match for β ~10. Hence, my values of α are generally lower than their ideal MHD counterpart by a ratio ~ β^1∕2∕3. This is mostlikely the signature of ambipolar diffusion in the disc.

In general, the values for α, υ, and ζ are compatible with the values found in the shearing-box literature for β ~10⁴−10⁵ and equivalent diffusivity regimes dominated by Λ_A = 1 (e.g. Bai 2013; Lesur et al. 2014), even though the dependence on β is slightly steeper than Bai (2013).

Interestingly, v_B > 0 for all the cases explored in the fiducial run. This indicates that despite the strong accretion in the disc mid-plane, the field is always diffusing outwards. This is qualitatively consistent with Bai & Stone (2017), but the diffusion speeds are reduced by an order of magnitude in my case. For β = 10⁴, I obtain v_B = 1.7 × 10⁻³, while Bai & Stone (2017) obtained¹ v_B = 4 × 10⁻². A careful examination of the ideal and non-ideal contributions to ${\mathcal{E}}_{\phi }$ (Fig. 4) shows that advection almost completely balances diffusion, diffusion winning only by a 1% excess. It is also worth noting that while the ideal and ambipolar EMFs vary strongly in the disc, their sum is approximately constant with θ, indicating that poloidal field lines are “moving” radially without much vertical deformation in the disc.

Using Eq. (10), the safe-similar scalings (Eq. (5)) and the definitions of Sect. 2.5, it is straightforward to show that the mass accretion rate is directly related to α and υ_± $$\frac{{\dot{M}}_{\text{acc}}}{2\pi \Sigma R^2{\Omega }_{\text{K}}}\equiv \dot{m}=\epsilon (\underset{{\dot{m}}_{\text{R}}}{\underbrace{\alpha \epsilon }}+\underset{{\dot{m}}_z}{\underbrace{2({\upsilon }_{+}+{\upsilon }_{-})}}),$$

where I have split Ṁ_acc in a contributionfrom the radial (ṁ_R) and vertical (ṁ_z) stresses. Even though these are usually referred to as “turbulent” and “wind” contributions (even when the flow is laminar), I emphasise that there is no turbulent stress. Hence both contributions are intrinsically due to the wind. Using the power-law scalings in Fig. 3, I find that $$\frac{{\dot{m}}_z}{{\dot{m}}_{\text{R}}}=1.5{\beta }^{0.26}.$$

Hence whenever β ≫ 1, the vertical stress is largely dominant and the radial transport of angular momentum can be neglected altogether, which is consistent with previous 3D numerical simulations (Béthune et al. 2017).

Neglecting ṁ_R, I can derive the mass accretion rate in a realistic disc using only the scaling for υ_±, which gives $$\begin{array}{ll}{\dot{M}}_{\text{acc}}=\hfill & 6.6\times {10}^{-8}\left(\frac{\Sigma }{10\text{g}{\text{cm}}^{-2}}\right){\left(\frac{R}{10\text{AU}}\right)}^{0.5}\hfill \\ \hfill & \times {\left(\frac{M}{M_{\odot }}\right)}^{0.5}{\left(\frac{\beta }{{10}^4}\right)}^{-0.78}{M}_{\odot }{\text{yr}}^{-1},\hfill \end{array}$$(17)

or substituting β with the field strength and the surface density, $$\begin{array}{ll}{\dot{M}}_{\text{acc}}=\hfill & 1.6\times {10}^{-8}{\left(\frac{\Sigma }{10\text{g}{\text{cm}}^{-2}}\right)}^{0.22}{\left(\frac{R}{10\text{AU}}\right)}^{2.08}{\left(\frac{M}{M_{\odot }}\right)}^{-0.28}\hfill \\ \hfill & \times {\left(\frac{\epsilon }{0.1}\right)}^{-0.78}{\left(\frac{B_z}{1\text{mG}}\right)}^{1.56}{M}_{\odot }{\text{yr}}^{-1}.\hfill \end{array}$$(18)

It is worth noting that the mass accretion rate does not depend strongly on the disc surface density, in contrast to the usual viscous disc model. Instead, I find that the field strength is the dominant control parameter in these solutions. More quantitatively, I find accretion rates compatible with expected rates from observations (~ 10⁻⁸–10⁻⁷ M_⊙ yr⁻¹) when β ∈ [10³, 10⁵] or B_z ~ a few mG. I caution that the observed accretion rates are typically inferred from the material that falls onto the central star, while here the mass accretion rate is computed in the disc at 10 AU. If the wind extracts a significant amount of material between the star and 10 AU, these accretion rates have no reason to match.

A common way to quantify the fraction of the accreted mass that is ejected is through the ejection efficiency parameter ξ, which is defined as $$\begin{array}{ll}\xi \hfill & \equiv \frac{2\pi R^2{[\overline{\rho v_z}]}_{-z_0}^{z_0}}{{\dot{M}}_{\text{acc}}}=\frac{\text{d}\log {\dot{M}}_{\text{acc}}}{\text{d}\log R}\text{insteadystate},\hfill \end{array}$$

where the last equality results from the continuity equation assuming ∂_t Σ = 0. It is easy to show that the ejection efficiency is directly connected to the transport coefficients: $$\begin{array}{ll}\xi \hfill & =\frac{{\zeta }_{+}+{\zeta }_{-}}{2\epsilon ({\upsilon }_{+}+{\upsilon }_{-})}\hfill \\ \hfill & =0.4{\beta }^{0.09},\hfill \end{array}$$

where the first equality assumes that ṁ_R is negligible,while the second uses the scaling found in the fiducial run (Fig. 3). This implies that ξ = O(1) for β ~10⁴ and that the accretion rate is approximately a linearly increasing function of radius in steady state in these models.

Finally, while most of the scalings presented here use β and therefore the vertical field strength as the main control parameter, these solutions are dominated by the azimuthal field component B_φ. Using the definition of υ, it easy to show that $$\begin{array}{ll}{\left|\frac{B_{\phi }}{B_z}\right|}_{z_0}\hfill & =\sqrt{\frac{\pi }{2}}\beta \upsilon \hfill \\ \hfill & \simeq 1.9{\beta }^{0.22},\hfill \end{array}$$(19)

where the last equality comes from the scaling of the fiducial run. Hence, for β = 10⁴, |B_φ | is more than 14 times larger than B_z at the disc surface, indicating that the field is strongly wrapped at this location. Finally, because υ₊ = −υ₋, the sign of B_ϕ is reversed at the top and bottom of the disc.

Fig. 2 Flow topology in the fiducial run (ambipolar diffusion only ΛA = 1). Top row: poloidal streamlines (white) and log of the sonic Mach number. Bottom row: poloidal field and log of density, normalised so that the mid-plane density at R = 1 is unity. The colour scales are identical in the columns. From left to right, the disc magnetisation increases: β = 105, 103, and 35. The green lines denotes critical lines of the flow: Alfvénic (plain) and fast magnetosonic (dot-dashed). The dashed green line represents the disc “surface” where the flow becomes ideal, which is arbitrarily located at z = 3.5h for all of the solutions.

Fig. 3 Transport coefficient as a function of the disc magnetisation for the fiducial run (ambipolar diffusion only ${\Lambda }_{\text{A}}{}_0=1$). Each blue dot is a steady-state solution obtained by continuation in the simulation. Because the flow is symmetric with respect to the mid-plane, only values computed from the disc top surface (ζ+ and υ+) are shown. The best fits using power laws are shown as a plain orange line.

Fig. 4 Zoom on the toroidal electromotive field of the disc, splitting the ideal MHD and ambipolar diffusion contributions for β = 105 in the fiducial case. The two electromotive fields balance almost exactly. The total EMF at θ = π∕2 is ${\mathcal{E}}_{\phi }/B_{z0}\simeq 6\times {10}^{-5}$ and is approximately constant throughout the disc. This gives vB = 6 × 10−4.

3.4 Magnetic torque and currents

Although it is possible to interpret accretion in terms of angular momentum budget (the disc angular momentum is transported into the wind), this gives us little information about how the disc accretes. This question can be addressed by considering the forces acting on the disc. In particular, because the outflow is magnetised, the Lorentz force F_L = J × B∕c is key. If F_L,φ < 0, the Lorentz force breaks down the disc rotation by creating a negative torque, driving accretion. The magnetic torque is clearly closely related to the poloidal current, $$F_{\text{L},\phi }=\frac{J_{\text{p}}\times B_{\text{p}}}{c}.$$

Because the poloidal field is more homogeneously distributed in the disc (Fig. 2), the magnetic torque is strongest in the regions of strong J_p. Hence the poloidal current is a direct indicator of the torque exerted by the wind on the disc. The distribution of poloidal currents in the fiducial β = 10⁵ case is shown in Fig. 5. I find a strongly focussed current sheet in the disc mid-plane, with currents directed outwards. This current induces a strong negative torque that in turn is responsible for accretion in the disc mid-plane. This current thencloses back in the disc surface and in the wind.

It might be surprising at first sight to find such a narrow current sheet in the system given that the disc is strongly diffusive (Λ_A = 1). It is well known, however, that ambipolar diffusion is a self-focusing diffusion process. It has been argued that in “free” reconnection sites, the current scales as z^−2∕3, where z is the distance from the null point (Zweibel & Brandenburg 1997). In the particular case presented here, I find that the current is even more focused with J_r ∝ z⁻¹ close to the mid-plane. This is expected because in the present case, I focus on a reconnection site that is “forced” by the outflow launched from the surface, hence the scaling need not be identical to that of Zweibel & Brandenburg (1997).

Fig. 5 Poloidal current in the fiducial simulation at β = 105. The current intensity is represented in colour. There is a strong outwards current in the disc mid-plane. The green lines represent the same characteristics as Fig. 2.

Fig. 6 Zoom on the radial current in the fiducial simulation at β = 105. The current sheet is strongly focused, with a scaling compatible with Jr ∝ z−1 (see text).

3.5 Wind properties

Winds in protoplanetary discs are often refereed to as magneto-thermal winds. This is because most (if not all) of the global simulations published so far used a “hot atmosphere”. These hot atmospheres are either achieved by a prescribed heating function (Béthune et al. 2017), or with a more complete treatment of the thermodynamics, including various sources of heating and cooling (Wang et al. 2019; Gressel et al. 2020). In any case, the energetics of these winds is partially dominated by thermal heating, hence their name.

In the fiducial solution presented here, the equation of state is locally isothermal. The temperature of the atmosphere therefore is that of the disc at the same cylindrical radius R. For this reason, my fiducial wind is not magneto-thermal, as I show below.

Because the wind is a steady-state ideal MHD flow, several conserved quantities can be derived from the equations of motion. I therefore considered a streamline passing through the mid-plane at R₀. I define Ω₀ as the Keplerian angular velocity at R₀: Ω₀ ≡ Ω_K(R₀). I then follow Blandford & Payne (1982), defining

• the mass-loading parameter

  $\kappa =\frac{4\pi \rho u_{\text{p}}{\Omega }_0{R}_0}{B_{\text{p}}{B}_0},$

  where v_p and B_p are the poloidal velocity and field strength, which quantifies the mass that is loaded onto the field lines.

• the rotation parameter

  $\omega \equiv \frac{\Omega }{{\Omega }_0}-\frac{\kappa B_0{B}_{\phi }}{4\pi \rho R{R}_0{\Omega }_0^2,}$

  which measures the rotation speed of magnetic surfaces.

• the magnetic lever arm

  $\lambda \equiv \frac{\Omega R^2}{{\Omega }_0{R}_0^2}-\frac{R{B}_{\phi }}{R_0\kappa B_0,}$

  which measures the angular momentum that is extracted by the wind.

• the Bernoulli invariant

  $\mathcal{B}=\underset{\text{Kinetic}}{\underbrace{\frac{v^2}{2{\Omega }_O^2{R}_0^2}}}\underset{\text{Gravitational}}{\underbrace{-\frac{R_0}{\sqrt{R^2+z^2}}}}\underset{\text{Magnetic}}{\underbrace{-\frac{\omega R{B}_{\phi }}{\kappa R_0{B}_0}}}\underset{\text{Thermal}}{\underbrace{+\frac{w}{{\Omega }_0^2{R}_0^2}}}$

  which measures the energy content of the flow. In this definition, I included the thermal contribution to the flow energetics as the work done by thermal pressure forces along the field line,

  $w\equiv {\displaystyle {\int }_s^{\infty }-}\text{d}\mathcal{l}\cdot \frac{\nabla P}{\rho }.$

  If the flow is adiabatic, this is simply the enthalpy. However, the locally isothermal approximation does not lead to an adiabatic flow, hence I have to use the integral form.

I computed the mass loading, rotation, and lever arm parameters along a magnetic field line for the fiducial simulation at β = 10⁵ (Fig. 7). As expected, these parameters are constant once above the non-ideal region of the disc. Surprisingly, the lever arm is very small, reaching only λ = 1.56. As is well known, cold MHD winds require λ > 3∕2, so this is just above this limit value. A close inspection of the Bernoulli invariant (Fig. 8) shows that thermal driving is indeed negligible in the energy budget of the outflow, confirming that the solution is a cold MHD wind, not a magneto-thermal wind. It also shows that $\mathcal{B}>0,$ indicating that the wind is free to escape the gravitational potential well up to infinity.

A systematic exploration of the κ, λ parameters as a function of β in the fiducial simulation gives Fig. 9. This figure can directly be compared to Fig. 2 in Blandford & Payne (1982). It is interesting to note that these new solutions at low magnetisation have a much larger κ and much smaller λ than the previous solutions of Blandford & Payne (1982). This trend has been observed by Jacquemin-Ide et al. (2019) in ideal MHD for β >1. The solutions I obtained here have even smaller λ, with an asymptote in the high β limit close to λ = 1.4. This is smaller than the λ = 3∕2 limit demonstrated by Blandford & Payne (1982). It should be pointed out, however, that this limit is computed for a wind emitted from z = 0, while here the wind is emitted from z = 5H, so the lever arm is physically allowed to be slightly smaller than 3∕2. This is confirmed by analysingthe Bernoulli invariant at β = 10⁷, which still shows that the magnetic energy dominates the energetics at the disc surface (not shown).

Fig. 7 MHD invariants plotted along one given field line for the fiducial run at β = 105. As expected, the invariants are approximately constant once in the ideal-MHD region lies above the disc. The asymptotic values are κ = 30.2, ω = 0.99, and λ = 1.56.

Fig. 8 Bernoulli invariant plotted along one given field line for the fiducial run at β = 105. The Bernoulli invariant becomes constant once in the ideal MHD region, as expected. The thermal contribution to the energetics isnegligible compared to the magnetic contribution, hence the wind is classified as a “cold” MHD wind.

Fig. 9 Dependence on κ and λ of the fiducial solutions as a function of the disc magnetisation β. As the field strength increases, λ increases, while the mass-loading parameter decreases.

4 Parameter space exploration

4.1 Ambipolar diffusion

In order to explore the effect of the disc diffusivity on the behaviour of the outflow, I varied the strength of ambipolar diffusion in the disc plane ${\Lambda }_{\text{A}}{}_0$ (Eq. (8)), keeping the same vertical profile. The resulting transport properties are shown in Fig. 10. The general trend is that as diffusion is increased, the mass-loss rate and angular momentum fluxes decrease, while magnetic flux transport increases. This overall suggests that wind-driven accretion tends to be reduced when ambipolar diffusion increases, as expected naively. It should be noted, however, that not all transport coefficients vary in the same proportion. Taking β =10⁴ as the reference case, I find that a decrease in ${\Lambda }_{\text{A}}{}_0$ by a factor 16 leads to adecrease in ζ, α, and υ by a factor 3.2, 4.0, and 2.25, respectively, while v_ψ increases by a factor 4.0. This therefore suggests a relatively weak sensitivity of the transport coefficient on ${\Lambda }_{\text{A}}{}_0$ (shallower than ${\Lambda }_{\text{A}}^{0.5}$), especially for the vertical angular momentum coefficient υ. The weak dependence of υ on the disc dissipation properties has been observed in shearing-box simulations (see Lesur 2020, Fig. 39 and related text) and is here confirmed in global geometry. Overall, this confirms that wind-driven angular momentum extraction is only moderately affected by the strength of ambipolar diffusion, and that the mass accretion rate (17) is approximately valid (up to a factor of a few) for the range of ambipolar diffusion considered here.

4.2 Ohmic diffusion

In addition to the strength of ambipolar diffusion, I also explored the effect of Ohmic diffusion on the flow topology. In protoplanetary discs, this leads to wind solutions that are valid in the Gammie (1996) historical “dead zone”, found around 1 AU in typical chemical models. To limit the computation time required by such models, I chose to focus on Rm₀ > 1. In the literature, Ohmic diffusion is frequently evaluated in terms of the Ohmic Elsasser number Λ_O, which depends on the field strength. Typically, Λ₀ = Rm₀β, so that Rm₀ ≳ 1 typically corresponds to the mid-plane value quotedat R ≲ 1 AU by several authors (Bai & Stone 2013, Fig. 1; Thi et al. 2019, Fig. 8). I finally point out that in all of these models, I kept ambipolar diffusion fixed to its fiducial value ${\Lambda }_{\text{A}}{}_0=1$. This means that Ohmic diffusion is added to the diffusivity tensor of the fiducial run.

As for ambipolar diffusion, I started with the transport coefficients as a function of Rm₀ (Fig. 11). Here, the angular momentum transport coefficients α and υ are barely affected by Ohmic diffusion, even when Rm₀ = 1. This suggests that the accretion rate in the Gammie (1996) dead zone probably is not so different from that observed in regions dominated by ambipolar diffusion. I note, however, that the field diffusion rate is increased by possibly several orders of magnitude. This is probably the most stringent effect of adding Ohmic diffusion to the system. While Rm₀ ≥ 10 solutions are top-down symmetric (therefore only ζ₊ and υ₊ are plotted), the Rm₀ = 1 solution breaks this top-down symmetry around β ≃ 10⁴. The direction of this symmetry breaking is random. Starting from slightly different initial condition, I could obtain ζ₊ > ζ₋ or the opposite, which shows that there is ‘no physically preferred direction’ for this symmetry breaking.

The fact that the transport coefficients are approximately unaffected by Ohmic diffusion should not lead to the incorrect impression that the solutions are physically identical. A careful examination of the solution at Rm₀ = 1 (Fig. 12) shows that the solution in the disc is very different from the fiducial solution. First, the current sheet that was initially localised in the mid-plane (Fig. 5) is now localised on the disc surfaces. This can be easily interpreted as a result of Ohmic diffusion. While ambipolar diffusion can have a focusing effect, Ohmic diffusion is really only a diffusive process. The current sheet therefore cannot form in the strongly diffusive mid-plane, and is therefore pushed away, at the disc surface.This change in localisation of the current sheet in turn affects the accretion flow because accretion is driven by the magnetic torque created by the poloidal current (see Sect. 3.4). This is more evidently seen in Fig. 13, where the accretion flow clearly moves from the disc mid-plane towards the disc surface as Rm₀ is decreased, thereby following the localisation of the current layer. Finally, poloidal current and accretion flows are both localised at the disc surface. However, because the current is set globally by the electrical circuit enforced by the wind, its intensity is not really affected by this change in localisation, implying that total torques and accretion rates are not dramatically modified, as indicated bytransport coefficients.

After understanding that the current sheet becomes divided into two parts and pushed towards the surface, this process might be imagined to be not be totally symmetric for some range of parameters. This is what happens in the Rm₀ = 1 solution around β = 10⁴. In this case, the current sheet is more pronounced in the southern hemisphere (Fig. 14), leading to an accretion stream that is morepronounced on the bottom side of the disc, and inclined field lines in the disc bulk. I stress that while accretion occurs in the southern side of the disc in this particular example, the wind is more pronounced on the northern side of the disc because |ζ₊ |> |ζ₋| and υ₊ > υ₋. Hence angular momentum and mass flow away on the side opposite to that of accretion. This behaviour was also observed by Béthune et al. (2017) in dissymmetric solutions.

5 Discussion and conclusions

I have used a finite-volume code to derive self-similar wind solutions applicable to the diffusive regime of protoplanetary discs. In contrast to previous self-similar solutions, no assumption was made on flux stationarity or top-down symmetry. I also circumvented the difficulties of full-blown 3D numerical simulations, which are often limited by their integration time and inner boundary conditions. Using this technic, it is possible to explore systematically a wide range of parameters at a reduced numerical cost.

Using this tool, I have presented a series of wind solutions ranging from β = 10⁸ to β = 35 that are valid in the non-ideal regions of protoplanetary discs R ≳ 1 AU. I showed that cold wind solutions (i.e. that do not require atmospheric heating) exist in this entire range of parameters. Some of these solutions are very similar to the solution found using 2.5D and 3D simulations, in the same range of parameters. These solutions exhibit several important properties, some of which can be confirmed observationally.

First and foremost, magnetised disc winds always exist and are unavoidable as soon as a large-scale magnetic field threads the disc. This statement is true even for β = 10⁸ fields that correspond to a few μG at 10 AU for typical surface densities. They give accretion rates Ṁ_acc (Eq. (18)), which essentially depend on the field strength, and which are broadly speaking compatible with observed accretion rates onto T Tauri stars provided that 10³ < β < 10⁵ or equivalently B_z ~ a few mG. Interestingly, Ṁ_acc depends only weakly on the disc surface density, in contrast to viscous disc models. However, the direct comparison between the accretion rate in the disc bulk (as measured by Eq. (18)) and the mass accretion measured onto the star is probably misleading. The latter is probably significantly smaller than the former because the ejection index (see below) is about one, hence the disc mass accretion rates could be significantly higher (possibly by an order of magnitude) than the rates quoted here.

The ejection index ξ, which essentially quantifies the ratio of ejected to accreted mass, depends only weakly on the magnetic field strength and is of the order of (but slightly lower than) unity for the range of magnetisation I explored. This means that the mass loss rate is comparable to the mass accretion rate at a given radius. Consequently, when steady state is assumed, the accretion rate is expected to increase significantly with radius. As expected from a high ξ wind (Lesur 2020, Eq. (11.22)), the wind lever arm is always relatively small, with typically λ < 2 for β >10⁴, without requiring any heating. Such a high ξ and low λ has been found in previous numerical work (Béthune et al. 2017; Bai 2017), and it was initially thought that atmospheric heating, included in all these models, was the main reason for this result, following the argument of Casse & Ferreira (2000b). I showed here that heating is not key to obtain these high ξ-low λ solutions, a result that was also recently reported by Jacquemin-Ide et al. (2019) in the context of fully ionised discs. The solutions presented here are locally isothermal, and I showed that thermal heating is negligible in the wind energy budget even for β =10⁸. This means that low λ and high ξ are quite clearly signatures of weakly magnetised (β ≫ 1) outflows.

The vertical field strength predicted to derive accretion rates compatible to typical T Tauri stars (about a mG at 10 AU, see Eq. (18)) are compatible with recent upper limits obtained from Zeeman measurements (e.g. Vlemmings et al. 2019). For field detection, it is quite clear that the toroidal is a better candidate than the vertical component because the former is expected to be much stronger (by a factor 10 or so at β ~10⁴) than the latter (Eq. (19)). The toroidal component is expected to change sign across the mid-plane, however, so that this detection is possible only for tracers that do not average out the field through the disc thickness.

For the disc microphysics, I find that the disc-averaged mass accretion and ejection rates depend only weakly on ambipolar and Ohmic diffusivities. This implies that very elaborated ionisation models in the disc are not required to obtain the right accretion or ejection bulk properties. However, the disc vertical structure does strongly depend on diffusion. Solutions with large Ohmic diffusion (i.e. valid closer to 1 AU) tend to exhibit accretion at the disc surface, while weaker Ohmic diffusion (R > 10 AU) shows accretion in the disc mid-plane. Some of the strong Ohmic diffusion solutions also exhibit top-down dissymmetry, where accretion mostly occurs on one side of the disc. An asymmetry like this was also found in simulations including Ohmic and ambipolar diffusion (Béthune et al. 2017; Gressel et al. 2020) in a similar range of parameters, so that this result has been reported before. These topology properties will likely have a strong effect on dust growth and planet migration theory, but again, they barely affect the vertically averaged transport properties.

Finally, in all of these models, the large-scale magnetic field is found to be transported outwards, in agreement with Bai & Stone (2017). However, the transport velocity is about an order of magnitude lower than that measured by Bai & Stone (2017) in similar conditions (ambipolar diffusion only). Moreover, I find that the transport rate varies significantly with the strength of Ohmic diffusion and to a lesser extent with ambipolar diffusion. Overall, this dependence might explain why simulations and models tend to disagree on the magnetic transport rate (Bai & Stone 2017; Leung & Ogilvie 2019; Gressel et al. 2020).

Clearly, the large-scale field transport is the main bottleneck of a complete theory of wind-driven accretion because the strength of the large-scale poloidal field is the main control parameter. This cannot be ignored for the secular evolution of discs because the field strength is expected to vary with time. While class II objects seem to require field strengths of a few mG at 10 s of AU to obtain the correct accretion rates (see above), core-collapse calculations tend to suggest a field strength of 100 mG when the disc forms (e.g. Masson et al. 2016). This means that the field strength must be reduced by at least one order of magnitude (more probably two) between class 0 and class 2, indicating that flux transport must be relativelyefficient. The problem of field transport therefore deserves more attention in the future if winds are to be used in secular evolution models.

Fig. 10 Transport coefficients as a function of the disc magnetisation and mid-plane Elsasser number ${\Lambda }_{\text{A}}{}_0$.

Fig. 11 Transport coefficients as a function of the disc magnetisation and mid-plane Reynolds number Rm0. These solutions also include ambipolar diffusion with ${\Lambda }_{\text{A}}{}_0=1$. The fiducial run is physically equivalent to Rm0 = ∞.

Fig. 12 Flow structure for the solution Rm0 = 1 at β = 105. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. Compared to the fiducial case, the current layer is now localised at the disc surface, leading to a surface accretion flow.

Fig. 13 Inward mass flux (top) and radial velocity Mach number (bottom) as a function of Ohmic diffusion strength at β = 105. The accretion flow moves towards the disc surface as Ohmic diffusion is increased, following the poloidal current sheet (see text).

Fig. 14 Flow structure for the solution Rm0 = 1 at β = 104. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. This solution exemplifies the dissymmetric solutions found around β = 104 at Rm0 = 1.

Acknowledgements

I am thankful to Jonatan Jacquemin-Ide, Etienne Martel, Jonathan Ferreira and Cornelis Dullemond for fruitful discussions regarding the physics of self-similar solutions. I acknowledge support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme (Grant agreement No. 815559 (MHDiscs)). All of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.

References

All Figures

	Fig. 1 Space-time diagram showing the evolution for the first 1000 Ω−1 of Bθ (top), Bφ (middle), and vθ (bottom) for the fiducial run. The system quickly reaches steady state. In order to continue the solution as a function of the field strength, the field is slowly increased every 300 Ω−1.
In the text

Fig. 2 Flow topology in the fiducial run (ambipolar diffusion only ΛA = 1). Top row: poloidal streamlines (white) and log of the sonic Mach number. Bottom row: poloidal field and log of density, normalised so that the mid-plane density at R = 1 is unity. The colour scales are identical in the columns. From left to right, the disc magnetisation increases: β = 105, 103, and 35. The green lines denotes critical lines of the flow: Alfvénic (plain) and fast magnetosonic (dot-dashed). The dashed green line represents the disc “surface” where the flow becomes ideal, which is arbitrarily located at z = 3.5h for all of the solutions.

In the text

Fig. 3 Transport coefficient as a function of the disc magnetisation for the fiducial run (ambipolar diffusion only ${\Lambda }_{\text{A}}{}_0=1$). Each blue dot is a steady-state solution obtained by continuation in the simulation. Because the flow is symmetric with respect to the mid-plane, only values computed from the disc top surface (ζ+ and υ+) are shown. The best fits using power laws are shown as a plain orange line.

In the text

	Fig. 4 Zoom on the toroidal electromotive field of the disc, splitting the ideal MHD and ambipolar diffusion contributions for β = 105 in the fiducial case. The two electromotive fields balance almost exactly. The total EMF at θ = π∕2 is ${\mathcal{E}}_{\phi }/B_{z0}\simeq 6\times {10}^{-5}$ and is approximately constant throughout the disc. This gives vB = 6 × 10−4.
In the text

	Fig. 5 Poloidal current in the fiducial simulation at β = 105. The current intensity is represented in colour. There is a strong outwards current in the disc mid-plane. The green lines represent the same characteristics as Fig. 2.
In the text

	Fig. 6 Zoom on the radial current in the fiducial simulation at β = 105. The current sheet is strongly focused, with a scaling compatible with Jr ∝ z−1 (see text).
In the text

	Fig. 7 MHD invariants plotted along one given field line for the fiducial run at β = 105. As expected, the invariants are approximately constant once in the ideal-MHD region lies above the disc. The asymptotic values are κ = 30.2, ω = 0.99, and λ = 1.56.
In the text

	Fig. 8 Bernoulli invariant plotted along one given field line for the fiducial run at β = 105. The Bernoulli invariant becomes constant once in the ideal MHD region, as expected. The thermal contribution to the energetics isnegligible compared to the magnetic contribution, hence the wind is classified as a “cold” MHD wind.
In the text

	Fig. 9 Dependence on κ and λ of the fiducial solutions as a function of the disc magnetisation β. As the field strength increases, λ increases, while the mass-loading parameter decreases.
In the text

	Fig. 10 Transport coefficients as a function of the disc magnetisation and mid-plane Elsasser number ${\Lambda }_{\text{A}}{}_0$.
In the text

	Fig. 11 Transport coefficients as a function of the disc magnetisation and mid-plane Reynolds number Rm0. These solutions also include ambipolar diffusion with ${\Lambda }_{\text{A}}{}_0=1$. The fiducial run is physically equivalent to Rm0 = ∞.
In the text

	Fig. 12 Flow structure for the solution Rm0 = 1 at β = 105. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. Compared to the fiducial case, the current layer is now localised at the disc surface, leading to a surface accretion flow.
In the text

	Fig. 13 Inward mass flux (top) and radial velocity Mach number (bottom) as a function of Ohmic diffusion strength at β = 105. The accretion flow moves towards the disc surface as Ohmic diffusion is increased, following the poloidal current sheet (see text).
In the text

	Fig. 14 Flow structure for the solution Rm0 = 1 at β = 104. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. This solution exemplifies the dissymmetric solutions found around β = 104 at Rm0 = 1.
In the text