Magnetic torques on T Tauri stars: accreting vs. non-accreting systems

Classical T Tauri stars (CTTs) magnetically interact with their surrounding disks, a process that is thought to regulate their rotational evolution. In this work, we compute torques acting onto the stellar surface of CTTs arising from different accreting (accretion funnels) and ejecting (stellar winds and magnetospheric ejections) flow components. Furthermore, we compare the magnetic braking due to stellar winds in two different systems: isolated and accreting stars. 2.5D magnetohydrodynamic, time-dependent, axisymmetric simulations are employed. For both systems the stellar wind is thermally driven. In the star-disk-interaction (SDI) simulations the accretion disk is Keplerian, viscous, and resistive. Two series of simulations are presented, one for each system. We find that in classical T Tauri systems the presence of magnetospheric ejections confines the stellar-wind expansion, resulting in an hourglass-shaped geometry of the outflow. In addition, the formation of the accretion columns modifies the amount of open magnetic flux exploited by the stellar wind. These effects have a strong impact on the stellar wind properties and we show that the stellar-wind braking is more efficient in the star-disk-interacting systems than in the isolated ones. We also derive torque scalings, over a wide range of magnetic field strengths, for each flow component in a SDI system that directly applies a torque on the stellar surface. In all the performed SDI simulations the stellar wind extracts less than 2% of the mass accretion rate and the disk is truncated up to 66% of the corotation radius. All the simulations show a net spin-up torque. In order to achieve a stellar-spin equilibrium we need either more massive stellar winds or disks being truncated closer to the corotation radius, which increases the torque efficiency by the magnetospheric ejections.

Despite being in a phase of stellar contraction and accretion, CTTs are observed to be slow rotators, with rotation periods of 8 days, which corresponds to about (or less than) 10% of their break-up limit (Herbst et al. 2007;Bouvier et al. 2014). In addition, rotation-period distributions from young stellar clusters indicate that CTTs rotate at a constant spin rate with time (Rebull et al. 2004;Irwin & Bouvier 2009;Gallet & Bouvier 2013;Gallet et al. 2019;Amard et al. 2016). The latter features suggest the presence of angular-momentum-loss mechanisms acting on CTTs, which prevent their stellar surfaces to spin up due to both accretion and contraction.
The general consensus is that the angular momentum evolution of CTTs is controlled by the interaction of the stellar magnetic field with the surrounding accretion disk and its environment. The star-disk magnetospheric interaction and the associated outflows have been extensively studied in the literature, using either semi-analytic models (e.g., Ghosh & Lamb 1979;Collier Cameron & Campbell 1993;Lovelace et al. 1995;Armitage & Clarke 1996;Agapitou & Papaloizou 2000;Ferreira et al. 2006;Matt & Pudritz 2005b;Mohanty & Shu 2008;Sauty et al. 2011, and references therein) or numerical simulations (e.g., Hayashi et al. 1996;Goodson et al. 1997;Miller & Stone 1997;Küker et al. 2003;Zanni & Ferreira 2009Čemeljić et al. 2013;Kulkarni & Romanova 2013;Romanova et al. 2013;Cemeljić 2019, and references therein). Within this scenario, several mechanisms have been proposed to explain the rotation periods of CTTs. One of the first star-disk interaction models applied to CTTs was based on the scenario originally proposed by Ghosh & Lamb (1979) for X-Ray pulsars. In this model, the stellar field maintains a connection with the disk outside the corota-Article number, page 1 of 18 arXiv:2009.00940v1 [astro-ph.SR] 2 Sep 2020 A&A proofs: manuscript no. Pantolmos_Zanni_Bouvier_v2_arxiv tion radius, where a Keplerian disk rotates slower than the star. As a consequence, the star can transfer angular momentum to the disk along the field lines connecting the two objects, with the disk rotation controlling the stellar rotation (i.e. the disk-locking mechanism, see e.g., Armitage & Clarke 1996;Agapitou & Papaloizou 2000;Matt & Pudritz 2005b;Zanni & Ferreira 2009). However, this mechanism requires an extended stellar magnetosphere that connects to the disk over a broad region beyond the corotation radius, which seems unlikely (see e.g., Matt & Pudritz 2005b). Different types of outflows, removing the excess stellar angular momentum from the star-disk system instead of transferring it back to the disk, have since been proposed as an alternative solution. Disk winds could effectively remove disk's angular momentum so that the magnetospheric accretion does not affect the stellar angular momentum evolution (e.g., Xwinds, Shu et al. 1994;Cai et al. 2008). Other authors proposed the idea of accretion-powered stellar winds (e.g., Matt & Pudritz 2005a in which magnetospheric accretion acts as an additional energy source to drive the stellar outflow, increasing the mass loss rate and consequently the spin-down efficiency of these winds. Different studies attempted to model the concurrent presence of stellar and disk winds (as a two-component outflow) with the main purpose of investigating the large-scale properties of protostellar jets (see e.g., Bogovalov & Tsinganos 2001;Fendt 2009;Matsakos et al. 2009). Another class of outflows exploits magnetospheric field lines that still connect the star to the disk through a quasi-periodic and unsteady process of inflation and reconnection of these magnetic surfaces. The stellar magnetic surfaces exploited by these outflows can interact with a disk magnetic field whose magnetic moment can be aligned (e.g., ReX-winds, Ferreira et al. 2000) or anti-aligned (e.g., conical winds or magnetospheric ejections; see Romanova et al. 2009;Zanni & Ferreira 2013) with respect to the stellar one. Since these ejections take advantage of magnetic field lines still connecting the star with the disk, they can exchange mass, energy and angular momentum with both of them. In particular the spin-down efficiency of this class of outflows is strongly increased when the disk is truncated close or beyond the disk corotation radius, so that they can efficiently tap the stellar rotational energy (propeller regime Romanova et al. 2005Romanova et al. , 2009Ustyugova et al. 2006;Zanni & Ferreira 2013). On the other hand, the different models of the propeller regime (e.g., Romanova et al. 2005Romanova et al. , 2009Ustyugova et al. 2006;Zanni & Ferreira 2013) tend to predict a strong accretion variability that does not seem to have any observational counterpart.
The aim of this study is to investigate the torques exerted onto a CTTs by different flow components of a star-disk interacting system. We will parametrize these external torques and provide formulae, which have been proven useful for studies that attempt to simulate the rotational evolution of late-type stars (e.g., Gallet & Bouvier 2013Gallet et al. 2019;Johnstone et al. 2015;Matt et al. 2015;Amard et al. 2016Amard et al. , 2019Sadeghi Ardestani et al. 2017;See et al. 2017;Garraffo et al. 2018). Our results will be based on magneto-hydrodynamic (MHD), timedependent, axisymmetric numerical simulations of an accretion disk interacting with the magnetosphere of a rotating star. The accretion disk is taken to be initially Keplerian, viscous, and resistive, modelled with an alpha prescription (Shakura & Sunyaev 1973). Furthermore, we introduce a new equation of state with a temperature-dependent polytropic index, which allows to simulate at the same time a quasi-isothermal stellar wind and an adiabatic disk. Our simulations will model the different flow components that directly apply a torque on the stellar surface: accretion funnel flows that increase the stellar angular momentum; mag-netized stellar winds that provide a spin-down torque; intermittent magnetospheric ejections (hereafter MEs, Zanni & Ferreira 2013), consequence of the differential rotation between the star and the inner disk, which can either spin up or spin down the stellar surface. In particular we will focus on regimes where the disk truncation radius does not exceed the corotation radius, providing a steady accretion flow. In addition, we present stellar wind simulations from isolated stars (e.g., Washimi & Shibata 1993;Keppens & Goedbloed 1999;Matt et al. 2012;Réville et al. 2015;Pantolmos & Matt 2017;Finley & Matt 2018), which are used to compare the stellar-wind torque efficiency in the two different systems (diskless stars vs. star-disk-interacting systems).
In Sect. 2.1 we discuss the numerical method, the initial and boundary conditions of the simulations. In Sect. 3.1, we present the global phenomenology of our numerical solutions, focusing on representative examples of the two systems (i.e., star-disk interaction and isolated stellar winds) simulated here. In Sect. 3.2, we present the the main results of this study. In particular, in Sect. 3.2.2, we derive new stellar-wind torque scalings appropriate for the classical T Tauri phase of stellar evolution, in Sect. 3.2.1 and Sect. 3.2.3, we present torque prescriptions due to accretion and magnetospheric ejections, respectively. In section 4 we compare our results with previous studies from the literature and finally, in section 5 we summarize the conclusions of this work. The details about the equation of state employed in this numerical work are provided in Appendix A and in Appendix B we derive the formulation of the stellar-wind torque discussed in Sect. 3.2.2.

MHD equations and numerical method
The models presented in this work are numerical solutions of the magneto-hydrodynamic (MHD) system of equations, including viscous and resistive effects. In Gaussian units these equations are: (1) The system of Eqs.
(1) consists of mass, momentum and energy conservation equations coupled to the induction equation to follow the evolution of the magnetic field. We indicate with ρ the mass density, P the plasma thermal pressure, B and υ the magnetic field and velocity vectors respectively and I the identity tensor. The total energy E is the sum of internal, kinetic and magnetic energy where the definition of the specific internal energy u(T ) as a function of temperature T is provided in Appendix A. For this work we have specifically developed a caloric equation of state of a calorically imperfect gas, i.e. whose specific heats and their ratio γ are temperature-dependent. In particular the plasma in our models will behave almost isothermally at high temperature and adiabatically at low temperature, so that we will be able to simulate at the same time quasi-isothermal hot stellar winds and cold adiabatic accretion disks. We set the equation of state so that γ = 1.05 for P/ρ > 0.1 GM * /R * and γ = 5/3 for P/ρ < 0.01 GM * /R * , where G is Newton's gravitational constant, M and R are the stellar mass and radius. The stellar gravitational acceleration is given by g = −(GM * /R 2 )R. The electric current is defined by the Ampère's law, J = ∇ × B/4π. We indicate with η m and ν m = η m /4π the magnetic resistivity and diffusivity respectively. The viscous stress tensor T is where η v and ν v = η v /ρ are the dynamic and the kinematic viscosity respectively. Notice that the viscous and magnetic diffusive terms in the right hand side of the total energy equation in the system of Eqs.
(1) correspond to the work of the viscous forces and the diffusion of the magnetic energy only. The dissipative viscous and Ohmic heating terms are not included to avoid, in particular, a runaway irreversible heating of the accretion disk, that would happen since no cooling radiative effects have been taken into account. In the absence of viscosity and resistivity, the system of Eqs. (1) reduces to the ideal MHD equations, which is the system that will be solved for the simulations of stellar winds from isolated stars.
Besides the system of Eqs.
(1), we also solve two passive scalar equations: where s is the specific entropy, whose definition for our newly implemented equation of state is provided in Appendix A, and Tr is a passive tracer. The entropy is used to monitor the dissipation and heating usually associated with the numerical integration of the total energy equation in system of Eqs. (1), e.g. by providing a maximum and minimum entropy value that can be attained during the computation, and it is used to compute the internal energy when the total energy equation provides an unphysical value of the latter. The passive scalar Tr is used to track the disk material and distinguish it from the coronal/stellar wind plasma. We employ a second-order Godunov method provided by the PLUTO code 1 (Mignone et al. 2007) to numerically solve the system of Eqs. (1)-(4). We use a mixture of linear and parabolic interpolation to perform the spatial reconstruction of the primitive variables. The approximate HLLD Riemann solver developed by Miyoshi et al. (2010) is employed to compute the intercell fluxes, exploiting its ability to subtract the contribution of a potential magnetic field (i.e. the initial unperturbed stellar magnetosphere) to compute the Lorentz forces. A second order Runge-Kutta scheme advances the MHD equations in time. The hyperbolic divergence cleaning method (Dedner et al. 2002) is used to control the ∇ · B = 0 condition for the magnetic field. The viscous and resistive terms, computed using a second-order finite difference approximation, have been integrated in time explicitly. We solved the MHD equations in a frame of reference co-rotating with the star. 1 PLUTO is freely available at http://plutocode.ph.unito.it All simulations have been carried out in 2.5 dimensions, i.e. in a two-dimensional computational domain with threedimensional vector fields, assuming axisymmetry around the stellar rotation axis. We solved the equations in a spherical system of coordinates (R, θ). From this point on we will use the capital letter R for the spherical radius and the lower case r = R sin θ to indicate the cylindrical one. Our computational domain covers a region R ∈ [1, 50.76]R * , where R * is the stellar radius, and θ ∈ [0, π]. We discretized the domain with 320 points in the radial direction using a logarithmic grid spacing (i.e. ∆R ∝ R) and with 256 points along θ with a uniform resolution, so that R∆θ ≈ ∆R, i.e. the cells are approximately square.

Initial and boundary conditions
In this work we will present simulations of two different systems: magnetized stars launching thermally-driven winds either (1) isolated or (2) interacting with an accretion disk. For the latter one, the initial conditions are made up of three parts: the accretion disk, the stellar corona, and the stellar magnetic field. For the isolated-stellar-wind simulation the initial conditions are the same, without the presence of a disk.
We set up a Keplerian accretion disk adopting an α parametrization (Shakura & Sunyaev 1973) for the viscosity. Neglecting the inertial terms, its thermal pressure P d and density ρ d can be determined by the vertical hydrostatic equilibrium, while the toroidal speed υ φd can be derived from the radial equilibrium. Assuming a polytropic condition (i.e., P d ∝ ρ γ d with γ = 5/3), we obtained where = c sd /υ K | θ=π/2 is the disk aspect ratio given by the ratio between the disk isothermal sound speed c sd = P d /ρ d and the Keplerian speed υ K = √ GM * /r evaluated at the disk midplane; ρ d0 and υ K * are the disk density and Keplerian speed at the disk midplane at R * . The accretion speed is computed solving the stationary angular momentum equation using an α parametrization for the kinematic viscosity ν v : where Ω K = GM * /r 3 is the Keplerian angular speed. Neglecting the T θφ component of the stress tensor Eq. (3) we obtained This equation shows that the inertial term due to the accretion flow in the radial momentum equation is of order O(α 2 v 4 ), while the thermal pressure gradient is of order O( 2 ), so that it was possible to neglect it when deriving the disk equilibrium Eq. (5). We neglected the T θφ component of the viscous stress tensor in order to avoid the backflow along the disk midplane, most likely unphysical, usually associated with the three-dimensional models of α accretion disks (see e.g., Regev & Gitelman 2002). Consistently, we set T θφ = 0 also during the time-dependent calculations.
We assume that the disk possesses, beside an alpha viscosity, also an anomalous magnetic diffusivity allowing the magnetic flux not to be perfectly frozen into the accretion flow. The idea behind these mechanisms is that some sort of instability, e.g., magneto-rotational or interchange, can trigger a large-scale turbulent transport of angular momentum (viscosity) and magnetic flux (resistivity). For both viscosity ν v and magnetic diffusivity ν m we use a customary alpha parametrization: akin to Eq. (6), where the isothermal sound speed c sν is now space and time dependent. In the outer part of the disk, where the disk structure remains approximately unchanged with respect to the initial disk structure, we fix the sound speed at its initial value provided by Eq. (5), while in the inner part of the disk down to the truncation region, we use the local sound speed value. We use a function F(r) that smoothly goes from zero for r < 0.5R t,i to one for r > 1.5R t,i (where R t,i is the position at which the initial disk solution is truncated, see below) to match the two c s values. We assume that, if the magnetic field is strong enough, the instabilities that trigger the anomalous alpha transport are suppressed. We therefore multiply the sound speed that defines the transport coefficients by an exponential function that goes to zero for µ = B 2 /8πP > 1. In practice, this term cancels the transport coefficients outside the disk, while determining their smooth transition to zero along the accretion funnels. Finally, we multiply the sound speed by a tracer Tr that is set to zero in the stellar atmosphere and to one inside the disk, in order to suppress the viscosity and resistivity in the stellar wind and the magnetic cavity. The full expression for c s in Eq. (8) is therefore: We initialized the stellar atmosphere surrounding the disk computing the thermal pressure and density profiles of a onedimensional, spherically symmetric, isentropic (according to the entropy s defined in Appendix A), transonic Parker-like wind model. This solution is defined by its density ρ * and sound speed c s * at the stellar surface. The poloidal speed is set to zero.
The stellar magnetosphere is modeled initially as a potential dipolar field aligned with stellar the rotation axis. Its two components are where B * is the magnetic field intensity at the stellar equator. The interface between the disk surface and the corona is placed at the position where the disk and coronal thermal pressures are equal. We compute the initial truncation radius R t,i by solving the following implicit equation in the variable R: where B d,θ=π/2 = B * (R * /R) 3 is the intensity of the initial magnetic field and P d,θ=π/2 is the disk thermal pressure both taken at the disk midplane; B + φ is the toroidal magnetic field at the disk surface while M s roughly corresponds to the sonic Mach number of the accretion flow induced by the large-scale magnetic field torque (see e.g. Eq. (3) in Combet & Ferreira 2008). For B + φ we take an estimate of the toroidal field induced at the disk surface by the star-disk differential rotation (see e.g. Collier Cameron & Campbell 1993): where R co = (GM * /Ω 2 * ) 1/3 is the corotation radius, where the disk Keplerian rotation equals the stellar angular speed Ω * . We take M s = 1.5: this means that at its inner radius the disk starts to accrete trans-sonically, which is a typical condition to form the accretion funnels (see e.g. Bessolaz et al. 2008), the initial disk structure Eq. (5) is deeply modified and the local viscous torque is completely negligible (notice, for example, that the typical sonic Mach number induced by the viscous torque is of the order M s = α v , see Eq. 7). Anyway we set an upper limit for the initial truncation radius R t,i < 0.8 R co . Finally, the magnetic surfaces initially threading the disk are set to rotate at the Keplerian angular speed calculated at the disk anchoring radius, while the magnetic surfaces inside R t,i are forced to corotate with the star.
On the rotation axis we assume axisymmetric boundary conditions. At the R = R * boundary we have to consider two types of conditions, one for a subsonic inflow (the stellar wind) and another for a supersonic outflow (the accretion columns). For the subsonic inflow condition we fix in the ghost zones the density and pressure profiles of the one-dimensional wind model used to initialize the stellar atmosphere. For the supersonic outflow conditions the density and pressure must be left free to adjust to the values in the accretion funnel: in this region we use a power-law extrapolation along the magnetic field lines for the density while the pressure is set assuming a constant entropy value along the magnetic surfaces. For intermediate situations (i.e. a subsonic outflow or a hydrostatic corona), we use the sonic Mach number calculated in the first row of cells of the domain to linearly interpolate between the two boundary conditions. Notice that, since the sonic Mach number of the flow can change with time, this boundary condition allows the stellar areas occupied by the accretion spots or the stellar winds to vary in time and adjust to evolution of the system. Clearly, the simulations of stellar winds from isolated stars employ the subsonic inflow conditions only, i.e. a fixed pressure and density profile. The boundary conditions for all the other quantities are the same in both cases. The radial component of the magnetic field is kept fixed, so as to conserve the total stellar flux. The θ component is left free to adjust using a linear extrapolation. The polodial speed υ p is set to be parallel to the poloidal magnetic field, using the conservation of the axisymmetric ideal MHD invariant k = ρυ p /B p along the field lines to determine its value. This condition ensures a smooth inflow for the stellar wind injection while the supersonic infalling accretion funnels are absorbed by the inner boundary without generating a shock. We follow Zanni & Ferreira (2009 to impose a boundary condition that ensures that the magnetic field is frozen into the rotating stellar surface, that is to say that the poloidal component of the electric field in the rotating frame of reference must be equal to zero, i.e. (υ φ − rΩ * )B p − υ p B φ = 0. In order to achieve that, we extract the radial derivative of the toroidal field B φ from the angular momentum conservation equation in the system of Eqs. (1): In order to use this derivative to linearly extrapolate the value of the toroidal field in the boundary zones, we compute a finitedifference approximation of the right-hand side of Eq. (13) calculated in the first row of cells in the computational domain adjacent to the inner boundary. We employ a first-order approximation for the spatial derivatives while for the local toroidal acceleration we use the expression where ∆t is the Alfvén crossing time of a grid cell at the inner domain. This condition ensures that the Lorentz force at the stellar boundary tries to force the magnetic surfaces to rotate at a rate Ω * on a timescale ∆t. Consistently, the boundary value of the toroidal speed is set to At the outer radial boundary a power-law extrapolation for density and pressure is used and all the other variables are linearly extrapolated. Particular attention is devoted to the boundary condition for the toroidal field in the region where the stellar wind exits the computational domain. We used an approach similar to the one employed at the inner radial boundary, only using a much longer timescale ∆t of the order of the Alfvén crossing time of the entire computational domain. This condition avoids artificial torques exerted on the star even when the matter crosses the outer boundary at sub-Alfvénic speeds.

Units and normalization
Simulations have been performed and results will be presented in dimensionless units. The stellar radius, R * is employed as the unit of length. Given the stellar mass, M * , the velocities can be expressed in units of the Keplerian speed at the stellar radius, υ K * = √ GM * /R * . The unit time is t 0 = R * /υ K * . Using the density of the stellar wind at the stellar surface, ρ * , as the reference density, magnetic fields are given in units of Mass-accretion/outflow rates and torques are expressed in units ofṀ 0 = ρ * R 2 * υ K * , and τ 0 = ρ * R 3 * υ 2 K * respectively. Using ρ * = 10 −12 g cm −3 , R * = 2R , M * = 0.7M as reference values, we obtain

Parameters of the study
Two sets of numerical simulations are presented in this work. The first set includes five simulations of isolated stellar winds (hereafter "ISW"), and the second one consists of five simulations of star-disk-interacting systems (hereafter "SDI").
Once the MHD equations and the initial conditions have been normalized, the ISW simulations depend on three dimensionless free parameters, the stellar rotation rate, given as the fraction of the break-up speed f * = R * Ω * /υ K * , the magnetic field intensity B * /B 0 and the wind sound speed at the stellar surface c s * /υ K * . The SDI simulations require four additional parameters to define the disk structure, its density ρ d0 /ρ * , the aspect ratio and the transport coefficients α v and α m . In both sets we fix the stellar rotation taking f * = 0.05, which corresponds to a stellar rotation period of and a corotation radius R co = 7.37 R * . The stellar wind sound speed at the stellar surface is assumed to be c s * = 0.35 υ K * , which corresponds to a specific enthalpy h * = 1.38 υ 2 K * . To define the initial disk structure in the SDI simulations we take in all cases ρ d0 = 100 ρ * , = 0.075 and α v = α m = 0.2. With this choice of parameters the initial disk accretion rate, determined by the viscous torque only, iṡ In the two sets of simulations we varied the stellar magnetic field strength only. The value of B * /B 0 for all the simulations presented in this work is listed in the 2nd column of Table 1.

Numerical Solutions of star-disk interaction
We ran all the SDI simulations for up to 30 stellar periods. After a short strong transient (2-3 stellar periods) the simulations reach a quasi-stationary state, where the integrated mass and angular momentum fluxes slowly vary in time around a mean value. All the SDI cases show a longer-term decrease of the accretion rate, with this effect being more prominent for higher magnetic field cases. This decrease leads the truncation radius getting closer to R co and consequently, simulations having B * > 3.25 B 0 start to display a strong variability most likely associated with a weak propeller regime. We decided to exclude these later highly variable phases from our analysis and focus on the steadily accreting stages. However, for all the SDI cases we based our analysis on timescales longer than 10 stellar periods.
An example of the SDI numerical solutions obtained in this study is shown in Fig. 1. Different groups of field lines can be distinguished: 1) polar magnetic field lines, anchored at the stellar surface that remain open over the entire simulation; 2) open field lines threading the disk beyond the corotation radius, with R co = 7.37 R * ; 3) closed field lines threading the disk in the vicinity of the truncation radius R t 4R * , steadily connecting the disk and the star, and maintaining approximately their geometry throughout the duration of the simulation; 4) closed field lines connecting the star and the disk beyond R t and within R co , periodically evolving through phases of inflation, reconnection,  Notes. (a) Isolated-stellar-wind (ISW) simulations. (b) Star-disk-interaction (SDI) simulations. In our simulations, a negative (positive) torque indicates angular momentum flowing towards (away from) the star (see also §3.1.1). Consequently, τ sw > 0, τ acc < 0, and τ me ≶ 0 (for the sign of τ me in each SDI case see also Fig. 10).  (colormaps) showing the temporal evolution of a star-disk-intercation simulation. The snapshots were taken after 5, 10, 20, and 30 stellar rotation periods, P * . The far-left panel illustrates the inner domain of this simulation, with the vertical white lines indicating the location of the truncation radius R t (dashed black core) and corotation radius R co (solid black core). Each plot includes magnetic field lines (white lines) and velocity-field vectors. The cyan lines delimit the stellar-wind flux tube and the red lines mark the Alfvén surface. and contraction; 5) stellar closed magnetic loops, located below R t , forming a dead zone (or magnetic cavity).
Each group of magnetic field lines is associated with a different type of plasma flow observed in each SDI solution. A conical magnetized stellar wind emerges along the open stellar field lines exerting a braking torque on the star. A disk-wind flows along the open field lines attached to the accretion disk, removing angular momentum from the disk and adding a torque to the viscous one to determine the accretion rate beyond corotation. Around R t , matter is lifted from the disk to form accretion funnel flows. Through this process the star and the disk directly exchange angular mometum. Finally, intermittent ejections propagate within the area neighboring the stellar-and disk-wind open magnetic surfaces, known as magnetospheric ejections (hereafter MEs; see e.g., Zanni & Ferreira 2013). Such outflows occur due to the differential rotation between the star and the disk, which results in the growth of toroidal field pressure. This process leads to the inflation of the field lines attached close to R t that eventually will reconnect, producing plasmoids that propagate ballistically outwards. Clearly, since they occur in a low plasma β = 8πP/B 2 region, where ν m = 0, these reconnection events are numerically driven and therefore depend both on the grid resolution and the numerical algorithm employed. On one hand we are rather confident that numerical effects do not have a strong impact on the launching mechanism of the MEs since, from the point of view of the disk, MEs are launched as magneto-centrifugal flows for which reconnection is not an acceleration driver. On the other hand, magnetic reconnection can modify some large-scale properties of the MEs. For example, numerically-driven magnetic reconnection can modify the stellar magnetic torque due to MEs, since the star can exchange angular momentum only with parts of the system to which it is magnetically connected and once the plasmoids disconnect they can not contribute to the stellar torque anymore. Numerical dissipation can also modify the position of the reconnection X-point and the periodicity of the reconnection events, which should be proportional to the beating frequency (i.e. the difference between the stellar and disk rotation frequency at R t ). Clearly, MEs are the part of our solutions which can be more sensible to numerical effects. While the MEs dynamical picture is physically sound, as confirmed by the number of independent studies presenting an analogous phenomenology (see e.g., Hayashi et al. 1996;Romanova et al. 2009), some quantitative aspects such as the mass and angular momentum stellar fluxes associated with them must be taken with caution. The main objective of this work is to quantify and model the contribution of all these type of flows to the stellar angular momentum. Therefore, we compute their global properties, in particular their mass fluxṀ, and angular momentum flux τ, using, respectivelẏ For the SDI simulations, these integrals are calculated at the stellar surface, separating it into the areas corresponding to the different flow components, using different subscripts for stellar winds (sw), magnetospheric ejections (me) and accretion (acc). We ensure that the sum of the mass and angular momentum fluxes computed for each flow component corresponds to the total flux crossing the inner boundary of our domain. It should be noted that the integrals computed inside the magnetic cavity do not contribute to the mass and angular momentum fluxes, at least in a time-averaged sense. We adopt a sign convention for integrals Eqs. (19), (20) so that a positive (negative) value ofṀ, τ denotes mass/angular momentum flowing away from (towards) the star. Finally, the unsigned magnetic flux is defined as: This quantity can be defined for the full stellar surface (Φ * ) or for the open flux carried by the stellar wind only (Φ open ).

Numerical solutions of isolated stellar winds
Each ISW simulation is stopped when it becomes quasistationary. For this set of simulations, a steady state is achieved after 2-3 stellar periods. In Fig. 2, we present an example of the ISW solutions obtained in this work. Such simulations are less dynamic, compared to the SDI solutions shown above, and only two regions can be identified in the plot: a stellar wind and a dead zone. Clearly, there is a main difference between the two stellarwind solutions from the two different systems (i.e., ISW vs. SDI) studied here. As it can be seen in Fig. 2, asymptotically the flow entirely opens the stellar magnetosphere, filling with plasma the whole domain. On the other hand, the SDI stellar-wind solution, illustrated in Fig. 1, is confined within a conical flux tube due to the presence of MEs. As we will show later, this difference in the expansion of the two stellar outflows has a significant effect on their magnetic torque efficiency. Analogously to SDI simulations, Eqs. (19) and (20) are used to determine the mass and angular momentum fluxes of the stellar wind, and Eq. (21) to compute both the open and total magnetic flux. Since the ISW simulations converge to a steady state, these integrals are also averaged in time and space, using spherical sections of the open flux tubes at different radii, in order to reduce noise and errors.
We recall here that Eq. (20) can be rewritten as where Λ is the total specific angular momentum carried away by the stellar wind, For axisymmetric, ideal MHD, steady-state flows, Λ is invariant along magnetic surfaces and, in the case of a trans-Alfvénic flow, is equal to where r A is the wind Alfvén radius, the radial distance at which the stellar outflow reaches the local Alfvén velocity (e.g., Weber & Davis 1967;Mestel 1968Mestel , 1999. Note that, r A represents the distance from the rotation axis or, in other words, the cylindrical Alfvén radius, r A = R A sin θ A , where R A is now the spherical Alfvén radius and θ A the angular distance of the Alfvén point from the rotation axis. The term r A /R * defines a dimensionless lever-arm that determines the efficiency of the braking torque acting on the star.

Magnetospheric accretion torque
The truncation, or magnetospheric, radius R t of a star-diskinteracting system is commonly parametrized as (see e.g., Bessolaz et al. 2008;Zanni & Ferreira 2009;Kulkarni & Romanova 2013), where R A is the Alfvén radius of a spherical free-fall collapse (e.g., Lamb et al. 1973;Elsner & Lamb 1977), given as and K A is a dimensionless constant that parametrizes the different geometry and dynamics of disk accretion compared to a free-fall collapse. For example, Bessolaz et al. (2008); Zanni & Ferreira (2009) showed that K A can be expressed as a function of the β = 8πP/B 2 parameter and the sonic Mach number of the A&A proofs: manuscript no. Pantolmos_Zanni_Bouvier_v2_arxiv accretion flow in the truncation region. Nevertheless, K A turns out to be always of order unity (see e.g., Kulkarni & Romanova 2013;Zanni & Ferreira 2013).
In the present work, we extract the position of the truncation radius from the simulations by taking the first closed magnetic surface that envelopes the accretion funnels and look for its intersection with the accretion disk. This is achieved by searching for the radial position of the entropy minimum along this magnetic field line, since the bulk of the disk is characterized by the minimum entropy in the whole domain. This criterion provides the radius at which the accretion flow of the disk starts to be deviated and uplifted to form the accretion columns (see Fig.  1). Another common approach (see e.g., Romanova et al. 2002;Kulkarni & Romanova 2013) is to look at the position where the magnetic energy equals the total (thermal plus kinetic) energy of the disk. This definition provides the position at which the midplane accretion flow of the disk is completely disrupted and reasonably provides a slightly smaller estimate of the truncation radius compared to our method.
As discussed in Sect. 3.1.1, in the five SDI simulation presented here, R t /R * slowly increases with time, with variations of ∼ 10% between the lowest and highest value. The time-averaged R t /R * , for all the cases of this work, is given in the 9th column of Table 1.
Defining an Υ-like parameter for the accreting flow, equation (25) can now be rewritten as The absolute values of the mass accretion rates and Υ acc for all the simulations of this study are listed in the 7th and 8th column of Table 1. The dependence of R t /R * on Υ acc is shown in Fig. 3. We fit the data with a power-law in the form of where K t , m t are dimensionless fitting constants. This expression reduces to Eq. (28) for m t = 2/7. The values of K t and m t , for the best fit (solid line) in the plot, are K t = 1.1452 and m t = 0.35, also given in Table 3. From Eq. (28), we get K A ≈ 0.56, which is in agreement with previous numerical studies (Long et al. 2005;Zanni & Ferreira 2009Kulkarni & Romanova 2013). The power index, m t , is found to be 20% higher than its theoretical value, m th t = 2/7 ≈ 0.286. This could be the consequence of the magnetosphere being compressed by the accreting matter, which leads to a slower decline of B * with R (see e.g., Fig. 8 in Zanni & Ferreira 2009). In addition, m t differs by almost a factor of two compared with the value of m t = 0.2 reported in Kulkarni & Romanova (2013). However, their simulations focused on nonaxisymmteric dipolar fields and perhaps further investigation on the dependence of K t and m t on the inclination of a stellar magnetosphere is required.
For a Keplerian disk, the specific angular momentum at R t is l = √ GM * R t . Then, the magnitude of the accretion spin-up torque, as angular momentum is transferred to the star, can be written as where K acc is a dimensionless constant whose value is related to the disk rotational profile in the truncation region. The negative sign indicates that the angular momentum is carried towards the stellar surface (see also §3.1.1). Different studies (Long et al. 2005;Kluźniak & Rappaport 2007;Zanni & Ferreira 2009 have shown that the disk is likely to become sub-Keplerian below the corotation radius. There are two processes torquing down the disk in this region. First, the small annulus of the disk that is threaded by the stellar magnetosphere and corresponds to the base of the accretion columns is subject to a magneticbraking torque by the stellar rotation trying to force the matter to corotate with the star. Through this mechanism, the disk directly transfers its angular momentum to the star so that, even if this process can modify the Keplerian rotation profile, it should not have a strong effect on the value of K acc . Second, the region of the disk which corresponds to the base of the MEs is subject to an external torque since MEs can directly remove angular momentum from the surface of the accretion disk. If this torque is strong enough, it can yield a sub-Keplerian rotation, with the excess angular momentum being ejected, instead of being transfered to the star. Therefore this process usually determines a value of K acc less than unity. At this point it should be noted that the torque provided by these mechanisms can explain the increase ofṀ acc with B * observed in our simulations. The magnetic torques are the dominant drivers of mass accretion in the star-disk-interaction region, as confirmed by the fact that the mass accretion rates measured in our models (see Table 1) are sensibly larger than the accretion rate determined by the viscous torque only, see Eq. (18). As a consequence, a higher B * leads to a stronger overall magnetic torque acting on the disk, which in turn leads to a higher accretion rate. On the other hand, it is not guaranteed that the torques driving accretion in the truncation region are matched by the accretion drivers in the outer disk (viscous and disk-wind torques). Therefore, the long-term decrease of the accretion rate and the corresponding increase of the truncation radius, observed in many of our simulations, is likely due to this mismatch.
The absolute value of the accretion torque τ acc versus the quantityṀ acc (GM * R t ) 1/2 is plotted in Fig. 4. In the figure, the best fit (solid line) provides K acc = 0.79 (see also Table 3). G. Pantolmos, C. Zanni, and J. Bouvier: Magnetic torques on T Tauri stars: accreting vs. non-accreting systems  . Normalized effective lever arm, r A /R * as a function of the wind magnetization, Υ for all the cases of this work. Each data point in the plot represents a single simulation. The blue and black colored points/fitting curves correspond to star-disk-interaction (SDI) and isolated-stellar-wind simulations (ISW) respectively. Data points with the same symbol, connected with a black dashed lines, correspond to numerical solutions having the same surface magnetic field strength, B * . By combining Eqs. (29) and (30), we obtain the accretion torque formulation

Stellar-wind torque
Combining Eqs. (19), (22) and (24), the torque due to a magnetized stellar wind (see e.g., Schatzman 1962;Weber & Davis 1967;Mestel 1968Mestel , 1999 can be written as    (32), provides the definition of the effective lever arm r A , which is calculated using Article number, page 9 of 18 A&A proofs: manuscript no. Pantolmos_Zanni_Bouvier_v2_arxiv S sw scales with radial distance as S sw ∝ R n , where n ≈ 3 indicates super-radial expansion, n ≈ 2, and n < 2 indicate radial and sub-radial expansion, respectively. In both panels the colors are the same as in Fig. 5.
The 6th column of Table 1 lists all the values of r A /R * for this study. Following Matt & Pudritz (2008), we look for a scaling of the lever arm r A with the wind magnetization, Υ, defined as where Φ * is the total unsigned surface magnetic flux and υ esc is the escape speed from the surface of the star.
For all the simulations of this study, parameters Υ and Υ open are tabulated in the 3rd and 5th column of Table 1, respectively. As explained in greater detail in Appendix B, simple powerlaws are usually expected and employed to parametrize the relation between r A /R * and Υ or Υ open . The dependence of r A /R * on the wind magnetization, Υ, is shown in Fig. 5. The following function is used to fit the points where K sw,s , m sw,s are dimensionless fitting constants. Two scaling laws are shown in the plot, which correspond to the two different sets of simulation studied here (i.e., ISW and SDI simulations). The values of K sw,s , m sw,s are given in Table 2 2 . Figure 5 shows that for a given value of Υ, a stellar wind originating from a star-disk-interacting system has a larger braking lever arm. The dependence of r A /R * on Υ open , for all the cases of the study, is presented in Fig. 6. The data points are fitted with a power-law function  Table 2 3 . For clarity, the fitting constants derived from the SDI simulations are also listed in Table 3. As in Fig. 5, for a given Υ open value the stellar wind in a star-disk-interacting system provides a larger magnetic lever arm, even if in this case the offset between the two curves is clearly smaller.
We try now to understand the reason behind the SDI and ISW differences. We verified that the presence of an accretion disk, the accretion columns and the MEs modify the properties of an isolated stellar wind in two main ways.
On one hand we observed that, in the SDI simulations, the amount of open flux that can be exploited by the stellar winds can change substantially. In Fig. 7, we plot the fractional open flux, given in the 4th column of Table 1, versus parameter Υ, for all the numerical solutions. In the plot, we identify two trends. with the Υ parameter in SDI systems is less straightforward. In our SDI cases, the increase of the Υ parameter corresponds also to an increase of the position of R t and, correspondingly, of the Υ acc parameter (see Table 1). A larger truncation radius tends to displace the accretion spots to higher latitudes and therefore to reduce the amount of fractional open flux. With our limited set of simulations it is not possible to determine which effect, the stellar wind push (i.e. the Υ value) or the position of the truncation radius and the accretion spots (i.e. the Υ acc value), is more important to determine the scaling of the open magnetic flux. In particular, since in our simulations we changed the dipolar field intensity B * only, the Υ and Υ acc parameters increase (decrease) at the same time, both contributing to the decrease (increase) of the fractional open flux. Clearly a wider parameter space exploration, with Υ and Υ acc varying independently, is required to address this issue.
On the other hand, we already mentioned that the presence of MEs confines the stellar wind inside an hourglass-shaped magnetic flux tube, instead of letting the wind expand freely as in the ISW cases. To better quantify this effect we plot in the right panel of Fig. 8 the area of the stellar wind flux tube S sw versus the radial distance from the star for two SDI and ISW cases with the same B * /B 0 value. Obviously, the radial dependence of S sw reflects the magnetic field topology since, due to magnetic flux conservation, S sw ∝ B −1 sw . In both cases the magnetic field close to the star (R 4R * ) tends to keep its dipolar potential topology, providing S sw ∝ R n , with n = 3. Note, however, that at the stellar surface the area of the stellar wind is slightly larger in the SDI case than in the ISW case, confirming that the SDI case is characterized by a larger open magnetic flux. It should be mentioned that the areal difference at the stellar surface between the two stellar-wind regions corresponds to a θ-coordinate difference of a few degrees on each hemisphere. Nevertheless, as shown in Fig. 7, such a small difference is capable to increase the open flux by ∼ 20% For low values of the wind magnetization (with Υ < 10 5 ), these angles can differ by up to 10 • on each hemisphere, which results in almost a factor of two increase of Φ open /Φ * for SDI cases. On a large scale (R > 10R * ), the ISW completely opens the magnetosphere, fills the entire domain and propagates almost radially, with n = 2. The stellar wind of the SDI cases remains confined within a smaller area that expands sub-radially, with n < 2.
The flux tube actually acts as a nozzle and strongly modifies the acceleration profile of the flow. Since we are dealing with mainly thermally driven winds (rotation and Lorentz forces are almost negligible), the faster expanding nozzle, the ISW one, should determine a faster decline of the thermal pressure and a faster acceleration of the wind. This is confirmed by the left panel of Fig. 8, where the radial profiles of the wind velocity, the slow-magnetosonic and Alfvén speeds of the two cases are plotted. These quantities are averaged in time and space (along the θ coordinate). Within a few stellar radii the three speed profiles are quite similar in both cases and the two outflows reach the slowmagnetosonic point around R ≈ 4R * . Confirming our hypothesis, on a larger scale the ISW case accelerates more rapidly and becomes noticeably faster than the SDI case. Moreover, since the magnetic field is compressed within a smaller area, the SDI Alfvén speed becomes larger than the ISW one. As a consequence, while the ISW becomes super-Alfvénic at a distance R ≈ 40R * , the stellar wind of the SDI case has not reached the Alfvén point within our computational domain. Even in SDI tests performed with a radial domain twice as large, the stellar wind did not reach the Alfvén surface. Notice that a large difference of the radial distance from the star of the Alfvén surface does not automatically reflects a comparable difference of average lever arm r A . In Appendix B we propose a simple relation between the cylindrical Alfvén radius r A and the average radial distance of the Alfvén surfaceR A , where θ oA is the opening angle of the Alfvén surface. Since for an ISW θ oA ≈ π/2, the lever arm difference is strongly reduced by the smaller opening angle of the SDI stellar winds. Following the approach adopted in Pantolmos & Matt (2017, see also Kawaler (1988; Tout & Pringle (1992); Matt & Pudritz (2008); Réville et al. (2015)), we tried to quantify the impact of these effects (the different amount of open flux and speed profile) and derived a simple analytic expression for the ratio of the lever arms r A in two SDI and ISW cases (see Appendix B for a full derivation): where υ sw,A is the average speed of the flow at the Alfvén surface and quantifies the acceleration efficiency of the stellar wind. According to Eq. (39), the larger lever arm displayed by an SDI stellar wind with the same magnetization Υ of an ISW (see Fig.  5), can be ascribed to the larger open flux and to the slower wind speed of the SDI case. If we compare two SDI and ISW cases with the same Υ open (see Fig. 6), the larger SDI lever arm should be determined by the slower speed (i.e. υ sw,A ) only. Note that this last point is hard to prove, since for the majority of our simulations, we cannot extract υ sw,A because the stellar outflows stay on average sub-Alfvénic within our computational box. Therefore we cannot verify that the offset in Fig. 6 is entirely due to differences in υ sw,A . More subtle geometrical effects could produce an additional scatter of the data points. In addition, the different m o exponent found in ISW and SDI systems could reflect the different wind acceleration of the two sets of simulations. As discussed in Appendix B, the flatter speed profile observed in SDI simulations (a smaller q value in Eq. (B.7)) should correspond to a larger m o exponent compared to ISWs, qualitatively consistent with our findings.
We evaluate the mass ejection and torque efficiencies of the stellar winds by plotting in Fig. 9 the time evolution ofṀ sw /Ṁ acc and τ sw /τ acc . The 10th column of Table 1 lists the time-averaged values of τ sw /τ acc . On average, the stellar outflow is able to extract between 0.2% (i.e., SDI case 5 with B * = 13 B 0 ) and 1% (i.e., SDI case 1 with B * = 1.625 B 0 ) of the mass accretion rate. In addition, the stellar wind braking torque ranges from about 20% (i.e., SDI case 5 with B * = 13 B 0 ) to 40% (i.e., SDI case 2 with B * = 3.25 B 0 ) of the accretion torque. These values are consistent with the results of Zanni & Ferreira (2009. We close this section by presenting the functional form of the torque scaling for stellar winds. By combining Eqs. (36) and (32) or Eqs. (37) and (32), we obtain where Eq.  Table 2, equation Eqs. (40) and (41) provide an estimate of the magnetic torque due to stellar winds over a wide range of surface magnetic field strengths, both for accreting and non-accreting stars.

Magnetospheric-ejections torque
As pointed out in Zanni & Ferreira (2013), MEs can also directly exchange angular momentum with the star. In that case, the direction of the angular-momentum flow depends on the differential rotation between the star and the MEs. If the plasma located at the cusp of the inflated field lines rotates faster than the stellar rotation, the star is subject to a spin-up torque. Similarly, the star experiences a spin-down torque if the matter rotates slower than the star. In order to parametrize the torque exerted by the MEs directly onto the star, we adopt the prescription introduced in Gallet et al. (2019), Fig. 11. τ me /τ acc versus time (given in P * ). Colors and linestyles have the same meaning as in Fig. 9. In the plot, a negative (positive) ratio indicates a spin-down (spin-up) torque due to magnetospheric ejections. where K me and K rot are dimensionless fitting coefficients. Equation (42) assumes that the launching region of the MEs is located close to R t , so that the torque depends on the local magnetic field strength and the differential rotation between the truncation region and the star (for a more detailed discussion see Sect. 2.3.2 in Gallet et al. 2019). In particular, the K rot factor takes into account the difference between the MEs rotation and the disk sub-Keplerian rotation around R t . The dependence of τ me on the term inside the curly brackets in the right-hand side of Eq. (42) is presented in Fig. 10. The best fit (solid line) in the plot has K me = 0.13 and K rot = 0.46 (see also Table 3). K rot is found to be less than unity, indicative of the sub-Keplerian rotation rate of the disk around R t , where the plasma is mass loaded along MEs field lines. In addition, K rot and K me differ by 20% and almost a factor of 2, respectively, compared to the values used in Gallet et al. (2019). These differences could be a result of numerical diffusion effects that appear in our SDI simulations. More specifically, at the boundary between the stellar wind and MEs region there are a few steadily open field lines anchored into the stellar surface along which the plasma accretes close to the star and is ejected at a larger distance, nevertheless exerting a spin-down torque onto the star. This effect is likely due to numerical diffusion of the accretion flow into the MEs and stellar wind regions. Despite being steadily open and providing a spin-down torque, we choose to classify these magnetic surfaces not as part of the stellar wind but as belonging to the MEs region, since the conditions at the stellar surface could in principle allow the MEs to accrete onto the stellar surface, while we strictly require that the stellar wind extracts mass from the star. This choice most likely tends to increase the spin-down efficiency of the MEs and, correspondingly, to slightly underestimate the stellar wind torques. As already discussed in Sect. 3.1.1, this is a further indication of the influence of the numerical subtleties on the quantitative properties of our solutions, and the MEs in particular.
Consistently with Eq. (42) the plot shows a change of sign of τ me , going from a negative (i.e., spin up) to a positive (i.e., spin down) value. This transition occurs at R t ≈ 0.6R co . The latter result agrees qualitatively with the findings of Zanni & Ferreira (2013).
The time-averaged values τ me as a fraction of τ acc are tabulated in the 11th column of Table 1 and the temporal evolution of the MEs torque efficiency, τ me /τ acc , is shown Fig. 11. As seen in Table 1, the efficiency of the MEs torque can vary between a 50% spin-up efficiency for the lowest field case and a 10% spin-down efficiency for the highest field case, changing sign in between the two. For the lower field cases (B * < 5B 0 or nominally B * < 500 G) MEs provide a significant contribution (between 20% and 50% of τ acc ) to the spin-up torque exerted onto the star. This regime requires relatively weak dipolar magnetic fields combined with a ratio R t /R co 0.4. Notice that this configuration is compatible with many observations of classical T Tauri stars (see e.g., Johnstone et al. 2014). The average spindown torque for the high field cases is almost negligible, with a maximum spin-down efficiency ≈ 10% in the highest field case (B * = 13B 0 , nominally corresponding to B * = 1.2 kG). On the other hand Fig. 11 shows that the spin-down efficiency of the MEs in the higher field cases is increasing in time, corresponding to a decrease of the accretion rate and the truncation radius moving closer to corotation. This result points to the fact that in order to maximize the spin-down efficiency of the MEs the disk should be truncated close to the corotation radius, possibly entering a propeller regime, as suggested in Zanni & Ferreira (2013).

Discussion
In this section, we discuss the astrophysical implications of our simulations and compare our findings with previous works from the literature. In addition, we refer to possible limitations of our models.
The angular momentum equation of a star rotating as a solid body is given bẏ In the right hand side of Eq. (43), the first term corresponds to the inverse of the characteristic spin-down/spin-up timescale t sdi = −J * /τ tot of the total exernal torque τ tot = τ acc + τ sw + τ me , where J * = k 2 M * R 2 * Ω * is the stellar angular momentum, with k 2 = 0.2 (i.e., mean radius of gyration of a fully convective star). The second term,Ṁ acc /M * ∼ t −1 acc , is linked to the mass-accretion timescale. This term considers the changes of the stellar moment of inertia due to mass accretion and gives t acc 10Myr, for a typicalṀ acc 10 −7 M yr −1 , so that it can be usually neglected. Finally, the third term, −Ṙ * /R * ∼ t −1 KH , refers to the change of the stellar moment of inertia due to the gravitational stellar contraction and is therefore associated with the Kelvin-Helmholtz timescale (see e.g. Bodenheimer 2011). From Collier Cameron & Campbell (1993; Matt et al. (2010), t KH can be written as where T e f f is the photospheric temperature. Equation (43) shows that, in order to achieve a steady stellar rotation (Ω * = 0), the total external torque should provide a net spin-down, a condition that is not verified in any of our simulations, see Table 1, with a characteristic timescale comparable to the Kelvin-Helmholtz one. It should be noted that stellar contraction on its own is able to decrease the rotation period of slow rotators with P * ≈ 8 days at 1Myr down to 2 days during the disk lifetime, taking a median value of 3 Myr, even if later studies propose a disk lifetime for slow rotators of 9 Myr (see e.g. Williams & Cieza 2011;Gallet et al. 2019). Besides, in all the SDI simulations discussed in this work, the star is subject to an external SDI torque that spins up the star providing a characteristic spin-up timescale varying between 0.4 and 1.5Myr (i.e. for SDI cases 5 and 1, respectively). Obviously, the SDI stellar torques in our simulation can only provide an additional spin-up torque to the stellar rotational evolution, they further shorten the spin-up timescale due to contraction and clearly cannot explain the presence of slow rotators with P * ≥ 8 days at about 10 Myr. A possible solution to this problem could be provided by more massive stellar winds, with a higher spin-down torque. We recall that in our SDI models the spin-down torque of a stellar wind ejecting less than 2% of the mass accretion rate corresponds to 20 − 40% of the spin-up accretion torque. Despite not being sufficient to provide an efficient enough spin-down torque, the important result presented in this paper is that, in the range of the parameter space explored, the stellar wind torque in accreting systems is more efficient than in isolated stars. We found that the presence of the accretion disk strongly influences the amount of open flux and the shape of the wind flux tube, two main factors that determine the wind magnetic lever-arm. We found that, for the same dipolar field intensity B * , the stellar wind torque in SDI systems is 1.3 to 2.3 times higher than in ISW systems. For fixed stellar parameters, using the SDI fits from Table 3, Eqs. (40) and (41) show that the stellar wind torque can be increased by a larger wind mass-loss rate. We can actually use these scalings to estimate the wind mass ejection efficiency required to balance at least the accretion spin-up torque in the SDI cases presented in this paper. Using Eq. (40) we find that an ejection efficiencyṀ sw /Ṁ acc 0.05 is needed while Eq. (41) provides an estimateṀ sw /Ṁ acc > 1. The large discrepancy between these two estimates clearly depends on the the different dependence of τ sw onṀ sw , with the torque increasing more rapidly withṀ sw in Eq. (40) than in Eq. (41). This difference is due to the fact that using Eq. (41), i.e. the scaling with the wind open magnetic flux, we suppose that when the wind-mass loss rate increases the open flux does not change, which could correspond to a situation in which the fractional open flux is solely determined by the position of the accretion spot and the truncation radius (i.e. the Υ acc value, see Eq. (29)). In Eq. (40), i.e. the scaling with the total unsigned flux, we implicitly assume that the raise of the wind mass-loss rate increases the fractional open flux, further amplifying the wind lever arm and the torque (see Eqs. (B.9) and (B.10) in Appendix B). This result clearly points to the fact that, in order to have a robust estimate of the wind mass-loss rate needed to provide an efficient enough spin-down torque, it is necessary to explore the parameter space more extensively to have a more quantitative estimate of the dependency of the fractional open flux on the parameters of the system.
Anyway, requiring a high wind mass-loss rate could present some issues. We recall that, in this work, we considered thermally-driven winds, that require the presence of a hot, 10 6 K corona. Observations show that T Tauri stars are Xray active (e.g., Güdel 2007, and references therein), which implies the presence of million-Kelvin coronae and coronal winds during this pre-main-sequence phase of stellar evolution (e.g., Schwadron & McComas 2008). However, Matt & Pudritz (2007) showed that thermally-driven winds should havė M sw 10 −11 M yr −1 for their energetics to be compatible with the X-ray activity of CTTs, which does not seem to be enough to provide an efficient spin-down torque. Therefore, it has been proposed that stellar winds could be driven by other MHD processes (e.g., dissipation of Alfvén waves, Decampli 1981;Scheurwater & Kuijpers 1988;Cranmer 2008Cranmer , 2009, possibly amplified by the impact of the accretion streams onto the stellar surface (Accretion Powered Stellar Winds, Matt & Pudritz 2005a), so as to produce more massive (withṀ sw /Ṁ acc ∼ 0.1) and colder (10 4 K) winds. The latter idea was supported by observations (Edwards et al. 2006), which suggested the presence of cool and massive stellar winds originating from T Tauri stars, and also confirmed by radiative-transfer calculations (Kurosawa et al. 2006). Cranmer (2008) tested the hypothesis of cold stellar winds in CTTs driven by Alfvén waves but found an ejection efficiency below 2% of the mass accretion rate. Despite the fact that the driving and energetics of CTTs stellar winds still remain an open question, our current findings suggest that the torque provided by magnetized stellar outflows can not single-handedly solve the angular momentum evolution problem of CTTs.
An additional spin-down torque could be provided by MEs. However, in this work, MEs provided either a spin-up torque or a weak spin-down torque. Zanni & Ferreira (2013) showed that the spin-down torque by MEs becomes more efficient when the truncation radius approaches the corotation radius. This is confirmed by our results, since we find an increase of the spin-down efficiency as the truncation radius increases. However truncating the disk close to corotation causes the transition to a (weak) propeller regime, characterized by a strong (an order of magnitude or more) and fast (of the order of one stellar period) variability of the mass accretion rate, which is not observed (Costigan et al. 2014;Venuti et al. 2014). The analysis presented in this paper focused on CTTs being in a steady accretion regime and therefore, we leave this investigation for future studies. It should be noted that the scenario of a propeller regime as a solution to explain the rotational evolution of CTTs is supported by both MHD simulations (Romanova et al. , 2009Ustyugova et al. 2006;Zanni & Ferreira 2013) and angular-momentum-evolution studies (Gallet et al. 2019).
It is also important to point out that, since in our models we varied the magnetic field strength B * only, the simulations displayed aṀ acc − B * correlation, withṀ acc increasing with B * , see discussion in §3.2.1, that has no observable counterpart. Obviously, in our simulations the accretion rate can be changed independently of the B * value, by varying for example the disk surface density and/or the α parameter(s).
Finally, our simulations do not consider complex field topologies. Observations of T Tauri stars show their magnetic fields to be multipolar (e.g., dipolar accompanied by strong octupolar components). Complex field geometries can influence the stellar-wind braking torque in two different ways. First, by affecting the location of the accretion hotspot (Mohanty & Shu 2008;Alencar et al. 2012), which can modify the area that ejects the stellar outflow. This effect could have an impact on the stellar-wind acceleration/speed and mass-loss rate, therefore changing the stellar wind torque. Second, by affecting the amount of open flux available to the wind, which has an impact on the length of the magnetic lever arm and consequently, on the braking torque. In particular, stellar winds with single high-order field topologies (e.g., quadrupoles, octupoles) exhibit decreased torque efficiency (for fixed wind energetics) since the magnetic field decays faster with distance (or the wind carries less open flux), which results in smaller magnetic lever arms (e.g., Réville et al. 2015;Garraffo et al. 2016). For mixed topolgies (e.g., superpositions of dipoles and quadrupoles/octupoles), simulations of isolated stellar winds demonstrate that the stellar-wind torque is mainly determined by the lowest order field topology (Finley & Matt 2017. However, supplementary analysis by See et al. (2019) showed that higher order fields (in mixed geometries) can have an effect on stellar-wind torques if the Alfvén surface is reached at a distance from the stellar surface where the higher-order-field strengths have not declined enough to be considered negligible. Therefore, more realistic/complex field topologies shall be considered in future SDI simulations to improve the accuracy of the current stellar-wind torque presciptions.

Conclusions
In this work, we presented 2.5D, MHD, time-dependent, axisymmetric numerical simulations of magnetized rotating stars interacting with their environment. We focused on two different types of numerical solutions and split our simulations in two sets. The first set included 5 numerical solutions of stellar winds coming from isolated stars (ISW), representing weak-line T Tauri and main-sequence stars. The second set included 5 star-diskinteraction simulations (SDI), representing classical T Tauri systems. In both sets we assumed dipolar field geometries, slow stellar rotation with a stellar spin rate corresponding to 5% of the break-up speed and thermally-driven stellar winds. In the SDI simulations the disk was taken to be viscous and resistive, using a standard alpha prescription, with fixed values for the α parameters and initial disk surface density. In each set, the five simulations are characterized by a different field strength. We provided parametrizations for all the external torques (due to stellar winds, magnetospheric ejections, and accretion funnel flows) exerted at the surface of a star mangnetically interacting with its accretion disk. We also compared the magnetic-braking efficiency of the stellar winds in the two distinct systems considered here (i.e., diskless stars and star-disk interacting systems).
The following points summarize the conclusions of this work.
1. In star-disk interacting systems we find power-law scalings of the stellar-wind Alfvén radius with the wind total Υ or open Υ open magnetization (e.g., Matt et al. 2012;Réville et al. 2015;Pantolmos & Matt 2017;Finley & Matt 2018) akin to those found in ISW systems (see Figs. 5,6,and Eqs. 34,35). Within the parameter space considered, our simulations showed that for a given value of Υ or Υ open stellar winds in SDI systems exhibit a larger magnetic lever arm and an overall more efficient spin-down torque compared to stellar winds from isolated stars. We found that the presence of the accretion funnels tends to increase the amount of open magnetic flux carried by the wind, while the magnetospheric ejections confine the stellar wind in a narrower nozzle that modifies the flow acceleration/velocity profile producing slower stellar-wind solutions. We verified that both these factors contribute to increase the lever arm. In our simulation sample the stellar wind torque is from 1.3 to 2.3 times stronger in SDI systems than it is in ISW systems with the same B * . For the relatively low wind mass-loss rates produced by our SDI simulations, corresponding to less than 2%

Appendix A: The equation of state
We assume the thermal equation of state of an ideal gas PV = Nk B T , where P is the pressure, V is the volume occupied by the gas, N is the number of molecules, k B is the Boltzmann constant and T is the temperature. In terms of the mass density ρ = Nm/V, P = ρk B T/m where m is the mean molecular mass.
To simplify the notation, assuming a constant mean molecular mass, we include k B and m in the definition of the "temperature" so that P = ρT . With this notation the temperature T has the c.g.s. dimensions of cm 2 s −2 , i.e. the square of a speed. We write the first law of thermodynamics: where U is the internal energy and S is the entropy, as: where u = U/M and s = S /M, with the gas mass M = Nm = ρV, are the specific, i.e. per unit mass, internal energy and entropy respectively. We also take the caloric equation of state for an ideal gas: where h = u + T is the specific enthalpy and C V and C p are the specific heat capacities at constant volume and pressure respectively. With our definition of temperature, C V and C p are adimensional and the Mayer's relation becomes C p − C V = 1. For a calorically perfect ideal gas C p and C V are constant and equal to C p = γ/ (γ − 1) and C V = 1/ (γ − 1) where γ = C p /C V is their ratio. We model the gas as thermally but not calorically perfect, i.e. its heat capacities can depend on temperature, and we assume a piecewise constant/linear dependency: where we take C p0 = γ 0 /(γ 0 − 1) and C p1 = γ 1 /(γ 1 − 1), with γ 0 and γ 1 specifying the ratio of the heat capacities for T < T 0 and T > T 1 respectively. Integrating Eq. (A.2) we can write the specific enthalpy as: for T ≤ T 0 0.5 C p1 − C p0 X 2 + C p0 X (T 1 − T 0 ) + C p0 T 0 for T 0 < T < T 1 C p1 T + 0.5 C p0 − C p1 (T 1 + T 0 ) for T ≥ where γ = C p /(C p − 1) = C p /C V . From the first law of thermodynamics Eq. (A.1) we can also derive the expression for the specific entropy s: where s 0 is an integration constant and f (T ) is: Since the specific entropy in the absence of irreversible processes (dissipative heating, cooling) obeys a scalar equation: it is possible to redefine the specific entropy as an arbitrary function of Eq. (A.5). Since the log and exp functions, which would be extensively used to derive entropy from temperature and temperature from entropy respectively, are computationally expensive, we use the simpler expression: In our SDI simulations we assume T 0 = 0.01 (GM /R ), γ 0 = 5/3, T 1 = 0.1 (GM /R ), γ 1 = 1.05.

Appendix B: Analytic scaling of the stellar-wind Alfvén radii
At the Alfvén surface of an MHD flow the wind speed, υ sw,A , equals the local Alfvén speed, υ A whereR A is the average spherical radius of the Alfvén surface and θ oA is the opening angle, measured from the rotation axis, of the outflow flux tube at the Alfvén surface. As discussed in this paper, for ISW cases this angle is about 90 • , since the flow completely opens the stellar magnetic field (see e.g. Fig. 2). For the SDI cases, θ oA is smaller, typically less than 45 • (see e.g. Fig. 1), due to the funnel-shaped geometry of the outflow. If we Article number, page 17 of 18 wherer A represents the average cylindrical radius of the Alfvén surface, considering only its geometrical properties. If we assume thatr A = r A , i.e. the effective Alfvén radius, this equation suggests that S A ∝ r A 2 independently of the opening angle of the wind, either fully open, as in an isolated star, or confined in a smaller funnel, as in SDI systems. Following for example Réville et al. (2015); Pantolmos & Matt (2017), we assume that υ sw,A /υ esc scales as a power law with S A