© ESO, 2013

1. Introduction

Accretion discs are found around several kinds of astrophysical objects: from young stars (protoplanetary discs) to supermassive black holes in active galactic nuclei. The dynamics of these discs is however poorly understood. It is known that, on average, matter moves inward resulting in the accretion of material onto the central object. Dynamically, one can show easily (Lynden-Bell & Pringle 1974) that this accretion is controlled by the loss of angular momentum of the falling gas. Hence, if one is to predict the accretion rate and survival time of discs, one needs to understand the angular momentum dynamics of these objects.

For accretion to happen on timescales compatible with observations, angular momentum must be transported rather efficiently. This has led to the α disc model (Shakura & Sunyaev 1973) in which angular momentum is transported radially outward in the disc by a prescribed “turbulent viscosity”. The origin of this turbulent viscosity has been debated for several decades and remains even today rather unclear. The magnetorotational instability (MRI: Balbus & Hawley 1991; Hawley et al. 1995) appears as the most promising way of producing turbulence in discs, though other processes might also be at work in some regions (see Armitage 2011, for a review). However, transporting angular momentum in discs is not the only way to cause accretion. Instead, one can suppose that angular momentum is “extracted” from a disc by a large-scale magnetic structure anchored in the disc midplane (Blandford & Payne 1982; Pudritz & Norman 1983). This sort of mechanism usually produces magnetically driven outflows, also known as disc winds.

It is most likely that these two mechanisms (angular momentum transport through turbulence and extraction by winds) actually coexist in astrophysical objects. However, this link is, from a theoretical point of view, poorly understood. The first reason is that stationary disc wind models require a strong magnetisation (plasma β ~ 1, Ferreira & Pelletier 1995) which quenches the MRI by magnetic tension effects (Balbus & Hawley 1991). At first sight, the conditions of existence for MRI turbulence and disc winds are therefore mutually exclusive. Moreover, a numerical computation including both a large-scale wind and small-scale turbulence in the disc midplane is technically very challenging: resolving turbulence in the bulk of the disc requires a large resolution in the disc midplane and fully 3D simulations, whereas computing the wind itself requires a very large computational domain. Several authors have therefore computed disc winds using a prescription to take into account the effects of disc turbulence (essentially a turbulent viscosity and resistivity), both analytically (Casse & Ferreira 2000) and numerically (Zanni et al. 2007). On the MRI side, most of our knowledge comes from shearing box simulations (Hawley et al. 1995; Stone et al. 1996; Longaretti & Lesur 2010), although several authors have also considered global configurations (Hawley 2000; Fromang et al. 2011; Flock et al. 2011). These global configurations are however limited to situations without any poloidal magnetic structure (or very weak mean poloidal magnetic fields) which precludes the production of outflows.

More recently the production of outflows from MRI turbulence in shearing boxes has been studied in the limit of weak poloidal fields, β = 10⁵ − 10⁴ (Suzuki & Inutsuka 2009). A magnetocentrifugal mechanism similar to Blandford & Payne (1982) has been identified in these simulations, with non-steady outflows starting about two scale heights above the midplane, both in the weak field regime β = 10⁴ (Fromang et al. 2013) and in the moderate field regime β = 10²–10³ (Bai & Stone 2013). However, the outflow mass loss rate was shown to depend strongly on the box vertical size and aspect ratio. Moreover, the outflow for such a weak mean field is dynamically insignificant for the disc as it extracts a very small amount of angular momentum. In the strong poloidal field regime (β ≲ 1) which is stable to MRI modes, Ogilvie (2012) studied the production of quasi 1D steady outflows in shearing boxes. This study demonstrated that some properties of global solutions were recovered in shearing boxes, although the procedure used in this work was not used to look in the parameter regime unstable to the MRI. More recently, a related set of simulations in 2D (axisymmetric) shearing boxes was presented by Moll (2012) who founds that the wind became striped and unsteady due to an unknown instability.

The aim of this paper is to study the potential link between the MRI and quasi-steady outflows which are found in global models and simulations. To this end, we study stratified shearing boxes threaded by a strong poloidal field (β ~ 10) which are not far from the steady solutions of Ogilvie (2012) but are lying in the MRI unstable regime. We first present the shearing box model and the numerical methods we have used in this investigation. We then look at the saturation of 1D MRI modes which naturally produce outflows. These “MRI-outflows” are compared to other solutions found in the literature, both local and global. We then study the 3D stability of these outflows and demonstrate that they are unstable on dynamical timescales. The consequences of this instability and its nature are briefly discussed. Finally, we summarise our findings and discuss their potential implications for astrophysical outflows and jets.

2. Local model

2.1. Equations

MRI-related turbulence and the shearing box approximation have been extensively described in the literature. However, since the physics we are looking for involves outflows and mass losses, we recall here briefly the basic equations for the shearing box model. Interested readers may also consult Hawley et al. (1995), Balbus (2003) and Regev & Umurhan (2008) for an extensive discussion of this approximation.

The shearing-box equations are found by considering a cartesian box centred at a fiducial radius R₀, rotating with the disc at constant angular velocity Ω = Ω_K(R₀) and having dimensions (L_x,L_y,L_z) with L_i ≪ R₀. We define r − R₀ → x, R₀φ → y similarly to Hawley et al. (1995). Furthermore, we assume the disc follows an isothermal equation of state with a constant sound speed c. In this approximation, the MHD equations read: $\begin{array}{ccc}& & \\ & & \\ & \mathrm{-}\mathrm{2}\mathit{\rho }\mathbf{\Omega }\mathbf{\times }{u}\mathrm{-}\mathit{\rho }\mathbf{\nabla }\mathit{\psi ,}& \\ & & \end{array}$where we have chosen the units so that μ₀ = 1 and ψ is a local expansion of the effective gravitational potential around R₀: $\mathit{\psi }\mathrm{=}\mathrm{-}\mathit{q}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{z}}^{\mathrm{2}}$(4)in which we have considered a Keplerian disc having a rotation profile Ω_K(r) ∝ r ^− q with q = 3/2.

One should note that the above set of equations admits a solution as a linear shear flow which is a local approximation to the Keplerian profile, u =  −qΩxe_y. In the following, we will consider perturbations v (not necessarily small) to this Keplerian profile so that u = v − qΩxe_y.

We consider a shearing box of size (L_x,L_y) = (16,16), the unit of length being defined by the disc scale height: H = c/Ω. The time unit is Ω^-1 and the velocity unit is c. We will assume ρ = 1 in the disc midplane at the hydrostatic equilibrium which sets our unit of mass. As shown below, we only consider the upper half of the disc, so that z ∈  [0,z_B]  where z_B is the altitude of the upper z boundary condition. Unless otherwise stated, we assume z_B = 6. Because the total vertical magnetic flux is conserved, we introduce the dimensionless magnetisation $\mathit{\mu }\mathrm{=}\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{\Sigma \Omega }\mathit{c}}\mathit{,}$(5)as a control parameter¹ of our simulations, where Σ = ^∫dV   ρ/(L_xL_y) is the equivalent surface density of the box.

In this study, we only consider the ideal MHD approximation presented above, in which we neglect all non-ideal effects such as viscosity and Ohmic resistivity. It is known that these non-ideal effects have a strong impact on the saturation level of MRI turbulence (Lesur & Longaretti 2007; Fromang et al. 2007; Longaretti & Lesur 2010; Simon et al. 2011; Fromang et al. 2013). However, much of this work is dedicated to the study of quasi-stationary 2D smooth outflows, for which viscosity and resistivity effects should be negligible. Some physical processes not discussed here, such as the saturation of 3D instabilities presented in Sect. 4, might depend on these non-ideal effects. These processes will be the subject of a separate paper.

2.2. Numerical model

2.2.1. Numerical method

To investigate the above system of equations, we use the PLUTO code (Mignone et al. 2007). PLUTO uses a finite volume method with a Godunov scheme to integrate the equations in their conservative form. MHD terms are computed using the constrained transport method of Evans & Hawley (1988) which enforces ∇·B = 0 at machine precision during the evolution of the physical system. In this procedure, the electromotive forces are reconstructed using the UCT_CONTACT scheme (Gardiner & Stone 2005) which have been shown to accurately reproduce analytical MRI growth rates in PLUTO (Flock et al. 2010). We use the Roe method to solve the Riemann problem at cell boundaries. This choice was dictated by the presence of strong numerical instabilities with the HLLD solver when the plasma beta β = 2ρc²/B² was too small (typically β < 1). Moreover, in order to stabilise the code, we switch to an HLL solver (more diffusive) and a minmod slope limiter whenever the magnetisation is very large (typically β < 10^-4) in 3D runs. Only a few cells have such a small β and the precise value of the threshold does not significantly change the outcome of the simulations. This however prevents the code from failing when the Alfvén speed becomes too large. Similar techniques have been used in numerical studies of supersonic interstellar turbulence (Lemaster & Stone 2009). All the simulations discussed in this paper are summarised in Table 1.

2.2.2. Boundary conditions

To reduce the computational costs of the simulations, we compute only the upper half of the disc. The lower half is deduced by symmetry: ρ( − z) = ρ(z); v_H( − z) = v_H(z); v_z(z) =  −v_z(z); B_H( − z) =  −B_H(z) and B_z( − z) = B_z(z) where H stands for the horizontal (x,y) component of a vector. It should be noted that this symmetry is a natural symmetry of the equations of motion. This implies that if the initial conditions satisfy this symmetry, the resulting solution will verify this symmetry at all times.

The boundary conditions we impose are shear-periodic in the radial direction and periodic in the azimuthal direction. The midplane symmetry described above is imposed at z = 0. The upper boundary condition (z = z_B) is the most delicate part of the setup. Unless otherwise stated, we consider modified outflow boundary conditions where we enforce a strictly vertical poloidal field at the boundary: $\begin{array}{ccc}& & {\mathit{\partial }}_{\mathit{z}}\mathit{\rho }\mathrm{\left(}{\mathit{z}}_{\mathrm{B}}\mathrm{\right)}\mathrm{=}{\mathit{\partial }}_{\mathit{z}}{{v}}_{\mathrm{H}}\mathrm{\left(}{\mathit{z}}_{\mathrm{B}}\mathrm{\right)}\mathrm{=}{\mathit{\partial }}_{\mathit{z}}{\mathit{B}}_{\mathit{y}}\mathrm{\left(}{\mathit{z}}_{\mathrm{B}}\mathrm{\right)}\mathrm{=}\mathrm{0}\\ & & {\mathit{B}}_{\mathit{x}}\mathrm{\left(}{\mathit{z}}_{\mathrm{B}}\mathrm{\right)}\mathrm{=}\mathrm{0.}\end{array}$Surprisingly, strict outflow boundary conditions (zero gradient for all fields) prevent MHD driven winds to be produced. This point is discussed more extensively in Sect. 3.4. Boundary conditions used for each run are specified in Table 1.

Because many of the simulations described here show an MHD-driven wind, a significant amount of mass is lost in our model. To mimic the mass inflow due to the material which would be accreted in a realistic disc including curvature effects, we have chosen to add mass in the midplane to compensate for wind losses. This is accomplished by adding mass at z ≲ z_inj at each numerical time step, modifying the mass conservation so that: ${\mathit{\partial }}_{\mathit{t}}\mathit{\rho }\mathrm{+}\mathbf{\nabla }\mathbf{·}\mathit{\rho }{v}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\exp \left[\mathrm{-}(\frac{\mathit{z}}{\mathrm{2}{\mathit{z}}_{\mathrm{inj}}}{)}^{\mathrm{2}}\right]\mathit{,}$(8)where ρ₀ is adjusted so that ^#ρ   dV is constant in time. This mass addition can be done in several ways. We have chosen to keep the velocity field of the gas constant while adding mass. This obviously implies that conservation of momentum and energy are broken in the injection region z ≲ z_inj. However, this choice has a little impact on the numerical results presented here: using a momentum conserving mass injection procedure did not modify significantly the outcome of the run 1DRef. In all the simulations discussed below have z_inj = 0.1, unless otherwise stated. Tests with z_inj = 0.05 have shown that none of the results we discuss thereafter are significantly affected.

2.2.3. Modified potential

The vertical hydrostatic equilibrium described by (2) leads to a gaussian vertical density profile: $\mathit{\rho }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\exp \left[\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}^{\mathrm{2}}}\right]\mathrm{·}$(9)Assuming a constant B_z in the box we get an Alfvén speed ${\mathit{V}}_{\mathrm{A}}\mathrm{=}{\mathit{B}}_{\mathit{z}}\mathit{/}\sqrt{\mathit{\rho }}\mathrm{\propto }\exp \mathrm{(}{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{4}{\mathit{H}}^{\mathrm{2}}\mathrm{)}$ which increases very steeply as z increases. Because of the CFL condition, this leads to very small timesteps which dramatically increase the computational time. To reduce the computational costs, several of our simulations were performed using the modified potential: ${\mathit{\psi }}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathit{q}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{[}\mathrm{1}\mathrm{-}\exp {\mathrm{(}\mathrm{-}{\mathit{z}}^{\mathrm{2}}\mathit{/}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{]}}\mathrm{·}$(10)This modified potential is roughly identical to the real potential for z < z₀ but is less steep for z > z₀ leading to a smoother density profile in the hydrostatic equilibrium and therefore smaller Alfvén velocities. It should be noticed that this problem arises only in the hydrostatic equilibrium without outflow. As we will show below, when an outflow is produced, the density profile is much smoother and we do not need a modified potential anymore.

We have used the above modified potential with z₀ = 4 in our 1D simulations to initiate the MHD wind in the tall box simulations (z_B ≥ 8). Once the quasi steady wind was formed, we relaxed back the potential to the original shearing box potential. Comparisons between the solutions obtained with and without the modified potential in the case z_B = 6 have shown no difference once the potential has been relaxed (runs 1DRef-1Dz6).

3. One-dimensional MRI outflows

In this section, we look at one-dimensional solutions of equations (1)−(3) i.e. (ρ(z,t), v(z,t), B(z,t)). This simplification allows us to isolate the physical mechanisms responsible for the MHD-driven outflows.

We initialise our box with a disc in vertical hydrostatic equilibrium (9). We add a mean vertical field B_z so that μ = 8 × 10^-2. To initialise the growth of MRI modes, we add a small perturbation B_x = 0.02sin(z) to the system. We show the temporal evolution of this run (1DRef) in Fig. 1.

3.1. From MRI modes to steady outflows

At it can be seen from Fig. 1, we first observe the development of a linear MRI mode in the simulation box (t = 4.0). These linear modes were described extensively by Latter et al. (2010) for stratified shearing boxes. The mode we observe in our particular setup is the n = 2 mode. This can be easily checked looking at the shape of the perturbation and comparing to the eigenmodes of Latter et al. (2010). Moreover, for μ = 8 × 10^-2, the n = 2 mode is the most unstable mode in the system (Fig. 2). Because the magnetic fluctuations are localised at z ~ H, the magnetic pressure tends to increase at that location, pushing toward the midplane the bulk of the disc and pushing up the disc atmosphere. At t = 7.0−8.0, the magnetic pressure is sufficiently large to push the atmosphere, creating a “bubble” of material. At this stage, the flow is no longer in a linear regime: the eigenmode is modified by nonlinear (e.g. magnetic pressure) terms. In the end, the system relaxes toward a quasi stationary state (t = 60). This implies that energy which is in injected by the MRI into the eigenmode is evacuated at the same rate in the outflow.

It should be pointed out that we have totally ignored secondary instabilities in this description. These secondary instabilities are usually x,y dependent modes which grow on the top of MRI eigenmodes (Goodman & Xu 1994; Latter et al. 2009; Pessah & Goodman 2009; Pessah 2010). It is widely believed that parasitic modes are responsible for the breakup of MRI modes into 3D turbulence, although the exact role they play is still controversial (Lesur & Longaretti 2011). Evidently in our 1D model, parasitic instabilities are inhibited, allowing the primary MRI mode to grow virtually for ever. To allow for the growth of parasitic modes, we have reproduced the same setup in 3D, adding 3D noise at moderate level ( ⟨ v² ⟩  ~ 10^-2) to the 1D perturbation described above. This systematically led to the production of an outflow before parasitic instabilities could do anything to the MRI eigenmodes.

Although surprising at first sight, this result can be understood using the phenomenology of parasitic modes. First, it should be noted that the most unstable parasitic modes are usually Kelvin-Helmholtz modes (Latter et al. 2009) growing on MRI modes. The maximum growth rate of Kelvin-Helmholtz modes can be estimated by the local vertical shear rate. If we consider a primary MRI mode of amplitude δv with a characteristic vertical size δl, then the maximum growth rate of the secondary mode is γ_S ~ δv/δl. For the secondary mode to have an impact on the MRI mode, we require γ_S > γ which implies δv > δlγ. Moreover, following Latter et al. (2010), we have δb ~ B_z0δv/ΩH where δb is the magnetic perturbation of the MRI mode and B_z0 is the mean vertical field. Therefore, the parasitic mode can destroy the MRI mode only if $\frac{\mathit{\delta b}}{{\mathit{B}}_{\mathit{z}\mathrm{0}}}\mathit{>}\frac{\mathit{\delta l}}{\mathit{H}}\frac{\mathit{\gamma }}{\mathrm{\Omega }}\mathrm{·}$(11)In our system, the MRI mode characteristic length δl, growth rate γ and the mean field amplitude are all of the order of 1 (in code units, see Sect. 2.1). This implies that parasitic modes will appear only when δb > 1. However, as can be seen in Fig. 1, the outflow starts when δb ≲ 1, i.e. before parasitic modes could act on the MRI mode. This result is due to the relatively large magnetisation used in these simulations (μ ≲ 1) compared to traditional MRI setups (μ ≲ 10^-3) . This large magnetisation implies the production of a large scale MRI mode whose growth rate is of the order of the orbital timescale. This makes the outflow more favourable compared to secondary instabilities. On the contrary, when the magnetisation is small, dominant MRI modes are found at smaller scale (large n, see Fig. 2). In this case, parasitic modes are favoured and a turbulent flow is obtained.

We should emphasize that the presence of an outflow does not mean that parasitic instabilities are totally absent from this picture. As we will show later (Sect. 4), solutions exhibiting an outflow are themself subject to parasitic instabilities. However, these instabilities have nothing in common with the traditional parasitic instabilities of MRI eigenmodes.

We have seen above that the evolution of a large-scale MRI mode in a strongly magnetised shearing box leads naturally to the production of a magnetically-driven outflow. This steady outflow is essentially one-dimensional and can be described by v(z),ρ(z),B(z). In the following we will concentrate on the structure of this outflow: phenomenology, critical points, boundary conditions and conserved quantities.

3.2. Outflow phenomenology

As we have shown above, the outflow is primarily produced by the magnetic pressure gradient. The magnetic pressure being maximum at z ~ 1.5, it pushes up the outflow at z > 2.0 but it also compresses the bulk of the disc. An alternative view of this effect can be obtained looking at currents. The outflow is in this case due to horizontal currents which are reversed at z ~ 1.5. We typically have J_x > 0 and J_y > 0 in the bulk of the disc whereas J_x < 0 and J_y < 0 in the atmosphere z > 2. It is important to note that the outflow acceleration can occur only if these currents are non-zero and change sign. This remark justifies the absence of any outflow with “zero gradient” boundary conditions (see Sect. 3.4).

It is of interest to put these outflow solutions in the context of the Blandford & Payne (1982) disc wind paradigm. In this model, the outflow is initiated by a magnetocentrifugal effect: the poloidal magnetic field lines are considered as rigid wires anchored in the bulk of the disc and fluid particles are allowed to drift along these wires. If the field lines are sufficiently inclined (typically more than 30° with respect to the vertical axis), then the particles are azimuthally accelerated by the anchored field lines. This leads to a centrifugal effect which ejects fluid particles along field lines. In this picture, the ejection is driven by an exchange of angular momentum: angular momentum is taken from the disc by the field line and it is then released in the ejected material.

To compare this mechanism to outflow solutions driven by MRI modes, we show in Fig. 3 the poloidal streamlines and field lines of such a solution. We first note that poloidal streamlines and fieldlines are not aligned. This property is allowed by the shearing box boundary conditions, but it physically means that magnetic flux is “accreted” toward the centre of the disc². In a more realistic model including curvature, such a state could not be sustained for a long time because magnetic flux would get accumulated at the disc centre, thereby modifying the disc properties (especially its rotational profile). This is the principal motivation for the presence of a strong “magnetic diffusivity” in disc wind models, either assuming the presence of small scale turbulence (Ferreira & Pelletier 1993a), ambipolar diffusion (Wardle & Königl 1993) or Hall and Ohm diffusion (Königl et al. 2010).

Despite this difference, we recover most of the phenomenological properties of the Blandford & Payne (1982) paradigm: field lines are inclined and drive an outflow which is inclined towards the same direction. This indicates that angular momentum is exchanged between the field and the flow. As we will see below (Sect. 3.5.3), angular momentum is effectively taken away from the disc by the field and then released to the ejected material. This effect is actually inevitable since the current configuration described previously inevitably leads to a positive magnetic torque  ∝ J_zB_x − J_xB_z in the outflow (and a negative one in the disc midplane). Because angular momentum is taken from the disc by the field, a strong radial flow is produced which explains the streamlines’ orientation for z < 1. Finally, we find an inclination angle of  ~40° for the poloidal magnetic field lines and  ~25° for the poloidal velocity field at the top of the box. This last value is very close to the critical value of Blandford & Payne (1982).

3.3. Critical points

In principle, it is possible to look systematically for a steady 1D solution of Eqs. (1) − (3). This is done writing the equations of motion in the form M·X = Y, where $\begin{array}{ccc}\mathbf{M}& \mathrm{=}& \left(\begin{array}{c}\\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \end{array}\right)\\ \mathbf{X}& \mathrm{=}& {\mathit{\partial }}_{\mathit{z}}\left(\begin{array}{c}\\ \\ \\ \\ \\ \\ \end{array}\right)\mathrm{;} \mathbf{Y}\mathrm{=}\left(\begin{array}{c}\\ \\ \\ \\ \\ \\ \end{array}\right)\end{array}$where M and Y are a matrix and a vector which do not contain any spatial derivative. To solve this nonlinear system, one then inverts the above system of equations to get a set of ordinary differential equations X = M^-1·Y. However, when a critical point is reached, M is singular and the system cannot be inverted anymore. In this case, the physical system needs to satisfy an extra condition in order to get through the critical point.

In the shearing box, we find three types of singular points:

• two slow magnetosonic points${\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{V}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}^{\mathrm{2}}\mathrm{-}\sqrt{{\mathrm{(}{\mathit{V}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{c}}^{\mathrm{2}}{\mathit{V}}_{\mathrm{A}\mathit{z}}^{\mathrm{2}}}\right]\mathit{,}$(14)

• two Alfvén points ${\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{=}{\mathit{V}}_{\mathrm{A}\mathit{z}}^{\mathrm{2}}\mathit{,}$(15)

• and two fast magnetosonic points ${\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{V}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}^{\mathrm{2}}\mathrm{+}\sqrt{{\mathrm{(}{\mathit{V}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{c}}^{\mathrm{2}}{\mathit{V}}_{\mathrm{A}\mathit{z}}^{\mathrm{2}}}\right]\mathit{,}$(16)

where ${\mathit{V}}_{\mathrm{A}}\mathrm{=}\sqrt{{{B}}^{\mathrm{2}}\mathit{/}\mathit{\rho }}$ and ${\mathit{V}}_{\mathrm{A}\mathit{z}}\mathrm{=}{\mathit{B}}_{\mathit{z}}\mathit{/}\sqrt{\mathit{\rho }}$. In the following, thanks to symmetries, we will only consider solutions with v_z > 0 so that only one critical point of each kind will be present.

We present the MHD wave speeds and flow speeds in Fig. 4. We find that the slow point is located around z = 0.52 and the Alfvén point is found at z = 2.47. The flow does not cross the fast point, however we find that the fast speed and the flow vertical speed tend to converge more rapidly close to the upper boundary. Note that a very similar behaviour was observed in a global self-similar solution close to the fast surface (Casse & Ferreira 2000, Fig. 1). This result indicates that the flow is still causally connected to the disc and therefore the boundary condition we impose at the top of the box still has an impact on the flow structure itself. This point will be discussed in the next section.

3.4. Influence of the upper boundary condition

We have seen above that the outflow is still physically connected to the disc since it is not super-fast. Moreover, it looks as if the fast magnetosonic point is located close to the imposed upper boundary condition. One might wonder what is the exact role played by this boundary condition. To investigate this issue, we have performed two kinds of test, either modifying the altitude z_B of the vertical boundary condition or modifying the nature of the upper boundary condition we apply.

First, varying the altitude of the vertical boundary condition z_B led to the results presented in Table 2. One finds that modifying the altitude of the boundary condition strongly modifies the outflow properties. In particular, we find that increasing the altitude of the boundary condition leads to a decrease in the flow ejection rate. This evolution is accompanied by a slow point moving closer to the midplane and an Alfvén point moving higher up in the atmosphere. This result clearly demonstrates that because the flow has not crossed the fast magnetosonic point, the solution we obtain is still constrained by the boundary condition we impose at z_B. In all the solutions described in Table 2, we have observed a “convergence” of the fast magnetosonic speed and the flow speed when approaching z_B, similarly to what we find in Fig. 4. This surprising result tends to indicate that our boundary conditions somehow force the fast point to be close to the upper boundary.

The fact that the slow point moves closer to the disc midplane when the mass loss rate decreases might look dubious to readers familiar with the phenomenology usually associated with the slow point. In particular, it is often said that the slow point “sets” the wind escape speed through the relation ρv_z = ρ₀(z_s)c where ρ₀ is taken as the hydrostatic density profile (Spruit 1996). This argument predicts that smaller mass loss rate should be associated to slow points located higher up in the atmosphere, which is exactly the opposite of what we find. This simple argument is however not valid in the present case since the density profile deviates significantly from the hydrostatic equilibrium. In particular, configurations 1Dz16 and 1Dz20 exhibit strongly compressed discs due to a large magnetic pressure gradient in the atmosphere which most probably prevents significant ejection from happening.

We have also tried to modify the nature of the upper boundary conditions. First, instead of imposing B_x(z_B) = 0, we have imposed a fixed angle to the poloidal field, i.e. B_x(z_B) = tan(θ)B_z with θ = 30°;   45° as in Ogilvie (2012). Surprisingly this did not modify significantly the outflow solution we obtained: the field values are modified by less than 5%. This can be explained by the fact that the outflow is super-Alfvénic when it reaches the top boundary. As we will see below (Sect. 3.6), sub-Alfvénic outflows are effectively very sensitive to the field configuration at the boundary, but super-Alfvénic outflows are not. We conclude from this that the inclination angle of the poloidal field line is set by the Alfvén point crossing condition. This result is corroborated by the constant inclination angle at the Alfvén surface found when one changes the box vertical size (Table 2)

We have finally tried to impose a zero gradient condition on B_x and B_y (classically called “outflow” boundary condition). This results in the suppression of the outflow solution. We observe instead a constant increase of the magnetic pressure in the atmosphere which results in a strong compression of the disc material in the midplane until the disc occupies one numerical grid cell. This result is similar to the low β simulations of Hawley et al. (1995) with mean vertical flux. This was to be expected since the outflow is driven by horizontal currents. Imposing a zero gradient condition means that no horizontal current is allowed at the boundary, blocking any potential outflow by cancelling the vertical component of the Lorentz force and the change of sign of the magnetic torque (Ferreira & Pelletier 1993b).

3.5. Conserved quantities

In the following, we assume all the quantities v, B and ρ depend only on z, as found in the steady ejection above. With this hypothesis, we reconstruct the conserved quantities used in global disc wind models (Blandford & Payne 1982; Pelletier & Pudritz 1992; Casse & Ferreira 2000).

3.5.1. Frozen field condition

We first note that under the above conditions, B_z and ρv_z are constant in the box thanks to flux and mass conservation. We therefore introduce a proportionality constant between these quantities: $\mathit{\rho }{\mathit{v}}_{\mathit{z}}\mathrm{=}\mathit{\alpha }{\mathit{B}}_{\mathit{z}}\mathit{.}$(17)Thanks to the x induction equation, we also have v_xB_z − v_zB_x = ℰ_y = const., which can be simply written as: ${\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}\mathrm{=}{\mathit{v}}_{\mathit{x}}\mathrm{-}{\mathit{v}}_{\mathit{z}}\frac{{\mathit{B}}_{\mathit{x}}}{{\mathit{B}}_{\mathit{z}}}$(18)where ${\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}$ is a constant. Physically, ${\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}$ is the advection speed of the poloidal magnetic field. In global models, this quantity is implicitly set to 0 in order to avoid flux accumulation at the disc inner edge and potential singularities at R = 0 (Chandrasekhar 1956; Mestel 1961; Pelletier & Pudritz 1992). This is not required here as we are using a shearing-box model. However, one should keep in mind that our model implicitly allows radial advection of magnetic field lines. The existence of this conserved quantity allows us to define a relation between the poloidal field and the poloidal velocity, namely: ${{v}}_{\mathrm{p}}\mathrm{=}\frac{\mathit{\alpha }}{\mathit{\rho }}{{B}}_{\mathrm{p}}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}{{e}}_{{x}}\mathit{,}$(19)where v_p and B_p are the poloidal (x,z) components of the velocity and magnetic fields. This equation shows a major difference between the global approach and the local approach. In the global approach, ${\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}\mathrm{=}\mathrm{0}$ (no advection of magnetic field lines), which implies that poloidal streamlines and magnetic field lines are aligned. This is not the case in our solutions, for which ${\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}\mathrm{=}{\mathit{v}}_{\mathit{x}}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{<}\mathrm{0}$.

3.5.2. Magnetic surface rotation

A relation, known as the magnetic surface rotation speed, can be deduced from the y component of the induction equation. Since B_y does not depend on x, we get after some algebra ${{B}}_{\mathrm{p}}\mathrm{·}\mathrm{\nabla }\left[\frac{\mathit{\alpha }}{\mathit{\rho }}{\mathit{B}}_{\mathit{y}}\mathrm{-}{\mathit{v}}_{\mathit{y}}\mathrm{+}\mathit{q}\mathrm{\Omega }\mathit{x}\right]\mathrm{=}\mathrm{0.}$(20)Note that the above relation is only valid along a poloidal magnetic field line. In contrast to Pelletier & Pudritz (1992), no equivalent relation is found along poloidal streamlines. This allows us to define a new conserved quantity ${\mathit{v}}_{\mathit{y}}^{\mathrm{\ast }}\mathrm{=}{\mathit{v}}_{\mathit{y}}\mathrm{-}\mathit{q}\mathrm{\Omega }\mathit{x}\mathrm{-}\frac{\mathit{\alpha }}{\mathit{\rho }}{\mathit{B}}_{\mathit{y}}$(21)which is conserved along a magnetic field line. It should be pointed out at this stage that ${\mathrm{ℰ}}_{\mathit{x}}\mathrm{=}{\mathit{u}}_{\mathit{y}}{\mathit{B}}_{\mathit{z}}\mathrm{-}{\mathit{u}}_{\mathit{z}}{\mathit{B}}_{\mathit{y}}\mathrm{=}{\mathit{B}}_{\mathit{z}}{\mathit{v}}_{\mathit{y}}^{\mathrm{\ast }}$ where ℰ_x is the x component of the electromotive force. Therefore, in the frame translating at ${\mathit{v}}_{\mathit{y}}^{\mathrm{\ast }}$, the x component of the electromotive force is zero, which justifies the naming of this invariant.

3.5.3. Angular momentum conservation

The angular momentum conservation equivalent is derived from the y component of the equation of motion (2). It should first be noted that this equation can be written: $\mathit{\rho }{{v}}_{\mathrm{p}}\mathrm{·}\mathrm{\nabla }\mathrm{ℒ}\mathrm{-}{{B}}_{\mathrm{p}}\mathrm{·}\mathrm{\nabla }{\mathit{B}}_{\mathit{y}}\mathrm{=}\mathrm{0}$(22)where ℒ = v_y + (2 − q)Ωx is the local equivalent of the specific angular momentum. Since B_y does not depend on x, we can rewrite the above equation using the flux freezing condition: $\mathit{\rho }{{v}}_{\mathrm{p}}\mathrm{·}\mathrm{\nabla }[\mathrm{ℒ}\mathrm{-}\frac{{\mathit{B}}_{\mathit{y}}}{\mathit{\alpha }}]\mathrm{=}\mathrm{0}\mathit{,}$(23)which indicates that $\mathit{f}\mathrm{=}\mathrm{ℒ}\mathrm{-}\frac{{\mathit{B}}_{\mathit{y}}}{\mathit{\alpha }}$ is conserved along streamlines.

The angular momentum along one streamline is shown in Fig. 5. This allows us to check that the angular momentum is effectively conserved in our simulations. Moreover, we find that most of the angular momentum is initially extracted from the disc by the toroidal field B_y. Higher above the disc (z ~ 2H), the angular momentum is exchanged between the toroidal field and the accelerated gas. This demonstrates that the magnetocentrifugal acceleration effect, present in the Blandford & Payne (1982) phenomenology, is also found in our outflow solutions.

3.5.4. Bernoulli invariant

The Bernoulli invariant is derived from Eq. (2) after a scalar multiplication by u. One gets: $\mathit{\rho }{{u}}_{\mathbf{p}}\mathbf{·}\mathbf{\nabla }\left[\frac{{\mathit{u}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{H}\mathrm{+}\mathit{\psi }\right]\mathrm{=}{J}\mathbf{\times }{B}\mathbf{·}{u}$(24)where H is the flow enthalpy. The right-hand side of this equation corresponds to the work done by magnetic forces. This term can be calculated as a function of the previously defined invariants which, after some algebra leads to: $\mathit{\rho }{{u}}_{\mathbf{p}}\mathbf{·}\mathbf{\nabla }[\frac{{\mathit{u}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{H}\mathrm{+}\mathit{\psi }\mathrm{-}\frac{{\mathit{B}}_{\mathit{y}}{\mathit{v}}_{\mathit{y}}^{\mathrm{\ast }}\mathrm{-}{\mathit{B}}_{\mathit{x}}{\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}}{\mathit{\alpha }}]\mathrm{=}\mathrm{-}\mathit{q}\mathrm{\Omega }\mathit{\rho }{\mathit{B}}_{\mathit{y}}{\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}\mathit{.}$(25)Because v_p and B_p are not strictly speaking collinear, Bernoulli’s equation is left with a non conservative term. However, as we will see later, this term is important only in the bulk of the disc, so that the flux defined above will be approximately invariant along a streamline. Comparing this invariant to the one defined in global geometry, we note the presence of a new term in our case, ${\mathit{B}}_{\mathit{x}}{\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}\mathit{/}\mathit{\alpha }$. As shown before, this term describes the energy associated with the field lines being advected.

We show the different terms involved in the Bernoulli invariant Eq. (25) in Fig. 6. One should note that the magnetic contribution is separated in two parts: (i) the conservative part $\mathrm{(}{\mathit{B}}_{\mathit{y}}{\mathit{v}}_{\mathit{y}}^{\mathrm{\ast }}\mathrm{-}{\mathit{B}}_{\mathit{x}}{\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}\mathrm{)}\mathit{/}\mathit{\alpha }$ and (ii) the non-conservative part ${}^{\mathrm{\int }}\mathrm{d}\mathit{l}\hspace{0.17em}\mathit{q}\mathrm{\Omega }\mathit{\rho }{\mathit{B}}_{\mathit{y}}{\mathit{v}}_{\mathit{x}}^{\mathrm{\ast }}$ where the integral is computed along a streamline.

We first find that the invariant is approximately conserved for z > 0.2. Initially (up to z ~ 1.5) the energy is stored in the conservative component of the magnetic field. The non-conservative part is constantly decreasing, demonstrating that this term is helping the outflow. Higher in the outflow, the magnetic energy is converted into kinetic and potential energy. Finally, we find that the thermal energy plays essentially no role in the ejection energetics. This demonstrates that the isothermal approximation used in this work does not inject a significant amount of energy in the outflow.

The fact that potential energy increases along a streamline might look at first surprising since in a typical Blandford & Payne (1982) situation, one would expect the potential energy (gravitational + centrifugal) to decrease along a field line. This is not the case here because the inclination angle of the outflowing streamlines is slightly less than 30° (see Sect. 3.2).

3.6. Magnetisation dependency

In the previous discussion, we kept a constant equivalent surface density Σ by artificially injecting mass in the disc midplane. This approximation, although partly motivated by the global disc structure, should be tested more extensively. To do so, we have performed simulations without mass injection. By definition these simulations cannot achieve a steady state. They are however representative of an extreme case in which no mass is coming from the outer disc.

It should first be pointed out that the outcome of these simulations depends strongly on the nature of the upper boundary condition. This is because as the disc mass is lost, the Alfvén point moves higher up in the atmosphere (see Fig. 7). At some point, the Alfvén point crosses the upper boundary and the ejection becomes sub-Alfvénic. When this happens in this simplified setting (no radial dependence), the poloidal field inclination is not set anymore by the Alfvén point crossing condition (Ogilvie 2012, see also Sect. 3.4). Instead, it should be set manually at the upper boundary. If we set B_x(z_B) = 0, the ejection stops as soon as the Alfvén point exits the simulation domain. This is expected since no ejection should be happening with strictly vertical poloidal field lines. On the contrary, if we impose a fixed inclination angle of 45° at the upper boundary condition, as in Ogilvie (2012), the outflow is kept. In the following we will consider the latter case.

Because the disc is losing mass, its magnetisation μ is increasing as a function of time (see Fig. 8). We note however, that the system appears to be approaching a steady state with μ ≃ 2.3. Looking at the correlation between the mass loss rate due to the wind Ṁ_W = (ρv_z)_{z    =    zB} and the magnetisation μ (Fig. 9), we find that the mass loss rate decreases steeply for μ > 1, which explains the quasi-steady state we observe. We note that this last result is very similar³ to Ogilvie (2012) (Fig. 4). This demonstrates a plausible causal transition between MRI-driven outflow solutions described here and the solutions found by Ogilvie (2012) in the MRI stable regime.

3.7. Comparison with a global outflow solution

Because the shearing box model is subject to several restrictive hypothesis, it is useful to compare the MRI driven outflows found here with global steady self-similar solutions of disc winds. The idea here is not to have a quantitive agreement between these solutions since they are computed in very different different configurations. Instead, we look for qualitative similarities which could indicate that the physical processes at work are similar. To this end, we have used the numerical procedure described by Ferreira & Pelletier (1995) and Casse & Ferreira (2000) to compute an isothermal “cold” solution with the following parameters⁴: ϵ = 0.1; α_m = 1.0; Pm = 0; μ = 0.7268; ξ = 9.92 × 10^-3. This last parameter is connected to the radial dependence of the accretion rate through the relation Ṁ_A(r) ∝ r^ξ. Although this solution depends explicitly on r and z through a self-similar ansatz, the radial dependence is much weaker than the vertical one up to a few scale heights. Therefore, we only show the z dependence of the resulting solution in Fig. 10. We have zoomed on the region from z = 0...6H to allow a comparison with shearing box solutions (Fig. 1), although the actual global solution extends up to 25H (fast magnetosonic point).

The direct comparison between this global solution and the shearing box solutions shows several differences: the velocities, and in particular v_z, are larger in the global solution. On the other hand, the horizontal magnetic field strength is weaker compared to the shearing box solution. However, we recover most of the physical properties found in global solutions in the shearing box solution: a strong magnetic shear close to the midplane accompanied by a magnetically compressed disc (compared to a purely hydrostatic equilibrium) and a strong accretion flow in the disc midplane (v_x ~ −2c).

It should be noted that the global solution presented here has a much weaker mass loss rate (ρv_z ~ 5 × 10^-3) compared to the shearing box solution (ρv_z ≃ 0.1). This can be partly due to the stronger mean vertical field required by the global solution (μ ~ 0.24), although this difference is not enough to explain entirely this discrepancy (Fig. 9). The fact that the mass loading is much smaller in the global solution explains the faster escape velocity found in this solution.

4. Three-dimensional solution and stability

4.1. Outflow evolution in 3D

To investigate the 3D evolution of the 1D outflow described in Sect. 3.1, we have perform the outflow simulations in 3D. The initial condition consists of the 1D initial perturbation described in Sect. 3 plus a random 3D noise added to v_x with max(|v_x|) = 0.02. We show in Fig. 11 the temporal evolution of simulation 3DRef for μ = 8 × 10^-2. These figures represent the evolution of the density, poloidal magnetic field inclination and vertical velocity averaged in the (x,y) plane in a (z,t) diagram. We first observe the presence of a strong outburst (t ~ 8) followed by a rather steady state during which the flow does not evolve rapidly. This state corresponds to the 1D solution described in Sect. 3 and is essentially a 1D outflow solution. However, after some time, this 1D outflow goes unstable (t ~ 40) and rapid variations in all quantities are observed. Surprisingly, the global structure of the 1D outflow is conserved: on time and horizontal averages, the typical v_z vertical and poloidal inclination profiles are consistent with the 1D steady solution (Fig. 1). We therefore produce a turbulent outflow driven by the Blandford & Payne (1982) magnetocentrifugal effect. We have additionally performed simulations similar to 3DRef removing the symmetry condition at z = 0. An unstable pattern having a similar growth rate as the one described here is found in these cases, demonstrating that this instability is neither driven nor limited by the symmetry condition we impose at z = 0.

We show in Fig. 12 snapshots of the same simulation taken at t = 20, t = 40 and t = 270. At t = 20, we confirm the 1D nature of the solution: no dependence can be seen in x or y. This also shows that the outflow is produced before parasitic modes get significantly excited. At t = 40 we see the development of a “secondary” instability of the outflow solution which produces a x-dependent structure. The physical processes responsible for this instability will be discussed in Sect. 4.2. When the “secondary” modes reach a significant amplitude, the structures break down into non-axisymmetric turbulent motions in which the main outflow properties are maintained (field line inclination, outflow speed). A typical example of such a state is shown at t = 270.

4.2. Outflow solution stability

We have shown above that the outflow solution is unstable to 3D perturbation, resulting in a “turbulent outflow” configuration. Several authors have discussed the possibility of having such an instability (Lubow et al. 1994; Lubow & Spruit 1995; Cao & Spruit 2002) although the applicability of these stability analyses to all outflow solutions is still uncertain (Königl & Wardle 1996; Königl 2004).

Instabilities comparable to the one described here were also found by Moll (2012). The evolution timescales and the global shape of unstable perturbations (“striped wind”) look similar in both cases, although the radial wavelength of Moll (2012) instability seems smaller than ours. It is therefore probable that all these instabilities have the same physical origin, though a precise comparison of the growing eigenmodes is required to confirm this point.

To analyse the instability observed in our outflow solutions, we have performed a simulation (3DLin) starting from the 1D solution corresponding to the final state of run 1DRef to which a small-amplitude (10^-3) 3D white noise was added. Since the growth phase of the instability implies only x-dependent modes, we use a Fourier decomposition in the x direction to characterise growing modes: $\begin{array}{ccc}& & \mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{+}\mathit{\delta \rho }\mathrm{\left(}\mathit{x,z,t}\mathrm{\right)}\\ & & \mathit{\delta \rho }\mathrm{\left(}\mathit{x,z,t}\mathrm{\right)}\mathrm{=}\mathrm{\Re }\left[\sum _{\mathrm{K}}{\mathit{\rho }}_{\mathrm{K}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\exp \mathrm{\left(}\mathrm{i}\mathit{kx}\mathrm{\right)}\right]\end{array}$and similarly for v and B. In this expansion, we have assumed the instability was growing on the top of the 1D solution ρ₀(z) given by the final state of run 1DRef. We first present the temporal evolution of the fluctuation $\sqrt{\mathrm{⟨}\mathit{\delta }{\mathit{\rho }}^{\mathrm{2}}\mathrm{⟩}}$ in Fig. 13. We find that perturbations grow exponentially with a growth rate γ = 0.33Ω up to t ≃ 40 where a saturation is reached. To further investigate the behaviour of this instability, we present the temporal evolution of  |ρ_K(z = 0.5)|  as a function of k in Fig. 14. As it can readily be seen, the growth is dominated by modes having n ≡ kL_x/2π = 4 which is consistent with the 4 “spots” observed in Fig. 12. The measured growth rate of this mode is γ = 0.33Ω which explains its fast appearance in 3D simulations once an outflow has formed. However, other neighbouring modes are growing as well (n = 3;5;6) although not as fast as the n = 4 mode. Finally once the n = 4 mode reaches large amplitudes (|ρ_K |  ~ 1), we note the sudden growth of the n = 8 and n = 9 modes which are the result of nonlinear interactions of the fastest growing modes n = 3;4;5;6.

To analyse the physics underlying the n = 4 mode, we present in Fig. 15 the density fluctuations corresponding to that eigenmode. We first note that the density fluctuation is highly localised in z around z ~ 0.5. In addition, this eigenmode has a tail whose inclination and shape closely follow that of the poloidal magnetic field (Fig. 3). The localisation of the eigenmode is rather surprising and requires some explanation. We first note that this localisation is much higher than the top of the mass injection region (z_inj = 0.1). However, comparing the relative vorticity component ω_y = ∂_zv_x of the 1D outflow (Fig. 16) to the vertical profile of the eigenmode (Fig. 17), we find that the density perturbation is localised close to a maximum of ω_y. This tends to suggest that this instability is somehow linked to the vertical ω_y profile of the outflow solution, and therefore to a kind of Kelvin-Helmholtz instability. This potential link is also consistent with the growth rate γ ≲ max(ω_y) and with the shape of the eigenmode and the vorticity profile around z = 0.5.

We would like to stress that these remarks are not a proof that this outflow instability is of the Kelvin-Helmholtz type. Among the effects we did not take into account in that analysis are the magnetic field, compressibility and the presence of a large v_z up in the atmosphere. However, if we assume that the source of the instability lies around z ~ 0.5 as suggested by the eigenmodes, then timescales, length scales and phenomenology match that of a Kelvin-Helmholtz instability. To ascertain these claims, a proper linear study taking into account all of the outflow properties is required, which is well beyond the scope of this paper.

5. Discussion

In this paper, we have shown that large scale MRI modes which are unstable when the disc magnetisation is moderately sub-thermal spontaneously produce super-Alfvénic outflows. The physical mechanism behind this outflow is a Blandford & Payne-like process where angular momentum is transferred to bending field line and then released to accelerated material. We demonstrated that this outflow is qualitatively similar to the outflow solutions found both in local boxes and in global self-similar geometry, making a clear connection between the MRI and the formation of disc winds. We have also shown that MRI outflows are unstable in 3D which could be a potential source of variability for disc winds and jets. These 3D instabilities could also be the origin of the turbulent resistivity α_m used in Ferreira & Pelletier (1995) and subsequent self-similar models.

However, the picture provided here is still incomplete. We first note that the simulations we have performed here only produce super-Alfvénic flows, whereas real escaping outflows should be faster than the fast magnetosonic speed. This problem seems to be linked to the shearing box geometry, as several other authors have noticed this in shearing boxes (Suzuki & Inutsuka 2009; Fromang et al. 2013; Bai & Stone 2013). Indeed, a shearing box does not allow for the opening of magnetic field lines which is expected in realistic global geometry. We believe that such an opening is required to get super-fast flows, which are therefore beyond the scope of our shearing box model. The fact that the flow is only super-Alfvénic indicates that the upper boundary condition still plays a role in determining the outflow structure, and indeed it does, at least for the mass-loss rate (Sect. 3.4).

We also emphasize that the shearing box model possesses a horizontal symmetry by the transformation (v_x,v_y) → (−v_x, − v_y), (B_x,B_y) → (−B_x, − B_y) and (x,y) → (−x, − y). This symmetry indicates that locally, there is no mathematical difference between a magneto-centrifugally accelerated wind where angular momentum is transferred from the disc to the jet and a magnetically decelerated accretion column (formally a magneto-centripetal wind) where the angular momentum is transferred from material falling radially inward to the disc. This symmetry is obviously broken once curvature terms are taken into account. The presence of this symmetry implies that shearing box simulations can spontaneously switch from one situation to the other, which is unexpected in realistic situations. Such sudden changes in the magnetic configuration were indeed observed in rare occasions in our 3D runs but also by Fromang et al. (2013).

Finally, we should point out that the presence of a non-zero toroidal electromotive force implies that magnetic field lines are accreted. As mentioned earlier, this situation is rather unrealistic in global geometry, although it is allowed in shearing boxes. More realistic configurations, probably including a sort of resistivity (either effective or molecular) are required to compensate for this effect.

All of these points indicate that shearing boxes are not very well suited to study globally the outflows produced by MRI turbulent accretion discs. In particular, little can be learned regarding outflow mass loss rate and velocities. We note however that our solutions can be qualitatively compared to global solution and several properties are recovered by the local model. Moreover, 3D instabilities can be studied much more easily in boxes than in global geometry where computational costs increase very rapidly.

Overall, our results tend to indicate a paradigm shift: up to now, the MRI driving “viscous” discs and disc winds at the origin of jets have often been considered as separate processes. Here we show that these two processes are actually intrinsically connected: outflows are a logical consequence of the MRI in strongly magnetised discs. Obviously, the next question is to understand what is driving the disc magnetisation. This could be due to local dynamos, large scale field redistribution through advection and diffusion, etc. To identify these processes, shearing boxes are clearly insufficient and global models, including both large-scale magnetic fields and turbulence, will be required.

---

Acknowledgments

G.L. thanks Sébastien Fromang for his comments and suggestions on the initial version of this manuscript, Céline Combet for providing the code used to compute the self-similar solutions presented in Sect. 3.7 and Andrea Mignone for his help and advices regarding the Pluto code. This work was granted access to the HPC resources of IDRIS under allocation x2012042231 made by GENCI (Grand Équipement National de Calcul Intensif). Some of the computations presented in this paper were performed using the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhone-Alpes region (GRANT CPER07_13 CIRA: http://www.ci-ra.org). G.L. acknowledges support by the European Community via contract PCIG09-GA-2011-294110.

References

All Tables

All Figures