© É. Martel and G. Lesur 2022

1 Introduction

Transition discs (TDs) are protoplanetary discs exhibiting a deficit of near-infrared 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 disc-less young stellar objects, hence their name. In this framework, TDs are the result of an inside-out dispersal process, which is usually believed to be a combination of viscous accretion, dust growth (Dullemond & Dominik 2005), giant planets (Marsh & Mahoney 1992), and photo-evaporation (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 ro-vibrational 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⁻⁹ M_⊙ yr⁻¹, it also predicts a large fraction of non-accreting 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 angular-momentum transport in the cavity. The most studied candidate for this is planet-disc 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 small-scale turbulence, which could be of hydrodynamic (vertical shear instability or VSI, Nelson et al. 2013) or magnetic (magneto-rotational 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) wind-driven models is mostly controlled by the strength of the large-scale 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 fast-accreting TD cavity.

The idea of having a magnetic wind-driven 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 self-consistent (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 wind-emitting outer disc, subject to realistic magnetic diffusion, and that the resulting configuration can be long-lived. Given the richness of the dynamics, we first concentrated on 2.5D models in this work, and we will discuss 3D models in a follow-up paper. The paper is divided as follows. We first introduce the model equations, physical quantities, and numerical setup. We then discuss an in-depth 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 non-ideal MHD equations

In the following, we place ourselves in the non-relativistic, 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 non-ideal 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)

and (6)

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³)^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 angular-momentum conservation equations (e.g. the last term of Eq. (16) in Zhu & Stone 2018). This will simplify the interpretation of angular-momentum 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₀ 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₀ 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 Godunov-type scheme and a second-order Runge-Kutta 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) = Ω₀ (R/R_int)^−3/2=R^−3/2. We chose 300 g cm⁻² as a unit for the surface density and express the accretion rate in M_⊚ yr⁻¹. We denote code units as c.u. hereafter. We use the subscript X₀ 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 so-called outer disc refers to the region where the disc is full and described by a standard proto-planetary 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 mid-plane (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 wave-absorbing 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 ρ₀ 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 self-similar stationary disc solutions (Jacquemin-Ide 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)

with (19)

where Σ₀(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₀ 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 theta-averaged ‘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⟩₁₀₀₀ being the time-averaged value of X during the last 1000 orbits at Rint, for example.

Fig. 1 Initial and time-averaged profile of Σ (top left panel), (top right), BZ,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 hint 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 high-resolution (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⁴, β_in = 1, Λ_A = 1 and R₀ = 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 striped-like 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⁻² are reached from time to time. After approximately 4000 orbits at the internal radius, both Σ and reach a quasi-stationary state.

Gaps and rings are detected in the outer part of the disc, in the spatio-temporal 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 = 12-18 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.

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, time-averaged on the last 1000 orbits at Rint. 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 time-averaged 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 elbow-shaped 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 time-averaged 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 elbow-like 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 Angular-momentum-flux streamlines

In order to deepen the analysis of the role of the magnetic structure, we concentrated on the time-averaged angular-momentum flux, defined by (24)

The poloidal flux lines associated with this angular-momentum 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 elbow-like shape for the angular-momentum 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 Time-averaged 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 elbow-shaped 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 mid-plane. This radius corresponds to the location of the basis of the elbow-shaped 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 well-defined transition between subsonic accretion outside the cavity with ⟨v_acc.⟩ ~ 10⁻³ ⟨c_s⟩ and transsonic accretion inside with ⟨v_acc.⟩ ~ ⟨c_s⟩.

Fig. 6 Time-averaged 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 angular-momentum 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 Time-averaged angular-momentum-flux streamlines over time-averaged 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 ρvr 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 cs. 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 mass-loss-rate parameter

The mass conservation equation is given by Eq. (27). Figure 10 shows the mass conservation for ɛ_int = 0.9 with time-averaged 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 ‘elbow-shaped structure’ in the poloidal streamlines).

In order to quantitatively account for the role of the wind, we constructed the mass-loss-rate 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 self-similar models (Lesur 2021b), we studied the values of ⟨ζ⟩ at zo = 6 h, which corresponds to ɛ_int = 0.6. The mass-loss-rate parameter is approximately constant in the external part of the disc around 6.2 × 10⁻⁵, while it peaks at 2.9 × 10⁻² 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 self-similar solutions, we show the self-similar scaling of the mass-loss-rate 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 self-similar.

The wind mass-loss-rate parameter is smaller than the self-similar scaling in the outer disc by a factor of a few. This discrepancy is probably due to the influence of the cavity mag-netosphere that compresses the disc magnetosphere, resulting in a deviation of ζ from the self-similar result. Moreover, it seems that the further we move outward, the closer we are to the self-similar values, indicating that we recover self-similar scalings far ‘enough’ from the cavity, as expected.

In the cavity, ζ is significantly weaker than expected from a naive extrapolation of self-similar scaling laws. This indicates that the mass-loss rate saturates at β ~ 1, a regime which was not explored by Lesur (2021b).

An alternative model to the self-similar one is used to describe ⟨ζ⟩ with greater accuracy. The self-similar 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 self-similar 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 time-averaged 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 self-similar 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 Angular-momentum conservation

In Fig. 12, we show the terms involved in the angular-momentum 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 mass-conservation 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⁻⁴, 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⁻⁴, 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 wind-emitting protoplanetary discs, with transport coefficients close to the ones found in self-similar wind models for β ~ 10⁴, 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.

Fig. 12 Angular momentum conservation multiplied by r−3,2 and time averaged. Full blue line is , red dot-dashed line is , green dashed line is , and purple dotted line is .

3.4 MHD wind

It is well known that steady-state 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 time-averaged quantities and listed in Table 3.

The mass-loading 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 large-scale 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 steady-state 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 sub-Keplerian in this region (Fig. 1). Quantitatively, we find ĸ_in = 2.2 × 10⁻², ĸ_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 super-Alfvé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 Time-averaged transport coefficients ⟨α⟩ and ⟨vw⟩ for ɛint = 0.9.

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 (semi-dashed, 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 space-time 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 cavity-edge expansion

As previously mentioned, the cavity edge moves slowly outwards during the simulation. Neglecting the impact of the wind in terms of mass-loss 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₀, 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 iso-contours of ψ describe the motion of the magnetic field lines in the disc plane. The spatio-temporal 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 large-scale 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⁻³ 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⁻³ 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⁻³ 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 Magnetic-field transport parameter vB as a function of radius. We note that the outer disc is transporting magnetic field lines inwards (vB < 0).

Fig. 17 Temporal evolution of Σ in dotted green, Ṁ in dashed blue, and Bz (vertically averaged) in black (full line) at R = 3 au for the fiducial simulation. ζ is calculated at ɛint = 0.3 and shown by a semi-dashed 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 time-averaged quantities and ignored fast variability. While our numerical solutions are quasi-steady 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 spatio-temporal 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 anti-parallel poloidal field lines meet at the elbow-shaped 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 mid-plane 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 self-similar models for which the ejection is continuous with a higher mass-loss 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 elbow-shaped structure on the time-averaged 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 Rint. 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 spatio-temporal 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 (Jacquemin-Ide et al. 2021) or non-ideal 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 time-averaged magnetic structure of the disc. The main features of the fiducial simulation are recovered, namely the elbow-shaped 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 ambipolar-dominated 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⁻⁷ M_⊙ yr⁻¹ 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 mass-loaded 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³ (run B3Bin0Am0) and β_out = 10⁵ (run B5Bin0Am0).

With regard to B5Bin0Am0, the spatio-temporal 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⁵) 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⁻⁴ and δ Σ ≈ 4.2 × 10⁻² (both in code units), so Eq. (38) gives c.u., where we chose R₀ ≈ 20. The simple model we used seems to overestimate the widening velocity of the cavity but still gives the correct order of magnitude.

The time-averaged 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₀Σ, 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 Time-averaged 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 Spatio-temporal 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 Rint. For B5Bin0Am0, we average on 1000 orbhs at Rint occurring between the 2 burst events seen in Fig. 21.

Fig. 23 Impact of internal initial magnetisation on the plasma parameter for βout = 104.

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⁴. 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³.

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 Spatio-temporal diagrams of Σ (first panel), Bz,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 double-sized cavity (R₀ = 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 elbow-shaped 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 steady-state 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, ambipolar-dominated protoplanetary disc (Lesur 2021b; Cui & Bai 2021). In particular, we find mass and angular-momentum 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 non-ideal 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 elbow-shaped magnetic surfaces at the transition radius. In any case, it points to the fact that magnetic field transport is a non-local 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.13 G (see Eq. (17), with β_out = 10⁴), 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 angular-momentum 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 jet-emitting 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 self-regulating the magnetic stresses. We find that most of the angular-momentum transport is due to the laminar stress (Appendix D) indicating that turbulent transport (possibly MRI-driven) 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 outer-disc 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 gas-depleted 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 wind-driven 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 non-axisymmetric 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 self-regulated with 0.1 ≲ ß ≲ 1 in the cavity, independently of the initial field strength.

As a result, the cavity is strongly magnetised and rotates at sub-Keplerian 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 non-ideal, non-relativistic 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 non-axisymmetric 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 non-axisymmetric 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 follow-up 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 long-term 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 non-linearity embedded in the ambipolar diffusivity (η ∝ B²). 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 X-ray luminosity of the star L_X (see their Fig. 2, panels 2 and 3) as well as the role of the temperature T₀ 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²⁹ erg s⁻¹) and 6 (where T₀ = 30 K). Therefore, we anticipate that an increase of two orders of magnitude for L_X would lead to Λ_Α > 10² 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 non-axisymmetric 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 Spatio-temporal diagrams for 〈Σ〉 and 〈β〉 for R20FID.

6 Conclusions

We performed 2.5D global numerical simulations of TDs in the context of non-ideal MHD with MHD wind launching. Our simulation design is initialised with a cavity in the gas surface density profile, and a power-law 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:

1. We modelled strongly accreting TDs that reach a quasi-steady 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;

2. The cavity itself is characterised by a strong sub-Keplerian rotation and a transsonic accretion velocity. These kinematic signatures could potentially be verified observationally;

3. The magnetic field is advected inwards in the outer disc, in contrast to full disc simulations. This points to the possible non-locality of large-scale field transport;

4. The cavity structure (density and field strength) is self-regulated. In particular, it is insensitive to a change in the initial internal magnetisation and is characterised by 0.1 ≲ β_nt ≲ 1;

5. 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 Rayleigh-Taylor instability might be responsible for this unsteadiness;

6. The physics of the cavity (accretion speed, wind lever arm, and mass loading) match previously published jet-emitting 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 Jacquemin-Ide, 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 2021-A0100402231 made by GENCI. A part of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.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_EX/ LAT_EXfile 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 protoplane-tary 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 non-ideal 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 ion-neutral 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 dimension-less 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 far-UV photon contribution that we modelled following Perez-Becker & 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 X-ray ionisation from the protostar by two bremsstrahlung-emitting corona (following Bai & Goodman 2009 and Igea & Glassgold 1999): (A.9)

with L_x,29 = L_x/10²⁹erg s⁻¹ and L_x, ζ₁, ζ₂, α, ß, N₁, N₂ 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 cosmic-ray ionisation rate following Umebayashi & Nakano (1980) so that: (A.10)

where ζ_CR,0 = 10⁻¹⁷ s⁻¹ 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 protoplan-etary 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² 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,₀ = 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 wind-driven 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 non-ideal and ideal MHD occurs (the non-ideal 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 cut-off is used for the η_Α profile so that if η_Α > η_Α,max, the value of η_Α is replaced by η_Α,max = 10 ε² 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 〈Σ〉₄₀₀₀ 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 surface-density 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⁻⁷ M_⊙ yr⁻¹ 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 re-accreted, leading to an increase of the mass-accretion 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⁵).

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 spatio-temporal 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 time-averaged on the first 4000 orbits at Rint. The blue lines are the reference runs (left panel: fiducial run, right panel: B5Bin0Am0), and the red-dashed mines are the corresponding runs without the relaxation.

Fig. B.2 Spatio-temporal diagrams of Σ (first panel), Bz,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² = 9/4 Ω², 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 time-averaged values of β with the criterion given in Eq. C.9. The value of ß_crit is below the time-averaged 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 Jacquemin-Ide 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; Jacquemin-Ide 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 Time-averaged 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 〈UW〉.

References

All Tables

All Figures

	Fig. 1 Initial and time-averaged profile of Σ (top left panel), (top right), BZ,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 hint 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, time-averaged on the last 1000 orbits at Rint. The vertical magnetic field is vertically averaged, and both profiles are given in arbitrary units.
In the text

	Fig. 5 Time-averaged 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 Time-averaged streamlines and density for the fiducial simulation. We note the peculiar shape of the streamlines around the transition radius.
In the text

	Fig. 7 Time-averaged angular-momentum-flux streamlines over time-averaged 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 ρvr 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 cs. 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 time-averaged 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 self-similar 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 dot-dashed line is , green dashed line is , and purple dotted line is .
In the text

	Fig. 13 Time-averaged transport coefficients ⟨α⟩ and ⟨vw⟩ 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 (semi-dashed, 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 Magnetic-field transport parameter vB as a function of radius. We note that the outer disc is transporting magnetic field lines inwards (vB < 0).
In the text

Fig. 17 Temporal evolution of Σ in dotted green, Ṁ in dashed blue, and Bz (vertically averaged) in black (full line) at R = 3 au for the fiducial simulation. ζ is calculated at ɛint = 0.3 and shown by a semi-dashed 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 Rint. The profile of is characterised by the presence of gaps in the external disc.
In the text

	Fig. 20 Time-averaged 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 Spatio-temporal 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 Rint. For B5Bin0Am0, we average on 1000 orbhs at Rint 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 = 104.
In the text

	Fig. 24 Spatio-temporal diagrams of Σ (first panel), Bz,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 Spatio-temporal 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 time-averaged on the first 4000 orbits at Rint. The blue lines are the reference runs (left panel: fiducial run, right panel: B5Bin0Am0), and the red-dashed mines are the corresponding runs without the relaxation.
In the text

	Fig. B.2 Spatio-temporal diagrams of Σ (first panel), Bz,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 Time-averaged 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 〈UW〉.
In the text