Issue 
A&A
Volume 667, November 2022



Article Number  A17  
Number of page(s)  25  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202142946  
Published online  01 November 2022 
Magnetised winds in transition discs
I. 2.5D global simulations^{★}
Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
email: etienne.martel@univgrenoblealpes.fr
Received:
17
December
2021
Accepted:
3
May
2022
Context. Protoplanetary discs are cold, dense, and weakly ionised environments that witness planetary formation. Among these discs, transition discs (TDs) are characterised by a wide cavity (up to tens of au) in the dust and gas distribution. Despite this lack of material, a considerable fraction of TDs are still strongly accreting onto their central star, possibly indicating that a mechanism is driving fast accretion in TD cavities.
Aims. The presence of radially extended ‘dead zones’ in protoplanetary discs has recently revived interest in magnetised disc winds (MDWs), where accretion is driven by a large magnetic field extracting angular momentum from the disc. We propose that TDs could be subject to similar disc winds, and that these could naturally explain the fastaccreting and longlived cavities inferred in TDs.
Methods. We present the results of the first 2.5D global numerical simulations of TDs harbouring MDWs using the PLUTO code. We imposed a cavity in the gas distribution with various density contrasts, and considered a powerlaw distribution for the largescale magnetic field strength. We assume the disc is weakly ionised and is therefore subject to ambipolar diffusion, as expected in this range of densities and temperatures.
Results. We find that our simulated TDs always reach a steady state with an inner cavity and an outer ‘standard’ disc. These models also maintain an approximately constant accretion rate through the entire structure, reaching 10^{−7} M_{⊚} yr^{−1} for typical surface density values. The MDW launched from the cavity is more magnetised and has a significantly larger lever arm (up to 10) than the MDW launched from the outer disc. The material in the cavity is accreted at sonic velocities, and the cavity itself is rotating at 70% of the Keplerian velocity due to the efficient magnetic braking imposed by the MDW. Overall, our cavity matches the dynamical properties of an inner jet emitting disc (JED) and of magnetically arrested discs (MADs) in blackhole physics. Finally, we observe that the cavity is subject to recurring accretion bursts that may be driven by a magnetic RayleighTaylor instability of the cavity edge.
Conclusions. Some strongly accreting TDs could be the result of magnetised wind sculpting protoplanetary discs. Kinematic diagnostics of the disc or the wind (orbital velocity, wind speeds, accretion velocities) could disentangle classic photoevaporation from MDW models.
Key words: accretion / accretion disks / protoplanetary disks / magnetohydrodynamics (MHD) / methods: numerical
© É. Martel and G. Lesur 2022
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.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Transition discs (TDs) are protoplanetary discs exhibiting a deficit of nearinfrared emission, indicating a significant drop in the abundance of small dust grains in the regions inside a few tens of au (Espaillat et al. 2014). These objects are believed to be the intermediate stage between ‘full‘ primordial T Tauri discs and discless young stellar objects, hence their name. In this framework, TDs are the result of an insideout dispersal process, which is usually believed to be a combination of viscous accretion, dust growth (Dullemond & Dominik 2005), giant planets (Marsh & Mahoney 1992), and photoevaporation (Clarke et al. 2001; Alexander et al. 2014).
Despite their cavities, a large fraction of TDs are accreting onto their protostars. While Najita et al. (2007) quoted a median accretion rate reduced by one order of magnitude in Taurus compared to ‘primordial’ discs, more recent studies find even stronger accretion rates. Fang et al. (2013) showed that accreting TDs have a median accretion rate similar to normal optically thick discs. Manara et al. (2014) finds that TDs accrete similarly to classic T Tauri stars and that there is no correlation between the accretion rate and the cavity size. The fact that TDs are accreting systems should not give the impression that their cavity is depleted only in dust grains: TDs also exhibit cavities in the gas distribution (Zhang et al. 2014), with gas surface density increasing with radius (Carmona et al. 2014). Probing rotational emission of CO, van der Marel et al. (2015, 2016) find a drop in gas surface density by two to four orders of magnitude, while the drop in dust surface density goes up to six orders of magnitude. Similar results hold in rovibrational CO lines, probing the cavity further in, leading to a gas drop of two to four orders of magnitude in the inner (<3 au) regions (Carmona et al. 2017).
The picture that emerges is that of discs with a drop in gas surface density by several orders of magnitude, which are accreting similarly to (or slightly less than) primordial discs. There can be only two explanations for this phenomenon: either accretion is due to a ‘hidden’ mass reservoir localised close to the star, and what we observe is the transient accretion of this reservoir, or gas somehow manages to penetrate the cavity with a much larger velocity than the usual viscous accretion velocity. In that case, one typically needs an accretion velocity of the order of the sound speed to reconcile the accretion rate with the drop in surface density (Wang & Goodman 2017).
In the first category of models, we find scenarios involving photoevaporation combined with an inner dead zone (Morishima 2012; Gárate et al. 2021). This inner dead zone, typically extending between 1 and 10 au, sets the radius of the mass reservoir and therefore the cavity inner edge. While it predicts a fraction of TDs with accretion rates of Ṁ ~ 10^{−9} M_{⊙} yr^{−1}, it also predicts a large fraction of nonaccreting TDs, which is not observed (Gárate et al. 2021). In addition, these models rely on the Ohmic dead zone model of Gammie (1996), while it is now understood that dead zones are much more extended radially because of ambipolar diffusion (Simon et al. 2013), casting doubts on the applicability of inner dead zone models. The second type of model requires a mechanism to boost angularmomentum transport in the cavity. The most studied candidate for this is planetdisc interaction with planets (typically more than 3) embedded in the inner cavity. This scenario, however, finds gaps that are not necessarily sufficiently ‘clean’ (Zhu et al. 2011) and predicts that multiple giant planet systems in resonance are much more common than observed (Dong & Dawson 2016).
It should be noted that all of these scenarios make the explicit assumption of viscous accretion, the viscosity being due to some kind of smallscale turbulence, which could be of hydrodynamic (vertical shear instability or VSI, Nelson et al. 2013) or magnetic (magnetorotational instability or MRI, Balbus & Hawley 1991) origin. It is, however, becoming clear that accretion in the regions outside of 1 au is probably partially driven by magnetic winds (Bai & Stone 2013; Lesur et al. 2014; Béthune et al. 2017). While the accretion rate of viscous models is proportional to the gas surface density, the accretion rate of magnetohydrodynamic (MHD) winddriven models is mostly controlled by the strength of the largescale magnetic field, and much less by the surface density (for instance, Lesur 2021b proposes Ṁ ∝ Σ^{0.2} B^{1.6}). Hence, if one carves a cavity in a disc without significantly modifying its magnetic field distribution, one could in principle create a population of accreting TDs not so different from classical T Tauri discs in terms of accretion rates. This kind of scenario is found in secular evolution models that include a realistic dependence of the wind stress on the surface density (e.g. Suzuki et al. 2016, see their Σdependent wind torque models). Hence, MHD winds could, in principle, generate and sustain a fastaccreting TD cavity.
The idea of having a magnetic winddriven cavity was first proposed by Combet & Ferreira (2008). In this work, the cavity (named the jet emitting disc, or JED) is diluted, accreting at sonic velocities, and it sustains an accretion rate similar to that of the outer disc. The same angle of attack was more recently tackled by Wang & Goodman (2017), who showed that the magnetic diffusion properties of TD cavities were reminiscent of the magnetic wind solutions of Wardle & Koenigl (1993), indicating that all of the conditions required for efficient magnetic wind launching were met in TD cavities. While this picture is promising to explain accreting TDs, no dynamical model exists connecting an outer ‘standard’ disc to an inner cavity accreting thanks to magnetised winds.
In this work, we present the first selfconsistent (under the standard MHD assumptions) numerical models of accreting TDs based on the MHD wind scenario. The model we propose does not enforce accretion (e.g., with an α parameter that would be added by hand). Accretion and the disc equilibrium are natural consequences of the first principles of MHD, in the sense that their origins lie within the magnetic stresses arising from the initial vertical magnetic field. Our aim is to demonstrate that a fast accreting cavity can connect to a standard windemitting outer disc, subject to realistic magnetic diffusion, and that the resulting configuration can be longlived. Given the richness of the dynamics, we first concentrated on 2.5D models in this work, and we will discuss 3D models in a followup paper. The paper is divided as follows. We first introduce the model equations, physical quantities, and numerical setup. We then discuss an indepth investigation of a fiducial model, which possesses a cavity with a drop of four orders of magnitude in gas surface density. We finally explore alternative models, varying the cavity depth and size, and the diffusion coefficients before concluding. We stress that we focus here on a proof of concept that such a TD configuration is sufficiently stable to be observable, but we do not discuss ‘how’ a primordial disc could have ended in such a configuration. This will be the subject of future work.
2 Physical and numerical setups
2.1 Physical model
2.1.1 Governing nonideal MHD equations
In the following, we place ourselves in the nonrelativistic, nonideal MHD regime and consider a thin, locally isothermal disc to follow the evolution of the gas. The mass and momentum conservation equations and the induction equation respectively read (1) (2) (3)
where ρ, P, u, and B are, respectively, the density, the thermal pressure, the plasma velocity, and magnetic field. Φ_{*}, = –GM_{*}/r is the gravitational potential due to the central star of mass M_{*}, G being the gravitational constant. To close this system of equations, we assume the plasma follows a nonideal Ohm’s law including ambipolar diffusion: (4)
where is a unit vector parallel to B, J is the electric current, c is the speed of light, and η_{A} is the ambipolar diffusivity. No turbulence was added in this model whatsoever. In addition to these equations, the plasma follows the Maxwell equations: (5)
We place ourselves in a spherical coordinate system (r,θ,ϕ) centred on the star. For convenience, we also introduce the cylindrical coordinates R = r sin θ, ϑ = φ and z = r cos θ.
Since we work in a thin disc, the azimuthal angular velocity Ω is expected to be close to the Keplerian angular velocity Ω_{K}(r) = (GM_{*}/r^{3})^{1/2}. It is therefore useful to introduce a deviation from the Keplerian velocity u, defined as (7)
with . We note that the latitudinal dependence of is somewhat arbitrary and need not be a particular equilibrium state. Here, our choice of ensures that our reference Keplerian velocity has constant specific angular momentum on spherical shells and eliminates surface terms that are otherwise present in angularmomentum conservation equations (e.g. the last term of Eq. (16) in Zhu & Stone 2018). This will simplify the interpretation of angularmomentum budgets later.
2.1.2 Equation of state and cooling function
As a simplification, we assume the flow follows an ideal equation of state and is approximately locally isothermal, that is T ≈ T_{eff}. (R) where T_{eff}. is a prescribed radial temperature profile. This is achieved by solving the energy equation (8)
where we defined a heating and cooling function (9)
where τ is the cooling time that is set to 0.1 time code unit (see below) and Γ = 1.0001 is the polytropic index of the gas. The target temperature profile is (10)
where T_{0} is the midplane temperature at the inner radius Rint. This choice of cooling function allows us to enforce a chosen temperature profile that mimics the real radiative equilibrium and avoid the development of the vertical shear instability (VSI, Nelson et al. 2013), which would appear in a strictly locally isothermal approximation.
Since the gas is ideal, we can define an isothermal sound speed . It can be shown that as a result of the vertical hydrostatic equilibrium, c_{s} and Ω_{K} are related to the vertical disc thickness h(R) through (11)
Assuming the disc is at thermal equilibrium (T = T_{eff.}(R)), we have c_{s} ∞ R^{−1/2}, and hence the disc aspect ratio ɛ ≡ h/R is constant. For the following, we chose T_{0} in (10) so that ε = 0.1.
2.2 Numerical method and parameters
2.2.1 Integration scheme
The simulations were performed using the PLUTO code (Mignone et al. 2007) that solves the MHD equations with a conservative Godunovtype scheme and a secondorder RungeKutta time stepping. We used a HLLD type Riemann solver to compute the intercell fluxes. In order to ensure the solenoidal constraint (6), we used the constrained transport approach (Kane 1966; Evans & Hawley 1988). The implementation of ambipolar diffusion in the PLUTO code follows that of Lesur et al. (2014) and Béthune et al. (2017).
2.2.2 Code units and notations
The internal radius is Rint = 1, which sets the length code unit and is chosen to be 1 au while R_{ext} = 50. The time code unit is , which is set to 1/2π years so that GM_{*} = 1, with M_{*} = 1 M_{⊙}, M_{⊙} being 1 solar mass. Therefore, Ω_{K}(R) = Ω_{0} (R/R_{int})^{−3/2}=R^{−3/2}. We chose 300 g cm^{−2} as a unit for the surface density and express the accretion rate in M_{⊚} yr^{−1}. We denote code units as c.u. hereafter. We use the subscript X_{0} to indicate that the quantity X is considered on the midplane (θ = π/2) and the subscript X_{p} when X is a poloidal quantity.
2.2.3 Dimensionless numbers and definitions
We used the plasma parameter β to quantify the disc magnetisation, which is defined from the midplane properties of the disc as (12)
When considering the initial state of a given simulation, we refer to the initial magnetisation inside the cavity as β_{in} and to the initial magnetisation in the external part of the disc as β_{out}. The second key parameter of this study is the strength of ambipolar diffusion, quantified with the Elsasser number (13)
where v_{A} = B/(4πρ)^{1/2} is the Alfvén speed. We refer the reader to Appendix A for detailed information on the justifications of the model we adopted for Λ_{A} and how we modelled its spatial dependencies in our simulations. These two dimensionless numbers are the main control parameters of our study.
The disc refers to the whole part of the simulation that covers r ∈[1; 50] and z/R ∈ [−0.3; 0.3]. The cavity is the region where the surface density is reduced by a given factor in the innermost part of the disc (i.e. from r = 1 to r = 10 in most of the models). The external part of the disc or socalled outer disc refers to the region where the disc is full and described by a standard protoplanetary disc (without a drop in the density profile) and which extends from r ≈ 10 to r = 50. Finally, we refer to the region defined by r ≤ 1 of our disc as ‘seed‘, which is at play in our simulations through the inner radial boundary condition.
2.2.4 Computational domain
The radial direction is divided into 320 cells that expand from the inner radius r ≡ R_{int} to the external one r ≡ R_{ext} and are uniformly meshed on a logarithmically shaped grid. The colatitude domain is mapped on a stretched grid near the poles (from θ = 0 to θ = 1.279 and from θ = 1.862 to θ = π, with 72 cells in each zone), while the grid is chosen to be uniform around the midplane (from θ = 1.279 to θ = 1.862 with 96 cells) for a total of 240, which increases the precision in the region of interest. The disc scale height h is then covered by 16 points in the case where ε is fixed as constant and equal to 0.1.
2.2.5 Boundary conditions
Outflow boundary conditions are used in the radial direction so that no matter can come from the inner radius. In addition, we added a waveabsorbing zone for radii r < 1.5, which dampens poloidal motions on an orbital timescale. We detail the impact of this procedure in Appendix B.
In these 2.5D simulations, axisymmetric conditions with respect to the polar axis are enough to handle the boundaries for the colatitude direction. With the aim of reducing the impact of the outer boundary conditions, we focused on radii lower than 30.
2.2.6 Initial condition, wind, and cavity
The initial temperature profile is the effective temperature profile given in (10). The initial states for the density and the azimuthal velocity v_{ϕ} = RΩ_{K} mimic those of Nelson et al. (2013) to account for the hydrostatic equilibrium, while v_{r} = v_{θ} = initially. These profiles read, without taking into account the cavity, (14) (15)
with ρ_{0} being the density at the internal radius. We chose q = −1 and p = −3/2 for the Eqs. (14) and (15), which is consistent with selfsimilar stationary disc solutions (JacqueminIde et al. 2021).
The initial vertical magnetic field follows a power law B_{z} ∞ r ^{(p}+^{q)/2} so that the plasma β parameter in the unperturbed disc is constant. To ensure that ∇·B = 0, we initialised the magnetic field using its vector potential A, defined so that B = ∇ × A. Following Zhu & Stone (2018), we chose (16)
where m = (p + q)/2 = −5/4. This results in apoloidal magnetic field that depends on the radius only: (17)
The initial strength of the magnetic field is controlled by β_{out}, so that .
To add a cavity and mimic a TD, we multiplied the density profile by a function f that depends on the radius only, so that (18)
where Σ_{0}(R) ∞ R^{p+1} is a standard surface density profile for a protoplanetary disc. The a, b, and c coefficients are defined as
where R_{0} is the radius of the cavity (in code units), n the number of cells on which the transition spans, and nδR the corresponding length in code units. We note that while the density profile exhibits an inner ‘hole’, the magnetic field distribution is kept as a power law (17). As a consequence, the initial magnetisation β(R) also exhibits a jump in the cavity since P ∞ Σ(R).
Therefore, β_{out}/β_{in} is equal to the contrast in the gas surface density. In short, the function f creates a cavity in Σ but does not affect B_{p}. As a result, we simulate a TD with a strongly magnetised cavity (β_{in} = 1). A typical radial profile of the quantities discussed above is shown in Fig. 1.
2.3 Integration and averages
Several integrations and averages are used throughout the text. In this manuscript, we use the following proxy for the vertical integration along θ: (20)
θ_{±} quantify the integration height as shown in Fig. 2 so that (21)
with hint being the integration height at radius R given by an integration effective aspect ratio ε_{int} ≡ h_{int}/R. We note that this integration ‘height‘ is not necessarily the disc thickness h. We introduce as (22)
which corresponds to a thetaaveraged ‘effective’ midplane β plasma parameter. It is defined so that it matches the midplane β parameter in a hydrostatic isothermal disc. This more general definition is needed when the disc midplane is displaced vertically such as inside the cavity (see Sect. 3.5.3). Finally, we added the time average defined by (23)
We ran the 2.5D simulations so that we reach 1000 orbits at R = 10, which means ≈31000 orbits at R_{int}. If not specified, time averages were calculated taking into account the whole simulation without the first 4000 orbits at R_{int} to suppress the transient state. Otherwise, we indicate our choice of notation when needed, ⟨X⟩_{1000} being the timeaveraged value of X during the last 1000 orbits at Rint, for example.
Fig. 1 Initial and timeaveraged profile of Σ (top left panel), (top right), B_{Z,0} (bottom left), and ν_{φ,ο} (bottom right), with respect to r. 
Fig. 2 Schematic view of the disc represented in orange. θ_{±} define the vertical integration surface and h_{int} is the integration scale height at a given radius r. 
2.4 Simulation table
All the simulations performed are listed in Table 1. The characteristic parameters are the external initial plasma parameter β_{out}, the internal initial plasma parameter β_{inı} and the initial ambipolar Elsasser number ʌ_{A},_{o}. Additionally, we performed a convergence test by running a highresolution (640 × 480) simulation similar to the fiducial one that exhibits profiles that differ by less than 8% in the cavity and by less than 1 % when considering the entire domain.
3 Fiducial simulation
We start this section by describing our fiducial simulation in detail (β_{out} = 10^{4}, β_{in} = 1, Λ_{A} = 1 and R_{0} = 10), before turning to an exploration of the parameter space.
3.1 Evolution of surface density and plasma magnetisation
We first look at the temporal evolution of the surface density (Fig. 3). We find that the cavity stands during the whole simulation as its radius remains close to its initial value. As we show in Sect. 3.5.1, the cavity tends to expand slightly. The cavity location, defined as the radius where the surface density equals half of its maximum value, is subject to a small variation of Δr/r = 10.3% over the duration of the simulation. While the external disc is relatively smooth with respect to time, the cavity is striped by temporal variations of Σ, which may suggest that matter is moving inside the cavity at relatively fast speeds. We study these stripes in depth in Sect. 3.5.3. A small accumulation of material is seen close to the inner radius at R ≤ 1.5. We refer the reader to Appendix B for a quantitative discussion on this accumulation.
Figure 3 also pictures the evolution of whose results are similar to the ones for Σ. Inside the cavity, exhibits a stripedlike pattern with an accumulation close to the internal radius. The edge of the cavity is not smooth at all but varies around its initial value of 10. Though stays on average around 1 in the cavity, some low values around 10^{−2} are reached from time to time. After approximately 4000 orbits at the internal radius, both Σ and reach a quasistationary state.
Gaps and rings are detected in the outer part of the disc, in the spatiotemporal diagram of both Σ and β (Fig. 3). We also emphasise that these structures are observed in all of our simulations (see Figs. 19, 21, and 25). Regarding the fiducial simulation, we detect two main gaps after the cavity edge and before R = 30. For better visibility, we show the surface density and the vertical magnetic field, time averaged on the last 1000 inner orbits and a focus in the region R = 1218 au (where the gaps are detected) in Fig. 4. Gaps are characterised by a drop of ~5% of the local surface density and their location is correlated with a sharp increase of the vertical magnetic field, which matches the secular wind instability described by Riols et al. (2020). These structures are enhanced in the simulation with a higher ambipolar Elsasser number, as can be seen in Fig. 20. In addition, we observe the merging of gaps on longer timescales (Fig. 19) similarly to Cui & Bai (2021). While of interest for the dynamics of the outer disc, we did not address the evolution of these rings and gaps any further and instead focused on the dynamics of the cavity.
Simulation information. B4BinOAmO is the fiducial simulation.
Fig. 3 Surface density (left panel) and plasma beta parameter β (right panel) as a function of R at midplane and time for the fiducial simulation. The cavity remains during the entire simulation and keeps a relatively strong magnetisation with . 
Fig. 4 Radial profiles of the surface density and of the vertical magnetic field in the midplane, timeaveraged on the last 1000 orbits at R_{int}. The vertical magnetic field is vertically averaged, and both profiles are given in arbitrary units. 
3.2 Disc structure
3.2.1 Magnetic structure
We show the timeaveraged magnetic field in Fig. 5. In the cavity, the poloidal magnetic field lines are pinched at the midplane, but they remain vertical in the outer disc. These two regions are separated by a transition zone located at the cavity edge that exhibits a magnetic loop. Inside this loop, the polarity of the azimuthal component is reversed, with Β_{φ} > 0 in the upper hemisphere close to the disc. The poloidal field lines present an elbowshaped structure above and below the transition with significant changes of direction at h_{int}/R ~ ±0.3, ±0.6, and ±0, 9.
3.2.2 Velocity streamlines
We show the timeaveraged density and streamlines in Fig. 6. The disc clearly appears around the midplane for R ≲ 10, while the depleted profile in ρ indicates the cavity for R < 10. We find that a wind is emitted from the cavity, with poloidal streamlines approximately parallel to magnetic field lines, as expected from ideal MHD. A closer inspection of the streamlines, however, shows that in the regions close to the transition radius R ≳ 8, matter is falling into the cavity. Figure 6 shows that this material is actually coming from the outer disc. It is originally ejected from this disc, before being deflected and accreted into the cavity, generating an elbowlike shape similar to the one found for magnetic field lines (Fig. 5). This accretion stream then stays localised close to the cavity midplane down to the inner radius of the simulation. In the outer disc, the motion of the gas is not as well organised, though it is approximately symmetric with respect to the midplane.
3.2.3 Angularmomentumflux streamlines
In order to deepen the analysis of the role of the magnetic structure, we concentrated on the timeaveraged angularmomentum flux, defined by (24)
The poloidal flux lines associated with this angularmomentum flux are shown in Fig. 7. It appears that angular momentum is extracted from the disc midplane and carried both radially and vertically in a relatively homogeneous manner. In particular, we note that there is no elbowlike shape for the angularmomentum flux, in contrast to the magnetic and velocity streamlines, indicating that the cavity + outer disc system has adapted its magnetic topology to transport angular momentum homogeneously.
Fig. 5 Timeaveraged poloidal magnetic field lines and toroidal field component ⟨Β_{φ}⟩ for the fiducial simulation. We note the peculiar field topology close to the truncation radius. 
3.3 Accretion theory
3.3.1 Accretion rate
The first step to study the accretion in the disc is to investigate the accretion rate Ṁ defined as (25)
The height over which pv_{r} is integrated has a direct influence on Ṁ, mostly because of the elbowshaped streamlines. It is then useful to change the thickness of the integration domain, which is controlled by the parameter ε_{int} = tan [(θ_{+} − θ_) /2]. Results are presented in Fig. 8 for three values of ε_{int}. For ε_{int} = 0.3 and around R = 10, the accretion rate is close to zero, indicating that the gas does not plunge directly in the cavity from the disc midplane. This radius corresponds to the location of the basis of the elbowshaped loop along which the gas is moving. Averaging higher above the disc allows us to cancel out this effect. Moving to ε_{int} = 0.6 and 0.9, the accretion rates in the disc and in the cavity eventually match by less than 50%, despite a jump of more than two orders of magnitude in Σ. This clearly indicates that the accreted material effectively ‘jumps’ above the transition radius, and that a steady state is reached with the whole system (cavity + outer disc) accreting at a constant rate.
The fact that the accretion rate is approximately constant while the surface density decreases by two orders of magnitude implies that the accretion speed should increase dramatically. This is clearly visible in Fig. 9, which shows the radial profile of the accretion speed v_{acc} for ε_{int} = 0.9, defined by (26)
This velocity profile exhibits a welldefined transition between subsonic accretion outside the cavity with ⟨v_{acc.}⟩ ~ 10^{−3} ⟨c_{s}⟩ and transsonic accretion inside with ⟨v_{acc.}⟩ ~ ⟨c_{s}⟩.
Fig. 6 Timeaveraged streamlines and density for the fiducial simulation. We note the peculiar shape of the streamlines around the transition radius. 
3.3.2 Governing equations for accretion
Accretion theory can be understood as the secular evolution of Ṁ and Σ. In systems driven by MHD processes, these two quantities are usually supplemented by the magnetic field B_{z} threading the disc. We apply the vertical integration procedure to the mass and angularmomentum conservation equations, which become (27) (28)
We defined W_{Rφ} and W_{θφ,} respectively, as the radial and surface stresses by (29)
We remind the reader that we use a peculiar definition of the velocity deviation υ so that no additional surface terms appear in Eq. (28). In order to take into consideration the role of the magnetic wind, we complete this set of equations by the vertical magnetic flux conservation (30)
Fig. 7 Timeaveraged angularmomentumflux streamlines over timeaveraged density for the fiducial simulation. Angular momentum leaves the disc midplane because of the wind. 
Fig. 8 Accretion rates for different integration height scales with respect to the radius inside the disc. The higher ρv_{r} is integrated the closer to a constant value Ṁ is in the cavity. The average value inside the cavity (from R = 1 to R = 10} is Ṁ = 1.4 ± 0.2 × 10^{−7} M_{⊙} yr^{−1}. 
Fig. 9 Accretion speed for ɛ_{int} = 0.9 in units of local sound speed c_{s}. The profile exhibits a clear transition between subsonic and transsonic accretion that occurs where the edge of the cavity is located. 
3.3.3 Mass conservation and masslossrate parameter
The mass conservation equation is given by Eq. (27). Figure 10 shows the mass conservation for ɛ_{int} = 0.9 with timeaveraged quantities. The first information is that inside the cavity, the time derivative of Σ is close to zero, meaning the simulation reaches a steady state up to r ~ 8. Closer to the cavity edge, we note that this same term is negative, which is linked to the slow expansion of the cavity; this is discussed later in Sect. 3.5.1.
The main contribution of the wind mass loss is located in the cavity at r < 5 and is completely compensated by the radial derivative of the accretion rate. Additionally, the ‘wind’ mass flux turns negative around the cavity edge, which is due to matter being accreted from the outer disc atmosphere (see the ‘elbowshaped structure’ in the poloidal streamlines).
In order to quantitatively account for the role of the wind, we constructed the masslossrate parameter ζ = ζ_{+} + ζ_ (Lesur 2021b), where ζ_{+} and ζ_ are defined by (31)
and the corresponding quantities are time averaged. The signs of ζ_{±} are chosen accordingly so that a positive value of ζ_{±} corresponds to matter leaving the surface at θ_{±}. Since ζ_{+} and are pretty much symmetric with respect to the midplane, we focus on ζ only. The results are illustrated in Fig. 11 where both ⟨ζ⟩ and –⟨ζ⟩ are shown. In order to compare with selfsimilar models (Lesur 2021b), we studied the values of ⟨ζ⟩ at zo = 6 h, which corresponds to ɛ_{int} = 0.6. The masslossrate parameter is approximately constant in the external part of the disc around 6.2 × 10^{−5}, while it peaks at 2.9 × 10^{−2} in the inner part. We find two zones where ⟨ζ⟩ < 0. One is close to the inner boundary and probably a boundary condition artefact, while the other extends from R ≈ 5 to R ≈ 17 au and is related to the material falling down on the disc around the transition zone, such a contribution being notably stronger for ɛ_{int} = 0.6.
To compare to selfsimilar solutions, we show the selfsimilar scaling of the masslossrate parameter with respect to ⟨β⟩ derived by Lesur (2021b), which reads ⟨ζ_{self}⟩ = 0.24 ⟨β⟩^{−0.69}. It comes as no surprise that this fit does not account for negative values of ⟨ζ⟩ since these are due to the transition radius, which is not selfsimilar.
The wind masslossrate parameter is smaller than the selfsimilar scaling in the outer disc by a factor of a few. This discrepancy is probably due to the influence of the cavity magnetosphere that compresses the disc magnetosphere, resulting in a deviation of ζ from the selfsimilar result. Moreover, it seems that the further we move outward, the closer we are to the selfsimilar values, indicating that we recover selfsimilar scalings far ‘enough’ from the cavity, as expected.
In the cavity, ζ is significantly weaker than expected from a naive extrapolation of selfsimilar scaling laws. This indicates that the massloss rate saturates at β ~ 1, a regime which was not explored by Lesur (2021b).
An alternative model to the selfsimilar one is used to describe ⟨ζ⟩ with greater accuracy. The selfsimilar fit is kept for the external parts of the disc . Another one is then calculated for the inner part only, , so that the final profile is given by (32)
We obtain a_{int} = −0.20 and ζ_{0 int} = 0.018. Such a model, with a_{ext} – a_{int} < 0, allows us to recover both of the previous regimes with a reasonably accurate depiction of the disc. The final profile exhibits a transition occurring at β_{t} ≈ 5, which is close to the lowest value of the ones used to build the selfsimilar fit in Lesur (2021b). The final curves are rendered in Fig. 11. The fit does not account for the negative values, but it properly catches both the inner and external parts of the disc.
Fig. 10 Mass conservation for ɛ_{int} = 0.9. The time derivative of Σ remains perfectly constant and equal to zero inside the cavity and only becomes moderately negative at the cavity edge R = 10. This suggests that the simulation indeed reaches a steady state for radii up to ≈8. The mass conservation is also correctly recovered. The three lines do not add up to zero because we used a moving average for better visibility, and the quantities are timeaveraged on a sample selection of output flies that do not contain all the time steps computed by the code. 
Fig. 11 ⟨ζ⟩ parameter for ɛ_{int} = 0.6. The selfsimilar fit shown here is for comparison only and was also obtained for ɛ_{int} = 0.6. It appears that it is coherent for the external disc, while it predicts a wind way too intense in the internal part, therefore, another model is used to describe ⟨ζ⟩ in the whole disc. 
3.3.4 Angularmomentum conservation
In Fig. 12, we show the terms involved in the angularmomentum conservation Eq. (28). These are time averaged and multiplied by r^{−3/2} for beugŗ readability.
The integration height is ɛ_{int} = 0.9 and chosen so that the influence of the cavity edge is diminished. In contrast to the massconservation equation, the time derivative is negligible. The surface stress (‘wind’) removes angular momentum from the whole disc with a major contribution right after the cavity at R ≈ 13. We also observe that the radial stress is always positive, except at the cavity edge.
Such a cancellation suggests that two accretion regimes are observed in the disc, which echoes the radial profile of both the accretion rate and speed. To characterise the radial stress term, we introduced the α parameter of Shakura & Sunyaev (1973). It must be noted that the origin of this stress is in no way solely linked to turbulence and is considerably driven by the laminar structure of the magnetic wind. Appendix D details the origin of the stress and sheds light on the turbulent versus laminar origin of α Nevertheless, the α parameter can still be used in this wind model, whose definition when time averaged is (33)
The corresponding profile is given in Fig. 13, where ɛ_{int} = 0.9. In the external part of the disc, ⟨α⟩ = 49 ± 5 × 10^{−4}, while it reaches a maximum value inside the cavity ⟨α⟩ = 13 + 5.
Following a similar procedure as the one for α, we define a dimensionless number associated with the surface stress component, v_{w}. As for ζ, we define v_{w,±}, which were chosen to be positive for angular momentum leaving the disc on both sides: (34)
We show the dependence of v_{w} on R in Fig. 13. In the external disc, ⟨v_{w}⟩ = 2.3 + 1.1 × 10^{−4}, while it rises up to ⟨v_{w}⟩ = 1.0±0.1 inside the cavity. The same observations are drawn for both ⟨α⟩ and ⟨v_{w}⟩ as for ⟨ζ⟩. Therefore, two separate regimes are at stake in the disc. The outer disc regime is typical of windemitting protoplanetary discs, with transport coefficients close to the ones found in selfsimilar wind models for β ~ 10^{4}, indicating that the dynamical properties of the outer disc are not perturbed by the presence of the cavity. On the contrary, the second regime describes the inner part of the disc with fast accretion and high values for α and v_{w}, which are both of the order of unity. Table 2 displays the transport coefficient values for all the simulations.
Transport coefficients for a subset of simulations.
Fig. 12 Angular momentum conservation multiplied by r^{−3,2} and time averaged. Full blue line is , red dotdashed line is , green dashed line is , and purple dotted line is . 
3.4 MHD wind
It is well known that steadystate MHD winds in ideal MHD can be characterised by a set of MHD invariants (Blandford & Payne 1982), which are conserved quantities along each poloidal field line (Fig. 5). In our axisymmetric simulations, a steady state is approximately achieved above the disc, in the ideal MHD region. Hence, we can measure these invariants on field lines attached in the cavity and in the outer disc.
In the following, we select a field line anchored in the disc midplane at R = R_{w}. The corresponding Keplerian angular velocity is Ω_{w}. while B_{w} is the poloidal magnetic field at the midplane. We then consider the following invariants, built on timeaveraged quantities and listed in Table 3.
The massloading parameter accounts for the quantity of matter that escapes the disc with the wind and is defined as (35)
The rotation parameter is given by (36)
The magnetic lever arm accounts for the angular momentum that is removed from the disc by the wind is defined as (37)
Of course, these invariants echo the transport coefficient definitions previously used to describe the disc, and one expects Κ ≈ β ζ/4 ɛ (Lesur 2021a).
To compute these invariants, we arbitrarily chose one field line in the cavity (referred to as ‘in’) leaving the midplane at R_{in} = 3.5 au and one in the external disc (referred to as ‘ext’) leaving the midplane at R_{ext} = 15 au (see the first panel of Fig. 14). We note that the disc thickness affects the MHD invariants since the physical foot points of the field lines are not located at the midplane but slightly above. Such a limitation especially concerns the field lines in the external disc, which are subject to a largescale oscillation close to the transition radius. Therefore, the calculated MHD invariants are subject to caution and we only draw general conclusions regarding the nature of the wind.
We show the invariants along the chosen field lines in Fig. 14. We find that all of the invariants remain reasonably constant once high enough above the disc, as expected from a steadystate ideal MHD flow. The wind launched from the cavity is different from the disc one. The cavity wind has a much weaker mass loading parameter and a much larger lever arm (by almost a factor 10). We also find that its rotation parameter differs significantly from 1, indicating that field lines are rotating at 80% of Ω_{ĸ} in the cavity. This point is probably related to the fact that the disc itself is subKeplerian in this region (Fig. 1). Quantitatively, we find ĸ_{in} = 2.2 × 10^{−2}, ĸ_{ext} = 2.5, λ_{in} = 23 and λ_{ext} = 3.2. These values are coherent with the transport coefficients computed in previous sections. We also note that the values of K and λ in the cavity match some of the historical solutions of Blandford & Payne (1982) (see their Fig. 2), which correspond to superAlfvénic and collimated outflows. These values are also consistent with the magnetic outflow solutions of Ferreira (1997) (see Fig. 3). Hence, the cavity we find quantitatively matches the inner JED proposed by Combet & Ferreira (2008).
Fig. 13 Timeaveraged transport coefficients ⟨α⟩ and ⟨v_{w}⟩ for ɛint = 0.9. 
MHD invariants for a subset of simulations, computed with timeaveraged quantities on the last 1000 orbits at R_{int}.
Fig. 14 General structure of the wind and MHD invariants. First panel shows the field lines in the internal and external disc. The grey dashed lines represent the surfaces at h/R = 0.3; 0.6 and 0.9. The three other plots display the MHD invariants for the internal field line (full, light blue) and the external one (semidashed, dark blue). The invariants are time averaged on the last 1000 orbits. 
3.5 Temporal evolution
We observe two kinds of time variability in the fiducial simulation: a secular variability responsible for the slow expansion of the cavity, and a short timescale variability responsible for the striped patterns observed in spacetime diagrams (Fig. 3). At this point, we started our exploration of time variability by focusing on the secular evolution, beginning with a discussion of the cavity expansion.
3.5.1 Slow cavityedge expansion
As previously mentioned, the cavity edge moves slowly outwards during the simulation. Neglecting the impact of the wind in terms of massloss rate at the cavity edge location, which is coherent with Fig. 10, and assuming piecewise constant accretion rates and surface densities across the cavity edge, one obtains (38)
where is the cavity edge ‘velocity’ and δṀ and δΣ are the jump in accretion rate and surface density at the cavity edge. By calculating Ṁ and Σ around R_{0}, we find while directly evaluating the cavity edge motion yields (both in eu.). Therefore, the cavity is expanding because of the slight mismatch in accretion rate observed in Fig. 8.
3.5.2 Magnetic field transport
To interpret the time evolution of the magnetic field, we studied the transport of magnetic flux inside the disc and defined a flux function ψ such that (39)
Assuming the total flux is constant with respect to time, the isocontours of ψ describe the motion of the magnetic field lines in the disc plane. The spatiotemporal diagram for ψ is shown in Fig. 15. The magnetic flux is advected slowly towards the star in the external disc, while it tends to diffuse outwards from the inner part of the disc to the cavity edge. The poloidal magnetic field lines in Fig. 5 show that ⟨B_{z,0}⟩ < 0 in the transition region (8 ≲ R ≲ 12) and ⟨B_{z,0}⟩ > 0 otherwise. This transition region is recovered in Fig. 15 as a region where ∂_{r}ψ < 0.
Overall, we observe that the negative field of the transition region is diffusing outwards, while the positive field of the outer disc is advected inwards. We therefore observe a reconnection of the largescale field around R ≈ 12, which progressively ‘eats’ the negative field of the transition region. In addition to this, we observe that field lines deep in the cavity also diffuse outwards.
To achieve a quantitative estimate of the field line advection speed, we first note that the evolution equations for ψ read (40)
Following Guilet & Ogilvie (2014), we rewrote these evolution equations as an advection equation for ψ: (41)
where we defined the ‘field advection velocity’: (42)
Eventually, we define a dimensionless advection parameter v_{B} = v_{ψ}/v_{K} which quantifies the advection speed (Bai & Stone 2017). In this framework, positive values of v_{B} imply an outward transport field, while negative values trace inward field transport.
We show the radial dependence of ⟨v_{b}⟩ in Fig. 16. In the external disc, we find that the magnetic field is advected inwards with a velocity of υ_{ψ} = −2.6 × 10^{−3} v_{K}. v_{B} changes its sign multiple times in the cavity, but it remains negative close to the cavity edge, between R ≈ 7 and R ≈ 11 au, where υ_{ψ} = +3.2 × 10^{−3} v_{K} Such a result is in accordance with Fig. 15 and indicates that field lines are converging at the transition radius with opposite vertical polarity. In the external parts, vr is negative and υ ψ = −2.6 × 10^{−3} v_{K} so that vertical magnetic field pointing upwards is advected. We note that this inwards advection of the outer disc field lines is in sharp contrast to other work that focused on ‘full’ discs (Bai & Stone 2017; Lesur 2021b). We revisit this discrepancy in the discussion.
Fig. 15 Flux function ψ(R, t) for the magnetic field, taking into account the flux at the surface of the seed and the radial flux. For radii larger than the one of the cavity, the field lines are advected toward the centre during the whole simulation. At R ≈ 13, the flux accumulates and exhibits a striped structure for smaller radii. 
Fig. 16 Magneticfield transport parameter v_{B} as a function of radius. We note that the outer disc is transporting magnetic field lines inwards (v_{B} < 0). 
Fig. 17 Temporal evolution of Σ in dotted green, Ṁ in dashed blue, and B_{z} (vertically averaged) in black (full line) at R = 3 au for the fiducial simulation. ζ is calculated at ɛ_{int} = 0.3 and shown by a semidashed red line with a logarithmic scale on the right of the panels. Apart from ζ, all the profiles are given in arbitrary units and divided by their maximum value reached during the timescale of the top panel. 
3.5.3 Fast variability of the cavity
Up to this point, we mostly considered timeaveraged quantities and ignored fast variability. While our numerical solutions are quasisteady if one looks at averages on 100s of orbits, they also exhibit a fast time variability (see the temporal stripes in Fig. 3), the origin of which requires clarification.
Figure 17 shows such a temporal evolution of Σ, Ṁ, and B_{z} at R = 3. These profiles encounter sharp fluctuations chaotically distributed over time. Therefore, the cavity is subject to bursts of matter that quickly falls onto the star (the typical width of a peak is ~5 orbits at R_{int}, which is still far larger than our temporal resolution). This variability explains the stripes seen in the spatiotemporal diagram (Fig. 3).
We focus on a few of these bursts in the bottom panels of Fig. 17, while instantaneous pictures of the density corresponding to the (b) panel are given in Fig. 18. For these bursts, we see that the local maximum values of B_{z}, Σ, and Ml are correlated. When an inflow of matter crosses the cavity, Σ peaks as well as Ṁ, which in turn increases ζ. In terms of temporal sequence, it seems that B_{z} increases slightly before Σ and Ṁ, which would indicate that B_{z} is the driver of these bursts, but we cannot be definitive on this sequence because of the lack of temporal resolution. Finally, we observe that ζ is always clearly delayed compared to the other quantities, indicating that the wind inside the cavity ejects more material once the bubble of material has passed.
For a more precise insight on accretion and temporal variability, we refer the reader to Fig. 18, which shows the density and poloidal magnetic field lines at different times. In the first panel, we see a filament of matter located above the disc that extends from (R = 10, Z = 5) to (R = 15, Z = 10). This structure is cut in two on the second panel, revealing two bubbles of matter, one being about to fall, while the other is about to be ejected and leave the disc in the wind. Concerning the filament and the bubble formation, we detect a current sheet at the location of the filaments, where the total magnetic field cancels out (Bφ, = 0 at the edge of the magnetic loop and B_{p} = 0, because two antiparallel poloidal field lines meet at the elbowshaped structure location). It is therefore a possibility that these structures form due to magnetic reconnection. Focusing on the falling material, we see it reaching the edge of the cavity on the third panel before crossing it on the next one. When the gas crosses the cavity, the disc oscillates locally above and below the midplane and is therefore highly dynamical. With a slight delay (last three panels), we see an outflow emerging from the cavity, and the wind density increases. Such an observation exhibits the link between wind and accretion (see Fig. 10). The ejection of gas from the cavity is not constant with respect to time and occurs occasionally with burst events for which ζ eventually peaks at 0.1. This explains why the effective value of ⟨ζ⟩ is lower than the one predicted by selfsimilar models for which the ejection is continuous with a higher massloss rate parameter.
Combining 6 and 18, we unveil a general scheme for feeding the cavity. First, the gas located inside the outer disc elevates from the midplane up to approximately two local disc heights and organises itself in a filamentary way. Then, bubbles of matter fall and cross the cavity, forming the elbowshaped structure on the timeaveraged profile.
3.5.4 Magnetic Rayleigh Taylor instability
To account for the formation and stability of the bubbles of matter at the cavity edge, we explored the possibility of having a magnetic Rayleigh Taylor instability (RTI; or interchange instability) in the cavity. The disc is geometrically thin inside the cavity and the density is relatively continuous radially. Under these conditions, we referred to the analyses of Spruit & Taam (1990), Spruit et al. (1995) and Stehle & Spruit (2001), which assume an infinitely thin disc. We reformulated the instability criterion of Spruit et al. (1995) (see their Eq. (59)) in terms of the plasma parameter in Appendix C. The resulting criterion C.9 states that a necessary condition for the occurrence of the RTI is . Figure 3 shows that is of the order of 0.1 in the cavity and rarely goes beyond this value, except for very short periods of time, for instance during the accretion 'bursts'. We conclude that the cavity β plasma parameter is too large to sustain the RTI on average, but we cannot exclude that it could be triggered in the rare excursions where the cavity reaches , as it does during some of the bursts.
Fig. 18 Density and magnetic field lines at different times showing the advection of a bubble of material (arrow) from the outer disc through the cavity. 
Fig. 19 Σ(R, t) and of B4Bin0Am1. The cavity stands throughout the simulation, though its edge falls down to R ≈ 2 au during the transient, before broadening up to R ≈ 4 au in a few thousands of orbits at R_{int}. The profile of is characterised by the presence of gaps in the external disc. 
4 Parameter space exploration
4.1 Ambipolar Diffusion
We checked the influence of Λ_{A,0} in the simulation B4Bin0Am1. This simulation is the same as the fiducial one, except for the initial value of Λ_{A} which is set to 10.
4.1.1 General structure of the disc and gaps
The spatiotemporal evolution of Σ and ß are shown in Fig. 19 for Λ_{A} = 10. During the transient state, the cavity edge falls down to R ≲ 2 au before expanding back up to R ≳ 4 au in a few thousands of orbits at R_{int}. Overall, the transient state lasts for a shorter period of time than in the fiducial run, and the cavity extension is smaller.
We observe the apparition of gaps in both the profiles of Σ and (Fig. 19) located in the external disc and broadening with time. Such structures are observed in numerous occasions in protoplanetary discs simulations either with ideal (JacqueminIde et al. 2021) or nonideal MHD (Béthune et al. 2017; Suriano et al. 2019; Riols et al. 2020; Cui & Bai 2021). We observe that gaps are associated with low ß regions and are localised relatively far from the disc's inner boundary. Some gaps merge with one another, meaning that only three of them remain after 15 000 orbits at R_{int}, similarly to what was found by Cui & Bai (2021). We reserve the study of the interaction between these gaps and the cavity for a future paper.
Figure 20 shows the flow and field topology for Λ_{A} = 10 as well as the timeaveraged magnetic structure of the disc. The main features of the fiducial simulation are recovered, namely the elbowshaped structure and the associated magnetic loop. These are, however, located closer to the star, the cavity radius being smaller in this simulation.
In contrast to the fiducial simulation, the outer disc is this time top down asymmetric, which has an impact on the shape of the elbow above and below the disc plane. The elbow is prominent above the disc but almost disappears below it, except for a small set of streamlines close to the cavity. The magnetic field lines exhibit a local slanted symmetry in the external disc at the gaps location. This is similar to the topology observed in ambipolardominated discs (Riols & Lesur 2018, 2019). The gaps seem to be characterised by small vortices in the (r, θ) plane, located at the disc surface at the corresponding radii, indicating a meridional circulation.
4.1.2 Transport coefficients and wind invariants
The accretion rate remains constant in the whole disc with a value close to 1.2 × 10^{−7} M_{⊙} yr^{−1} and an accretion velocity that is still subsonic in the outer disc and peaks up to 2 c_{s} at the internal radius. Therefore, the accretion picture is identical to the one for the fiducial run with an internal transsonic regime connecting through the cavity edge to a weakly magnetised wind.
Regarding the wind, we obtain a highly massloaded field line in the external disc that removes little angular momentum (λ_{ext} = 1.5 and K_{ext} = 5.0) and a lighter one in the internal disc that carries a massive load of angular momentum (λ_{in} = 4.9 and K_{in} = 0.24). We note that the disc wind is overall less magnetised and more massive, while the general picture of the fiducial run remains. The rotational invariant contrast is higher than in the fiducial simulation, its internal value being three times lower and the external one three times higher.
4.2 Influence of the initial plasma parameter
We studied the impact of the plasma parameter by varying both its internal β_{in} and external β_{out} initial value.
4.2.1 Role of the external initial plasma parameter
We explored how the outer disc magnetisation impacts the general properties of the system. We varied the initial value of ß between β_{out} = 10^{3} (run B3Bin0Am0) and β_{out} = 10^{5} (run B5Bin0Am0).
With regard to B5Bin0Am0, the spatiotemporal evolutions of Σ and are shown in the left panels of Fig. 21. Right at the beginning of the simulation, a burst of matter appears in the cavity, which is subsequently refilled. Its radius then remains fixed at ~4 au until other bursts happen at ~17 400 and ~27 000 orbits at R_{int}. Such local events do not dramatically change the general properties of the disc, which is similar to the fiducial one overall.
The bursts of matter (at ~17 400 and ~27 000 orbits at R_{int}, assuming the first one is due to the initial transient) give the illusion that some gas might be created inside the cavity, challenging mass conservation. These bursts are actually due to gas accumulating at the boundary of an accretion 'barrier'. We refer the reader to Appendix B for a more detailed description of these bursts. For now, we point out that these bursts highlight a limitation of our model regarding the implementation of the inner boundary conditions, but they only occur in the weakly magnetised (β_{out} = 10^{5}) simulations. Lastly, we add that estimating for this simulation is too difficult, since the cavity edge barely moves during the entire simulation.
In contrast to the run B5Bin0Am0, the cavity in the simulation B3Bin0Am0 quickly expands up to R ≈ 15 au and keeps growing throughout the simulation, faster than in the fiducial run (see the right panels of Fig. 21). We estimate its velocity as c.u., which is about three times faster than the fiducial run. We obtain δṀ ≈ 2.3 × 10^{−4} and δ Σ ≈ 4.2 × 10^{−2} (both in code units), so Eq. (38) gives c.u., where we chose R_{0} ≈ 20. The simple model we used seems to overestimate the widening velocity of the cavity but still gives the correct order of magnitude.
The timeaveraged surface density from the fiducial run, B3Bin0Am0, and B5Bin0Am0 are shown in Fig. 22, which demonstrates that the size of the cavity is ruled by the initial external plasma parameter. The lower β_{out} is, the wider the cavity gets when the disc reaches a steady state. On the contrary, the plasma parameter inside the cavity does not depend on its external structure and converges to β_{in} ≲ 1 in all of these simulations (Sect. 4.2.2 tackles this observation in depth). Once the transient state is gone, we note that the cavity expands faster for lower β_{out}. This can be understood using Eq. (38), which can be recast as (43)
where we defined the accretion velocities v_{acc} ≡ Ṁ/2π R_{0}Σ, and assumed Σ_{out} ≫ Σ_{in}. The expansion speed is then controlled by the term in parenthesis, since the accretion velocity in the cavity is always sonic (see 4.2.2). It is well known that the accretion velocity in the outer 'standard' disc is a decreasing function of β. Writing v_{acc}. ∝ β^{−σ} with σ > 0, Lesur (2021b) proposes σ = 0.78 and Bai & Stone (2013) σ = 0.66, which indicates that 0 < σ < 1. Assuming that a value exists for which , we obtain the following scaling: (44)
where we used the fact that in our setup. The relation (44) shows that for , we have approximately , indicating that the cavity expansion speed should increase as β_{out} decreases, which is precisely what we observe for B3Bin0Am0. For , we obtain, on the contrary, , showing a change of sign (hence a contraction of the cavity), albeit with a reduced speed. This regime might correspond to B5Bin0Am0, indicating that .
Fig. 20 Timeaveraged structure of the disc for B4Bin0Am1. Left panel: poloidal streamlines and density. Right panel: magnetic structure of the disc with magnetic poloidal field lines and 〈Β_{ϕ}〉. 
Fig. 21 Spatiotemporal diagrams for Σ and for B5Bin0Am0 and B3Bin0Am0. The cavity expands more in B3Bin0Am0 but shrinks in B5Bin0Am0. Bursts of matter occur in B5Bin0Am0 (see main text for more explanations). 
Fig. 22 Impact of initial external magnetisation on the surface density. We average the profile on the last 1000 orbrts at R_{int}. For B5Bin0Am0, we average on 1000 orbhs at R_{int} occurring between the 2 burst events seen in Fig. 21. 
Fig. 23 Impact of internal initial magnetisation on the plasma parameter for β_{out} = 10^{4}. 
4.2.2 Role of the internal initial plasma parameter
To study the impact of β_{in}, we ran a set of simulations that cover all the possible initial gaps β_{in}/β_{out} where logβ_{out} ∈ {3, 4, 5} and . We compared each result to the one obtained with β_{in} = 1 and the corresponding value of β_{out} A striking result is the fact that the disc's inner structure does not depend on β_{in}. No matter which β_{in} we initially choose, a transition occurs in the cavity in order to impose β_{in} ≈ 1. Interestingly, this threshold value is the one required to achieve transsonic accretion, which is mentioned in Wang & Goodman (2017). We itlustrate this statement with Fig. 23 for the particular case of β_{out} = 10^{4}. We focus on the transient state of B4Bin3Am0 in Fig. 24. The transition is due to matter leaving the cavity because of the fast accretion at stake after a sharp increase of the magnetic field (and therefore a decrease of β). This reorganisation of the cavity is a consequence of a rapid advection of magnetic flux from the cavity onto the seed, which initially has low magnetisation because of our initial setup. Due to the total magnetic flux conservation, there is a shortage of magnetic flux inside the cavity, up until the inner seed reaches a state where its magnetisation is almost constant. The magnetic field then accumulates at the inner boundary and β decreases accordingly so that accretion is enhanced. At this point, matter leaves the cavity as it is accreted onto the star. It is then clear that the cavity converges towards the same overall structure as the fiducial simulation one.
We note that taking β_{in} equal to ß_{out} would simulate a full disc with no cavity. Hence, a threshold should exist regarding the value of β_{in;} above which no cavity is able to form. Considering Fig. 23, it seems that this threshold is ≳10^{3}.
From these observations, we deduce that the cavity is regulated by the value of the plasma parameter, which must take a value close to 1 . The reason for this regulation is not entirely clear, and we add a word of caution regarding the role of the inner radial boundary condition, especially with respect to the magnetic field transport at R_{int}. We discuss this influence in Appendix B. Referring to Sects. 3.5.3 and 3.5.4, we suggest that the RTI may be responsible for this regulation, but a dedicated study would be required to ascertain this claim.
Fig. 24 Spatiotemporal diagrams of Σ (first panel), B_{z,0}, the vertical magnetic field at the midplane, and ψ, the flux function defined in Eq. (39). These profiles focus on the first orbits of the run B4Bin3Am0. 
4.3 Zoom with a larger cavity radius
We performed a simulation with a doublesized cavity (R_{0} = 20 au) in order to check the impact of the cavity size. The simulation was integrated for 1000 orbits at R = 10 au so that it reaches 355 orbits at R = 20 au The general observations are confirmed, such as the elbowshaped structure, the magnetic loop, the magnetic field advection in the outer disc, as well as the conclusions regarding the accretion. While the cavity size is identical to B3Bin0Am0, the behaviour of the disc is exactly the same as the fiducial one (see Fig. 25), indicating that β_{out} is the main parameter regulating the cavity expansion. This means that the global picture where two types of discs are connected is robust and not linked to limitations in the cavity size or artefacts due to the inner boundary condition.
5 Discussion and comparison with previous work
We modelled TDs sustained by MHD winds by performing 2.5D global simulations. This model acts as a proof of concept, showing that steadystate discs with both a cavity and a wind can be obtained. The resulting simulated discs are characterised by two different zones with contrasted dynamics.
First, our 'outer disc' behaves like a standard, weakly magnetised, ambipolardominated protoplanetary disc (Lesur 2021b; Cui & Bai 2021). In particular, we find mass and angularmomentum transport coefficients, wind properties, and accretion rates comparable to those found in the literature for 'full' discs. We also find weak gaps, which are characteristic of nonideal MHD discs (Riols & Lesur 2019; Riols et al. 2020). However, the magnetic field transport in the outer disc differs from previous studies: we find that magnetic field lines are advected inwards in the outer disc, in contrast to measurements in full discs which always show outwards transport (Bai & Stone 2017; Gressel et al. 2020; Lesur 2021b). This discrepancy is likely due to the fact that the field lines in the cavity are more collimated (i.e. less opened), which results in a lower pressure on the magnetic surfaces in the outer disc, but it is possibly also connected to the peculiar elbowshaped magnetic surfaces at the transition radius. In any case, it points to the fact that magnetic field transport is a nonlocal phenomenon: it depends on the global disc structure.
In contrast to the outer disc, the cavity (or inner disc) is strongly magnetised (β ≈ 1) because of its low surface density. We emphasise here that the absolute magnetic field strength in the cavity is not stronger than standard protoplanetary disc models. In practice, and given our set of units, we have B_{0} ≈ 0.13 G (see Eq. (17), with β_{out} = 10^{4}), so, initially, B_{z} ≈ 1.25 mG at R = 42 au in our simulations, which is of the same order of magnitude as the upper limit of B_{z}(R = 42 au) = 0.8 mG found in Vlemmings et al. (2019), for example. Hence, while the cavity is strongly magnetised, its field strength is compatible with observational constraints.
Compared to the outer disc, the mass and angularmomentum transport coefficients in the cavity are all of the order of unity, resulting in transsonic accretion velocities and faster wind with large lever arms (λ ≳ 10). Overall, this picture quantitatively matches the inner jetemitting disc proposed by Combet & Ferreira (2008). Interestingly, in all of our models, the cavity manages to reach an accretion rate close to the outer disc one by selfregulating the magnetic stresses. We find that most of the angularmomentum transport is due to the laminar stress (Appendix D) indicating that turbulent transport (possibly MRIdriven) is unimportant in the cavity. This is not surprising since our discs are dominated by ambipolar diffusion, which mostly suppresses MRI turbulence (Bai 2011).
We find a significant deviation of the rotation profile in the cavity as a result of the strong magnetic stress due to the wind and typical rotation velocities of the order of 70–80% of the Keplerian velocity. This fact, combined with the transsonic accretion, implies that the kinematics of these cavities have singular observational signatures. Fast accretion kinematics have been observed in some TDs (Rosenfeld et al. 2014), but we note that these signatures might also be due to a warped circumbinary disc (Casassus et al. 2015).
As a result of the stress balance mentioned above, we obtain accreting cavities that survive thousands of orbits and that are slowly expanding or contracting, depending on the outerdisc magnetisation. This result suggests that a cavity could be carved spontaneously if the magnetisation of the outer disc is high enough. There are already hints of such a process in global simulations; for instance, Cui & Bai (2021) show a gasdepleted cavity forming in the inner profile of Σ (see their Fig. 5, first row and first column panel). While this is by no means proof, since the boundary conditions are probably unrealistic, it shows that the secular evolution of winddriven discs should be investigated systematically to check whether or not cavities could spontaneously form in these models.
The temporal analysis of the disc reveals the appearance of dynamical structures. In particular, we highlight the formation of gas filaments above the disc surface that end up forming two bubbles of gas each, one being ejected, while the other one falls down onto the cavity before crossing it. At some point, the falling matter has to cross the poloidal magnetic field lines at the magnetic field loop location, resembling, to some extent, the magnetospheric accretion observed in young stars (Bouvier et al. 2007, 2020b,a; Pouilly et al. 2020) and magnetospheric ejection events (Zanni & Ferreira 2013; Čemeljić et al. 2013). However, there is no magnetosphere in our simulations, so the magnetic topology is quite different from that of magnetospheric interaction.
By analogy with magnetospheric accretion, we checked whether the time variability seen in our simulations could be due to a magnetic RTI. We have studied two criteria for the RTI, in the form of a radial interchange of poloidal field lines (see Sect. 3.5.4 and Appendix C). We found, however, that the RTI requires magnetisations stronger than the ones found in our simulations, ruling out the RTI in the form we have assumed. It is, however, still possible that another branch of this instability is present. It is also possible that the nonaxisymmetric version of the RTI could be triggered in 3D simulations. We therefore defer this study to a future publication.
On longer timescales, averaging out the fast variability, the magnetic field strength appears to be selfregulated with 0.1 ≲ ß ≲ 1 in the cavity, independently of the initial field strength.
As a result, the cavity is strongly magnetised and rotates at subKeplerian velocities, indicating a substantial magnetic support against gravity in this region. In essence, the regime of our cavity is similar to the magnetically arrested disc (MAD) proposed by Narayan et al. (2003) in the context of black hole accretion discs. McKinney et al. (2012) showed that MADs could be regulated by magnetic RTI, leading to magnetically chocked accretion flows (MCAF). The MAD model is also associated with the formation of plasmoids by reconnection events (Ripperda et al. 2022). These features are recovered in our models of TDs, despite the fact that we used Newtonian dynamics (MADs are usually found in GRMHD simulations) and the presence of a strong ambipolar diffusivity in our models. Hence, our models could be interpreted as nonideal, nonrelativistic models of MADs.
The time variability of the cavity is likely to be related to the axisymmetric approximation used in this work since it suppresses nonaxisymmetric instabilities, which seem to play a key role in MADs simulations (e.g. McKinney et al. 2012; Liska et al. 2022). Additionally, we note that the question of nonaxisymmetric hydrodynamical instabilities such as the Rossby wave instability (RWI) (Lovelace et al. 1999; Li et al. 2000) at the cavity edge is still open to debate in a magnetised environment (Bajer & Mizerski 2013). We will address these points using full 3D simulations in a followup paper.
Regarding the caveats of our simulations, we remark that the inner radial boundary is probably the most stringent caveat of our numerical model. In particular, we found that this inner boundary condition sometimes expels some poloidal magnetic flux, resulting in the bursts seen in Fig. 21. However, the weakly magnetised simulations (such as B5Bin0Am0) are the only ones exhibiting these events, and once the transient state is over, all the simulations reach comparable steady states. So, the inner boundary condition is likely not affecting the longterm evolution of our models. Future models should nevertheless try to include either an inner turbulent disc, or possibly the magnetospheric interaction with the central star.
A possible limitation of our model concerns the role of the MRI. Our simulated discs are dominated by ambipolar diffusion, and as such, subject to MRI quenching by the nonlinearity embedded in the ambipolar diffusivity (η ∝ B^{2}). This saturation is different from the saturation by 3D turbulence observed in the ideal MHD regime. It is suggested that the MRI saturates in very similar ways in 3D and 2D under strong ambipolar diffusion (see e.g. Béthune et al. 2017; Cui & Bai 2021). This is also confirmed by our own 3D simulations, which will be published in a forthcoming paper. Hence, the fact that our simulations are 2.5D has a very limited impact on the turbulent transport one may observe.
We note that our simulations used a simplified treatment of thermodynamics and ionisation chemistry. More numerically involved models, such as that of Wang & Goodman (2017), use a refined computation of the ionisation fraction and Λ_{Α} inside the cavity of a TD, including several chemical species. This work highlights, in particular, the influence of the Xray luminosity of the star L_{X} (see their Fig. 2, panels 2 and 3) as well as the role of the temperature T_{0} at 16 au (Fig. 2, panels 6 and 7). Regarding our profile of Λ_{Α} ≈ 1–10, our work is similar to their models 2 (with L_{X} = 10^{29} erg s^{−1}) and 6 (where T_{0} = 30 K). Therefore, we anticipate that an increase of two orders of magnitude for L_{X} would lead to Λ_{Α} > 10^{2} in most of our cavity. Such a change would greatly alter the dynamical regime of the cavity since MRI would then play a significant role see Appendix A and Blaes & Balbus (1994) and Bai (2011). However, the role of the temperature is less straightforward and seems to have little impact on Λ_{Α}.
Additionally, dust plays a significant role in the work of Wang & Goodman (2017) regarding the ionisation of the disc. As a matter of fact, only their models with dust reach low values of The effect of dust in TDs is a major subject that is not addressed in our work. Dust can modify the ionisation fraction but also create peculiar structures at the interface between the disc and cavity. We mention in particular the interplay between dust and the radiation pressure, which is known to create nonaxisymmetric structures at the cavity edge (Bi & Fung 2022) or an inner rim with an accumulation of matter due to photophoresis (Cuello et al. 2016).
Fig. 25 Spatiotemporal diagrams for 〈Σ〉 and 〈β〉 for R20FID. 
6 Conclusions
We performed 2.5D global numerical simulations of TDs in the context of nonideal MHD with MHD wind launching. Our simulation design is initialised with a cavity in the gas surface density profile, and a powerlaw distribution for the vertical magnetic field strength, resulting in a strongly magnetised cavity surrounded by a standard weakly magnetised disc.
The main results are summarised in the following points:
We modelled strongly accreting TDs that reach a quasisteady state that lasts for at least thousands of years. The accretion rate inside the cavity connects smoothly to the accretion rate in the external part of the disc;
The cavity itself is characterised by a strong subKeplerian rotation and a transsonic accretion velocity. These kinematic signatures could potentially be verified observationally;
The magnetic field is advected inwards in the outer disc, in contrast to full disc simulations. This points to the possible nonlocality of largescale field transport;
The cavity structure (density and field strength) is selfregulated. In particular, it is insensitive to a change in the initial internal magnetisation and is characterised by 0.1 ≲ β_{nt} ≲ 1;
The temporal analysis of the cavity dynamics highlights the formation and accretion of bubbles of gas above the disc which cross the cavity at sonic speeds. The magnetic RayleighTaylor instability might be responsible for this unsteadiness;
The physics of the cavity (accretion speed, wind lever arm, and mass loading) match previously published jetemitting disc solutions (Ferreira 1997; Combet & Ferreira 2008). The presence of a strong radial magnetic support and possible regulation by the RTI is also reminiscent of MADs in black hole physics (Narayan et al. 2003; McKinney et al. 2012). These resemblances suggest that TDs could be an instance of MADs applied to protoplanetary discs. Transition discs with strong accretion rates and arbitrarily large cavities can be achieved by magnetic winds emitted from the cavity. This model is promising and should be tested with observations.
Acknowledgements
The authors would like to thank the anonymous referee for constructive comments that have greatly improved the quality of this work. They wish to thank Jonatan JacqueminIde, Andres Carmona, Antoine Riols, Ileyk El Mellah and Jonathan Ferreira for fruitful discussions and comments. This work is supported by the European Research Council (ERC) European Union Horizon 2020 research and innovation programme (Grant agreement no. 815559 (MHDiscs)). This work was granted access to the HPC resources of TGCC under the allocation 2021A0100402231 made by GENCI. A part of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univgrenoblealpes.fr), which is supported by Grenoble research communities. This work makes use of matplotlib (Hunter 2007) for graphics, NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020) and Pickle (Van Rossum 2020). This article has been typeset from a T_{E}X/ LAT_{E}Xfile prepared by the authors.
Appendix A Ambipolar diffusivity for a transition disc  a simple model
The aim of this appendix is to model the ambipolar diffusivity spatial dependence in both a TD and a standard protoplanetary disc (i.e. without cavity). The general procedure to reach such a result follows and adapts the main calculation steps that are presented in Combet et al. (2010). As assumed in Eq. 4, only the ambipolar diffusivity does appear in the MHD equations, which we assume is the dominant nonideal effect in the regime of discs we used at R ≥ 10 (Riols et al. 2020; Simon et al. 2015). Therefore, the only momentum exchange that occurs between particles happens exclusively between ions and neutrals. In a plasma made of molecular ions, electrons and neutrals, the ambipolar diffusivity is given by (A.1)
(Wardle 2007), where ρ_{n} and ρ_{i} are, respectively, the density of the neutrals (the gas so ρ_{n} = ρ) and of the ions and γ_{in} = 〈σv〉_{in}/(m_{n} + m_{i}) with 〈σv〉_{in} is the ionneutral collision rate whose value is (A.2)
(Bai 2011), with m_{H} being the atomic mass and μ = 2,34 m_{H} the mean molecular weight. Introducing the ionisation fraction ξ = ρ_{i}/ρ_{n}, one obtains (A.3)
Ambipolar diffusion is usually evaluated with the dimensionless ambipolar Elsasser number Λ_{A} defined in Eq. 13. To obtain this number, we have to evaluate the ionisation fraction. Let us consider a simple chemical lattice with no metals nor grains: (A.4) (A.5)
with ζ_{i} being the ionisation rate and δ the dissociative recombination rate. Following Fromang et al. (2002), we take (A.6)
In this toy model, we then have (A.7)
(Lesur et al. 2014), where ξ_{FUV} accounts for the farUV photon contribution that we modelled following PerezBecker & Chiang (2011) as (A.8)
with Σ_{*} the column density computed from the star to the point of interest.
To calculate ζ_{i}, we considered several ionisation sources. We modelled the Xray ionisation from the protostar by two bremsstrahlungemitting corona (following Bai & Goodman 2009 and Igea & Glassgold 1999): (A.9)
with L_{x,29} = L_{x}/10^{29}erg s^{−1} and L_{x}, ζ_{1}, ζ_{2}, α, ß, N_{1}, N_{2} being the numerical values defined in Bai & Goodman (2009), while N_{H1} and N_{H2} are the column densities of hydrogen vertically computed above and below the calculation point.
We modelled the cosmicray ionisation rate following Umebayashi & Nakano (1980) so that: (A.10)
where ζ_{CR,0} = 10^{−17} s^{−1} and is the matter column density above and below the point of interest.
Fig. A.1 Ambipolar Elsasser number Λ_{A} in a standard protoplanetary disc (top panel) and in a TD (bottom panel). In spite of these 2 profiles being slightly different, no major changes occur from one kind of disc to another around the midplane. 
Lastly, we added the radioactive decay that is assumed constant (Umebayashi & Nakano 2008) with a ionisation rate given by (A.11)
Combining the Eqs. A.9, A.10, and A.11, we obtain ζ_{i} = ζ_{X} + ζ_{CR} + ζ_{rad}., paving the way to finally reach Λ_{A} using Eqs. 13, A.6, and A.7. We note that due to the dependency of η_{A} and v_{A} on the norm of the magnetic field, this latter is cancelled out and does not need to be computed to obtain ΛA. The previous calculations can be performed either for a standard protoplanetary disc or for a TD. The only thing that needs to be changed to account for such discs is the surface density profile, where Eq. 18 allows us to consider (or not) the effects of the cavity.
The results of such calculations are displayed in Fig. A.1, which represents the spatial dependency of AA in both a standard protoplanetary disc and a TD. Though these two profiles look different at first glance, a deeper investigation reveals that the values taken by AA in the discs remain quite close to unity in both cases, while the general trend of Λ_{A} in a standard protoplanetary disc is recovered even in the case of a TD (Thi et al. 2019).
Moreover, Λ_{A} remains fairly below the critical value Λ_{a,crit}. = 10^{2} with or without a cavity. Λ_{A} must stay below Λ_{a,crit}. so that the MRI effects are negligible (Blaes & Balbus 1994; Bai 2011). Therefore, assuming a characteristic value of Λ_{A},_{0} = 1 captures, with a reasonable level of accuracy, the physics of ambipolar diffusion, and the cavity does not alter the ambipolar Elsasser number profile. The results we achieve from this simple toy model are to be compared to the more detailed work of Wang & Goodman (2017), in which many chemical species were taken into account to compute the ambipolar Elsasser number inside the cavity of a winddriven TD.
Following Lesur (2021b) and Thi et al. (2019), we implemented the profile of Λ_{A} so that (A.12)
where λ is a parameter that controls the height where a transition between nonideal and ideal MHD occurs (the nonideal MHD part being the inside of the disc) and is chosen as constant and equal to 3 h. Λ_{a,0} remains a free parameter (see Section 1 for more details). Additionally, a cutoff is used for the η_{Α} profile so that if η_{Α} > η_{Α,max}, the value of η_{Α} is replaced by η_{Α,max} = 10 ε^{2} in code units, such a choice being reflected in the Λ_{A} profile with Eq. 13.
Appendix B Poloidal velocity relaxation and inner boundary condition
We aim to address the influence of the poloidal velocity relaxation on our results to test our control on the inner boundary condition. Two additional simulations were conducted with the same setup as B4Bin0Am0 (fiducial run) and B5Bin0Am0, respectively, but without the relaxation procedure. The results are given in Fig. B.1, where we show the surface density 〈Σ〉_{4000} time averaged on the first 4000 orbits at the internal radius (when the differences are enhanced), with a focus on the innermost radii. We highlight that these differences do not increase for t > 4000 orbits at R_{int}. For B5Bin0Am0, the right panel of Fig. B.1 suggests that the relaxation procedure influences how the initial burst is evacuated since we detect differences between the surfacedensity profiles at R > 1.5. However, releasing this inner constraint reduces the inner peak of the profile of Σ, but it does not prevent the initial accumulation of matter from appearing. In particular, the bursts of matter seen in Fig. 21 are not due to this condition (and are probably due to the inner boundary condition; see next paragraph). For the fiducial simulation, we estimate differences of 15% until R = 2, 7% until R = 10, and less than 2% until R = 50, and we conclude that the slight accumulation described in Section 3.1 is due to this procedure, contrary to the occurrence of bursts as seen in Fig. 21.
Regarding the bursts of B5Bin0Am0 (see Fig. 21), we focus on one of them in Fig. B.2. The first panel displays the spatiotemporal diagram of the surface density on which the burst is clearly detected at 17435 orbits at R_{int} and localised by the red dashed line. The accumulation of matter is correlated with a decrease of the vertical magnetic field at the midplane (second panel of Fig. B.2). This magnetic field is not lost but is expelled outwards, as is evident from the magnetic flux function (third panels of Fig. B.2). Such a shortage of magnetic field leads to an increase of β and blocks accretion (we recall that the accretion speed is v_{acc}. ∝ β^{−σ} with σ > 0). As a result, Ṁ falls from 0,25 down to 0,1 × 10^{−7} M_{⊙} yr^{−1} in the region between the inner radial boundary and the burst, and matter piles up in the cavity. This episode ends when the magnetic flux is eventually reaccreted, leading to an increase of the massaccretion rate and the disappearance of the density excess in the cavity. At some point, the magnetic flux is advected back onto the seed until it saturates so that B_{z} can accumulate again close to the inner boundary condition before accretion is enhanced back to normal. The reason why such magnetic flux evacuates from the seed from time to time remains unclear, and these occurrences close to the inner boundary suggest that these might be a boundary condition artefact. However, we mention that the total magnetisation of the seed eventually saturates with a roughly constant value, so a sharp increase of magnetic field (as is the case for this burst; see the middle panel of Fig. B.2, a few orbits before the location of the red dashed line) could force the seed to lose magnetic flux to ensure its conservation. We finish by adding that these bursts are only detected for the weakly magnetised simulations (the ones with ß_{out} = 10^{5}).
Regarding the inner radial boundary condition for the magnetic field, we tried several configurations (an outflow condition, which is the one we eventually chose, and a perfect conductor). Both of these conditions lead to the same steady states.
We also ran a simulation with a stronger magnetic field close to the inner boundary condition, but no significant changes were noticed. The additional magnetic field was chosen so that the magnetisation of the seed is set close to its saturation value in the fiducial run. However, in any case, the same transient state occurs and leaves the stage to a similar steady state (the magnetisation of the seed reaches the same saturation value and the same stripes are observed in the spatiotemporal diagram of ψ).
We therefore conclude that our setup is robust regarding the initial state and the boundary conditions. The inner boundary still plays a role because of its magnetisation and the fact that only a given amount of magnetic field can be advected. This probably leads to the burst events seen in simulation B5Bin0Am0.
Fig. B.1 Surface density timeaveraged on the first 4000 orbits at R_{int}. The blue lines are the reference runs (left panel: fiducial run, right panel: B5Bin0Am0), and the reddashed mines are the corresponding runs without the relaxation. 
Fig. B.2 Spatiotemporal diagrams of Σ (first panel), B_{z,0} the vertical magnetic field at the midplane and ψ the flux function defined in Eq. 39, for simulation B5Bin0Am0. These profiles focus on the second burst detected in the left panels of Fig. 21. The red dashed line marks the beginning of the burst when detected using Σ. 
Appendix C Interchange instability criterion calculations
We express the instability criterion for the interchange instability (or RTI) calculated in Spruit et al. (1995) (Eq. 59) in terms of the plasma parameter. This criterion reads (C.1)
where S is the shear that we approximate with S^{2} = 9/4 Ω^{2}, and g_{m} is (C.2)
is the radial component of the magnetic field at the disc surface. Let us rewrite the previous expression in terms of ß, q (defined with ) and δ (defined as δ = −d ln Σ/d ln R): (C.3) (C.4)
where X' denotes the derivative of X with respect to R. With , we obtain (C.5)
Therefore, the instability criterion becomes (C.6)
ε being constant in the disc as well as β inside the cavity. Ω_{κ} varies as R^{−3/2} and Σ as R^{−δ} so that (C.7)
By taking , the RTI can be triggered when (C.8)
If we now assume that δ = q = 1 for simplicity, we finally obtain (C.9)
where ε = 0.1.
Fig. C.1 Interchange instability criteria. The red dotted line shows the critical value of β_{crit}, while the green solid line is obtained with Eq. C.8. 
Figure C.1 compares the timeaveraged values of β with the criterion given in Eq. C.9. The value of ß_{crit} is below the timeaveraged values of . Though this simple analysis makes it difficult to make definitive conclusions on this subject, it seems that the interchange instability is not triggered inside the cavity.
Appendix D Laminar transport coefficients
In order to discuss the role of the MRI, we must highlight the impact of the laminar stress and its contribution to the transport coefficients. In this article, we focus on the total stresses, defined in Eq. 29. To compare the turbulent effects, we decompose the stresses with a turbulent and a laminar part. In this prospect, we introduce the deviation to the temporal mean such that (D.1)
Focusing on W_{rϕ}, we expand the magnetic term as (D.2)
Concerning the turbulent stresses, we referred to JacqueminIde et al. (2021) (see their Appendix A) as we only computed the laminar ones and compared the laminar transport coefficients to the ones studied in the article. Therefore, we adopted the following definition for the laminar radial stress: (D.3)
and for the laminar surface stress (D.4)
These definitions are coherent with previous works (Béthune et al. 2017; Mishra et al. 2020; JacqueminIde et al. 2021). Hence, the laminar transport coefficients are given by
while we define their turbulent counterparts as
The results are shown in Fig. D.1. The laminar contribution is the major one for 〈u_{w}〉 in the whole disc, so we only show its laminar contribution with respect to the full coefficient, as they take essentially the same values Nevertheless, despite the laminar term being high for 〈α〉, a strong turbulent term is at stake, especially in the external part of the disc where it is dominant Inside the cavity, 〈α〉 is fairly distributed between the laminar and turbulent contributions However, we recall that the wind may act also on the turbulent component of α) since the magnetic field also appears in Eq. D.2. We finally conclude that the MRI is probably acting on the disc outer parts in the 〈α〉 coefficient, while the surface stress embodied by 〈U_{W}〉 is definitely dominated by the laminar part and due to the wind.
Fig. D.1 Timeaveraged transport coefficients and their laminar and turbulent contributions. We give the laminar and turbulent contributions for (α) and the total profile with its laminar contribution for 〈U_{W}〉. 
References
 Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475 [Google Scholar]
 Bai, X.N. 2011, ApJ, 739, 50 [CrossRef] [Google Scholar]
 Bai, X.N., & Goodman, J. 2009, ApJ, 701, 737 [NASA ADS] [CrossRef] [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]
 Bajer, K., & Mizerski, K. 2013, Phys. Rev. Lett., 110, 104503 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bi, J., & Fung, J. 2022, ApJ, 928, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Bouvier, J., Alencar, S. H. P., Boutelier, T., et al. 2007, A&A, 463, 1017 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bouvier, J., Alecian, E., Alencar, S. H. P., et al. 2020a, A&A, 643, A99 [EDP Sciences] [Google Scholar]
 Bouvier, J., Perraut, K., Bouquin, J.B. L., et al. 2020b, A&A, 636, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carmona, A., Pinte, C., Thi, W. F., et al. 2014, A&A, 567, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carmona, A., Thi, W. F., Kamp, I., et al. 2017, A&A, 598, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Casassus, S., Marino, S., Pérez, S., et al. 2015, ApJ, 811, 92 [Google Scholar]
 Čemeljić, M., Shang, H., & Chiang, T.Y. 2013, ApJ, 768, 5 [CrossRef] [Google Scholar]
 Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Combet, C., & Ferreira, J. 2008, A&A, 479, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Combet, C., Ferreira, J., & Casse, F. 2010, A&A, 519, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cuello, N., Gonzalez, J.F., & Pignatale, F. C. 2016, MNRAS, 458, 2140 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, C., & Bai, X.N. 2021, MNRAS, 507, 1106 [NASA ADS] [CrossRef] [Google Scholar]
 Dong, R., & Dawson, R. 2016, ApJ, 825, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [CrossRef] [EDP Sciences] [Google Scholar]
 Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 497 [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Fang, M., Kim, J. S., van Boekel, R., et al. 2013, ApJS, 207, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Ferreira, J. 1997, A&A, 319, 340 [Google Scholar]
 Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [Google Scholar]
 Gárate, M., Delage, T. N., Stadler, J., et al. 2021, A&A, 655, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [Google Scholar]
 Guilet, J., & Ogilvie, G. 2014, MNRAS, 441, 852 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848 [Google Scholar]
 JacqueminIde, J., Lesur, G., & Ferreira, J. 2021, A&A, 647, A192 [EDP Sciences] [Google Scholar]
 Kane, Y., 1966, IEEE Trans. Antennas Propag., 14, 302 [CrossRef] [Google Scholar]
 Lesur, G. R. J. 2021a, J. Plasma Phys., 87, 205870101 [CrossRef] [Google Scholar]
 Lesur, G. R. J. 2021b, A&A, 650, A35 [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]
 Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] [Google Scholar]
 Liska, M. T. P., Musoke, G., Tchekhovskoy, A., Porth, O., & Beloborodov, A. M. 2022, ApJl, 935, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marsh, K. A., & Mahoney, M. J. 1992, ApJ, 395, L115 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mishra, B., Begelman, M. C., Armitage, P. J., & Simon, J. B. 2020, MNRAS, 492, 1855 [NASA ADS] [CrossRef] [Google Scholar]
 Morishima, R. 2012, MNRAS, 420, 2851 [NASA ADS] [CrossRef] [Google Scholar]
 Najita, J. R., Strom, S. E., & Muzerolle, J. 2007, MNRAS, 378, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69 [NASA ADS] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
 PerezBecker, D., & Chiang, E. 2011, ApJ, 735, 8 [Google Scholar]
 Pouilly, K., Bouvier, J., Alecian, E., et al. 2020, A&A, 642, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riols, A., & Lesur, G. 2018, A&A, 617, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riols, A., & Lesur, G. 2019, A&A, 625, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJ, 924, L32 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenfeld, K. A., Chiang, E., & Andrews, S. M. 2014, ApJ, 782, 62 [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Simon, J. B., Bai, X.N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013, ApJ, 764, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C., & Taam, R. E. 1990, A&A, 229, 475 [NASA ADS] [Google Scholar]
 Spruit, H. C., Stehle, R., & Papaloizou, J. C. B. 1995, MNRAS, 275, 1223 [NASA ADS] [CrossRef] [Google Scholar]
 Stehle, R., & Spruit, H. C. 2001, MNRAS, 323, 587 [NASA ADS] [CrossRef] [Google Scholar]
 Suriano, S. S., Li, Z.Y., Krasnopolsky, R., Suzuki, T. K., & Shang, H. 2019, MNRAS, 484, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thi, W. F., Lesur, G., Woitke, P., et al. 2019, A&A, 632, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Umebayashi, T., & Nakano, T. 1980, PASJ, 32, 405 [NASA ADS] [Google Scholar]
 Umebayashi, T., & Nakano, T. 2008, ApJ, 690, 69 [Google Scholar]
 van der Marel, N., van Dishoeck, E. F., Bruderer, S., Pérez, L., & Isella, A. 2015, A&A, 579, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Rossum, G. 2020, The Python Library Reference, release 3.8.2 (Python Software Foundation) [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 [CrossRef] [Google Scholar]
 Vlemmings, W. H. T., Lankhaar, B., Cazzoletti, P., et al. 2019, A&A, 624, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, L., & Goodman, J. J. 2017, ApJ, 835, 59 [Google Scholar]
 Wardle, M. 2007, Astrophys. Space Sci., 311, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218 [Google Scholar]
 Zanni, C., & Ferreira, J. 2013, A&A, 550, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhang, K., Isella, A., Carpenter, J. M., & Blake, G. A. 2014, ApJ, 791, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34 [Google Scholar]
 Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
MHD invariants for a subset of simulations, computed with timeaveraged quantities on the last 1000 orbits at R_{int}.
All Figures
Fig. 1 Initial and timeaveraged profile of Σ (top left panel), (top right), B_{Z,0} (bottom left), and ν_{φ,ο} (bottom right), with respect to r. 

In the text 
Fig. 2 Schematic view of the disc represented in orange. θ_{±} define the vertical integration surface and h_{int} is the integration scale height at a given radius r. 

In the text 
Fig. 3 Surface density (left panel) and plasma beta parameter β (right panel) as a function of R at midplane and time for the fiducial simulation. The cavity remains during the entire simulation and keeps a relatively strong magnetisation with . 

In the text 
Fig. 4 Radial profiles of the surface density and of the vertical magnetic field in the midplane, timeaveraged on the last 1000 orbits at R_{int}. The vertical magnetic field is vertically averaged, and both profiles are given in arbitrary units. 

In the text 
Fig. 5 Timeaveraged poloidal magnetic field lines and toroidal field component ⟨Β_{φ}⟩ for the fiducial simulation. We note the peculiar field topology close to the truncation radius. 

In the text 
Fig. 6 Timeaveraged streamlines and density for the fiducial simulation. We note the peculiar shape of the streamlines around the transition radius. 

In the text 
Fig. 7 Timeaveraged angularmomentumflux streamlines over timeaveraged density for the fiducial simulation. Angular momentum leaves the disc midplane because of the wind. 

In the text 
Fig. 8 Accretion rates for different integration height scales with respect to the radius inside the disc. The higher ρv_{r} is integrated the closer to a constant value Ṁ is in the cavity. The average value inside the cavity (from R = 1 to R = 10} is Ṁ = 1.4 ± 0.2 × 10^{−7} M_{⊙} yr^{−1}. 

In the text 
Fig. 9 Accretion speed for ɛ_{int} = 0.9 in units of local sound speed c_{s}. The profile exhibits a clear transition between subsonic and transsonic accretion that occurs where the edge of the cavity is located. 

In the text 
Fig. 10 Mass conservation for ɛ_{int} = 0.9. The time derivative of Σ remains perfectly constant and equal to zero inside the cavity and only becomes moderately negative at the cavity edge R = 10. This suggests that the simulation indeed reaches a steady state for radii up to ≈8. The mass conservation is also correctly recovered. The three lines do not add up to zero because we used a moving average for better visibility, and the quantities are timeaveraged on a sample selection of output flies that do not contain all the time steps computed by the code. 

In the text 
Fig. 11 ⟨ζ⟩ parameter for ɛ_{int} = 0.6. The selfsimilar fit shown here is for comparison only and was also obtained for ɛ_{int} = 0.6. It appears that it is coherent for the external disc, while it predicts a wind way too intense in the internal part, therefore, another model is used to describe ⟨ζ⟩ in the whole disc. 

In the text 
Fig. 12 Angular momentum conservation multiplied by r^{−3,2} and time averaged. Full blue line is , red dotdashed line is , green dashed line is , and purple dotted line is . 

In the text 
Fig. 13 Timeaveraged transport coefficients ⟨α⟩ and ⟨v_{w}⟩ for ɛint = 0.9. 

In the text 
Fig. 14 General structure of the wind and MHD invariants. First panel shows the field lines in the internal and external disc. The grey dashed lines represent the surfaces at h/R = 0.3; 0.6 and 0.9. The three other plots display the MHD invariants for the internal field line (full, light blue) and the external one (semidashed, dark blue). The invariants are time averaged on the last 1000 orbits. 

In the text 
Fig. 15 Flux function ψ(R, t) for the magnetic field, taking into account the flux at the surface of the seed and the radial flux. For radii larger than the one of the cavity, the field lines are advected toward the centre during the whole simulation. At R ≈ 13, the flux accumulates and exhibits a striped structure for smaller radii. 

In the text 
Fig. 16 Magneticfield transport parameter v_{B} as a function of radius. We note that the outer disc is transporting magnetic field lines inwards (v_{B} < 0). 

In the text 
Fig. 17 Temporal evolution of Σ in dotted green, Ṁ in dashed blue, and B_{z} (vertically averaged) in black (full line) at R = 3 au for the fiducial simulation. ζ is calculated at ɛ_{int} = 0.3 and shown by a semidashed red line with a logarithmic scale on the right of the panels. Apart from ζ, all the profiles are given in arbitrary units and divided by their maximum value reached during the timescale of the top panel. 

In the text 
Fig. 18 Density and magnetic field lines at different times showing the advection of a bubble of material (arrow) from the outer disc through the cavity. 

In the text 
Fig. 19 Σ(R, t) and of B4Bin0Am1. The cavity stands throughout the simulation, though its edge falls down to R ≈ 2 au during the transient, before broadening up to R ≈ 4 au in a few thousands of orbits at R_{int}. The profile of is characterised by the presence of gaps in the external disc. 

In the text 
Fig. 20 Timeaveraged structure of the disc for B4Bin0Am1. Left panel: poloidal streamlines and density. Right panel: magnetic structure of the disc with magnetic poloidal field lines and 〈Β_{ϕ}〉. 

In the text 
Fig. 21 Spatiotemporal diagrams for Σ and for B5Bin0Am0 and B3Bin0Am0. The cavity expands more in B3Bin0Am0 but shrinks in B5Bin0Am0. Bursts of matter occur in B5Bin0Am0 (see main text for more explanations). 

In the text 
Fig. 22 Impact of initial external magnetisation on the surface density. We average the profile on the last 1000 orbrts at R_{int}. For B5Bin0Am0, we average on 1000 orbhs at R_{int} occurring between the 2 burst events seen in Fig. 21. 

In the text 
Fig. 23 Impact of internal initial magnetisation on the plasma parameter for β_{out} = 10^{4}. 

In the text 
Fig. 24 Spatiotemporal diagrams of Σ (first panel), B_{z,0}, the vertical magnetic field at the midplane, and ψ, the flux function defined in Eq. (39). These profiles focus on the first orbits of the run B4Bin3Am0. 

In the text 
Fig. 25 Spatiotemporal diagrams for 〈Σ〉 and 〈β〉 for R20FID. 

In the text 
Fig. A.1 Ambipolar Elsasser number Λ_{A} in a standard protoplanetary disc (top panel) and in a TD (bottom panel). In spite of these 2 profiles being slightly different, no major changes occur from one kind of disc to another around the midplane. 

In the text 
Fig. B.1 Surface density timeaveraged on the first 4000 orbits at R_{int}. The blue lines are the reference runs (left panel: fiducial run, right panel: B5Bin0Am0), and the reddashed mines are the corresponding runs without the relaxation. 

In the text 
Fig. B.2 Spatiotemporal diagrams of Σ (first panel), B_{z,0} the vertical magnetic field at the midplane and ψ the flux function defined in Eq. 39, for simulation B5Bin0Am0. These profiles focus on the second burst detected in the left panels of Fig. 21. The red dashed line marks the beginning of the burst when detected using Σ. 

In the text 
Fig. C.1 Interchange instability criteria. The red dotted line shows the critical value of β_{crit}, while the green solid line is obtained with Eq. C.8. 

In the text 
Fig. D.1 Timeaveraged transport coefficients and their laminar and turbulent contributions. We give the laminar and turbulent contributions for (α) and the total profile with its laminar contribution for 〈U_{W}〉. 

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.