Open Access
Issue
A&A
Volume 625, May 2019
Article Number A108
Number of page(s) 18
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201834813
Published online 22 May 2019

© A. Riols and G. Lesur 2019

Licence Creative Commons
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 radio-interferometer ALMA and the new generation of instruments like SPHERE at the Very Large Telescope have imaged a variety of structures in proto-planetary 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 dust-drift-driven 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 magneto-hydrodynamics (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 magneto-rotational 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 large-scale 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, self-organise into large-scale 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 non-ideal 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 “anti-diffusion” 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 non-ideal 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 “viscous-type” 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 wind-driven 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 “mass-flux” (or “stripe” instabilities) seen by Moll (2012) and Lesur et al. (2013) in the highly magnetised (MRI-stable) 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 wind-driven 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 well-defined growth rates corresponding to those predicted by the theory. We have also explored the non-linear 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 self-organisation

2.1. Zonal flows and rings occurrence in MHD simulations

Self-organisation 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 sub-Keplerian and super-Keplerian 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 sub-structure. Figure 1 summarises ten years of research and numerical effort to identify the conditions under which zonal flows can develop.

thumbnail 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).

Open with DEXTER

In 3D unstratified simulations, zonal structures seem absent in the ideal or ambipolar cases but arise in the Hall-dominated regime via an anti-diffusive 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 “viscous-type” 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 pile-up 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 self-organised 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 non-ideal 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 large-scale 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 non-axisymmetric 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 large-scale 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 beta-parameter in the midplane)

with , is stronger in the gaps and therefore anti-correlated 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).

thumbnail 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.

Open with DEXTER

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 vx within the midplane (|z|≲H) is opposed to the formation of the ring. In particular, it is anti-correlated 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 wind-driven instability is probably connected with the formation of these plumes, which organise the flow into radial sub-structures.

thumbnail 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.

Open with DEXTER

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 = αcsH. 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,

(1)

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 large-scale 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 Bz in the density minima or gaps (green patches). Because matter and magnetic flux are transported radially at the same speed δvR, 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.

thumbnail 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 Bz 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.

Open with DEXTER

Quantitatively, an instability occurs whenever the surplus of mass launched in the wind is larger than the mass flowing radially. If we note ra the typical separation between the rings and , the instability criterion becomes

Using relation (1) and the transport equation for Bz, 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ϕ >  vK in the regions where vR <  0 and vϕ <  vK where vR >  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 re-distribution of angular momentum. Finally, we do not yet know how the instability behaves if the disc is subject to magnetic diffusion and non-ideal 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 & Lynden-Bell 1965). This corresponds to a Cartesian patch of the disc, centred at r = R0, where the Keplerian rotation is approximated locally by a linear shear flow plus a uniform rotation rate Ω = Ω ez. 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 −zd and zd, where zd is some arbitrary altitude, and a mass-weighted 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

(2)

We also note the difference between the field at z = zd and z = −zd. The compressible, inviscid and isothermal equations of motion in the horisontal plane, integrated azimuthally and vertically between −zd and zd are

(3)

(4)

(5)

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 ⟨Bz⟩:

(6)

where ℰy = vzBx − vxBz 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 Wij. In the limit of highly sub-sonic 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 Wxz ≃ 0. In the azimuthal momentum equation, however, the stress is comparable to others terms. By using the decomposition of Eq. (2), we have

(7)

where

can be identified respectively as a turbulent and laminar radial transport coefficient. The term related to the vertical stress Wyz can be written as

(8)

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 = ανcsΩ). 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 sub-script “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 (z-dependent) wind solutions studied by Lesur et al. (2013) and Riols et al. (2016). The solenoidal condition gives immediately δBz = 0, which means that . We can then define a constant magnetisation of the equilibrium, as

Solutions can be either symmetric about the midplane with Bx = By = 0 at z = 0 or anti-symmetric with ∂zBx = ∂zBy = 0.

  1. In the symmetric case, we have ⟨Bx⟩ = ⟨By⟩ = 0 but . Horizontal components satisfy but . Therefore there is a net vertical stress through the disc, which is responsible for a mean accreting flow.

  2. In the second case ⟨Bx⟩ = ⟨By⟩≠0 but . We define Bx0 and By0 the mean radial and toroidal magnetic field throughout the midplane. Departures to the vertical average satisfy but . It is straightforward to check that Wyz0 = 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 anti-symmetric solutions have been shown to naturally emerge from turbulent shearing box simulations (Lesur et al. 2014; Bai 2015), and even from global simulations when non-ideal effects are included (Béthune et al. 2017; Bai 2017). For that reason, we focus mainly on the anti-symmetric solutions in the following sections.

3.6. Linearisation around equilibrium

We use the sub-script “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 kxH ≲ 1 and assume σ ≪ Ω. To simplify the problem, we made the following further assumptions:

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

  2. The vertical component of the Reynolds stress tensor as well as the turbulent vertical stress αW = 0 are neglected

  3. The stress is purely turbulent αν0 ≫ αŁ0

  4. In Eq. (6), the vertical stretching of radial field ⟨uz⟩⟨Bx⟩ and the vertical diffusion of Bz are neglected.

  5. We have assumed that . This is true if the radial velocity perturbation does not vary too much between -zd and zd.

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 kx to H and σ to Ω and define η = (η + η0)/(H2Ω). Linearisation of Eqs. (3)–(6) around equilibrium leads to

(9)

(10)

(11)

(12)

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

(13)

Solutions of the problem are found by setting D = 0. It is straightforward to show that growth rates follow the dispersion relation

(14)

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 non-ideal 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 kx ≫ 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 MRI-driven 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 kx for three different resistivities (η ≃ 0.002, 0.01, 0.06). These values are here purely ad-hoc and just to illustrate the dependence of σ(kx) 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 kxH ≲ 0.5H which correspond to rings of size ra ≳ 10 H. The instability occurs on long timescale 1/σ larger than 150 Ω−1 (20 orbits).

thumbnail Fig. 5.

Growth rate of the wind-driven 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.

Open with DEXTER

3.8. Eigenmodes

The linearised system admits solutions of the form:

(15)

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 anti-correlated with the zonal flow . The vertical field and the magnetisation are anti-correlated 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

Shearing-box simulations were run with the PLUTO code (Mignone et al. 2007), a finite-volume 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 shear-periodic 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 Bx = 0 and By = 0 at z = ±Lz/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),

(16)

where is the mass injection rate computed at each time step to maintain the total mass constant in the box, and zi 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 Lx = Ly = 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 small-scale 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 sub-structures, 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.

thumbnail Fig. 6.

Evolution of the transport coefficient αν, αL and the wind mass loss efficiency ζ in the 2D simulation without explicit diffusion (μeq = 10−3)

Open with DEXTER

thumbnail Fig. 7.

Top: space-time 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, time-evolution of , normalised projection of Bz onto the prominent axisymmetric Fourier mode (with largest amplitude). This mode corresponds respectively to kx = 6, 4 and 2kx0 for μeq = 10−4,  10−3 and 10−2, with kx0 = 2π/Lx the fundamental box radial wavenumber.

Open with DEXTER

To understand whether rings form via a linear instability or a more complicated non-linear process, we explored the dynamics of axisymmetric modes in Fourier space. The procedure is simple: at each time-step, 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 simulation1. 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 kx = 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 quasi-exponentially, 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 non-linear 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.

thumbnail 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.

Open with DEXTER

4.3. Non-ideal 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 η/(ΩH2) = 0.01. This corresponds to a magnetic Reynolds number Rm = ΩH2/η = 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 kx = 3kx0 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 vertically-averaged Bz. Again three or four structures seem to emerge from the initial turbulent phase. The growth rate associated with the kx = 4kx0 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.

thumbnail Fig. 9.

Top panel: space-time diagrams showing the surface density Σ(x, t) in the 2D ambipolar simulation (Ammid = 1, μeq = 10−3 and zd = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially. Bottom panel: evolution of the Bz component projected onto the kx = 4kx0 mode.

Open with DEXTER

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 ηAJ × (B × B)/B2 merely acts as a diffusion on the Bz 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 quasi-laminar discs with αL0 ≳ αν0 (see Appendix B), typically those obtained in non-ideal 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 quasi-steady 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 large-scale 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 kx = 6 × 2π/Lx and normalised with respect to the mean (kx = 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:

(17)

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).

thumbnail 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 kx = 6kx0, averaged within z ± 2H and normalised to its mean.

Open with DEXTER

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 ∂⟨Bz⟩/∂x and indeed acts to diffuse the large scale structure. The ratio between the two terms is ηt/(ΩH2) ≈ 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:

(18)

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 large-scale 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. Non-linear saturation and ring/gap contrast

As suggested by Fig. 7, the ring instability saturates and enters a non-linear 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 non-linear saturation leads to a drop in the mean stress (Fig. 6) and the production of a strong mean azimuthal field By, 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 non-linear 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 Bz ∼ Bz0). 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 non-linear regime can be then estimated by setting in Eq. (15). We obtain in the ideal limit:

(19)

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 kx. 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 kx = 2, 4, 5, 5kx0 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.

thumbnail 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.

Open with DEXTER

4.6. Ring separation

So far, we have simply measured the ring separation as the most prominent kx 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 kx = 5kx0 = 1.5 H−1 and 6kx0 = 1.8 H−1, while the linear theory gives a maximum growth rate at kx ≃ 0.45 H−1. Actually, this is not in contradiction with the theory since σ appears to be quite flat with kx in this regime. The reason is that turbulent dissipation is low (Pmt = 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 Bz associated with the initial turbulence peaks around kx = 6kx0 and beyond for μeq = 10−4, while it peaks at much lower kx ≲ 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 kx = 6kx0 mode and the large-scale mode which, in principles, should grow fastest (kx ≃ 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:

(20)

where P0 is the mean equilibrium pressure profile and kR the local radial wave number of the ring perturbation. A typical disc model, which matches disc observations (Andrews et al. 2009), has

(21)

The criterion for dust accumulation becomes:

(22)

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 non-ideal 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 Ammid = 1. These simulations are initialised with random perturbations and are computed with Lx = Ly = 20 H, Lz = 8 H and resolution NX = NY = 256, NZ = 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 non-axisymmetric) 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 kx (unlike the 2D case for which it was quite flat.) The maximum growth rate σ ≃ 0.0039 is obtained for kxH = 0.2 and stability is reached for kxH ≳ 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.

thumbnail Fig. 12.

Top: space-time 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 = ΩH2/η = 100); right panel: ambipolar diffusion (Ammid = 1). The density is integrated within z ± 1.5H. Bottom panels: for each run, time-evolution of turbulent αν, laminar αL and ζ the mass loss efficiency.

Open with DEXTER

thumbnail Fig. 13.

Left: time-evolution of , the normalised projection of Bz onto the prominent axisymmetric Fourier mode (kx = 2π/Lx) 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.

Open with DEXTER

In non-ideal 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 −⟨Bx⟩⟨By⟩. Azimuthal structures are very faint and have little impact on the unstable dynamics. We did indeed check that turbulent diffusion of Bz 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 Lx = 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 large-scale 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 “no-stress” 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 ζ = (ρuz)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 non-ideal 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 T-Tauri 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 large-scale 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 low-velocity 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 (photo-evaporation 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 long-lived rings of size ≫H, which might remain stable during a large fraction of the disc lifetime. Such wind-driven 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 wind-driven 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 “radial-drift” 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 back-react 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 self-gravity 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 X-ray 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.


1

This is defined as the mode with maximum amplitude, which corresponds to the kx = 4kx0 = 8π/Lx 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 kx = 5 and kx = 6kx0 components.

Acknowledgments

This work acknowledges funding from The French ANR under contracts ANR-17-ERC2-0007 (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.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07-13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

References

Appendix A: Linearisation of total stress in the azimuthal momentum equation

In this appendix, we describe how we linearised the total stress in the y-momentum 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

(A.1)

(A.2)

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 Wyz and the x derivative of the averaged radial stress gives:

(A.3)

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 magneto-hydrostatic equilibrium and the magnetic field is essentially toroidal:

(A.4)

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 −zd and zd, we obtain immediately that any perturbation of By is symmetric with respect to the midplane (if the background equilibrium By is also symmetric). We then integrated Eq. (A.4) between zi and z, with zi chosen to have . This altitude is always defined since By is a decreasing function of z and is located in general below the Alfven point since By is close to 0 at this altitude. By averaging in z and assuming ρ[zi]≈0, we obtain

(A.5)

If we assume that the density profile is not too far from a Gaussian distribution of width Hm = hH with h >  1, then we have:

(A.6)

and we find

(A.7)

We linearised this relation and obtain the following relation for the dimensionless toroidal perturbation

(A.8)

is the “toroidal magnetisation”. The total linearised stress can be then written

(A.9)

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 magneto-geostrophic 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:

(B.1)

with λ ≈ 0.5, a parameter introduced in Appendix A and the toroidal magnetisation of the equilibrium. The growth rates follow the dispersion relation

(B.2)

with

(B.3)

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 −BxBy 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 Bz0 = 0.044 (μeq = 10−3). We added some random noise at t = 0 and let the system evolves until it reaches a new equilibrium. In-between, an MRI mode is triggered and lead to the launch of a magneto-centrifugal 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.

thumbnail Fig. C.1.

Top: some quantities related to the 1D wind equilibrium. Bottom: amplitude (in logarithm) of the unstable Fourier mode kx = π/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.

Open with DEXTER

We then ran a 2D axisymmetric simulation, starting from the wind solution and add a small perturbation of the form exp(ikxx) with kx = 2π/Lx. Here we use a box size Lx = 20H and resolution (NX, NZ)=(256, 512). Figure C.1 (bottom) shows the time-evolution of the axisymetric perturbation and , projected onto the Fourier component kx = 2π/Lx. 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, anti-correlated 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 anti-correlated with and directed towards the gap. The outflow and mass loss rate perturbations ρ0z 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 By structures, near the midplane. Such behaviour results from the magneto-hydrostatic 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 kx.

We finally achieved the same simulation for different kx. 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 kx is not well reproduced by our model. The reason is the omission of the EMF term uzbx in the induction equation, which clearly enhances the instability at large kx. We measured the ratio uzbx/uxbz for each kx and include in the model a correction factor co(kx) in the linearised induction equation. With such a correction, we obtain a very good match between the numerical and theoretical growth rates at large kx (plain green curve in Fig. C.3).

thumbnail Fig. C.2.

Shape, in the poloidal plane, of the unstable axisymmetric mode kx = π/10 around the 1D wind equilibrium (μeq = 10−3).

Open with DEXTER

thumbnail Fig. C.3.

Growth rate of axismmyetric modes measured in the laminar-wind simulations for different kx (red diamond makers). The blue-dashed 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 uzBx term in the induction equation (this term is simply estimated by calculating the ratio uzBx/uxBz associated with the unstable mode in simulations.

Open with DEXTER

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 re-filled 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 vx, 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 –vxBz, as expected from the theory. The role of the vertical EMF component vzBx is more ambiguous: it seems to reinforce Bz 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 non-linear regime, unlike the others. In any case, this EMF term can be comparable to −vxBz and thus, our assumption (4) in Sect. 3.6 is not necessarily satisfied. The ambipolar term ηAJy (with Jy 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.

thumbnail 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).

Open with DEXTER

thumbnail Fig. D.2.

Spacetime diagrams showing the vertically averaged Bz in the 2D ambipolar simulation (Ammid = 1, μeq = 10−3 and zd = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially.

Open with DEXTER

All Figures

thumbnail 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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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 Bz 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.

Open with DEXTER
In the text
thumbnail Fig. 5.

Growth rate of the wind-driven 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.

Open with DEXTER
In the text
thumbnail Fig. 6.

Evolution of the transport coefficient αν, αL and the wind mass loss efficiency ζ in the 2D simulation without explicit diffusion (μeq = 10−3)

Open with DEXTER
In the text
thumbnail Fig. 7.

Top: space-time 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, time-evolution of , normalised projection of Bz onto the prominent axisymmetric Fourier mode (with largest amplitude). This mode corresponds respectively to kx = 6, 4 and 2kx0 for μeq = 10−4,  10−3 and 10−2, with kx0 = 2π/Lx the fundamental box radial wavenumber.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 9.

Top panel: space-time diagrams showing the surface density Σ(x, t) in the 2D ambipolar simulation (Ammid = 1, μeq = 10−3 and zd = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially. Bottom panel: evolution of the Bz component projected onto the kx = 4kx0 mode.

Open with DEXTER
In the text
thumbnail 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 kx = 6kx0, averaged within z ± 2H and normalised to its mean.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 12.

Top: space-time 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 = ΩH2/η = 100); right panel: ambipolar diffusion (Ammid = 1). The density is integrated within z ± 1.5H. Bottom panels: for each run, time-evolution of turbulent αν, laminar αL and ζ the mass loss efficiency.

Open with DEXTER
In the text
thumbnail Fig. 13.

Left: time-evolution of , the normalised projection of Bz onto the prominent axisymmetric Fourier mode (kx = 2π/Lx) 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.

Open with DEXTER
In the text
thumbnail Fig. C.1.

Top: some quantities related to the 1D wind equilibrium. Bottom: amplitude (in logarithm) of the unstable Fourier mode kx = π/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.

Open with DEXTER
In the text
thumbnail Fig. C.2.

Shape, in the poloidal plane, of the unstable axisymmetric mode kx = π/10 around the 1D wind equilibrium (μeq = 10−3).

Open with DEXTER
In the text
thumbnail Fig. C.3.

Growth rate of axismmyetric modes measured in the laminar-wind simulations for different kx (red diamond makers). The blue-dashed 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 uzBx term in the induction equation (this term is simply estimated by calculating the ratio uzBx/uxBz associated with the unstable mode in simulations.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail Fig. D.2.

Spacetime diagrams showing the vertically averaged Bz in the 2D ambipolar simulation (Ammid = 1, μeq = 10−3 and zd = 2H). The two vertical dashed lines delimit the “linear” phase, during which zonal modes, in Fourier space, grow exponentially.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.