{'fraction_non_alphanumeric': 0.034114888628370456, 'fraction_numerical': 0.02043767096522079, 'mean_word_length': 4.192141214662518, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnetised Winds in Transition Discs I: 2.5D Global Simulations

Protoplanetary discs (PPDs) are cold, dense and weakly ionised environments that witness the planetary formation. Among these discs, transition discs (TDs) are characterised by a wide cavity in the dust and gas distribution. Despite this lack of material, many TDs strongly accrete onto their central star, possibly indicating that a mechanism is driving fast accretion in TDs cavities. The presence of radially extended 'dead zones' in PPDs recently revived the interest in magnetised disc winds (MDWs), where accretion is driven by a large magnetic field. We propose that TDs could be subject to similar winds, explaining their fast-accreting and long-lived cavities. We present the results of the first 2.5D global numerical simulations of TDs harbouring MDWs using the PLUTO code. We impose a cavity in the gas distribution and consider a power law distribution for the large-scale magnetic field strength. We assume weakly ionised discs subject to ambipolar diffusion, as expected in this range of densities and temperatures. Our simulated TDs reach a steady state with an inner cavity and an outer 'standard' disc. The accretion rates remain approximately constant through the entire discs, reaching $10^{-7}~M_\odot/\text{yrs}^{-1}$ for typical surface density values. The MDW launched from the cavity is more magnetised with a much larger lever arm than the MDW launched from the outer disc. The material in the cavity is accreted at sonic velocities, and the cavity itself is rotating at 70% of the Keplerian velocity due to the efficient magnetic braking imposed by the MDW. Overall, our cavity matches the dynamical properties of an inner jet emitting disc (JED) and of magnetically arrested discs (MADs) in black hole physics. Finally, kinematic diagnostics (wind speeds, orbital and accretion velocities) could disentangle classical photo-evaporation from MDW models.


Introduction
Transition discs 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 10s of a.u. (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 fractions of TDs are accreting onto their protostars. While Najita et al. (2007) quotes 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) shows that accreting TDs have a median accretion rate similar to normal optically thick discs. Manara et al. (2014) finds that TDs accrete like classical 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. (2015van der Marel et al. ( , 2016 find a drop in gas surface density by 2 to 4 orders of magnitude, while the drop in dust surface density goes up to 6 orders of magnitude. Similar results hold in ro-vibrational CO lines, probing the cavity further in, leading to a gas drop of 2 to 4 orders of magnitude in the inner (< 3 a.u.) regions (Carmona et al. 2017).
Hence, 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 explanation to 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 manage to get through 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 to an inner dead zone (Morishima 2012;Gárate et al. 2021). This inner dead zone, typically extending between 1 and 10 a.u., sets the radius of the mass reservoir, Article number, page 1 of 25 arXiv:2205.02126v1 [astro-ph.EP] 4 May 2022 A&A proofs: manuscript no. Transition_Discs_VF_arxiv and therefore the cavity inner edge. While it predicts a fractions of TDs with accretion ratesṀ ∼ 10 −9 M .yrs −1 , 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 zones models. The second type of models requires some 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 which are not necessarily sufficiently 'clean' (Zhu et al. 2011), and predicts that multiple giant planet systems in resonance should be 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 a.u. is probably partially driven by magnetic winds 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 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 modifying significantly 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 jet emitting disc, or JED) is diluted, accreting at sonic velocities, and sustains accretion rates 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, there exists no dynamical model 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 (for example 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 concentrate on 2.5 D models in this first article, and we will discuss 3 D models in a follow up paper. The paper is divided as follows: we first introduce the models equa-tions, physical quantities and numerical setup. We then investigate in depth a fiducial model, which possesses a cavity with a drop of 4 order in 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.

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 where ρ, P, u and B are respectively the density, the thermal pressure and 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: whereb is a unit vector parallel to B and J the electric current, c is the speed of light and η A is the ambipolar diffusivity. No turbulence is added in this model whatsoever. In addition to these equations, the plasma follows the Maxwell's equations and 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) = (G M /r 3 ) 1/2 . It is therefore useful to introduce a deviation from the Keplerian velocity v defined as withΩ(r) ≡ Ω K (r)/ sin 2 θ. We note that the latitudinal dependence ofΩ is somewhat arbitrary, and need not be a particular equilibrium state. Here, our choice ofΩ(r) ensures that our reference Keplerian velocity has constant specific angular momentum on spherical shells and eliminates surface terms which are otherwise present in angular momentum conservation equations (e.g. the last term of equation (16) in Zhu & Stone 2018). This will simplify the interpretation of angular momentum budgets later.

Equation of state and cooling function
As a simplification, we assume the flow follows an ideal equation of state, and is approximately locally isothermal, i.e. T ≈ T eff. (R) where T eff. is a prescribed radial temperature profile. This is achieved solving the energy equation where we have defined a heating/cooling function where τ is the cooling time that equals 0.1 time code unit (see below) and Γ = 1.0001 is the polytropic index of the gas. The target temperature profile is where T 0 is the midplane temperature at the inner radius R int . This choice of cooling function allows us to enforce a chosen temperature profile which 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 c 2 s = P/ρ. 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 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. In the following, we choose T 0 in (10) so that ε = 0.1 .

Integration scheme
The simulations are 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 use a HLLD type Riemann solver to compute the intercell fluxes. In order to ensure the solenoidal constraint (6), we use the constrained transport approach (Kane Yee 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).

Code units and notations
The internal radius is R int = 1, which sets the length code unit, and is chosen to be 1 a.u. while R ext = 50. The time code unit is Ω −1 0 ≡ Ω K (R int ) −1 = 1 which is set to 1/2π years so that G M * = 1 with M * = 1 M , M being 1 solar mass. Therefore Ω K (R) = Ω 0 (R/R int ) −3/2 = R −3/2 . We choose as a unit for the surface density 300 g cm −2 and express the accretion rate in M .yrs −1 . We denote by 'c.u.' the use of code units. We use the subscript X 0 to indicate that the quantity X is considered on the midplane (θ = π/2) and the subscript X p when X is a poloidal quantity.

Dimensionless numbers and definitions
We use the plasma parameter β to quantify the disc magnetisation, defined from the midplane properties of the disc as 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 where v A = B/(4π ρ) 1/2 is the Alfvén speed. We refer to the appendix A for detailed information on the justifications of the model we adopt for Λ A and how we model 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 protoplanetary disc (without a drop in the density profile) and which extends from r ≈ 10 to r = 50.
Finally, we call 'seed' the region defined by r ≤ 1 of our disc, which is at play in our simulations through the inner radial boundary condition.

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 that are uniformly meshed on a logarithmically shaped grid. The colatitude domain is mapped on a stretched grid near the poles (from θ = 0 to θ = 1.279 and from θ = 1.862 to θ = π, with 72 cells in each zone) while the grid is chosen to be uniform around the midplane (from θ = 1.279 to θ = 1.862 with 96 cells) for a total of 240 that 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 constant and equal to 0.1.

Boundary conditions
Outflow boundary conditions are used in the radial direction so that no matter can come from the inner radius. In addition, we add a wave absorbing zone for radii r < 1.5 which damps poloidal motions on an orbital timescale. We detail the impact of this procedure in appendix B.
In these 2.5 D 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 will focus on radii lower than 30.

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 Article number, page 3 of 25 A&A proofs: manuscript no. Transition_Discs_VF_arxiv velocity v ϕ = R Ω K mimic Nelson et al. (2013) to account for the hydrostatic equilibrium, while v r = v θ = 0 initially. These profiles read, without taking into account the cavity yet with ρ 0 being the density at the internal radius. We choose q = −1 and p = −3/2 for the equations (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 initialise the magnetic field using its vector potential A defined so that B = ∇ × A. Following Zhu & Stone (2018), we choose where m = (p + q)/2 = −5/4. This results in a poloidal magnetic field which depends on the radius only The initial strength of the magnetic field is controlled by β out , so that B 0 ∝ β −1/2 out . To add a cavity and mimic a transition disc, we multiply the density profile by a function f that depends on the radius only so that where Σ 0 (R) ∝ R p+1 is a standard surface density profile for protoplanetary disc. The a, b, c coefficients are defined as where R 0 is the radius of the cavity (in code units), n the number of cells on which the transition spans and n δR the corresponding length in code units. Note that while the density profile exhibit 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 transition disc with a strongly magnetised cavity (β in = 1). A typical radial profile of the quantities discussed above are shown in Fig. 1.

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 θ θ ± quantify the integration height as shown in Fig. 2 so that with h int the integration height at radius R given by an integration effective aspect ratio ε int ≡ h int /R. Note that this integration 'height' is not necessarily the disc thickness h. We introduce β as 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 section 3.5.3). Finally, we add the time-average defined by We run the 2.5 D simulations so that we reach 1000 orbits at R = 10 which means ≈ 31000 orbits at R int . If not specified, time-averages are calculated taking into account the whole simulation without the first 4000 orbits at R int to suppress the transient state. Otherwise, we indicate our choice of notation when needed X 1000 being the time-averaged value of X during the last 1000 orbits at R int for example.

Simulations table
All the simulations performed are listed in table 2.4. The characteristic parameters are the external initial plasma parameter β out , the internal initial plasma parameter β in and the initial ambipolar Elsasser number Λ A, 0 . Additionally, we perform 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.

Fiducial simulation
We start by describing in details our fiducial simulation (β out = 10 4 , β in = 1, Λ A = 1 and R 0 = 10), before turning to an exploration of the parameter space.

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 it will be shown in section 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 ∆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 Σ that may suggest that matter is moving inside the cavity at relatively fast speeds. We study in depth these stripes in section 3.5.3. A small accumulation of material is seen close to the inner radius at R ≤ 1.5. We refer to the appendix B for a quantitative discussion on this accumulation. Figure 3 also pictures the evolution of β whose results are similar to the ones for Σ. Inside the cavity, β exhibits a stripedlike pattern with an accumulation close to the internal radius. The edge of the cavity is not smooth at all but varies around its initial value of 10. Though β stays on average around 1 in the cavity, some low values around 10 −2 are reached from time to time. After approximately 4000 orbits at the internal radius, both Σ and β reach a quasi-stationary state.
tions (see Fig. 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 a.u. (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 ambipo-  lar Elsasser number as it 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 do not address the evolution of these rings and gaps any further and instead focus on the dynamics of the cavity.

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 which exhibits a magnetic loop. Inside this loop, the polarity of the azimuthal component is reversed, with B ϕ > 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.

Velocity stream lines
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.

Angular Momentum stream lines
In order to deeper the analysis of the role of the magnetic structure, we concentrate on the time-averaged angular momentum flux, defined by The poloidal flux lines associated to 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.

Accretion rate
The first step to study the accretion in the disc is to investigate the accretion rateṀ defined aṡ The height over which ρv r is integrated has a direct influence oṅ M mostly because of the elbow-shaped stream lines. 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 3 values of ε int . For ε int = 0.3 and around R = 10, the accretion rate is close to zero indicating that the gas does not plunge directly in the cavity from the disc midplane. This radius corresponds to the location of the basis of the Article number, page 7 of 25 A&A proofs: manuscript no. Transition_Discs_VF_arxiv Fig. 8. Accretion rate for different integration height scales with respect to the radius inside the disc. The higher ρv r is integrated the closer to a constant valueṀ is in the cavity. The average value inside the cavity 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 This velocity profile exhibits a well-defined transition between subsonic accretion outside the cavity with v acc. ∼ 10 −3 c s and transsonic accretion inside with v acc. ∼ c s .

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 where we have defined W rϕ and W θϕ respectively the radial and surface stresses by We recall that we use a peculiar definition of the velocity deviation v 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

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, as it will be discussed later in section 3.5.1. The main contribution of the wind mass loss is located in the cavity at R < 5 and is completely compensated by the radial derivative of the accretion rate. Additionally, the 'wind' mass flux turns negative around the cavity edge, which is due to matter being accreted from the outer disc atmosphere (see the 'elbowshaped structure' in the poloidal streamlines).
In order to quantitatively account for the role of the wind, we construct the mass loss rate parameter ζ = ζ + +ζ − (Lesur 2021b), where ζ + and ζ − are defined by where 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 study the values of ζ at z 0 = 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 −5 , while it peaks at 2.9 × 10 −2 in the inner part. We find two zones where ζ < 0. One is close to the inner boundary and  Fig. 10. Mass conservation for ε int = 0.9. The time derivative of Σ remains perfectly constant and equal to zero inside the cavity and only gets 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 use a moving average for better visibility and the quantities are time-averaged on a sample selection of output files that do not contain all the timesteps computed by the code.
probably a boundary condition artefact, while the other extends from R ≈ 5 to R ≈ 17 a.u. 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 selfsimilar 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 by essence.
The wind mass loss rate parameter is smaller than the selfsimilar scaling in the outer disc by a factor of a few. This discrepancy is probably due to the influence of the cavity magnetosphere that compresses the disk magnetosphere, resulting in a deviation of ζ from the self-similar result. Moreover, it seems that the further we move outward, the closer we get to the selfsimilar 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 has not been 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 ζ ext = ζ 0, ext β a ext , with a ext = −0.69 and ζ 0, ext = 0.24. Another one is then calculated for the inner part only, ζ int = ζ 0, int β a int , with a int < 0, so that the final profile is given by We get a int = −0.20 and ζ 0 int = 0.018. Such a model, with a ext − a int < 0 allows to recover both the 2 previous regimes with a reasonably accurate depiction of the disc. The final profile exhibits a transition occurring at β t ≈ 5 which is close to  Fig. 11. ζ parameter for ε int = 0.6. The self similar fit shown here is for comparison only and was obtained also 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.
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 properly catches both the inner and external parts of the disc.

Angular momentum conservation
We show in Fig. 12 the terms involved in the angular momentum conservation equation (28), time-averaged and multiplied by r −3/2 for better 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 2 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 introduce Shakura & Sunyaev (1973) α parameter. It must be noted that the origin of this stress is in no way solely linked to turbulence and considerably driven by the laminar structure of the magnetic wind. The appendix D details the origin of the stress and sheds light on the turbulent vs. laminar origin of α. Nevertheless, the α parameter can still be used in this wind model whose definition when time-averaged is The corresponding profile is given in Fig. 13, where ε int = 0.9. In the external part of the disc, α = 49 ± 5 × 10 −4 while it reaches a maximum value inside the cavity α = 13 ± 5. Following a similar procedure as the one for α, we define a dimensionless number associated to the surface stress component, υ W . As for ζ, we define υ W, ± which are chosen to be positive for angular momentum leaving the disc on both sides: We show the dependence of υ W on R in Fig. 13. In the external disc, υ W = 2.3 ± 1.1 × 10 −4 while it rises up to υ W = 1.0 ± 0.1 inside the cavity. The same observations as for ζ are drawn for both α and υ W . Therefore, two separated 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 4 , indicating that the dynamical properties of the outer disc are not perturbed by the presence of the cavity. On the contrary, the second regime describes the inner part of the disc with fast accretion and high values for α and υ W , which are both of the order of unity. Table 2 displays the transport coefficients values for all the simulations.

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 lines (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 fields lines attached in the cavity and in the outer disc.
In the following, we select a field line anchored in the disc midplane at R ≡ R w . The corresponding Keplerian angular velocity is Ω w while B w is the poloidal magnetic field at the midplane. We then consider the following invariants, built on timeaveraged quantities and listed in table 3 -The mass loading parameter which accounts for the quantity of matter that escapes the disc with the wind -The rotation parameter -The magnetic lever arm that accounts for the angular momentum that is removed from the disc by the wind Of course, these invariants echoes the transport coefficients definitions previously used to describe the disc and one expects κ ≈ β ζ/4 ε and λ ≈ 1 + ε υ W /ζ (Lesur 2021a).
To compute these invariants, we arbitrarily choose one field line in the cavity (referred to as 'in') leaving the midplane at R in = 5 a.u. and one in the external disc (referred to as 'ext') leaving the midplane at R ext = 15 a.u (see the first panel of Fig. 14). 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 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 Ω K 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 −2 , κ ext = 2.5, λ in = 23 and λ ext = 3.2. These values are coherent with the transport coefficients computed in previous sections. We also note that the values of κ and λ in the cavity match some of the historical solutions of Blandford & Payne (1982) (see their figure 2), which correspond to super-Alfvénic and collimated outflows. These values are also consistent with the magnetic outflow solutions of Ferreira (1997) (see figure 3). Hence, the cavity we find quantitatively matches the inner JED proposed by Combet & Ferreira (2008).

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). We start here our exploration of time variability by focusing on the secular evolution, beginning with a discussion of the cavity expansion.

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 getsṘ whereṘ 0 is the cavity edge 'velocity' and δṀ and δΣ are the jump in accretion rate and surface density at the cavity edge. By calculatingṀ and Σ around R 0 , we findṘ 0 = 1.8 × 10 −5 while evaluating directly the cavity edge motionṘ 0 yieldsṘ 0 = 1.4 × 10 −5 (both in c.u.). Therefore, the cavity is expanding because of the slight mismatch in accretion rate observed in Fig. 8.

Magnetic field transport
To interpret the time evolution of the magnetic field, we study the transport of magnetic flux inside the disc and define a flux function ψ such that ψ(r, t) = R 2 int π/2 0 B r (R int , θ, t) sin θ dθ − r r=R int r B θ (r, π/2, t) dr.
Article number, page 11 of 25 A&A proofs: manuscript no. Transition_Discs_VF_arxiv 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 on 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 get a quantitative estimate of the field line advection speed, we first note that the evolution equations for ψ read Following Guilet & Ogilvie (2014), we rewrite these evolution equations as an advection equation for ψ where we have defined the 'field advection velocity' Eventually, we define a dimensionless advection parameter ν B = v ψ /v K which quantifies the advection speed (Bai & Stone 2017). In this framework, positive values of ν B implies an outward transport field while negative values trace inward field transport. We show the radial dependence of ν B in Fig. 16. In the external disc we find that the magnetic field is advected inwards with a velocity v ψ = −2.6 × 10 −3 v K . ν B changes its sign multiple times in the cavity, but remains negative close to the cavity edge, -ν B Fig. 16. The magnetic field transport parameter ν B as a function of radius. Note that the outer disc is transporting magnetic field lines inwards (ν B < 0).
between R ≈ 7 and R ≈ 11 a.u. where v ψ = +3.2 × 10 −3 v K . Such a result is in accordance with Fig. 15 and indicates that field lines are converging at the transition radius with opposite vertical polarity. In the external parts, ν B is negative and v ψ = −2.6 × 10 −3 v K so that vertical magnetic field pointing upwards is advected. We note that this inwards advection of the outer disc field lines is in sharp contrast to other work which focused on 'full' discs (Bai & Stone 2017;Lesur 2021b). We will come back to this discrepancy in the discussion.

Fast variability of the cavity
Up to this point, we have mostly considered time-averaged quantities, and ignored fast variability. While our numerical solution 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) whose origin ought to be clarified. Figure 17 shows such a temporal evolution of Σ,Ṁ and B z at R = 3. These profiles encounter sharp fluctuations over time, chaotically distributed. 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Ṁ are correlated. When an inflow of matter crosses the cavity, Σ peaks as well aṡ M which 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 to Fig. 18 that shows the density and poloidal magnetic field lines at different times. On 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 to leave the disc in the wind. Concerning the filament as well as the bubbles formation, we detect a current sheet localised at the location of the filaments, where the total magnetic field cancels (B ϕ = 0 at the edge of the magnetic loop and B p = 0 because two antiparallel poloidal field lines meet at the elbow shape structure location). It is therefore a possibility that these structures form due to magnetic reconnection. Focusing on the falling material, we see it reaching the edge of the cavity on the third panel before crossing it on the next one. When the gas crosses the cavity, the disc oscillates locally above and below the midplane and is therefore highly dynamical. With a slight delay (last three panels), we see an outflow emerging from the cavity and the wind density increases. Such an observation exhibits the link between wind and accretion (see Fig. 10). The ejection of gas from the cavity is not constant with respect to time and occurs occasionally with burst events for which ζ eventually peaks at 0.1. This explain 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 2 local disc height 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.

Magnetic Rayleigh Taylor instability
To account for the formation and stability of the bubbles of matter at the cavity edge, we explore the possibility of having a mag-netic 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 refer to the analysis of Spruit & Taam (1990); Spruit et al. (1995); Stehle & Spruit (2001), which assume an infinitely thin disc. We reformulate the instability criterion of Spruit et al. (1995) (see their equation 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 β < β crit. 0.0355. Figure 3 shows that β is of the order of 0.1 in the cavity and rarely go 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 reach β < 0.1, as during some of the bursts.

Ambipolar Diffusion
We check the influence of Λ A, 0 in the simulation B4Bin0Am1, which is the same as the fiducial one except for the initial value of Λ A which is set to 10.

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 a.u. before expanding back up to R 4 a.u. in a few thousands of orbits at R int . Overall, the transient state lasts A&A proofs: manuscript no. Transition_Discs_VF_arxiv 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 inner boundary. Some gaps merge with one another, so that only 3 of them remain after 15000 orbits at R int , similarly to Cui & Bai (2021). We deserve the study of the interaction between these gaps and the cavity to 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 disappear below, except for a small set of stream lines 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. 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.

Transport coefficients and wind invariants
The accretion rate remains constant in the whole disc with a value close to 1.2 × 10 −7 M .yrs −1 and an accretion velocity that is still subsonic in the outer disc and peaks up to 2 c s at the internal radius. Therefore, the accretion picture is identical to the one for the fiducial run with an internal transsonic regime connecting through the cavity edge to a weakly magnetised wind.
Regarding the wind, we obtain a highly mass loaded field line in the external disc that removes little angular momentum (λ ext = 1.5 and κ ext = 5.0) and a lighter one in the internal disc that carries a massive load of angular momentum (λ in = 4.9 and κ 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 3 times lower and the external 3 times higher.

Influence of the initial plasma parameter
We study the impact of the plasma parameter varying both its internal β in and external β out initial value.

Role of the external initial plasma parameter
We explore here how the outer disc magnetisation impacts the general properties of the system. We vary the initial value of β between β out = 10 3 (run B3Bin0Am0) and β out = 10 5 (run B5Bin0Am0).
General observations for B5Bin0Am0: The spatiotemporal evolutions of Σ and β are shown in the left panels of Fig. 21. Right at the beginning of the simulation, a burst of matter appears in the cavity which is subsequently refilled. Its radius then remains fixed at ∼ 4 a.u. 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 overall similar to the fiducial one.
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 to the 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 only occur in the weakly magnetised (β out = 10 5 ) simulations.
Lastly, we add that estimatingṘ 0 for this simulation is too difficult, since the cavity edge barely move during the entire simulation.
General observations for B3Bin0Am0: In contrast to the run B5Bin0Am0, the cavity quickly expands up to R ≈ 15 a.u. and keeps growing during all the simulation, faster than in the fiducial run (see the right panels of Fig. 21). We estimate its velocity asṘ 0 ≈ 3.0 × 10 −5 c.u., which is about 3 times faster than the fiducial run. We get δṀ ≈ 2.3 × 10 −4 and δ Σ ≈ 4.2 × 10 −2 (both in code units) so that Eq. 38 givesṘ 0 ≈ 4.2 × 10 −5 c.u., where we choose R 0 ≈ 20. The simple model we use seems to overestimate the widening velocity of the cavity but still gives a correct order of magnitude.
The time-averaged surface density from the fiducial run, B3Bin0Am0 and B5Bin0Am0 are shown in Fig. 22, which shows 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 (the subsection 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 aṡ where we have defined the accretion velocities v acc. ≡ M/2π R 0 Σ, and we have 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 there exists a valueβ for whichṘ 0 = 0, we get the scalinġ where we have used the fact that Σ in /Σ out = β in /β out = β −1 out in our setup. The relation (44) shows that for β out β , we have approx-imatelyṘ 0 ≈ v acc., in β −1 out , indicating that the cavity expansion speed should increase as β out gets lower, which is precisely what we observe for B3Bin0Am0. For β out β we get on the contrarẏ R 0 ≈ −v acc., inβ σ−1 β −σ out , showing a change of sign (hence a contraction of the cavity), albeit with a reduced speed. This regime might correspond to B5Bin0Am0, indicating thatβ 10 4 .

Role of the internal initial plasma parameter
To study the impact of β in , we run a set of simulations that covers all the possible initial gaps β in /β out where log β out ∈ {3, 4, 5} and log β in ∈ 0; log β out . We compare each result to the one obtained with β in = 1 and the corresponding value of β out . A striking result is the fact that the disc 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 get transsonic accretion as it is mentioned in Wang & Goodman (2017). We illustrate this statement with Fig. 23 for the particular case of β out = 10 4 . We focus on the transient state of B4Bin3Am0 in Fig. 24. The transition is due to matter leaving the cavity because of the fast Fig. 19. Σ(R, t) and β(R, t) of B4Bin0Am1. The cavity stands during the whole simulation though its edge falls down to R ≈ 2 a.u. during the transient, before broadening up to R ≈ 4 a.u. in a few thousands of orbits at R int . The profile of β is characterised by the presence of gaps in the external disc. NameṀ (10 −7 M .yrs −1 ) ζ in (10 −2 ) ζ ext (10 −5 ) α in α ext (10 −3 ) υ W, in υ W, ext (10 −4 ) B4Bin0Am0 Name  Table 3. MHD invariants for a subset of simulations, computed with time-averaged quantities on the last 1000 orbits at R int .
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 a 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 there should exist a threshold regarding the value of β in above which no cavity is able to form. Considering Fig. 23, it seems that this threshold is 10 3 .
From these observations, we deduce that the cavity is regulated by the value of the plasma parameter which must take a value close to 1. The reason for this regulation is not entirely clear and we add a word of caution regarding the role of the inner radial boundary condition, especially with respect to the magnetic field transport at R int . We discuss this influence in appendix B.
Referring to the sections 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.

Zoom with a larger cavity radius
We perform a simulation with a double-sized cavity (R 0 = 20 a.u.) in order to check the impact of the cavity size. The simulation was integrated for 1000 orbits at R = 10 a.u. so that it reaches 355 orbits at R = 20 a.u. 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 (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.

Discussion and comparison with previous work
We model transition discs sustained by MHD winds by performing 2.5 D 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 weaklymagnetised ambipolar-dominated protoplanetary discs (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 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 ≈ 0.13 G (see Eq. 17, with β out = 10 4 ) so that initially, B z ≈ 1.25 mG at R = 42 a.u. in our simulations, which is of the same order of magnitude as the upper limit of B z (R = 42 a.u.) = 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 matches quantitatively 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 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, we refer to the text for more explanations. Fig. 22. Impact of the initial external magnetisation on the surface density. We average the profile on the last 1000 orbits at R int . For B5Bin0Am0, we average on 1000 orbits at R int occurring between the 2 burst events seen in Fig. 21. 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 ac- Fig. 23. Impact of the internal initial magnetisation on the plasma parameter for β out = 10 4 . cretion, implies that the kinematics of these cavities have singular observational signatures. Fast accretion kinematics have been observed in some transition discs (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 or orbits and which 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 figure 5, first row and first column panel). While this is by no mean a 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 2 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, recalling to some extent the magnetospheric accretion observed in young stars (Bouvier et al. 2007;Pouilly et al. 2020;Bouvier et al. 2020b,a) 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 have checked whether the time variability seen in our simulations could be due to a magnetic RTI. We have studied 2 criteria for the RTI, in the form of a radial interchange of poloidal field lines (see section 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 on 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) shown that MADs could be regulated by magnetic RTI leading to magnetically chocked accretion flows (MCAF). The MAD model is also associated to the formation of plasmoids by reconnection events (Ripperda et al. 2022). These features are recovered in our models of transition discs, despite the fact that we have 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 nonaxisymmetric hydrodynamical instabilities such as the Rossby Wave Instability (RWI) (Lovelace et al. 1999;Li et al. 2000) at the cavity edge is still open to debate in a magnetised environment (Bajer & Mizerski 2013). We will address these points using full 3 D 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 is sometime expelling 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 either include an inner turbulent disc, or possibly the magnetospheric interaction with the central star.
A possible limitation of our model one could raise 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 2 ). This saturation is different from the saturation by 3D turbulence observed in the ideal MHD regime. It is suggested that the MRI saturates in very similar ways in 3D and 2D under strong ambipolar diffusion (see e.g. Béthune et al. 2017;Cui & Bai 2021). This is also confirmed by our own 3D simulations which will be published in a forthcoming paper. Hence, the fact that our simulations are 2.5D have a very limited impact on the turbulent transport one may observe.
Note that our simulations used a simplified treatment of thermodynamics and ionisation chemistry. More numerically involved models, such as Wang & Goodman (2017), use a refined computation of the ionisation fraction and Λ A 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 0 at 16 a.u. (Fig. 2, panels 6 and 7). Regarding our profile of Λ A ≈ 1 − 10, our work is similar to their models 2 (with L X = 10 29 erg s −1 ) and 6 (where T 0 = 30 K). Therefore, we anticipate that an increase of 2 orders of magnitude for L X would lead to Λ A > 10 2 in most of our cavity. Such a change would greatly alter the dynamical regime of the cavity since MRI would then play a significant role (see the appendix A and Blaes & Balbus 1994;Bai 2011). However, the role of the temperature is less straightforward and seems to have a little impact on Λ A .
Additionally, dust plays a significant role in Wang & Goodman (2017) regarding the ionisation of the disc. As a matter of fact, only their models with dust reach low values of Λ A . The effect of dust in transition discs 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, that is known to create nonaxisymmetric structures at the cavity edge (Bi & Fung 2022) or an inner rim with an accumulation of matter due to photophoresis (Cuello et al. 2016).

Conclusions
We performed 2.5 D global numerical simulations of transition discs 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 have modelled strongly accreting transition discs that reach a quasi steady state that last for at least tens of kyrs. 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 selfregulated. In particular, it is insensitive to a change in the initial internal magnetisation and is characterised by 0.1 β int 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 transition discs could be an instance of MADs applied to protoplanetary discs.
The aim of this appendix is to model the ambipolar diffusivity spatial dependence in both a transition disc and a standard protoplanetary disc (i.e. without cavity). The general procedure to reach such a result follows and adapts the main calculation steps that are presented in Combet et al. (2010). As assumed in Eq. 4, only the ambipolar diffusivity does appear in the MHD equations, which we assume is the dominant non-ideal effect in the regime of discs we use at R ≥ 10 ( Riols et al. 2020;Simon et al. 2015). Therefore, the only momentum exchange that occurs between particles happens only between ions and neutrals. In a plasma made of molecular ions, electrons and neutrals, the ambipolar diffusivity is given by (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 the ion-neutral collision rate whose value is (Bai 2011) σv in = 2,0 × 10 −9 m H µ with ζ i the ionisation rate and δ the dissociative recombination rate. Following Fromang et al. (2002), we take δ = 3 × 10 −6 T −1/2 cm 3 s −1 . (A.6) In this toy model we then have (Lesur et al. 2014) where ξ FUV accounts for the far UV photons contribution that we model following Perez-Becker & Chiang (2011) as with Σ the column density computed from the star to the point of interest.
• Cosmic-ray ionisation following (Umebayashi & Nakano 1980) ζ CR = ζ CR,0 e −Σ col. /96 g cm −2 s −1 , (A.10) where ζ CR,0 = 10 −17 s −1 and Σ col. is the matter column density above and below the point of interest. • Radioactive decay is assumed constant (Umebayashi & Nakano 2008) ζ rad. = 10 −19 s −1 . (A.11) Combining the equations A.9, A.10 and A.11, we obtain ζ i = ζ X + ζ CR + ζ rad. , paving the way to finally get Λ A using equations 13, A.6 and A.7. Note that due to the dependency of η A and v A on the norm of the magnetic field, this latter cancels and does not need to be computed to get Λ A . The previous calculations can be performed either for a standard protoplanetary disc or for a transition disc. The only thing that needs to be changed to account for such discs is the surface density profile, where Eq. 18 allows to consider or not the effects of the cavity.
The results of such calculations are displayed in Fig. A.1 that represents the spatial dependency of Λ A in both a standard protoplanetary disc and a transition disc. Though these 2 profiles look different at first glance, a deeper investigation reveals that the values taken by Λ A in the discs remain pretty much close to unity in both cases, while the general trend of Λ A in a standard protoplanetary disc is recovered even in the case of a transition disc (Thi et al. 2019). Moreover, Λ A remains fairly below the critical value Λ A, crit. = 10 2 with or without a cavity. Λ A must stay below Λ A, crit. so that the MRI effects are negligible (Blaes & Balbus 1994;Bai 2011). Therefore, assuming a characteristic value of Λ A, 0 = 1 captures within a reasonable accuracy the physics of ambipolar diffusion and the cavity does not alter the ambipolar Elsasser number profile. The results we get from this simple toy model are to be compared to the more detailed work of Wang & Goodman (2017) where many chemical species are taken into account to compute the ambipolar Elsasser number inside the cavity of a wind-driven transition disc.
Following Lesur (2021b) and Thi et al. (2019), we implement the profile of Λ A so that 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 constant and equal to 3 h. Λ A, 0 remains a free parameter (see 2.4 for more details). Additionally, a cutoff is used for the η A profile so that if η A > η A, max , the value of η A is replaced by η A, max ≡ 10 ε 2 in code units, such a choice being reflected on the Λ A profile with Eq. 13.
≡ 2 S 2 , (C.1) where S is the shear that we approximate with S 2 = 9/4 Ω 2 and g m is B + R is the radial component of the magnetic field at the disc surface. Let us rewrite the previous expression in terms of β, q (defined with B + R = q B z ) and δ (defined as δ = −d ln Σ/d ln R).
where X denotes the derivative of X with respect to R. With P = c 2 s ρ = (h Ω K ) 2 Σ/(  Therefore, the instability criterion becomes ε being constant in the disc as well as β inside the cavity. Ω K varies as R −3/2 and Σ as R −δ so that S 2 < 4 ε Ω 2 K √ 2π β q (−δ + 1 + δ/2) (C.7) By taking S 2 /Ω 2 K = 9/4, the RTI can be triggered when If we now assume that δ = q = 1 for simplicity, we finally get β < 8 ε 9 √ 2π ≈ 0.355 ε = 0.0355 ≡ β crit. , (C.9) where ε = 0.1. Figure C.1 compares the time-averaged values of β with the criterion given in Eq. C.9. The value of β crit. is anyhow below the time-averaged values of β. Though this simple analysis makes it difficult to be definitive 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 Article number, page 24 of 25 Étienne Martel and Geoffroy Lesur : Magnetised Winds in Transition Discs I: 2.5 D Global Simulations stresses with a turbulent and a laminar part. In this prospect, we introduce the deviation to the temporal mean such that δX ≡ X − X .
(D. The results are shown in Fig. D.1. The laminar contribution is the major one for υ W in the whole disc so that 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 on the turbulent component of α too 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 υ W is definitely dominated by its the laminar part and due to the wind.