# Magnetised winds in transition discs. I. 2.5D global simulations

Étienne Martel, Geoffroy Lesur

► To cite this version:

Étienne Martel, Geoffroy Lesur. Magnetised winds in transition discs. I. 2.5D global simulations. Astronomy and Astrophysics - A&A, 2022, 667, 10.1051/0004-6361/202142946. insu-03860267

HAL Id: insu-03860267

https://insu.hal.science/insu-03860267

Submitted on 19 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Magnetised winds in transition discs

## I. 2.5D global simulations*

Étienne Martel and Geoffroy Lesur

Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
e-mail: etienne.martel@univ-grenoble-alpes.fr

Received 17 December 2021 / Accepted 3 May 2022

## ABSTRACT

**Context.** Protoplanetary discs are cold, dense, and weakly ionised environments that witness planetary formation. Among these discs, transition discs (TDs) are characterised by a wide cavity (up to tens of au) in the dust and gas distribution. Despite this lack of material, a considerable fraction of TDs are still strongly accreting onto their central star, possibly indicating that a mechanism is driving fast accretion in TD cavities.

**Aims.** The presence of radially extended ‘dead zones’ in protoplanetary discs has recently revived interest in magnetised disc winds (MDWs), where accretion is driven by a large magnetic field extracting angular momentum from the disc. We propose that TDs could be subject to similar disc winds, and that these could naturally explain the fast-accreting and long-lived cavities inferred in TDs.

**Methods.** We present the results of the first 2.5D global numerical simulations of TDs harbouring MDWs using the PLUTO code. We imposed a cavity in the gas distribution with various density contrasts, and considered a power-law distribution for the large-scale magnetic field strength. We assume the disc is weakly ionised and is therefore subject to ambipolar diffusion, as expected in this range of densities and temperatures.

**Results.** We find that our simulated TDs always reach a steady state with an inner cavity and an outer ‘standard’ disc. These models also maintain an approximately constant accretion rate through the entire structure, reaching $10^{-7} M_{\odot} \text{yr}^{-1}$ for typical surface density values. The MDW launched from the cavity is more magnetised and has a significantly larger lever arm (up to 10) than the MDW launched from the outer disc. The material in the cavity is accreted at sonic velocities, and the cavity itself is rotating at 70% of the Keplerian velocity due to the efficient magnetic braking imposed by the MDW. Overall, our cavity matches the dynamical properties of an inner jet emitting disc (JED) and of magnetically arrested discs (MADs) in black-hole physics. Finally, we observe that the cavity is subject to recurring accretion bursts that may be driven by a magnetic Rayleigh-Taylor instability of the cavity edge.

**Conclusions.** Some strongly accreting TDs could be the result of magnetised wind sculpting protoplanetary discs. Kinematic diagnostics of the disc or the wind (orbital velocity, wind speeds, accretion velocities) could disentangle classic photo-evaporation from MDW models.

**Key words.** accretion, accretion disks – protoplanetary disks – magnetohydrodynamics (MHD) – methods: numerical

## 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

<footnote>* The data underlying this article will be shared on reasonable request to the corresponding author.</footnote>
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 $\dot{M} \sim 10^{-9} M_{\odot} yr^{-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 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 $\dot{M} \propto \Sigma^{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 $\Sigma$-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 $\alpha$ 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, non-ideal 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

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \quad (1) $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla P - \rho \nabla \Phi_* + \frac{\mathbf{J} \times \mathbf{B}}{c}, \quad (2) $$

$$ \partial_t \mathbf{B} = -\nabla \times \mathcal{E}, \quad (3) $$

where $\rho$, $P$, $\mathbf{u}$, and $\mathbf{B}$ are, respectively, the density, the thermal pressure, the plasma velocity, and magnetic field. $\Phi_* = -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:

$$ \mathcal{E} = -\mathbf{u} \times \mathbf{B} - \frac{4\pi}{c} \eta_A \mathbf{J} \times \hat{\mathbf{b}} \times \hat{\mathbf{b}}, \quad (4) $$

where $\hat{\mathbf{b}}$ is a unit vector parallel to $\mathbf{B}$, $\mathbf{J}$ is the electric current, $c$ is the speed of light, and $\eta_A$ is the ambipolar diffusivity. No turbulence was added in this model whatsoever. In addition to these equations, the plasma follows the Maxwell equations:

$$ \nabla \cdot \mathbf{B} = 0 \quad (5) $$

and

$$ \mathbf{J} = \frac{c}{4\pi} \nabla \times \mathbf{B}. \quad (6) $$

We place ourselves in a spherical coordinate system $(r, \theta, \phi)$ centred on the star. For convenience, we also introduce the cylindrical coordinates $R = r \sin\theta$, $\vartheta = \phi$ and $z = r \cos\theta$.

Since we work in a thin disc, the azimuthal angular velocity $\Omega$ is expected to be close to the Keplerian angular velocity $\Omega_K(r) = (G M_*/r^3)^{1/2}$. It is therefore useful to introduce a deviation from the Keplerian velocity $v$, defined as

$$ v = \mathbf{u} - r \sin\theta \tilde{\Omega}(r) e_{\varphi}, \quad (7) $$
with $\tilde{\Omega}(r) \equiv \Omega_K(r)/\sin^2\theta$. We note that the latitudinal dependence of $\tilde{\Omega}$ is somewhat arbitrary and need not be a particular equilibrium state. Here, our choice of $\tilde{\Omega}(r)$ 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 \approx T_{\text{eff.}}(R)$ where $T_{\text{eff.}}$ is a prescribed radial temperature profile. This is achieved by solving the energy equation

$$ \partial_t P + \mathbf{u} \cdot \nabla P + \Gamma P \nabla \cdot \mathbf{u} = \Lambda, \quad (8) $$

where we defined a heating and cooling function

$$ \Lambda = \frac{P T - T_{\text{eff.}}}{T \tau}, \quad (9) $$

where $\tau$ is the cooling time that is set to 0.1 time code unit (see below) and $\Gamma = 1.0001$ is the polytropic index of the gas. The target temperature profile is

$$ T_{\text{eff.}}(R) = T_0 \left( \frac{R}{R_{\text{int}}} \right)^{-1}, \quad (10) $$

where $T_0$ is the midplane temperature at the inner radius $R_{\text{int}}$. 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 $c_s^2 \equiv P/\rho$. It can be shown that as a result of the vertical hydrostatic equilibrium, $c_s$ and $\Omega_K$ are related to the vertical disc thickness $h(R)$ through

$$ h(R) = c_s(R)/\Omega_K(R). \quad (11) $$

Assuming the disc is at thermal equilibrium ($T = T_{\text{eff.}}(R)$), we have $c_s \propto R^{-1/2}$, and hence the disc aspect ratio $\epsilon \equiv h/R$ is constant. For the following, we chose $T_0$ in (10) so that $\epsilon = 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 $R_{\text{int}} = 1$, which sets the length code unit and is chosen to be 1 au while $R_{\text{ext}} = 50$. The time code unit is $\Omega_0^{-1} \equiv \Omega_K(R_{\text{int}})^{-1} = 1$, which is set to $1/2\pi$ years so that $G M_* = 1$, with $M_* = 1 M_\odot$, $M_\odot$ being 1 solar mass. Therefore, $\Omega_K(R) = \Omega_0 (R/R_{\text{int}})^{-3/2} = R^{-3/2}$. We chose 300 g cm$^{-2}$ as a unit for the surface density and express the accretion rate in $M_\odot \text{yr}^{-1}$. We denote code units as c.u. hereafter. We use the subscript $X_0$ to indicate that the quantity $X$ is considered on the midplane ($\theta = \pi/2$) and the subscript $X_p$ when $X$ is a poloidal quantity.

### 2.2.3. Dimensionless numbers and definitions

We used the plasma parameter $\beta$ to quantify the disc magnetisation, which is defined from the midplane properties of the disc as

$$ \beta = \frac{8\pi P_0}{B_{p,0}^2}. \quad (12) $$

When considering the initial state of a given simulation, we refer to the initial magnetisation inside the cavity as $\beta_{\text{in}}$ and to the initial magnetisation in the external part of the disc as $\beta_{\text{out}}$. The second key parameter of this study is the strength of ambipolar diffusion, quantified with the Elsasser number

$$ \Lambda_A \equiv \frac{v_A^2}{\Omega_K \eta_A}, \quad (13) $$

where $v_A = B/(4\pi\rho)^{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 $\Lambda_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 \in [1; 50]$ and $z/R \in [-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 \approx 10$ to $r = 50$. Finally, we refer to the region defined by $r \le 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 \equiv R_{\text{int}}$ to the external one $r \equiv R_{\text{ext}}$ and are uniformly meshed on a logarithmically shaped grid. The colatitude domain is mapped on a stretched grid near the poles (from $\theta = 0$ to $\theta = 1.279$ and from $\theta = 1.862$ to $\theta = \pi$, with 72 cells in each zone), while the grid is chosen to be uniform around the mid-plane (from $\theta = 1.279$ to $\theta = 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 $\epsilon$ 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_θ = 0$ initially. These profiles read, without taking into account the cavity,

$$ \rho(R, z) = \rho_0 \left(\frac{R}{R_{\text{int}}}\right)^p \exp\left[\left(\frac{\Omega_K(R) R^3}{c_s(R)}\right)^2 \left(\frac{1}{\sqrt{R^2 + z^2}} - \frac{1}{R}\right)\right], \quad (14) $$

$$ v(R, z) = v_K(R) \left[ (p+q) \left( \frac{c_s(R)}{\Omega_K(R) R^2} \right)^2 + (1+q) - \frac{q R}{\sqrt{R^2 + z^2}} \right]^{1/2}, \quad (15) $$

with $ρ_0$ being the density at the internal radius. We chose $q = -1$ and $p = -3/2$ for the Eqs. (14) and (15), which is consistent with self-similar stationary disc solutions (Jacquemin-Ide et al. 2021).

The initial vertical magnetic field follows a power law $B_z \propto 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

$$ A_φ = \begin{cases} \frac{1}{2} B_0 R & \text{if } R \le R_{\text{int}} \\ B_0 \frac{R_{\text{int}}^2}{R} \left( \frac{1}{2} - \frac{1}{m+2} \right) + R \left( \frac{R}{R_{\text{int}}} \right)^m \frac{1}{(m+2)} & \text{if } R > R_{\text{int}} \end{cases}, \quad (16) $$

where $m = (p+q)/2 = -5/4$. This results in a poloidal magnetic field that depends on the radius only:

$$ \mathbf{B} = B_0 \left( \frac{R}{R_{\text{int}}} \right)^m \mathbf{e}_z. \quad (17) $$

The initial strength of the magnetic field is controlled by $β_{out}$, so that $B_0 \propto β_{out}^{-1/2}$.

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

$$ \Sigma(R) = f(R) \times \Sigma_0(R), \quad (18) $$

with

$$ f(R) = a \left( 1 - c \tanh \left[ b \left( 1 - \frac{R}{R_0} \right) \right] \right), \quad (19) $$

where $Σ_0(R) \propto R^{p+1}$ is a standard surface density profile for a protoplanetary disc. The $a$, $b$, and $c$ coefficients are

defined as

$$ \left\{ \begin{aligned} b &= \frac{2}{n} \left( \frac{\delta R}{R_0} \right)^{-1} \\ a &= \frac{\beta_{\text{in}}/\beta_{\text{out}} + \tanh(b)}{1 + \tanh(b)}, \\ c &= \frac{1 - \beta_{\text{in}}/\beta_{\text{out}}}{\beta_{\text{in}}/\beta_{\text{out}} + \tanh(b)} \end{aligned} \right. $$

where $R_0$ is the radius of the cavity (in code units), $n$ the number of cells on which the transition spans, and $n \delta 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 \propto Σ(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 $θ$:

$$ \bar{X}(r, t) = r \int_{\theta_{+}}^{\theta_{-}} X(\mathbf{r}, t) \sin \theta \, d\theta. \quad (20) $$

$θ_{±}$ quantify the integration height as shown in Fig. 2 so that

$$ \frac{\theta_{-} - \theta_{+}}{2} = \arctan \left( \frac{h_{\text{int}}}{R} \right) = \arctan \varepsilon_{\text{int}}, \quad (21) $$

with $h_{int}$ being the integration height at radius $R$ given by an integration effective aspect ratio $ε_{int} \equiv h_{int}/R$. We note that this integration ‘height’ is not necessarily the disc thickness $h$. We introduce $\bar{\beta}$ as

$$ \bar{\beta} \equiv \frac{8\pi \Sigma \bar{c}_{s,0}^2}{\sqrt{2\pi R \varepsilon \left(\bar{B}_r^2 + \bar{B}_{\theta}^2\right)}}, \quad (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

$$ \langle X(\mathbf{r}) \rangle = \frac{1}{T} \int_{t_0}^{t_0+T} X(\mathbf{r}, t) \, dt. \quad (23) $$

We ran the 2.5D simulations so that we reach 1000 orbits at $R = 10$, which means ≈31 000 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, $\langle X \rangle_{1000}$ being the time-averaged value of $X$ during the last 1000 orbits at $R_{int}$, for example.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Initial and time-averaged profile of $\Sigma$ (top left panel), $\bar{\beta}$ (top right), $B_{z,0}$ (bottom left), and $v_{\varphi,0}$ (bottom right), with respect to R.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic view of the disc represented in orange. $\theta_{\pm}$ define the vertical integration surface and $h_{\text{int}}$ is the integration scale height at a given radius R.</figcaption></figure>

## 2.4. Simulation table

All the simulations performed are listed in Table 1. The characteristic parameters are the external initial plasma parameter $\beta_{\text{out}}$, the internal initial plasma parameter $\beta_{\text{in}}$, and the initial ambipolar Elsasser number $\Lambda_{A,0}$. 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 ($\beta_{\text{out}} = 10^4$, $\beta_{\text{in}} = 1$, $\Lambda_A = 1$ and $R_0 = 10$), before turning to an exploration of the parameter space.

### 3.1. Evolution of surface density and plasma magnetisation

We first look at the temporal evolution of the surface density (Fig. 3). We find that the cavity stands during the whole simulation as its radius remains close to its initial value. As we show in Sect. 3.5.1, the cavity tends to expand slightly. The cavity

<table><caption>Table 1. Simulation information. B4Bin0Am0 is the fiducial simulation.</caption><thead><tr><th>Name</th><th>β<sub>out</sub></th><th>β<sub>in</sub></th><th>Λ<sub>A,0</sub></th><th>R<sub>0</sub> (au)</th></tr></thead><tbody><tr><td>B4Bin0Am0</td><td>10<sup>4</sup></td><td>1</td><td>1</td><td>10</td></tr><tr><td>B3Bin0Am0</td><td>10<sup>3</sup></td><td>1</td><td>1</td><td>10</td></tr><tr><td>B5Bin0Am0</td><td>10<sup>5</sup></td><td>1</td><td>1</td><td>10</td></tr><tr><td>B4Bin0Am1</td><td>10<sup>4</sup></td><td>1</td><td>10</td><td>10</td></tr><tr><td>B4Bin1Am0</td><td>10<sup>4</sup></td><td>10</td><td>1</td><td>10</td></tr><tr><td>B4Bin2Am0</td><td>10<sup>4</sup></td><td>10<sup>2</sup></td><td>1</td><td>10</td></tr><tr><td>B4Bin3Am0</td><td>10<sup>4</sup></td><td>10<sup>3</sup></td><td>1</td><td>10</td></tr><tr><td>B5Bin1Am0</td><td>10<sup>5</sup></td><td>10</td><td>1</td><td>10</td></tr><tr><td>B5Bin2Am0</td><td>10<sup>5</sup></td><td>10<sup>2</sup></td><td>1</td><td>10</td></tr><tr><td>B5Bin3Am0</td><td>10<sup>5</sup></td><td>10<sup>3</sup></td><td>1</td><td>10</td></tr><tr><td>B5Bin4Am0</td><td>10<sup>5</sup></td><td>10<sup>4</sup></td><td>1</td><td>10</td></tr><tr><td>B3Bin1Am0</td><td>10<sup>3</sup></td><td>10<sup>1</sup></td><td>1</td><td>10</td></tr><tr><td>B3Bin2Am0</td><td>10<sup>3</sup></td><td>10<sup>2</sup></td><td>1</td><td>10</td></tr><tr><td>R20FID</td><td>10<sup>4</sup></td><td>1</td><td>1</td><td>20</td></tr></tbody></table>

**Notes.** B4Bin0Am1 quantifies the influence of $\Lambda_{A,0}$ while B5Bin0Am0 and B3Bin0Am0 are the reference runs for $\beta_{\text{out}} = 10^5$ and $\beta_{\text{out}} = 10^3$. All the runs with $\text{Bin} \neq 0$ in their label explore the role of the initial value of $\beta$ at $R_{\text{int}}$. R20FID is the same simulation as the fiducial one, with a cavity that is twice as big. The bold font indicates the fiducial simulation values.

location, defined as the radius where the surface density equals half of its maximum value, is subject to a small variation of $\Delta 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 $\Sigma$, 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 \le 1.5$. We refer the reader to Appendix B for a quantitative discussion on this accumulation.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Surface density (left panel) and plasma beta parameter $\bar{\beta}$ (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 $\bar{\beta} \sim 1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

Figure 3 also pictures the evolution of $\bar{\beta}$, whose results are similar to the ones for $\Sigma$. Inside the cavity, $\bar{\beta}$ 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 $\bar{\beta}$ 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 $\Sigma$ and $\bar{\beta}$ reach a quasi-stationary state.

Gaps and rings are detected in the outer part of the disc, in the spatio-temporal diagram of both $\Sigma$ and $\beta$ (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.

## 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 $B_\varphi > 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 \approx \pm 0.3, \pm 0.6$, and $\pm 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 \ge 10$, while the depleted profile in $\rho$ 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 \ge 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

$$ \mathcal{L}_p = r \sin\theta \langle \rho u_p u_\varphi \rangle - r \sin\theta \langle B_p B_\varphi \rangle. \quad (24) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5. Time-averaged poloidal magnetic field lines and toroidal field component $\langle B_φ \rangle$ for the fiducial simulation. We note the peculiar field topology close to the truncation radius.</figcaption></figure>

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.

### 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 $\dot{M}$ defined as

$$ \dot{M}(R,t) = -2\pi R \overline{\rho v_r}. \qquad (25) $$

The height over which $\rho v_r$ is integrated has a direct influence on $\dot{M}$, 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 $\epsilon_{int} = \tan[(\theta_+ - \theta_-)/2]$. Results are presented in Fig. 8 for three values of $\epsilon_{int}$. For $\epsilon_{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

<figure><img src="image_7.png" /><figcaption>Fig. 6. Time-averaged streamlines and density for the fiducial simulation. We note the peculiar shape of the streamlines around the transition radius.</figcaption></figure>

higher above the disc allows us to cancel out this effect. Moving to $\epsilon_{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 $\Sigma$. 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 $\epsilon_{int} = 0.9$, defined by

$$ \langle v_{\text{acc.}} \rangle = \frac{\langle \dot{M} \rangle}{2\pi R \langle \Sigma \rangle}. \qquad (26) $$

This velocity profile exhibits a well-defined transition between subsonic accretion outside the cavity with $\langle v_{acc.} \rangle \sim 10^{-3} \langle c_s \rangle$ and transsonic accretion inside with $\langle v_{acc.} \rangle \sim \langle c_s \rangle$.

#### 3.3.2. Governing equations for accretion

Accretion theory can be understood as the secular evolution of $\dot{M}$ and $\Sigma$. In systems driven by MHD processes, these two quantities are usually supplemented by the magnetic field $B_z$ threading the
<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Accretion rates for different integration height scales with respect to the radius inside the disc. The higher $\rho v_r$ is integrated the closer to a constant value $\dot{M}$ is in the cavity. The average value inside the cavity (from $R=1$ to $R=10$) is $\dot{M} = 1.4 \pm 0.2 \times 10^{-7} M_{\odot} \text{yr}^{-1}$.</figcaption></figure>

disc. We apply the vertical integration procedure to the mass and angular-momentum conservation equations, which become

$$ \partial_t \Sigma - \frac{1}{2\pi r} \partial_r \dot{M} = -[\sin\theta \rho v_\theta]_{\theta_+}^{\theta_-} \quad (27) $$

<figure><img src="image_10.png" /><figcaption>Fig. 9. Accretion speed for $\epsilon_{int} = 0.9$ in units of local sound speed $c_s$. The profile exhibits a clear transition between subsonic and transsonic accretion that occurs where the edge of the cavity is located.</figcaption></figure>

$$ \partial_t (\overline{r \sin\theta \rho v_\varphi}) - \frac{1}{2\pi r} \dot{M} \partial_r (r^2 \sin^2\theta \tilde{\Omega}(r)) = -\frac{1}{r} \partial_r (r^2 W_{r\varphi}) - W_{\theta\varphi}. \quad (28) $$

We defined $W_{r\varphi}$ and $W_{\theta\varphi}$, respectively, as the radial and surface stresses by

$$ \left\{ \begin{aligned} W_{r\varphi} &\equiv \overline{\rho \sin\theta v_r v_\varphi} - \sin\theta \frac{B_r B_\varphi}{4\pi} \\ W_{\theta\varphi} &\equiv \left[ r \sin^2\theta \left( \rho v_\theta v_\varphi - \frac{B_\theta B_\varphi}{4\pi} \right) \right]_{\theta_+}^{\theta_-} \end{aligned} \right. \quad (29) $$

We remind the reader 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

$$ \partial_t B_{\theta,0} = \frac{1}{r} \partial_r (r \mathcal{E}_{\varphi,0}). \quad (30) $$

### 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 $\epsilon_{int} = 0.9$ with time-averaged quantities. The first information is that inside the cavity, the time derivative of $\Sigma$ is close to zero, meaning the simulation reaches a steady state up to $R \approx 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 $\zeta = \zeta_+ + \zeta_-$ (Lesur 2021b), where $\zeta_+$ and $\zeta_-$ are defined by

$$ \langle \zeta_{\pm} \rangle = \pm \frac{\langle \rho v_z \rangle (\theta_{\pm})}{\langle \Sigma \rangle \Omega_{\mathrm{K}}} = \pm \frac{\langle \rho v_r \cos\theta \rangle - \langle \rho v_{\theta} \sin\theta \rangle}{\langle \Sigma \rangle \Omega_{\mathrm{K}}} \quad (31) $$
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Mass conservation for $\varepsilon_{\text{int}} = 0.9$. The time derivative of $\Sigma$ 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 files that do not contain all the time steps computed by the code.</figcaption></figure>

and the corresponding quantities are time averaged. The signs of $\zeta_{\pm}$ are chosen accordingly so that a positive value of $\zeta_{\pm}$ corresponds to matter leaving the surface at $\theta_{\pm}$. Since $\zeta_{+}$ and $\zeta_{-}$ are pretty much symmetric with respect to the midplane, we focus on $\zeta$ only. The results are illustrated in Fig. 11 where both $\langle\zeta\rangle$ and $-\langle\zeta\rangle$ are shown. In order to compare with self-similar models (Lesur 2021b), we studied the values of $\langle\zeta\rangle$ at $z_0 = 6h$, which corresponds to $\varepsilon_{\text{int}} = 0.6$. The mass-loss-rate parameter is approximately constant in the external part of the disc around $6.2 \times 10^{-5}$, while it peaks at $2.9 \times 10^{-2}$ in the inner part. We find two zones where $\langle\zeta\rangle < 0$. One is close to the inner boundary and probably a boundary condition artefact, while the other extends from $R \approx 5$ to $R \approx 17$ au and is related to the material falling down on the disc around the transition zone, such a contribution being notably stronger for $\varepsilon_{\text{int}} = 0.6$.

To compare to self-similar solutions, we show the self-similar scaling of the mass-loss-rate parameter with respect to $\langle\beta\rangle$ derived by Lesur (2021b), which reads $\langle\zeta_{\text{self}}\rangle = 0.24 \langle\beta\rangle^{-0.69}$. It comes as no surprise that this fit does not account for negative values of $\langle\zeta\rangle$ 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 magnetosphere that compresses the disc magnetosphere, resulting in a deviation of $\zeta$ 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, $\zeta$ is significantly weaker than expected from a naive extrapolation of self-similar scaling laws. This indicates that the mass-loss rate saturates at $\beta \sim 1$, a regime which was not explored by Lesur (2021b).

An alternative model to the self-similar one is used to describe $\langle\zeta\rangle$ with greater accuracy. The self-similar fit is kept for the external parts of the disc $\zeta_{\text{ext}} = \zeta_{0,\text{ext}} \langle\beta\rangle^{d_{\text{ext}}}$, with $a_{\text{ext}} = -0.69$ and $\zeta_{0,\text{ext}} = 0.24$. Another one is then calculated for the inner

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** $\langle\zeta\rangle$ parameter for $\varepsilon_{\text{int}} = 0.6$. The self-similar fit shown here is for comparison only and was also obtained for $\varepsilon_{\text{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 $\langle\zeta\rangle$ in the whole disc.</figcaption></figure>

part only, $\zeta_{\text{int}} = \zeta_{0,\text{int}} \langle\beta\rangle^{d_{\text{int}}}$, with $a_{\text{int}} < 0$, so that the final profile is given by

$$ \langle\zeta\rangle_{\text{fit}} = \frac{\zeta_{0,\text{ext}} \langle\beta\rangle^{d_{\text{ext}}}}{1 + \frac{\zeta_{0,\text{ext}}}{\zeta_{0,\text{int}}} \langle\beta\rangle^{a_{\text{ext}} - a_{\text{int}}}}. \quad (32) $$

We obtain $a_{\text{int}} = -0.20$ and $\zeta_{0,\text{int}} = 0.018$. Such a model, with $a_{\text{ext}} - a_{\text{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 $\beta_t \approx 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.

### 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 better readability.

The integration height is $\varepsilon_{\text{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 \approx 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 $\alpha$ 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 $\alpha$. Nevertheless, the $\alpha$ parameter can still be used in this wind model, whose definition when time averaged is

$$ \langle\alpha\rangle = \frac{\langle W_{r\varphi}\rangle}{\langle\bar{P}\rangle}. \quad (33) $$
<table><caption>**Table 2.** Transport coefficients for a subset of simulations.</caption><thead><tr><th>Name</th><th>M (10<sup>-7</sup> M<sub>☉</sub> yr<sup>-1</sup>)</th><th>ζ<sub>in</sub> (10<sup>-2</sup>)</th><th>ζ<sub>ext</sub> (10<sup>-5</sup>)</th><th>α<sub>in</sub></th><th>α<sub>ext</sub> (10<sup>-3</sup>)</th><th>ν<sub>W,in</sub></th><th>ν<sub>W,ext</sub> (10<sup>-4</sup>)</th></tr></thead><tbody><tr><td>B4Bin0Am0</td><td>1.4</td><td>2.9</td><td>6.2</td><td>13</td><td>4.9</td><td>1.0</td><td>2.3</td></tr><tr><td>B3Bin0Am0</td><td>5.1</td><td>3.8</td><td>14</td><td>19</td><td>23</td><td>1.4</td><td>31</td></tr><tr><td>B5Bin0Am0</td><td>0.27</td><td>2.1</td><td>4.5</td><td>6.6</td><td>1.0</td><td>0.18</td><td>0.19</td></tr><tr><td>B4Bin0Am1</td><td>1.2</td><td>-1.8</td><td>10</td><td>2.8</td><td>16</td><td>0.16</td><td>2.7</td></tr><tr><td>R20FID</td><td>1.1</td><td>4.5</td><td>5.8</td><td>15</td><td>2.4</td><td>1.2</td><td>10</td></tr></tbody></table>

**Notes.** The accretion rate is calculated inside the cavity. The bold font indicates the fiducial simulation values.

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Angular momentum conservation multiplied by $r^{-3/2}$ and time averaged. Full blue line is $\langle \partial_r(r \sin\theta \rho v_\phi) \rangle r^{-3/2}$, red dot-dashed line is $\langle \partial_r(r^2 W_{r\phi}) \rangle r^{-3/2}$, green dashed line is $\langle W_{\theta\phi} \rangle r^{-3/2}$, and purple dotted line is $-\langle \frac{1}{2\pi r} \dot{M} \partial_r(r^2 \tilde{\Omega}) \rangle r^{-3/2}$.</figcaption></figure>

The corresponding profile is given in Fig. 13, where $\epsilon_{int} = 0.9$. In the external part of the disc, $\langle \alpha \rangle = 49 \pm 5 \times 10^{-4}$, while it reaches a maximum value inside the cavity $\langle \alpha \rangle = 13 \pm 5$.

Following a similar procedure as the one for $\alpha$, we define a dimensionless number associated with the surface stress component, $\nu_W$. As for $\zeta$, we define $\nu_{W,\pm}$, which were chosen to be positive for angular momentum leaving the disc on both sides:

$$ \langle \nu_W \rangle = \langle \nu_{W,+} \rangle + \langle \nu_{W,-} \rangle = \frac{\langle W_{\theta\phi} \rangle}{r \langle P_0 \rangle}. \quad (34) $$

We show the dependence of $\nu_W$ on $R$ in Fig. 13. In the external disc, $\langle \nu_W \rangle = 2.3 \pm 1.1 \times 10^{-4}$, while it rises up to $\langle \nu_W \rangle = 1.0 \pm 0.1$ inside the cavity. The same observations are drawn for both $\langle \alpha \rangle$ and $\langle \nu_W \rangle$ as for $\langle \zeta \rangle$. 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 $\beta \sim 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 $\alpha$ and $\nu_W$, which are both of the order of unity. Table 2 displays the transport coefficient values for all the simulations.

### 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

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Time-averaged transport coefficients $\langle \alpha \rangle$ and $\langle \nu_W \rangle$ for $\epsilon_{int} = 0.9$.</figcaption></figure>

<table><caption>**Table 3.** MHD invariants for a subset of simulations, computed with time-averaged quantities on the last 1000 orbits at $R_{int}$.</caption><thead><tr><th>Name</th><th>λ<sub>in</sub></th><th>λ<sub>ext</sub></th><th>κ<sub>in</sub></th><th>κ<sub>ext</sub></th><th>ω<sub>in</sub></th><th>ω<sub>ext</sub></th></tr></thead><tbody><tr><td>B4Bin0Am0</td><td>23</td><td>3.2</td><td>2.2 × 10<sup>-2</sup></td><td>2.5</td><td>0.67</td><td>1.2</td></tr><tr><td>B3Bin0Am0</td><td>185</td><td>1.3</td><td>3.1 × 10<sup>-3</sup></td><td>9.9</td><td>0.64</td><td>0.93</td></tr><tr><td>B5Bin0Am0</td><td>4.4</td><td>1.2</td><td>1.3</td><td>18</td><td>0.69</td><td>1.1</td></tr><tr><td>B4Bin0Am1</td><td>4.9</td><td>1.5</td><td>0.24</td><td>5.0</td><td>0.23</td><td>3.6</td></tr><tr><td>R20FID</td><td>26</td><td>2.1</td><td>1.8 × 10<sup>-2</sup></td><td>2.8</td><td>0.52</td><td>1.4</td></tr></tbody></table>

**Notes.** The bold font indicates the fiducial simulation values.

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 \equiv R_W$. The corresponding Keplerian angular velocity is $\Omega_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

$$ \kappa \equiv 4\pi \frac{\rho v_p \Omega_W R_W}{B_p B_W}. \quad (35) $$

The rotation parameter is given by

$$ \omega \equiv \frac{\Omega}{\Omega_W} - \frac{\kappa B_W B_p}{4\pi \rho R R_W \Omega_W^2}. \quad (36) $$
<figure><img src="image_16.png" /><figcaption>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.</figcaption></figure>

The magnetic lever arm accounts for the angular momentum that is removed from the disc by the wind is defined as

$$ \lambda \equiv \frac{\Omega R^2}{\Omega_w R_w^2} - \frac{R B_{\varphi}}{R_w B_w \kappa}. \qquad (37) $$

Of course, these invariants echo the transport coefficient definitions previously used to describe the disc, and one expects $\kappa \approx \beta\zeta/4\epsilon$ and $\lambda \approx 1 + \epsilon v_w/\zeta$ (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 $\Omega_{\kappa}$ 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 $\kappa_{in} = 2.2 \times 10^{-2}$, $\kappa_{ext} = 2.5$, $\lambda_{in} = 23$ and $\lambda_{ext} = 3.2$. These values are coherent with the transport coefficients computed in previous sections. We also note that the values of $\kappa$ and $\lambda$ 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).
### 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

$$ \dot{R}_0 = -\frac{1}{2\pi R_0} \frac{\delta \dot{M}}{\delta \Sigma}, \qquad (38) $$

where $\dot{R}_0$ is the cavity edge ‘velocity’ and $\delta \dot{M}$ and $\delta \Sigma$ are the jump in accretion rate and surface density at the cavity edge. By calculating $\dot{M}$ and $\Sigma$ around $R_0$, we find $\dot{R}_0 = 1.8 \times 10^{-5}$ while directly evaluating the cavity edge motion $\dot{R}_0$ yields $\dot{R}_0 = 1.4 \times 10^{-5}$ (both in c.u.). 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 $\psi$ such that

$$ \psi(r, t) = R_{\text{int}}^2 \int_{0}^{\pi/2} B_r(R_{\text{int}}, \theta, t) \sin\theta \, d\theta - \int_{r=R_{\text{int}}}^{r} r B_{\theta}(r, \pi/2, t) \, dr. \qquad (39) $$

Assuming the total flux is constant with respect to time, the iso-contours of $\psi$ describe the motion of the magnetic field lines in the disc plane. The spatio-temporal diagram for $\psi$ 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 $\langle B_{\theta,0} \rangle < 0$ in the transition region ($8 \le R \le 12$) and $\langle B_{\theta,0} \rangle > 0$ otherwise. This transition region is recovered in Fig. 15 as a region where $\partial_r \psi < 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 \approx 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 $\psi$ read

$$ \begin{cases} \partial_t \psi(R, t) = -R \mathcal{E}_{\varphi,0}(R, t) \\ \partial_R \psi(R, t) = -R B_{\theta,0}(R, t) \end{cases} \qquad (40) $$

Following Guilet & Ogilvie (2014), we rewrote these evolution equations as an advection equation for $\psi$:

$$ \partial_t \psi + v_\psi \partial_R \psi = 0, \qquad (41) $$

<figure><img src="image_17.png" /><figcaption>Fig. 15. Flux function $\psi(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 \approx 13$, the flux accumulates and exhibits a striped structure for smaller radii.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Magnetic-field transport parameter $v_B$ as a function of radius. We note that the outer disc is transporting magnetic field lines inwards ($v_B < 0$).</figcaption></figure>

where we defined the ‘field advection velocity’:

$$ v_{\psi} = -\frac{\mathcal{E}_{\varphi,0}(R, t)}{B_{\theta,0}(R, t)}. \qquad (42) $$

Eventually, we define a dimensionless advection parameter $v_B = v_\psi/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 $\langle v_B \rangle$ in Fig. 16. In the external disc, we find that the magnetic field is advected inwards with a velocity of $v_\psi = -2.6 \times 10^{-3} v_K$. $v_B$ changes its sign multiple times in the cavity, but it remains negative close to the cavity edge, between $R \approx 7$ and $R \approx 11$ au, where $v_\psi = +3.2 \times 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, $v_B$ is negative and $v_\psi = -2.6 \times 10^{-3} v_K$ so that vertical magnetic field pointing upwards is advected. We note that this inwards advection of the outer disc field lines is in sharp contrast to other work that focused
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Temporal evolution of $\Sigma$ in dotted green, $\dot{M}$ in dashed blue, and $B_z$ (vertically averaged) in black (full line) at $R = 3$ au for the fiducial simulation. $\zeta$ is calculated at $\varepsilon_{\text{int}} = 0.3$ and shown by a semi-dashed red line with a logarithmic scale on the right of the panels. Apart from $\zeta$, all the profiles are given in arbitrary units and divided by their maximum value reached during the timescale of the top panel.</figcaption></figure>

on ‘full’ discs (Bai & Stone 2017; Lesur 2021b). We revisit this discrepancy in the discussion.

### 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 $\Sigma$, $\dot{M}$, 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_{\text{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$, $\Sigma$, and $\dot{M}$ are correlated. When an inflow of matter crosses the cavity, $\Sigma$ peaks as well as $\dot{M}$, which in turn increases $\zeta$. In terms of temporal sequence, it seems that $B_z$ increases slightly before $\Sigma$ and $\dot{M}$, 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 $\zeta$ 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_{\varphi} = 0$ at the edge of the magnetic loop and $B_{\rho} = 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 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 $\zeta$ eventually peaks at 0.1. This explains why the effective value of $\langle\zeta\rangle$ 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
<figure><img src="image_20.png" /><figcaption>**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.</figcaption></figure>

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 $\bar{\beta} < \beta_{crit.} \approx 0.0355$. Figure 3 shows that $\bar{\beta}$ 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 $\beta$ 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 $\bar{\beta} < 0.1$, as it does during some of the bursts.
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Σ(R, t) and $\bar{\beta}(R, t)$ of B4Bin0Am1. The cavity stands throughout the simulation, though its edge falls down to R ≈ 2 au during the transient, before broadening up to R ≈ 4 au in a few thousands of orbits at R<sub>int</sub>. The profile of $\bar{\beta}$ is characterised by the presence of gaps in the external disc.</figcaption></figure>

# 4. Parameter space exploration

## 4.1. Ambipolar Diffusion

We checked the influence of $\Lambda_{A,0}$ in the simulation B4Bin0Am1. This simulation is the same as the fiducial one, except for the initial value of $\Lambda_A$, which is set to 10.

### 4.1.1. General structure of the disc and gaps

The spatio-temporal evolution of Σ and $\beta$ are shown in Fig. 19 for $\Lambda_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<sub>int</sub>. 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 $\bar{\beta}$ (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 $\beta$ 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<sub>int</sub>, 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 $\Lambda_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 \times 10^{-7} M_{\odot} \text{yr}^{-1}$ and an accretion velocity that is still subsonic in the outer disc and peaks up to $2 c_s$ at the internal radius. Therefore, the accretion picture is identical to the one for the fiducial run with an internal transsonic regime connecting through the cavity edge to a weakly magnetised wind.

Regarding the wind, we obtain a highly mass-loaded field line in the external disc that removes little angular momentum ($\lambda_{\text{ext}} = 1.5$ and $\kappa_{\text{ext}} = 5.0$) and a lighter one in the internal disc that carries a massive load of angular momentum ($\lambda_{\text{in}} = 4.9$ and $\kappa_{\text{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 $\beta_{\text{in}}$ and external $\beta_{\text{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 $\beta$ between $\beta_{\text{out}} = 10^3$ (run B3Bin0Am0) and $\beta_{\text{out}} = 10^5$ (run B5Bin0Am0).

With regard to B5Bin0Am0, the spatio-temporal evolutions of Σ and $\bar{\beta}$ 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<sub>int</sub>. 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<sub>int</sub>, assuming the first one is due to the initial transient) give the illusion that some gas might be created inside the cavity, challenging
<figure><img src="image_22.png" /><figcaption>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 $\langle B_{\phi} \rangle$.</figcaption></figure>

mass conservation. These bursts are actually due to gas accu- mulating 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 lim- itation of our model regarding the implementation of the inner boundary conditions, but they only occur in the weakly magne- tised ($\beta_{out} = 10^5$) simulations. Lastly, we add that estimating $\dot{R}_0$ 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 simula- tion B3Bin0Am0 quickly expands up to $R \approx 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 $\dot{R}_0 \approx 3.0 \times 10^{-5}$ c.u., which is about three times faster than the fiducial run. We obtain $\delta \dot{M} \approx 2.3 \times 10^{-4}$ and $\delta \Sigma \approx 4.2 \times 10^{-2}$ (both in code units), so Eq. (38) gives $\dot{R}_0 \approx 4.2 \times 10^{-5}$ c.u., where we chose $R_0 \approx 20$. The simple model we used seems to overesti- mate 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 $\beta_{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 exter- nal structure and converges to $\beta_{in} \lesssim 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 $\beta_{out}$. This can be understood using Eq. (38), which can be recast as

$$ \dot{R}_0 = v_{\text{acc., in}} \left( \frac{\Sigma_{\text{in}}}{\Sigma_{\text{out}}} - \frac{v_{\text{acc., out}}}{v_{\text{acc., in}}} \right), \quad (43) $$

where we defined the accretion velocities $v_{acc.} \equiv \dot{M}/2\pi R_0 \Sigma$, and assumed $\Sigma_{out} \gg \Sigma_{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 $\beta$. Writing $v_{acc.} \propto \beta^{-\sigma}$ with $\sigma > 0$, Lesur (2021b) proposes $\sigma = 0.78$ and Bai & Stone (2013) $\sigma = 0.66$, which indicates that $0 < \sigma < 1$. Assuming that a value $\tilde{\beta}$ exists for which $\dot{R}_0 = 0$, we obtain the following scaling:

$$ \dot{R}_0 = v_{\text{acc., in}} \frac{1}{\beta_{\text{out}}} \left[ 1 - \left( \frac{\beta_{\text{out}}}{\tilde{\beta}} \right)^{1-\sigma} \right], \quad (44) $$
<figure><img src="image_23.png" /><figcaption>**Fig. 21.** Spatio-temporal diagrams for $\Sigma$ and $\bar{\beta}$ 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).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 22.** Impact of initial external magnetisation on the surface density. We average the profile on the last 1000 orbits at $R_{\text{int}}$. For B5Bin0Am0, we average on 1000 orbits at $R_{\text{int}}$ occurring between the 2 burst events seen in Fig. 21.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 23.** Impact of internal initial magnetisation on the plasma parameter for $\beta_{\text{out}} = 10^4$.</figcaption></figure>

where we used the fact that $\Sigma_{\text{in}}/\Sigma_{\text{out}} = \beta_{\text{in}}/\beta_{\text{out}} = \beta_{\text{out}}^{-1}$ in our setup. The relation (44) shows that for $\beta_{\text{out}} \ll \bar{\beta}$, we have approximately $\dot{R}_0 \approx v_{\text{acc.}}$ in $\beta_{\text{out}}^{-1}$, indicating that the cavity expansion speed should increase as $\beta_{\text{out}}$ decreases, which is precisely what we observe for B3Bin0Am0. For $\beta_{\text{out}} \gg \bar{\beta}$, we obtain, on the contrary, $\dot{R}_0 \approx -v_{\text{acc.}}$ in $\beta^{\sigma-1}\beta_{\text{out}}^{-\sigma}$, showing a change of sign (hence a contraction of the cavity), albeit with a reduced speed. This regime might correspond to B5Bin0Am0, indicating that $\bar{\beta} \approx 10^4$.

### 4.2.2. Role of the internal initial plasma parameter

To study the impact of $\beta_{\text{in}}$, we ran a set of simulations that cover all the possible initial gaps $\beta_{\text{in}}/\beta_{\text{out}}$ where $\log\beta_{\text{out}} \in \{3, 4, 5\}$ and $\log\beta_{\text{in}} \in [0; \log\beta_{\text{out}}]$. We compared each result to the one obtained with $\beta_{\text{in}} = 1$ and the corresponding value of $\beta_{\text{out}}$. A striking result is the fact that the disc’s inner structure does not depend on $\beta_{\text{in}}$. No matter which $\beta_{\text{in}}$ we initially choose, a transition occurs in the cavity in order to impose $\beta_{\text{in}} \approx 1$. Interestingly, this threshold value is the one required to achieve transsonic accretion, which is mentioned in Wang & Goodman (2017). We illustrate this statement with Fig. 23 for the particular case of $\beta_{\text{out}} = 10^4$. We focus on the transient state of B4Bin3Am0 in
<figure><img src="image_26.png" /><figcaption>**Fig. 24.** Spatio-temporal diagrams of Σ (first panel), $B_z$, the vertical magnetic field at the midplane, and $\psi$, the flux function defined in Eq. (39). These profiles focus on the first orbits of the run B4Bin3Am0.</figcaption></figure>

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 $\beta$). 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 $\beta$ 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 $\beta_{in}$ equal to $\beta_{out}$ would simulate a full disc with no cavity. Hence, a threshold should exist regarding the value of $\beta_{in}$, above which no cavity is able to form. Considering Fig. 23, it seems that this threshold is $\ge 10^3$.

From these observations, we deduce that the cavity is regulated by the value of the plasma parameter, which must take a value close to 1. The reason for this regulation is not entirely clear, and we add a word of caution regarding the role of the inner radial boundary condition, especially with respect to the magnetic field transport at $R_{int}$. We discuss this influence in Appendix B. Referring to Sects. 3.5.3 and 3.5.4, we suggest that the RTI may be responsible for this regulation, but a dedicated study would be required to ascertain this claim.

### 4.3. Zoom with a larger cavity radius

We performed a simulation with a double-sized cavity ($R_0 = 20$ au) in order to check the impact of the cavity size. The simulation was integrated for 1000 orbits at $R = 10$ au so that it reaches 355 orbits at $R = 20$ au. The general observations are confirmed, such as the 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 $\beta_{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 ($\beta \approx 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 \approx 0.13$ G (see Eq. (17), with $\beta_{out} = 10^4$), so, initially, $B_z \approx 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 ($\lambda \ge 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
<figure><img src="image_27.png" /><figcaption>**Fig. 25.** Spatio-temporal diagrams for $\langle\Sigma\rangle$ and $\langle\beta\rangle$ for R20FID.</figcaption></figure>

accretion, implies that the kinematics of these cavities have sin-
gular 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 sim-
ulations; 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 inves-
tigated 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 sim-
ulations, 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 sup-
port 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 asso-
ciated with the formation of plasmons by reconnection events
(Ripperda et al. 2022). These features are recovered in our mod-
els of TDs, despite the fact that we used Newtonian dynamics
(MADs are usually found in GRMHD simulations) and the pres-
ence 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 sup-
presses 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 bound-
ary 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 diffu-
sion, and as such, subject to MRI quenching by the non-linearity
embedded in the ambipolar diffusivity ($\eta \propto 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 forth-
coming 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 Λ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<sub>X</sub> (see their Fig. 2, panels 2 and 3) as well as the role of
the temperature T<sub>0</sub> at 16 au (Fig. 2, panels 6 and 7). Regarding
our profile of ΛA ≈ 1–10, our work is similar to their models 2
(with L<sub>X</sub> = 10<sup>29</sup> erg s<sup>-1</sup>) and 6 (where T<sub>0</sub> = 30 K). Therefore,
we anticipate that an increase of two orders of magnitude for
L<sub>X</sub> would lead to ΛA > 10<sup>2</sup> 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 tem-
perature is less straightforward and seems to have little impact
on ΛA.

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 Λ<sub>A</sub>. 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).

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 mag-
netic 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 ≲ β<sub>int</sub> ≲ 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, Andrés 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 `TEX`/`LATEX` file prepared by the authors.

References

Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in *Protostars and Planets VI*, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475

Bai, X.-N. 2011, ApJ, 739, 50

Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

Bai, X.-N., & Stone, J. M. 2017, ApJ, 836, 46

Bajer, K., & Mizerski, K. 2013, Phys. Rev. Lett., 110, 104503

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Bi, J., & Fung, J. 2022, ApJ, 928, 74

Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Bouvier, J., Alencar, S. H. P., Boutelier, T., et al. 2007, A&A, 463, 1017

Bouvier, J., Alecian, E., Alencar, S. H. P., et al. 2020a, A&A, 643, A99

Bouvier, J., Perraut, K., Bouquin, J.-B. L., et al. 2020b, A&A, 636, A108

Carmona, A., Pinte, C., Thi, W. F., et al. 2014, A&A, 567, A51

Carmona, A., Thi, W. F., Kamp, I., et al. 2017, A&A, 598, A118

Casassus, S., Marino, S., Pérez, S., et al. 2015, ApJ, 811, 92

Ćemeljić, M., Shang, H., & Chiang, T.-Y. 2013, ApJ, 768, 5

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Combet, C., & Ferreira, J. 2008, A&A, 479, 481

Combet, C., Ferreira, J., & Casse, F. 2010, A&A, 519, A108

Cuello, N., Gonzalez, J.-F., & Pignatale, F. C. 2016, MNRAS, 458, 2140

Cui, C., & Bai, X.-N. 2021, MNRAS, 507, 1106

Dong, R., & Dawson, R. 2016, ApJ, 825, 77

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 497

Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659

Fang, M., Kim, J. S., van Boekel, R., et al. 2013, ApJS, 207, 5

Ferreira, J. 1997, A&A, 319, 340

Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18

Gammie, C. F. 1996, ApJ, 457, 355

Gárate, M., Delage, T. N., Stadler, J., et al. 2021, A&A, 655, A18

Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126

Guilet, J., & Ogilvie, G. 2014, MNRAS, 441, 852

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848

Jacquemin-Ide, J., Lesur, G., & Ferreira, J. 2021, A&A, 647, A192

Kane, Y., 1966, IEEE Trans. Antennas Propag., 14, 302

Lesur, G. R. J. 2021a, J. Plasma Phys., 87, 205870101

Lesur, G. R. J. 2021b, A&A, 650, A35

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023

Liska, M. T. P., Musoke, G., Tchekhovskoy, A., Porth, O., & Beloborodov, A. M. 2022, ApJ, 935, L1
Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805

Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18

Marsh, K. A., & Mahoney, M. J. 1992, ApJ, 395, L115

McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Mishra, B., Begelman, M. C., Armitage, P. J., & Simon, J. B. 2020, MNRAS, 492, 1855

Morishima, R. 2012, MNRAS, 420, 2851

Najita, J. R., Strom, S. E., & Muzerolle, J. 2007, MNRAS, 378, 369

Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Perez-Becker, D., & Chiang, E. 2011, ApJ, 735, 8

Pouilly, K., Bouvier, J., Alecian, E., et al. 2020, A&A, 642, A99

Riols, A., & Lesur, G. 2018, A&A, 617, A117

Riols, A., & Lesur, G. 2019, A&A, 625, A108

Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95

Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJ, 924, L32

Rosenfeld, K. A., Chiang, E., & Andrews, S. M. 2014, ApJ, 782, 62

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013, ApJ, 764, 66

Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117

Spruit, H. C., & Taam, R. E. 1990, A&A, 229, 475

Spruit, H. C., Stehle, R., & Papaloizou, J. C. B. 1995, MNRAS, 275, 1223

Stehle, R., & Spruit, H. C. 2001, MNRAS, 323, 587

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., Suzuki, T. K., & Shang, H. 2019, MNRAS, 484, 107

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74

Thi, W. F., Lesur, G., Woitke, P., et al. 2019, A&A, 632, A44

Umebayashi, T., & Nakano, T. 1980, PASJ, 32, 405

Umebayashi, T., & Nakano, T. 2008, ApJ, 690, 69

van der Marel, N., van Dishoeck, E. F., Bruderer, S., Pérez, L., & Isella, A. 2015, A&A, 579, A106

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

Van Rossum, G. 2020, The Python Library Reference, release 3.8.2 (Python Software Foundation)

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

Vlemmings, W. H. T., Lankhaar, B., Cazzoletti, P., et al. 2019, A&A, 624, A7

Wang, L., & Goodman, J. J. 2017, ApJ, 835, 59

Wardle, M. 2007, Astrophys. Space Sci., 311, 35

Wardle, M., & Koenig, A. 1993, ApJ, 410, 218

Zanni, C., & Ferreira, J. 2013, A&A, 550, A99

Zhang, K., Isella, A., Carpenter, J. M., & Blake, G. A. 2014, ApJ, 791, 42

Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34

Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47
## Appendix A: Ambipolar diffusivity for a transition disc - a simple model

The aim of this appendix is to model the ambipolar diffusivity spatial dependence in both a TD and a standard protoplanetary disc (i.e. without cavity). The general procedure to reach such a result follows and adapts the main calculation steps that are presented in Combet et al. (2010). As assumed in Eq. 4, only the ambipolar diffusivity does appear in the MHD equations, which we assume is the dominant non-ideal effect in the regime of discs we used at $R \ge 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

$$ \eta_A = \frac{B^2}{4\pi \gamma_{in} \rho_n \rho_i} \qquad (A.1) $$

(Wardle 2007), where $\rho_n$ and $\rho_i$ are, respectively, the density of the neutrals (the gas so $\rho_n = \rho$) and of the ions and $\gamma_{in} = \langle \sigma v \rangle_{in}/(m_n + m_i)$ with $\langle \sigma v \rangle_{in}$ is the ion-neutral collision rate whose value is

$$ \langle \sigma v \rangle_{in} = 2,0 \times 10^{-9} \left( \frac{m_H}{\mu} \right)^{1/2} \mathrm{cm}^3 \mathrm{s}^{-1} \qquad (A.2) $$

(Bai 2011), with $m_H$ being the atomic mass and $\mu = 2,34 m_H$ the mean molecular weight. Introducing the ionisation fraction $\xi = \rho_i/\rho_n$, one obtains

$$ \eta_A = 1,6 \times 10^{16} \left( \frac{\xi}{1 \times 10^{-13}} \right)^{-1} \left( \frac{B}{1G} \right)^2 \left( \frac{\rho}{1 \times 10^{14} \mathrm{cm}^{-1}} \right)^{-2} \mathrm{cm}^2 \mathrm{s}^{-1} \qquad (A.3) $$

Ambipolar diffusion is usually evaluated with the dimensionless ambipolar Elsasser number $\Lambda_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:

$$ m + \text{ionising radiation} \rightarrow m^+ + e^- \qquad \zeta_i, \qquad (A.4) $$

$$ m^+ + e^- \rightarrow m \qquad \delta, \qquad (A.5) $$

with $\zeta_i$ being the ionisation rate and $\delta$ the dissociative recombination rate. Following Fromang et al. (2002), we take

$$ \delta = 3 \times 10^{-6} T^{-1/2} \mathrm{cm}^3 \mathrm{s}^{-1}. \qquad (A.6) $$

In this toy model, we then have

$$ \xi = \sqrt{\frac{\zeta_i}{\delta \rho}} + \xi_{FUV} \qquad (A.7) $$

(Lesur et al. 2014), where $\xi_{FUV}$ accounts for the far-UV photon contribution that we modelled following Perez-Becker & Chiang (2011) as

$$ \xi_{FUV} = 2 \times 10^{-5} \exp \left[ - \left( \frac{\Sigma_{\star}}{0,03 \mathrm{g} \mathrm{cm}^{-2}} \right)^4 \right], \qquad (A.8) $$

with $\Sigma_{\star}$ the column density computed from the star to the point of interest.

To calculate $\zeta_i$, we considered several ionisation sources. We modelled the X-ray ionisation from the protostar by two

<figure><img src="image_28.png" /><figcaption>Fig. A.1. Ambipolar Elsasser number $\Lambda_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.</figcaption></figure>

bremsstrahlung-emitting corona (following Bai & Goodman 2009 and Igea & Glassgold 1999):

$$ \zeta_X = L_{X,29} \left(\frac{R}{1 \mathrm{au}}\right)^{-2,2} \left[ \zeta_1 \left( e^{-(N_{H1}/N_1)^\alpha} + e^{-(N_{H2}/N_1)^\alpha} \right) + \zeta_2 \left( e^{-(N_{H1}/N_2)^\beta} + e^{-(N_{H2}/N_2)^\beta} \right) \right], \qquad (A.9) $$

with $L_{X,29} \equiv L_X/10^{29} \mathrm{erg} \mathrm{s}^{-1}$ and $L_X$, $\zeta_1$, $\zeta_2$, $\alpha$, $\beta$, $N_1$, $N_2$ being the numerical values defined in Bai & Goodman (2009), while $N_{H1}$ and $N_{H2}$ are the column densities of hydrogen vertically computed above and below the calculation point.

We modelled the cosmic-ray ionisation rate following Umebayashi & Nakano (1980) so that:

$$ \zeta_{CR} = \zeta_{CR,0} e^{-\Sigma_{col.}/96 \mathrm{g} \mathrm{cm}^{-2} \mathrm{s}^{-1}}, \qquad (A.10) $$

where $\zeta_{CR,0} = 10^{-17} \mathrm{s}^{-1}$ and $\Sigma_{col.}$ is the matter column density above and below the point of interest.

Lastly, we added the radioactive decay that is assumed constant (Umebayashi & Nakano 2008) with a ionisation rate given by

$$ \zeta_{\mathrm{rad.}} = 10^{-19} \mathrm{s}^{-1}. \qquad (A.11) $$

Combining the Eqs. A.9, A.10, and A.11, we obtain $\zeta_i = \zeta_X + \zeta_{CR} + \zeta_{\mathrm{rad.}}$, paving the way to finally reach $\Lambda_A$ using Eqs. 13, A.6, and A.7. We note that due to the dependency of $\eta_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 $\Lambda_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 $\Lambda_A$ 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 $\Lambda_A$ in the discs remain quite close to unity in both cases, while the general trend of $\Lambda_A$ in a standard protoplanetary disc is recovered even in the case of a TD (Thi et al. 2019).
Moreover, $\Lambda_A$ remains fairly below the critical value $\Lambda_{A, \text{crit.}} = 10^2$ with or without a cavity. $\Lambda_A$ must stay below $\Lambda_{A, \text{crit.}}$ so that the MRI effects are negligible (Blaes & Balbus 1994; Bai 2011). Therefore, assuming a characteristic value of $\Lambda_{A,0} = 1$ captures, with a reasonable level of accuracy, the physics of ambipolar diffusion, and the cavity does not alter the ambipolar Elsasser number profile. The results we achieve from this simple toy model are to be compared to the more detailed work of Wang & Goodman (2017), in which many chemical species were taken into account to compute the ambipolar Elsasser number inside the cavity of a wind-driven TD.

Following Lesur (2021b) and Thi et al. (2019), we implemented the profile of $\Lambda_A$ so that

$$ \Lambda_A(z, R) = \Lambda_{A,0} \exp\left(\frac{z}{\lambda h}\right)^4, \qquad (A.12) $$

where $\lambda$ 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 $3h$. $\Lambda_{A,0}$ remains a free parameter (see Section 1 for more details). Additionally, a cut-off is used for the $\eta_A$ profile so that if $\eta_A > \eta_{A, \text{max}}$, the value of $\eta_A$ is replaced by $\eta_{A, \text{max}} \equiv 10 \epsilon^2$ in code units, such a choice being reflected in the $\Lambda_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 $\langle\Sigma\rangle_{4000}$ time averaged on the first 4000 orbits at the internal radius (when the differences are enhanced), with a focus on the innermost radii. We highlight that these differences do not increase for $t > 4000$ orbits at $R_{\text{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 $\Sigma$, 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 spatio-temporal diagram of the surface density on which the burst is clearly detected at 17435 orbits at $R_{\text{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 $\beta$ and blocks accretion (we recall that the accretion speed is $v_{\text{acc.}} \propto \beta^{-\sigma}$ with $\sigma > 0$). As a result, $\dot{M}$ falls from 0,25

down to $0,1 \times 10^{-7} M_{\odot} \text{yr}^{-1}$ in the region between the inner radial boundary and the burst, and matter piles up in the cavity. This episode ends when the magnetic flux is eventually 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 $\beta_{\text{out}} = 10^5$).

Regarding the inner radial boundary condition for the magnetic field, we tried several configurations (an outflow condition, which is the one we eventually chose, and a perfect conductor). Both of these conditions lead to the same steady states.

We also ran a simulation with a stronger magnetic field close to the inner boundary condition, but no significant changes were noticed. The additional magnetic field was chosen so that the magnetisation of the seed is set close to its saturation value in the fiducial run. However, in any case, the same transient state occurs and leaves the stage to a similar steady state (the magnetisation of the seed reaches the same saturation value and the same stripes are observed in the spatio-temporal diagram of $\psi$).

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.

## 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

$$ g_m \partial_R \ln \frac{\Sigma}{B_z} > 2 \left( r \frac{d\Omega}{dr} \right)^2 \equiv 2 S^2, \qquad (C.1) $$

where $S$ is the shear that we approximate with $S^2 = 9/4 \Omega^2$, and $g_m$ is

$$ g_m \equiv \frac{B_R^+ B_z}{2\pi \Sigma}. \qquad (C.2) $$

$B_R^+$ is the radial component of the magnetic field at the disc surface. Let us rewrite the previous expression in terms of $\beta$, $q$ (defined with $B_R^+ = q B_z$) and $\delta$ (defined as $\delta = -d \ln \Sigma / d \ln R$):

$$ \frac{B_R^+ B_z}{2\pi \Sigma} \partial_R \ln \frac{\Sigma}{B_z} = \frac{B_R^+ B_z \Sigma'}{2\pi \Sigma} - \frac{B_R^+ B_z B_z'}{2\pi \Sigma B_z}, \qquad (C.3) $$

$$ = \frac{q B_z^2}{2\pi \Sigma R} - \frac{\delta}{4\pi \Sigma} (B_z^2)', \qquad (C.4) $$
<figure><img src="image_29.png" /><figcaption>**Fig. B.1.** Surface density time-averaged on the first 4000 orbits at $R_{\text{int}}$. 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.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. B.2.** Spatio-temporal diagrams of $\Sigma$ (first panel), $B_{z,0}$ the vertical magnetic field at the midplane and $\psi$ 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 $\Sigma$.</figcaption></figure>

where $X'$ denotes the derivative of $X$ with respect to $R$. With $P = c_s^2 \rho = (h \Omega_K)^2 \Sigma / (\sqrt{2\pi} h)$, we obtain

$$ \beta = \frac{4 \sqrt{2\pi} R \epsilon \Omega_K^2 \Sigma}{B_z^2}. \qquad (C.5) $$

Therefore, the instability criterion becomes

$$ S^2 < -\frac{4 \epsilon \Omega_K^2 q \delta}{\sqrt{2\pi} \beta} - \frac{q}{4\pi \Sigma} \partial_R \left( \frac{4 \sqrt{2\pi} R \epsilon \Omega_K^2 \Sigma}{\beta} \right), \qquad (C.6) $$

$\epsilon$ being constant in the disc as well as $\beta$ inside the cavity. $\Omega_K$ varies as $R^{-3/2}$ and $\Sigma$ as $R^{-\delta}$ so that

$$ S^2 < \frac{4 \epsilon \Omega_K^2}{\sqrt{2\pi} \beta} q (-\delta + 1 + \delta/2). \qquad (C.7) $$

By taking $S^2/\Omega_K^2 = 9/4$, the RTI can be triggered when

$$ \beta < \frac{16 \epsilon}{9 \sqrt{2\pi}} q \left(1 - \frac{\delta}{2}\right). \qquad (C.8) $$

If we now assume that $\delta = q = 1$ for simplicity, we finally obtain

$$ \beta < \frac{8 \epsilon}{9 \sqrt{2\pi}} \approx 0.355 \epsilon = 0.0355 \equiv \beta_{\text{crit}}, \qquad (C.9) $$

where $\epsilon = 0.1$.

<figure><img src="image_31.png" /><figcaption>**Fig. C.1.** Interchange instability criteria. The red dotted line shows the critical value of $\beta_{\text{crit.}}$, while the green solid line is obtained with Eq. C.8.</figcaption></figure>

Figure C.1 compares the time-averaged values of $\bar{\beta}$ with the criterion given in Eq. C.9. The value of $\beta_{\text{crit.}}$ is below the time-averaged values of $\bar{\beta}$. 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

$$ \delta X \equiv X - \langle X \rangle. \qquad (D.1) $$

Focusing on $W_{r\varphi}$, we expand the magnetic term as

$$ \overline{\langle B_r B_\varphi \rangle} = \overline{\langle B_r \rangle \langle B_\varphi \rangle} + \overline{\langle \delta B_r \rangle \langle \delta B_\varphi \rangle}. \qquad (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:

$$ \langle W_{r\varphi}^{\text{lamin.}} \rangle = -\frac{1}{4\pi} \overline{\sin\theta \langle B_r \rangle \langle B_\varphi \rangle}, \qquad (D.3) $$

and for the laminar surface stress

$$ \langle W_{\theta\varphi}^{\text{lamin.}} \rangle = -r \left[ \sin^2 \theta \frac{\langle B_{\theta} \rangle \langle B_{\varphi} \rangle}{4\pi} \right]_{\theta_{+}}^{\theta_{-}}. \qquad (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

$$ \begin{cases} \langle \alpha^{\text{lamin.}} \rangle = \frac{\langle W_{r\varphi}^{\text{lamin.}} \rangle}{\langle \bar{P} \rangle}, \\ \langle \nu_{W}^{\text{lamin.}} \rangle = \frac{\langle W_{\theta\varphi}^{\text{lamin.}} \rangle}{r \langle P_0 \rangle} \end{cases} $$

while we define their turbulent counterparts as

$$ \begin{cases} \langle \alpha^{\text{turb.}} \rangle = \langle \alpha \rangle - \langle \alpha^{\text{lamin.}} \rangle \\ \langle \nu_{W}^{\text{turb.}} \rangle = \langle \nu_{W} \rangle - \langle \nu_{W}^{\text{lamin.}} \rangle \end{cases}. $$

The results are shown in Fig. D.1. The laminar contribution is the major one for $\langle \nu_W \rangle$ 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 $\langle \alpha \rangle$, a strong turbulent term is at stake, especially in the external part of the disc where it is dominant. Inside the cavity, $\langle \alpha \rangle$ is fairly distributed between the laminar and turbulent contributions. However, we recall that the wind may act also on the turbulent component of $\langle \alpha \rangle$ 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 $\langle \alpha \rangle$ coefficient, while the surface stress embodied by $\langle \nu_W \rangle$ is definitely dominated by the laminar part and due to the wind.

<figure><img src="image_32.png" /><figcaption>Fig. D.1. Time-averaged transport coefficients and their laminar and turbulent contributions. We give the laminar and turbulent contributions for $\langle \alpha \rangle$ and the total profile with its laminar contribution for $\langle \nu_W \rangle$.</figcaption></figure>