Ring formation and dust dynamics in wind-driven protoplanetary discs: global simulations

Large-scale vertical magnetic fields are believed to play a key role in the evolution of protoplanetary discs. Associated with non-ideal effects, such as ambipolar diffusion, they are known to launch a wind that could drive accretion in the outer part of the disc ($R>1$ AU). They also potentially lead to self-organisation of the disc into large-scale axisymmetric structures, similar to the rings recently imaged by sub-millimetre or near-infrared instruments (ALMA and SPHERE). The aim of this paper is to investigate the mechanism behind the formation of these gaseous rings, but also to understand the dust dynamics and its emission in discs threaded by a large-scale magnetic field. To this end, we performed global magneto-hydrodynamics (MHD) axisymmetric simulations with ambipolar diffusion using a modified version of the PLUTO code. We explored different magnetisations with the midplane $\beta$ parameter ranging from $10^5$ to $10^3$ and included dust grains -- treated in the fluid approximation -- ranging from $100 \mu$m to 1 cm in size. We first show that the gaseous rings (associated with zonal flows) are tightly linked to the existence of MHD winds. Secondly, we find that millimetre-size dust is highly sedimented, with a typical scale height of 1 AU at $R=100$ AU for $\beta=10^4$, compatible with recent ALMA observations. We also show that these grains concentrate into pressure maxima associated with zonal flows, leading to the formation of dusty rings. Using the radiative transfer code MCFOST, we computed the dust emission and make predictions on the ring-gap contrast and the spectral index that one might observe with interferometers like ALMA.


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 −10 to 10 −7 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(Flock et al. , 2013Bai & 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;Lesur et al. 2014;Bai 2015).
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. 2015Flaherty et al. , 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 to 70 % of the observed discs (Garufi et al. 2018). Some examples of discs showing rings are HL tau (ALMA Partnership et al. 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 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 largescale 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. 2016Béthune et al. , 2017Suriano et al. 2018a) 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. 2018a,b), 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 Section 2, we describe the model and review the main characteristics of non-ideal physics and dust-gas interaction. In Section 3, we present the numerical methods used to simulate the dust-gas dynamics in a global model. In Section 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 Section 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 Section 6, we produce synthetic images of our disc simulation in ALMA bands by using MCFOST. Finally, in Section 7, we conclude and discuss the potential implications of our work for protoplanetary disc evolution and planet formation.

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.

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 P = ρc 2 s , where c s (R, z) is a fixed sound-speed profile (see Section 2.1.4 for temperature prescription). The evolution of its density ρ, and total velocity field v follows: We note that the unit of time is taken as the inner Ω −1 and the unit of length corresponds to R = 10AU. The unit of mass is such that ρ = 1 in the midplane at R = 10AU. The magnetic field B, which appears in the Lorentz force (∇ × B) × B, is governed by the non-ideal induction equation, 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 Section 2.1.3.

Surface density profile
We assume a disc model with surface density following 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 −1 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 Section 4.2).

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 −9 cm 3 s −1 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 v A = B/ √ ρ, 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 2 . 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 where Am 0 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 0 1. (Thi et al. 2018). We note that to convert numerical values of Σ i into g/cm 2 , we used the disc model of Eq. (5).

Temperature and hot corona
Complex thermo-chemical models of protoplanetary discs (Kamp & Dullemond 2004;Thi et al. 2018) 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 coronato-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. with 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 Ω −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 ρ d k and velocity v d k . 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: Article number, page 3 of 17 with γ g−>d k 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 τ s k 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 τ s k , we consider in this study that dust particles are spherical and small enough 1 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: A useful dimensionless quantity to parametrise the coupling between dust and gas is the Stokes number: St = Ωτ s k .

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 ρ = ρ 0 = 1, c s = c s 0 = 0.05 and Ω = Ω 0 = 1, which we label St 0 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 Σ = ρ 0 H √ 2π, where ρ 0 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 −3 , we obtain the following conversion in the midplane: 1 We note that in principle, dust particles could be porous and have a more complex geometry. This would not change the dynamics described here, but only the relation between the Stokes number and the physical size of the particles (11).

Initial equilibrium
For the initial hydrostatic disc equilibrium, we use the same prescription as Zhu & Stone (2018): where c s = P/ρ is the sound speed that is constant on the cylindrical radius: 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 s 0 /v K 0 . 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 s 0 /v K 0 = 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 0 = 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 :

Averages
We define a vertical integral of a quantity Q between angles θ − and θ + as and, respectively, a radial integral and temporal average:

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: where Σ = ρ is the gas surface density defined on spherical geometry andσ w = ρu θ sin θ θ + θ − > 0 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 θ + : and the vertical stress difference between the two angles: the conservation of angular momentum in spherical coordinates reads: with Ω(r) = GM/r 3 , 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: Multiplying Eq. (21) by r and integrating in radius, we obtain the following conservation equation: Accretion onto the star can be either due to a radial transport of angular momentum, via the radial stress T rφ , or due to angular momentum extracted by a wind through the disc surface, via the vertical stress T θφ . An important quantity that characterises the radial transport efficiency is the parameter α (or Shakura & Sunyaev parameter) which is: where R rφ = ρu φ u r sin θ is the Reynolds stress and M rφ = −B φ B r sin θ is the Maxwell stress. The coefficient α can be decomposed into a sum of a laminar part: plus a turbulent part: α ν = α − α L .

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.

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.

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 K 0 < 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 −14 in code units. At each time step, we also relax the temperature towards its initial value over a characteristic timescale of 0.1Ω −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   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) 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 Article number, page 6 of 17 A. Riols, G. Lesur and F. Menard: Ring formation and dust dynamics in PP discs 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. the inner boundary, we relax the poloidal velocity toward zero for R < 15AU within a characteristic timescale equal to 0.5Ω −1.

General properties of the outflow
In this section, we study a global axisymmetric and pure gaseous simulation with β = 10 3 and Am 0 = 1. This simulation is integrated over 1000 T 0 , where T 0 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 −8 M /yrs. One particularity is that the wind is highly dissymmetric around the midplane, with a more massive outflow in the northern hemisphere (σ w 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 = 50AU. 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 ( 50AU), 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 60AU. Another important feature that we discuss later in  Fig. 5. Surface density Σ defined as ρ (top) and vertically-integrated magnetic field B z (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. Fig. 6. Snapshot of the simulation for β = 10 3 at t = 350T 0 , zoomed in the inner part of the disc. The colour map represents the density, while the white lines are the poloidal magnetic field lines. Section 4.3 is the formation of rings, which are visible in Fig. 1 (top) and appear after 200 inner orbits in the simulation.

Angular momentum budget
We find that the accretion rate in the disc, averaged in time, fol-lowsṀ a = 9 × 10 −8 (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 T θφ 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 z B φ . Because the B z and B φ field keep the same sign within the disc, the vertical flux of angular momentum M zφ = −B z B φ is unidirectional in the disc bulk. The angular momentum is extracted from the disc southern surface where the gradient ∂ z M zφ is maximum (i.e. where B φ changes its sign) and is transported towards the northern hemisphere.
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 0 . The radial transport of angular momentum is mainly due to the Maxwell stress and settles around α M 10 −2 , 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 α M L 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: where k z is the vertical wavenumber and V A is the vertical Alfven speed. The marginal case corresponds to k z = k z min 1/H, which gives the condition: In our simulations, β ranges from ten in the gaps to 10 4 -10 6 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 A z 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.

Ring structures and origin
As already mentioned in Section 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 In Fig. 7, we also show the instantaneous poloidal streamlines around the ring structures at t = 350T 0 . 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 12 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 Fig. 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 1/r ∂/∂r(rv r ), 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σ w 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.
One major difference in the global configuration is that there is a net accretion flow 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 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 1 r ∂ ∂r (−ru θ B r ) plays little role in the evolution of B z . 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 ( 1 r ∂ ∂r (ru r B θ ) > 0 to the left of the B z maxima), while the ambipolar term seems to diffuse the field outward ( 1 r ∂ ∂r (rη A J ⊥ ) > 0 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 0 , this effect cannot be directly checked in the simulations.
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. 2018b). 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.

Initialisation
We now explore the dust dynamics in the non-ideal windy MHD flows presented in Section 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 0 (at z = 0 and R = 10 AU) from 0.001 to 0.1 (see Section 2.2.2 for the conversion between sizes and Stokes numbers). In total, we ran four different simulations: the two first at β = 10 5 and β = 10 4 both contain 0.1, 0.3, 1, and 3 mm grains (respectively, St 0 of 0.001, 0.003, 0.01 and 0.03). The third simulation is for β = 10 3 and contains 0.1, 0.3, 3, and 10 mm grains (respectively, St 0 of 0.001, 0.003, 0.03 and 0.1). The last simulation at β = 10 3 contains a unique species of 1 mm in size (St 0 =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 with the initial H d /H g = 0.5. ρ d 0 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.

Settling and dust scale height
At the beginning of each simulation, large dust particles fall towards the midplane within a time proportional to Ω −1 /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 Section 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 1mm particle (St 0 = 0.01), the typical settling time is about 100 Ω −1 0 at R=10 AU and 2000 Ω −1 0 at R=200 AU (the disc's outer radius). This remains smaller than the computation time ( 4000 Ω −1 0 or 650 inner orbits). However, for smaller particles, typically with St 0 = 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 3 and two different particle sizes corresponding to a = 100µm (St 0 = 0.001) and a = 3mm (St 0 = 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 0 , β, R) as the altitude z, in such a way that 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 4 and a 1mm particle, we find that the dust layer size is 0.26H g (0.13 AU) at R = 10AU, and 0.18H g (0.9 AU) at R = 100AU. 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 0 ). 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: We note that for the largest particle size (St 0 = 0.1 for β = 10 3 or St 0 = 0.03 for β = 10 5 ), 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 1985;Dubrulle et al. 1995) in which the equilibrium distribution in the vertical direction results from the balance between Top panels shows the gas pressure as a function of radius R for β = 10 3 (left) and β = 10 4 (right). The central and bottom panels show the 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. 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 3 and for small particles. In our case, we do not see such a deviation. Moreover, for β = 10 4 and β = 10 3 , 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 z rms ∼ 0.07c s ) compared to local simulations (where v z rms ∼ 0.04c s , see Fig. 2 of Riols & Lesur (2018) for β = 10 3 ).

Concentration in pressure maxima
As we showed in Section 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 3 , both large and small grains concentrate into pressure maxima (with position delimited by the red dashed vertical lines). However, 3mm grains, the contrast in density, at some radii, can be 200 to 500 times larger than the contrast measured for 0.1 mm grains. Although the density contrast is pretty high for a = 3 mm, the dust-to-gas ratio remains bounded and smaller than 0.04. We also note that the radial dust structures are thinner in the case of larger particles. For β = 10 4 , 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 4 is about 2.5 H at R 20 AU, while it is about 4 H for β = 10 3 . By fitting the radial density profiles with gaussian functions, we find that the half-width of the dust rings at R = 20AU is ∼ 0.3 − 0.4 AU (0.3-0.4 H) for β = 10 3 and a = 3 mm, while it is ∼ 0.15 AU for β = 10 4 . For comparison, the pressure bumps in the gas have a typical half-width of 1.2 H.
We note that for β = 10 5 , a few rings form, but they are much fainter than in the case β = 10 4 or β = 10 3 , 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.

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 −2 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 no feedback on the gas (since the dust-to-gas ratio remains low 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 3 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 15km. 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.
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.
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.

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 3 , the accretion rate is about 10 −7 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 expectṡ M acc (R in ) 10 −8 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 4 . 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;. 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 to 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 3 . 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 (1mm) for β = 10 3 . 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 et al. 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. (2018a) 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 0 ) 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 4 and β = 10 3 .
We remind the reader that the the choice of cooling timescale adopted in our simulations prevents the development of the vertical shear instability (VSI). Recent radiation hydrodynamic simulations (Stoll & Kley 2014Flock 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 (2019) 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 3 ).
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 5 or β = 10 4 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 3 disc is