Issue 
A&A
Volume 625, May 2019



Article Number  A108  
Number of page(s)  18  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201834813  
Published online  22 May 2019 
Spontaneous ring formation in windemitting accretion discs
Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France
email: antoine.riols@univgrenoblealpes.fr
Received:
10
December
2018
Accepted:
14
April
2019
Rings and gaps have been observed in a wide range of protoplanetary discs, from young systems like HLTau to older discs like TW Hydra. Recent disc simulations have shown that magnetohydrodynamic (MHD) turbulence (in both the ideal or nonideal regime) can lead to the formation of rings and be an alternative to the embedded planets scenario. In this paper, we have investigated the way in which these ring form in this context and seek a generic formation process, taking into account the various dissipative regimes and magnetisations probed by the past simulations. We identify the existence of a linear and secular instability, driven by MHD winds, and giving birth to rings of gas that have a width larger than the disc scale height. We show that the linear theory is able to make reliable predictions regarding the growth rates, the contrast and spacing between ring and gap, by comparing these predictions to a series of 2D (axisymmetric) and 3D MHD numerical simulations. In addition, we demonstrate that these rings can act as dust traps provided that the disc is sufficiently magnetised, with plasma beta lower than 10^{4}. Given its robustness, the process identified in this paper could have important implications, not only for protoplanetary discs but also for a wide range of accreting systems threaded by largescale magnetic fields.
Key words: accretion, accretion disks / protoplanetary disks / magnetohydrodynamics (MHD) / instabilities / turbulence
© A. Riols and G. Lesur 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The radiointerferometer ALMA and the new generation of instruments like SPHERE at the Very Large Telescope have imaged a variety of structures in protoplanetary discs around young stars (Garufi et al. 2017). One of the most striking features are the concentric rings (or gaps), observed in many discs: HL tau (ALMA Partnership 2015), TW Hydra (Andrews et al. 2016) or the disc around Herbig Ae star HD 163296 (Isella et al. 2016). These structures may influence the disc evolution and could be a privileged location of dust accumulation (Pinilla et al. 2012), a key step towards planetary cores formation
One important challenge in accretion discs theory is to understand the origin of these rings. The scenario commonly invoked is the presence of embedded planets forming and opening gaps (Kley & Nelson 2012; Baruteau et al. 2014; Dong et al. 2015). Although the planet hypothesis is attractive and seems consistent with recent simulations (Dipierro et al. 2015), it challenges the planet formation theory, in particular for young systems like HLTau (<1 Myr old). The core accretion model at distance of a few tens of AU indeed requires more than a million years to form planets (Helled & Bodenheimer 2014). A large number of alternative mechanisms has been suggested, such as dustdriftdriven viscous ring instability (Wünsch et al. 2005; Dullemond & Penzlin 2018), snow lines (Okuzumi et al. 2016), dead zones (Flock et al. 2015) or secular gravitational instabilities in the dust (Takahashi & Inutsuka 2014). One recent and appealing scenario is the formation of concentric rings and gaps by magnetohydrodynamics (MHD) processes in the disc.
Since the early 90s, it has been accepted that MHD processes are ubiquitous and crucial in the evolution of accreting systems. Magnetised discs, if sufficiently ionised, are indeed prone to the magnetorotational instability (MRI, Balbus & Hawley 1991; Hawley et al. 1995), leading to turbulence and angular momentum transport. When threaded by a mean vertical field, the disc may also evacuate a significant part of their angular momentum and energy through largescale winds intimately connected to the MRI (Lesur et al. 2013; Fromang et al. 2013). Even in poorly ionised regions (r ≳ 0.1−1 AU) subject to non ideal effects (ambipolar diffusion and Hall effect), accretion can operate via MHD winds, despite the absence of vigorous MRI turbulence (Bai & Stone 2013; Lesur et al. 2014; Bai 2015; Béthune et al. 2017).
Several simulations in the local and global configuration, including a mean vertical field, have brought evidence that MHD flows and their winds, selforganise into largescale axisymmetric structures or “zonal flows” associated with rings of matter (Kunz & Lesur 2013; Bai 2015; Béthune et al. 2016, 2017). These features appear predominantly in the presence of nonideal effects, but were also noticed in MRI simulations without any explicit diffusion (Steinacker & Papaloizou 2002; Bai & Stone 2014; Suriano et al. 2018). Recent works attempted to explain their origin, though without any persuasive outcome. Bai & Stone (2014) proposed that rings form through an “antidiffusion” associated with the anisotropy of MRI turbulence. However, their result appears in conflict with most of the simulations that measured turbulent magnetic diffusivities (Guan & Gammie 2009; Fromang & Stone 2009; Lesur & Longaretti 2009). Using global simulations, Suriano et al. (2018, 2019) suggest that the structures are formed via reconnection of pinched poloidal field lines in the midplane current layer. Nevertheless, there is a lack of evidence that this mechanism is generic and works in all magnetic configurations. It requires a particular symmetry, with a poloidal field bending in the midplane, which is not the geometry observed in many nonideal simulations.
In this paper, we bring evidence that the process forming rings and gaps in MHD simulations (ideal, resistive, or ambipolar) is generic and supported by a local wind instability. The instability requires a mean mass ejection and a radial transport of vertical magnetic flux, whose origin can be either the “α” viscosity or the zonal flow itself. In turbulent discs, the criterion for instability is that the mass loss rate increases faster than the stress with the disc vertical magnetisation. In the presence of turbulence, the mechanism is reminiscent of a “viscoustype” instability, similar to that imagined by Lightman & Eardley (1974), except that the mass is free to escape the disc vertically. In a sense, it is also similar to the winddriven instability proposed by Lubow et al. (1994) and Cao & Spruit (2002). Unlike the latter, however, it does not rely on the assumptions that angular momentum is removed by the wind magnetic torque, neither that the mass loss rate increases with the poloidal field inclination, which are arguable assumptions (see Konigl & Wardle 1996). We believe that the mechanism described in this paper is closely related to the “massflux” (or “stripe” instabilities) seen by Moll (2012) and Lesur et al. (2013) in the highly magnetised (MRIstable) regime.
This paper is organised as follows: in Sect. 2, we present the main characteristics of zonal flows (or rings) in MHD simulations and show that, counterintuitively, they do not form through a radial transport of matter, but appear as a consequence of gaps emptied by vertical outflows. In Sect. 3, we suggest the existence of a winddriven instability and describe how we calculated its growth rate theoretically. In Sect. 4, we perform 2D MHD simulations (with or without explicit diffusion) for a wide range of magnetisations to test the existence and properties of the instability. We show in particular that axisymmetric modes projected into the Fourier space grow exponentially, with welldefined growth rates corresponding to those predicted by the theory. We have also explored the nonlinear saturation of the instability and attempt to predict the contrast between rings and gaps and their radial separation. Finally, in Sect. 6, we discuss the potential implications of our work beyond the scope of stellar discs.
2. Phenomenology of selforganisation
2.1. Zonal flows and rings occurrence in MHD simulations
Selforganisation of the gas into “zonal flows” is found in various MHD simulations of discs threaded by a net vertical field. The term zonal flows refers to the succession of subKeplerian and superKeplerian bands, associated with large scale rings of matter, as revealed by the global simulations of Béthune et al. (2017). In the following, the terms ring or zonal flow will be used to refer to the same substructure. Figure 1 summarises ten years of research and numerical effort to identify the conditions under which zonal flows can develop.
Fig. 1.
Ring occurrence in MHD simulations with net vertical magnetic flux. The orange check mark means that zonal flow exist, but their persistence within the turbulent flow and their convergence with box size are uncertain at current time. (1) Suriano et al. (2017), (2) Bai & Stone (2014), (3) our own simulations, (4) Hawley (2001), (5) Kunz & Lesur (2013), (6) Béthune et al. (2016), (7) Krapp et al. (2018), (8) Bai (2015), (9) Riols & Lesur (2018), (10) Suriano et al. (2018, 2019). 
In 3D unstratified simulations, zonal structures seem absent in the ideal or ambipolar cases but arise in the Halldominated regime via an antidiffusive process (Kunz & Lesur 2013; Béthune et al. 2016; Krapp et al. 2018). Early global simulations without explicit diffusion and neglecting the vertical gravity (Hawley 2001; Steinacker & Papaloizou 2002) have yet reported gaps and rings, whose origin were originally attributed to a “viscoustype” instability. However, the instability criterion is generally not fulfilled in simulations (see Sect. 3.7), and these structures may be created through a boundary effect (perhaps artificial) relying on the pileup of matter at the disc inner edge. Our own local unstratified simulations in the ideal limit or with ambipolar diffusion (not shown here) indicate that selforganised structures vanish as the horizontal box size is extended. Such behaviour is not obtained in stratified simulations, for which the separation of rings converges with the radial box size (see Sect. 5).
Stratified simulations (2D or 3D) show a radically different result:, ambipolar diffusion (with or without ohmic diffusion, see Simon & Armitage 2014; Bai 2015; Riols & Lesur 2018; Suriano et al. 2018) and plasma combining the three nonideal effects (Ohmic, Hall, and ambipolar, see Bai 2015) favour the emergence of zonal flows. However, in stratified discs, it is still uncertain whether the Hall effect alone can trigger their formation (Kunz & Lesur 2013; Lesur et al. 2014). Although ambipolar diffusion is believed to enhance the process of rings formation, local and global simulations without explicit diffusion (Bai & Stone 2014; Suriano et al. 2019) or with pure ohmic resistivity (Suriano et al. 2017, see also Sect. 4) also exhibit an efficient production of largescale rings. This is particularly obvious in the 2D axisymmetric case when fields are independent of the azimuthal coordinate. Perhaps, the 3D ideal case is a matter of discussion, since the structures obtained are generally more difficult to identify and fill the box entirely in local simulations. Actually, we see in Sect. 5 that rings of intermediate size undeniably form in 3D, but are efficiently diffused by vigorous nonaxisymmetric MRI turbulence.
We note that axisymmetric structures have been also discovered in zero net flux simulations (Johansen et al. 2009; Simon et al. 2012) but are transient, and probably emerge from stochastic processes. If we put aside the Hall effect, whose role in rings formation is probably of different nature (Kunz & Lesur 2013), and seek for a generic mechanism, it comes naturally from Fig. 1 that the rings formation process is ideal in essence and requires a largescale poloidal field with a vertical stratification.
2.2. Characteristic features
We remind the reader of the generic characteristics of zonal flows found in various MHD simulations (Hall free regime). Figure 2 shows the typical flow and magnetic field topology around such structure, obtained in the 3D simulation of Riols & Lesur (2018) with ambipolar diffusion. The quantities are averaged in time and in the azimuthal direction, for a magnetisation μ = 10^{−3} and ambipolar coefficient Am = 1 in the midplane (see Riols & Lesur 2018, for further detail about the numerical setup). In the top panel we can identify a strong and coherent windy plume that emanates from the density minimum (gap). The structure is inclined and connected to a large scale roll in the midplane. The bottom panel shows that the poloidal magnetic field is concentrated radially within the gap, while outside the net vertical flux is close to zero. It follows that the magnetisation (defined as the inverse of the plasma betaparameter in the midplane)
with , is stronger in the gaps and therefore anticorrelated with the rings (see also Suriano et al. 2018, for instance). The poloidal field lines follow the plume and clearly drive the outflow in this region. The existence of such localised and inclined magnetic shell is reminiscent of other MHD simulations (ideal or not, see Bai & Stone 2014; Bai 2015) and their origin seems clearly connected to the formation of the ring. We note that the inclined outflow and the spontaneous breaking of vertical symmetry is found in various simulations, including global models (Fig. 25 of Béthune et al. 2017; Gressel et al. 2015). Thus they cannot be considered as glitches of the shearing box. We stress, however, that local models do not correctly reproduce the outflow behaviour at large distance from the disc midplane. Actually, if the star is located on the left in Fig. 2, the upper field lines must bend at some point, further out in the atmosphere (see Fig. 25 of Béthune et al. 2017).
Fig. 2.
Example of a ring and gap structure and its surrounding outflow topology in the 3D shearing box simulation of Riols & Lesur (2018) with ambipolar diffusion. Top panel: gas density (colourmap in log scale) and streamlines (cyan lines) in the poloidal plane (x and z are respectively the radial and vertical coordinates). Bottom panel: vertical magnetic field (colourmap) and poloidal magnetic field lines (in white) showing the inclined “plume” structure from where most of the mass is extracted. The thickness of the lines is proportional to the field intensity. 
2.3. The key role of wind plumes
Another important and recurrent feature of the ring structure is related to the radial velocity profile. Figure 3 shows that the averaged v_{x} within the midplane (z≲H) is opposed to the formation of the ring. In particular, it is anticorrelated with the azimuthal velocity perturbation (zonal flow) and phased with an angle of π/2 with respect to the density maximum. During the ring formation, the matter is thus radially dragged towards the gap in average. Physically, this is expected since the turbulent stress acts like an effective viscosity that tends to diffuse any density bulge in the disc. Therefore, the only way to produce a ring is through vertical transport of matter. We find that the strong wind plume identified in Fig. 2 indeed clear out the gas in the gap regions. In other words, the ring structure does not result from matter being radially concentrated but from material in the gaps being depleted. We note that a similar result is obtained for other magnetisations and in larger box simulations with more than one ring. A detailed mass budget is provided in Appendix D for such simulation. In the next section, we show that a winddriven instability is probably connected with the formation of these plumes, which organise the flow into radial substructures.
Fig. 3.
Radial velocity computed from the 3D simulation of Fig. 2. It is averaged in the azimuthal direction, in z within ±1.5H and during the first 50 orbits (∼300 Ω^{−1}), corresponding to the growth of the ring structure. 
3. A unified theory based on a wind instability
3.1. Naive picture of the instability
Here we set out the basic principles of a wind instability, that may lead to the growth of ring structures in accretion discs. To start with, consider a classical “α” disc with constant and uniform viscosity ν_{t} = αc_{s}H. Assume initially a small axisymmetric and radial perturbation of the surface density δΣ. In the local approximation, then the viscous stress generates a radial flow associated with angular momentum transport,
which acts to diffuse the initial density ring perturbation. This is a classical and trivial result of the viscous disc theory. In the absence of winds, the perturbation decays and the disc remains stable. If now we include a wind, associated with a largescale poloidal magnetic field that removes the mass and angular momentum, the outcome may be radically different and potentially lead to an instability. In Fig. 4 we sketch out the main dynamical ingredients of such instability. If the magnetic diffusivity is not too high, the radial flow generated by the viscous stress concentrates the vertical magnetic flux B_{z} in the density minima or gaps (green patches). Because matter and magnetic flux are transported radially at the same speed δv_{R}, we expect the magnetisation to increase within the density minima. In such circumstances, the wind that was originally uniformly distributed in R, will adjust to this new configuration. A key hypothesis is that the vertical mass flux of the wind (proportional to the surface density) increases with μ. Therefore, if the magnetised gaps eject more material than viscosity can bring in, the initial perturbation is reinforced.
Fig. 4.
Sketch illustrating the development of rings in accretion discs. The black and blue arrows represent the radial flow and the wind, respectively. The green patches correspond to density minima where the vertical field B_{z} and the magnetisation μ grow. In this representation, the poloidal field lines (// to the wind) do not bend in the midplane but cross the disc with some angle. This configuration is actually not intuitive but is encountered in various disc simulations around the midplane. 
Quantitatively, an instability occurs whenever the surplus of mass launched in the wind is larger than the mass flowing radially. If we note r_{a} the typical separation between the rings and , the instability criterion becomes
Using relation (1) and the transport equation for B_{z}, it is possible to show that the instability criterion is simply p > 0. In sum, the instability is driven by the combination of an outflow and a radial flux transport, and appears optimal in the ideal limit. Rings do not originate from a radial concentration of matter but result from the vertical depletion of their surrounding gaps. The instability is of same nature as Lubow et al. (1994) and Cao & Spruit (2002) but does not require a wind torque (this point is discussed in Sect. 6). We note that the radial pressure gradient associated with the rings has to be balanced by the Coriolis force. The geostrophic equilibrium naturally gives birth to a zonal flow with v_{ϕ} > v_{K} in the regions where v_{R} < 0 and v_{ϕ} < v_{K} where v_{R} > 0 (see Fig. 4).
Although the instability mechanism relies on simple physical arguments, it needs to be demonstrated rigorously through a linear analysis of the MHD equations. Moreover there are many caveats to the simple picture described here: first the viscosity or transport coefficient α generally depends on μ, which is a widely accepted result based on turbulent MHD simulations. This effect may reduce the strength of the instability. Second, it is unclear what happens if angular momentum is free to flow along the poloidal magnetic field line. The presence of a mean toroidal field and a vertical stress could in principle have important consequences on the redistribution of angular momentum. Finally, we do not yet know how the instability behaves if the disc is subject to magnetic diffusion and nonideal effects. In the next section, we derive a general instability criterion taking into account several of these effects.
3.2. Averaged equations in the local framework
To simplify the problem, we have used the local shearing sheet framework (Goldreich & LyndenBell 1965). This corresponds to a Cartesian patch of the disc, centred at r = R_{0}, where the Keplerian rotation is approximated locally by a linear shear flow plus a uniform rotation rate Ω = Ω e_{z}. We note (x, y, z) the radial, azimuthal and vertical directions respectively.
To analyse the radial disc structure, we azimuthally and vertically integrated the equations of motion and therefore neglected the vertical dependence of the flow. For that purpose, we introduced two average procedures in the plane (y, z): a standard average between −z_{d} and z_{d}, where z_{d} is some arbitrary altitude, and a massweighted average so that for any field φ, we have:
where ρ is the fluid density and the surface density. Each field can be decomposed into a sum of a mean component (depending on x only) and a fluctuation
We also note the difference between the field at z = z_{d} and z = −z_{d}. The compressible, inviscid and isothermal equations of motion in the horisontal plane, integrated azimuthally and vertically between −z_{d} and z_{d} are
where is the mass loss rate, the stress tensor integrated in the vertical direction and its boundary value. We assumed that mass is replenished locally at a constant rate (associated with the radial accretion in a global view, which is absent in the local shearing box framework).
These equations of motions are coupled with the induction equation that describes the evolution of the magnetic field B. We are particularly interested in the evolution of the vertical magnetic flux ⟨B_{z}⟩:
where ℰ_{y} = v_{z}B_{x} − v_{x}B_{z} is the toroidal electromotive force and η is assumed to be a constant and uniform ohmic resistivity. An additional constraint is the solenoidal condition which gives
3.3. Stress tensor and electromotive force
We first develop the terms related to the stress tensor and W_{ij}. In the limit of highly subsonic fluctuations, it is straightforward to show that the radial stress and the vertical stress in the radial momentum equations are negligible compared to thermal pressure. Thus we can assume that and W_{xz} ≃ 0. In the azimuthal momentum equation, however, the stress is comparable to others terms. By using the decomposition of Eq. (2), we have
where
can be identified respectively as a turbulent and laminar radial transport coefficient. The term related to the vertical stress W_{yz} can be written as
where is the turbulent vertical transport. Finally the electromotive force in Eq. (6) can also be decomposed into a laminar and a turbulent part, the latter being assumed to behave as an effective magnetic diffusivity η_{t}:
3.4. Power laws for mass loss rate and turbulent coefficients
To close the system of equations we need to relate the mass loss efficiency and the turbulent coefficient α_{ν}, α_{W}, η_{t} to the integrated disc quantities. It is reasonable to assume that these coefficients depend principally on the main dimensionless parameter of the disc, namely the vertical magnetisation μ. We suppose that this dependence can be captured by a simple power law:
and similar relations for α_{W} and η_{t}. μ_{eq}, ζ_{0} and α_{ν0} correspond to values of a given equilibrium (see next section). Numerical simulations are actually a suitable tool to probe and test these scaling laws. The dependence of α_{ν} on the magnetisation has been explored in various ideal simulations of MRI turbulence, with or without magnetic diffusion and thermal effects (Hawley et al. 1995; Simon et al. 2013; Salvesen et al. 2016; Scepi et al. 2018). In all cases, there is some consensus that
The vertical transport due to a wind has also been measured in simulations (Bai & Stone 2013; Fromang et al. 2013; Zhu & Stone 2018; Scepi et al. 2018) but is generally supplied by a large scale toroidal field instead of turbulent fields. Properly characterising α_{W} would require us to measure the laminar and turbulent contribution of the stress separately, which has never been done in the literature. The turbulent diffusivity has been calculated in 3D unstratified shearing box simulations (Guan & Gammie 2009; Fromang & Stone 2009; Lesur & Longaretti 2009) and a fair assumption is to consider η_{t} between 0.2 and 0.5ν_{t} (where ν_{t} = α_{ν}c_{s}Ω). However, we see in Sect. 4 that such a ratio can actually be much weaker in 2D. Finally, evaluating the mass loss rate is probably the hardest part since it depends on the vertical box size and the nature of boundary conditions (Fromang et al. 2013). Simulations of 1D laminar winds predict p ≈ 0.6−0.7 (see Fig. 4 of Riols et al. 2016) while full 3D turbulent simulations suggest p ≃ 1; see Fig. 5 of Suzuki & Inutsuka (2009) and Fig. 4 of Scepi et al. (2018).
3.5. Local equilibrium solutions
We use a subscript “0” to indicate the equilibrium solutions (independent of time and x) of the system of equations derived in Sect. 3.2. The local disc equilibria are obtained by setting ∂⋅/∂t = 0 and ∂⋅/∂x = 0 in Eqs. (3)–(6). We note that in absence of turbulence, these equilibria correspond to the vertical average of the 1D (zdependent) wind solutions studied by Lesur et al. (2013) and Riols et al. (2016). The solenoidal condition gives immediately δB_{z} = 0, which means that . We can then define a constant magnetisation of the equilibrium, as
Solutions can be either symmetric about the midplane with B_{x} = B_{y} = 0 at z = 0 or antisymmetric with ∂_{z}B_{x} = ∂_{z}B_{y} = 0.

In the symmetric case, we have ⟨B_{x}⟩ = ⟨B_{y}⟩ = 0 but . Horizontal components satisfy but . Therefore there is a net vertical stress through the disc, which is responsible for a mean accreting flow.

In the second case ⟨B_{x}⟩ = ⟨B_{y}⟩≠0 but . We define B_{x0} and B_{y0} the mean radial and toroidal magnetic field throughout the midplane. Departures to the vertical average satisfy but . It is straightforward to check that W_{yz0} = 0 when such symmetry is enforced.
Historically, the first class of solutions (1) were considered to be the most intuitive and representative of a disc structure; the reason being that magnetic field lines do not bend outside the midplane. However antisymmetric solutions have been shown to naturally emerge from turbulent shearing box simulations (Lesur et al. 2014; Bai 2015), and even from global simulations when nonideal effects are included (Béthune et al. 2017; Bai 2017). For that reason, we focus mainly on the antisymmetric solutions in the following sections.
3.6. Linearisation around equilibrium
We use the subscript “0” to indicate equilibrium solution of Eqs. (3)–(6) enforcing the second class of symmetry. We remind the reader that such symmetry implies that at equilibrium. To study the stability of these solutions, we introduce small normalised axisymmetric perturbations of the form around the equilibrium flow:
We can perform a similar decomposition for the magnetisation
as well as the turbulent coefficients and mass loss rate whose first order perturbation is proportional to . To be consistent with the formalism of Sect. 3.2, the radial scale of perturbations is supposed to be much larger than the turbulent lengthscale. Because turbulent eddies are, in principle, limited to the disc scale height, we have considered modes with k_{x}H ≲ 1 and assume σ ≪ Ω. To simplify the problem, we made the following further assumptions:

Perturbations are in a geostrophic equilibrium (pressure gradient balances Coriolis force in the x direction). This is satisfied if σ ≪ Ω.

The vertical component of the Reynolds stress tensor as well as the turbulent vertical stress α_{W} = 0 are neglected

The stress is purely turbulent α_{ν0} ≫ α_{Ł0}

In Eq. (6), the vertical stretching of radial field ⟨u_{z}⟩⟨B_{x}⟩ and the vertical diffusion of B_{z} are neglected.

We have assumed that . This is true if the radial velocity perturbation does not vary too much between z_{d} and z_{d}.
Some of these assumptions are tested in simulations (see Appendices C and D). A more general case, including a laminar stress and a departure from the geostrophic equilibrium due to a strong toroidal field, is treated in Appendix B. By using these hypothesis, we show that the linearised momentum transport due to the total stress at leading order simply reduced to
We note that the perturbation of vertical stress cancels out with the term – in the azimuthal momentum equation (see Appendix A). We normalised k_{x} to H and σ to Ω and define η^{⋆} = (η + η_{0})/(H^{2}Ω). Linearisation of Eqs. (3)–(6) around equilibrium leads to
3.7. Stability criterion and growth rates
The linear system of Eqs. (9)–(12) can be simplified and cast into a (3 × 3) matrix problem whose determinant is
Solutions of the problem are found by setting D = 0. It is straightforward to show that growth rates follow the dispersion relation
with
There are actually two different regimes, depending on the strength of the turbulent stress. The first one corresponds to B < 0, i.e α_{ν0} ≪ ζ_{0} and ). In that regime, solutions are always unstable (ℜ(σ) > 0). The radial flow that concentrates magnetic flux is directly produced by the inertial term , via the geostrophic equilibrium and the conservation of angular momentum. However, this particular regime, though quite exotic and allowing the instability for α = 0, is never encountered in simulations. The second regime corresponds to B > 0 or α_{ν0} ≫ ζ_{0}. In that case, the instability occurs if C < 0, which corresponds in the ideal limit (η^{⋆} = 0), to
If p < 1 and q < 1, the nonideal term in C is always positive. Therefore it contributes to weaken the growth rate. We note that the hypothetical case q > 1 and p = ζ_{0} = 0 (no winds) can in principle lead to an instability in presence of finite resistivity but is excluded by the simulations. It corresponds to the usual, but hypothetical “viscous” instability, like imagined by Lightman & Eardley (1974).
In the ideal limit (η^{⋆} = 0), it is straightforward to show that the optimal growth rate is obtained in the limit k_{x} ≫ 1. Under the condition α_{ν0} ≫ ζ_{0}, which is verified in most of simulations, it can be demonstrated that the optimal growth rate, at leading order, is independent of α_{ν0} and proportional to ζ_{0}:
For typical MRIdriven turbulence with μ_{eq} = 10^{−3}, realistic values of the parameters are α_{ν} ≃ 0.3, p ≃ 1, q ≃ 0.5 and ζ_{0} ≃ 0.01. Using these values, we find in the ideal limit σ = 0.00666. To illustrate the effect of magnetic diffusivity, we plot in Fig. 5 the growth rate of the instability σ as a function of k_{x} for three different resistivities (η^{⋆} ≃ 0.002, 0.01, 0.06). These values are here purely adhoc and just to illustrate the dependence of σ(k_{x}) on η^{⋆}. Realistic values, directly inferred from 2D and 3D simulations, are used in Sects. 4 and 5. In all cases, the maximum growth rate is obtained for k_{x}H ≲ 0.5H which correspond to rings of size r_{a} ≳ 10 H. The instability occurs on long timescale 1/σ larger than 150 Ω^{−1} (20 orbits).
Fig. 5.
Growth rate of the winddriven instability as a function of the radial wavenumber for different magnetic diffusivities. Parameters of the model have been chosen in agreement with numerical simulations for fiducial μ_{eq} = 10^{−3}: α_{ν} ≃ 0.3, p ≃ 1, q ≃ 0.5, ζ_{0} ≃ 0.01. 
3.8. Eigenmodes
The linearised system admits solutions of the form:
with a small number compared to 1. In the turbulent regime with α_{ν} ≫ ζ, and if p > q, it is straightforward to show that E > 0. The radial velocity perturbation is out of phase (with an angle of π/2) with respect to the density maximum and anticorrelated with the zonal flow . The vertical field and the magnetisation are anticorrelated with the density maximum. This configuration is exactly that depicted in Sect. 2.
4. 2D axisymmetric simulations
To test our instability model, we needed to confront the theoretical predictions of Sect. 3 with MHD numerical simulations exhibiting rings structures. As a preliminary check, we performed in Appendix C a stability analysis around an initial laminar wind equilibrium, using 2D shearing box simulations. We find that axisymmetric perturbations undergo clean exponential growth with rates and eigenfunctions compatible with our theory. In this section, we explore the case of a turbulent disc, by using 2D simulations initialised with a net vertical field and random noise. In particular, we check whether the numerical growth rate and spacing of axisymmetric modes, as well as their physical behaviour, are compatible with the linear theory.
4.1. Numerical setup
Shearingbox simulations were run with the PLUTO code (Mignone et al. 2007), a finitevolume method with a Godunov scheme that integrates the compressible MHD equations in their conservative form. The fluxes are computed with the HLLD Riemann solver for runs without ambipolar diffusion and with HLL otherwise. The gas is isothermal and inviscid (no explicit viscosity). Boundary conditions are shearperiodic in x and periodic in y. In the vertical direction, we used standard outflow boundary conditions for the velocity field and imposed hydrostatic balance in the ghost cells for the density. In this way, we significantly reduced the excitation of artificial waves near the boundary. For the magnetic field, we used the “vertical field” or open boundary conditions with B_{x} = 0 and B_{y} = 0 at z = ±L_{z}/2. Because the instability we seek occurs on a long timescale, it is important to maintain a constant disc surface density Σ. Otherwise, the wind would empty the disc before the instability reaches any saturated state. For that purpose, we regularly injected mass near the midplane at each numerical time step. The source term in the mass conservation equation is similar to the one prescribed by Lesur et al. (2013),
where is the mass injection rate computed at each time step to maintain the total mass constant in the box, and z_{i} is a free parameter that corresponds to the altitude below which mass is replenished. We note that σ_{i} is uniformly distributed in x and y and thus independent of local density variations in the box.
For most of the simulations, we chose a large horizontal box size L_{x} = L_{y} = 20 H to be able to capture the largest rings. In z, the box spans −4 H to 4 H. We adopted a resolution of 256 points in the horizontal directions and 128 points in the vertical direction. In the ideal regime, such a resolution is insufficient to resolve the smallscale MRI turbulence properly, but enough to obtain the right axisymmetric dynamics. We checked that doubling resolution (512 × 256) does not change the results in terms of growth rate and spacing.
Finally, for simulations with ambipolar diffusion, we use exactly the same setup as Riols & Lesur (2018) where the ambipolar Elsasser number Am is one in the midplane, and increases abruptly above a certain height corresponding to the FUV ionisation layer (see Sect. 2.4 of the paper for further detail of the prescription and its physical motivation).
4.2. Simulations in the ideal limit and numerical growth rates
We first run a series of 2D axisymmetric simulations in the ideal limit (without explicit resistivity or ambipolar diffusion) by varying the vertical magnetisation μ_{eq} from 10^{−5} to 0.03. All runs are initialised with a hydrostatic equilibrium in density and a weak random noise in velocity. Figure 6 shows the evolution of the turbulent and laminar transport coefficients α_{ν}, α_{L} and the wind loss efficiency ζ for the case μ_{eq} = 10^{−3}. To complement the analysis, the evolution of column density Σ(x) is illustrated in Fig. 7 (top panels) for three different magnetisation μ_{eq} = 10^{−4}, 10^{−3} and 10^{−2}. Each quantity is integrated between −1.5 and 1.5 H. In all cases, we identify three distinct phases associated with: (a) the development of vigorous MHD turbulence and the launching of a wind during the first ∼50−150 Ω^{−1}, (b) the growth of zonal flows and density rings, and c) the saturation of zonal flows which is accompanied by a severe drop in the turbulent stress and a modest drop in the wind loss efficiency ζ (see Fig. 6). We note that during the initial turbulent phase, α_{ν} ≫ α_{L}, while after the saturation of substructures, the radial stress is predominantly laminar α_{ν} ≪ α_{L}. Figure 7 shows that the timescale associated with rings formation seems to decrease with μ_{eq} while their spacing and strength seems to increase with μ_{eq}.
Fig. 6.
Evolution of the transport coefficient α_{ν}, α_{L} and the wind mass loss efficiency ζ in the 2D simulation without explicit diffusion (μ_{eq} = 10^{−3}) 
Fig. 7.
Top: spacetime diagram showing the column density as a function of time and x, in 2D axisymmetric turbulent simulations, for different magnetisation (left to right, μ_{eq} = 10^{−4}, 10^{−3} and 10^{−2}). The density is integrated within z ± 1.5H and runs were computed in the ideal limit (without explicit magnetic diffusion). Bottom panels: for each run, timeevolution of , normalised projection of B_{z} onto the prominent axisymmetric Fourier mode (with largest amplitude). This mode corresponds respectively to k_{x} = 6, 4 and 2k_{x0} for μ_{eq} = 10^{−4}, 10^{−3} and 10^{−2}, with k_{x0} = 2π/L_{x} the fundamental box radial wavenumber. 
To understand whether rings form via a linear instability or a more complicated nonlinear process, we explored the dynamics of axisymmetric modes in Fourier space. The procedure is simple: at each timestep, we perform a 1D FFT (along the x direction) of the vertically averaged quantities. We then project each field (or any physical quantity) onto the dominant axisymmetric mode that grows during the simulation^{1}. In this way, we keep the interesting dynamics related to the prominent ring and filter out part of the dynamics associated with the turbulent flow. For a given field φ, we note such projection, normalised with respect to its mean (the k_{x} = 0 mode in Fourier space).
The bottom panels of Fig. 7 show the evolution of for three different magnetisations. In all cases, the prominent ring mode begins to grow quasiexponentially, indicating that a linear instability is at work. We note that the projection procedure is necessary to obtain a clean exponential growth at the beginning of simulations. After a few tens of orbits, the growth stops and saturates around one. This is an indication that a nonlinear regime is reached.
During the linear phase, we measure the growth rates associated with the prominent mode and report them in Fig. 8 for different μ_{eq}. For μ_{eq} > 10^{−4}, growth rates increase with the magnetisation and vary from 0.002 to 0.05 Ω; the dependence of σ with μ_{eq} follows a power law with index 0.7−0.8. We checked that σ is not changed when doubling the resolution of simulations. It is worth noting that such growth rates are ∼10−20 smaller than those characterising the MRI at a similar scale.
Fig. 8.
Optimal growth rate of the instability as a function of the magnetisation μ. The red diamond markers are the growth rates of the most unstable Fourier mode measured in 2D ideal simulations. The orange and green crosses are those measured in resistive and ambipolar runs. The purple solid line correspond to the model prediction, using the scaling relations (17) and (18). The blue and yellow dashed lines use the same model but with constant η^{⋆} = 0.01 in both cases and constant α_{ν} = 0.1 for the last case. 
4.3. Nonideal case
We first investigated the effect of ohmic resistivity on the growth of axisymmetric structures. For that, we ran a 2D simulation with μ_{eq} = 10^{−3} and explicit η/(ΩH^{2}) = 0.01. This corresponds to a magnetic Reynolds number Rm = ΩH^{2}/η = 100. During the first tens of orbits, a turbulent state develops but with a much weaker strength and transport than in the ideal case (for comparison, α ≃ 0.05). Such a result is expected since the MRI is quenched by the resistivity. However, the wind mass loss efficiency ζ is of the same order of magnitude. Three rings, associated with a mode k_{x} = 3k_{x0} develop in the box and grow at a rate σ ≃ 0.011Ω. The number of rings, their properties, and the growth rate are actually comparable to those obtained in the ideal case.
We then studied the effect of ambipolar diffusion by running a simulation with the same μ_{eq} = 10^{−3} and Am = 1 in the midplane. Figure 9 shows the evolution of surface density and verticallyaveraged B_{z}. Again three or four structures seem to emerge from the initial turbulent phase. The growth rate associated with the k_{x} = 4k_{x0} mode in Fourier space is σ ≃ 0.0085, very similar to the value obtained in the ideal limit. Thus, ambipolar diffusion does not seem to alter the instability mechanism, although it weakens or even suppresses the initial turbulent state.
Fig. 9.
Top panel: spacetime diagrams showing the surface density Σ(x, t) in the 2D ambipolar simulation (Am_{mid} = 1, μ_{eq} = 10^{−3} and z_{d} = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially. Bottom panel: evolution of the B_{z} component projected onto the k_{x} = 4k_{x0} mode. 
We note that contrary to the ideal case, the turbulent transport is either comparable to or smaller than the laminar transport α_{L}. To go further in the analysis and check that the instability identified numerically is of same nature as that described in Sect. 3, we inspected in the ambipolar run the different flux and source terms in Eqs. (3)–(6). For ease of reading, this analysis is presented in Appendix D. We show in particular that the ambipolar term η_{A}J × (B × B)/B^{2} merely acts as a diffusion on the B_{z} field.
4.4. Parameters p and q and confrontation with the model
In this section, we describe our investigation of whether the dependence of numerical growth rates σ(μ_{eq}) can be predicted by the simple model exposed in Sect. 3. We have focussed particularly on the ideal simulations, for which α_{L0} ≪ α_{ν0} during the linear phase. The linear theory can actually be generalised to quasilaminar discs with α_{L0} ≳ α_{ν0} (see Appendix B), typically those obtained in nonideal simulations. The instability in this regime is conceptually not different and the theoretical growth rates are comparable to those obtained with a pure turbulent stress.
The model of Sect. 3 depends on three main parameters p = d log(α_{ν})/d log(μ), q = d log(ζ)/d log(μ) and the diffusivity η_{⋆}. To evaluate these coefficients, it is suitable to measure them directly from numerical simulations. A naive way is to run a series of turbulent simulations varying μ and measuring α_{ν}, ζ and η_{⋆}. Though simple, this method is complicated to accomplish in practise. Indeed, as suggested by Fig. 6, the quasisteady turbulent state obtained in the early stage of simulations is short (< 100 Ω^{−1}) and rapidly affected by the zonal structures. This turns to be particularly critical at large μ_{eq}, for which measures of transport coefficient and mass loss rate are not statistically meaningful and can be highly inaccurate.
A possible way to circumvent this issue is to infer directly the scaling relations from a thorough inspection of the largescale ring perturbations themselves. The idea is to compute the perturbed coefficients and associated with the dominant axisymmetric mode and see how it correlates in time with the perturbed magnetisation . In addition to giving p and q, it provides a great opportunity to check the linear relations conjectured in Sect. 3.6, which are fundamental requisites of the instability.
Figure 10 shows the two ratios and for the case μ_{eq} = 10^{−4}. Here , and are the vertically averaged perturbations of transport, mass loss efficiency and magnetisation projected onto the Fourier mode k_{x} = 6 × 2π/L_{x} and normalised with respect to the mean (k_{x} = 0 mode). We note that the latter is averaged in time during the growth phase. Although these ratios are highly fluctuating (this is particularly true for p), and seem linearly correlated to the magnetisation , at least statistically. Most importantly, we check that p > q. An interesting, but unexpected, result is that the ratios seem to keep a similar value in the saturation regime (t ≳ 500Ω^{−1}). By averaging over time, we find q ≃ 0.4 and p ≃ 0.9 in agreement with the scaling relations obtained in 3D fully turbulent simulations (Scepi et al. 2018). We performed the same calculation for μ_{eq} = 10^{−3} and found p ≃ 0.8 and q ≃ 0.55.
In summary, we adopted the following scaling laws in our model:
with the constants calibrated to fit with the simulation data at intermediate μ = 10^{−3} (values are taken from Fig. 6). Simulations in this regime generally provide a more accurate measure of α_{ν0} and ζ_{0} since the turbulent phase settles longer. The dependence of α is very close to that obtained in past 3D simulations (see Eq. (20) of Salvesen et al. 2016).
Fig. 10.
Coefficients p and q measured from the 2D ideal simulation with μ_{eq} = 10^{−4}. The red (respectively blue) line shows the ratio between the perturbation of ζ (respectively α) and the perturbation associated with magnetisation, as a function of time. Each quantity (ζ, α and μ) is projected onto the mode k_{x} = 6k_{x0}, averaged within z ± 2H and normalised to its mean. 
To estimate the turbulent magnetic diffusivity, we have used a slightly different method. Instead of projecting the flow into the Fourier space, we calculated the averaged term directly in real space during the linear growth of the mode. We checked that it correlates quite well in space with ∂⟨B_{z}⟩/∂x and indeed acts to diffuse the large scale structure. The ratio between the two terms is η_{t}/(ΩH^{2}) ≈ 0.003 for μ_{eq} = 10^{−4} and η_{t} ≈ 0.008 for μ_{eq} = 10^{−3}. This corresponds to rather low turbulent diffusivities, compared to values given by Lesur & Longaretti (2009) in the unstratified 3D case. The turbulent Prandlt number Pm = α_{ν}/η_{t} is of order 25 and we checked that it depends little on resolution. Therefore, we consider:
These scaling laws are then used as inputs of our model (Sect. 3.7), and allow us to compute the optimal growth rate as a function of μ (using Eq. (14)). The result is superimposed in Fig. 8 (purple curve). The model fits quite well with the numerical values obtained in Sect. 4.2, and in particular reproduces the slope of σ(μ) for intermediate μ. We stress that such a result does not rely exclusively on the linear correlations measured in Fig. 10, since growth rates depend also on the absolute values of ζ. To show that the model is robust, we also plot the theoretical growth rates obtained for a constant magnetic diffusivity (blue curve) and constant α and η^{⋆} (yellow dashed curve). In both cases, there are few differences with the prescription η_{t} = 0.04 α_{ν}. This is unsurprising since maximum growth rates depend little on α and η^{⋆} (for the range of diffusivities considered here, see Fig. 5).
We finally comment on the extreme cases, corresponding to the lowest and highest magnetisations in Fig. 8. In these regimes, the numerical values slightly deviate from the model. For μ ≳ 10^{−1.5}, measures are perhaps inaccurate since growth rates associated with largescale modes become comparable to those associated with the MRI phase. For μ ≲ 10^{−4}, the small discrepancy is probably due to a sudden drop in α and η_{t}. The fact that the turbulence dies out at low μ is expected since the 2D box does not sustain a dynamo in the limit μ → 0. As turbulent dissipation is weakened, zonal flows are slightly enhanced. This effect, however, remains marginal and only holds for 10^{−5} ≲ μ ≲ 10^{−4}. Below μ = 10^{−5}, the vertical flux of mass associated with the wind drops abruptly, and we checked numerically that zonal flows vanish.
4.5. Nonlinear saturation and ring/gap contrast
As suggested by Fig. 7, the ring instability saturates and enters a nonlinear regime, once the azimuthal structures reach a significant amplitude. During this phase, the structures stop growing but prevail in the flow and remain stable for hundreds of orbits. Their nonlinear saturation leads to a drop in the mean stress (Fig. 6) and the production of a strong mean azimuthal field B_{y}, altering significantly the initial hydrostatic density profile. But it is still not clear how this saturation occurs and what determines the final amplitude (or density contrast) of the rings. By examining the density contrast in the case μ = 10^{−4} or μ = 10^{−3} (upper panels of Fig. 7), we find that the gaps in the nonlinear regime still contain a large amount of gas. For instance, the surface density perturbation associated with the prominent Fourier mode, settles towards for μ_{eq} = 10^{−4} and for μ_{eq} = 10^{−3}. Therefore, saturation does not occur because the material in the gaps has been emptied. The lower panels of Fig. 7 show instead that the instability saturates when (i.e perturbation of B_{z} ∼ B_{z0}). In other words, the instability stops when there is no more vertical flux to drag in. The density contrast between the gaps and the rings in the nonlinear regime can be then estimated by setting in Eq. (15). We obtain in the ideal limit:
where the last equality is obtained in the limit . Figure 11 shows the theoretical density contrast for different magnetisations (plain lines) as a function of k_{x}. This is calculated using the same parameters and scaling relations as in Sect. 4.4. The result is that ΔΣ/Σ_{0} depends mainly on the radial separation, which is a function of the magnetisation. We report on the same figure the density contrast of leading axisymmetric modes (respectively k_{x} = 2, 4, 5, 5k_{x0} measured in simulations for μ_{eq} = 10^{−4}, 10^{−3.5}, 10^{−3} and 10^{−2.5}. Numerically, there is a close relationship between the density contrast and the rings separation, which appears consistent with the theoretical calculations.
Fig. 11.
Contrast between rigs and gaps versus radial separation. The plain curves are estimated from the theory (Eq. (19)) while the diamond markers are points measured from simulations. The blue, orange, green, and red colours correspond respectively to μ_{eq} = 10^{−4}, 10^{−3.5}, 10^{−3}, and 10^{−2.5}. Dotted lines delimit the region above which dust can be concentrated within the rings, assuming H/R = 0.05 (black) and H/R = 0.1 (purple). See Sect. 4.7 for further detail. 
4.6. Ring separation
So far, we have simply measured the ring separation as the most prominent k_{x} in simulations, but we also want to know if it is possible to predict this quantity from the linear theory. For large magnetisation (μ_{eq} = 10^{−2} and 10^{−2.5}), this corresponds roughly to the radial scale that maximises the theoretical growth rate. However, for smaller magnetisation, this is no more the case. For instance, the prominent modes in simulation with μ_{eq} = 10^{−4} corresponds to k_{x} = 5k_{x0} = 1.5 H^{−1} and 6k_{x0} = 1.8 H^{−1}, while the linear theory gives a maximum growth rate at k_{x} ≃ 0.45 H^{−1}. Actually, this is not in contradiction with the theory since σ appears to be quite flat with k_{x} in this regime. The reason is that turbulent dissipation is low (Pm_{t} = 25) and permits small scales axisymmetric modes to be amplified with significant growth rates. Modes then compete with each other and it becomes extremely difficult to predict the spacing from linear theory. The initial amplitude will be probably a decisive factor in the selection of the dominant mode(s). Such amplitude depends on the spectrum of the initial turbulence phase, which is a priori difficult to predict. We find (not shown here) that the spectrum of B_{z} associated with the initial turbulence peaks around k_{x} = 6k_{x0} and beyond for μ_{eq} = 10^{−4}, while it peaks at much lower k_{x} ≲ 0.3 in the case μ_{eq} = 10^{−2}. For μ_{eq} = 10^{−4}, there is a difference of roughly a factor of 25 in initial energy between the k_{x} = 6k_{x0} mode and the largescale mode which, in principles, should grow fastest (k_{x} ≃ 0.45 H^{−1}). Therefore, it is likely that the ring spacing in our simulations is partly ruled by the initial turbulent conditions. Whether the spectrum of this turbulence is universal or depends on the historic and detailed physics of the disc remains an open question. We note however that all simulations, even in 3D (see Sect. 5), indicate a similar trend that the ring spacing increases with magnetisation. We finally emphasise that it does not depend too much on horizontal box size (see a comparison with small box Riols & Lesur 2018) or resolution.
4.7. Criterion for dust concentration
With the result of Sect. 4.5, we can infer a minimum magnetisation μ for which dust can accumulate into the gaps. A simple criterion is that the local pressure gradient in the disc have to be positive:
where P_{0} is the mean equilibrium pressure profile and k_{R} the local radial wave number of the ring perturbation. A typical disc model, which matches disc observations (Andrews et al. 2009), has
The criterion for dust accumulation becomes:
The black and purple dotted lines in Fig. 11 delimit the region above which a radial concentration of dust is possible, respectively for H/R = 0.05 and H/R = 0.1. We predict that magnetisations greater than 10^{−4} or 10^{−3.5} will lead to the formation of dusty rings. A more accurate estimation would require to take into account the radial diffusion of particles by the turbulence, which depends on the grain size and the nature of nonideal effects, but this is beyond the scope of this paper.
5. 3D simulations
We extended our analysis to 3D simulations and show that ring structures exhibit properties that are again compatible with the instability model of Sect. 3.
5.1. Dependence on diffusive processes
First, we performed three different simulations for μ_{eq} = 10^{−3}, respectively with no explicit diffusion, ohmic resistivity η = 0.01 and ambipolar diffusion Am_{mid} = 1. These simulations are initialised with random perturbations and are computed with L_{x} = L_{y} = 20 H, L_{z} = 8 H and resolution N_{X} = N_{Y} = 256, N_{Z} = 128. The top panels of Fig. 12 show that in all cases, zonal flows are produced, with growth rates and size that can substantially vary from on case to another. To help the analysis, we show in the lower panels the evolution of the mean transport coefficients α_{ν}, α_{L}, and mean mass loss efficiency ζ. In the ideal limit (simulation without explicit diffusion), these quantities are almost identical to those obtained in 2D (see Fig. 6 for comparison). However, the ring structure is much wider and takes longer time to form. The growth rate, measured from the left panel of Fig. 13 is indeed σ_{3D} ≃ 0.0034 ≃ σ_{2D}/2.5. The main difference is that the turbulent magnetic diffusivity η_{t} in 3D is considerably enhanced by the azimuthal (or nonaxisymmetric) MRI dynamics. We indeed measure η ≃ 0.12 which corresponds to , in agreement with other 3D numeric simulations (Lesur & Longaretti 2009). Therefore, according to the linear theory, the growth rate is reduced and drops rapidly with k_{x} (unlike the 2D case for which it was quite flat.) The maximum growth rate σ ≃ 0.0039 is obtained for k_{x}H = 0.2 and stability is reached for k_{x}H ≳ 0.42. This explains why we obtain a single ring in the 3D case, while three or four rings were obtained in 2D. We note that, unlike the 2D case, α does not drop significantly after saturation of the instability and vigorous turbulent motions persist along with the zonal structure.
Fig. 12.
Top: spacetime diagram showing the column density as a function of time and x, in 3D turbulent simulations with μ_{eq} = 10^{−3}. Left panel: no explicit diffusion; centre panel: ohmic diffusion (Rm = ΩH^{2}/η = 100); right panel: ambipolar diffusion (Am_{mid} = 1). The density is integrated within z ± 1.5H. Bottom panels: for each run, timeevolution of turbulent α_{ν}, laminar α_{L} and ζ the mass loss efficiency. 
Fig. 13.
Left: timeevolution of , the normalised projection of B_{z} onto the prominent axisymmetric Fourier mode (k_{x} = 2π/L_{x}) in the ideal simulation with μ_{eq} = 10^{−3}. Right: growth rate of the linear instability as a function of the magnetisation μ. The red diamond markers are measures from 3D ideal simulations. The green and yellow crosses are those measured in resistive and ambipolar runs. The plain lines correspond to model predictions, using similar scaling relations as (17), but with η_{t} = 0.5α_{ν} in the ideal case. 
In nonideal simulations, Fig. 12 (bottom panels) indicates that the wind mass loss efficiency ζ ≃ 0.015 is comparable to that in the ideal case. However, the turbulent stress is weak or even reduced to zero in the ambipolar case. The radial transport of angular momentum is provided essentially by the laminar component of the stress −⟨B_{x}⟩⟨B_{y}⟩. Azimuthal structures are very faint and have little impact on the unstable dynamics. We did indeed check that turbulent diffusion of B_{z} structures is negligible compared to the explicit one. As a consequence, growth rates are much larger than in the ideal case. We measure respectively σ = 0.018 and σ = 0.021 in the ohmic and ambipolar simulations. Using an extension of the linear theory (see Appendix B) in the limit α_{L} ≫ α_{ν}, we find that the theoretical growth rates match with the numerical values.
We note that σ is twice as big as the 2D axisymmetric case. We attribute this difference to the fact that ζ is on average twice as large as in 2D. Finally, with ambipolar diffusion, the number of rings is identical between the 2D and 3D case. With pure ohmic diffusion, the 3D simulation exhibits a single ring (instead of 3 in 2D). It is possible that the stronger laminar stress in 3D (a factor ≳2) favours the emergence of larger scale structures, thought initial conditions may also play a secondary role in the final shape.
5.2. Dependence on magnetisation μ
Finally, we ran a set of 3D ideal simulations (without explicit diffusion) by varying μ. We scanned a large range of magnetisations from μ_{eq} = 10^{−4} to 10^{−1}. For the largest magnetisations ( μ = 10^{−2} and 0.08), we used a box of size L_{x} = 30 and 40 H respectively. Two or three rings are obtained for μ = 10^{−4} while a single ring is obtained for intermediate magnetisations. In the extreme case (μ ∼ 0.08), three or four rings form at the very early phase, but rapidly merge into a single ring. The growth rates of each prominent modes are shown in the right panel of Fig. 13 (red diamonds). We superimposed the theoretical growth rate σ(μ) (purple line) obtained by using the linear theory of Sect. 3. Here we assumed the same scaling relations for α and ζ as in 2D (Eq. (17)) but with η_{t} = 0.5α_{ν}. The theory reproduces quite well the numerical values, except maybe in the large magnetisation regime (μ ∼ 0.08). In the same figure, we also plot the numerical growth rates obtained in the resistive run (yellow marker) and those measured in ambipolar runs (green markers). The green dashed line accounts for the growth rate computed with an extension of the linear theory, in the regime α_{L} ≫ α_{ν} (see Appendix B) assuming ζ = 7 μ^{0.85} and α_{ν} = 0; α_{L} = 4μ^{0.45}.
6. Discussion and conclusion
In summary, we showed that the ring and gap structures obtained in simulations of magnetised discs are formed via a linear instability. The process is generic and works in all geometrical configurations (2D or 3D), with or without turbulence, and in various diffusive regimes (ideal, ohmic or ambipolar). This instability is driven by a magnetic wind, associated with a largescale poloidal field, and relies on the assumption that the mass loss rate increases locally with the magnetisation. The process works as follows:

A small radial perturbation of density generates a radial flow directed towards the gaps, via the turbulent stress if α ≠ 0.

The magnetic field, initially uniform, is radially transported towards the gaps.

The excess of poloidal flux induces a stronger wind and an excess of ejected mass in the gaps.

The initial density perturbation is thus reinforced.
We showed that in theory, the instability exists in the “nostress” regime (α = 0). In that case, the radial flow results simply from angular momentum conservation. Indeed, assuming a geostrophic balance and an initial density perturbation, an azimuthal or zonal flow is produced. Therefore, the inertial term in the azimuthal direction produces a radial flow opposed to the initial density ring. However in practice the radial transport of matter and magnetic field is mainly (or even totally) induced by the viscous stress. Growth rates are typically of the order of the mass loss rate efficiency ζ = (ρu_{z})_{∣z = H}/(ΣΩ^{−1}). In Sects. 4 and 5, we brought several evidence that such linear instability is at work in our MHD simulations. We showed in particular that axisymmetric modes undergo exponential growth in the early phase of the simulations, at rates compatible with the linear theory. A remarkable result is that growth rates increase rapidly with μ, and can be of order ≃0.1 Ω for μ = 0.1. This is a direct consequence of the wind reinforcement at large μ. Although the instability can theoretically exist in the limit of small wavenumbers, the rings separation is set by nonideal processes in the disc and is typically ≃10 H for realistic diffusion coefficients. The final ring density contrast can be also estimated from a simple saturation predictor. Our model suggests that in typical TTauri discs with Σ ∝ R^{−1}, the dust accumulates into the gas rings for magnetisation μ ≳ 10^{−4}. We finally provided detailed diagnostics of the axisymmetric flow, which are all consistent with the instability mechanism described above.
Unlike the mechanism suggested by Lubow et al. (1994) and Cao & Spruit (2002), the magnetic torque exerted by the wind is unnecessary here. We note that our linear analysis has been carried out around antisymmetric equilibria, for which the mean vertical stress and accretion flow are zero. Given its robustness, we think that the inclusion of a mean wind torque in our model will not dramatically affect the process, though an extension of the linear theory needs to be developed in this configuration. One possible effect is that the ring structures will be advected by the mean accretion flow and therefore slowly drift towards the central object. Global simulations with ambipolar diffusion have shown that ring structures persist despite the recurrent change in the largescale wind geometry and vertical symmetries during the disc evolution (Béthune et al. 2017; Suriano et al. 2019). Further investigations will nevertheless be necessary to characterise the instability mechanism in the global configuration. Another caveat is related to the effect of boundary conditions and mass replenishment. Indeed, the mass loss rate is known to depend on the location of the vertical boundary (Fromang et al. 2013), as the outflow never crosses the fast magnetosonic point. The growth rate is then impacted too, according to our calculations. We then expect that the dependence of growth rates on μ is still correct, but that the absolute values of σ might depend on the external environment of the disc, but also on the topology of magnetic fields and wind launching conditions, which in reality can differ from the simulations and from object to objects. Regarding the mass replenishment procedure, it takes place on a long timescale, similar to the instability timescale. However, we think that this does not affect the whole process since it is independent on the radial coordinate. We also conducted simulations without replenishment and found that rings still form within the same timescale, despite the progressive loss of mass experienced by the disc. This indicates that mass replenishment does not play any role in the ring formation mechanism.
We think that the process described in this paper could have important implications for discs around young stars. Observationally, there are indications that some discs may emit a wind through the action of a large scale poloidal field. Indeed recent studies have reported lowvelocity molecular outflow consistent with MHD winds in Class 0 and I discs (Launhardt et al. 2009; Bjerkeli et al. 2016; Tabone et al. 2017; Hirota et al. 2017; Louvet et al. 2018), although alternative processes are not excluded (photoevaporation winds or jet “cocoon”). In theory, the magnetic flux from the parental core during the disc formation could be sufficient to provide the degree of disc magnetisation required for the launching of a wind. However an important unknown is how the poloidal magnetic flux is radially transported during the disc lifetime. If such transport is efficient, as is suggested by recent works (Guilet & Ogilvie 2013; Zhu & Stone 2018), the wind instability process described in this paper may be possible during the early stages of the disc evolution (Class 0 to Class II). As shown in Sect. 4, the instability form longlived rings of size ≫H, which might remain stable during a large fraction of the disc lifetime. Such winddriven instability could potentially explain the gaps imaged by ALMA or SPHERE in young discs like HL Tau or eventually leads to cavities like those observed in transitional discs (an alternative is the transonic winddriven accretion recently proposed by Wang & Goodman 2017).
Recent MHD simulations with similar setup (Riols & Lesur 2018) have shown that the pressure maxima associated with the rings can efficiently trap dust grains of intermediate size (mm to dm). This could be directly relevant for planet formation theory, as it may help to overcome the “radialdrift” barrier (growth of grains to pebbles) and the “fragmentation” barrier (growth from planetesimals to planets) during the early phases when dust abounds (Birnstiel et al. 2016). Indeed, the confinement of solids in pressure maxima can stop the radial migration of grains, accelerate their growth, and prevent their fragmentation (Gonzalez et al. 2017). Further work is needed to understand how fast grain can grow and how the dust backreact on the gas rings. Another interesting avenue of research is to investigate the interaction between these rings and the disc gravitational instability, which is expected in class 0 and class I discs (Tobin et al. 2013; Mann et al. 2015). It is indeed unclear whether selfgravity could enhance or destroy the zonal MHD flows.
Finally, the linear instability identified in this paper could be applied to other accreting systems, such as active galactic nuclei, dwarf novae or Xray binaries. Indeed some of these objects show signatures of strong jets, potentially driven by large scale poloidal fields via a magnetocentrifugal effect. The spontaneous formation of rings via a linear instability could have important consequences in their dynamics and their variability.
This is defined as the mode with maximum amplitude, which corresponds to the k_{x} = 4k_{x0} = 8π/L_{x} mode at μ_{eq} = 10^{−3}. At low magnetisation (μ = 10^{−4}) however zonal flow cannot be properly reduced to a single mode. The peak of the spectrum oscillate between the k_{x} = 5 and k_{x} = 6k_{x0} components.
Acknowledgments
This work acknowledges funding from The French ANR under contracts ANR17ERC20007 (MHDiscs). This work was granted access to the HPC resources of IDRIS under the allocation A0040402231 made by GENCI (Grand Equipment National de Calcul Intensif). Part of this work was performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujfgrenoble.fr), which is supported by the RhôneAlpes region (GRANT CPER0713 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
References
 ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N. 2015, ApJ, 798, 84 [Google Scholar]
 Bai, X.N. 2017, ApJ, 845, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2014, ApJ, 796, 31 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, Protostars and Planets VI, 667 [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2016, A&A, 589, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Bjerkeli, P., van der Wiel, M. H. D., Harsono, D., Ramsey, J. P., & Jørgensen, J. K. 2016, Nature, 540, 406 [NASA ADS] [CrossRef] [Google Scholar]
 Cao, X., & Spruit, H. C. 2002, A&A, 385, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Penzlin, A. B. T. 2018, A&A, 609, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., & Stone, J. M. 2009, A&A, 507, 19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garufi, A., Benisty, M., Stolker, T., et al. 2017, The Messenger, 169, 32 [NASA ADS] [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, J.F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
 Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Guan, X., & Gammie, C. F. 2009, ApJ, 697, 1901 [NASA ADS] [CrossRef] [Google Scholar]
 Guilet, J., & Ogilvie, G. I. 2013, MNRAS, 430, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F. 2001, ApJ, 554, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Helled, R., & Bodenheimer, P. 2014, ApJ, 789, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Hirota, T., Machida, M. N., Matsushita, Y., et al. 2017, Nat. Astron., 1, 0146 [Google Scholar]
 Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Konigl, A., & Wardle, M. 1996, MNRAS, 279, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Krapp, L., Gressel, O., BenítezLlambay, P., et al. 2018, ApJ, 865, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 [NASA ADS] [CrossRef] [Google Scholar]
 Launhardt, R., Pavlyuchenkov, Y., Gueth, F., et al. 2009, A&A, 494, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesur, G., & Longaretti, P.Y. 2009, A&A, 504, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 268, 1010 [NASA ADS] [CrossRef] [Google Scholar]
 Mann, R. K., Andrews, S. M., Eisner, J. A., et al. 2015, ApJ, 802, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Moll, R. 2012, A&A, 548, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Okuzumi, S., Momose, M., Sirono, S.I., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riols, A., & Lesur, G. 2018, A&A, 617, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riols, A., Ogilvie, G. I., Latter, H., & Ross, J. P. 2016, MNRAS, 463, 3096 [NASA ADS] [CrossRef] [Google Scholar]
 Salvesen, G., Armitage, P. J., Simon, J. B., & Begelman, M. C. 2016, MNRAS, 460, 3488 [NASA ADS] [CrossRef] [Google Scholar]
 Scepi, N., Lesur, G., Dubus, G., & Flock, M. 2018, A&A, 620, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, J. B., & Armitage, P. J. 2014, ApJ, 784, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Bai, X.N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013, ApJ, 775, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Steinacker, A., & Papaloizou, J. C. B. 2002, ApJ, 571, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Suriano, S. S., Li, Z.Y., Krasnopolsky, R., & Shang, H. 2017, MNRAS, 468, 3850 [NASA ADS] [CrossRef] [Google Scholar]
 Suriano, S. S., Li, Z.Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239 [NASA ADS] [CrossRef] [Google Scholar]
 Suriano, S. S., Li, Z.Y., Krasnopolsky, R., Suzuki, T. K., & Shang, H. 2019, MNRAS, 484, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Suzuki, T. K., & Inutsuka, S.I. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Takahashi, S. Z., & Inutsuka, S.I. 2014, ApJ, 794, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Tobin, J. J., Chandler, C. J., Wilner, D. J., et al. 2013, ApJ, 779, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., & Goodman, J. J. 2017, ApJ, 835, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Wünsch, R., Klahr, H., & Różyczka, M. 2005, MNRAS, 362, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Linearisation of total stress in the azimuthal momentum equation
In this appendix, we describe how we linearised the total stress in the ymomentum equation (radial plus vertical) around a given disc equilibrium. For that we used assumptions (2) and (3) of Sect. 3.6. These assumptions are checked numerically for unstable modes around a laminar disc equilibrium (see Appendix C). Using the formalism of Sect. 3 and Eqs. (7) and (8), the linearised quantities are
Using the solenoidal condition, integrated azimuthally and vertically, we obtain that . Therefore, the term associated with the vertical stress perturbation compensates exactly the term in the laminar radial stress. Adding W_{yz} and the x derivative of the averaged radial stress gives:
with . We needed then to relate the perturbation to the variables of the reduced model or û. For that, we assumed that vertically, the disc is in a magnetohydrostatic equilibrium and the magnetic field is essentially toroidal:
This is true below the Alfven point, where inertial terms of the wind are small compared to magnetic pressure. If we integrate Eq. (A.4) between −z_{d} and z_{d}, we obtain immediately that any perturbation of B_{y} is symmetric with respect to the midplane (if the background equilibrium B_{y} is also symmetric). We then integrated Eq. (A.4) between z_{i} and z, with z_{i} chosen to have . This altitude is always defined since B_{y} is a decreasing function of z and is located in general below the Alfven point since B_{y} is close to 0 at this altitude. By averaging in z and assuming ρ[z_{i}]≈0, we obtain
If we assume that the density profile is not too far from a Gaussian distribution of width H_{m} = hH with h > 1, then we have:
and we find
We linearised this relation and obtain the following relation for the dimensionless toroidal perturbation
is the “toroidal magnetisation”. The total linearised stress can be then written
Appendix B: Extension of the linear theory to laminar discs or discs with strong toroidal field
In this appendix, we have relaxed assumptions (1) and (3) in Section 3.6 and extended the linear analysis to a more general configuration. The geostrophic balance is replaced by a magnetogeostrophic balance, into which we included the effect of toroidal magnetic pressure. The radial stress can be purely laminar or a mix of both turbulent and laminar components. The following calculation applies in particular to discs with ambipolar diffusion. The linearised system of equations is:
with λ ≈ 0.5, a parameter introduced in Appendix A and the toroidal magnetisation of the equilibrium. The growth rates follow the dispersion relation
with
We note that the case with a pure laminar stress is equivalent to that of a pure turbulent stress, but with q = 0. Thus, the instability criterion for laminar disc is less restrictive. In practise, given the fact that q = 0 and that turbulent diffusion is absent, growth rates can be enhanced (by a factor ≈2) compared to the turbulent case. A strong toroidal field can also help the instability, by affecting the geostrophic equilibrium and enhancing the zonal flows.
Appendix C: Stability of laminar disc with initial 1D uniform wind
To rigorously demonstrate the existence of an instability, similar to that predicted in Sect. 3, we use a very simple setup which consists of an initial uniform and laminar wind solution (independent on x). In that case, the radial stress is provided by the mean −B_{x}B_{y} relative to the laminar solution.
The first step is to compute a wind equilibrium in the local frame. For that, we employ the same technique as Lesur et al. (2013): we run a 1D shearing box simulation with an initial hydrostatic disc threaded by a vertical field of strength B_{z0} = 0.044 (μ_{eq} = 10^{−3}). We added some random noise at t = 0 and let the system evolves until it reaches a new equilibrium. Inbetween, an MRI mode is triggered and lead to the launch of a magnetocentrifugal wind. Figure C.1 (top) shows the equilibrium obtained in a box that spans −15 to 15H in the vertical direction. This class of solution is characterised by a very strong toroidal field and develops only if the system is initialised with a very clean 1D perturbation, so it is likely that it is never encountered in nature. The wind is however more realistic with mass loss efficiency ζ_{0} ≃ 0.01.
Fig. C.1.
Top: some quantities related to the 1D wind equilibrium. Bottom: amplitude (in logarithm) of the unstable Fourier mode k_{x} = π/10 as a function of time. The background equilibrium is a laminar 1D wind computed for μ_{eq} = 10^{−3}. Blue and green lines are, respectively, the surface density and vertical magnetic field perturbation. The blue dashed lines delimits the linear growth phase. 
We then ran a 2D axisymmetric simulation, starting from the wind solution and add a small perturbation of the form exp(ik_{x}x) with k_{x} = 2π/L_{x}. Here we use a box size L_{x} = 20H and resolution (N_{X}, N_{Z})=(256, 512). Figure C.1 (bottom) shows the timeevolution of the axisymetric perturbation and , projected onto the Fourier component k_{x} = 2π/L_{x}. Clearly, the mode is amplified exponentially during the first 200 Ω^{−1} with growth rate σ ≃ 0.0074. We show in Fig. C.2 the shape of the unstable mode in the poloidal (x,z) plane. The density field develops a ring structure while the magnetic perturbation, anticorrelated with ρ, forms an inclined shell similar to that obtained in 3D turbulent simulations (see Fig. 2 for comparison). The radial velocity of the mode is anticorrelated with and directed towards the gap. The outflow and mass loss rate perturbations ρ_{0}û_{z} are positive inside the gap and negative outside. All of these properties are indications that the mode is triggered by the same instability as described in Sect. 3. One other interesting result is the apparent correlation between the density and the B_{y} structures, near the midplane. Such behaviour results from the magnetohydrostatic equilibrium in the vertical direction (very well checked) and actually confirms our calculation in Appendix A. In particular we find when quantities are integrated vertically up to the Alfven point (instead of λ ≃ 0.5 in Appendix A). This parameter does not seem to depend on the integration boundary neither on k_{x}.
We finally achieved the same simulation for different k_{x}. The growth rates obtained for seven different radial wavenumbers are shown in Fig. C.3 (red diamond markers). The blue dashed curve is the model prediction, assuming α_{L} = 0.4, ζ_{0} = 0.01, λ ≃ 0.2, λμ_{0y} = 3.2 and p = 0.5. We note that we use the extended model of Appendix B to compute the growth rates. The value of p is different from turbulent simulations and has been obtained by varying the magnetisation and measuring ζ of different 1D wind solutions. Although the maximum growth rate is in agreement with the numerical data, the behaviour at large k_{x} is not well reproduced by our model. The reason is the omission of the EMF term u_{z}b_{x} in the induction equation, which clearly enhances the instability at large k_{x}. We measured the ratio u_{z}b_{x}/u_{x}b_{z} for each k_{x} and include in the model a correction factor c_{o}(k_{x}) in the linearised induction equation. With such a correction, we obtain a very good match between the numerical and theoretical growth rates at large k_{x} (plain green curve in Fig. C.3).
Fig. C.2.
Shape, in the poloidal plane, of the unstable axisymmetric mode k_{x} = π/10 around the 1D wind equilibrium (μ_{eq} = 10^{−3}). 
Fig. C.3.
Growth rate of axismmyetric modes measured in the laminarwind simulations for different k_{x} (red diamond makers). The bluedashed curve shows the theoretical growth rate predicted from the model extension in Appendix B (assuming a pure laminar stress). The green plain curve shows the same model but taking into account the u_{z}B_{x} term in the induction equation (this term is simply estimated by calculating the ratio u_{z}B_{x}/u_{x}B_{z} associated with the unstable mode in simulations. 
Appendix D: Check of model assumptions in simulations
To check the assumptions made in Sect. 3.6 and validate the linear theory of Sect. 3, we examine in details the mass, momentum, and magnetic budget in the 2D ambipolar run with μ = 10^{−3}. Figure D.1 shows the different flux and source terms in Eqs. (3)–(6), averaged during the growth phase (delimited by the dashed lines in Fig. D.2). The terms are calculated in real space (and thus comprise all the axisymmetric modes) and are integrated within ±1H. The first panel shows that the gaps are depleted by the vertical mass flux while they are refilled by the radial mass flux. It can be noted that the first term is slightly larger than the second. This ensures that the instability is driven by the outflow. The second panel clearly demonstrates that perturbations are in a geostrophic balance (Assumption 1 in Sect. 3.6), the Coriolis force being equilibrated by the radial thermal pressure gradient. In the azimuthal direction (third panel), both the radial flux of turbulent stress and the inertial term contribute to the production of v_{x}, directed towards the gaps. As suggested by the linear equations in Sect. 3.6, the inertial term due to the zonal flow can be an additional source of radial transport and even substitute for the α viscosity.
In the final panel we show that the vertical field in the gaps is mainly produced via the EMF term associated with radial motion –v_{x}B_{z}, as expected from the theory. The role of the vertical EMF component v_{z}B_{x} is more ambiguous: it seems to reinforce B_{z} in the nascent gaps at x ≃ −8 H and x = 0 but has a negative effect in the central gap at x = −5. Actually, Fig. D.2 shows that this structure has already evolved into a nonlinear regime, unlike the others. In any case, this EMF term can be comparable to −v_{x}B_{z} and thus, our assumption (4) in Sect. 3.6 is not necessarily satisfied. The ambipolar term η_{A}J_{⊥y} (with J_{⊥y} the y projection of the current perpendicular to B), is always negative within the gaps and positive outside. Therefore it tends to diffuse the magnetic field out of the gaps. This last result is clear evidence that the process forming rings (or gaps) is supported by the ideal terms and not by the ambipolar diffusion. In conclusion we showed that the force balance assumed in the linear regime (Sect. 3.6) is compatible with numerical data.
Fig. D.1.
Radial profiles of the flux and source terms in Eqs. (3)–(6), integrated vertically within z ± H and averaged in time during the linearphase for μ_{eq} = 10^{−3} and Am = 1. The linear phase is delimited by the vertical dashed lines in Fig. 9 (60 < t < 200 Ω^{−1}). 
Fig. D.2.
Spacetime diagrams showing the vertically averaged B_{z} in the 2D ambipolar simulation (Am_{mid} = 1, μ_{eq} = 10^{−3} and z_{d} = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially. 
All Figures
Fig. 1.
Ring occurrence in MHD simulations with net vertical magnetic flux. The orange check mark means that zonal flow exist, but their persistence within the turbulent flow and their convergence with box size are uncertain at current time. (1) Suriano et al. (2017), (2) Bai & Stone (2014), (3) our own simulations, (4) Hawley (2001), (5) Kunz & Lesur (2013), (6) Béthune et al. (2016), (7) Krapp et al. (2018), (8) Bai (2015), (9) Riols & Lesur (2018), (10) Suriano et al. (2018, 2019). 

In the text 
Fig. 2.
Example of a ring and gap structure and its surrounding outflow topology in the 3D shearing box simulation of Riols & Lesur (2018) with ambipolar diffusion. Top panel: gas density (colourmap in log scale) and streamlines (cyan lines) in the poloidal plane (x and z are respectively the radial and vertical coordinates). Bottom panel: vertical magnetic field (colourmap) and poloidal magnetic field lines (in white) showing the inclined “plume” structure from where most of the mass is extracted. The thickness of the lines is proportional to the field intensity. 

In the text 
Fig. 3.
Radial velocity computed from the 3D simulation of Fig. 2. It is averaged in the azimuthal direction, in z within ±1.5H and during the first 50 orbits (∼300 Ω^{−1}), corresponding to the growth of the ring structure. 

In the text 
Fig. 4.
Sketch illustrating the development of rings in accretion discs. The black and blue arrows represent the radial flow and the wind, respectively. The green patches correspond to density minima where the vertical field B_{z} and the magnetisation μ grow. In this representation, the poloidal field lines (// to the wind) do not bend in the midplane but cross the disc with some angle. This configuration is actually not intuitive but is encountered in various disc simulations around the midplane. 

In the text 
Fig. 5.
Growth rate of the winddriven instability as a function of the radial wavenumber for different magnetic diffusivities. Parameters of the model have been chosen in agreement with numerical simulations for fiducial μ_{eq} = 10^{−3}: α_{ν} ≃ 0.3, p ≃ 1, q ≃ 0.5, ζ_{0} ≃ 0.01. 

In the text 
Fig. 6.
Evolution of the transport coefficient α_{ν}, α_{L} and the wind mass loss efficiency ζ in the 2D simulation without explicit diffusion (μ_{eq} = 10^{−3}) 

In the text 
Fig. 7.
Top: spacetime diagram showing the column density as a function of time and x, in 2D axisymmetric turbulent simulations, for different magnetisation (left to right, μ_{eq} = 10^{−4}, 10^{−3} and 10^{−2}). The density is integrated within z ± 1.5H and runs were computed in the ideal limit (without explicit magnetic diffusion). Bottom panels: for each run, timeevolution of , normalised projection of B_{z} onto the prominent axisymmetric Fourier mode (with largest amplitude). This mode corresponds respectively to k_{x} = 6, 4 and 2k_{x0} for μ_{eq} = 10^{−4}, 10^{−3} and 10^{−2}, with k_{x0} = 2π/L_{x} the fundamental box radial wavenumber. 

In the text 
Fig. 8.
Optimal growth rate of the instability as a function of the magnetisation μ. The red diamond markers are the growth rates of the most unstable Fourier mode measured in 2D ideal simulations. The orange and green crosses are those measured in resistive and ambipolar runs. The purple solid line correspond to the model prediction, using the scaling relations (17) and (18). The blue and yellow dashed lines use the same model but with constant η^{⋆} = 0.01 in both cases and constant α_{ν} = 0.1 for the last case. 

In the text 
Fig. 9.
Top panel: spacetime diagrams showing the surface density Σ(x, t) in the 2D ambipolar simulation (Am_{mid} = 1, μ_{eq} = 10^{−3} and z_{d} = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially. Bottom panel: evolution of the B_{z} component projected onto the k_{x} = 4k_{x0} mode. 

In the text 
Fig. 10.
Coefficients p and q measured from the 2D ideal simulation with μ_{eq} = 10^{−4}. The red (respectively blue) line shows the ratio between the perturbation of ζ (respectively α) and the perturbation associated with magnetisation, as a function of time. Each quantity (ζ, α and μ) is projected onto the mode k_{x} = 6k_{x0}, averaged within z ± 2H and normalised to its mean. 

In the text 
Fig. 11.
Contrast between rigs and gaps versus radial separation. The plain curves are estimated from the theory (Eq. (19)) while the diamond markers are points measured from simulations. The blue, orange, green, and red colours correspond respectively to μ_{eq} = 10^{−4}, 10^{−3.5}, 10^{−3}, and 10^{−2.5}. Dotted lines delimit the region above which dust can be concentrated within the rings, assuming H/R = 0.05 (black) and H/R = 0.1 (purple). See Sect. 4.7 for further detail. 

In the text 
Fig. 12.
Top: spacetime diagram showing the column density as a function of time and x, in 3D turbulent simulations with μ_{eq} = 10^{−3}. Left panel: no explicit diffusion; centre panel: ohmic diffusion (Rm = ΩH^{2}/η = 100); right panel: ambipolar diffusion (Am_{mid} = 1). The density is integrated within z ± 1.5H. Bottom panels: for each run, timeevolution of turbulent α_{ν}, laminar α_{L} and ζ the mass loss efficiency. 

In the text 
Fig. 13.
Left: timeevolution of , the normalised projection of B_{z} onto the prominent axisymmetric Fourier mode (k_{x} = 2π/L_{x}) in the ideal simulation with μ_{eq} = 10^{−3}. Right: growth rate of the linear instability as a function of the magnetisation μ. The red diamond markers are measures from 3D ideal simulations. The green and yellow crosses are those measured in resistive and ambipolar runs. The plain lines correspond to model predictions, using similar scaling relations as (17), but with η_{t} = 0.5α_{ν} in the ideal case. 

In the text 
Fig. C.1.
Top: some quantities related to the 1D wind equilibrium. Bottom: amplitude (in logarithm) of the unstable Fourier mode k_{x} = π/10 as a function of time. The background equilibrium is a laminar 1D wind computed for μ_{eq} = 10^{−3}. Blue and green lines are, respectively, the surface density and vertical magnetic field perturbation. The blue dashed lines delimits the linear growth phase. 

In the text 
Fig. C.2.
Shape, in the poloidal plane, of the unstable axisymmetric mode k_{x} = π/10 around the 1D wind equilibrium (μ_{eq} = 10^{−3}). 

In the text 
Fig. C.3.
Growth rate of axismmyetric modes measured in the laminarwind simulations for different k_{x} (red diamond makers). The bluedashed curve shows the theoretical growth rate predicted from the model extension in Appendix B (assuming a pure laminar stress). The green plain curve shows the same model but taking into account the u_{z}B_{x} term in the induction equation (this term is simply estimated by calculating the ratio u_{z}B_{x}/u_{x}B_{z} associated with the unstable mode in simulations. 

In the text 
Fig. D.1.
Radial profiles of the flux and source terms in Eqs. (3)–(6), integrated vertically within z ± H and averaged in time during the linearphase for μ_{eq} = 10^{−3} and Am = 1. The linear phase is delimited by the vertical dashed lines in Fig. 9 (60 < t < 200 Ω^{−1}). 

In the text 
Fig. D.2.
Spacetime diagrams showing the vertically averaged B_{z} in the 2D ambipolar simulation (Am_{mid} = 1, μ_{eq} = 10^{−3} and z_{d} = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially. 

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