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

# Magnetic Braking of Accreting T Tauri Stars: Effects of Mass Accretion Rate, Rotation, and Dipolar Field Strength

Lewis G. Ireland¹, Claudio Zanni², Sean P. Matt¹, and George Pantolmos³

¹Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK; L.G.Ireland@exeter.ac.uk

²INAF—Osservatorio Astrofisico di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, Italy

³Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

Received 2020 July 27; revised 2020 November 2; accepted 2020 November 4; published 2020 December 29

## Abstract

The rotational evolution of an accreting pre-main-sequence star is influenced by its magnetic interaction with its surrounding circumstellar disk. Using the PLUTO code, we perform 2.5D magnetohydrodynamic, axisymmetric, time-dependent simulations of star-disk interaction—with an initial dipolar magnetic field structure, and a viscous and resistive accretion disk—in order to model the three mechanisms that contribute to the net stellar torque: accretion flow, stellar wind, and magnetospheric ejections (periodic inflation and reconnection events). We investigate how changes in the stellar magnetic field strength, rotation rate, and mass accretion rate (changing the initial disk density) affect the net stellar torque. All simulations are in a net spin-up regime. We fit semi-analytic functions for the three stellar torque contributions, allowing for the prediction of the net stellar torque for our parameter regime, as well as the possibility of investigating spin evolution using 1D stellar evolution codes. The presence of an accretion disk appears to increase the efficiency of stellar torques compared to isolated stars, for cases with outflow rates much smaller than accretion rates, because the star-disk interaction opens more of the stellar magnetic flux compared to that from isolated stars. In our parameter regime, a stellar wind with a mass-loss rate of ≈1% of the mass accretion rate is capable of extracting ≲50% of the accreting angular momentum. These simulations suggest that achieving spin equilibrium in a representative T Tauri case within our parameter regime, e.g., BP Tau, would require a wind mass-loss rate of ≈25% of the mass accretion rate.

Unified Astronomy Thesaurus concepts: Stellar winds (1636); Stellar evolution (1599); Stellar magnetic fields (1610); Stellar rotation (1629); Magnetohydrodynamical simulations (1966); Pre-main-sequence stars (1290); Stellar accretion (1578); Classical T Tauri stars (252); Stellar accretion disks (1579)

## 1. Introduction

Pre-main-sequence (PMS) solar-like stars, such as classical T Tauri stars (CTTS), exhibit rotational evolution that still lacks a comprehensive theoretical explanation. Many of these protostellar objects are observed to have a wide range of rotation rates, and a majority of these appear to be rotating at speeds much lower than their break-up velocity (see, e.g., Bouvier et al. 1993; Rebull et al. 2004; Herbst et al. 2007). CTTS are observed to still actively accrete material from a circumstellar accretion disk, exhibiting mass accretion rates of ~$10^{-10}-10^{-8} M_{\odot} \text{ yr}^{-1}$ (Edwards et al. 1994; Gullbring et al. 1998; Hartmann et al. 1998; Ingleby et al. 2014; Alcalá et al. 2017). Combined with the fact that these objects are still undergoing gravitational contraction, one would expect these protostars to be spinning up; however, observations have suggested solar-like stars in the PMS phase appear to have a fairly constant distribution of spin rates for a few million years (Bouvier et al. 1993; Edwards et al. 1993; Rebull et al. 2004; Irwin & Bouvier 2009; Gallet et al. 2019). Therefore, to explain the rotational evolution we observe, we require an additional mechanism that efficiently removes angular momentum from a star during its PMS phase.

CTTS are magnetically active objects, exhibiting ~kG strength multipolar magnetic fields (see, e.g., Johns-Krull 2007; Gregory et al. 2008; Johnstone et al. 2014; Donati et al. 2019, 2020), and it has been theorized that magnetic interaction between a star and its disk allows for the removal of angular momentum, allowing accreting stars to achieve slow rotation rates (see, e.g., Camenzind 1990; Koenigl 1991). The classical Ghosh & Lamb (1979) model proposes that a spin-down torque

is exerted onto the stellar surface from magnetic field lines connecting the star and the accretion disk beyond the corotation radius, where the disk is rotating more slowly than the stellar rotation rate. However, this assumption that the magnetic field is predominantly closed and connected to a large proportion of the disk is not consistent with our modern picture. It has been shown both analytically and via magnetohydrodynamical (MHD) simulations that the “twisting” of magnetic field lines as a result of differential rotation between the star and the disk inflates and opens the magnetosphere at mid-latitudes (see Matt & Pudritz 2005a, and references therein). Combined with the dilution of the poloidal field intensity (Agapitou & Papaloizou 2000; Zanni & Ferreira 2009), the efficiency of the spin-down torque associated with the extended magnetosphere in the Ghosh & Lamb (1979) picture is vastly reduced.

Other proposed solutions consider the presence of outflows, allowing angular momentum to be removed from the system. Magnetized stellar winds extract angular momentum from the stellar surface via open magnetic field lines, providing a spin-down torque contribution (Matt & Pudritz 2005b). Matt & Pudritz (2008a) show that, for spin equilibrium, where the spin-up accretion torque balances with the spin-down stellar wind torque, the stellar wind mass-loss rate is likely ~10% of the mass accretion rate. However, this is unlikely to be achievable if the outflows are predominantly powered by the star's rotational, thermal, and magnetic energies alone; therefore, an additional energy source that is not of thermal origin would be required (Matt & Pudritz 2007). The idea of accretion-powered stellar winds was proposed by Matt & Pudritz (2005b). However, possible mechanisms behind the power transfer have
only been explored in a few papers (Cranmer 2008, 2009), and the energy available has strict limits (Matt & Pudritz 2008a; Zanni & Ferreira 2011).

The opening and reconnecting of closed magnetospheric field lines, coined magnetospheric ejections (MEs), is another ejection mechanism; it is a result of the build-up of toroidal magnetic field pressure due to the star-disk differential rotation (Zanni & Ferreira 2013). Performing 2.5D axisymmetric star-disk interaction simulations, Zanni & Ferreira (2013) found that, similar to a magnetic slingshot, MEs are powered by the star-disk differential rotation and can extract angular momentum from both the disk and the star. In fact, they find it possible for a weak stellar wind torque, with a mass outflow rate ~1% of the mass accretion rate, to balance the spin-up accretion torque. In general, MEs take a net amount of angular momentum (from the disk) away from the system, reducing the accretion torque. As they are magnetically connected (intermittently) to both the star and the disk, a fraction of the extracted angular momentum can be given to the star. This happens when the disk is truncated close enough to the star such that the disk's orbital motion is faster than the stellar rotation. In that case, the magnetic field lines anchored in the disk rotate faster than the star, making them twist forward, with respect to the stellar rotation, at the stellar surface (but they twist backward at the disk surface, with respect to the disk's rotation). This imparts a spin-up torque on the star (and spin-down on the disk). In this case, the MEs can also add mass back onto the star; most of the stellar ME mass flux comes from the disk and can be sensibly larger than 1% of the mass accretion rate (see Zanni & Ferreira 2013). Alternatively, when the disk is truncated far enough from the star such that the inner disk's angular rotation is slower than the star, the magnetic field lines twist backward, with respect to the stellar rotation, at the stellar surface, imparting a spin-down torque on the star. Therefore, the ME stellar torque contribution is primarily a function of whether the MEs anchor to the disk in a place where the disk rotates faster or slower than the star. Thus, MEs may be important in the spin evolution of PMS stars.

The relevance of outflows powered by magnetospheric accretion, either of stellar or disk origin, on the spin evolution of CTTS seems to be supported by the fact that, on average, nonaccreting or faintly accreting weak-lined T Tauri stars (WTTS) rotate faster than CTTS (see, e.g., Venuti et al. 2017). While the presence of strong magnetospheric winds in CTTS is supported by the observation of P-Cygni line profiles (see, e.g., Edwards et al. 2006) that can be reproduced either by stellar wind or inner disk wind models (Kurosawa et al. 2011), there are no reliable detections of winds in WTTS.

We attempt to characterize the angular momentum problem over a parameter regime that may represent stars over a wide range of ages, and to further explore the ME mechanism and their interaction with stellar winds, using SDI simulations to extract torque prescriptions that can be used in secular spin-evolution models. Over the last decade, stellar wind torque prescriptions have been proposed (see, e.g., Matt & Pudritz 2008b; Ud-Doula et al. 2009; Matt et al. 2012; Réville et al. 2015; Finley & Matt 2017, 2018; Pantolmos & Matt 2017) using isolated stellar wind (ISW) MHD simulations. The first spin-evolution calculation to use stellar torque prescriptions derived from SDI simulations was Gallet et al. (2019), using simulations from Zanni & Ferreira (2013), which they adopted

in their spin evolution models. They found that a strong stellar wind with a mass outflow rate corresponding to ≈10% of the accretion rate, or a disk with a truncation radius close to its corotation radius, was required for an efficient spin-down torque configuration. They propose an accretion torque formulation that accounts for a fraction of the analytical value due to the extraction of angular momentum from the disk via MEs, and a stellar ME torque prescription that depends on the torque exerted by magnetic field lines coupled to the MEs, as a result of the differential rotation between the star and the accretion disk. However, their prescription was only based on four simulations, and employs coefficients that appear to be guided by these simulations, rather than explicitly fitted. Furthermore, they use a stellar wind torque prescription calibrated using ISW simulations only (Matt et al. 2012).

In this paper, we make further improvements to this torque prescription, to better represent a wide parameter regime of SDI simulations and capture the interplay between stellar winds and SDI, in a dynamically self-consistent way. We introduce a new stellar wind torque formulation that encapsulates the effects of SDI, illustrating that the stellar wind (open) magnetic flux is strongly dependent on disk properties, predominantly the truncation radius. We propose an accretion torque formulation that depends on how the accretion disk deviates from Keplerian rotation for each simulation. We also propose an improved stellar ME torque prescription that again considers how strongly sub-Keplerian the disk is, but also accounts for the perturbation of the magnetosphere from simple dipolar scaling due to the existence of an accretion disk.

We perform 38 2.5D axisymmetric SDI simulations, changing the initial disk density, the rotation rate, and the magnetic field strength, in order to capture how the net stellar torque, due to the stellar wind, accretion, and MEs, depends on both stellar and accretion disk properties. Our simulations also adopt an equation of state with a variable polytropic index, allowing us to simultaneously treat the cooler disk gas with a more physically meaningful $\gamma = 5/3$, while also allowing for a Parker stellar wind with $\gamma = 1.05$. We therefore model the stellar winds as thermally driven outflows, by setting suitable density and temperature profiles. Despite its simplicity, this approach has two main caveats. First, it has been shown that, in order to be compatible with the emission properties and the accretion power available in CTTS, a thermally driven wind must possess a very low mass-loss rate (Matt & Pudritz 2007). From this point of view, our driving mechanism is most likely unrealistic, but since our goal is to evaluate the extraction of angular momentum via the stellar wind, the torque does not much depend on the precise mechanism that drives the wind. Second, this approach does not allow us to self-consistently model the energetic coupling between the wind and the accretion flow, which is usually invoked to produce wind mass-loss rates that are a substantial fraction of the accretion rate (Matt & Pudritz 2005b; Cranmer 2008, 2009; Zanni & Ferreira 2011). In our simulations, the mass-loss rate and wind power can be set independently of the accretion luminosity, by adjusting the wind boundary conditions. Our setup results in wind mass-loss rates typically on the order of 1% of the accretion rate, which is consistent with the mass flux predicted by models of stellar winds powered by the dissipation of MHD turbulence and waves excited by the inhomogeneous impact of the accretion funnels (Cranmer 2008, 2009).
In Section 2, we introduce the numerical simulation setup, the initial and boundary conditions, normalizations, and the simulation parameter study. In Section 3, we list our SDI simulations, present a qualitative description of their behavior and how it depends on the parameters in our study, and introduce the global quantities used in our investigation, such as the mass flow rate, torque, and unsigned magnetic flux. In Section 4, we introduce and fit semi-analytic formulations for the truncation radius and the torques from accretion, MEs, and stellar wind, which allows for the net stellar torque to be calculated. In Section 5, we discuss our findings and summarize our conclusions.

## 2. Numerical Method

### 2.1. Numerical Setup

The simulations in this paper are solved numerically using the following MHD equations, including resistive and viscous effects:

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) &= 0, \\
\frac{\partial (\rho \mathbf{v})}{\partial t} + \nabla \cdot \left[ \rho \mathbf{v} \mathbf{v} + \left( p + \frac{\mathbf{B} \cdot \mathbf{B}}{8\pi} \right) \mathbf{I} - \frac{\mathbf{B}\mathbf{B}}{4\pi} - \tau \right] &= \rho \mathbf{g}, \\
\frac{\partial E}{\partial t} + \nabla \cdot \left[ \left( E + p + \frac{\mathbf{B} \cdot \mathbf{B}}{8\pi} \right) \mathbf{v} - \frac{(\mathbf{v} \cdot \mathbf{B})\mathbf{B}}{4\pi} \right] &= \rho \mathbf{g} \cdot \mathbf{v} + \nabla \tau \cdot \mathbf{v} - \frac{\mathbf{B}}{4\pi} \cdot (\nabla \times \eta_m \mathbf{J}), \\
\frac{\partial \mathbf{B}}{\partial t} + \nabla \times (\mathbf{B} \times \mathbf{v} + \eta_m \mathbf{J}) &= 0.
\end{aligned}
\quad (1) $$

The expressions in Equation (1) represent the mass-continuity, momentum, energy, and magnetic induction equation, respectively. Here, $\rho$ is the mass density, $\mathbf{v}$ is the plasma flow velocity, $p$ is the plasma thermal pressure, $\mathbf{B}$ is the magnetic field, $\mathbf{I}$ is the identity matrix, $\mathbf{g} = -(GM_*/R^2)\hat{\mathbf{R}}$ is the gravitational acceleration (where $G$ is the gravitational constant, $M_*$ is the stellar mass, and $R$ is the spherical radius), $\mathbf{J} = \nabla \times \mathbf{B}/4\pi$ is the electric current, and $\eta_m$ is the magnetic resistivity. The magnetic diffusivity is defined as $\nu_m = \eta_m/4\pi$. The total energy density $E$ is defined as

$$ E = \rho u + \rho \frac{\mathbf{v} \cdot \mathbf{v}}{2} + \frac{\mathbf{B} \cdot \mathbf{B}}{8\pi}, \quad (2) $$

where $u(T)$ is specific internal energy as a function of temperature $T$, determined via our new adopted caloric equation of state for a calorically imperfect gas, where the specific heats and adiabatic index $\gamma$ are temperature-dependent; this is described in detail in the appendix of Pantolmos et al. (2020). The hot stellar wind plasma will behave near-isothermally ($\gamma = 1.05$), and the cold accretion disk will behave adiabatically ($\gamma = 5/3$). The viscous stress tensor $\tau$ is defined as

$$ \tau = \eta_v \left[ (\nabla \mathbf{v}) + (\nabla \mathbf{v})^T - \frac{2}{3} (\nabla \cdot \mathbf{v}) \mathbf{I} \right], \quad (3) $$

where $\eta_v$ is the dynamic viscosity. The kinematic viscosity is defined as $\nu_v = \eta_v/\rho$. The viscous and magnetic diffusive terms on the right-hand side of the total energy equation (in the

system of Equation (1)) correspond to the work of the viscous forces and diffusion of the magnetic energy. Dissipative viscous and ohmic heating terms are omitted from the total energy equation to avoid potential runaway heating of the accretion disk, as a result of the lack of cooling radiative effects in our equation (Zanni & Ferreira 2009).

We also solve the following scalar equations:

$$ 
\begin{aligned}
\frac{\partial \rho s}{\partial t} + \nabla \cdot (\rho s \mathbf{v}) &= 0, \\
\frac{\partial \rho T r}{\partial t} + \nabla \cdot (\rho T r \mathbf{v}) &= 0,
\end{aligned}
\quad (4) $$

where $s$ is the specific entropy (defined using the new equation of state), and $T$ is a passive tracer. The specific entropy can be used to monitor dissipation and heating from numerical integration of the total energy equation. $T r$ is used to distinguish between disk material and the stellar wind plasma.

To numerically solve Equations (1)–(4), we employ a second-order Godunov method from the PLUTO code<fnref n="4" /> (Mignone et al. 2007). The MHD equations are solved in the rotating frame of the star. We achieve spatial reconstruction of primitive variables using a combination of linear and parabolic interpolation. To compute intercell fluxes, we use the approximate HLLD Riemann solver (Miyoshi et al. 2010), computing only the deviation from the initial magnetic field. To advance the MHD equations in time, we employ a second-order Runge-Kutta scheme. We use the hyperbolic divergence cleaning method (Dedner et al. 2002) to control the divergence-free constraint $\nabla \cdot \mathbf{B} = 0$. A second-order finite difference approximation, integrated in time explicitly, is adopted to compute viscous and resistive terms.

The simulations in this paper are 2.5D (2D domain with 3D vector components), and they use the $(R, \theta)$ spherical coordinate system. We assume axisymmetry about the stellar rotation axis. Hereafter, we define the spherical radius using $R$ and the cylindrical radius using $r = R \sin \theta$. We use a logarithmic grid in the $R$ direction, covering [1.50.756521]$R_*$ across 320 grid cells, where $R_*$ is the stellar radius, which ensures that the cell sizes satisfy $\Delta R \sim R\Delta\theta$. We use a uniform grid in the $\theta$ direction, covering [0, $\pi$] across 256 grid cells.

### 2.2. Initial Conditions

The computational domain is initialized with a stellar corona, a viscous accretion disk, and a magnetic field configuration.

The stellar corona is initialized using density and pressure profiles from a 1D spherically symmetric, isentropic, transonic Parker wind solution. This solution is defined by its base density $\rho_*$ and sound speed at the stellar surface $c_{s,*}$. The poloidal speed is set to zero.

The Keplerian accretion disk is initialized adopting an $\alpha$ parameterization Shakura & Sunyaev (1973) for the viscosity. Neglecting inertial terms, the thermal disk pressure $p_d$ and density $\rho_d$ are determined via vertical hydrostatic equilibrium, and the disk toroidal velocity $v_{\phi,d}$ is derived via radial equilibrium. Assuming the thermal disk pressure and density are related via the polytropic relation $p_d \propto \rho_d^\alpha$, where $\gamma = 5/3$,

<footnote n="4">⁴ PLUTO is freely available to download at http://plutocode.ph.unito.it.</footnote>
we write

$$
\begin{align}
\rho_d &= \rho_{d,\star} \left\{ \frac{2}{5\epsilon^2} \left[ \frac{R_*}{R} - \left(1 - \frac{5\epsilon^2}{2}\right) \frac{R_*}{r} \right] \right\}^{3/2}, \nonumber \\
p_d &= \rho_{d,\star} \epsilon^2 v_{K,\star}^2 \left( \frac{\rho_d}{\rho_{d,\star}} \right)^{5/3}, \nonumber \\
v_{\phi,d} &= \left(1 - \frac{5\epsilon^2}{2}\right)^{1/2} \left( \frac{GM_*}{r} \right)^{1/2}, \tag{5}
\end{align}
$$

where $\epsilon = c_{s,d}/v_K|_{\theta=\pi/2}$ is the disk aspect ratio evaluated at the disk midplane, $c_{s,d} = (p_d/\rho_d)^{1/2}$ is the disk isothermal sound speed, $v_K = (GM_*/r)^{1/2}$ is the Keplerian velocity, and $\rho_{d,\star}$ and $v_{K,\star}$ are the disk density and Keplerian velocity, respectively, evaluated on the midplane at $R_*$. The kinematic viscosity of the disk, using an $\alpha$ parameterization (Shakura & Sunyaev 1973), is expressed as

$$ \nu_v = \frac{2}{3} \alpha_v \frac{c_{s,d}^2}{\Omega_K}, \quad (6) $$

where $\Omega_K = (GM_*/r^3)^{1/2}$ is the Keplerian angular velocity, and the magnetic diffusivity of the disk is expressed as a function of the kinematic viscosity:

$$ \nu_m = \frac{3}{2} \alpha_m \frac{\nu_v}{\alpha_v}, \quad (7) $$

where $\alpha_m$ is the anomalous diffusivity coefficient. Thus, the magnetic Prandtl number $Pr_m = \nu_v/\nu_m \equiv 2\alpha_v/3\alpha_m$. The radial accretion velocity is expressed as

$$ v_{R,d} = -\alpha_v \frac{c_{s,d}^2}{v_K} \sin \theta. \quad (8) $$

In this expression, we neglect the $\tau_{\theta\phi}$ component of the stress tensor in Equation (3), which avoids the likely unphysical backflow along the disk midplane, usually associated with 3D models of $\alpha$ accretion disks (see, e.g., Regev & Gitelman 2002; Zanni & Ferreira 2013). Equation (8) demonstrates that the inertial term due to radial accretion flow in the radial momentum equation is of order $\mathcal{O}(\alpha_v \epsilon^4)$, whereas the thermal pressure gradient is of order $\mathcal{O}(\epsilon^2)$; therefore, the inertial term can be neglected when deriving the disk equilibrium solution (Equation (5)). The passive tracer $Tr$ is initialized to one inside the disk and zero in the stellar corona.

Even though we use an $\alpha$ model to compute our initial disk structure, the thermal and density profiles of our initial disk are different from a Shakura & Sunyaev $\alpha$ disk or the illuminated $\alpha$ disk modeled by D’Alessio et al. (1998), since we made a simplifying adiabatic assumption to compute the disk thermal structure. Our disk is intended to be only a dynamical model that provides a source of mass and angular momentum accreting at a suitable rate.

The magnetospheric field is initially a purely dipolar configuration, which is aligned to the rotation axis. The $R$ and $\theta$ magnetic field components are

$$ B_R(R, \theta) = B_* \left(\frac{R_*}{R}\right)^3 \cos \theta \quad (9) $$

and

$$ B_{\theta}(R, \theta) = \frac{1}{2} B_{*} \left( \frac{R_{*}}{R} \right)^{3} \sin \theta, \quad (10) $$

respectively, where $B_*$ is the stellar surface polar magnetic field strength.

The disk is initially truncated, such that it is excluded from the region where the coronal magnetic stresses exceed the disk thermal pressure. We determine an initial truncation radius $R_{t,i}$, which corresponds to the radius at which the following implicit equation in the $R$ direction is true:

$$ m_s = \left| \frac{B_{\phi}^{+} B_{d,\theta=\pi/2}}{2\pi p_{d,\theta=\pi/2}} \right|, \quad (11) $$

where $m_s$ is the sonic Mach number of the accretion flow due to the large-scale magnetic field torque, and $B_{\phi}^+$ is the theoretical toroidal field at the disk surface due to star-disk differential rotation, estimated as

$$ B_{\phi}^{+} = \frac{B_{d,\theta=\pi/2}}{\alpha_m \epsilon} \left[ \left( \frac{R}{R_{\text{co}}} \right)^{3/2} - 1 \right], \quad (12) $$

where $R_{\text{co}} = (GM_*/\Omega_*^2)^{1/3}$ is the Keplerian corotation radius. Here, the disk’s Keplerian rotation equals the stellar angular velocity $\Omega_*$, and $B_{d,\theta=\pi/2} = B_*(R/R_*)^{-3}$ and $p_{d,\theta=\pi/2}$, respectively, are the initial magnetic field strength and the disk thermal pressure, both at the disk midplane. We set $m_s = 1.5$, resulting in transonic accretion at the inner disk radius; a typical condition for accretion funnels to form (see, e.g., Bessolaz et al. 2008). The upper limit is set as $R_{t,i} < 0.8R_{\text{co}}$. Once $R_{t,i}$ is determined, we set the disk pressure and density equal to zero (but initialize with the stellar corona solution) inside that radius.

Due to star-disk interaction, the initial disk structure is strongly modified; for $R < R_{t,i}$, we use the local sound speed to define the kinematic viscosity (Equation (6)) and magnetic diffusivity (Equation (7)), which is smoothly matched to the initial disk sound speed for $R > R_{t,i}$, where disk structure is largely unaffected. We multiply the viscosity and diffusivity by the passive scalar $Tr$, so that only the accretion disk is viscous and resistive. For $B^2/8\pi p > 1$, the viscosity and diffusivity are set to zero, to ensure ideal MHD conditions in regions of strong magnetization, such as the accretion funnels and stellar winds. Magnetic surfaces that initially thread the disk are imposed to rotate at Keplerian angular velocity, calculated at the disk anchoring radius, while those within $R_{t,i}$ are imposed to corotate with the star.

## 2.3. Boundary Conditions

The computational domain is enclosed by a boundary and augmented with ghost cells in order to allow boundary conditions to be implemented. An “axisymmetric” inner and outer boundary for the $\theta$ coordinate is employed, which symmetrizes variables and flips the signs of both normal and $\phi$ components of vector fields across the boundary. At the inner boundary for the $R$ coordinate, there are two conditions: for subsonic flow into the computational domain (out from the stellar boundary, i.e., stellar wind), we impose fixed density and thermal pressure profiles from the Parker wind solution in the ghost cells; for supersonic flow from the computational
domain (into the stellar boundary, i.e., accretion columns), density and thermal pressure are free to adjust to the values determined by the accretion funnel. In the latter, a power-law extrapolation along magnetic field lines is employed for the density and the thermal pressure is set assuming constant entropy along magnetic surfaces. Intermediate regions, i.e., subsonic outflow or a hydrostatic corona, use the sonic Mach number calculated in the first row of cells in the domain to linearly interpolate between the two boundary conditions. The R magnetic field component is kept fixed to conserve the total stellar flux, whereas the θ component is free to change using a linear extrapolation. We also impose the poloidal velocity v_p to be parallel to the poloidal magnetic field, along with the continuity of the axisymmetric MHD invariant k = 4πρv_p/B_p along field lines, to ensure smooth inflow of stellar wind and supersonic infall via accretion funnels onto the stellar surface, without any shock generation. Following Zanni & Ferreira (2009), we model the stellar surface as a perfect conductor that rotates at an angular velocity Ω_*, and we impose a boundary condition on B_φ that forces magnetic field lines to rotate at the stellar rotation rate. The radial derivative of B_φ, which is used to do the boundary extrapolation, is derived by replacing the local acceleration in the angular momentum equation with the following:

$$ \rho \frac{\partial r \Omega}{\partial t} = \rho \frac{r \Omega_* + v_p B_{\phi} / B_p - r \Omega}{\Delta t}, \quad (13) $$

where Δt is the Alfvén crossing time of one grid cell close to the stellar surface. The stellar rotation at the inner boundary is set as

$$ v_{\phi} = r\Omega + v_p \frac{B_{\phi}}{B_p}. \quad (14) $$

We impose a power-law extrapolation for density and thermal pressure, and a linear extrapolation for all other variables, at the outer boundary for the R coordinate. For the toroidal magnetic field along open field lines anchored to the stellar surface, we impose a boundary condition similar to the approach adopted at the inner boundary, in order to avoid artificial torques being exerted on the star due to sub-Alfvénic matter crossing the outer boundary.

## 2.4. Units and Normalization

Simulations in this paper are performed in dimensionless units. Here, we list normalization factors required to express our results in physical units. Length is expressed in units of the stellar radius R*, density is expressed in units of the initial coronal density at the stellar surface ρ*, and velocities are expressed in units of Keplerian velocity at the stellar surface v_{K,*} = (GM_*/R_*)^{1/2}, given the stellar mass M*. Hence, time is expressed in units of t_0 = R_*/v_{K,*}. Magnetic field strength is expressed in units of B_0 = (4πρ_* v_{K,*}^2)^{1/2}. Mass accretion/ejection rates are expressed in units of Ṁ_0 = ρ_* R_*^2 v_{K,*}, and torque is expressed in units of J̇_0 = ρ_* R_*^3 v_{K,*}^2.

In order to make direct comparisons with young stars, we adopt R* = 2 R_⊙, M* = 0.5 M_⊙, and ρ* = 10⁻¹² g cm⁻³, and

express the normalizations as

$$ \begin{align} v_{K,*} &= 218.38 \left( \frac{M_*}{0.5 M_\odot} \right)^{1/2} \left( \frac{R_*}{2 R_\odot} \right)^{-1/2} \text{km s}^{-1} \nonumber \\ B_0 &= 77.41 \left( \frac{\rho_*}{10^{-12} \text{ g cm}^{-3}} \right)^{1/2} \left( \frac{M_*}{0.5 M_\odot} \right)^{1/2} \left( \frac{R_*}{2 R_\odot} \right)^{-1/2} \text{G} \nonumber \\ t_0 &= 0.074 \left( \frac{M_*}{0.5 M_\odot} \right)^{-1/2} \left( \frac{R_*}{2 R_\odot} \right)^{3/2} \text{days} \nonumber \\ \dot{M}_0 &= 6.71 \times 10^{-9} \left( \frac{\rho_*}{10^{-12} \text{ g cm}^{-3}} \right) \left( \frac{M_*}{0.5 M_\odot} \right)^{1/2} \left( \frac{R_*}{2 R_\odot} \right)^{3/2} \text{M}_\odot \text{yr}^{-1} \nonumber \\ \dot{J}_0 &= 1.29 \times 10^{36} \left( \frac{\rho_*}{10^{-12} \text{ g cm}^{-3}} \right) \left( \frac{M_*}{0.5 M_\odot} \right) \left( \frac{R_*}{2 R_\odot} \right)^2 \text{erg.} \tag{15} \end{align} $$

## 2.5. Simulation Parameters

In our parameter study, we vary the following input parameters:

1. The initial disk density, ρ_{d,*};

2. Surface polar magnetic field strength B*, controlled by the parameter v_A/v_{esc}, i.e., the ratio of the surface polar Alfvén velocity v_A = B_*/(4πρ_*)^{1/2} and the stellar escape velocity v_{esc} = (2GM_*/R_*)^{1/2};

3. Stellar break-up fraction f = Ω_*R_*/v_{K,*}, where Ω_* is the stellar rotation rate.

The stellar rotation period can be written as

$$ P_* = 9.30 \left(\frac{f}{0.05}\right)^{-1} \left(\frac{M_*}{0.5 M_\odot}\right)^{-1/2} \left(\frac{R_*}{2 R_\odot}\right)^{3/2} \text{days.} \quad (16) $$

We fix the disk thermal aspect ratio at ϵ = 0.075, the viscous and resistive transport coefficients at α_v = 0.2 and α_m = 0.2, respectively, and the stellar wind sound speed at the stellar surface at c_{s,*} = 0.35v_{K,*}. The initial accretion rate of the disk, determined solely by the viscous torque, is

$$ \dot{M}_{\text{acc},i} \approx 0.12 \left(\frac{\alpha_{v}}{0.2}\right) \left(\frac{\rho_{d,*}}{100}\right) \left(\frac{\epsilon}{0.075}\right)^3 \dot{M}_0. \quad (17) $$

## 3. Simulations of Star-Disk Interaction

We run 38 simulations in total, changing the simulation parameters described in Section 2.5. Each simulation is run for a period that corresponds to roughly 20 stellar rotation periods for a break-up fraction f = 0.05 (t_0 = 2513.2742), therefore fewer stellar rotation periods are covered as f decreases. However, we still find sensible time averages over a quasi-steady state for our lowest f simulations. Typical simulation times in our study correspond to a period on the order of weeks or months, where changes in the stellar rotation rate would be negligible, thus using a fixed f for each simulation is justified. The input parameters and outputted variables (in simulation units) for all simulations can be found in Table 1; for convenience, normalization factors for converting these parameters to physical units are shown in Table 2 (see Section 2.4). Negative mass fluxes or torques represent an inflow of mass or angular momentum onto the star, whereas positive quantities represent an outflow from the star.
<table><caption>Table 1
Models and Variable Input Parameters (Left of Vertical Line), and Outputed Global Variables (Right)</caption><thead><tr><th>Model</th><th>f</th><th>ρ<sub>l,★</sub><br>ρ<sub>★</sub></th><th>v<sub>A</sub><br>v<sub>esc</sub></th><th>B<sub>★</sub><br>B<sub>0</sub></th><th>M<sub>wind</sub><br>M<sub>0</sub><br>[10<sup>-2</sup>]</th><th>M<sub>ME,★</sub><br>M<sub>0</sub><br>[10<sup>-2</sup>]</th><th>M<sub>acc</sub><br>M<sub>0</sub></th><th>J<sub>wind</sub><br>J<sub>0</sub></th><th>J<sub>ME,★</sub><br>J<sub>0</sub></th><th>J<sub>acc</sub><br>J<sub>0</sub></th><th>Φ<sub>wind</sub><br>Φ<sub>0</sub></th><th>Φ<sub>ME,★</sub><br>Φ<sub>0</sub></th><th>Φ<sub>acc</sub><br>Φ<sub>0</sub></th><th>⟨r<sub>A</sub>⟩<br>R<sub>★</sub></th><th>R<sub>l</sub><br>R<sub>★</sub></th><th>Υ<sub>wind</sub><br>[10<sup>4</sup>]</th><th>Υ<sub>acc</sub></th></tr></thead><tbody><tr><td>1</td><td>0.001</td><td>100</td><td>4.60</td><td>12.8</td><td>0.308</td><td>-0.131</td><td>-0.327</td><td>0.0381</td><td>-0.194</td><td>-0.438</td><td>6.77</td><td>4.81</td><td>69.6</td><td>35.04</td><td>5.47</td><td>1.05</td><td>352</td></tr><tr><td>2</td><td>0.001</td><td>200</td><td>4.60</td><td>12.8</td><td>0.401</td><td>-0.744</td><td>-0.627</td><td>0.0643</td><td>-0.329</td><td>-0.814</td><td>9.17</td><td>5.99</td><td>66.0</td><td>40.0</td><td>4.46</td><td>1.48</td><td>184</td></tr><tr><td>3</td><td>0.001</td><td>300</td><td>4.60</td><td>12.8</td><td>0.438</td><td>-0.952</td><td>-1.03</td><td>0.0946</td><td>-0.541</td><td>-1.21</td><td>11.7</td><td>7.03</td><td>62.4</td><td>46.5</td><td>3.78</td><td>2.23</td><td>112</td></tr><tr><td>4</td><td>0.001</td><td>400</td><td>4.60</td><td>12.8</td><td>0.456</td><td>-1.61</td><td>-1.48</td><td>0.128</td><td>-0.730</td><td>-1.68</td><td>14.1</td><td>7.99</td><td>59.1</td><td>52.6</td><td>3.33</td><td>3.06</td><td>78.1</td></tr><tr><td>5</td><td>0.0025</td><td>100</td><td>4.60</td><td>12.8</td><td>0.368</td><td>-0.177</td><td>-0.334</td><td>0.126</td><td>-0.211</td><td>-0.451</td><td>7.92</td><td>3.85</td><td>69.4</td><td>37.0</td><td>5.42</td><td>1.21</td><td>345</td></tr><tr><td>6</td><td>0.0025</td><td>200</td><td>4.60</td><td>12.8</td><td>0.462</td><td>-0.841</td><td>-0.625</td><td>0.206</td><td>-0.304</td><td>-0.820</td><td>10.4</td><td>4.66</td><td>66.1</td><td>42.1</td><td>4.47</td><td>1.67</td><td>184</td></tr><tr><td>7</td><td>0.0025</td><td>300</td><td>4.60</td><td>12.8</td><td>0.495</td><td>-0.986</td><td>-1.09</td><td>0.315</td><td>-0.526</td><td>-1.29</td><td>13.5</td><td>5.65</td><td>62.0</td><td>50.4</td><td>3.73</td><td>2.60</td><td>106</td></tr><tr><td>8</td><td>0.0025</td><td>400</td><td>4.60</td><td>12.8</td><td>0.517</td><td>-1.96</td><td>-1.46</td><td>0.398</td><td>-0.729</td><td>-1.67</td><td>15.8</td><td>6.37</td><td>59.0</td><td>55.5</td><td>3.33</td><td>3.43</td><td>78.7</td></tr><tr><td>9</td><td>0.005</td><td>100</td><td>4.60</td><td>12.8</td><td>0.395</td><td>-0.160</td><td>-0.38</td><td>0.297</td><td>-0.214</td><td>-0.498</td><td>8.44</td><td>3.64</td><td>69.1</td><td>38.8</td><td>5.21</td><td>1.27</td><td>303</td></tr><tr><td>10</td><td>0.005</td><td>200</td><td>4.60</td><td>12.8</td><td>0.497</td><td>-0.628</td><td>-0.711</td><td>0.51</td><td>-0.326</td><td>-0.909</td><td>11.5</td><td>4.31</td><td>65.4</td><td>45.3</td><td>4.31</td><td>1.89</td><td>162</td></tr><tr><td>11</td><td>0.005</td><td>300</td><td>2.30</td><td>6.38</td><td>0.521</td><td>-4.55</td><td>-1.30</td><td>0.711</td><td>-0.643</td><td>-1.26</td><td>14.3</td><td>4.49</td><td>21.8</td><td>52.3</td><td>2.12</td><td>2.77</td><td>22.1</td></tr><tr><td>12</td><td>0.005</td><td>300</td><td>4.60</td><td>12.8</td><td>0.533</td><td>-1.11</td><td>-1.11</td><td>0.737</td><td>-0.510</td><td>-1.33</td><td>14.5</td><td>5.00</td><td>61.7</td><td>52.6</td><td>3.72</td><td>2.77</td><td>103</td></tr><tr><td>13</td><td>0.005</td><td>300</td><td>9.20</td><td>25.5</td><td>0.365</td><td>-0.300</td><td>-1.06</td><td>0.579</td><td>-0.569</td><td>-1.58</td><td>15.1</td><td>6.44</td><td>141</td><td>56.4</td><td>5.89</td><td>4.42</td><td>434</td></tr><tr><td>14</td><td>0.005</td><td>400</td><td>4.60</td><td>12.8</td><td>0.559</td><td>-1.90</td><td>-1.56</td><td>0.939</td><td>-0.697</td><td>-1.79</td><td>17.2</td><td>5.51</td><td>58.5</td><td>58.0</td><td>3.29</td><td>3.73</td><td>73.8</td></tr><tr><td>15</td><td>0.01</td><td>100</td><td>4.60</td><td>12.8</td><td>0.425</td><td>-0.434</td><td>-0.380</td><td>0.663</td><td>-0.197</td><td>-0.533</td><td>8.99</td><td>3.44</td><td>68.8</td><td>39.5</td><td>5.19</td><td>1.34</td><td>303</td></tr><tr><td>16</td><td>0.01</td><td>200</td><td>2.30</td><td>6.38</td><td>0.516</td><td>-2.46</td><td>-0.810</td><td>1.10</td><td>-0.407</td><td>-0.827</td><td>11.5</td><td>3.34</td><td>25.8</td><td>46.1</td><td>2.54</td><td>1.80</td><td>35.5</td></tr><tr><td>17</td><td>0.01</td><td>200</td><td>4.60</td><td>12.8</td><td>0.525</td><td>-0.597</td><td>-0.763</td><td>1.12</td><td>-0.327</td><td>-0.971</td><td>12.2</td><td>4.19</td><td>64.8</td><td>46.3</td><td>4.23</td><td>2.01</td><td>151</td></tr><tr><td>18</td><td>0.01</td><td>200</td><td>9.20</td><td>25.5</td><td>0.325</td><td>-0.231</td><td>-0.791</td><td>1.17</td><td>-0.419</td><td>-1.28</td><td>13.5</td><td>5.55</td><td>143</td><td>59.9</td><td>6.40</td><td>3.93</td><td>582</td></tr><tr><td>19</td><td>0.01</td><td>300</td><td>4.60</td><td>12.8</td><td>0.559</td><td>-1.22</td><td>-1.21</td><td>1.58</td><td>-0.507</td><td>-1.45</td><td>15.3</td><td>4.73</td><td>61.1</td><td>53.2</td><td>3.61</td><td>2.97</td><td>95.5</td></tr><tr><td>20</td><td>0.01</td><td>400</td><td>4.60</td><td>12.8</td><td>0.580</td><td>-2.20</td><td>-1.61</td><td>2.09</td><td>-0.651</td><td>-1.84</td><td>17.8</td><td>5.18</td><td>58.2</td><td>60.0</td><td>3.25</td><td>3.87</td><td>71.5</td></tr><tr><td>21</td><td>0.025</td><td>100</td><td>4.60</td><td>12.8</td><td>0.535</td><td>-0.302</td><td>-0.583</td><td>1.96</td><td>-0.210</td><td>-0.797</td><td>11.1</td><td>3.53</td><td>66.6</td><td>38.3</td><td>4.64</td><td>1.62</td><td>198</td></tr><tr><td>22</td><td>0.025</td><td>200</td><td>4.60</td><td>12.8</td><td>0.610</td><td>-0.985</td><td>-1.09</td><td>3.18</td><td>-0.409</td><td>-1.32</td><td>14.9</td><td>4.39</td><td>61.9</td><td>45.5</td><td>3.74</td><td>2.56</td><td>105</td></tr><tr><td>23</td><td>0.025</td><td>300</td><td>4.60</td><td>12.8</td><td>0.647</td><td>-1.55</td><td>-1.51</td><td>4.27</td><td>-0.547</td><td>-1.79</td><td>17.8</td><td>4.75</td><td>58.7</td><td>51.4</td><td>3.34</td><td>3.45</td><td>76.5</td></tr><tr><td>24</td><td>0.025</td><td>300</td><td>9.20</td><td>25.5</td><td>0.477</td><td>-0.759</td><td>-1.80</td><td>4.49</td><td>-0.724</td><td>-2.61</td><td>19.9</td><td>7.19</td><td>135</td><td>61.3</td><td>5.01</td><td>5.85</td><td>256</td></tr><tr><td>25</td><td>0.025</td><td>400</td><td>3.45</td><td>9.57</td><td>0.669</td><td>-5.32</td><td>-2.02</td><td>5.22</td><td>-0.795</td><td>-2.16</td><td>20.2</td><td>4.67</td><td>36.1</td><td>55.9</td><td>2.41</td><td>4.29</td><td>32.1</td></tr><tr><td>26</td><td>0.025</td><td>400</td><td>4.60</td><td>12.8</td><td>0.676</td><td>-2.46</td><td>-1.91</td><td>5.44</td><td>-0.687</td><td>-2.22</td><td>20.4</td><td>5.09</td><td>55.7</td><td>56.7</td><td>3.03</td><td>4.33</td><td>60.2</td></tr><tr><td>27</td><td>0.05</td><td>100</td><td>2.30</td><td>6.38</td><td>0.822</td><td>-1.40</td><td>-0.656</td><td>3.56</td><td>-0.216</td><td>-0.712</td><td>11.5</td><td>2.46</td><td>26.6</td><td>29.4</td><td>2.73</td><td>1.14</td><td>43.9</td></tr><tr><td>28</td><td>0.05</td><td>100</td><td>3.45</td><td>9.57</td><td>0.814</td><td>-0.802</td><td>-0.787</td><td>4.49</td><td>-0.206</td><td>-0.973</td><td>13.1</td><td>3.10</td><td>44.7</td><td>33.2</td><td>3.41</td><td>1.49</td><td>82.3</td></tr><tr><td>29</td><td>0.05</td><td>100</td><td>4.60</td><td>12.8</td><td>0.774</td><td>-0.499</td><td>-0.955</td><td>5.50</td><td>-0.244</td><td>-1.28</td><td>15.0</td><td>3.91</td><td>62.2</td><td>37.7</td><td>3.87</td><td>2.07</td><td>121</td></tr><tr><td>30</td><td>0.05</td><td>200</td><td>3.45</td><td>9.57</td><td>0.848</td><td>-2.15</td><td>-1.34</td><td>6.68</td><td>-0.401</td><td>-1.50</td><td>16.8</td><td>3.67</td><td>40.5</td><td>39.7</td><td>2.83</td><td>2.34</td><td>48.4</td></tr><tr><td>31</td><td>0.05</td><td>200</td><td>4.60</td><td>12.8</td><td>0.837</td><td>-1.27</td><td>-1.55</td><td>7.83</td><td>-0.436</td><td>-1.88</td><td>18.5</td><td>4.48</td><td>58.2</td><td>43.3</td><td>3.32</td><td>2.88</td><td>74.4</td></tr><tr><td>32</td><td>0.05</td><td>300</td><td>3.45</td><td>9.57</td><td>0.897</td><td>-4.76</td><td>-2.01</td><td>9.15</td><td>-0.666</td><td>-2.20</td><td>20.7</td><td>4.11</td><td>36.0</td><td>45.2</td><td>2.41</td><td>3.39</td><td>32.2</td></tr><tr><td>33</td><td>0.05</td><td>300</td><td>4.60</td><td>12.8</td><td>0.880</td><td>-2.28</td><td>-2.13</td><td>9.97</td><td>-0.625</td><td>-2.45</td><td>21.6</td><td>4.81</td><td>54.8</td><td>47.5</td><td>2.96</td><td>3.75</td><td>54.1</td></tr><tr><td>34</td><td>0.05</td><td>400</td><td>3.45</td><td>9.57</td><td>0.943</td><td>-8.19</td><td>-3.21</td><td>12.1</td><td>-1.11</td><td>-3.30</td><td>24.9</td><td>4.97</td><td>31.0</td><td>50.7</td><td>2.09</td><td>4.64</td><td>20.2</td></tr><tr><td>35</td><td>0.05</td><td>400</td><td>4.60</td><td>12.8</td><td>0.931</td><td>-4.34</td><td>-2.96</td><td>12.6</td><td>-0.867</td><td>-3.32</td><td>25.4</td><td>5.45</td><td>50.4</td><td>52.0</td><td>2.60</td><td>4.89</td><td>39.0</td></tr><tr><td>36</td><td>0.0625</td><td>200</td><td>4.60</td><td>12.8</td><td>0.999</td><td>-1.82</td><td>-1.79</td><td>10.4</td><td>-0.447</td><td>-2.17</td><td>20.2</td><td>4.79</td><td>56.2</td><td>40.9</td><td>3.12</td><td>2.87</td><td>64.5</td></tr><tr><td>37</td><td>0.0625</td><td>300</td><td>4.60</td><td>12.8</td><td>1.02</td><td>-2.49</td><td>-2.35</td><td>12.6</td><td>-0.629</td><td>-2.71</td><td>22.7</td><td>5.03</td><td>53.4</td><td>44.4</td><td>2.84</td><td>3.56</td><td>48.9</td></tr><tr><td>38</td><td>0.0625</td><td>400</td><td>4.60</td><td>12.8</td><td>1.09</td><td>-5.56</td><td>-3.40</td><td>16.0</td><td>-0.868</td><td>-3.79</td><td>26.8</td><td>6.00</td><td>48.3</td><td>48.3</td><td>2.46</td><td>4.66</td><td>33.9</td></tr></tbody></table>

Note.

<footnote>a $B_{★}$ is not a fundamental input parameter, but it is simply derived from $v_A/v_{esc}$ and tabulated here for convenience.</footnote>
<table><caption>**Table 2**
Normalization of Each Parameter Type in Table 1</caption><thead><tr><td>Parameter</td><td>Normalization</td></tr></thead><tbody><tr><td>B<sub>0</sub></td><td>(4πρ<sub>*</sub>v<sup>2</sup><sub>K,*</sub>)<sup>1/2</sup></td></tr><tr><td>M&#x0307;<sub>0</sub></td><td>ρ<sub>*</sub>R<sup>2</sup><sub>vK,*</sub></td></tr><tr><td>j&#x0307;<sub>0</sub></td><td>ρ<sub>*</sub>R<sup>3</sup><sub>vK,*</sub></td></tr><tr><td>Φ<sub>0</sub></td><td>(4πρ<sub>*</sub>R<sup>4</sup><sub>vK,*</sub>)<sup>1/2</sup></td></tr></tbody></table>

### 3.1. Morphology of Star–Disk Interaction Systems

In this section, we qualitatively analyze the properties of the SDI system, focusing on the dynamical processes that are a result of this interaction, namely, the stellar wind, the accretion flow, and the time-dependent MEs. Zanni & Ferreira (2013) have a detailed analysis of this latter phenomenon, which we briefly describe below. Due to the time-dependent nature of these simulations, they can only achieve a quasi-steady or quasi-periodic state; however, as the focus of this paper is to investigate global properties and the long-term evolution of such a system, any periodic transient effects have been neglected by using time averages.

Figure 1 shows the density distribution for the computational domain of a representative SDI simulation (model 17), taken at a time corresponding to three stellar rotation periods (middle of time-averaged domain). The top panel of Figure 2 also shows a zoomed-in region of this Figure 1. The magnetic field (white) is strongly coupled to the accretion disk. We dissect the system into three distinct regions. Region (1) shows the ejection of the stellar wind occurring along these open magnetic field lines anchored to the stellar surface, providing the star with a spin-down torque contribution.

Region (2) contains the MEs, exhibiting the inflation and reconnection of magnetospheric field lines. These occur due to the differential rotation between the star and the disk, resulting in a strong toroidal field that inflates the magnetosphere and opens field lines close to the truncation radius (labeled $R_t$). This causes ejections to detach from the magnetosphere and propagate outward in between the open magnetic field lines and those connected to the disk, bounded by the red field lines (a dotted red line is used to connect asymmetric ME regions about the equatorial plane, and should not be interpreted as a magnetic field line).

At the inner disk edge, the magnetic field lines connect the disk to the star, and act as accretion funnels that allow the exchange of angular momentum between the disk and the stellar surface to occur at the truncation radius, which is shown in region (3). Region (3) also contains the “dead zone” region, which consists of corotating coronal gas below the accretion flow, around the magnetic equator. The incoming angular momentum from accretion provides the star with a spin-up torque contribution. Although our system is initialized with north–south symmetry (other than the magnetic polarity reversal, characteristic of a dipole field), our simulation setup does not assume or enforce any north–south symmetry during the evolution. We observe some models to establish funnel flows that are approximately symmetric, some to be dominated in the north or south hemispheric flows, and some that switch between regimes one or more times during the simulation. This can also produce hemispheric asymmetry in the MEs. For our example in Figure 1, the dominant accretion channel occurs in the northern hemisphere. The determination of mass-loss rates,

torques, and unsigned magnetic fluxes associated with each of these mechanisms is discussed in Section 3.2.

Our criterion for defining the boundary between regions (1) and (2) selects the magnetic field lines closest to the midplane (one for each hemisphere) that is attached to the star and reaches the outer domain, and where both $v_R > 0$ (ejection) and $B_R B_\phi < 0$ (spin-down torque). These field lines also correspond to the boundary between region (2) and the outer disk domain. Our criterion for defining the boundary between regions (2) and (3) selects the outermost closed magnetic field line that has a maximum R coordinate no more than 5% larger than its midplane radius in at least one of the hemispheres (in Figure 1, this criteria is met in the upper hemisphere, where the dominant accretion channel occurs), and where flow is supersonic at its connection with the stellar surface. This avoids misclassifying field lines that are still closed but are strongly deformed due to MEs (and in practice, these strongly deformed field lines open and close intermittently, while the lines we select typically remain closed for the duration). The truncation radius is then taken to be the radial coordinate corresponding to the minimum specific entropy (which locates the densest region of the disk material) along this last closed field line (Pantolmos et al. 2020).

We plot the Keplerian corotation radius $R_{co} = f^{-2/3}R_*$, which represents the radius at which a Keplerian disk rotates at the same speed as the star. For such a disk, if $r < R_{co}$, the magnetic

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Logarithmic density color map showing the entire domain of model 17, at $t = 3P_*$, where $P_*$ is the stellar rotation period. Magnetic field lines (white) and velocity vectors normalized by $v_{K,*}$ (green arrows) are included. We illustrate the three regions that contribute to the net stellar torque: (1) the stellar wind, (2) MEs, and (3) accretion flow. Region (2) is bounded by the red field lines, and regions (2) and (3) are separated by the cyan field line. Vertical white dashed lines (from left to right) show the position of the truncation radius, $R_t$, the magnetic corotation radius, $R_{co,m}$, and the Keplerian corotation radius, $R_{co}$, respectively. A zoomed-in region of this plot is shown in Figure 2.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Logarithmic density color map showing the inner domain of (a) model 17 (fiducial), (b) model 20 (twice the disk density), (c) model 16 (half the magnetic field strength), and (d) model 31 (five times the rotation rate). For each simulation, magnetic field lines (white) are included. Stellar wind and ME regions are separated by red field lines. ME and accretion regions are separated by cyan field line. Vertical white dashed lines (from left to right) show the position of the truncation radius, $R_t$, and the magnetic corotation radius, $R_{co,m}$, respectively.</figcaption></figure>

field is twisted such that field lines “lead” the stellar rotation, thus magnetic torques threading this region spin up the star; conversely, magnetic torques threading regions where $r > R_{co}$ spin down the star. However, MEs extract angular momentum from the accretion disk, resulting in both a smaller accretion torque being exerted onto the stellar surface and the disk rotating slower than Keplerian speed. Thus, in general, for $r < R_{co}$, the star–disk differential rotation is reduced, resulting in a weaker

magnetic “twisting” effect and reduced magnetic torques, relative to Keplerian rotation. Thus, we plot the magnetic corotation radius, $R_{co,m}$, defined such that the magnetic twist becomes zero for $R_t = R_{co,m}$. The technical criteria of how we determine $R_{co,m}$ in our simulations are discussed in Section 4.2.

The value of $R_t$ must be less than $R_{co}$ for disk material to lose angular momentum and accrete onto the star (Matt & Pudritz 2005a); for a disk truncated outside $R_{co}$, i.e., in the
“propeller” regime, accretion is not possible (see, e.g., Illarionov & Sunyaev 1975) or intermittent (see, e.g., Romanova et al. 2005; D’Angelo & Spruit 2010). How the MEs exchange angular momentum with the star depends on the star–disk differential rotation in the truncation region. For a system where $R_t \gtrsim R_{co,m}$, the MEs extract angular momentum from a region of the disk rotating slower than the star, and therefore the MEs exert a spin-down torque contribution on the star. However, all of our simulations reside in the regime where $R_t \lesssim R_{co,m}$, so the MEs extract angular momentum from a region of the disk rotating faster than the star, and we observe that the MEs exert a time-averaged spin-up torque contribution on the star. The region $R_{co,m} \lesssim R_t \lesssim R_{co}$ corresponds to a situation where it is possible that the MEs exert a spin-down torque on the stellar surface while the disk is still sub-Keplerian, and thus in principle can still be steadily accreting.

In Figure 2, we plot a density distribution of the inner domain for four SDI models, taken at a time corresponding to the middle of their time-averaged domain, to demonstrate the qualitative effects of changing the initial disk density, $\rho_{d,\star}$, the magnetic field strength, $B_*$, and the stellar break-up fraction, $f$. These include (a) our representative model (model 17), (b) a disk with double the initial disk density (model 20), (c) a star with half the magnetic field strength (model 16), and (d) a star rotating five times faster (model 31). All three latter cases essentially increase the ratio between the disk pressure and magnetic stress, which decreases $R_t$ compared to our fiducial model. Generally, increasing $\rho_{d,\star}$ increases the mass accretion rate and hence the disk pressure; increasing $f$ decreases the differential rotation between the star and the disk, thus decreasing the magnetic stress felt by the disk; and decreasing $B_*$ directly decreases the magnetic stress. We find that $R_{co,m}$ decreases with $f$, but also has a small inverse proportionality on $R_t$; we demonstrate this quantitatively in Section 4.2. For all cases, we include a zoom panel at identical coordinates, to demonstrate the relative positions of the field lines where accretion and ME mechanisms deposit angular momentum (cyan and red lines, respectively). For cases (b)–(d), a smaller truncation radius results in these field lines connecting to the star at lower latitudes, opening up a larger area on the stellar surface for the stellar wind to be ejected.

## 3.2. Global Flow Quantities and Efficiencies

The main aim of this paper is to investigate how the stellar wind, accretion flow, and MEs interact with the star. Thus, we investigate global flow properties, such as the mass flow rate, torque, and unsigned magnetic flux. Defining a surface S perpendicular to the poloidal flow at a given radius, these quantities are written as

$$ \dot{M} = \int_S \rho v_p \cdot dS = 2\pi R^2 \int_{\theta_1}^{\theta_2} \rho v_R \sin\theta \, d\theta, \quad (18) $$

$$ \begin{aligned} j &= \int_S r \left( \rho v_{\phi} v_p - \frac{B_{\phi} B_p}{4\pi} \right) \cdot dS \\ &= 2\pi R^3 \int_{\theta_1}^{\theta_2} \left( \rho v_{\phi} v_R - \frac{B_{\phi} B_R}{4\pi} \right) \sin^2 \theta \, d\theta, \end{aligned} \quad (19) $$

and

$$ \Phi = \int_S |\mathbf{B} \cdot dS| = 2\pi R^2 \int_{\theta_1}^{\theta_2} |B_R \sin\theta| \, d\theta, \quad (20) $$

respectively, where $v_R$ and $B_R$ are the radial velocity and magnetic field strength, respectively, and $B_{\phi}$ is the toroidal magnetic field strength. These global quantities can be split into additive contributions from the stellar wind, MEs, and accretion, which we hereafter distinguish with subscripts “wind,” “ME,\*,” and “acc,”<fnref n="5" /> respectively. The two angles $\theta_1$ and $\theta_2$ represent the area enclosing regions that correspond to each component.

For our simulations, we wish to investigate the effects of these global quantities at the stellar surface, which can be determined, in theory, by evaluating Equations (18)–(20) at $R = R_\star$. To avoid numerical effects near the inner boundary, we instead take the median value of mass flow rates and torques between the stellar surface and 2/3 of the truncation radius. However, we are still able to calculate the unsigned stellar magnetic flux by evaluating Equation (20) at the stellar surface, as the radial magnetic field strength component is kept fixed at its initial dipolar configuration. The total unsigned stellar magnetic flux can also be expressed analytically:

$$ \Phi_{\star} = \Phi_{\text{wind}} + \Phi_{\text{acc}} + \Phi_{\text{ME},\star} \equiv \alpha \pi R_{\star}^2 B_{\star}, \quad (21) $$

where $\alpha = 2$ for our dipolar surface configuration.

In Figure 3, using model 17 as a representative example, we plot the temporal evolution of the mass flow rates, torques, and unsigned magnetic fluxes (normalized by $\Phi_{\star}$) for the stellar wind, MEs, accretion, and the sum of these contributions. At the beginning of the simulation, there is a large transient event as the simulation adjusts from the initial state to a quasi-steady state, beginning around $t \approx P_\star$ for the case shown. In this relaxed state, quasi-periodic structure can be seen in the mass flows and torques, which are signatures of the MEs (green). Furthermore, it is clear that stellar wind (blue) global properties are slowly increasing with time, further demonstrating the quasi-steady configuration of these simulations. The unsigned accretion magnetic flux (orange) also decreases with time, conserving the unsigned stellar magnetic flux. This indicates that the truncation radius is slowly decreasing as the system evolves. Regardless, we take a time-averaged value for these global quantities (indicated by the gray shaded region) when simulations appear to become quasi-steady-state, which we take to be the latter half of the total timescale for each simulation. This method is applied for all simulations, and all values are listed in Table 1. In spite of the slow drift of the simulated quantities during the period of our time averages, we show in Section 4 that the global properties obey strict relationships and allow us still to derive robust relationships between global system properties.

Our time average mitigates any transient features, such as the periodic behavior of the MEs, and provides global quantities characteristic of the long-term evolution of the system. In the time-averaged domain, a significant proportion of the mass flow rate comes from material accreting onto the stellar surface. The mass-loss rate from the stellar wind and the mass inflow rate from the MEs are comparable. A majority of the net stellar torque is also associated with the spin-up contribution from the accretion flow, with a further contribution from the MEs (~30%), and a small spin-down contribution due to the stellar winds (~10%); thus, this system is in a net spin-up configuration. In these cases, the MEs accrete some of their material onto the star. This hierarchy is also present in the normalized unsigned magnetic flux; however, the high value

<footnote n="5">⁵ Global quantities with the “acc” subscript also encapsulate the “dead zone” region.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Evolution of the stellar (from top to bottom) mass flow rates, torques, and normalized unsigned magnetic fluxes, for the stellar wind (blue), MEs (green), accretion (orange), and the net value (black), as a function of time in units of stellar period, $P_{\star}$. Respectively colored dotted lines indicate time-averaged values that are taken over a period indicated by the gray shaded region. This figure refers to model 17.</figcaption></figure>

for accretion is indicative of the inclusion of the dead zone region. Therefore, a majority of the unsigned stellar magnetic flux is closed, and only a small fraction is open in the stellar wind, participating in the intermittent opening and closing of the MEs, or carrying the accretion flow.

To illustrate the extent of the explored parameter space, Figure 4 shows $\dot{J}_{\star}$ as a function of $f$ for all simulations, where the symbol shapes, colors, and borders encode the variations of each parameter. All of our simulations exist in the net spin-up torque regime. In general, a decrease in any one of these three input parameters decreases the spin-up torque experienced by the star. In Section 4, we show how each component of the net stellar torque depends on these input parameters.

# 4. Torque Formulations for Simulations of Star–Disk Interaction

## 4.1. Disk Truncation Radius

As disk material spirals inward and gets close enough to the star, the stresses arising from the stellar magnetosphere disrupt fluid motions, resulting in a truncation of the disk at some

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Net stellar torque, $J_{\star}$, as a function of the stellar break-up fraction, $f$, for all simulations. To differentiate between different input parameters, we use differing symbol styles for $f$, marker colors for $\rho_{d,\star}$, and marker edge colors for $B_{\star}$. Gray shaded region represents the net spin-down torque regime, while the transparent region represents the net spin-up torque regime; all of our simulations exist in the latter.</figcaption></figure>

height above the stellar surface, $R_t$. In our study, we determine the time-averaged truncation radius over the latter half of each simulation, where simulations appear to be in a quasi-steady state; these values (normalized by the stellar radius) can be found in Table 1. The truncation radius for a dipolar magnetic field structure is typically parameterized as (see, e.g., Romanova et al. 2002; Matt & Pudritz 2005a; Bessolaz et al. 2008; Zanni & Ferreira 2009)

$$ \frac{R_t}{R_*} = K_t \Upsilon_{\text{acc}}^{m_t}, \qquad (22) $$

where

$$ \Upsilon_{\text{acc}} = \frac{B_{\star}^2 R_{\star}^2}{4\pi |\dot{M}_{\text{acc}}| v_{\text{esc}}} \qquad (23) $$

is a “disk magnetization” parameter, which effectively represents the ratio of the accretion flow’s magnetic field and kinetic energies,<fnref n="6" /> $|\dot{M}_{\text{acc}}|$ is the absolute value of the mass accretion rate, and $v_{\text{esc}} = (2GM_{\star}/R_{\star})^{1/2}$ is the escape velocity from the stellar surface.<fnref n="7" /> Values for $\Upsilon_{\text{acc}}$ in our simulations can be found in Table 1.

In Figure 5, we plot the normalized truncation radius, $R_t/R_*$, as a function of $\Upsilon_{\text{acc}}$, which shows that all simulations fit the simple power-law formulation of Equation (22). Increasing $B_{\star}$ results in the disk being truncated further from the star, due to the increased magnetic stresses counteracting the pressure of the disk. The truncation radius is also inversely proportional to the mass accretion rate. There is no explicit dependence of $R_t$ on the stellar rotation rate or disk density, but we find that the mass accretion rate in our simulations does depend on these two parameters. We present our torque formulations here and throughout in terms of the global mass accretion rate onto the

<footnote n="6">⁶ Pantolmos et al. (2020) define $B_{\star}$ as the equatorial magnetic field strength, i.e., $B_{\star,\text{Pantolmos}} \equiv B_{\star}/2$, therefore $K_{t,\text{Pantolmos}} = 4m_t K_t$.</footnote>

<footnote n="7">⁷ A full derivation, including both the thermal and ram pressure, would introduce a factor of $[m_s/(1+m_s^2)]^{m_t}$ in Equation (22) (as shown in Appendix A), where $m_s$ is the sonic Mach number, which we neglect because we show that a constant $K_t$ fits the data well.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Normalized truncation radius, $R_t/R_*$, as a function of the disk magnetization parameter, $\Upsilon_{\text{acc}}$, for all simulations, where the symbols have the same meaning as Figure 4. Black line shows the best fit to Equation (22), giving $R_t/R_* = 0.756\Upsilon_{\text{acc}}^{0.340}$.</figcaption></figure>

star, so that the formulations are independent of some of the details of our accretion disk structure. In Appendix A, we show how the accretion rate can be predicted for our particular disk setup, and as a function of the fundamental input parameters.

The best-fit dimensionless parameters to Equation (22) are $K_t = 0.756$ and $m_t = 0.340$. This gives a steeper scaling than the value of $m_t = 2/7$ given by several previous analytic formulations (e.g., Elsner & Lamb 1977; Romanova et al. 2002; Matt & Pudritz 2005a; Bessolaz et al. 2008). Our simulations fit a steeper power law, (at least in part) because the magnetosphere is significantly perturbed relative to the initial, simple scaling of magnetic field strength with $R^{-3}$, which is discussed further in Section 4.3. For convenience, all of the best-fit dimensionless parameters presented in Section 4 can be found in Table 3.

## 4.2. Accretion Torque

The accretion torque accounts for the angular momentum carried by the disk material at $R_t$, which falls onto the star. It can be given by

$$ \dot{J}_{\text{acc}} = -|\dot{M}_{\text{acc}}|\Omega(R_t)R_t^2, \quad (24) $$

where $\Omega(R_t)$ is the rotation rate at the truncation radius. The analytical value of the accretion torque can be determined by assuming Keplerian rotation, i.e., $\Omega(R_t) \equiv \Omega_K(R_t) = (GM_*/R_t^3)^{1/2}$, giving $\dot{J}_{\text{acc}} = -|\dot{M}_{\text{acc}}|(GM_\star R_t)^{1/2}$. However, Zanni & Ferreira (2013) noted that the accretion disk is not in Keplerian rotation, due to the MEs extracting a fraction of the angular momentum from the surface of the disk. Therefore, we introduce the following functional fit for the accretion torque as a function of $R_t/R_*$:

$$ \dot{J}_{\text{acc}} = -K_{\text{acc}}|\dot{M}_{\text{acc}}|(GM_{\star}R_{\star})^{1/2}\left(\frac{R_t}{R_{\star}}\right)^{(1/2)+m_{\text{acc}}}, \quad (25) $$

where $K_{\text{acc}}$ and $m_{\text{acc}}$ are best-fit dimensionless parameters. In Figure 6, we plot the accretion torque as a function of the right-hand side of Equation (25), where the best-fit dimensionless

<table><caption>**Table 3**

Best-fit Dimensionless Parameters from Scaling Laws in Section 4</caption><thead><tr><th>Formulation</th><th>Parameter</th><th>Value</th><th>Equations</th></tr></thead><tbody><tr><td rowspan="2">R<sub>t</sub>/R<sub>*</sub></td><td>K<sub>t</sub></td><td>0.756</td><td rowspan="2">(22)</td></tr><tr><td>m<sub>t</sub></td><td>0.340</td></tr><tr><td rowspan="2">J<sub>acc</sub></td><td>K<sub>acc</sub></td><td>0.775</td><td rowspan="2">(25)</td></tr><tr><td>m<sub>acc</sub></td><td>-0.147</td></tr><tr><td rowspan="2">B<sub>θ</sub>(R<sub>t</sub>)</td><td>K<sub>B</sub></td><td>0.236</td><td rowspan="2">(31)</td></tr><tr><td>m<sub>B</sub></td><td>-2.54</td></tr><tr><td rowspan="2">J<sub>ME,*</sub></td><td>K<sub>ME,*</sub></td><td>0.00772</td><td rowspan="2">(33)</td></tr><tr><td>m<sub>ME,*</sub></td><td>-2.54</td></tr><tr><td rowspan="3">⟨r<sub>λ</sub>⟩/R<sub>*</sub></td><td>K<sub>λ,1</sub></td><td>1.13</td><td rowspan="3">(37)</td></tr><tr><td>K<sub>λ,2</sub></td><td>0.0356</td></tr><tr><td>m<sub>λ</sub></td><td>0.373</td></tr><tr><td rowspan="3">Φ<sub>wind</sub>/Φ<sub>*</sub></td><td>K<sub>Φ,1</sub></td><td>1.36</td><td rowspan="3">(42)</td></tr><tr><td>m<sub>Φ,1</sub></td><td>-1.34</td></tr><tr><td>m<sub>Φ,2</sub></td><td>0.0614</td></tr><tr><td>˙M<sub>acc</sub></td><td>K<sub>M</sub></td><td>0.00221</td><td>(A12)</td></tr></tbody></table>

parameters $K_{\text{acc}} = 0.775$ and $m_{\text{acc}} = -0.147$, giving a truncation radius power-law index of $(1/2) + m_{\text{acc}} = 0.353$ for our best-fit parameters. The accretion torque in our formulation predicts a lower value compared to the analytical formulation. For our parameter regime, one can predict the accretion torque via $\rho_{d,\star}$, using the mass accretion rate parameterization in Appendix A.

If the interpretation of Zanni & Ferreira (2013) is correct, the deviation of our formulation from the analytic one should be related to the disk rotation rate. Taken at face value, Equations (24) and (25) imply that the actual rotation of the disk at $R_t$ is

$$ \begin{aligned} \Omega(R_t) &= K_{\text{acc}} \left(\frac{R_t}{R_*}\right)^{m_{\text{acc}}} \Omega_K(R_t) \\ &= K_{\text{acc}} \left(\frac{GM_*}{R_*^3}\right)^{1/2} \left(\frac{R_t}{R_*}\right)^{-(3/2)+m_{\text{acc}}}. \end{aligned} \quad (26) $$

To test whether this is true, we computed $\Omega(R_t)$ in our simulations using the time-averaged local value, and we plot these as a function of truncation radius in Figure 7. The red line shows the expectation for purely Keplerian rotation ($K_{\text{acc}} = 1$ and $m_{\text{acc}} = 0$), which indicates that the simulated disks rotate more slowly than Keplerian, and their rotation has a slightly different radial dependence. The blue line shows the prediction of Equation (26), using the values of $K_{\text{acc}}$ and $m_{\text{acc}}$ that were fit to Equation (25). The approximate match between the blue line and data corroborates that the divergence of our accretion torque formulation from the analytical solution is primarily a result of the accretion disk rotating at sub-Keplerian velocity in the truncation region. The match of the prediction to the simulation values is not perfect, which may be due to the difficulty of defining a single value of $\Omega(R_t)$ in the simulations or may suggest that other mechanisms contribute a small additional adjustment to the accretion torque.

In a Keplerian disk, the corotation radius, $R_{\text{co}}$, is defined by $\Omega_* = \Omega_K(R_{\text{co}})$. For our simulations, we can now define the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The (absolute) accretion torque, $|\dot{J}_{acc}|$, as a function of the parameterization in Equation (25) for all simulations. Black line shows $y = x$, illustrating the goodness of fit. Symbols have the same meaning as in Figure 4, except that they are colored according to their truncation radius, in this and all following figures.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Rotation rate of disk material at the truncation radius, $\Omega(R_t)$, as a function of $R_t/R_*$ for all simulations. Red line shows the Keplerian solution, indicating that all simulation cases exhibit sub-Keplerian rotation. Blue line shows the prediction of Equation (26), but using the fit coefficients for the accretion torque of Equation (25). The approximate match of the blue line to the data indicates that sub-Keplerian rotation is the primary reason that accretion torques in the SDI are less than the analytic predictions.</figcaption></figure>

magnetic corotation radius as

$$R_{\text{co,m}} = \left[ K_{\text{acc}} \left( \frac{R_t}{R_*} \right)^{m_{\text{acc}}} \right]^{2/3} R_{\text{co}}. \quad (27)$$

It is clear that the location where the disk angular rotation rate equals that of the star is always smaller than the value for a Keplerian disk by a factor of $K_{\text{acc}}(R_t/R_*)^{m_{\text{acc}}}$. This value of $R_{\text{co,m}}$ is shown in Figures 1 and 2.

### 4.3. Stellar Magnetospheric Ejection Torque

The MEs exert a torque on the disk and star, as well as carrying some angular momentum out of the system in an outflow. This is a complex and highly dynamic region. Our

approach for finding a formulation for the torque on the star, due to MEs, is thus guided by the physics, but there are necessarily some heuristic aspects in the following formulations. We start with an analysis of the ME torque on the surface of the disk. The differential magnetic torque, due to the exertion of magnetic field lines threading an annulus of radial width $dR$ in a disk, is given as (see, e.g., Armitage & Clarke 1996; Matt & Pudritz 2005a)

$$\frac{d\dot{J}_{\text{mag}}}{dR} = qB_{\theta}^{2}R^{2}, \quad (28)$$

where $q$ is a measure of the twist of the field, expressed as the ratio of the toroidal to the vertical field at the disk surface, $B_{\theta}^+$ and $B_{\theta}$, respectively. Next, we will assume that some fraction of the torque exerted by the MEs on the disk is also exerted on the star. Recognizing that the size of the ME region ($\Delta R$) is small and bounded on one side by $R_t$, we can write the torque exerted on the star due to MEs as

$$\dot{J}_{\text{ME},*} \propto q(R_t) B_{\theta}(R_t)^2 R_t^2 \Delta R. \quad (29)$$

Thus, we need to include formulations for the magnetic twist ($q$), vertical field strength ($B_{\theta}$), and size of the ME region connected in the disk ($\Delta R$).

The field twist at the disk surface can be expressed in terms of the stellar and local rotation rates, to account for the differential rotation between the MEs and the star (see, e.g., Livio & Pringle 1992; Armitage & Clarke 1996; Matt & Pudritz 2005a; Gallet et al. 2019):

$$q(R_t) = \left. -\frac{B_{\theta}^{+}}{B_{\theta}} \right|_{R_t} \propto \left[ \frac{\Omega_{*}}{\Omega(R_t)} - 1 \right] \propto \left[ \left( \frac{R_t}{R_{\text{co,m}}} \right)^{3/2} - 1 \right]. \quad (30)$$

The polar magnetic field strength evaluated at the truncation radius can be written as

$$B_{\theta}(R_t) = K_B B_* \left( \frac{R_t}{R_*} \right)^{m_B}, \quad (31)$$

where $K_B$ and $m_B$ are best-fit dimensionless parameters. Figure 8 shows $B_{\theta}(R_t)/B_*$ (computed using time averages in our simulations) as a function of $R_t/R_*$. Red line shows the expectation for a simple dipole ($K_B = 0.5$, $m_B = -3$, which also corresponds to our initial conditions), indicating that the SDI perturbs the field in the simulations to have a slightly different radial dependence. The black line shows the best fit, giving $K_B = 0.236$ and $m_B = -2.54$. Finally, we assume the size of the ME connected region scales simply with $R_t/R_*$:

$$\Delta R/R_* \propto (R_t/R_*)^{m_{\Delta R}}, \quad (32)$$

where the index $m_{\Delta R}$ can be inferred from our final formulation below.

Combining Equations (29)–(32), we express the stellar ME torque as

$$\dot{J}_{\text{ME},*} = K_{\text{ME},*} \left[ \left( \frac{R_t}{R_{\text{co,m}}} \right)^{3/2} - 1 \right] B_*^2 R_*^3 \left( \frac{R_t}{R_*} \right)^{m_{\text{ME},*}}, \quad (33)$$

where $K_{\text{ME},*}$ and $m_{\text{ME},*}$ are the best-fit dimensionless parameters. It is clear that $m_{\text{ME},*}$ is equivalent to $2 + 2m_B + m_{\Delta R}$, and the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Vertical magnetic field strength in the disk midplane, $B_\theta(R_t)/B_\star$, as a function of $R_t/R_\star$ for all simulations. Red line shows a simple dipolar magnetic field ($B_\theta(R_t)/B_\star = 0.5 (R_t/R_\star)^{-3}$), indicating the field perturbations in the SDI results in a shallower radial dependence. Black line shows the fit to Equation (31), giving $B_\theta(R_t)/B_\star = 0.236 (R_t/R_\star)^{-2.54}$.</figcaption></figure>

constant $K_{ME,\star}$ combines all of the proportionality constants above, as well as the fraction of the ME torque that is applied to the star.<fnref n="8" /> Figure 9 shows the stellar ME torque as a function of the right-hand side of Equation (33), where the best-fit dimensionless parameters are $K_{ME,\star} = 0.00772$ and $m_{ME,\star} = -2.54$. In all of our simulations, the MEs add angular momentum to the star. Combining the parameterization of the truncation radius (Equations (22)–(23)), we can equivalently express the stellar ME torque in an alternative form, in terms of $R_t/R_\star$, $f$, and $|\dot{M}_{acc}|$:

$$ \dot{J}_{ME,\star} = \frac{4\sqrt{2}\pi K_{ME,\star}}{K_t^{1/m_t}} |\dot{M}_{acc}| (GM_\star R_\star)^{1/2} \times \left[ \frac{f}{K_{acc}} \left(\frac{R_t}{R_\star}\right)^{(3/2)-m_{acc}} - 1 \right] \left(\frac{R_t}{R_\star}\right)^{m_{ME,\star}+(1/m_t)} , \quad (34) $$

giving a positive truncation radius power-law index of $m_{ME,\star} + (1/m_t) = 0.401$ for our best-fit parameters. For our parameter regime, one can predict the stellar ME torque via $\rho_{d,\star}$, using the mass accretion rate parameterization in Appendix A.

### 4.4. Stellar Wind Torque

The stellar wind torque is parameterized as

$$ \dot{J}_{wind} = \dot{M}_{wind} \Omega_{\star} \langle r_{A} \rangle^{2}, \quad (35) $$

where $\langle r_A \rangle$ is the torque-average Alfvén radius, which is actually defined by Equation (35) (following Washimi & Shibata 1993; Matt & Pudritz 2008b). Thus, in our simulations,

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The (absolute) stellar ME torque, $|\dot{J}_{ME,\star}|$, as a function of the parameterization in Equation (33) for all simulations. Black line shows $y = x$, illustrating the goodness of fit.</figcaption></figure>

we compute $\langle r_A \rangle$ using

$$ \frac{\langle r_A \rangle}{R_\star} = \left( \frac{\dot{J}_{wind}}{\dot{M}_{wind} \Omega_\star R_\star^2} \right)^{1/2} . \quad (36) $$

This quantity represents an “effective magnetic lever arm” that quantifies the efficiency of the wind spin-down torque (see, e.g., Weber & Davis 1967). Values for $\langle r_A \rangle/R_\star$ in our simulations can be found in Table 1.

A torque formulation for stellar wind simulations was developed by Réville et al. (2015); it relates the Alfvén radius to the unsigned stellar wind (or open) magnetic flux. Using Equations (9) and (20), the open flux, $\Phi_{wind}$, can be determined by integrating between opening angles enclosing the stellar wind region at the stellar surface, i.e., evaluating Equation (20) at $R = R_\star$ between $0 \le \theta \le \theta_{wind,a}$ and $\theta_{wind,b} \le \theta \le \pi$, where $\theta_{wind,a}$ and $\theta_{wind,b}$ are the stellar wind anchoring angles for the northern and southern hemisphere, respectively. Réville et al. (2015) show that $\langle r_A \rangle/R_\star$ is related to $\Phi_{wind}$ via a single power law, independent of magnetic geometry:

$$ \frac{\langle r_A \rangle}{R_\star} = K_{A,1} \left\{ \frac{\Upsilon_{wind}}{[1 + (f/K_{A,2})^2]^{1/2}} \right\}^{m_A}, \quad (37) $$

where $K_{A,1}$, $K_{A,2}$, and $m_A$ are best-fit dimensionless parameters, and

$$ \Upsilon_{wind} = \frac{\Phi_{wind}^2}{4\pi R_{\star}^2 \dot{M}_{wind} v_{esc}} \quad (38) $$

is the “magnetization parameter” of the stellar wind, based on the open magnetic flux. Values for $\Upsilon_{wind}$ in our simulations can be found in Table 1. Figure 10 shows $\langle r_A \rangle/R_\star$ as a function of $\Upsilon_{wind}/[1 + (f/K_{A,2})^2]^{1/2}$, where the best-fit dimensionless parameters are $K_{A,1} = 1.13$, $K_{A,2} = 0.0356$, and $m_w = 0.373$. Increasing $\Phi_{wind}$ results in a increased Alfvén radius, giving a more efficient torque on the star (i.e., more torque for a given $\dot{M}_{wind} \Omega_\star$). The Alfvén radius is also inversely proportional to the stellar wind mass-loss rate, which itself is dependent on the coronal temperature or density, and $\Phi_{wind}$.

<footnote n="8">⁸ Note that the specific formulation/interpretation here is somewhat heuristic. An alternative interpretation, for example, could assume a constant $\Delta R$ and instead that there is some $R_t$-dependence in the fraction of the ME torque applied to the star.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 10. Normalized Alfvén radius in the stellar wind, $\langle r_A \rangle / R_*$, as a function of the wind open flux magnetization parameter, $\Upsilon_{wind}/[1 + (f/K_{A,2})^2]^{1/2}$, for all simulations. Black line shows the fit to Equation (37), giving $\langle r_A \rangle / R_* = 1.13 \{\Upsilon_{wind}/[1 + (f/0.0356)^2]^{1/2}\}^{0.373}$.</figcaption></figure>

Equation (38) can also be written as

$$ \Upsilon_{\text{wind}} = \left( \frac{\Phi_{\text{wind}}}{\Phi_{\star}} \right)^2 \Upsilon_{\star}, \qquad (39) $$

where

$$ \Upsilon_{\star} = \frac{\Phi_{\star}^2}{4\pi R_{\star}^2 \dot{M}_{\text{wind}} v_{\text{esc}}}. \qquad (40) $$

Therefore, in order to predict the stellar wind torque as a function of the properties at the stellar surface, given by $\Upsilon_{\star}$, we must have a formulation for the open magnetic flux, $\Phi_{wind}$.

In the SDI setup, the differential rotation between the star and disk opens magnetic surfaces at latitudes higher than where the field connects to $R_t$ (e.g., Uzdensky et al. 2002). The flux opened in this way includes the stellar wind flux (region (1), Figure 1), $\Phi_{wind}$, and also the flux participating in MEs (region (2)), $\Phi_{ME,\star}$, the latter of which opens and closes quasi-periodically. In our simulations, we find the sum of these fluxes to scale as

$$ \frac{\Phi_{\text{wind}} + \Phi_{\text{ME},\star}}{\Phi_{\star}} = K_{\Phi} \left( \frac{R_t}{R_{\star}} \right)^{m_{\Phi}}, \qquad (41) $$

where the best-fit dimensionless parameters are $K_{\Phi} = 1.26$ and $m_{\Phi} = -1.27$.<fnref n="9" />

Due to the complexity of MEs, it is not obvious how to predict and remove their contribution in Equation (41), although we expect that $\Phi_{ME,\star}$ depends on the magnetic twist and thus the spin rate of the star. In practice, we tried several ideas and found the best agreement with the simple formulation

$$ \frac{\Phi_{\text{wind}}}{\Phi_{\star}} = K_{\Phi,1} \left( \frac{R_t}{R_{\star}} \right)^{m_{\Phi,1}} f^{m_{\Phi,2}}. \qquad (42) $$

Figure 11 shows $\Phi_{wind}/\Phi_{\star}$ as a function of the right-hand side of Equation (42), showing our simulations to follow this two-variable, power-law scaling with the best-fit parameters

<figure><img src="image_12.png" /><figcaption>Figure 11. Fraction of open flux in stellar wind, $\Phi_{wind}/\Phi_{\star}$, as a function of the parameterization in Equation (42) for all simulations. Black line shows $y = x$, indicating goodness of fit.</figcaption></figure>

$K_{\Phi,1} = 1.36$, $m_{\Phi,1} = -1.34$, and $m_{\Phi,2} = 0.0614$. Qualitatively, this suggests that, for a fixed $R_t$, the fraction of flux in the MEs increases for slower stellar rotation (smaller $f$), which also corresponds to more star-disk differential rotation in the truncation vicinity. Therefore, for our SDI simulations, we find that $\Phi_{wind}/\Phi_{\star}$ strongly depends on the position of the truncation radius, and is independent of the stellar wind mass-loss rate, demonstrating that the accretion disk dominates over the wind dynamics. The fact that the wind’s open magnetic flux is independent of $\dot{M}_{wind}$ (and in fact depends on $\dot{M}_{acc}$, via $R_t$) represents a significant difference from the dynamics of ISWs (i.e., in systems without disks), and the implications of this are further discussed in Section 5.4.

Combining Equations (35), (37), (39), and (42), we express the stellar wind torque as

$$ \dot{J}_{\text{wind}} = K_{\text{A},1}^2 K_{\Phi,1}^{4m_{\text{A}}} \dot{M}_{\text{wind}} (GM_{\star})^{1/2} f^{1+4m_{\text{A}}m_{\Phi,2}} \left(\frac{R_t}{R_{\star}}\right)^{4m_{\text{A}}m_{\Phi,1}} \\ \times \left\{ \frac{\Upsilon_{\star}}{[1 + (f/K_{\text{A},2})^2]^{1/2}} \right\}^{2m_{\text{A}}}. \qquad (43) $$

Using Equations (20), (22), and (23), we can equivalently express the stellar wind torque (Equation (43)) in terms of $R_t/R_*$ and the mass ejection efficiency $\dot{M}_{wind}/|\dot{M}_{acc}|$:

$$ \dot{J}_{\text{wind}} = \frac{K_{\text{A},1}^2 (\alpha \pi K_{\Phi,1})^{4m_{\text{A}}}}{K_t^{2m_{\text{A}}/m_t}} |\dot{M}_{\text{acc}}|(GM_{\star}R_{\star})^{1/2} \left(\frac{\dot{M}_{\text{wind}}}{|\dot{M}_{\text{acc}}|}\right)^{1-2m_{\text{A}}} \\ \times f^{1+4m_{\text{A}}m_{\Phi,2}} \left[1 + \left(\frac{f}{K_{\text{A},2}}\right)^2\right]^{-m_{\text{A}}} \left(\frac{R_t}{R_{\star}}\right)^{4m_{\text{A}}m_{\Phi,1}+(2m_{\text{A}}/m_t)} , \qquad (44) $$

where $\alpha = 2$ for our dipolar configuration, giving a positive truncation radius power-law index of $4m_{A}m_{\Phi,1} + (2m_{A}/m_t) = 0.194$ for our best-fit parameters. Figure 12 shows the stellar wind torque as a function of the right-hand side of

<footnote n="9">⁹ If the field were a simple dipole, the analytic solution would be $K_{\Phi} = 1$ and $m_{\Phi} = -1$, so our best-fit values reflect the perturbation of the magnetospheric geometry by the SDI.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Stellar wind torque, $\dot{J}_{\text{wind}}$, as a function of the parameterization in Equation (44) for all simulations. Black line shows $y = x$, indicating goodness of fit.</figcaption></figure>

Equation (44), showing excellent correspondence with our stellar wind torque formulation.

# 5. Discussion and Conclusions

In this work, we perform 38 2.5D axisymmetric MHD SDI simulations to investigate how the net torque experienced by the star is impacted by variations in the stellar rotation rate (expressed as the stellar break-up fraction $f$), the stellar surface polar magnetic field strength ($B_*$), and the mass accretion rate (by changing the initial accretion disk density ($\rho_{d,\star}$)). All of our simulations are in a net spin-up regime. In this parameter study, we are motivated to further investigate the transport of angular momentum to and from the star via MEs during its PMS phase, and to produce improved torque formulations for the three mechanisms influencing the net stellar torque, namely, the stellar wind, the accretion, and the MEs. We create torque formulations for each, allowing for the investigation of the spin evolution of PMS objects using 1D stellar evolution codes. In Appendix B, we write a guide for implementing our torque formulations in a stellar evolution code, and discuss the range of validity when using our simulations.

## 5.1. Truncation Radius

We express a parameterization of the truncation radius with the “disk magnetization” parameter, $\Upsilon_{\text{acc}}$, where the proportionality constant $K_t = 0.756$ and the power-law index $m_t = 0.340$, which scales with $B_*$ and inversely with the mass accretion rate. This is a steeper scaling compared to the analytical power-law index ($m_t = 2/7$), due to the fact that the accretion disk perturbs the magnetosphere. Using an identical setup, Pantolmos et al. (2020) find $K_t = 0.704$ (corrected for polar magnetic field strength) and $m_t = 0.350$ for their $B_*$ parameter study, showing close agreement with our determined best-fit dimensionless parameters. Gallet et al. (2019) assume a dipolar magnetosphere and use $K_t = 0.766$ and $m_t = 2/7$ for their $R_t$ parameterization (corrected for polar magnetic field strength), thus, we predict higher truncation radii than those predicted by Gallet et al. (2019) for our parameter range. Kulkarni & Romanova (2013) found a much smaller $m_t = 0.2$ in their 3D simulations of magnetospheric accretion, for models

that have a dipolar magnetic field tilted 30° about the rotational axis. Therefore, further investigation into magnetic field inclination and its influence on $m_t$ would be required to make a direct comparison.

In Appendix A, we parameterize the mass accretion rate in terms of our three input parameters, allowing for the prediction of the truncation radius, accretion torque, and stellar ME torque, as a function of $\rho_{d,\star}$ rather than $\dot{M}_{\text{acc}}$.

## 5.2. Accretion Torque

For the accretion torque, we propose a formulation that considers the sub-Keplerianity of the accretion disk, due to the extraction of a fraction of angular momentum from the disk caused by MEs. We use a sub-Keplerian correction term $K_{\text{acc}}(R_t/R_*)^{m_{\text{acc}}} (<1)$ to represent the accretion torque in terms of the fraction of the analytical accretion torque (which assumes the accretion disk is rotating at Keplerian velocity), i.e., $\dot{J}_{\text{acc}} = K_{\text{acc}}(R_t/R_*)^{m_{\text{acc}}} \dot{J}_{\text{acc,K}}$. From this, we define a magnetic corotation radius $R_{\text{co,m}} < R_{\text{co}}$, where the magnetic twist becomes zero at $R_t = R_{\text{co,m}}$. For our simulations, $0.590 \le K_{\text{acc}}(R_t/R_*)^{m_{\text{acc}}} \le 0.695$, which is higher than the sub-Keplerian coefficient of $K_{\text{acc}} \approx 0.4$, with no additional dependence on $R_t$, that was used in the spin-evolution calculations of Gallet et al. (2019); we comment on the choice of coefficients in Gallet et al. (2019) in Section 5.5. Therefore, we predict a higher spin-up accretion torque for a given SDI case in our parameter regime. Pantolmos et al. (2020) adopted a fixed $K_{\text{acc}} = 0.790$ with no additional dependence on $R_t$, predicting an even higher accretion torque than our formulation; however, this value was determined using a smaller set of simulations, solely changing $B_*$.

## 5.3. Stellar ME Torque

Our stellar ME torque formulation offers additional self-consistency to the prescription introduced in Gallet et al. (2019) and adopted by Pantolmos et al. (2020). First, Gallet et al. (2019) use the term $K_{\text{rot}}$ to account for the sub-Keplerian correction to the magnetic field twist (or the star-disk differential rotation). However, our $K_{\text{acc}}(R_t/R_*)^{m_{\text{acc}}}$ parameter can be used interchangeably to represent both the accretion torque and field twist corrections. Their value for $K_{\text{rot}} = 0.7$ is at the upper limit of our $K_{\text{acc}}(R_t/R_*)^{m_{\text{acc}}}$ range, corresponding to our low $R_t/R_* \sim 2$ cases. Gallet et al. (2019) also assume a dipolar magnetic field when expressing the polar magnetic field strength at the truncation radius, whereas we employ the determined radial dependence (at the truncation radius) of the perturbed magnetic field structure. Finally, Gallet et al. (2019) assume a linear dependence between the radial width $\Delta R$ and the truncation radius, whereas we allow the dependence on the truncation radius to be fitted. Their formulation also assumes a relatively higher $K_{\text{ME},\star} = 0.0525$ (corrected for dipolar magnetic field strength) compared to our best-fit $K_{\text{ME},\star} = 0.0072$. Overall, we predict a lower spin-up stellar ME torque for a given SDI case in our parameter regime.

All simulations demonstrate a spin-up stellar ME torque, with our truncation radii well within their respective magnetic corotation radii. We observe the largest scatter between simulation and derived values of the stellar ME torque compared to other torque formulations, which is likely due to the periodicity of the stellar ME torque affecting our
time-averaged value, or from added complexity not considered
in our stellar ME torque parameterization.

5.4. Stellar Wind Torque

We determine a stellar wind torque formulation through the
parameterization of the effective Alfvén radius, based on the
open flux formulation employed by Réville et al. (2015). For a
known amount of open flux, our simulations follow a
$\langle r_A \rangle - \Upsilon_{\text{wind}} f$ relationship (Equation (37)) with a power-law
index of $m_A = 0.373$. ISW simulations with an initial dipolar
configuration have been shown to produce power-law indices
of $0.310 \le m_A \le 0.371$ (see, e.g., Réville et al. 2015; Finley &
Matt 2017; Pantolmos & Matt 2017). The power-law index
depends on the wind acceleration, which can vary from the
wind channel being modified by the presence of MEs, but also
varies between numerical studies, e.g., due to differences in
coronal temperature. A direct comparison of ISW and SDI
simulations in Pantolmos et al. (2020), exploring the $B_*$
parameter space, produces power-law indices of $m_A = 0.355$
and $m_A = 0.439$, respectively. Qualitatively, this is a result of
the accretion disk modifying the stellar wind geometry and
acceleration; a more detailed comparison between SDI and
ISW systems can be found in Pantolmos et al. (2020). The
deviations in their $m_A$ value from ours are likely within the
systematic uncertainties of these simulations. Generally, we
expect our values to be more precise because they are based on
a higher number of simulations and varied parameters. In
addition, the rotational dependence in Equation (37) introduces
added complexity, which may slightly modify the power-law
index, compared to the formulation of Pantolmos et al. (2020).

In order to compute the stellar wind torques in practice, one
must be able to predict the open flux based on the magnetic
properties at the stellar surface (e.g., B<sub>*</sub>) and other system
properties. In ISW simulations, the amount of open flux
depends on the wind mass-loss rate; specifically, Φ<sub>wind</sub>/Φ<sub>*</sub>
inversely scales with Υ<sub>*</sub> (see, e.g., Réville et al. 2015; Finley &
Matt 2017, 2018; Pantolmos & Matt 2017). By contrast, our
simulations show that, in the presence of an accretion disk (and
in the regime of M&#x0307;<sub>wind</sub>/|M&#x0307;<sub>acc</sub>| &le; 1%), the open flux depends
primarily on R<sub>t</sub> (and thus M&#x0307;<sub>acc</sub>) and results in an increase of the
open flux, compared to an ISW case. The apparent scaling of
the open flux with Υ<sub>*</sub> in the SDI simulations of Pantolmos et al.
(2020) is a consequence of a relationship between the mass
accretion rate and wind mass-loss rate, due to the fact that they
only vary B<sub>*</sub>; solely increasing B<sub>*</sub> causes both mass fluxes to
increase. The independence of the open flux on M&#x0307;<sub>wind</sub> (and its
dependence instead on M&#x0307;<sub>acc</sub>) represents a significant difference
from the dynamics of ISWs. Once the amount of open flux is
determined, the stellar wind flows along all flux that is available,
extracting angular momentum. The somewhat surprising result
of this is that the stellar wind torque in SDI systems depends
upon the accretion rate onto the star, in addition to the usual
stellar wind parameters. As shown in Section 4.4, for SDI
cases, J&#x0307;<sub>wind</sub> &prop; f<sup>1.09</sup> M&#x0307;<sub>wind</sub><sup>0.254</sup> |M&#x0307;<sub>acc</sub>|<sup>0.680</sup> B<sub>*</sub><sup>0.132</sup> (using best-fit para-
meters from Table 3 and neglecting the wind centrifugal term<sup>10</sup>),
whereas for ISW cases, J&#x0307;<sub>wind</sub> ISW &prop; fM&#x0307;<sub>wind</sub><sup>0.56</sup> B<sub>*</sub><sup>0.88</sup> (e.g., Matt
et al. 2012). For fixed stellar wind parameters, the stellar wind
torque in the SDI increases when M&#x0307;<sub>acc</sub> increases, due to a
decrease in R<sub>t</sub> and thus increase in open magnetic flux. The

dependence of the wind torque on magnetic field strength is
also much weaker than for ISWs because, in SDI systems, the
fractional open flux decreases almost linearly with an increase
in B<sub>*</sub> (due to the sensitivity of fractional open flux on R<sub>t</sub> and R<sub>*</sub>
on B<sub>*</sub>), such that the total amount of open flux is nearly
independent of B<sub>*</sub>.

Since our simulations exhibit more open magnetic flux than
ISW cases with the same wind parameters, their stellar wind
torques are more efficient (i.e., they have more torque for a
given B<sub>*</sub> and M&#x0307;<sub>wind</sub>) than ISWs. This could not be the case for
high mass-loss rates and/or large truncation radii. For example,
comparing our SDI formulation to the ISW formulation in
Finley & Matt (2017) with f = 0.05, for R<sub>t</sub>/R<sub>*</sub> &ge; 3 and
M&#x0307;<sub>wind</sub>/|M&#x0307;<sub>acc</sub>| &ge; 10%, the torque of a stellar wind in an SDI
system becomes less efficient than an ISW with the same mass-
loss rate and magnetic field. This is mainly due to the fact that,
in these conditions, according to the derived scalings, a more
massive stellar wind could open more magnetic flux than the
SDI. However, our formulations are based on simulations with
M&#x0307;<sub>wind</sub>/|M&#x0307;<sub>acc</sub>| &le; 1%, so it is not clear if our formulations can be
extrapolated to ~10%. It is possible that, when M&#x0307;<sub>wind</sub>/|M&#x0307;<sub>acc</sub>| is
large enough, the stellar wind could again contribute to the
opening of the magnetic flux in SDI systems as well, but when
this transition occurs remains unknown. Thus, until the regime
of high M&#x0307;<sub>wind</sub>/|M&#x0307;<sub>acc</sub>| is explored with simulations, caution
should be used when extrapolating our formulations. For
clarity, we note that we are solely interested in how stellar
wind torques scale in the presence of SDI; comparing to the
ISW case gives us an indication of potential limits to our SDI
scalings, when the wind becomes “important enough” relative
to the disk (such that the stellar wind torque starts to scale like
an ISW).

For the stellar wind torque, Gallet et al. (2019) use the Matt
et al. (2012) ISW prescription with a 30% higher value for K<sub>1</sub>.
However, their formulation is not sensitive to the accretion rate,
resulting in increasing underestimations of the stellar wind
torque as M&#x0307;<sub>acc</sub> increases. Therefore, we predict a higher spin-
down stellar wind torque for a given SDI case in our parameter
regime.

5.5. Spin-equilibrium Solution

We use our torque formulation to explore outside the
parameter regime studied and predict equilibrium conditions
for our simulations, defined as the conditions where the net
stellar torque is zero, i.e., *J*<sub>wind</sub> + *J*<sub>acc</sub> + *J*<sub>ME,*</sub> = 0. We give
the formulation for the equilibrium state in Appendix C, where
we show that the equilibrium stellar spin rate (as a fraction of
break-up), *f*<sub>eq</sub>, can be written in terms of the truncation radius
and the mass ejection efficiency alone. In Figure 13, we plot *f*<sub>eq</sub>
as a function of *R*<sub>t</sub>/*R*<sub>*</sub>. The symbols show our simulation
values, but plotted at the spin rate the star would need to have
in order for the net torque in each simulation to be zero. The
blue and orange solid lines show continuous solutions for
*M&#x0307;<sub>wind</sub>/|M&#x0307;<sub>acc</sub>| = 0.001 and 0.01 (using Equation (C13)),
encapsulating the majority of our simulations. The apparent
approximate inverse relationship with truncation radius illus-
trates that spin-equilibrium can be reached at lower spin rates
as the effective magnetization (relative to the mass accretion
rate) increases.

The lower black line shows $f_{\text{ME}} = K_{\text{acc}}(R_t/R_*)^{-(3/2)+m_{\text{acc}}}$,
which illustrates the transition from a spin-up (below line) and
spin-down (above line) stellar ME torque contribution for these

<footnote n="10"><sup>10</sup> The wind centrifugal term is [1 + (f/K<sub>A,2</sub>)<sup>2</sup>]<sup>-m<sub>A</sub></sup>, and generally is nearly
equal to unity.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Illustration of $f_{\text{eq}}$ as a function of $R_t/R_*$ for all simulations, and for $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}| = 0.001$ (blue line) and 0.01 (orange line), using Appendix C. We plot $f_{\text{ME}}$ (black solid line) to differentiate between spin-up and spin-down stellar ME torque regions, and $f_{\text{prop}}$ (black solid line) to indicate the transition into the propeller regime. Comparisons are made with spin-equilibrium solutions from Long et al. (2005) (cyan dashed line), and Gallet et al. (2019) for $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}| = 0.001$ (green dashed line) and 0.01 (red dashed line).</figcaption></figure>

spin-equilibrium cases. This $\dot{J}_{\text{ME},*} = 0$ condition corresponds to a situation in which the disk is still sub-Keplerian so that, in principle, it can still steadily accrete. We also plot $f_{\text{prop}} = (R_t/R_*)^{-3/2}$, which illustrates the transition into the propeller regime (above line), which corresponds to a situation in which the stellar centrifugal barrier hinders accretion by triggering super-Keplerian rotation.

Our formulation predicts that, in order for our simulated cases to have a net zero torque, a majority of our simulations must be in the propeller regime. In these cases, the stellar ME torque (Equation (33)) in fact provides a spin-down torque onto the star, predicting a slower spin rate than if MEs were neglected. In Figure 13, the $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}| = 0.01$ solution (orange line) sits below the propeller regime transition (black line) when $R_t/R_* < 4$. Therefore, our spin-equilibrium solutions suggest that it is possible to be in a spin-down stellar ME torque configuration without being in the propeller regime. Since our spin-equilibrium solutions either exist in the propeller regime or have a spin-down stellar ME torque configuration, it is not possible to actually verify whether our formulation is valid in such a regime at present. It is possible that the open flux could saturate for $R_t \gtrsim R_{\text{co}}$, since it is unlikely to have steady accretion beyond this limit; even if accretion becomes intermittent, it is necessary that $R_t < R_{\text{co}}$, even just temporarily.

The $J_*=0$ solution using the torque formulation in Gallet et al. (2019) (labeled G19) predicts slower values of $f_{\text{eq}}$ compared to our simulations (except for very small $R_t$). Overall, our formulation predicts a lower total spin-up torque (accretion plus MEs) and a higher total spin-down torque (stellar wind), giving a lower net spin-up stellar torque compared to the formulation in Gallet et al. (2019). Therefore, for $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}| = 0.01$, the lower $f_{\text{eq}}$ values provided by Gallet et al. (2019) are mostly due to the lower accretion torque. For $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}| = 0.001$, the more efficient spin-down torque from the MEs in the Gallet et al. (2019) formulation contributes significantly to decrease the $f_{\text{eq}}$ value. The solution in Long et al. (2005) (labeled L05), where $R_{\text{co}} = 1.4R_t$, predicts slower

values of $f_{\text{eq}}$ compared to our simulations. This also demonstrates that the Gallet et al. (2019) solution is consistent with the $R_t \propto R_{\text{co}}$ scaling proposed by this work.

As an alternative to computing $f_{\text{eq}}$, we can also fix $f$ and $R_t$ at the value of each simulation and then solve for which value of mass ejection efficiency ($\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}|$) gives a net zero torque. Doing this, we find our formulation to suggest that $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}| \gtrsim 30\%$ is required to achieve spin-equilibrium in our simulated cases, with the lowest mass fraction corresponding to our most rapidly rotating cases ($f=0.0625$). Using $R_*, M_*, B_*, P_*$, and $\dot{M}_{\text{acc}}$ observations of BP Tau as a representative CTTS case (see Johnstone et al. 2014; Gallet et al. 2019), our formulation suggests $\dot{M}_{\text{wind}}/|\dot{M}_{\text{acc}}|$ would need to be $\approx 25\%$ for this star to be in spin-equilibrium. We are skeptical of such high mass-loss rates in CTTS; assuming that these high mass-loss rates are not possible, the implication of our current torque formulation is that BP Tau could in fact be spinning up. However, this result is a prediction obtained solely via extrapolation of our current torque formulation, which may not be compatible within other parameter regimes, e.g., for stars in the propeller regime. The validity of our formulation outside of the low mass-loss regime must be investigated in future work.

## 5.6. Limitations

Due to the complexity of SDI systems, various limitations of our study are listed hereafter. First, we treat the accretion disk, and its ability to transfer angular momentum, using simple viscosity and diffusion formulations, with constant anomalous coefficients for each. These coefficients control the viscous torque and thus the mass accretion rate. We have not varied these coefficients or related physics, so it is not yet clear whether or how these would affect the truncation radius or the open flux for a given accretion rate (which is how our formulations are phrased).

All of our simulations exist in the net spin-up regime. Furthermore, MEs contribute solely spin-up torque to the star for all simulations. Therefore, we were unable to explore the possibility of MEs helping the star reach an equilibrium or net spin-down situation. Our torque formulations predict that increasing the magnetic field strength or decreasing the mass accretion rate (extrapolating beyond the parameter regime simulated) could push simulations into the propeller regime, where MEs could produce a spin-down contribution to the star and provide solutions that approach or even reach a net spin-down configuration. Alternatively (or in addition), more massive stellar winds, i.e., higher stellar wind mass-loss rates, may enable the transition into this net spin-down regime, and could be achieved by increasing the coronal temperature. There are also limitations to using a simple polytropic Parker wind to model our stellar winds, so their properties should be taken to be hypothetical. We are still unsure what mass-loss rates are realistic for these stars, as well as what drives these stellar winds. The outputted mass-loss rates should not be taken as a true prediction of these simulations, but they are a result of unconstrained parameters, such as the coronal density, temperature, and adiabatic index. Finally, there may be other mechanisms that are not captured by our simulations.

A proper modeling of the thermal structure of the disk should take into account a plethora of effects, such as internal heating, stellar illumination and disk geometry, opacity, emissivity, dust content, and dust-to-gas ratio (see, e.g., Dullemond & Monnier 2010), that we neglected and replaced
with a simple adiabatic assumption. It is widely accepted that nonideal MHD effects play a fundamental role in determining the accretion dynamics and the disk structure (see, e.g., Lesur et al. 2014; Béthune et al. 2017), but we used a simpler alpha parameterization to control the gas–magnetic field coupling.

The assumption of axisymmetry neglects 3D effects that may result in deviations of the open flux, the amount of flux participating in the MEs, and the size of the accretion flow, i.e., for a fixed accretion rate, there may be differences in the density within the funnel flow, as well as the truncation radius, and thus the flow dynamics. To improve upon this, future studies must consider 3D SDI simulations. It is possible that the formulation developed in this paper can be applied to even multipolar cases if the dipolar component indeed determines the position of the truncation radius, i.e., if the dipole is stronger than any other multipole in the truncation region and the multipoles do not introduce any closed loops at latitudes higher than the accretion spots, as the stellar flux available to accelerate a wind at latitudes higher than the accretion spots would be equivalent to that in our purely dipolar cases. However, to test this, future studies should consider more realistic magnetic field topologies, i.e., nonaxisymmetry and additional multipolar components, and investigate the effectiveness of our torque formulations for these increasingly complex cases.

L.G.I. and S.P.M. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 682393; AWESOME Stars: Accretion, Winds, and Evolution of Spins and Magnetism of Stars; http://empslocal.ex.ac.uk/AWESOME Stars). The authors would like to acknowledge the use of the University of Exeter High-Performance Computing (HPC) facility in carrying out this work. G.P. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 742095; SPIDI: Star–Planet–Inner Disk Interactions); https://www.spidi-eu.org. We thank Andrea Mignone and others for the development and maintenance of the PLUTO code. Figures within this work are produced using the Python package Matplotlib (Hunter 2007).

Software: Matplotlib (Hunter 2007), PLUTO (Mignone et al. 2007, 2012).

## Appendix A

### Parameterization of the Mass Accretion Rate

Our formulation calculates the stellar torques for a specified or known mass accretion rate, so the formulations should be somewhat independent of the details of the physical processes and properties that determine $\dot{M}_{acc}$ in the disk. In our simulations, $\dot{M}_{acc}$ is a result of the particular initial disk profile, input parameters, and the physics included. In this appendix, we show how to predict $\dot{M}_{acc}$ from our simulation input parameters, although this may only apply for the initial disk conditions and physics that we have included in our simulations. The parameterization of $\dot{M}_{acc}$ allows for the truncation radius, and hence the accretion and stellar ME torque formulations, to be determined a priori from $\rho_d$.

### A.1. Poloidal Equilibrium Derivation of the Truncation Radius

The truncation radius can be derived from the poloidal equilibrium of the accretion disk and the stellar magnetic field. With $\rho$, $p$, and $v_R$ being the disk midplane density, thermal pressure and accretion speed, respectively, $c_s = (p/\rho)^{1/2}$ being the equatorial isothermal sound speed, and $m_s = |v_R/c_s|$ being the sonic Mach number of the accretion flow, this is expressed as

$$ \frac{B_{\theta}^{2}}{8\pi} = \frac{p + \rho v_{R}^{2}}{\beta_{\text{tot}}} = \frac{p(1 + m_{s}^{2})}{\beta_{\text{tot}}}, \qquad (A1) $$

where $B_\theta$ is the poloidal field at the disk midplane and the parameter $\beta_{tot}$ is of order unity. The usual plasma $\beta = 8\pi p/B_\theta^2$ parameter is therefore

$$ \beta = \frac{\beta_{\text{tot}}}{1 + m_s^2}. \qquad (A2) $$

Notice that here we are considering a region of the disk that is still accreting along the midplane. Our definition of the truncation radius corresponds to the position where the accretion flow starts to deviate and uplift to form the accretion columns, not to the radius at which the accretion flow is disrupted. In other words, our definition of the truncation radius corresponds to the intersection with the magnetic surface that envelops the accretion columns, not the one onto which the accretion columns lean on.

Taking the thermal height scale of the disk $H = Rc_s/v_K = c_s/\Omega_K$ where $v_K = (GM_*/R)^{1/2}$ and $\Omega_K = v_K/R$ are the Keplerian toroidal and angular speeds, respectively, one can express the disk density and pressure as a function of the mass accretion rate $\dot{M}_{acc}$:

$$ \begin{aligned} \rho &= \frac{\dot{M}_{\text{acc}}}{4\pi H R v_R}, \\ p &= \rho c_s^2 = \frac{\dot{M}_{\text{acc}} v_K}{4\pi R^2 m_s^2}. \end{aligned} \qquad (A3) $$

Therefore, Equation (A1) can be rewritten as

$$ \frac{2\dot{M}_{\text{acc}} v_K}{R^2 B_{\theta}^2} = \frac{\beta_{\text{tot}} m_s}{1 + m_s^2}, \qquad (A4) $$

demonstrating that, while the truncation can be seen as a poloidal pressure equilibrium, it also depends on the Mach number and accretion speed, which are determined by the angular momentum transport that drives accretion.

### A.2. Angular Momentum Conservation

Using the $\phi$ component of the momentum equation (Equation (1)), we write the stationary angular momentum conservation equation, considering the torque due to a large-scale magnetic field and a viscous torque:

$$ \rho v_R \frac{d}{dR} (R^2 \Omega) = \frac{R B_{\theta} d B_{\phi}}{4\pi dz} + \frac{1}{R} \frac{d}{dR} (R^2 \tau_{R\phi}), \qquad (A5) $$

where $B_\theta$ and $B_\phi$ are the vertical and the toroidal components of the magnetic field strength, respectively, and $\tau_{R\phi} = \rho v_R R (d\Omega/dr)$ is the $R\phi$ component of the viscous stress tensor (see Equation (3)). By assuming an $\alpha$ parameterization (Shakura & Sunyaev 1973) for
the kinematic viscosity, $\nu_v = \alpha_v c_s^2 / \Omega_K = \alpha_v c_s^2 R / \nu_K$. In a Keplerian disk, $\tau_{R\phi} = -\alpha_v p$. We define the dimensionless quantity $Q_{rot} = v_K^{-1} d(R^2\Omega)/dR$ (where $Q_{rot} = 1/2$ for a Keplerian disk). Therefore, Equation (A5) becomes

$$Q_{\text{rot}} \rho v_R v_K = \frac{R B_0}{4\pi} \frac{d B_{\phi}}{dz} - \frac{\alpha_v}{R} \frac{d}{dR} (R^2 p). \quad (A6)$$

Assuming $p = c_s^2 \rho$, $dB_{\phi}/dz \approx B_{\phi}^+/H$, $q = -B_{\phi}^+/B_0$, ($\alpha_v/R$) $d(R^2p)/dR \approx \alpha_v$, and $\epsilon = H/R = c_s/v_K$, and $m_s = |v_R|/c_s$, we can rewrite Equation (A5) as

$$m_s = \left| \frac{2q}{Q_{\text{rot}}\beta} + \alpha_v \epsilon \right|, \quad (A7)$$

where $\beta = 8\pi p / B_0^2$. While viscous accretion is usually strongly subsonic (typically $\alpha_v \epsilon \ll 1$), an equipartition field, with both $q$ and $\beta$ of order unity, tends to yield a sonic accretion velocity. Since these are the typical conditions in the truncation region, it is safe to neglect the viscous torque:

$$m_s \approx \frac{2|q|}{Q_{\text{rot}}\beta}. \quad (A8)$$

### A.3. Mass Accretion Rate Scaling

If we assume the typical SDI parameterization for $q$:

$$q = K_q \left( \frac{\Omega_*}{\Omega} - 1 \right), \quad (A9)$$

where $K_q$ is a best-fit dimensionless parameter, and neglect the viscous torque, we can combine Equations (A2), (A4), and (A8) in order to obtain the following set of equations:

$$\frac{2\dot{M}_{\text{acc}}v_K}{R^2 B_0^2} = \frac{\beta_{\text{tot}}m_s}{1 + m_s^2}, \\ \frac{K_q B_0^2}{4\pi Q_{\text{rot}} p} \left| \frac{\Omega_*}{\Omega} - 1 \right| = m_s. \quad (A10)$$

Here, $\beta_{\text{tot}}$ is plausibly of order unity, since the magnetic field should be able to stop the disk accretion flow. Sonic accretion, together with a magnetic field that can balance both the ram and thermal disk pressure, seems to represent $m_s$ of order unity as an optimal condition to truncate the disk. Sonic accretion, along with an equipartition poloidal field, also corresponds to a field twist $q \lesssim 1$, which is compatible with the maximum twist allowed to avoid the inflation of the closed magnetosphere. We notice also that the factor $m_s/(1+m_s^2)$ has a maximum at $m_s$, meaning a $m_s = 1$ condition roughly corresponds to the outermost position where the conditions for truncation are met.

The first expression in Equation (A10) provides the familiar relation for the truncation radius $R_t = R_t(\dot{M}_{\text{acc}}, B_*)$ (see Equation (22)), while the second provides an implicit relation for the truncation radius $R_t = R_t(Q_{\text{rot}}(R_t), \Omega_*, \Omega(R_t))$. At the truncation radius, we assume that the thermal disk pressure scales with the initial disk pressure profile (Equation (5)), i.e., $p(R_t) \propto \rho_{d,*} v_{\text{esc}}^2 (R_t/R_*)^{-5/2}$ (for fixed disk aspect ratio $\epsilon$), $B_\theta(R_t) \approx B_\theta(R_t) = K_B B_\star (R_t/R_*)^{m_B}$ (Equation (31)), and set $\Omega_*/\Omega(R_t) = (R_t/R_{\text{co,m}})^{3/2} = (f/K_{\text{acc}})(R_t/R_*)^{(3/2)-m_{\text{acc}}}$

<figure><img src="image_15.png" /><figcaption>Figure 14. Outputed $|\dot{M}_{\text{acc}}|$ as a function of the parameterization in Equation (A12) for all simulations. Black line shows $y = x$, indicating goodness of fit, giving $K_M = 0.00221$.</figcaption></figure>

(Equation (27)), giving

$$\frac{R_t}{R_*} = \frac{m_s' Q_{\text{rot}}}{2K_q K_B^2} \left\{ \left| \frac{f}{K_{\text{acc}}} \left( \frac{R_t}{R_*} \right)^{(3/2)-m_{\text{acc}}} - 1 \right| \left( \frac{B_*^2}{4\pi \rho_{d,*} v_{\text{esc}}^2} \right) \right\}^{-2/(5+4m_B)}. \quad (A11)$$

The $m_s'$ factor does not correspond to the actual Mach number in the evolved solution, due to the assumptions made above. Using our $R_t$ parameterization (Equations (22)–(23)), a parameterization of the mass accretion rate for our simulations can be written as

$$|\dot{M}_{\text{acc}}| = K_M \frac{B_*^2 R_*^2}{4\pi v_{\text{esc}}} \left\{ \left| \frac{K_t^{(3/2)-m_{\text{acc}}} f}{K_{\text{acc}}} \left( \frac{B_*^2 R_*^2}{4\pi |\dot{M}_{\text{acc}}| v_{\text{esc}}} \right) \right|^{m_t[(3/2)-m_{\text{acc}}]} - 1 \left| \frac{-1}{4\pi \rho_{d,*} v_{\text{esc}}^2} \right|^{-2/[m_t(5+4m_B)]} \right\}, \quad (A12)$$

where $K_M = [m_s' Q_{\text{rot}}/(2K_q K_B^2 K_t)]^{-1/m_t}$ is a best-fit dimensionless parameter. This is a transcendental equation that must be solved iteratively. In Figure 14, we plot the outputed mass accretion rate as a function of the theoretical mass accretion rate calculated using Equation (A12), and show excellent correspondence between these with best-fit dimensionless parameter $K_M = 0.00221$.

## Appendix B

### Using Our Formulation to Calculate Net Torque

For clarity, here is a simple guide to using the torque formulations derived in this study.

1. Specify the following input parameters: stellar mass ($M_*$), stellar radius ($R_*$), surface polar magnetic field strength ($B_*$), stellar break-up fraction ($f$), mass accretion rate ($\dot{M}_{\text{acc}}$), and stellar wind mass-loss rate ($\dot{M}_{\text{wind}}$).
Alternatively, one can use Appendix A to specify the initial stellar surface disk density ($\rho_{d,\star}$) rather than the mass accretion rate.

2. Calculate the truncation radius ($R_t$), using Equations (22)–(23).

3. Calculate the accretion torque ($\dot{J}_{acc}$), stellar ME torque ($\dot{J}_{ME,\star}$), and stellar wind torque ($\dot{J}_{wind}$) using Equations (25), (34), and (44), respectively.

4. Sum each torque contribution to calculate the net stellar torque ($\dot{J}_* = \dot{J}_{acc} + \dot{J}_{ME,\star} + \dot{J}_{wind}$).

SDI simulations in this study, and in general, are performed over several dynamical timescales, but evolutionary timescales of these stars are on the order of millions of years. Therefore, to predict the spin evolution of these stars, one can determine the net stellar torque for a star with given global properties, using torque formulations derived from these SDI simulations, and then evolve this over time using a stellar evolution code. However, it is important to consider the validity of our torque formulation in regimes not considered by the parameter study; for example, if the truncation radius becomes smaller than the stellar radius, i.e., $R_t < R_\star$ (outside the magnetic accretion regime), we recommend enforcing a limit of $R_t \ge R_\star$. On the other hand, if $R_t$ exceeds the magnetic corotation radius ($R_t > R_{co,m}$, i.e., the net spin-up stellar torque regime), our physically motivated formulations can be extrapolated and used here, but doing so will be uncertain because real systems may transition to different regimes in parameter ranges far outside what we have explored. The situation is even more unclear for the propeller regime (where $R_t > R_{co}$), where the accretion flow is expected to be significantly disrupted.

## Appendix C
Spin-equilibrium Solution for Our Formulation

By setting the sum of the accretion (Equation (25)), stellar ME (Equation (34)), and stellar wind (Equation (44)) torques equal to zero, i.e., $\dot{J}_{wind} + \dot{J}_{acc} + \dot{J}_{ME,\star} = 0$, one can express the following transcendental equation (using best-fit dimensionless parameters in Table 3):

$$58.2 \left( \frac{\dot{M}_{\text{wind}}}{|\dot{M}_{\text{acc}}|} \right)^{0.254} f_{\text{eq}}^{1.09} \left[ 1 + \left( \frac{f_{\text{eq}}}{0.0356} \right)^2 \right]^{-0.373} \left( \frac{R_t}{R_\star} \right)^{0.191} \\ - 0.775 \left( \frac{R_t}{R_\star} \right)^{0.353} + 0.312 \left[ \frac{f_{\text{eq}}}{0.775} \left( \frac{R_t}{R_\star} \right)^{1.65} - 1 \right] \\ \times \left( \frac{R_t}{R_\star} \right)^{0.402} = 0, \qquad (\text{C1})$$

where $f_{eq}$ is the spin rate (as a fraction of break-up) for the star to have a net zero torque, i.e., the “equilibrium spin rate.” To calculate $f_{eq}$, for a given $R_t$ and $\dot{M}_{wind}/|\dot{M}_{acc}|$, this expression must be solved iteratively. Alternatively, if the spin rate is specified/fixed, one can instead solve for the predicted equilibrium mass ejection efficiency for a given $R_t$ (or vice versa).

## ORCID iDs

Lewis G. Ireland https://orcid.org/0000-0002-8833-1204

Claudio Zanni https://orcid.org/0000-0003-0204-8190

Sean P. Matt https://orcid.org/0000-0001-9590-2274

George Pantolmos https://orcid.org/0000-0001-7788-3727

## References

Agapitou, V., & Papaloizou, J. C. B. 2000, MNRAS, 317, 273

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20

Armitage, P. J., & Clarke, C. J. 1996, MNRAS, 280, 458

Bessolaz, N., Zanni, C., Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 155

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

Bouvier, J., Cabrit, S., Fernandez, M., Martin, E. L., & Matthews, J. M. 1993, A&A, 272, 176

Camenzind, M. 1990, RvMA, 3, 234

Cranmer, S. R. 2008, ApJ, 689, 316

Cranmer, S. R. 2009, ApJ, 706, 824

D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411

D'Angelo, C. R., & Spruit, H. C. 2010, MNRAS, 406, 1208

Dedner, A., Kemm, F., Kröner, D., et al. 2002, JCoPh, 175, 645

Donati, J. F., Bouvier, J., Alencar, S. H., et al. 2019, MNRAS, 483, L1

Donati, J. F., Bouvier, J., Alencar, S. H., et al. 2020, MNRAS, 491, 5660

Dullemond, C. P., & Monnier, J. D. 2010, ARA&A, 48, 205

Edwards, S., Fischer, W., Hillenbrand, L., & Kwan, J. 2006, ApJ, 646, 319

Edwards, S., Hartigan, P., Ghandour, L., & Andrulis, C. 1994, AJ, 108, 1056

Edwards, S., Strom, S. E., Hartigan, P., et al. 1993, AJ, 106, 372

Elsner, R. F., & Lamb, F. K. 1977, ApJ, 215, 897

Finley, A. J., & Matt, S. P. 2017, ApJ, 845, 46

Finley, A. J., & Matt, S. P. 2018, ApJ, 854, 78

Gallet, F., Zanni, C., & Amard, L. 2019, A&A, 632, A6

Ghosh, P., & Lamb, F. K. 1979, ApJ, 234, 296

Gregory, S. G., Matt, S. P., Donati, J. F., & Jardine, M. 2008, MNRAS, 389, 1839

Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Herbst, W., Eislöffel, J., Mundt, R., & Scholz, A. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 297

Hunter, J. D. 2007, CSE, 9, 90

Illarionov, A. F., & Sunyaev, R. A. 1975, A&A, 39, 185

Ingleby, L., Calvet, N., Hernández, J., et al. 2014, ApJ, 790, 47

Irwin, J., & Bouvier, J. 2009, in IAU Symp. 258, The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse (Cambridge: Cambridge Univ. Press), 363

Johns-Krull, C. M. 2007, ApJ, 664, 975

Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J. F., & Hussain, G. 2014, MNRAS, 437, 3202

Koenig, A. 1991, ApJL, 370, L39

Kulkarni, A. K., & Romanova, M. M. 2013, MNRAS, 433, 3048

Kurosawa, R., Romanova, M. M., & Harries, T. J. 2011, MNRAS, 416, 2623

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

Livio, M., & Pringle, J. E. 1992, MNRAS, 259, 23P

Long, M., Romanova, M. M., & Lovelace, R. V. E. 2005, ApJ, 634, 1214

Matt, S., & Pudritz, R. E. 2005a, MNRAS, 356, 167

Matt, S., & Pudritz, R. E. 2005b, ApJL, 632, L135

Matt, S., & Pudritz, R. E. 2007, in IAU Symp. 243, Star-Disk Interaction in Young Stars, ed. J. Bouvier & I. Appenzelle (Cambridge: Cambridge Univ. Press), 299

Matt, S., & Pudritz, R. E. 2008a, ApJ, 681, 391

Matt, S., & Pudritz, R. E. 2008b, ApJ, 678, 1109

Matt, S. P., MacGregor, K. B., Pinsonneault, M. H., & Greene, T. P. 2012, ApJL, 754, L26

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

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Miyoshi, T., Terada, N., Matsumoto, Y., et al. 2010, ITPS, 38, 2236

Pantolmos, G., & Matt, S. P. 2017, ApJ, 849, 83

Pantolmos, G., Zanni, C., & Bouvier, J. 2020, A&A, 643, 129

Rebull, L. M., Wolff, S. C., & Strom, S. E. 2004, AJ, 127, 1029

Regev, O., & Gitelman, L. 2002, A&A, 396, 623

Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015, ApJ, 798, 116

Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2002, ApJ, 578, 420

Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2005, ApJL, 635, L165
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Ud-Doula, A., Owocki, S. P., & Townsend, R. H. D. 2009, MNRAS, 392, 1022

Uzdensky, D. A., Königl, A., & Litwin, C. 2002, ApJ, 565, 1205

Venuti, L., Bouvier, J., Cody, A. M., et al. 2017, A&A, 599, A23

Washimi, H., & Shibata, S. 1993, MNRAS, 262, 936

Weber, E. J., & Davis, L., Jr 1967, ApJ, 148, 217

Zanni, C., & Ferreira, J. 2009, A&A, 508, 1117

Zanni, C., & Ferreira, J. 2011, ApJL, 727, L22

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