© A. Riols et al. 2020

1. Introduction

Protoplanetary discs, the nurseries of planets, are known to dissipate over a few million years. Measurements of UV excess emission indicate that these discs are accreting with rates ranging from 10⁻¹⁰ to 10⁻⁷ solar mass per year onto the central object (Venuti et al. 2014). The accretion process plays a central role in the formation of stars, but also in discs and planets evolution. However, its origin is still actively debated. For several decades, accretion was believed to result from turbulent motions, which act like an effective viscosity in the discs (Shakura & Sunyaev 1973) transporting the angular momentum outwards. The most commonly invoked excitation mechanism of turbulence in accretion discs is the magneto-rotational instability (MRI, Balbus & Hawley 1991; Hawley et al. 1995). Several MHD simulations in the ideal limit showed that the MRI could provide the desired accretion rate and transport efficiency (Flock et al. 2011, 2013; Bai & Stone 2013). However, due to their optical thickness and low temperatures, protoplanetary discs are weakly ionised in their midplane beyond 0.1–1 au (Gammie 1996). As a result, non-ideal MHD effects, such as ohmic diffusion, Hall effect, and ambipolar diffusion dominate at these radii and tend to suppress any form of MRI-driven turbulence (Fleming et al. 2000; Sano & Stone 2002; Wardle & Salmeron 2012; Bai 2013, 2015; Lesur et al. 2014).

Recent observations at sub-millimetric wavelength with ALMA seem to confirm that the level of turbulence in protoplanetary discs is indeed very low. Firstly, there is some evidence that the settling of millimetre dust is quite strong in some discs, like HLTau (Pinte et al. 2016), which excludes the presence of vigorous turbulent motions in these discs. Secondly, the measurements of gas velocity dispersion in CO (e.g. Flaherty et al. 2015, 2017) indicate that the departure from the Keplerian velocity in a number of discs is much weaker than what was theoretically predicted from ideal MRI simulations. Yet, these discs are known to accrete, which suggests that the transport of angular momentum is essentially driven by laminar processes.

Quite apart from the accretion problem, observations have revealed the presence of many sub-structures in several systems, both in the dust continuum (e.g. with ALMA), and in the scattered light emission (e.g with SPHERE). In particular, one of the most remarkable features is the concentric rings (or gaps) that appear in 60–70% of the observed discs (Garufi et al. 2018). Some examples of discs showing rings are 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 are possibly altering the long-term disc evolution and could be privileged locations of dust accumulation (Pinilla et al. 2012), a key step towards planetary core formation. How these rings form remains an open question, though a number of physical mechanisms have been proposed so far. The most commonly invoked scenario is the presence of planets opening gaps in the disc (Dong et al. 2015; Dipierro et al. 2015). Other scenarios involve snow lines (Okuzumi et al. 2016), dust-drift-driven viscous ring instability (Wünsch et al. 2005; Dullemond & Penzlin 2018), MHD turbulence with dead zones (Flock et al. 2015), zonal flows (Johansen et al. 2009), or secular gravitational instabilities in the dust (Takahashi & Inutsuka 2014).

In summary, the fundamental question is: are these rings the result of planet formation (as the planet-induced scenario suggests) or the driver (as envisioned by purely gas-induced mechanisms)? We note that answering this question is more difficult than it seems, since in both cases, one eventually expects to see planets in the core of the gaps. Hence, finding planets in these gaps cannot be the key signature allowing one to distinguish between these scenarios. Instead, one should rely on specific signatures of planet-free rings and gaps, which can be tested on the sky.

It has been recently suggested that magnetic winds in discs could play a crucial role in both the accretion and the ring formation processes. Molecular outflows consistent with magnetic winds have been observed in several Class 0 and I discs (Launhardt et al. 2009; Lumbreras & Zapata 2014; Bjerkeli et al. 2016; Tabone et al. 2017; Hirota et al. 2017; Louvet et al. 2018), although their origin is still debated. Simulations including a large-scale poloidal field and non-ideal effects (Hall effect, ambipolar diffusion, mainly) have shown that discs are able to launch a wind from their surface removing angular momentum and driving accretion in the process (Bai 2013; Lesur et al. 2014; Simon et al. 2015; Gressel & Pessah 2015; Béthune et al. 2017). The accretion is made possible by the laminar torque associated with the large-scale magnetic field supporting the wind. The typical accretion rates measured in these simulations are compatible with the observations. Moreover, large-scale axisymmetric structures or ‘zonal flows’ associated with rings of matter (Kunz & Lesur 2013; Bai 2015; Béthune et al. 2016, 2017; Suriano et al. 2018) form spontaneously in these non-ideal simulations. The ring formation has been attributed to several different MHD processes: self-organisation by ambipolar diffusion (Béthune et al. 2017), reconnection of poloidal fields lines (Suriano et al. 2018, 2019), or anti-diffusive processes (Bai 2014). More recently however, Riols & Lesur (2019) argued that their formation could result from a more generic process involving a wind instability and operating in both ideal and non-ideal MHD.

One challenging problem is to make a bridge between these non-ideal simulations and observations. Since current sub-millimetre or near-infrared observations by ALMA and SPHERE have mostly traced the dust emission so far (since the gas is more difficult to detect), it is crucial to understand the dynamics of dust grains in non-ideal MHD flows. In particular, there is concern surrounding the impact of gaseous rings and a large-scale wind on their spatial distribution. Several simulations by Zhu et al. (2015), Xu et al. (2017) have combined the dust dynamics with ambipolar diffusion but focused on a large particle regime (≳mm). More recently, Riols & Lesur (2018) studied the settling and dynamics of smaller grains in ambipolar-dominated flow, but their simulations were restricted to a local shearing box and do not capture the global structure of the wind flow. In this paper, we aim to revisit the simulations of Riols & Lesur (2018) but in the global configuration. Our ambition is also to go one step further and produce synthetic images of the disc by computing the dust emission. In this way, we aim to provide predictions on the typical dust settling, ring-gap contrast and spectral index that one might observe.

For that purpose, we performed global axisymmetric simulations of stratified discs with a net vertical magnetic field, including ambipolar diffusion and dust. For that, we used a customised version of the PLUTO code, based on PLUTO v4.3. To simplify the problem, we considered a portion of the disc spanning from 10 to 200 AU where we could neglect the Hall effect (Simon et al. 2015). The dust population is approximated as a pressure-less multi-fluid made up of different sized particles, from one hundred micrometres to one centimetre, relevant for radio-observations. The back reaction of the dust on the gas is taken into account, but we neglected coagulation or fragmentation processes. To produce synthetic images of the disc, we used the radiative transfer code MCFOST. As a preliminary step, the radiative transfer calculation was not performed during the hydrodynamic steps, but was based on a fixed snapshot (the hydrodynamic calculation was post-processed with MCFOST).

The paper is organised as follows: in Sect. 2, we describe the model and review the main characteristics of non-ideal physics and dust-gas interaction. In Sect. 3, we present the numerical methods used to simulate the dust-gas dynamics in a global model. In Sect. 4, we perform numerical simulations without dust in order to characterise the properties of the disc prone to ambipolar diffusion (wind, angular momentum transport and rings). In particular, we investigate the origin of the rings and their long-term evolution. In Sect. 5, we explore the dynamics of the dust in the ambipolar-dominated flows. We study, in particular, the sedimentation of the grains and show that pressure maxima associated with rings accumulate the solids into thin structures with high density contrast. In Sect. 6, we produce synthetic images of our disc simulation in ALMA bands by using MCFOST. Finally, in Sect. 7, we conclude and discuss the potential implications of our work for protoplanetary disc evolution and planet formation.

2. Theoretical framework

We model the gas and dust of a protoplanetary disc by using a global framework, in which quantities are assumed to be axisymmetric (2.5D framework). We note (r, θ, ϕ) the coordinates related to the spherical frame (θ being the co-latitude), and (R, ϕ,z) those related to the cylindrical frame. The disc extends from an inner cylindrical radius R_in = 10 AU to an outer limit R_out = 200 AU. We adopt a multi-fluid approximation in which the ionized gas and the dust interact and exchange momentum through drag forces.

2.1. Gas phase

2.1.1. Equations of motion

The gas is coupled to a magnetic field B and is assumed to be inviscid and locally isothermal on cylindrical radii, its pressure P and density ρ related by , where c_s(R, z) is a fixed sound-speed profile (see Sect. 2.1.4 for temperature prescription). The evolution of its density ρ, and total velocity field v follows:

(1)

(2)

where the total velocity field can be decomposed into a mean Keplerian flow plus a perturbation u:

(3)

The term g = −GM/r²e_r corresponds to the gravitational field exerted by the central star, where e_r is the unit vector pointing in the radial direction in the spherical frame. The last term in the momentum Eq. (2) contains the acceleration γ_{d − > g} exerted by the dust drag force on a gas parcel (detailed in Sect. 2.2.1). We note that the unit of time is taken as the inner Ω⁻¹ and the unit of length corresponds to R = 10 AU. The unit of mass is such that ρ = 1 in the midplane at R = 10 AU.

The magnetic field B, which appears in the Lorentz force (∇×B)×B, is governed by the non-ideal induction equation,

(4)

The term ∇ × [η_A(J×e_b)×e_b] corresponds to ambipolar diffusion, with J = ∇ × B the current density, e_b = B/‖B‖ the unit vector parallel to the field line, and η_A the ambipolar diffusivity. To simplify the problem, we assume that the diffusion of magnetic field is only due to ion-neutral collisions (ambipolar diffusion) and is assumed to be the dominant non-ideal effect in the regions considered in this paper, that is R ≥ 10 AU (see justifications in Simon et al. 2015). Therefore, the Hall effect and ohmic diffusion are neglected. The diffusivity η_A is not uniform, and its vertical profile is prescribed in Sect. 2.1.3.

2.1.2. Surface density profile

We assume a disc model with surface density following

(5)

This corresponds to a rather flat density profile that has been measured in some discs like HD 163296 (Williams & McPartland 2016). We note that the power-law index of the surface density inferred from observations ranges from −1.1 to 0.2 (Kwon et al. 2015). This choice is also motivated by theoretical considerations. Indeed, most (if not all) of the MHD models published to date assume Sigma ∝R⁻¹ with H/R = cst. This choice of density profile turns out to be peculiar in the usual alpha disc theory (Shakura & Sunyaev 1973), since it predicts no accretion, even if α is non-zero. This was indeed observed by Béthune et al. (2017) and Bai (2017). We therefore choose a shallower disc profile, so that even if no wind is present, turbulence alone could drive a measurable accretion rate (see Sect. 4.2).

2.1.3. Ambipolar profile

We consider here that the only charged particles are the ions and electrons. The ambipolar diffusivity η_A is given by:

(Wardle 2007), where γ_i = ⟨σv⟩_i/(m_n + m_i), and ⟨σv⟩_i = 1.3 × 10⁻⁹ cm³ s⁻¹ is the ion-neutral collision rate. ρ_n and ρ_i are, respectively, the neutral and ion mass density; m_n and m_i their individual mass. By introducing the Alfvén speed , we also define the dimensionless ambipolar Elsässer number

To calculate ρ_i, one generally needs to take into account the various ionisation sources (X rays, cosmic rays, radioactive decay and FUV from the star) as well as the dissociative recombination mechanisms. This requires the computation of the complex chemistry occurring in the gas phase and on grain surfaces. To simplify the problem as much as possible, we use a crude model that captures the essential physics of disc ionisation, which is the effect of FUV in the disc corona. In this model, Am is constant and of order unity in the midplane, and it increases abruptly up to a certain height z_io. The height z_io corresponds to the base of the ionisation layer, above which the FUV can penetrate. If we note Σ_i(z) as the horizontally averaged column density, integrated from the vertical boundary, we assume that FUV are blocked for Σ_i(z) > Σ_ic = 0.0005 g cm⁻². This corresponds to z_io ≃ 3.5H, where H is the gas disc scale height. To avoid any discontinuity at z = z_io, we use a smooth function to make the transition between the midplane regions and the ionised layer, so that the Ambipolar Elsasser number is

(6)

where Am₀ is the constant midplane value, Σ⁺ is the column density integrated from the top boundary, and Σ⁻ the column density integrated from the bottom boundary. This diffusivity profile is designed to mimic the conductivity obtained from complete thermochemical models of protoplanetary discs including dust grains, with Am₀ ≃ 1. (Thi et al. 2019). We note that to convert numerical values of Σ_i into g cm⁻², we used the disc model of Eq. (5).

2.1.4. Temperature and hot corona

Complex thermo-chemical models of protoplanetary discs (Kamp & Dullemond 2004; Thi et al. 2019) suggest that there is a sharp transition in temperature between the disc midplane and its ionised corona. The warm and optically thin corona can reach temperatures of a few thousand Kelvin, while the disc temperature remains below 100 K. The corresponding ratio of corona to midplane sound speed ranges typically from three to 10. To reproduce this temperature structure, we force the gas to be warmer above a certain height z_io ≃ 3.5H. This is chosen to be the same altitude as the one at which we impose a jump in ambipolar Elssaser number (Am). We take a constant corona-to-midplane temperature ratio at all radii equal to six. Again to avoid any discontinuity, we use a smooth function similar to that of Eq. (6) to make the transition between the two regions.

(7)

with

(8)

We note that inside the two distinct regions (disc and corona), we force the flow to maintain a constant temperature on cylindrical radii. For that, we evolve the energy equation with rapid cooling to keep the gas close to being locally isothermal (with a cooling time proportional to Ω⁻¹).

2.2. Dust phase

2.2.1. Equations of motion

The dust is composed of a mixture of different species, characterising different grain sizes. Each species, labelled by the subscript k, is described by a pressureless fluid, with a given density ρ_dk and velocity v_dk. We assume in this paper that dust grains remain electrically neutral so that they do not feel the magnetic field. The equations of motion for each species are:

(9)

(10)

with γ_{g − > dk} the drag acceleration induced by the gas on a dust particle. We obtain, by conservation of total angular momentum:

The acceleration associated with the drag and acting on a single particle is given by:

where τ_sk is the stopping time that corresponds to the timescale on which frictional drag causes a significant change in the momentum of the dust grain. This is a direct measure of the coupling between dust particles and gas. To evaluate τ_sk, we consider in this study that dust particles are spherical and small enough¹ to be in the Epstein regime (Weidenschilling 1977). For a spherical particle of size a_k and internal density ρ_s (which should not be confused with the fluid density ρ_d), the stopping time is:

(11)

A useful dimensionless quantity to parametrise the coupling between dust and gas is the Stokes number:

2.2.2. Converting particle size to Stokes number

In this paper, we keep the grain size a_k of each dust component fixed in space and time. To make an easier comparison with the literature, we also label each dust component by its Stokes number measured for ρ = ρ₀ = 1, c_s = c_s0 = 0.05 and Ω = Ω₀ = 1, which we label St₀ and refer to as the “reference Stokes number”. We note, however, that the local Stokes number (which controls the dynamics) depends on the position in the disc for one species. Therefore, we need to find a correspondence between the grain size and the Stokes number. For that, we assume that in the limit of small magnetisation, the disc is in hydrostatic equilibrium and its column (or surface) density is where ρ₀ is the midplane density and H = c_s/Ω is the disc scale height. Thus, combining these different relations, we obtain:

in the midplane. To estimate the Stokes number as a function of grain size, we assume that the disc surface density follows Eq. (5). By considering ρ_s = 2.5 g cm⁻³, we obtain the following conversion in the midplane:

(12)

2.3. Initial equilibrium

For the initial hydrostatic disc equilibrium, we use the same prescription as Zhu & Stone (2018):

(13)

(14)

where is the sound speed that is constant on the cylindrical radius:

(15)

The subscript “0” denotes quantities at the inner radius in the midplane. Taking Eq. (13) and developing at first order in z/R, we find that the disc is in a classical hydrostatic equilibrium with density varying on a characteristic vertical scale H = c_s/Ω. Considering the fact that Ω ∝ R^−3/2 and c_s given by Eq. (15), the ratio H/R is constant with radius and is given by c_s0/v_K0. Protoplanetary discs are geometrically thin, typically with H/R between 0.03 and 0.2 (Bitsch et al. 2013). In this paper, we thus choose H/R = c_s0/v_K0 = 0.05. A constant H/R is probably not the most encountered configuration in observed systems, although it is in line with some objects like HK Tau B or Oph163131.

We note that our initialisation is problematic at the pole, since the density, azimuthal velocity and sound speed become infinite at that location. To avoid this, we use R = max(R, R_min) in the above equations with R_min = 10 AU. The magnetic field is initialised in a way similar to that in Zhu & Stone (2018). Initially, there is a pure vertical field component, which has the following dependence on R:

with m = −2.5/2 and R_min = R₀ = 10 AU. Again the prescription at R <  R_min is to avoid the divergence of the vertical field at the pole. With such dependence, it is possible to define a constant magnetisation or plasma β parameter in the midplane for R >  R_min:

2.4. Averages

We define a vertical integral of a quantity Q between angles θ₋ and θ₊ as

(16)

and, respectively, a radial integral and temporal average:

(17)

2.5. Conservation of mass and angular momentum in spherical coordinates

Using the average defined in Eq. (16), the mass conservation equation then reads, in spherical coordinates:

(18)

where is the gas surface density defined on spherical geometry and is the mass loss rate measured between angles θ₋ and θ₊, which mark the location of the disc surface. If we note the vertically integrated radial stress between θ₋ and θ₊:

(19)

and the vertical stress difference between the two angles:

(20)

the conservation of angular momentum in spherical coordinates reads:

(21)

with , where we assumed that u ≪ ΩR, which allows us to neglect the time derivative term. We define, respectively, the accretion rate within the disc and the mass outflow rate as:

(22)

Multiplying Eq. (21) by r and integrating in radius, we obtain the following conservation equation:

(23)

Accretion onto the star can be either due to a radial transport of angular momentum, via the radial stress , or due to angular momentum extracted by a wind through the disc surface, via the vertical stress 𝒯_θϕ. An important quantity that characterises the radial transport efficiency is the parameter α (or Shakura & Sunyaev parameter) which is:

(24)

where is the Reynolds stress and is the Maxwell stress. The coefficient α can be decomposed into a sum of a laminar part:

(25)

plus a turbulent part: α_ν = α − α_L.

3. Numerical setup

3.1. Numerical code

To simulate the dynamics of the gas and dust, we use a modified version of the Godunov-based PLUTO code (Mignone et al. 2007). This code solves the approximate Riemann problem at each cell interface and is written so that it ensures conservation of mass and angular momentum. Inter-cell fluxes are computed with a HLLD solver for the gas and an HLL solver for the dust. We use a piecewise linear space reconstruction to estimate the Godunov fluxes at cell interfaces, with the monotonised central limiter. Physical quantities are evolved in time with a second order explicit Runge-Kutta scheme. Parabolic terms, including ambipolar diffusion, are also solved explicitly and introduce a strong constraint on the CFL condition. We note that we do not use the FARGO orbital advection scheme, as the CFL is limited by the parabolic terms due to ambipolar diffusion.

3.2. Grid and resolution

In radius, we use a logarithmic grid with 512 points, spanning from 10 AU to 200 AU. The increment δr is then proportional to the local radius r. In θ, the grid is split into three different regions. The first one is the midplane region from θ = 1.28 to θ = 1.86 (corrsponding to ±16° on either side of the midplane), which has 96 uniformly distributed cells. This corresponds to a vertical segment spanning ±5 scale height, with a resolution of nine points per H. The second and third regions are symmetric around the mid-plane and correspond, respectively, to θ = [0, 1.28] and θ = [1.86, π]. In each case, the grid is stretched and contain 64 points. Our grid resolution in r and θ is enough to capture most of the dynamics and the structures encountered in ambipolar-dominated discs.

3.3. Boundary and internal conditions

At the inner radius, we impose outflow boundary conditions for all the quantities except for the azimuthal velocity, on which we imposed the initial keplerian velocity given in Eq. (14). We also make sure that v_r <  0 to prevent any matter from entering the disc. At the outer radius, we apply similar outflow condition, but we do not force the keplerian equilibrium. In the θ direction, we implement a specific boundary condition allowing our simulations to extend all the way to the pole (from θ = 0 to θ = π). For that boundary, we use a similar implementation to Zhu & Stone (2018; see their Appendix).

To avoid small time steps, we limit the Alfvén velocity in the whole domain to eight times the Keplerian velocity at the inner disc radius : v_A/v_K0 <  8. This is done by artificially increasing the density in the cells where this condition is violated. In general, this affects the cells in the disc’s corona rather than the cells in the bulk of the disc. In addition to that and in order to stabilise the code, we also introduce a floor in the density, independent of the floor in Alfvén velocity. The condition reads:

We impose a density floor for the dust as well so that ρ_d never decrease below 10⁻¹⁴ in code units. At each time step, we also relax the temperature towards its initial value over a characteristic timescale of 0.1Ω⁻¹, so that both the disc and the corona remain locally isothermal during the course of the simulation (each of them relaxing toward a different temperature or sound speed). We note that this artificial relaxation is similar to a rapid cooling of the flow, which prevents vertical shear instability. Finally to avoid spurious and fast accretion of the field lines near the inner boundary, we relax the poloidal velocity toward zero for R <  15 AU within a characteristic timescale equal to 0.5Ω^−1.

4. Ambipolar simulation without dust

4.1. General properties of the outflow

In this section, we study a global axisymmetric and pure gaseous simulation with β = 10³ and Am₀ = 1. This simulation is integrated over 1000 T₀, where T₀ is the orbital period at the inner radius. In Fig. 1 (left), we show the time-averaged density and streamlines in the inner region between R = 10 and R = 100 AU. The upper (respectively lower) base of the disc is defined by the surface at angle θ₋ = π/2 − θ_d (respectively θ₊ = π/2 + θ_d ) with θ_d ≃ 10°. They are delimited by the dashed black lines in Fig. 1 and correspond to the frontiers at z ≃ 3.5H between the ambipolar-dominated region and the fully ionised atmosphere. This angle corresponds also to the jump in temperature, meaning the separation between the cold disc and the hot corona. We clearly see that a wind is emitted and accelerated from the base of the disc at θ ≃ π/2 ± θ_d, with supersonic speed higher in the atmosphere. The total mass outflow rate, averaged in time and integrated between the inner and outer radius is Ṁ_w(r_out) ≃ 6 × 10⁻⁸ M_⊙ yrs⁻¹. One particularity is that the wind is highly dissymmetric around the midplane, with a more massive outflow in the northern hemisphere ( being four times larger than in the southern hemisphere). This is similar to the configuration depicted in Figs. 21 and 22 of Béthune et al. (2017). Figure 2 (top) shows the vertical profile of the radial and vertical velocities, integrated in time and in radius between R = 25 and R = 50 AU. We find that the outflow in the northern hemisphere, which carries more mass, appears slower than the wind in the southern hemisphere. Figure 1 (right panel) shows the topology of the averaged magnetic field threading the disc, while Fig. 2 (bottom) displays the vertical profile of the radial and toroidal components. In the inner regions (≲50 AU), the field lines cross the midplane with some inclination and bend outside in the southern hemisphere, around z ≃ 3.5H. A remarkable feature is that the toroidal magnetic field keeps its polarity and is almost constant within the midplane. There is a transition between a negative toroidal field and a positive field at R ≃ 60 AU. Another important feature that we discuss later in Sect. 4.3 is the formation of rings, which are visible in Fig. 1 (top) and appear after 200 inner orbits in the simulation.

Fig. 1. Left panel: average gas density and streamlines in the poloidal plane. Right panel: average toroidal field and poloidal field lines. All quantities are averaged between t = 250T0 and t = 1000T0.

Fig. 2. Top panel: vertical profiles of the radial and vertical velocity. Bottom panel: vertical profiles of the radial and toroidal magnetic fields. Quantities are averaged between R = 25 AU and R = 50 AU and in time between t = 250T0 and t = 1000T0.

4.2. Angular momentum budget

We find that the accretion rate in the disc, averaged in time, follows Ṁ_a = 9 × 10⁻⁸(R/10 AU)^0.6 M_⊙ yrs⁻¹. This is roughly five times larger than the cumulative mass outflow rate at 100 AU. The dependence on R of the mass accretion rate is expected in this situation as a result of mass conservation, and is similar to the one found in previous models (Béthune et al. 2017). It is a known signature of weakly magnetised magneto-thermal winds. To characterise the origin of this accretion flow, in Fig. 3 we show the two different terms appearing in the angular momentum budget (Eq. (23)). We find that the vertical transport of angular momentum, associated with the wind stress 𝒯_θϕ is the main driver of accretion. This result is similar to that of Béthune et al. (2017), although in our case the radial transport is not completely negligible, probably because the surface density profile we use here is different from that of Béthune et al. (2017). We checked that the vertical flux of angular momentum is dominated by its Maxwell component −B_zB_ϕ. Because the B_z and B_ϕ field keep the same sign within the disc, the vertical flux of angular momentum ℳ_zϕ = −B_zB_ϕ is unidirectional in the disc bulk. The angular momentum is extracted from the disc southern surface where the gradient ∂_zℳ_zϕ is maximum (i.e. where B_ϕ changes its sign) and is transported towards the northern hemisphere.

Fig. 3. Time evolution of the radial (blue) and vertical (orange) torques, integrated radially and vertically, appearing in the right hand side of the angular momentum equation (Eq. (23)).

In Fig. 4, we also show the radial profile of the vertically averaged radial transport coefficients α_R and α_M. This is averaged between 250 and 1000 T₀. The radial transport of angular momentum is mainly due to the Maxwell stress and settles around α_M ≃ 10⁻², a value that has also been obtained in 3D shearing box ambipolar simulations (see Fig. 12 of Riols & Lesur 2019). Figure 4 shows that the laminar component of the Maxwell stress α_ML contributes to a large fraction of the total stress (see Eq. (25)), like in the box simulations. This term is actually associated with the laminar magnetic field that supports the wind structure. We find, however, that the turbulent stress is not negligible in the simulation, raising the question of the origin of such turbulence. Actually, despite the strong ambipolar diffusion, it can be seen that the disc is still linearly unstable to magnetic perturbations (a form of damped MRI). The criterion for instability, given by Kunz & Balbus (2004), for radially independent fluctuations and an unstratified disc is:

(26)

where k_z is the vertical wavenumber and V_A is the vertical Alfven speed. The marginal case corresponds to k_z = k_zmin ≃ 1/H, which gives the condition:

(27)

Fig. 4. Profile in R of the radial transport coefficient αR and αM (defined in Eq. (24)) averaged in time between 250 and 1000 orbits. The dashed line shows the laminar contribution to these coefficients.

In our simulations, β ranges from ten in the gaps to 10⁴–10⁶ in the ring regions. The condition of instability is marginally satisfied in the gap regions and largely verified in the ring regions with Am_c ranging from 0.008 to 0.0008. However, this criterion assumes that the field is purely vertical. The condition (27) can be generalised to a disc with a horizontal component of the magnetic field. In that case, the critical Am_c has to be multiplied by a factor V_A/V_Az where V_A is the total Alfven speed. In simulations, we typically find that the critical Am, taking into account this factor, is slightly less than one in the ring regions, allowing large-scale unstable modes to develop.

4.3. Ring structures and origin

As already mentioned in Sect. 4.1, the disc develops radial density variations, meaning rings and gaps, associated with zonal flows. These self-organised structures are not new in these types of simulations and were actually found in a variety of non-ideal MHD flows threaded by a net vertical field: Hall-dominated turbulence (Kunz & Lesur 2013), ambipolar-dominated flows (with or without ohmic diffusion, see Simon & Armitage 2014; Bai 2015), and flows combining all non-ideal effects (Bai 2015; Béthune et al. 2017). To look into these structures in more detail, in Fig. 5 (top) we show the time-averaged surface density as a function of radius. The two first ring-gap structures in the inner disc are the most prominent ones and reach amplitudes ΔΣ ≃ 150% of the mean background density Σ₀. The following ring-gap structures are fainter and their amplitudes are about 50% to 60% of the background density. Their typical radial size is about five times the disc scale height. Within the simulation time (∼1000 orbits), these structures remain almost steady. Figure 6 shows the configuration of the poloidal magnetic field lines in the inner disc at t = 350T₀, while Fig. 5 (bottom) shows the time and vertically averaged B_z as a function of radius. We can clearly see that the vertical magnetic field is concentrated into very thin shells at the location of the gaps (minima of Σ). Outside these shells, the net vertical flux is close to zero. This behaviour is reminiscent of other MHD simulations (ideal or not, see Bai & Stone 2014; Bai 2015; Suriano et al. 2018; Riols & Lesur 2018) and the origin of the magnetic shell seems clearly connected to the formation of the ring.

Fig. 5. Surface density Σ defined as (top) and vertically-integrated magnetic field (bottom) as a function of radius R. Quantities are averaged in time between 250 and 1000 orbits, and vertically between θ = θ− and θ = θ+ (z ≃ ±3.5H). The dashed red lines are the location of the density maxima.

In Fig. 7, we also show the instantaneous poloidal streamlines around the ring structures at t = 350T₀. For z ≲ 3.5H (white dashed lines), the flow is mainly radial (accreting in the southern hemisphere and decreting in the northern hemisphere) and seems to follow the rings, moving in a wave-like manner. In particular, it plunges toward the midplane at the gap location. This is line with the poloidal velcocity measured by Teague et al. (2019) in ¹²CO lines in the disc surface around HD 163296. Above this altitude, the flow is mainly vertical (pointing outward) and dominated by the wind component.

One important question is to understand the origin of these structures (rings + magnetic shell) when ambipolar diffusion is the dominant non-ideal MHD effect. Recently, Riols & Lesur (2019) proposed that the rings appear spontaneously in wind-emitting discs and are the result of a linear and secular instability, driven by the MHD wind. The process works as follows: a small initial radial perturbation of density generates a radial flow directed towards the gaps, because of the radial viscous stress. The magnetic field, initially uniform, is radially transported towards the gaps (by the radial flow) and the excess of poloidal flux induces a more efficient mass ejection in the gaps, which reinforces the initial density perturbation. In this instability model, the gas is then depleted vertically by wind plumes rather than being radially accumulated into rings. This process, however, has only been shown in a local shearing box model, and its persistence in global configurations remains to be demonstrated. In particular, global models include a net vertical stress and a net radial accretion velocity that were not present in the local framework.

We then explore how these rings form and evolve in the global configuration. Firstly, in the snapshots of Figs. 6 and 7, it seems obvious, in particular in the southern hemisphere, that the wind is not uniform and homogeneous in the radius. Actually, it is composed of wind plumes that emanate from the regions where the poloidal field is maximum (gap regions). We checked that streamlines connected to the gap regions have the largest vertical (outward) velocity component. In addition to that, in Fig. 8 (top) we show the radial divergence of the flow , averaged in time and in the vertical direction in the disc, as a function of the radius. We see that the divergence is positive within the rings and negative within the gaps. This implies that on average there is a converging radial flow oriented toward the gaps in the disc. In the bottom panel of Fig. 8, we also show that the mass loss rate computed at z = ±5H (in the wind region) is maximum at the gap location, at least for the first three gaps, so the wind removes material from the regions where the vertical magnetic flux is concentrated. All these elements indicate that the mechanism studied by Riols & Lesur (2019) in the local framework also works in the global configuration.

Fig. 6. Snapshot of the simulation for β = 103 at t = 350T0, zoomed in the inner part of the disc. The colour map represents the density, while the white lines are the poloidal magnetic field lines.

Fig. 7. Snapshot of the simulation for β = 103 at t = 350T0, zoomed in on the inner part of the disc. The colour map represents the density, while the black lines are the poloidal streamlines. The dashed white lines delimit the z = ±3.5H surface.

Fig. 8. Top panel: divergence of the radial flow integrated vertically between θ = θ− and θ = θ+ (z ≃ ±3.5H) as a function of radius R. Bottom panel: mass loss rate as a function of radius R (computed at z ≃ ±5H). Quantities are averaged in time between 400 and 650 orbits. The dashed red lines are at the location of the density maxima.

One major difference in the global configuration is that there is a net accretion flow (with z_d = 3.5H), which varies from 0.002 v_K0 at R = 15 AU to 0.01 v_K0 at R = 100 AU. This accretion flow is located mainly at the disc’s southern surface (see Figs. 2 and 7) and should transport the magnetic flux radially over several AU. However, we do not see any motion of the magnetic structures during the course of the simulation; they rather stay at a fixed position. To better understand this, we examine the induction equation for

(28)

Figure 9 shows the radial profile of each term in the right-hand side of the induction equation. We see that the term related to the vertical stretching of the field plays little role in the evolution of . The latter is rather controlled by the two other terms, meaning radial transport by u_r and ambipolar diffusion. Clearly, Fig. 9 shows that these two terms cancel each other out, the first one tends to transport the flux inward ( to the left of the B_z maxima), while the ambipolar term seems to diffuse the field outward ( to the right of the B_z maxima). Therefore, ambipolar diffusion plays a key role in the magnetic flux transport and seems to prevent the magnetic shells associated with the ring structures from being accreted inwards. We note that the mass is mainly accreted at the disc surface, so the density-weighted radial velocity remains small. Hence, rings are expected to be advected at the mass accretion speed, which is much smaller than the gas velocity at the disc surface. Because our models are integrated only for 1000 T₀, this effect cannot be directly checked in the simulations.

Fig. 9. Radial profiles of the ideal and non-ideal terms in the induction equation for Bθ ≃ −Bz (Eq. (28)), integrated vertically between θ = θ− and θ = θ+ (z ≃ ±3.5H) and averaged in time between t = 300 and t = 500T0. The red dashed lines correspond to the local maxima of Bz (averaged vertically and in the same time window).

Finally, we remind the reader that our simulations are axisymmetric and are thus unable to capture non-axisymmetric effects such as the Rossby wave instability (RWI), which might be important in the evolution of the pressure maxima. We of course checked that minima of potential vorticity (∇×v)_z/Σ are located in the ring regions and could in principle excite the RWI (Lovelace et al. 1999). Nevertheless, the axisymmetry and stability of the radial structures seems to be preserved in 3D non-axisymmetric simulations (Béthune et al. 2017; Suriano et al. 2019). One possibility is that the RWI is stabilised by the turbulent stress or by the strong toroidal field (Yu & Li 2009; Gholipour & Nejad-Asghar 2015). Whether these stabilising effects are important or not in 3D simulations remains, however, to be demonstrated.

5. Dust dynamics

5.1. Initialisation

We now explore the dust dynamics in the non-ideal windy MHD flows presented in Sect. 4. We performed simulations for different β and different grain sizes ranging from a = 100 μm up to a = 1 cm. This corresponds to Stokes number St₀ (at z = 0 and R = 10 AU) from 0.001 to 0.1 (see Sect. 2.2.2 for the conversion between sizes and Stokes numbers). In total, we ran four different simulations: the two first at β = 10⁵ and β = 10⁴ both contain 0.1, 0.3, 1, and 3 mm grains (respectively, St₀ of 0.001, 0.003, 0.01 and 0.03). The third simulation is for β = 10³ and contains 0.1, 0.3, 3, and 10 mm grains (respectively, St₀ of 0.001, 0.003, 0.03 and 0.1). The last simulation at β = 10³ contains a unique species of 1 mm in size (St₀ = 0.01).

The initial conditions for the dust are similar for all β and Stokes numbers. The dust velocity is initially Keplerian with zero perturbation. The density at t = 0 is

(29)

with the initial H_d/H_g = 0.5. ρ_d0 being fixed so that the ratio of surface densities Σ_d/Σ is 0.0025. In most of the simulations, we simultaneously integrate four different species so that the total dust-to-gas density ratio is 0.01. We note that for simplicity we use the same initial mass distribution for each species, which is not likely to be the case in real protoplanetary discs. However, since the gas-to-dust ratio remains small when summed over all sizes, this approximation has little effect on the results, as we can later re-normalise the dust size distribution to whatever is needed.

5.2. Settling and dust scale height

At the beginning of each simulation, large dust particles fall towards the midplane within a time proportional to Ω⁻¹/St (Dullemond & Dominik 2004). This is the settling phase. As soon as turbulent diffusion and mixing become significant, dust particles stop settling towards the midplane, and their mean vertical distribution reach an equilibrium. As we showed in Sect. 2.2.2, the midplane Stokes number St depends on R^1/2, so the settling time is proportional to R within the disc midplane. For a 1 mm particle (St₀ = 0.01), the typical settling time is about at R = 10 AU and at R = 200 AU (the disc’s outer radius). This remains smaller than the computation time ( or 650 inner orbits). However, for smaller particles, typically with St₀ = 0.001, the settling time is longer than the simulation time at large radii, so the vertical dust distribution at these locations might not be in equilibrium and may still be evolving slowly in time.

Figure 10 shows the typical dust density distribution in the disc averaged between 400 and 650 inner orbits, for β = 10³ and two different particle sizes corresponding to a = 100 μm (St₀ = 0.001) and a = 3 mm (St₀ = 0.03). Clearly, the dust layer is thinner for the highest Stokes number, which is physically expected, since large particles are more sensitive to gravitational settling. From these density maps, it is possible to calculate the typical dust scale height as a function of the radius R. To be quantitative, we define the dust scale height H_d (St₀, β, R) as the altitude z, in such a way that

(30)

Fig. 10. Dust density distribution for grains of size 100 μm (left) and 3 mm (right), for β = 103. This is averaged in time between 400 and 650 inner orbits

We find that the dust to gas scale height ratio H_d/H_g is a decreasing function of radius R. This is due to the fact that the effective Stokes number increases with R^1/2 (see Eq. (12)). We note that the power law depends on the choice of the disc’s surface density profile and its dependence on R. For instance, for β = 10⁴ and a 1 mm particle, we find that the dust layer size is 0.26H_g (0.13 AU) at R = 10 AU, and 0.18H_g (0.9 AU) at R = 100 AU.

For each β and particle size, we divide the radial domain into ten bins (or parcels) and calculate the average ratio H_d/H_g and average Stokes number within each parcel of the disc. In Fig. 11, we show the resulting H_d/H_g as a function of the average ‘effective’ Stokes number, where each colour accounts for a different particle size (or fixed St₀). As expected, we find that the ratio H_d/H_g decreases with the Stokes number, and we end up with the following scaling relations:

(31)

(32)

(33)

Fig. 11. Dust to gas scale height ratio as a function of the effective Stokes number for β = 105 (left), β = 104 (centre) and β = 103 (right). Each colour accounts for a given particle size (purple: 1 cm, blue: 3 mm, yellow: 1 mm, red: 0.3 mm and green: 0.1 mm.)

We note that for the largest particle size (St₀ = 0.1 for β = 10³ or St₀ = 0.03 for β = 10⁵), we lack resolution since the dust is contained within two grid cells around the midplane. In that case, the ratio H_d/H_g is probably overestimated. This could explain why the blue points in the left panel of Fig. 11 lie above the expected scaling law (dashed line).

The dependence on St^−1/2 is actually reminiscent of many other disc simulations coupling dust and gas turbulence (Fromang & Papaloizou 2006; Zhu et al. 2015; Riols & Lesur 2018). It can be obtained using the classical diffusion theory (Morfill et al. 1985; Dubrulle et al. 1995) in which the equilibrium distribution in the vertical direction results from the balance between the gravitational settling and the diffusion due to the turbulence. In that theory, the turbulence is assumed to be homogeneous and isotropic, and can be described through a uniform, effective diffusion coefficient. We note that the local simulations of Riols & Lesur (2018) found a deviation from the simple power law St^−1/2 at β = 10³ and for small particles. In our case, we do not see such a deviation. Moreover, for β = 10⁴ and β = 10³, the ratio H_d/H_g is two times larger than those measured in the local shearing box (see comparison with Fig. 7 of Riols & Lesur 2018). This is due to a larger vertical rms turbulent velocity measured in the midplane in global simulations (v_zrms ∼ 0.07c_s) compared to local simulations (where v_zrms ∼ 0.04c_s, see Fig. 2 of Riols & Lesur 2018 for β = 10³).

5.3. Concentration in pressure maxima

As we showed in Sect. 4.3, rings of matter arise spontaneously in wind emitting discs, and they correspond to pressure maxima in the locally isothermal case. These maxima are then the locations for dust accumulation and could play an important role in the formation of planetesimals. The concentration of dust is, however, only possible if the drift velocity, proportional to St and the local pressure gradient, is larger than the velocity associated with the radial turbulent diffusion. We then expect concentration of dust into radial bands to be more efficient for sufficiently large values of the Stokes number. In Fig. 12, we show the radial profiles of the gas pressure and dust surface density Σ_d, averaged in time, for two different plasma beta parameters and two different grain sizes. For β = 10³, both large and small grains concentrate into pressure maxima (with position delimited by the red dashed vertical lines). However, 3 mm grains, the contrast in density, at some radii, can be 200–500 times larger than the contrast measured for 0.1 mm grains. Although the density contrast is pretty high for a = 3 mm, the total dust-to-gas ratio remains bounded and smaller than 0.3 in most of the rings. We also note that the radial dust structures are thinner in the case of larger particles. For β = 10⁴, the results are similar, except that the pressure maxima are less pronounced and the dust is less concentrated on average. This is particularly obvious for 0.1 mm grains, for which the contrast in density never exceeds ten. Also, the typical separation between the dust rings in the case β = 10⁴ is about 2.5H at R ≃ 20 AU, while it is about 4H for β = 10³. By fitting the radial density profiles with gaussian functions, we find that the half-width of the dust rings at R = 20 AU is ∼0.3−0.4 AU (0.3–0.4H) for β = 10³ and a = 3 mm, while it is ∼0.15 AU for β = 10⁴. For comparison, the pressure bumps in the gas have a typical half-width of 1.2H.

Fig. 12. Top panels: gas pressure as a function of radius R for β = 103 (left) and β = 104 (right). Central and bottom panels: dust surface density for grains of 3 mm and 100 μm, respectively. Quantities are averaged vertically between θ = θ− and θ = θ+ (z ± 3.5H) and in time between 400 and 650 inner orbits.

We note that for β = 10⁵, a few rings form, but they are much fainter than in the case β = 10⁴ or β = 10³, and therefore we do not see any clear radial structures in the dust. We cannot rule out that more prominent structures could form if the simulation were run for a longer time period.

6. Radiative transfer with MCFOST

The aim of this Section is to produce synthetic images derived from our simulations in order to predict the observable properties of ambipolar-dominated discs threaded by a net vertical field. For that we use the 3D continuum and line radiative transfer code MCFOST based on the Monte Carlo method (Pinte et al. 2006).

Assuming a given gas and dust density distribution, the code first computes the temperature structure assuming that the dust is in radiative equilibrium with the local radiation field. This is done by propagating photon packets, originally emitted by the central star (assumed to be of solar-type), through a combination of scattering, absorption and re-emission events until they exit the computation grid. During this step, we use 1 280 000 photon packets in total. It is assumed that the dust opacities are fixed and do not depend on the temperature during this stage. Observables’ quantities (images) are then obtained via a ray-tracing method, which calculates the intensities by integrating the radiative transfer equation along rays perpendicular to each pixel of the image.

The disc extends from an inner cylindrical radius R_in = 10 AU to an outer limit R_out = 200 AU, with a total gaseous mass equal to 5 × 10⁻² solar mass. Dust grains are defined as homogeneous and spherical particles (Mie theory) with sizes distributed according to the power law dn(a)/da = a^p with p = −3.5. During the radiative transfer calculation, we consider 100 different sizes of grains ranging from 0.03 μm to 1 cm. We note that our PLUTO simulations only contain four dust species of 0.1, 0.3, 3 and 10 mm, which all possess the same initial mass. Thus, to construct the distribution of each grain population, MCFOST automatically interpolates the axisymmetric PLUTO density fields between the input sizes, and then normalises them to make sure that dn(a)/da = a^p with the additional constraint that the gas-to-dust ratio is 100. Such a re-normalisation procedure can be achieved because the dust has almost relatively weak feedback on the gas (since the dust-to-gas ratio remains relatively low (< 0.3) in all simulations). We note that during the interpolation, we assume that the smallest grains of 0.03 μm follow the gas distribution exactly.

In the upper panels of Fig. 13, we show the image of our disc for β = 10³ calculated at a given wavelength (λ = 1 mm), corresponding to ALMA band 7, in a face-on (left) and edge-on (right) configuration. To take into account the resolution of instruments like ALMA, in the lower panels we show the same images convolved with a Gaussian kernel with a full width half maximum (FWHM) of 3 AU. This resolution is typically expected for an object located at 140 pc and for a radio-telescope baseline of 15 km. In the original image, between ten and 50 AU we see a series of eight concentric rings that correspond to the location of the pressure maxima in Fig. 12 where the millimetre dust is trapped. In the convolved image, these rings are reduced to one bright inner ring and a second fainter ring (see also top panel of Fig. 14). A wide gap is obtained at 50 AU, followed by other bright rings in the outer regions (in particular one at R = 70 AU and another at R = 100 AU in the convolved image). The typical contrast at 50 AU in the convolved image is 7–8, while the separation between the gaps or the rings is close to 20 AU. In the edge-on configuration, the disc thickness is about 5 AU at R = 200 AU, which corresponds to 0.5 H_g. For λ = 1 mm, most of the thermal emission comes from grains of λ/2π ≃ 160 μm, which are indeed contained in a layer of ≈0.5H_g at R = 200 AU according to Fig. 11. We note that the emission is lower in the midplane of the disc than in the layers above (see Appendix A showing the vertical profile of the flux density). This is due to the smaller temperature in the midplane, resulting from the strong sedimentation of mm to cm grains.

Fig. 13. Top: flux density νFν in W m−2 pixel−1 calculated by MCFOST, with the disc viewed face-on (left panel) and edge-on (right panel) for λ = 1 mm. This is calculated for t = 500T0 and β = 103. Bottom: same images but convolved with a Gaussian kernel with a full width half maximum (FWHM) of 3 AU. The white dashed lines in the right panels shows the surface z = Hg for R = 200 AU.

We also show, in the bottom panel of Fig. 14, the spectral index α_s = d log F(ν)/d log ν as a function of the radius obtained by comparing the image at λ = 1 mm with the image at λ = 3 mm. The spectral index ranges between 2.5 and 3.5, suggesting that the disc deviates slightly from a black body (whose spectral index is 2 in the Rayleigh-Jeans regime). The spectral index is found to be at a maximum in the gap regions and at a minimum within the rings, in agreement with the fact that the surface density, and therefore the optical thickness, is lower in the gaps.

Fig. 14. Top panel: convolved flux density νFν in W m−2 pixel−1 as a function of radius. Bottom panel: spectral index αs = d log F(ν)/d log ν as a function of radius measured between λ = 1 mm (ALMA band 7) and λ = 3 mm (ALMA band 3).

Finally, we computed the scattered light emission from micron-size dust at the mid-infrared wavelength λ = 1.65 microns, which corresponds to band H of the SPHERE instrument. The synthetic images of the discs viewed face-on and edge-on are shown in the top panels of Fig. 15. In the face-on view, we distinguish several ring structures, although they remain much fainter than in the ALMA bands. In the edge on view, the emission comes mainly from the surface z ≃ 2H. Since the disc is optically thick, the infrared emission from the midplane region is totally absorbed, allowing only the detection of scattered light from the outer edges of the upper and lower surfaces of the disc. To take into account the beam of the instrument, we convolve the original images (top panels) with a Gaussian kernel of FWHM of 10 AU (bottom panels). We find that rings are absent in the convolved images, since their actual size is smaller than the beam size.

Fig. 15. Top: flux density νFν in W m−2 pixel−1 calculated by MCFOST, with the disc viewed face-on (left panel) and edge-on (right panel) for λ = 1.65 microns. This is calculated for t = 500T0 and β = 103. Bottom: same images but convolved with a Gaussian kernel with a full width half maximum (FWHM) of 10 AU. The white dashed lines in the right panels shows the surface z = Hg for R = 200 AU.

7. Conclusion and discussion

In summary, we investigated the gas and dust dynamics of protoplanetary discs threaded by a large-scale magnetic field and subject to ambipolar diffusion. Using radially and vertically stratified global axisymmetric MHD simulations, we showed that accretion occurs mainly at the disc surface and is mainly driven by the laminar torque associated with large-scale MHD winds. For plasma β of 10³, the accretion rate is about 10⁻⁷ M_⊙ yr⁻¹ at R = 10 AU. We note that the mass loss rate represents a significant fraction of the mass accretion rate, so that d log Ṁ_acc/d log R ≃ 0.6, typical of magneto-thermal winds. Hence, if we assume winds exist down to the disc’s inner radius R_in ≃ 0.1 AU, one expects Ṁ_acc(R_in) ≲ 10⁻⁸ M_⊙ yr⁻¹, which is fully compatible with observed accretion rates onto YSOs. Within the midplane, the disc displays gaseous ring-gap structures associated with zonal flows, with a segregation between vertical magnetic flux and matter. These structures are stable during the simulation and are likely to be triggered by the same mechanism described in Riols & Lesur (2019), which involves a wind instability.

We showed that centimetre to millimetre dust particles are highly sedimented in the midplane with a typical scaleheight H_d/H_g ≲ 0.25 for β = 10⁴. For R = 100 AU, this corresponds to a dust layer of size ≲1 AU. This result fits the observational constraint of Pinte et al. (2016) in HL Tau (ALMA) particularly well and is in contrast with ideal MRI turbulence simulations for which the dust-to-gas scale height ratio is much larger than that inferred from observations (Johansen et al. 2006; Fromang et al. 2006).

Another result is that the dust particles are trapped within the pressure maxima and form thin axisymmetric structures with high contrast (for mm to cm particles) and typical separation of 2.5–5 H_g at R = 20 AU, wider separation being obtained in stronger field situations. The width of the millimetre dust rings is about 0.3−0.4H for β = 10³. This is comparable with the width of dust ring structures observed by the DSHARP survey (Dullemond & Penzlin 2018) in objects like AS 209 (although the comparison is done at different radii, and we implicitly assume that the ratio width over H does not depend on radius). We estimate that the typical dust mass contained inside a single ring at R ≃ 20 AU is about 5 Earth mass. Using the radiative transfer code MCFOST, we were able to produce synthetic images of the disc in the ALMA band (1 mm) for β = 10³. These images, when convolved with a Gaussian kernel that accounts for the instrumental resolution, show concentric rings and gaps with typical contrast ≲10 and separation of 20 AU. The spectral index within the gaps goes up to 3.5, while it remains below 3 within the rings, similarly to observed spectral indices in HL-Tau (ALMA Partnership 2015). At a mid-infrared wavelength (scattered light emission), we found that instruments like SPHERE would be unable to detect the rings, due to their insufficient resolution.

We note that our simulations are limited to the region between 10 and 200 AU, but we expect the instability giving rise to the ring-gap structures to still work in the inner and densest regions where ohmic diffusion predominates (see Riols & Lesur 2019). This could explain, for instance, the dust rings observed in the TW Hydra disc on scales smaller than 10 AU (Andrews et al. 2016).

It is relevant to compare our work with recent simulations of ambipolar-dominated flow. In particular, Suriano et al. (2018) showed that the segregation of vertical magnetic flux and matter is due to the reconnection of highly pinched poloidal fields. This requires a current sheet in the midplane and a poloidal field that bends within it. In our case, we only see pinched poloidal fields at the very beginning of the simulations (t <  100T₀) when the rings are not yet formed. At later times, the field takes the configuration depicted in Fig. 6, similar to what was obtained by Béthune et al. (2017). This configuration shows a break of symmetry, inappropriate regarding reconnection in the midplane. The cause of the breaking of symmetry in our case is yet unknown, but appears to be a robust result independent of the initial condition. Another comparison can be made with the recent local box simulations of Riols & Lesur (2018) who computed the dust and gas motions around a fixed radius R = 30 AU. We found that the properties of the gaseous and dusty rings (separation and depth) are comparable in the local and global frameworks. However the settling of the dust is more pronounced in the shearing box in particular for β = 10⁴ and β = 10³.

We remind the reader that the choice of cooling timescale adopted in our simulations prevents the development of the vertical shear instability (VSI). Recent radiation hydrodynamic simulations (Stoll & Kley 2014, 2016; Flock et al. 2017) suggest that the VSI should actually be present in the outer regions of protoplanetary discs and induce substantial diffusion of dust particles. We argue, however, that the physics of the VSI strongly depends on the thermodynamics, and that the flux-limiting diffusion used by the previous authors might not be a good approximation away from the midplane where the VSI triggers the strongest motions. Moreover, a recent study by Cui & Bai (2020) investigating the interplay between MHD winds in PP discs with the VSI shows that the hydrodynamic instability is weakened by MHD effects in the moderate magnetisation regime (β = 10³).

This work is also connected to planet formation theory. The formation of gaseous rings may indeed be a way to overcome the ‘radial-drift’ barrier (growth of grains to pebbles) and the ‘fragmentation’ barrier (growth from planetesimals to planets) (Birnstiel et al. 2016). In particular, the trapping of grains in pressure maxima can hinder their radial migration, accelerate their growth, and prevent their fragmentation (Gonzalez et al. 2017). Gaseous rings associated with zonal flows appear to be an alternative to other trapping mechanisms such as vortices, whose stability and origin remain uncertain (Lesur & Papaloizou 2009). An interesting avenue of research would be to better characterise the grain growth inside these ring structures. Also, the strong sedimentation in the midplane, especially for β = 10⁵ or β = 10⁴ could also reinforce the so-called pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012), and therefore speed up the growth of massive cores.

Our results suggest that self-organisation is unavoidable in weakly ionised discs threaded by a large-scale poloidal field, provided that the field is strong enough. The question of the strength of the field threading these discs is therefore of central importance. Unfortunately, there is no direct measurement of the field strength at these distances. This is a rather difficult task since the field strength corresponding to a β = 10³ disc is approximately 15 mG at 10 AU, which is a relatively weak field. On the theory side, the most recent core-collapse models (Masson et al. 2016) suggest that discs form with a constant large-scale field of 100 mG, meaning stronger but still comparable to our β = 10³ disc model. As the disc evolves, the field strength is also expected to evolve, but no theory is yet able to predict whether the field should diffuse away, or on the opposite accumulates towards the centre (see e.g. Guilet & Ogilvie 2013; Bai 2017). In any case, this suggests that β ≃ 10² − 10⁴ is well within the large-scale field strength expected for protoplanetary discs, and therefore that wind-driven self-organisation is most probably present in these objects.

In the end, the main question is to unambiguously identify which mechanism is responsible for the structures we see, whether it is induced by embedded planets, wind-driven self-organisation, or other secular gas instabilities. As pointed out in the introduction, finding planets in the gaps is not sufficient, since all of the mechanisms will accelerate planet formation in the rings thanks to the presence of local pressure maxima that trap growing grains. The presence of a deviation from the Keplerian rotation profile is not really useful either, since this is expected as a result of the radial geostrophic equilibrium of the disc, which is always satisfied, with or without planets. More specific signatures are therefore desperately needed. A first potential signature would be the eccentricity of the gap/ring structures, which is expected to be zero in the case of wind-driven self-organisation, but to exceed 0.1 if a planet of mass ≳0.1 Jupiter mass opens the gap (Zhang et al. 2018). Such eccentricity might be difficult to measure with current facilities but feasible with the next generation of instruments. In the particular case of wind-driven self-organisation, the kinematics of the gas is peculiar, since the gas tends to radially converge towards the gap, before being ejected. At the disc surface, which is typically tracked by CO observations, the gas is typically flowing radially, “bouncing” onto the rings (Fig. 7). This is the kind of kinematics that can be tested on the sky (see e.g. Teague et al. 2019).

Our global simulations also predict that the typical separation between the rings increases with the disc magnetisation, a result that has been demonstrated, as well in the shearing box (Riols & Lesur 2019). We also expected the accretion rate and mass loss rate in the wind to increase with the magnetisation. If the rings are formed via MHD winds, we would then expect a correlation between the ring separation and the accretion rate. This relation could be tested in the current and future observations.

Ideally, one would also want to look for the presence of magnetic field concentration in the gaps. We expect a contrast in field intensity reaching several tens, implying a field strength of a few hundreds mG at 10 AU in the gaps. This might not be achievable with current instruments, but such a detection would undoubtedly demonstrate the reality of this kind of mechanism.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 815559 (MHDiscs)). This work was performed using the Froggy platform of the CIMENT HPC infrastructure (https://ciment.ujf-grenoble.fr).

References

Appendix A: Visibility in ALMA bands

In Fig. A.1, we show the normalised vertical profiles of the flux density νF_ν computed in different ALMA bands with MCFOST in the case of a disc viewed edge-on. These profiles are calculated at the centre of the image (x = 0) and are normalised to their maximum (longer wavelengths have lower emissions). We see that the width of the profiles decreases with the wavelength, since a longer wavelength traces a larger, more sedimented particle size. In the midplane region at z = 0, we observe a gap in the emission corresponding to a darkening of the disc. The depth of the gap increases at lower wavelength. The darkening is due to large sedimented grains (of cm size) absorbing the light and obscuring the disc midplane. They also contribute to lower the temperature in the midplane, and thus reducing the emission from this region.

Fig. A.1. Vertical profiles of the flux density νFν for different wavelengths corresponding to different ALMA bands, computed with MCFOST in case of a disc viewed edge-on. The flux density is normalised to its maximum.

All Figures

	Fig. 1. Left panel: average gas density and streamlines in the poloidal plane. Right panel: average toroidal field and poloidal field lines. All quantities are averaged between t = 250T0 and t = 1000T0.
In the text

	Fig. 2. Top panel: vertical profiles of the radial and vertical velocity. Bottom panel: vertical profiles of the radial and toroidal magnetic fields. Quantities are averaged between R = 25 AU and R = 50 AU and in time between t = 250T0 and t = 1000T0.
In the text

	Fig. 3. Time evolution of the radial (blue) and vertical (orange) torques, integrated radially and vertically, appearing in the right hand side of the angular momentum equation (Eq. (23)).
In the text

	Fig. 4. Profile in R of the radial transport coefficient αR and αM (defined in Eq. (24)) averaged in time between 250 and 1000 orbits. The dashed line shows the laminar contribution to these coefficients.
In the text

	Fig. 5. Surface density Σ defined as (top) and vertically-integrated magnetic field (bottom) as a function of radius R. Quantities are averaged in time between 250 and 1000 orbits, and vertically between θ = θ− and θ = θ+ (z ≃ ±3.5H). The dashed red lines are the location of the density maxima.
In the text

	Fig. 6. Snapshot of the simulation for β = 103 at t = 350T0, zoomed in the inner part of the disc. The colour map represents the density, while the white lines are the poloidal magnetic field lines.
In the text

	Fig. 7. Snapshot of the simulation for β = 103 at t = 350T0, zoomed in on the inner part of the disc. The colour map represents the density, while the black lines are the poloidal streamlines. The dashed white lines delimit the z = ±3.5H surface.
In the text

	Fig. 8. Top panel: divergence of the radial flow integrated vertically between θ = θ− and θ = θ+ (z ≃ ±3.5H) as a function of radius R. Bottom panel: mass loss rate as a function of radius R (computed at z ≃ ±5H). Quantities are averaged in time between 400 and 650 orbits. The dashed red lines are at the location of the density maxima.
In the text

	Fig. 9. Radial profiles of the ideal and non-ideal terms in the induction equation for Bθ ≃ −Bz (Eq. (28)), integrated vertically between θ = θ− and θ = θ+ (z ≃ ±3.5H) and averaged in time between t = 300 and t = 500T0. The red dashed lines correspond to the local maxima of Bz (averaged vertically and in the same time window).
In the text

	Fig. 10. Dust density distribution for grains of size 100 μm (left) and 3 mm (right), for β = 103. This is averaged in time between 400 and 650 inner orbits
In the text

	Fig. 11. Dust to gas scale height ratio as a function of the effective Stokes number for β = 105 (left), β = 104 (centre) and β = 103 (right). Each colour accounts for a given particle size (purple: 1 cm, blue: 3 mm, yellow: 1 mm, red: 0.3 mm and green: 0.1 mm.)
In the text

	Fig. 12. Top panels: gas pressure as a function of radius R for β = 103 (left) and β = 104 (right). Central and bottom panels: dust surface density for grains of 3 mm and 100 μm, respectively. Quantities are averaged vertically between θ = θ− and θ = θ+ (z ± 3.5H) and in time between 400 and 650 inner orbits.
In the text

	Fig. 13. Top: flux density νFν in W m−2 pixel−1 calculated by MCFOST, with the disc viewed face-on (left panel) and edge-on (right panel) for λ = 1 mm. This is calculated for t = 500T0 and β = 103. Bottom: same images but convolved with a Gaussian kernel with a full width half maximum (FWHM) of 3 AU. The white dashed lines in the right panels shows the surface z = Hg for R = 200 AU.
In the text

	Fig. 14. Top panel: convolved flux density νFν in W m−2 pixel−1 as a function of radius. Bottom panel: spectral index αs = d log F(ν)/d log ν as a function of radius measured between λ = 1 mm (ALMA band 7) and λ = 3 mm (ALMA band 3).
In the text

	Fig. 15. Top: flux density νFν in W m−2 pixel−1 calculated by MCFOST, with the disc viewed face-on (left panel) and edge-on (right panel) for λ = 1.65 microns. This is calculated for t = 500T0 and β = 103. Bottom: same images but convolved with a Gaussian kernel with a full width half maximum (FWHM) of 10 AU. The white dashed lines in the right panels shows the surface z = Hg for R = 200 AU.
In the text

	Fig. A.1. Vertical profiles of the flux density νFν for different wavelengths corresponding to different ALMA bands, computed with MCFOST in case of a disc viewed edge-on. The flux density is normalised to its maximum.
In the text