Issue 
A&A
Volume 650, June 2021



Article Number  A35  
Number of page(s)  13  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202040109  
Published online  02 June 2021 
Systematic description of winddriven protoplanetary discs^{★}
Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble,
France
email: geoffroy.lesur@univgrenoblealpes.fr
Received:
10
December
2020
Accepted:
24
January
2021
Aims. Planetforming discs are believed to be very weakly turbulent in the regions outside of 1 AU. For this reason, it is now believed that magnetised winds could be the dominant mechanism driving accretion in these systems. However, currently, no selfconsistent approach can describe discs that are subject to a magnetised wind in a way similar to the α disc model. In this article, I explore the parameter space of winddriven protoplanetary discs in a systematic manner and present scaling laws that can be used in reduced models in a similar way to α disc models.
Methods. I computed a series of selfsimilar wind solutions, assuming that the disc is dominated by ambipolar and Ohmic diffusion. These solution were obtained by searching for stationary solutions in the finitevolume code PLUTO using a relaxation method and continuation.
Results. Selfsimilar solutions are obtained for values of plasma β ranging from 10^{2} to 10^{8} for several Ohmic and ambipolar diffusion strengths. Mass accretion rates of about 10^{−8} M_{⊙} yr^{−1} are obtained for the poloidal field strength β = O(10^{4}) or equivalently, 1 mG at 10 AU. In addition, the ejection efficiency is always close to 1, implying that wind massloss rate can be higher than the inner mass accretion rate when the windemitting region is large. The resulting magnetic lever arms are typically lower than 2, possibly reaching 1.5 in the weakest field cases. Remarkably, the mean transport properties (accretion rate and massloss rate) mostly depend on the field strength and much less on the disc diffusivities or surface density. The disc internal structure is nevertheless strongly affected by Ohmic resistivity, strongly resistive discs being subject to accretion at the surface while ambipolar only models lead to midplane accretion. Finally, I provide a complete set of scaling laws and semianalytical wind solutions, which can be used to fit and interpret observations.
Conclusions. Magnetised winds are unavoidable in protoplanetary discs as soon as they are embedded in an ambient poloidal magnetic field. Very detailed disc microphysics are not always needed to describe them, and simplified models such as selfsimilar solutions can capture most of the physics seen in full 3D simulations. The remaining difficulty to set up a complete theory of winddriven accretion lies in the transport of the largescale field, which remains poorly constrained and is not well understood.
Key words: magnetohydrodynamics (MHD) / protoplanetary disks
The selfsimilar solutions presented in this paper are available for download on github: https://github.com/glesur/PPDwind
© G. R. J. Lesur 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Protoplanetary discs are relatively cold and dense objects, which typically last for a few million years around young stellar objects (YSOs). However, it is well known that these discs are primarily ‘accretion discs’, in which matter is slowly falling onto the central star. This accretion rate has now been measured in dozens of objects, and it typically lies in the range of 10^{−10}–10^{−7} M_{⊙} yr^{−1}. Accretion in astrophysical discs is usually thought to be due to a magnetohydrodynamic (MHD) instability, the magnetorotational instability (MRI, Balbus & Hawley 1991), which transports angular momentum outwards and mass inwards. It was quickly realised that applying the MRI in the context of protoplanetary discs is a difficult task because of the dramatically low ionisation fraction expected in these objects. The inclusion of Ohmic (Gammie 1996) and ambipolar diffusions (PerezBecker & Chiang 2011) led to the conclusion that MRI cannot account for the observed accretion rates in these discs by two to three orders of magnitude.
The lack of a proper mechanism to trigger turbulence and angular momentum transport in protoplanetary discs revived an old idea that was mostly abandoned since the finding of the MRI: magnetised winds. While the disc bulk might be too diffusive to sustain MHD turbulence, it is still weakly coupled to the ambient largescale field originating from the parent cloud of the YSO. This ambient field can in principle be enough to trigger a magnetic breaking of the disc, leading to mass accretion and the formation of a wind. This idea was originally proposed in the context of YSOs by Wardle & Koenigl (1993). It was later revived through numerical simulations, first in local shearing box setups (Bai & Stone 2013; Simon et al. 2013; Lesur et al. 2014) and then in global simulations (Gressel et al. 2015, 2020; Béthune et al. 2017; Bai 2017; Wang et al. 2019). All of these simulations include Ohmic and ambipolar diffusion, while some of them also include the Hall effect, which is thought to be important in the densest regions of the disc (around 1 AU). The strength of these diffusions is in turn computed from ionisation and chemical models of various complexity. In addition, almost all of these models include some sort of heating of the disc atmosphere, designed to mimic nonthermal radiation heating that is expected in these systems (e.g. Thi et al. 2019). Each published model is of very high complexity, and only a handful of simulations can be performed to explore a (very) limited subspace of the parameter space.
In this context, the purpose of this work is to take a step back and simplify the physics in order to explore the coupling between protoplanetary discs and magnetised winds in a more or less systematic manner. This work essentially follows the approach initiated by Ferreira (1997) on globally selfsimilar models that was later extended to viscous (Casse & Ferreira 2000a) and weakly magnetised (JacqueminIde et al. 2019) discs. I use a relatively different technic to derive the wind solutions, however, which I present in the next section, along with the physical model. I then focus on a fiducial series of wind solutions and explore their physical properties. Finally, I vary the disc diffusivities to offer an overview of the more exotic configurations before concluding.
2 Model
In the following, I consider a magnetised disc orbiting a central object of mass M. The disc is only subject to the gravitational pull of the central object (selfgravity being neglected) and to the Lorentz force due toelectrical currents. The disc is assumed to be weakly ionised, hence I consider a generalised Ohm law that includes Ohmic and ambipolar diffusivities. These diffusivities are prescribed in Sect. 2.3.
2.1 Equations
In the following, I use either spherical (r, θ, φ) or cylindrical (R, φ, z) coordinates, depending on the context. The cylindrical radius is written in upper case, and the spherical radius is lower case. I solve the nonideal MHD equations, here in spherical coordinates, (1) (2) (3) (4)
where I havedefined the isothermal sound speed c_{s}(r), assuming the flow was locally isothermal, the plasma current J = c∇ × B∕4π, the gravitational constant G the magnetic field direction , the electromotive field and the Ohmic, Hall, and ambipolar diffusivities η_{O}, η_{H}, and η_{A}. The diffusivities are also functions of space and magnetic field strength. I stress that ‘no turbulence is assumed’ in this model, only molecular magnetic diffusivities that result from the low ionisation fraction of the plasma. If accretion occurs, it is only the result of torques that are selfconsistently computed in the model.
2.2 Selfsimilar ansatz
In the context of protoplanetary discs, the strong diffusivities are known to almost suppress all nonaxisymmetric structures (Béthune et al. 2017; Bai 2017) in the regions outside 1 AU. For this reason, I assume that the flow is 2.5D: I conserve three components for v and B but neglect φ derivatives.
In order to simplify the problem even more, I assume that the flow is globally selfsimilar. This implies that with increasing distance from the central object, the flow “looks” the same. This approach has several advantages: first it avoids issues with the inner boundary conditions, which are usually problematic in global numerical models. Second, it allows us to systematically explore the parameter space at exquisite resolution with limited numerical cost because the problem essentially becomes a 1D problem. I follow Ferreira & Pelletier (1993) and define the selfsimilar scaling for any field Q as (5)
where γ_{Q} is the selfsimilar exponent, and is a 1D function that completely determines the flow. The scaling of the gravitational force as 1∕r^{2} imposes the selfsimilar scaling to all of the other components, that is,
Similarly to the velocity field, the sound speed c_{s} is proportional to r^{−1∕2}, and is a prescribed function of θ. The selfsimilar scaling is here written in spherical coordinates, but can be transformed into cylindrical (R, z, φ) coordinates by noting that r = Rsin(θ) and z∕R = tan^{−1}(θ),
where .
In the following, it will be useful to define the Keplerian velocity and the Keplerian angular frequency Ω_{K} ≡ v_{K}∕R. I also use the disc geometrical thickness h ≡ c_{s}(R, z = 0)∕Ω_{K}(R). By construction from the selfsimilar scaling, h∕R ≡ ε is a constant of the problem.
2.3 Diffusivities
The diffusivities η are usually computed from complex thermochemical networks (e.g. Thi et al. 2019). It is customary to use dimensionless numbers to quantify the diffusivities: the Ohmic Reynolds number Rm and the ambipolar Elsasser number Λ_{A}, defined as (6) (7)
these dimensionless numbers are more useful than the more traditionally used Elsasser numbers because in most cases, they do not depend on the magnetic field strength (Wardle & Ng 1999). Hence, they only depend on the gas properties (density, temperature, and composition).
The use of the selfsimilar ansatz implies that Rm and Λ_{A} must be functions of z∕R (or θ) only. This is not true in full chemical models, except for ambipolar diffusion, which is known to be of the order of unity across a wide range of scales (see e.g. Thi et al. 2019, Fig. 8, middle panel). In the following, I therefore mostly focus on models dominated by ambipolar diffusion, and introduce Ohmic diffusion subsequently.
Following the Thi et al. (2019) models, I prescribe the ambipolar diffusion profile to be (8)
where λ is a free parameter quantifying the thickness of the nonideal layer in units of h.
In models including Ohmic diffusion, I use the following diffusivity profile:
This profile,and in particular the dependence on the density ratio, is chosen to be consistent with the ambipolar diffusion profile, assuming the plasma is made of two types of charged species (such as electrons and molecular ions).
2.4 Numerical method
In contrast to commonly employed selfsimilar approaches (Casse & Ferreira 2000a; JacqueminIde et al. 2019) where stationary equations are solved with shooting methods through the critical points of the flow, here I solved the timedependent Eqs. (1)–(3) using the code PLUTO, a finitevolume, shockcapturing scheme (Mignone et al. 2007). I used a spherical domain, with the shape of a shell with only one grid point in the radial direction and 2048 points distributed homogeneously in the θ direction. I chose the grid to be centred on r = r_{0} = 1, and the shell extended from θ = 0.15 to θ = π − 0.15. The radial boundary conditions were set enforcing the selfsimilar relations described in Sect. 2.2. For the boundary in the θ direction, I used standard outflow boundary conditions.
The initial condition was a disc in hydrostatic equilibrium threaded by a largescale vertical (z) magnetic field whose initial value was set to have β = 10^{5} in the disc midplane. This initial condition is strongly unstable and the disc very quickly (in less than 10 Ω^{−1}) launched an outflow before reaching a quasi steadystate.
Given that the code is timedependent and my choice of boundary conditions, the disc mass is not necessarily conserved. A fraction of the accreted material can be lost in the wind, leading to a slow decrease of the disc mass. In previous models, this has been taken into account by slightly adjusting the power exponent γ_{ρ} to ensure that mass is constant in the disc. This is not possible in our case because the solution is dynamically evolving: γ_{ρ} would need to be changed as a function of time. Instead, I therefore chose to renormalise the mass at each time step to keep the total mass inside the domain constant.
For similar reasons, the total magnetic flux threading the disc can evolve on secular timescales in this approach. This is not allowed in commonly used selfsimilar approaches, which usually assume that the total toroidal electromotive field is null. Here, I chose not toenforce such a constraint. This allowed me to measure the transport of magnetic flux, but it also implies that the magnetic flux evolves with time. As for the density, I therefore multiplied the each component of the field by a fixed factor at each time step, adjusted so that the magnetic flux threading the disc was equal to the desired value.
2.5 Accretion theory and diagnostics
Several diagnostics can be derived from selfsimilar solution. The most useful diagnostics are related to accretion theory, that is, how the disc surface density and mass accretion rate evolve with space and time. I therefore developed the mass, angular momentum, and magnetic flux conservation equations as (9) (10) (11)
where the overline denotes an azimuthal and ensemble average (in the case of a selfsimilar stationary solution, this average is not strictly needed), is the disc surface density, z_{0} is the disc surface, is the mass accretion rate, T_{xϕ} = ρv_{x}(v_{φ} − v_{K}) − B_{x}B_{φ}∕4π is the stress φ component, and B_{z0} = B_{z}(z = 0) and are the vertical magnetic field strength and azimuthal EMF in the disc plane, respectively. It should be noted that in principle, accretion theory only require the first two relations. However, it is now well established that in a winddriven disc, stresses and massloss rates are also a function of the mean poloidal field strength. In this context, accretion theory therefore has to be supplemented by Eq. (11).
Although Eqs. (9)–(11) fully describe the secular evolution of any disc at hand, it should be realised that the right handside terms are a priori unknown. This closure problem is well known in the disc community, and is usually solved by using the α disc paradigm when only the T_{rφ} term is present. In the case of a winddriven disc, four terms are actually present. Following the α disc idea, I therefore define four parameters, (12) (13) (14) (15)
where I have assumed that the disc was isothermal to define α. With these definitions, a positive v_{B} implies that the field is ‘transported outwards’. The massloss rate ζ is comparable to the definition of Scepi et al. (2018), up to a factor of the order of unity. In numerical applications, I use z_{0} = 6h.
It is also customary to measure the magnetic field strength as a function of the plasma β parameter. In this manuscript, I define the plasma β on the strength of the mean vertical field threading the disc midplane, (16)
The transport coefficients ζ, α, υ, and v_{B} being dimensionless coefficients, they are expected to depend only on dimensionless numbers. We therefore expect them to depend on the magnetic field strength β, but also on the diffusivity of the disc Rm and Λ_{A}. In principle, when all of these dependencies are known, accretion theory is said to be complete and the evolution of a winddriven protoplanetary disc can be computed in a way similar to the historical α disc. Therefore, determining these dependencies is the main objective of this article.
Fig. 1 Spacetime diagram showing the evolution for the first 1000 Ω^{−1} of B_{θ} (top), B_{φ} (middle), and v_{θ} (bottom) for the fiducial run. The system quickly reaches steady state. In order to continue the solution as a function of the field strength, the field is slowly increased every 300 Ω^{−1}. 
3 Fiducial simulation
3.1 Time evolution
In the following, I define the fiducial run as a simulation with ambipolar diffusion only, Λ_{A} = 1 and λ = 3, so that the nonideal MHD zone extends up to 3h above the disc midplane. The disc aspect ratio was fixed to ε = 0.1 and the disc and atmosphere were assumed to be locally isothermal, that is, c_{s} (R, Z) = c_{s}(R, z = 0). As for all the simulations, the fiducial run was started with a midplane β =10^{5}. I have also run this simulation at halfresolution (i.e. 1024 points in θ) to confirm the convergence of the solution. The two solutions differ by less than 1%, so I am confident that the full resolution simulations presented here are perfectly resolved.
While most of the solutions presented in this article are quasisteady state, I briefly discuss the spacetime evolution of the simulations used to compute the solutions. This spacetime evolution is shown in Fig. 1 for the fiducial run. The system clearly rapidly relaxes into steady state. During the first 100 Ω^{−1}, I applied a linear damping to the equations of motion to damp epicyclic oscillations more rapidly than the system would naturally do. These oscillations are due to the sudden launching of the wind, and are a spurious result of my windfree initial conditions (we used a similar procedure in the shearing box Lesur et al. 2013).
After the first 100 Ω^{−1}, the system reaches a quasisteady state, from which I can compute the required information: angular momentum transport, mass flux, etc. In order to continue the solution as a function of β, I slowly increased the total field strength between 300 Ω^{−1} and 400 Ω^{−1}, and let the system reach its new equilibrium for the next 200 Ω^{−1}. I then repeated this procedure until the code stopped because of singular points in the solution. In doing so, I reached β =35 in the fiducial run in steps of 10% decrease in β. In the fiducial setup, I also continued the solutions from β = 10^{5} up to β = 10^{8} in steps of 10% increase to study the domain in which the solutions obtained by this procedure are valid. Hence, the domain explored for the fiducial setup is wider than for the other solutions (Sect. 4). In the following, I average the quasisteady flow obtained for the last 100 Ω^{−1} of “relaxation” for each magnetisation and use the result to compute the flow properties.
3.2 Flow topology
The 2D topology of the flow can be deduced from the 1D “shell” solutions of the simulation by reconstructing the 2D fields using the selfsimilar scalings (5). Using this procedure, I computed the flow topology for the two extreme cases (β = 10^{5} and β = 35) and for an intermediate case (β = 10^{3}) shown in Fig. 2. In this fiducial simulation, the flow was symmetric with respect to the disc midplane, as seen here. The strongly magnetised solution produces a supersonic accretion flow, while the weaker magnetised case tends to produce subsonic accretion (v_{p}∕c_{s} < 0.1). In the weakfield case, the field lines are almost straight, indicating that the decoupling between the field and the flow is efficient, thanks to ambipolar diffusion. As the field becomes stronger, however, the field lines are clearly pinched around the midplane. This indicates an efficient field dragging by the accretion flow.
The wind was super fastmagnetosonic in all cases. The fast surface was closer to the disc surface when the magnetisation became stronger, while the Alfvén surface tended to stay at the same altitude (around 4.5h). Generally speaking, stronger magnetisation leads to faster and more massive winds (for an identical disc). Interestingly, while the shape of the stream and field lines differ significantly in the disc, their aspect is very similar in the wind region starting above the disc surface. In all of these cases, the flow showed signs of collimation towards the axis.
3.3 Transport properties
The transport properties defined in Sect. 2.5 can be evaluated as a function of the magnetisation for the fiducial run. This gives the dependencies shown in Fig. 3. The resulting transport coefficients are remarkably close to power laws of β (except for v_{B}, which shows a shallower dependence on β than a power law). Several important trends can be deduced from this result.
First, it should be pointed out that because the solutions are all steady state and laminar, the α measured here is a ‘pure laminar stress’. It is not due to any form of turbulence because it is just the result of the largescale wind that transports angular momentum both radially and vertically. Nevertheless, it can be compared to what is found in turbulent discs. The dependence of α on β^{−1} is steeper than the one usually found in idealMHD shearingbox simulations, which usually give α ∝β^{−1∕2} (e.g. Salvesen et al. 2016; Scepi et al. 2018). However, the ideal and nonideal values for α match for β ~10. Hence, my values of α are generally lower than their ideal MHD counterpart by a ratio ~ β^{1∕2}∕3. This is mostlikely the signature of ambipolar diffusion in the disc.
In general, the values for α, υ, and ζ are compatible with the values found in the shearingbox literature for β ~10^{4}−10^{5} and equivalent diffusivity regimes dominated by Λ_{A} = 1 (e.g. Bai 2013; Lesur et al. 2014), even though the dependence on β is slightly steeper than Bai (2013).
Interestingly, v_{B} > 0 for all the cases explored in the fiducial run. This indicates that despite the strong accretion in the disc midplane, the field is always diffusing outwards. This is qualitatively consistent with Bai & Stone (2017), but the diffusion speeds are reduced by an order of magnitude in my case. For β = 10^{4}, I obtain v_{B} = 1.7 × 10^{−3}, while Bai & Stone (2017) obtained^{1} v_{B} = 4 × 10^{−2}. A careful examination of the ideal and nonideal contributions to (Fig. 4) shows that advection almost completely balances diffusion, diffusion winning only by a 1% excess. It is also worth noting that while the ideal and ambipolar EMFs vary strongly in the disc, their sum is approximately constant with θ, indicating that poloidal field lines are “moving” radially without much vertical deformation in the disc.
Using Eq. (10), the safesimilar scalings (Eq. (5)) and the definitions of Sect. 2.5, it is straightforward to show that the mass accretion rate is directly related to α and υ_{±}
where I have split Ṁ_{acc} in a contributionfrom the radial (ṁ_{R}) and vertical (ṁ_{z}) stresses. Even though these are usually referred to as “turbulent” and “wind” contributions (even when the flow is laminar), I emphasise that there is no turbulent stress. Hence both contributions are intrinsically due to the wind. Using the powerlaw scalings in Fig. 3, I find that
Hence whenever β ≫ 1, the vertical stress is largely dominant and the radial transport of angular momentum can be neglected altogether, which is consistent with previous 3D numerical simulations (Béthune et al. 2017).
Neglecting ṁ_{R}, I can derive the mass accretion rate in a realistic disc using only the scaling for υ_{±}, which gives (17)
or substituting β with the field strength and the surface density, (18)
It is worth noting that the mass accretion rate does not depend strongly on the disc surface density, in contrast to the usual viscous disc model. Instead, I find that the field strength is the dominant control parameter in these solutions. More quantitatively, I find accretion rates compatible with expected rates from observations (~ 10^{−8}–10^{−7} M_{⊙} yr^{−1}) when β ∈ [10^{3}, 10^{5}] or B_{z} ~ a few mG. I caution that the observed accretion rates are typically inferred from the material that falls onto the central star, while here the mass accretion rate is computed in the disc at 10 AU. If the wind extracts a significant amount of material between the star and 10 AU, these accretion rates have no reason to match.
A common way to quantify the fraction of the accreted mass that is ejected is through the ejection efficiency parameter ξ, which is defined as
where the last equality results from the continuity equation assuming ∂_{t} Σ = 0. It is easy to show that the ejection efficiency is directly connected to the transport coefficients:
where the first equality assumes that ṁ_{R} is negligible,while the second uses the scaling found in the fiducial run (Fig. 3). This implies that ξ = O(1) for β ~10^{4} and that the accretion rate is approximately a linearly increasing function of radius in steady state in these models.
Finally, while most of the scalings presented here use β and therefore the vertical field strength as the main control parameter, these solutions are dominated by the azimuthal field component B_{φ}. Using the definition of υ, it easy to show that (19)
where the last equality comes from the scaling of the fiducial run. Hence, for β = 10^{4}, B_{φ}  is more than 14 times larger than B_{z} at the disc surface, indicating that the field is strongly wrapped at this location. Finally, because υ_{+} = −υ_{−}, the sign of B_{ϕ} is reversed at the top and bottom of the disc.
Fig. 2 Flow topology in the fiducial run (ambipolar diffusion only Λ_{A} = 1). Top row: poloidal streamlines (white) and log of the sonic Mach number. Bottom row: poloidal field and log of density, normalised so that the midplane density at R = 1 is unity. The colour scales are identical in the columns. From left to right, the disc magnetisation increases: β = 10^{5}, 10^{3}, and 35. The green lines denotes critical lines of the flow: Alfvénic (plain) and fast magnetosonic (dotdashed). The dashed green line represents the disc “surface” where the flow becomes ideal, which is arbitrarily located at z = 3.5h for all of the solutions. 
Fig. 3 Transport coefficient as a function of the disc magnetisation for the fiducial run (ambipolar diffusion only ). Each blue dot is a steadystate solution obtained by continuation in the simulation. Because the flow is symmetric with respect to the midplane, only values computed from the disc top surface (ζ_{+} and υ_{+}) are shown. The best fits using power laws are shown as a plain orange line. 
Fig. 4 Zoom on the toroidal electromotive field of the disc, splitting the ideal MHD and ambipolar diffusion contributions for β = 10^{5} in the fiducial case. The two electromotive fields balance almost exactly. The total EMF at θ = π∕2 is and is approximately constant throughout the disc. This gives v_{B} = 6 × 10^{−4}. 
3.4 Magnetic torque and currents
Although it is possible to interpret accretion in terms of angular momentum budget (the disc angular momentum is transported into the wind), this gives us little information about how the disc accretes. This question can be addressed by considering the forces acting on the disc. In particular, because the outflow is magnetised, the Lorentz force F_{L} = J × B∕c is key. If F_{L,φ} < 0, the Lorentz force breaks down the disc rotation by creating a negative torque, driving accretion. The magnetic torque is clearly closely related to the poloidal current,
Because the poloidal field is more homogeneously distributed in the disc (Fig. 2), the magnetic torque is strongest in the regions of strong J_{p}. Hence the poloidal current is a direct indicator of the torque exerted by the wind on the disc. The distribution of poloidal currents in the fiducial β = 10^{5} case is shown in Fig. 5. I find a strongly focussed current sheet in the disc midplane, with currents directed outwards. This current induces a strong negative torque that in turn is responsible for accretion in the disc midplane. This current thencloses back in the disc surface and in the wind.
It might be surprising at first sight to find such a narrow current sheet in the system given that the disc is strongly diffusive (Λ_{A} = 1). It is well known, however, that ambipolar diffusion is a selffocusing diffusion process. It has been argued that in “free” reconnection sites, the current scales as z^{−2∕3}, where z is the distance from the null point (Zweibel & Brandenburg 1997). In the particular case presented here, I find that the current is even more focused with J_{r} ∝ z^{−1} close to the midplane. This is expected because in the present case, I focus on a reconnection site that is “forced” by the outflow launched from the surface, hence the scaling need not be identical to that of Zweibel & Brandenburg (1997).
Fig. 5 Poloidal current in the fiducial simulation at β = 10^{5}. The current intensity is represented in colour. There is a strong outwards current in the disc midplane. The green lines represent the same characteristics as Fig. 2. 
Fig. 6 Zoom on the radial current in the fiducial simulation at β = 10^{5}. The current sheet is strongly focused, with a scaling compatible with J_{r} ∝ z^{−1} (see text). 
3.5 Wind properties
Winds in protoplanetary discs are often refereed to as magnetothermal winds. This is because most (if not all) of the global simulations published so far used a “hot atmosphere”. These hot atmospheres are either achieved by a prescribed heating function (Béthune et al. 2017), or with a more complete treatment of the thermodynamics, including various sources of heating and cooling (Wang et al. 2019; Gressel et al. 2020). In any case, the energetics of these winds is partially dominated by thermal heating, hence their name.
In the fiducial solution presented here, the equation of state is locally isothermal. The temperature of the atmosphere therefore is that of the disc at the same cylindrical radius R. For this reason, my fiducial wind is not magnetothermal, as I show below.
Because the wind is a steadystate ideal MHD flow, several conserved quantities can be derived from the equations of motion. I therefore considered a streamline passing through the midplane at R_{0}. I define Ω_{0} as the Keplerian angular velocity at R_{0}: Ω_{0} ≡ Ω_{K}(R_{0}). I then follow Blandford & Payne (1982), defining

the massloading parameter
where v_{p} and B_{p} are the poloidal velocity and field strength, which quantifies the mass that is loaded onto the field lines.

the rotation parameter
which measures the rotation speed of magnetic surfaces.

the magnetic lever arm
which measures the angular momentum that is extracted by the wind.

the Bernoulli invariant
which measures the energy content of the flow. In this definition, I included the thermal contribution to the flow energetics as the work done by thermal pressure forces along the field line,
If the flow is adiabatic, this is simply the enthalpy. However, the locally isothermal approximation does not lead to an adiabatic flow, hence I have to use the integral form.
I computed the mass loading, rotation, and lever arm parameters along a magnetic field line for the fiducial simulation at β = 10^{5} (Fig. 7). As expected, these parameters are constant once above the nonideal region of the disc. Surprisingly, the lever arm is very small, reaching only λ = 1.56. As is well known, cold MHD winds require λ > 3∕2, so this is just above this limit value. A close inspection of the Bernoulli invariant (Fig. 8) shows that thermal driving is indeed negligible in the energy budget of the outflow, confirming that the solution is a cold MHD wind, not a magnetothermal wind. It also shows that indicating that the wind is free to escape the gravitational potential well up to infinity.
A systematic exploration of the κ, λ parameters as a function of β in the fiducial simulation gives Fig. 9. This figure can directly be compared to Fig. 2 in Blandford & Payne (1982). It is interesting to note that these new solutions at low magnetisation have a much larger κ and much smaller λ than the previous solutions of Blandford & Payne (1982). This trend has been observed by JacqueminIde et al. (2019) in ideal MHD for β >1. The solutions I obtained here have even smaller λ, with an asymptote in the high β limit close to λ = 1.4. This is smaller than the λ = 3∕2 limit demonstrated by Blandford & Payne (1982). It should be pointed out, however, that this limit is computed for a wind emitted from z = 0, while here the wind is emitted from z = 5H, so the lever arm is physically allowed to be slightly smaller than 3∕2. This is confirmed by analysingthe Bernoulli invariant at β = 10^{7}, which still shows that the magnetic energy dominates the energetics at the disc surface (not shown).
Fig. 7 MHD invariants plotted along one given field line for the fiducial run at β = 10^{5}. As expected, the invariants are approximately constant once in the idealMHD region lies above the disc. The asymptotic values are κ = 30.2, ω = 0.99, and λ = 1.56. 
Fig. 8 Bernoulli invariant plotted along one given field line for the fiducial run at β = 10^{5}. The Bernoulli invariant becomes constant once in the ideal MHD region, as expected. The thermal contribution to the energetics isnegligible compared to the magnetic contribution, hence the wind is classified as a “cold” MHD wind. 
Fig. 9 Dependence on κ and λ of the fiducial solutions as a function of the disc magnetisation β. As the field strength increases, λ increases, while the massloading parameter decreases. 
4 Parameter space exploration
4.1 Ambipolar diffusion
In order to explore the effect of the disc diffusivity on the behaviour of the outflow, I varied the strength of ambipolar diffusion in the disc plane (Eq. (8)), keeping the same vertical profile. The resulting transport properties are shown in Fig. 10. The general trend is that as diffusion is increased, the massloss rate and angular momentum fluxes decrease, while magnetic flux transport increases. This overall suggests that winddriven accretion tends to be reduced when ambipolar diffusion increases, as expected naively. It should be noted, however, that not all transport coefficients vary in the same proportion. Taking β =10^{4} as the reference case, I find that a decrease in by a factor 16 leads to adecrease in ζ, α, and υ by a factor 3.2, 4.0, and 2.25, respectively, while v_{ψ} increases by a factor 4.0. This therefore suggests a relatively weak sensitivity of the transport coefficient on (shallower than ), especially for the vertical angular momentum coefficient υ. The weak dependence of υ on the disc dissipation properties has been observed in shearingbox simulations (see Lesur 2020, Fig. 39 and related text) and is here confirmed in global geometry. Overall, this confirms that winddriven angular momentum extraction is only moderately affected by the strength of ambipolar diffusion, and that the mass accretion rate (17) is approximately valid (up to a factor of a few) for the range of ambipolar diffusion considered here.
4.2 Ohmic diffusion
In addition to the strength of ambipolar diffusion, I also explored the effect of Ohmic diffusion on the flow topology. In protoplanetary discs, this leads to wind solutions that are valid in the Gammie (1996) historical “dead zone”, found around 1 AU in typical chemical models. To limit the computation time required by such models, I chose to focus on Rm_{0} > 1. In the literature, Ohmic diffusion is frequently evaluated in terms of the Ohmic Elsasser number Λ_{O}, which depends on the field strength. Typically, Λ_{0} = Rm_{0}β, so that Rm_{0} ≳ 1 typically corresponds to the midplane value quotedat R ≲ 1 AU by several authors (Bai & Stone 2013, Fig. 1; Thi et al. 2019, Fig. 8). I finally point out that in all of these models, I kept ambipolar diffusion fixed to its fiducial value . This means that Ohmic diffusion is added to the diffusivity tensor of the fiducial run.
As for ambipolar diffusion, I started with the transport coefficients as a function of Rm_{0} (Fig. 11). Here, the angular momentum transport coefficients α and υ are barely affected by Ohmic diffusion, even when Rm_{0} = 1. This suggests that the accretion rate in the Gammie (1996) dead zone probably is not so different from that observed in regions dominated by ambipolar diffusion. I note, however, that the field diffusion rate is increased by possibly several orders of magnitude. This is probably the most stringent effect of adding Ohmic diffusion to the system. While Rm_{0} ≥ 10 solutions are topdown symmetric (therefore only ζ_{+} and υ_{+} are plotted), the Rm_{0} = 1 solution breaks this topdown symmetry around β ≃ 10^{4}. The direction of this symmetry breaking is random. Starting from slightly different initial condition, I could obtain ζ_{+} > ζ_{−} or the opposite, which shows that there is ‘no physically preferred direction’ for this symmetry breaking.
The fact that the transport coefficients are approximately unaffected by Ohmic diffusion should not lead to the incorrect impression that the solutions are physically identical. A careful examination of the solution at Rm_{0} = 1 (Fig. 12) shows that the solution in the disc is very different from the fiducial solution. First, the current sheet that was initially localised in the midplane (Fig. 5) is now localised on the disc surfaces. This can be easily interpreted as a result of Ohmic diffusion. While ambipolar diffusion can have a focusing effect, Ohmic diffusion is really only a diffusive process. The current sheet therefore cannot form in the strongly diffusive midplane, and is therefore pushed away, at the disc surface.This change in localisation of the current sheet in turn affects the accretion flow because accretion is driven by the magnetic torque created by the poloidal current (see Sect. 3.4). This is more evidently seen in Fig. 13, where the accretion flow clearly moves from the disc midplane towards the disc surface as Rm_{0} is decreased, thereby following the localisation of the current layer. Finally, poloidal current and accretion flows are both localised at the disc surface. However, because the current is set globally by the electrical circuit enforced by the wind, its intensity is not really affected by this change in localisation, implying that total torques and accretion rates are not dramatically modified, as indicated bytransport coefficients.
After understanding that the current sheet becomes divided into two parts and pushed towards the surface, this process might be imagined to be not be totally symmetric for some range of parameters. This is what happens in the Rm_{0} = 1 solution around β = 10^{4}. In this case, the current sheet is more pronounced in the southern hemisphere (Fig. 14), leading to an accretion stream that is morepronounced on the bottom side of the disc, and inclined field lines in the disc bulk. I stress that while accretion occurs in the southern side of the disc in this particular example, the wind is more pronounced on the northern side of the disc because ζ_{+} > ζ_{−} and υ_{+} > υ_{−}. Hence angular momentum and mass flow away on the side opposite to that of accretion. This behaviour was also observed by Béthune et al. (2017) in dissymmetric solutions.
5 Discussion and conclusions
I have used a finitevolume code to derive selfsimilar wind solutions applicable to the diffusive regime of protoplanetary discs. In contrast to previous selfsimilar solutions, no assumption was made on flux stationarity or topdown symmetry. I also circumvented the difficulties of fullblown 3D numerical simulations, which are often limited by their integration time and inner boundary conditions. Using this technic, it is possible to explore systematically a wide range of parameters at a reduced numerical cost.
Using this tool, I have presented a series of wind solutions ranging from β = 10^{8} to β = 35 that are valid in the nonideal regions of protoplanetary discs R ≳ 1 AU. I showed that cold wind solutions (i.e. that do not require atmospheric heating) exist in this entire range of parameters. Some of these solutions are very similar to the solution found using 2.5D and 3D simulations, in the same range of parameters. These solutions exhibit several important properties, some of which can be confirmed observationally.
First and foremost, magnetised disc winds always exist and are unavoidable as soon as a largescale magnetic field threads the disc. This statement is true even for β = 10^{8} fields that correspond to a few μG at 10 AU for typical surface densities. They give accretion rates Ṁ_{acc} (Eq. (18)), which essentially depend on the field strength, and which are broadly speaking compatible with observed accretion rates onto T Tauri stars provided that 10^{3} < β < 10^{5} or equivalently B_{z} ~ a few mG. Interestingly, Ṁ_{acc} depends only weakly on the disc surface density, in contrast to viscous disc models. However, the direct comparison between the accretion rate in the disc bulk (as measured by Eq. (18)) and the mass accretion measured onto the star is probably misleading. The latter is probably significantly smaller than the former because the ejection index (see below) is about one, hence the disc mass accretion rates could be significantly higher (possibly by an order of magnitude) than the rates quoted here.
The ejection index ξ, which essentially quantifies the ratio of ejected to accreted mass, depends only weakly on the magnetic field strength and is of the order of (but slightly lower than) unity for the range of magnetisation I explored. This means that the mass loss rate is comparable to the mass accretion rate at a given radius. Consequently, when steady state is assumed, the accretion rate is expected to increase significantly with radius. As expected from a high ξ wind (Lesur 2020, Eq. (11.22)), the wind lever arm is always relatively small, with typically λ < 2 for β >10^{4}, without requiring any heating. Such a high ξ and low λ has been found in previous numerical work (Béthune et al. 2017; Bai 2017), and it was initially thought that atmospheric heating, included in all these models, was the main reason for this result, following the argument of Casse & Ferreira (2000b). I showed here that heating is not key to obtain these high ξlow λ solutions, a result that was also recently reported by JacqueminIde et al. (2019) in the context of fully ionised discs. The solutions presented here are locally isothermal, and I showed that thermal heating is negligible in the wind energy budget even for β =10^{8}. This means that low λ and high ξ are quite clearly signatures of weakly magnetised (β ≫ 1) outflows.
The vertical field strength predicted to derive accretion rates compatible to typical T Tauri stars (about a mG at 10 AU, see Eq. (18)) are compatible with recent upper limits obtained from Zeeman measurements (e.g. Vlemmings et al. 2019). For field detection, it is quite clear that the toroidal is a better candidate than the vertical component because the former is expected to be much stronger (by a factor 10 or so at β ~10^{4}) than the latter (Eq. (19)). The toroidal component is expected to change sign across the midplane, however, so that this detection is possible only for tracers that do not average out the field through the disc thickness.
For the disc microphysics, I find that the discaveraged mass accretion and ejection rates depend only weakly on ambipolar and Ohmic diffusivities. This implies that very elaborated ionisation models in the disc are not required to obtain the right accretion or ejection bulk properties. However, the disc vertical structure does strongly depend on diffusion. Solutions with large Ohmic diffusion (i.e. valid closer to 1 AU) tend to exhibit accretion at the disc surface, while weaker Ohmic diffusion (R > 10 AU) shows accretion in the disc midplane. Some of the strong Ohmic diffusion solutions also exhibit topdown dissymmetry, where accretion mostly occurs on one side of the disc. An asymmetry like this was also found in simulations including Ohmic and ambipolar diffusion (Béthune et al. 2017; Gressel et al. 2020) in a similar range of parameters, so that this result has been reported before. These topology properties will likely have a strong effect on dust growth and planet migration theory, but again, they barely affect the vertically averaged transport properties.
Finally, in all of these models, the largescale magnetic field is found to be transported outwards, in agreement with Bai & Stone (2017). However, the transport velocity is about an order of magnitude lower than that measured by Bai & Stone (2017) in similar conditions (ambipolar diffusion only). Moreover, I find that the transport rate varies significantly with the strength of Ohmic diffusion and to a lesser extent with ambipolar diffusion. Overall, this dependence might explain why simulations and models tend to disagree on the magnetic transport rate (Bai & Stone 2017; Leung & Ogilvie 2019; Gressel et al. 2020).
Clearly, the largescale field transport is the main bottleneck of a complete theory of winddriven accretion because the strength of the largescale poloidal field is the main control parameter. This cannot be ignored for the secular evolution of discs because the field strength is expected to vary with time. While class II objects seem to require field strengths of a few mG at 10 s of AU to obtain the correct accretion rates (see above), corecollapse calculations tend to suggest a field strength of 100 mG when the disc forms (e.g. Masson et al. 2016). This means that the field strength must be reduced by at least one order of magnitude (more probably two) between class 0 and class 2, indicating that flux transport must be relativelyefficient. The problem of field transport therefore deserves more attention in the future if winds are to be used in secular evolution models.
Fig. 10 Transport coefficients as a function of the disc magnetisation and midplane Elsasser number . 
Fig. 11 Transport coefficients as a function of the disc magnetisation and midplane Reynolds number Rm_{0}. These solutions also include ambipolar diffusion with . The fiducial run is physically equivalent to Rm_{0} = ∞. 
Fig. 12 Flow structure for the solution Rm_{0} = 1 at β = 10^{5}. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. Compared to the fiducial case, the current layer is now localised at the disc surface, leading to a surface accretion flow. 
Fig. 13 Inward mass flux (top) and radial velocity Mach number (bottom) as a function of Ohmic diffusion strength at β = 10^{5}. The accretion flow moves towards the disc surface as Ohmic diffusion is increased, following the poloidal current sheet (see text). 
Fig. 14 Flow structure for the solution Rm_{0} = 1 at β = 10^{4}. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. This solution exemplifies the dissymmetric solutions found around β = 10^{4} at Rm_{0} = 1. 
Acknowledgements
I am thankful to Jonatan JacqueminIde, Etienne Martel, Jonathan Ferreira and Cornelis Dullemond for fruitful discussions regarding the physics of selfsimilar solutions. I acknowledge support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme (Grant agreement No. 815559 (MHDiscs)). All of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univgrenoblealpes.fr), which is supported by Grenoble research communities.
References
 Bai, X.N. 2013, ApJ, 772, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N. 2017, ApJ, 845, 75 [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 769, 76 [Google Scholar]
 Bai, X.N., & Stone, J. M. 2017, ApJ, 836, 46 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Casse, F., & Ferreira, J. 2000a, A&A, 353, 1115 [NASA ADS] [Google Scholar]
 Casse, F., & Ferreira, J. 2000b, A&A, 361, 1178 [NASA ADS] [Google Scholar]
 Ferreira, J. 1997, A&A, 319, 340 [Google Scholar]
 Ferreira, J., & Pelletier, G. 1993, A&A, 276, 625 [NASA ADS] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [CrossRef] [Google Scholar]
 JacqueminIde, J., Ferreira, J., & Lesur, G. 2019, MNRAS, 490, 3112 [CrossRef] [Google Scholar]
 Lesur, G. 2020, ArXiv eprints [arXiv:2007.15967] [Google Scholar]
 Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leung, P. K. C., & Ogilvie, G. I. 2019, MNRAS, 487, 5155 [Google Scholar]
 Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 PerezBecker, D., & Chiang, E. 2011, ApJ, 727, 2 [Google Scholar]
 Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857 [NASA ADS] [CrossRef] [Google Scholar]
 Scepi, N., Lesur, G., Dubus, G., & Flock, M. 2018, A&A, 620, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, J. B., Bai, X.N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013, ApJ, 775, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Thi, W. F., Lesur, G., Woitke, P., et al. 2019, A&A, 632, A44 [CrossRef] [EDP Sciences] [Google Scholar]
 Vlemmings, W. H. T., Lankhaar, B., Cazzoletti, P., et al. 2019, A&A, 624, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, L., Bai, X.N., & Goodman, J. 2019, ApJ, 874, 90 [Google Scholar]
 Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218 [Google Scholar]
 Wardle, M., & Ng, C. 1999, MNRAS, 303, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Zweibel, E. G., & Brandenburg, A. 1997, ApJ, 478, 563 [Google Scholar]
I measured v_{B} in units of c_{s}, while Bai & Stone (2017) quoted v_{B} in units of v_{K}, hence a factor of ε should be added when the two results are compared.
All Figures
Fig. 1 Spacetime diagram showing the evolution for the first 1000 Ω^{−1} of B_{θ} (top), B_{φ} (middle), and v_{θ} (bottom) for the fiducial run. The system quickly reaches steady state. In order to continue the solution as a function of the field strength, the field is slowly increased every 300 Ω^{−1}. 

In the text 
Fig. 2 Flow topology in the fiducial run (ambipolar diffusion only Λ_{A} = 1). Top row: poloidal streamlines (white) and log of the sonic Mach number. Bottom row: poloidal field and log of density, normalised so that the midplane density at R = 1 is unity. The colour scales are identical in the columns. From left to right, the disc magnetisation increases: β = 10^{5}, 10^{3}, and 35. The green lines denotes critical lines of the flow: Alfvénic (plain) and fast magnetosonic (dotdashed). The dashed green line represents the disc “surface” where the flow becomes ideal, which is arbitrarily located at z = 3.5h for all of the solutions. 

In the text 
Fig. 3 Transport coefficient as a function of the disc magnetisation for the fiducial run (ambipolar diffusion only ). Each blue dot is a steadystate solution obtained by continuation in the simulation. Because the flow is symmetric with respect to the midplane, only values computed from the disc top surface (ζ_{+} and υ_{+}) are shown. The best fits using power laws are shown as a plain orange line. 

In the text 
Fig. 4 Zoom on the toroidal electromotive field of the disc, splitting the ideal MHD and ambipolar diffusion contributions for β = 10^{5} in the fiducial case. The two electromotive fields balance almost exactly. The total EMF at θ = π∕2 is and is approximately constant throughout the disc. This gives v_{B} = 6 × 10^{−4}. 

In the text 
Fig. 5 Poloidal current in the fiducial simulation at β = 10^{5}. The current intensity is represented in colour. There is a strong outwards current in the disc midplane. The green lines represent the same characteristics as Fig. 2. 

In the text 
Fig. 6 Zoom on the radial current in the fiducial simulation at β = 10^{5}. The current sheet is strongly focused, with a scaling compatible with J_{r} ∝ z^{−1} (see text). 

In the text 
Fig. 7 MHD invariants plotted along one given field line for the fiducial run at β = 10^{5}. As expected, the invariants are approximately constant once in the idealMHD region lies above the disc. The asymptotic values are κ = 30.2, ω = 0.99, and λ = 1.56. 

In the text 
Fig. 8 Bernoulli invariant plotted along one given field line for the fiducial run at β = 10^{5}. The Bernoulli invariant becomes constant once in the ideal MHD region, as expected. The thermal contribution to the energetics isnegligible compared to the magnetic contribution, hence the wind is classified as a “cold” MHD wind. 

In the text 
Fig. 9 Dependence on κ and λ of the fiducial solutions as a function of the disc magnetisation β. As the field strength increases, λ increases, while the massloading parameter decreases. 

In the text 
Fig. 10 Transport coefficients as a function of the disc magnetisation and midplane Elsasser number . 

In the text 
Fig. 11 Transport coefficients as a function of the disc magnetisation and midplane Reynolds number Rm_{0}. These solutions also include ambipolar diffusion with . The fiducial run is physically equivalent to Rm_{0} = ∞. 

In the text 
Fig. 12 Flow structure for the solution Rm_{0} = 1 at β = 10^{5}. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. Compared to the fiducial case, the current layer is now localised at the disc surface, leading to a surface accretion flow. 

In the text 
Fig. 13 Inward mass flux (top) and radial velocity Mach number (bottom) as a function of Ohmic diffusion strength at β = 10^{5}. The accretion flow moves towards the disc surface as Ohmic diffusion is increased, following the poloidal current sheet (see text). 

In the text 
Fig. 14 Flow structure for the solution Rm_{0} = 1 at β = 10^{4}. Left: streamlines and sonic Mach number. Middle: field lines and density map. Right: poloidal current lines and current density. This solution exemplifies the dissymmetric solutions found around β = 10^{4} at Rm_{0} = 1. 

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