© G. Pantolmos et al. 2020

1. Introduction

Classical T Tauri stars (CTTs) are young stellar objects (a few Myr old) with M_* ≲ 2 M_⊙ that are surrounded by accretion disks (see e.g., the review by Hartmann et al. 2016). These stars are magnetically active, exhibiting multipolar fields that are typically ∼kG strong (e.g., Johns-Krull 2007; Donati et al. 2008, 2019, 2020; Gregory et al. 2008; Johnstone et al. 2014) and show signatures of mass accretion, with Ṁ_acc varying from 10⁻⁸ to 10⁻¹⁰ M_⊙ yr⁻¹ (e.g., Gullbring et al. 1998; Hartmann et al. 1998; Herczeg & Hillenbrand 2008; Ingleby et al. 2014; Venuti et al. 2014; Alcalá et al. 2017). The measured stellar magnetic fields are strong enough to disrupt the disk and channel the accretion flow into funnels that impact the stellar surface at near free-fall speed, forming hot accretion spots.

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). These features suggest the presence of angular-momentum-loss mechanisms acting on CTTs, which prevent their stellar surfaces from spinning 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 2005a; 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, 2013; Čemeljić et al. 2013; Kulkarni & Romanova 2013; Romanova et al. 2013; Čemeljić 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 (SDI) 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 corotation radius, where a Keplerian disk rotates more slowly 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 2005a; 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 2005a). Different types of outflows, which remove 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 the disk’s angular momentum so that the magnetospheric accretion does not affect the stellar angular momentum evolution (e.g., X-winds, Shu et al. 1994; Cai et al. 2008). Other authors have proposed the idea of accretion-powered stellar winds (e.g., Matt & Pudritz 2005b, 2007), where 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 have 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 magnetic moment. Since these ejections take advantage of magnetic field lines that still connect 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 to or beyond the disk corotation radius, so that they can efficiently tap the stellar rotational energy (propeller regime, Romanova et al. 2005, 2009; Ustyugova et al. 2006; Zanni & Ferreira 2013). On the other hand, the different models of the propeller regime (e.g., Romanova et al. 2005, 2009; Ustyugova 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 the different flow components of an SDI system. We will parameterize these external torques and provide formulae, which have proven useful for studies that have attempted to simulate the rotational evolution of late-type stars (e.g., Gallet & Bouvier 2013, 2015; Gallet et al. 2019; Johnstone et al. 2015; Matt et al. 2015; Amard et al. 2016, 2019; Sadeghi Ardestani et al. 2017; See et al. 2017; Garraffo et al. 2018). Our results will be based on magnetohydrodynamic (MHD), time-dependent, 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, and is modeled with an alpha prescription (Shakura & Sunyaev 1973). Furthermore, we introduce a new equation of state with a temperature-dependent polytropic index, which allows us to simulate a quasi-isothermal stellar wind and an adiabatic disk simultaneously. Our simulations will model the different flow components that apply a torque directly to the stellar surface: accretion funnel flows that increase the stellar angular momentum; magnetized stellar winds that provide a spin-down torque; and intermittent magnetospheric ejections (MEs; Zanni & Ferreira 2013), a 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 and thus provides 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 versus SDI systems).

In Sect. 2, we discuss the numerical method and the initial and boundary conditions of the simulations. In Sect. 3, we present the main results of this study. In Sect. 4, we discuss the astrophysical consequences of our results, and, finally, in Sect. 5 we summarize our conclusions. The details regarding the equation of state employed in this numerical work are given in Appendix A. In Appendix B, we provide the analytic derivation of the formulation for the stellar-wind torques.

2. Numerical setup

2.1. MHD equations and numerical method

The models presented in this work are numerical solutions of the MHD system of equations, including viscous and resistive effects. In Gaussian units, these equations are:

(1)

The system of Eq. (1) consists of mass, momentum, and energy conservation equations coupled with the induction equation to follow the evolution of the magnetic field. We indicate the mass density with ρ, the plasma thermal pressure with P, the magnetic field and velocity vectors with B and υ, respectively, and the identity tensor with I. The total energy E is the sum of internal, kinetic, and magnetic energy:

(2)

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, meaning a gas whose specific heats and γ (i.e., the ratio of the specific heats) are temperature-dependent. In particular, the plasma in our models will behave almost isothermally at high temperatures and adiabatically at low temperatures so that we will be able to simulate quasi-isothermal hot stellar winds and cold adiabatic accretion disks simultaneously. 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 and M_⋆ and R_⋆ are the stellar mass and radius, respectively. The stellar gravitational acceleration is given by . The electric current is defined by Ampère’s law, J = ∇ × B/4π. We indicate the magnetic resistivity and diffusivity with η_m and ν_m = η_m/4π, respectively. The viscous stress tensor 𝒯 is

(3)

where η_v and ν_v = η_v/ρ are the dynamic and the kinematic viscosities, respectively.

It should be noticed that the viscous and magnetic diffusive terms on the right-hand side of the total energy equation in the system of Eq. (1) only correspond to the work of the viscous forces and the diffusion of the magnetic energy. The dissipative viscous and Ohmic heating terms are not included to avoid, in particular, a runaway irreversible heating of the accretion disk, which would happen since no cooling radiative effects have been taken into account. In the absence of viscosity and resistivity, the system of Eq. (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 Eq. (1), we also solve two passive scalar equations:

(4)

where s is the specific entropy, whose definition for our newly implemented equation of state is provided in Appendix A, and T​r 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 Eq. (1), for example by providing a maximum and minimum entropy value that can be attained during the computation. The entropy equation is also used to compute the thermal pressure when the total energy equation provides unphysical negative values of the internal energy. The passive scalar T​r is used to track the disk material and distinguish it from the coronal and stellar-wind plasma.

We employed a second-order Godunov method provided by the PLUTO code¹ (Mignone et al. 2007) to numerically solve the system of Eqs. (1)–(4). We used a mixture of linear and parabolic interpolation to perform the spatial reconstruction of the primitive variables. We employed the approximate HLLD Riemann solver developed by Miyoshi et al. (2010) 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. We used a second-order Runge-Kutta scheme to advance the MHD equations in time. The hyperbolic divergence cleaning method (Dedner et al. 2002) was used to control the ∇ ⋅ B = 0 condition for the magnetic field. The viscous and resistive terms, computed using a second-order finite difference approximation, were integrated in time explicitly. We solved the MHD equations in a frame of reference corotating with the star.

All simulations were carried out in 2.5 dimensions, that is to say, in a two-dimensional computational domain with three-dimensional 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 radius. 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).

2.2. Initial and boundary conditions

In this work we will present simulations of two different systems: magnetized stars launching thermally driven winds that either (1) are isolated or (2) interact with an accretion disk. For the latter, the initial conditions are made up of three parts: the accretion disk, the stellar corona, and the stellar magnetic field. For the isolated stellar winds (ISWs), the initial conditions are the same but without the presence of a disk.

We set up a Keplerian accretion disk adopting an α parametrization (Shakura & Sunyaev 1973) for the viscosity. If we neglect 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., with γ = 5/3), we obtained

(5)

where ϵ = c_sd/υ_K|_{θ = π/2} is the disk aspect ratio given by the ratio between the disk isothermal sound speed and the Keplerian speed evaluated at the disk midplane; ρ_d0 and υ_K* are the disk density and Keplerian speed at the disk midplane at R_*. The accretion speed was computed by solving the stationary angular momentum equation using an α parametrization for the kinematic viscosity ν_v:

(6)

where is the Keplerian angular speed. Neglecting the 𝒯_θϕ component in the stress tensor Eq. (3), we obtained

(7)

This equation shows that the inertial term due to the accretion flow in the radial momentum equation is of the order of , and it can be neglected when deriving the disk equilibrium Eq. (5) since the thermal pressure gradient is of the order of 𝒪(ϵ²). We neglected the 𝒯_θϕ component of the viscous stress tensor in order to avoid the backflow along the disk midplane, most likely unphysical, that is usually associated with the three-dimensional models of α accretion disks (see e.g., Regev & Gitelman 2002). We consistently set 𝒯_θϕ = 0 in the time-dependent calculations as well.

We assumed that the disk possesses, in addition to an alpha viscosity, an anomalous magnetic diffusivity, allowing the magnetic flux to not be perfectly frozen in the accretion flow. The idea behind these mechanisms is that some sort of instability, for example 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 used a customary alpha parametrization,

(8)

that is 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 fixed the sound speed at its initial value provided by Eq. (5); in the inner part of the disk down to the truncation region, we used the local sound speed value. We used a function F(r) that goes smoothly 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 assumed that if the magnetic field is strong enough, the instabilities that trigger the anomalous alpha transport are suppressed. We therefore multiplied the sound speed that defines the transport coefficients by an exponential function that goes to zero for μ = B²/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 multiplied the sound speed by a tracer T​r 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:

(9)

We initialized the stellar atmosphere surrounding the disk with the thermal pressure and density profiles of a one-dimensional, 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 was initially modeled as a potential dipolar field aligned with the stellar rotation axis. Its two components are

(10)

where B_* is the magnetic field intensity at the stellar equator.

The interface between the disk surface and the corona was placed at the position where the disk and coronal thermal pressures are equal. We computed the initial truncation radius R_t, i by solving the following implicit equation in the variable R:

(11)

where: B_{d, θ = π/2} = B_*(R_*/R)³ is the intensity of the initial magnetic field and P_{d, θ = π/2} is the disk thermal pressure (both taken at the disk midplane); is the toroidal magnetic field at the disk surface; and 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 , we took an estimate of the toroidal field induced at the disk surface by the star-disk differential rotation (see e.g., Collier Cameron & Campbell 1993):

(12)

where is the corotation radius (i.e., the radial distance at which the disk Keplerian rotation equals the stellar angular speed Ω_*). We took M_s = 1.5: This means that at its inner radius, the disk starts to accrete transsonically, which is a typical condition for forming the accretion funnels (see e.g., Bessolaz et al. 2008); the initial disk structure from Eq. (5) is deeply modified and the local viscous torque is completely negligible (we notice, for example, that the typical sonic Mach number induced by the viscous torque is of the order of M_s = α_vϵ, see Eq. (7)). 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.

Along the rotation axis, we assumed axisymmetric boundary conditions. At the R = R_* boundary, we had 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 fixed the density and pressure profiles of the one-dimensional wind model used to initialize the stellar atmosphere in the ghost zones. 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 used a power-law extrapolation along the magnetic field lines for the density while the pressure was set assuming a constant entropy value along the magnetic surfaces. For intermediate situations (i.e., a subsonic outflow or a hydrostatic corona), we used the sonic Mach number calculated in the first row of cells of the domain to linearly interpolate between the two boundary conditions. We 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 with time and adjust to the evolution of the system. It is clear that the simulations of stellar winds from isolated stars only employ the subsonic inflow conditions (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 was kept fixed so as to conserve the total stellar flux. The θ component was left free to adjust using a linear extrapolation. The poloidal speed υ_p was 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 followed Zanni & Ferreira (2009, 2013) to impose a boundary condition that ensures that the magnetic field is frozen in 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 − υ_pB_ϕ = 0). In order to achieve that, we extracted the radial derivative of the toroidal field B_ϕ from the angular momentum conservation equation in the system of Eq. (1):

(13)

In order to use this derivative to linearly extrapolate the value of the toroidal field in the boundary zones, we computed a finite-difference approximation of the right-hand side of Eq. (13), which was calculated in the first row of cells in the computational domain adjacent to the inner boundary. We employed a first-order approximation for the spatial derivatives, while for the local toroidal acceleration we used the expression

(14)

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. We consistently set the boundary value of the toroidal speed to

(15)

At the outer radial boundary, a power-law extrapolation for density and pressure was used and all the other variables were linearly extrapolated. Particular attention was 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 that employed at the inner radial boundary, only using a much longer timescale Δt that was taken to be 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.

2.3. Units and normalization

Simulations were 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, . The unit time is t₀ = 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 and outflow rates are expressed in units of . Torques are expressed in units of . Using ρ_* = 10⁻¹² g cm⁻³, R_* = 2 R_⊙, and M_* = 0.7 M_⊙ as reference values, we obtain

(16)

2.4. Parameters of the study

Two sets of numerical simulations are presented in this work. The first set includes five simulations of ISWs, and the second consists of five simulations of SDI systems.

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₀; 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 fixed the stellar rotation taking f_* = 0.05, which corresponds to a stellar rotation period of

(17)

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 . To define the initial disk structure in the SDI simulations, we took ρ_d0 = 100 ρ_*, ϵ = 0.075, and α_v = α_m = 0.2 in all cases. With this choice of parameters, the initial disk accretion rate, determined by the viscous torque only, is

(18)

In the two sets of simulations, we only varied the stellar magnetic field strength. The values of B_*/B₀ for all the simulations presented in this work are listed in the second column of Table 1.

3. Results

In this section, we summarize the outcome of our numerical models. In Sect. 3.1, we discuss the general phenomenology of our solutions by focusing on representative examples of the two systems (i.e., SDI and ISW). In Sect. 3.2, we present the torque scalings. In particular, the accretion torque is discussed in Sect. 3.2.1. In Sect. 3.2.2, we compare the stellar-wind torque in accreting and non-accreting systems. In Sect. 3.2.3, we examine the torque due to MEs.

3.1. MHD simulations

3.1.1. Numerical solutions of star-disk interaction

We ran all the SDI simulations for up to 30 stellar periods. After a short strong transient (two to three stellar periods), the simulations reach a quasi-stationary state, where the integrated mass and angular momentum fluxes slowly vary with time around a mean value. All the SDI cases show a longer-term decrease in the accretion rate, with this effect being more prominent for cases with stronger magnetic fields. This decrease leads to the truncation radius getting closer to R_co and, consequently, simulations that have B_* >  3.25 B₀ start to display a strong variability that is 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 ten 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 their approximate 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, and contraction; (5) stellar closed magnetic loops, located below R_t, forming a dead zone (or magnetic cavity).

Fig. 1. Logarithmic normalized density (colormaps) showing the temporal evolution of an SDI simulation. The snapshots were taken after five, ten, 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 Rt (dashed black core) and corotation radius Rco (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.

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 providing, together with viscosity, the torque that determines 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 momentum. Finally, intermittent ejections propagate within the area neighboring the stellar- and disk-wind open magnetic surfaces, known as 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, which will eventually reconnect and produce plasmoids that propagate ballistically outward. It is clear that, since they occur in a low plasma β = 8πP/B² 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 exerted by MEs since the star can only exchange angular momentum with parts of the system to which it is magnetically connected, and once the plasmoids disconnect, they can no longer contribute to the stellar torque. 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 frequencies at R_t). Clearly, MEs are the part of our solutions that can be most sensitive 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 must be taken with caution, such as the mass and angular momentum stellar fluxes associated with them.

The main objective of this work is to quantify and model the contribution of all these types of flows to the stellar angular momentum. Therefore, we computed their global properties, in particular their mass flux Ṁ and angular momentum flux τ, using, respectively,

(19)

(20)

For the SDI simulations, these integrals were calculated at the stellar surface, separating it into the areas corresponding to the different flow components and using different subscripts for stellar winds (sw), MEs (me), and accretion (acc). We ensured that the sum of the mass and angular momentum fluxes computed for each flow component corresponded 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 adopted a sign convention for integral Eqs. (19) and (20) so that a positive (negative) value of Ṁ, τ denotes mass and angular momentum flowing away from (toward) the star. Finally, the unsigned magnetic flux is defined as:

(21)

This quantity can be defined for the full stellar surface (Φ_*) or for the open flux carried by the stellar wind only (Φ_open).

3.1.2. Numerical solutions of isolated stellar winds

Each ISW simulation was stopped when it became quasi-stationary. For this set of simulations, a steady state was achieved after two to three 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 major difference between the two stellar-wind solutions from the two different systems (i.e., ISW versus SDI) studied here. As can be seen in Fig. 2, the flow entirely opens the stellar magnetosphere asymptotically, filling the whole domain with plasma. 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 later show, this difference in the expansion of the two stellar outflows has a significant effect on their magnetic torque efficiency.

Fig. 2. Same as the full-domain panels of Fig. 1 for a simulation of an ISW. The image was taken at three stellar periods.

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) is used 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

(22)

where Λ is the total specific angular momentum carried away by the stellar wind,

(23)

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

(24)

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 1968, 1999). We note that r_A represents the distance from the rotation axis or, in other words, the cylindrical Alfvén radius, r_A = R_Asin θ_A, where R_A is now the spherical Alfvén radius and θ_A is 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.

3.2. Scalings of magnetic torques

3.2.1. Magnetospheric accretion torque

The truncation, or magnetospheric, radius R_t of an SDI system is commonly parameterized as (see e.g., Bessolaz et al. 2008; Zanni & Ferreira 2009; Kulkarni & Romanova 2013)

(25)

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

(26)

and K_A is a dimensionless constant that parameterizes the different geometry and dynamics of disk accretion compared to a free-fall collapse. For example, Bessolaz et al. (2008) and Zanni & Ferreira (2009) showed that K_A can be expressed as a function of the β = 8πP/B² parameter and the sonic Mach number of the accretion flow in the truncation region. Nevertheless, K_A turns out to always be of order unity (see e.g., Kulkarni & Romanova 2013; Zanni & Ferreira 2013).

In the present work, we extracted the position of the truncation radius from the simulations by taking the first closed magnetic surface that envelops the accretion funnels and looking for its intersection with the accretion disk. This was 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, R_t/R_* slowly increases with time in the five SDI simulations presented here, with variations of ∼10% between the lowest and highest values. The time-averaged R_t/R_* for all cases in this work is given in the ninth column of Table 1.

Defining an Υ-like parameter for the accreting flow,

(27)

Equation (25) can now be rewritten as

(28)

The absolute values of the mass accretion rates and Υ_acc for all the simulations of this study are listed in the seventh and eight columns 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

(29)

Fig. 3. Normalized truncation radius, Rt/R*, as a function of parameter Υacc. Each symbol (or data point) in the plot represents a single SDI simulation. The solid line shows the fitting function (29), with Kt = 1.1452 and mt = 0.35.

where K_t and 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, which is 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 2009, 2013; Kulkarni & Romanova 2013). The power index, m_t, is found to be 20% higher than its theoretical value, . 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 non-axisymmetric dipolar fields, and further investigation on the dependence of K_t and m_t on the inclination of a stellar magnetosphere may be required.

For a Keplerian disk, the specific angular momentum at R_t is . The magnitude of the accretion spin-up torque, as angular momentum is transferred to the star, can therefore be written as

(30)

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 toward the stellar surface (see also Sect. 3.1.1). Different studies (Long et al. 2005; Kluźniak & Rappaport 2007; Zanni & Ferreira 2009, 2013) 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 magnetic-braking torque by the stellar rotation trying to force the matter to corotate with the star. Through this mechanism, the disk transfers its angular momentum directly 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 that 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 transferred to the star. Therefore, this process usually determines a value of K_acc that is 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 SDI region, as confirmed by the fact that the mass accretion rates measured in our models (see Table 1) are significantly larger than the accretion rates 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 in the accretion rate and the corresponding increase in 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).

Fig. 4. Normalized accretion torque (absolute values), τacc, versus quantity Ṁacc(GM*Rt)1/2. Symbols are the same as in Fig. 3. The solid line corresponds to the best fit of the data, using Eq. (30) with Kacc = 0.79.

By combining Eqs. (29) and (30), we obtain the accretion torque formulation

(31)

3.2.2. 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 1968, 1999) can be written as

(32)

where the positive sign in Eq. (32) denotes angular momentum carried away from the star (see also Sect. 3.1.1). In multidimensional solutions, the flow becomes super-Alfvénic along a surface (see e.g., Fig. 2), and therefore ⟨r_A⟩ represents the wind average lever arm in Eq. (32) (see e.g., Washimi & Shibata 1993; Matt et al. 2012; Réville et al. 2015; Pantolmos & Matt 2017; Finley & Matt 2018). In fact, Eq. (32) provides the definition of the effective lever arm ⟨r_A⟩, which is calculated using

(33)

The sixth column of Table 1 lists all the values of ⟨r_A⟩/R_* for this study.

Following Matt & Pudritz (2008), we looked for a scaling of the lever arm ⟨r_A⟩ with the wind magnetization, Υ, defined as

(34)

where Φ_* is the total unsigned surface magnetic flux and υ_esc is the escape speed from the surface of the star. The parameter Υ_open, which depends on the open unsigned magnetic flux Φ_open carried by the stellar wind (see e.g., Réville et al. 2015), is defined as

(35)

Parameters Υ and Υ_open for all the simulations of this study are tabulated in the third and fifth columns of Table 1, respectively.

As explained in greater detail in Appendix B, simple power laws are usually expected and employed to parameterize 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:

(36)

Fig. 5. Normalized effective lever arm, ⟨rA⟩/R*, as a function of the wind magnetization, Υ, for all the cases in this work. As in Fig. 3, each data point is a single simulation. The blue and black colored points (or fitting curves) correspond to SDI and ISW simulations, respectively. Data points with the same symbol, connected with dashed black lines, correspond to numerical solutions that have the same surface magnetic field strength, B*.

where K_sw, s and m_sw, s are dimensionless fitting constants. Two scaling laws are shown in the plot, which correspond to the two different sets of simulations studied here (i.e., ISW and SDI simulations). The values of K_sw, s and m_sw, s are given in Table 2². Figure 5 shows that for a given value of Υ, a stellar wind originating from an SDI system has a larger braking lever arm.

The dependence of ⟨r_A⟩/R_* on Υ_open for all the cases in the study is presented in Fig. 6. The data points are fitted with a power-law function

(37)

Fig. 6. ⟨rA⟩/R* versus parameter Υopen for all simulations in this study. Colors, symbols, and line styles have the same meaning as in Fig. 5.

where K_sw, o and m_sw, o are dimensionless fitting constants. The values of K_sw, o and m_sw, o of the two fits are tabulated in Table 2³. For clarity, the fitting constants derived from the SDI simulations are also listed in Table 3. As in Fig. 5, the stellar wind in an SDI system provides a larger magnetic lever arm for a given Υ_open value, even if the offset between the two curves is clearly smaller in this case.

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 ISW in two main ways.

On one hand, we observed that the amount of open flux that can be exploited by the stellar winds can change substantially in the SDI simulations. In Fig. 7, we plot the fractional open flux, given in the fourth column of Table 1, versus parameter Υ for all the numerical solutions. In the plot, we identify two trends. First, for both sets of simulations, Φ_open/Φ_* decreases as a function of Υ, and second, for a given value of Υ, SDI simulations tend to produce a larger amount of fractional open flux. In the case of ISWs, the decrease in the fractional open flux with Υ can be understood if we look at the Υ parameter as the ratio between the magnetic field energy density and the kinetic energy density of the flow (see e.g., ud-Doula & Owocki 2002; Ud-Doula et al. 2009). The larger the magnetic energy with respect to the stellar-wind push, the harder it is to open up the closed field structure, resulting in a smaller open flux fraction. While in ISW cases the amount of open flux is determined by the stellar wind push only, in SDI systems it also depends on the position of the accretion spots, since the star-disk differential rotation has the tendency to open up all the magnetic surfaces that are not mass-loaded by the columns⁴. While this effect is likely responsible for the larger open flux measured in SDI simulations, the interpretation of the decrease in the open flux with the Υ parameter in SDI systems is less straightforward. In our SDI cases, the increase in the Υ parameter also corresponds to an increase in the position of R_t and, accordingly, to an increase in the Υ_acc parameter (see Table 1). A larger truncation radius tends to displace the accretion spots to higher latitudes and therefore reduces 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 for determining the scaling of the open magnetic flux. In particular, since in our simulations we only changed the dipolar field intensity B_*, the Υ and Υ_acc parameters increase (decrease) at the same time, both contributing to the decrease (increase) of the fractional open flux. It is clear that a wider parameter space exploration, with Υ and Υ_acc varying independently, is required to address this issue.

Fig. 7. Fractional open flux Φopen/Φ* versus Υ. Colors, symbols, and line styles have the same meaning as in Fig. 5.

On the other hand, we have 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 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₀ value in the right-hand panel of Fig. 8. Obviously, the radial dependence of S_sw reflects the magnetic field topology since, due to magnetic flux conservation, . In both cases, the magnetic field close to the star (R ≲ 4R_*) tends to keep its dipolar potential topology, providing S_sw ∝ Rⁿ, with n = 3. We 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 of increasing the open flux by ∼20% For low values of the wind magnetization (with Υ <  10⁵), these angles can differ by up to 10° on each hemisphere, which results in an increase in Φ_open/Φ_* by nearly a factor of two 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.

Fig. 8. Left: profiles of normalized Alfvén (dashed), slow-magnetosonic (dash-dotted), and stellar-wind-poloidal (solid) speeds as a function of radius for two cases presented in this work. The velocities are averaged over the θ coordinate and in time. Right: surface area of the stellar-wind flux tube in normalized units versus R/R*. Ssw scales with radial distance as Ssw ∝ Rn, where n ≈ 3 indicates super-radial expansion, and 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 flux tube actually acts as a nozzle and strongly modifies the acceleration profile of the flow. Since we are dealing with stellar winds whose main driver is the thermal pressure (rotation and Lorentz forces are almost negligible), the faster-expanding nozzle, the ISW nozzle, should determine a faster decline of the thermal pressure and a faster acceleration of the wind. This is confirmed by the left-hand panel of Fig. 8, where the radial profiles of the wind velocity and 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 slow-magnetosonic point around R ≈ 4R_*. Confirming our hypothesis, the ISW case accelerates more rapidly on a larger scale 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 that of the ISW. As a consequence, while the ISW becomes super-Alfvénic at a distance R ≈ 40 R_*, the stellar wind of the SDI case does not reach 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. We notice that a large difference in the radial distance from the star of the Alfvén surface does not automatically reflect 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 surface ,

(38)

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 amounts of open fluxes and the differences in the speed profiles) 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):

(39)

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. We 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 sub-Alfvénic on average 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 exponents 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 the time evolution of Ṁ_sw/Ṁ_acc and τ_sw/τ_acc in Fig. 9. The tenth 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₀) and 1% (i.e., SDI case 1 with B_* = 1.625 B₀) 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₀) to 40% (i.e., SDI case 2 with B_* = 3.25 B₀) of the accretion torque. These values are consistent with the results in Zanni & Ferreira (2009, 2013).

Fig. 9. Absolute values of the mass ejection efficiency Ṁsw/Ṁacc and torque efficiency τsw/τacc as a function of time expressed in units of the stellar rotation period. The straight lines in the plot represent the time-averaged values for the two cases exhibiting the lowest and highest values of these two quantities.

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

(40)

(41)

where Eq. (40) employs the total unsigned surface magnetic flux Φ_*, while Eq. (41) uses the open magnetic flux Φ_open. Using the corresponding values of K_sw, s, K_sw, o, and m_s, m_o given in Table 2, 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.

3.2.3. Magnetospheric-ejection 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 more slowly than the star.

In order to parameterize the torque exerted by the MEs directly onto the star, we adopted the prescription introduced in Gallet et al. (2019),

(42)

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 on 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). We find K_rot 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 ME field lines. In addition, K_rot and K_me differ by 20% and almost a factor of two, 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 ME regions, there are a few steadily open field lines anchored onto the stellar surface, along which the plasma accretes close to the star and is ejected at a larger distance; these open field lines nevertheless exert a spin-down torque onto the star. This effect is likely due to numerical diffusion of the accretion flow into the ME and stellar-wind regions. Despite being steadily open and providing a spin-down torque, we chose to classify these magnetic surfaces not as part of the stellar wind but as belonging to the ME region, since the conditions at the stellar surface could in principle allow the MEs to accrete onto the stellar surface, while we strictly required 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.

Fig. 10. Normalized τme versus quantity in curly brackets in Eq. (42). Symbols are the same as in Fig. 3. In the plot, a negative (positive) τme spins up (down) the stellar rotation. Using Eq. (42), the best fit to the points gives Kme = 0.13 and Krot = 0.46.

Consistent 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. These results agree 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 ME torque efficiency, τ_me/τ_acc, is shown in Fig. 11. As seen in Table 1, the efficiency of the ME 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₀ or nominally B_* <  500 G), MEs provide a significant contribution (between 20% and 50% of τ_acc) to the spin-up torque exerted on the star. This regime requires relatively weak dipolar magnetic fields combined with a ratio R_t/R_co ≲ 0.4. We notice that this configuration is compatible with many observations of CTTs (see e.g., Johnstone et al. 2014). The average spin-down torque for the high field cases is almost negligible, with a maximum spin-down efficiency ≈10% in the highest field case (B_* = 13B₀, 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 increases with time, corresponding to a decrease in 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).

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 MEs.

4. Discussion

In this section, we present the astrophysical implications of our simulations and compare our findings with previous works from the literature. In addition, we discuss the possible limitations of our models.

The angular momentum equation of a star rotating as a solid body is given by

(43)

On the right-hand side of Eq. (43), the first term corresponds to the inverse of the characteristic spin-down (or spin-up) timescale t_sdi = −J_*/τ_tot of the total external torque τ_tot = τ_acc + τ_sw + τ_me, where is the stellar angular momentum, with k² = 0.2 (i.e., the mean radius of gyration of a fully convective star). The second term, , is linked to the mass accretion timescale. This term considers the changes of the stellar moment of inertia due to mass accretion. Therefore, Ṁ_acc/M_* can usually be neglected because t_acc ≳ 10 Myr for a typical Ṁ_acc ≲ 10⁻⁷ M_⊙ yr⁻¹. Finally, the third term, , 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) and Matt et al. (2010), t_KH can be written as

(44)

where T_eff is the photospheric temperature. Equation (43) shows that, in order to achieve a steady stellar rotation (), 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 timescale. 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 two days during the disk lifetime, taking a median value of 3 Myr; however, more recent studies have proposed a disk lifetime for slow rotators of 9 Myr (see e.g., Williams & Cieza 2011; Gallet et al. 2019). Furthermore, 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.5 Myr (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 find 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 also estimate 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 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 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; this 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 rise of the wind mass-loss rate increases the fractional open flux, further amplifying the torque (see Eqs. (B.9) and (B.10)). 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.

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⁶ K) corona. Observations show that T Tauri stars are X-ray 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) have shown that thermally driven winds should have Ṁ_sw ≲ 10⁻¹¹ M_⊙ yr⁻¹ 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., the dissipation of Alfvén waves, Decampli 1981; Scheurwater & Kuijpers 1988; Cranmer 2008, 2009), possibly amplified by the impact of the accretion streams onto the stellar surface (accretion-powered stellar winds, Matt & Pudritz 2005b), so as to produce more massive (with Ṁ_sw/Ṁ_acc ∼ 0.1) and colder (10⁴ K) winds. This idea was supported by observations (Edwards et al. 2006) that suggest the presence of cool and massive stellar winds originating from T Tauri stars and was 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 open questions, our current findings suggest that the torque provided by magnetized stellar outflows cannot 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 in 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 we therefore 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. 2005, 2009; Ustyugova 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 only varied the magnetic field strength B_*, the simulations displayed an Ṁ_acc − B_* correlation, with Ṁ_acc increasing with B_* (see discussion in Sect. 3.2.1), which 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, it can affect 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 and 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 topologies (e.g., superpositions of dipoles and quadrupoles and/or octupoles), simulations of ISWs demonstrate that the stellar-wind torque is mainly determined by the lowest-order field topology (Finley & Matt 2017, 2018). However, supplementary analysis by See et al. (2019) shows 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, complex and more realistic field topologies should be considered in future SDI simulations to improve the accuracy of the current stellar-wind torque prescriptions.

5. Conclusions

In this work, we presented 2.5D, MHD, time-dependent, and axisymmetric numerical simulations of magnetized rotating stars interacting with their environment. We focused on two different types of numerical solutions and split our simulations into two sets. The first set included five numerical solutions of stellar winds coming from isolated stars (ISW), representing weak-line T Tauri and main-sequence stars. The second set included five SDI simulations, 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, MEs, and accretion funnel flows) exerted at the surface of a star magnetically 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 SDI systems).

The following points summarize the conclusions of this work.

1. In SDI 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 find that the presence of the accretion funnels tends to increase the amount of open magnetic flux carried by the wind and the MEs confine the stellar wind to a narrower nozzle that modifies the flow acceleration and velocity profiles, producing slower stellar-wind solutions. We verified that both these factors contribute to increasing the lever arm. In our simulation sample, the stellar-wind torque is 1.3 to 2.3 times stronger in SDI systems than 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% of the mass accretion rate, the stellar wind is able to extract between 20 and 40% of the accreting angular momentum.

2. In the SDI simulations, for the range of field strengths and mass accretion rates considered in this work, the disk is truncated by up to 66% of the corotation radius. We obtained a simple power-law scaling of the truncation radius R_t with the accretion-related magnetization parameter Υ_acc (see Eq. (28)), and we further verified that the spin-up torque due to accretion is linearly proportional to the mass accretion rate times the disk specific angular momentum at the truncation region (see Eq. (31)).

3. Following Gallet et al. (2019), we modeled the magnetic torque exerted by MEs on the star as a differential rotation effect between the star and the MEs, and we derived a torque parametrization over a range of stellar field strengths (see Eq. (42)). In our set of simulations, MEs provide either a spin-up torque, up to ∼50% of the accretion torque, or a weak spin-down torque, up to ∼10% of the accretion torque, depending on the relative position of the truncation radius R_t with respect to the corotation radius R_co. The transition from a spin-up to a spin-down ME torque occurs at R_t ≈ 0.6R_co.

4. In all the SDI cases analyzed, the star is subject to a net spin-up torque. The latter result yields a spin-up timescale due to the sum of all external torques of about 1 Myr, which is even shorter than the spin-up timescale due to stellar contraction (see Sect. 4 for more details). Therefore, the stars simulated in our study would have the tendency to spin-up rather rapidly, and, without the presence of a more efficient spin-down mechanism, they could not explain the existence of slow rotators (with P_* ≳ 8 days) after about 10 Myr. We argue that a stellar wind with a mass-loss rate higher than the one considered in our models could increase the spin-down efficiency. However, a more extensive exploration of the parameter space is necessary in order to provide robust estimates on the stellar-wind ejection efficiency needed to counteract the spin-up torques. Finally, our simulations show that, as the truncation radius gets closer to corotation and with the system moving toward a propeller regime, the spin-down efficiency of the MEs increases, further contributing to the solution of the rotational evolution problem of CTTs. This result agrees with the conclusions of other previous studies (Romanova et al. 2005, 2009; Ustyugova et al. 2006; Zanni & Ferreira 2013; Gallet et al. 2019).

Acknowledgments

The authors thank Catherine Dougados, Sean Matt, Lewis Ireland, and Florian Gallet for helpful discussions during the production of the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 742095; SPIDI: Star-Planets-Inner Disk-Interactions)’; https://readymag.com/gk/spidi/. The numerical simulations presented in this paper were performed with the Dahu supercomputer of the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. All the figures within this work were produced using the Python-library Matplotlib (Hunter 2007).

References

Appendix A: The equation of state

We assumed the thermal equation of state of an ideal gas PV = Nk_BT, 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_BT/m, where m is the mean molecular mass. To simplify the notation, and 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 dimensions of cm² s⁻² (i.e., the square of a speed) in Gaussian units. We write the first law of thermodynamics,

where U is the internal energy and S is the entropy, as:

(A.1)

where u = U/M and s = S/M are the specific (i.e., per unit mass) internal energy and entropy, respectively, with the gas mass M = Nm = ρV. We also take the caloric equation of state for an ideal gas:

(A.2)

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 modeled the gas as thermally, but not calorically, perfect (i.e., its heat capacities can depend on temperature), and we assumed a piecewise constant or linear dependency:

(A.3)

where we take C_p0 = γ₀/(γ₀ − 1) and C_p1 = γ₁/(γ₁ − 1), with γ₀ and γ₁ specifying the ratios of the heat capacities for T <  T₀ and T >  T₁, respectively. Integrating Eq. (A.2), we can write the specific enthalpy as:

(A.4)

where X = (T − T₀)/(T₁ − T₀). The specific internal energy is defined as u = h − T. Using Eqs. (A.1) and (A.2), we can derive the sound speed c_s:

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:

(A.5)

where s₀ is an integration constant and f(T) is:

(A.6)

Since the specific entropy in the absence of irreversible processes (dissipative heating and 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 used the simpler expression:

(A.7)

In our SDI simulations, we assumed T₀ = 0.01 (GM_⋆/R_⋆), γ₀ = 5/3, T₁ = 0.1 (GM_⋆/R_⋆), and γ₁ = 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:

(B.1)

where ρ_sw, A and B_A are the local plasma density and magnetic field strength, respectively. Approximating the stellar-wind open magnetic flux to be Φ_open ∼ B_AS_A, where S_A is the area of the Alfvén surface, Eq. (B.1) can be manipulated to yield

(B.2)

Using Eqs. (34) and (35) to define the magnetizations Υ and Υ_open, Eq. (B.2) can be rewritten as

(B.3)

We approximate S_A as

(B.4)

where 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, this angle is about 90° for ISW cases, 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 approximate the average cylindrical radius of the Alfvén surface as

(B.5)

we obtain, from Eq. (B.4),

(B.6)

where represents the average cylindrical radius of the Alfvén surface, considering only its geometrical properties. If we assume that (i.e., the effective Alfvén radius), this equation suggests that S_A ∝ ⟨r_A⟩² 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) and Pantolmos & Matt (2017), we assume that υ_sw, A/υ_esc scales as a power law with S_A:

(B.7)

where we suppose that the speed of the plasma at the Alfvén surface is proportional to the specific energy injected in the wind, expressed as a function of the escape speed, times the expansion rate of the flow. The exponent q is determined by the local acceleration of the wind at S_A. Combining Eqs. (B.3), (B.6), and (B.7), we can recover the scaling Eq. (37):

(B.8)

where the exponent m_o = 1/(2 + q) should depend on the speed profile of the wind at the Alfvén surface. If we further assume a power-law scaling for Φ_open/Φ_* with parameter Υ,

(B.9)

we can also recover the scaling Eq. (36):

(B.10)

where m_s = m_o(2n + 1). We can see that while Eq. (B.8) depends on the shape of the wind flux tube and the speed profile at the Alfvén surface, Eq. (B.10) also depends on the fractional open flux dependence on Υ. We notice that, while the scaling Eq. (B.9) represents a reasonable assumption for ISWs, the stellar open flux could also depend on other factors in SDI systems, most notably the mass accretion rate and the Υ_acc accretion parameter.

In order to directly compare the Alfvén radii of stellar winds in isolated and accreting stars without taking into account the exact speed profile Eq. (B.7) or the scaling of the fractional open flux Eq. (B.9), we can combine Eqs. (B.3) and (B.6) to get

(B.11)

All Tables

All Figures

Fig. 1. Logarithmic normalized density (colormaps) showing the temporal evolution of an SDI simulation. The snapshots were taken after five, ten, 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 Rt (dashed black core) and corotation radius Rco (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.

In the text

	Fig. 2. Same as the full-domain panels of Fig. 1 for a simulation of an ISW. The image was taken at three stellar periods.
In the text

	Fig. 3. Normalized truncation radius, Rt/R*, as a function of parameter Υacc. Each symbol (or data point) in the plot represents a single SDI simulation. The solid line shows the fitting function (29), with Kt = 1.1452 and mt = 0.35.
In the text

	Fig. 4. Normalized accretion torque (absolute values), τacc, versus quantity Ṁacc(GM*Rt)1/2. Symbols are the same as in Fig. 3. The solid line corresponds to the best fit of the data, using Eq. (30) with Kacc = 0.79.
In the text

Fig. 5. Normalized effective lever arm, ⟨rA⟩/R*, as a function of the wind magnetization, Υ, for all the cases in this work. As in Fig. 3, each data point is a single simulation. The blue and black colored points (or fitting curves) correspond to SDI and ISW simulations, respectively. Data points with the same symbol, connected with dashed black lines, correspond to numerical solutions that have the same surface magnetic field strength, B*.

In the text

	Fig. 6. ⟨rA⟩/R* versus parameter Υopen for all simulations in this study. Colors, symbols, and line styles have the same meaning as in Fig. 5.
In the text

	Fig. 7. Fractional open flux Φopen/Φ* versus Υ. Colors, symbols, and line styles have the same meaning as in Fig. 5.
In the text

Fig. 8. Left: profiles of normalized Alfvén (dashed), slow-magnetosonic (dash-dotted), and stellar-wind-poloidal (solid) speeds as a function of radius for two cases presented in this work. The velocities are averaged over the θ coordinate and in time. Right: surface area of the stellar-wind flux tube in normalized units versus R/R*. Ssw scales with radial distance as Ssw ∝ Rn, where n ≈ 3 indicates super-radial expansion, and n ≈ 2 and n <  2 indicate radial and sub-radial expansion, respectively. In both panels, the colors are the same as in Fig. 5.

In the text

	Fig. 9. Absolute values of the mass ejection efficiency Ṁsw/Ṁacc and torque efficiency τsw/τacc as a function of time expressed in units of the stellar rotation period. The straight lines in the plot represent the time-averaged values for the two cases exhibiting the lowest and highest values of these two quantities.
In the text

	Fig. 10. Normalized τme versus quantity in curly brackets in Eq. (42). Symbols are the same as in Fig. 3. In the plot, a negative (positive) τme spins up (down) the stellar rotation. Using Eq. (42), the best fit to the points gives Kme = 0.13 and Krot = 0.46.
In the text

	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 MEs.
In the text