© L. Gehrig et al. 2022

1 Introduction

Observations of young pre-main sequence (PMS) stars show that the stellar rotation periods range between ~1−10 days, which is significantly lower than the break up speed (e.g., Irwin & Bouvier 2009; Gallet & Bouvier 2013; Serna et al. 2021). Usually, these stars are surrounded by a protoplanetary disk (e.g., Collier Cameron & Campbell 1993; Hartmann et al. 1998) and are expected to spin up due to the accretion of disk material containing high angular momentum during the PMS phase, which usually lasts ≲10 Myr (e.g., Mamajek 2009). Additionally, the contraction of the young star should reduce the rotational period during the first several mega-year. Throughout the PMS phase, observations of stellar rotation lead to the assumption that the stellar period remains approximately constant during the disk phase (e.g., Edwards et al. 1993; Bouvier et al. 1993; Irwin & Bouvier 2009; Gallet & Bouvier 2013), implying the interaction between the star and the disk results in efficient breaking mechanisms that remove angular momentum from the star (e.g., Serna et al. 2021).

Early studies by Koenigl (1991), based on the model introduced by Ghosh & Lamb (1979), describe angular momentum transfer via the star-disk interaction taking into account strong stellar magnetic fields. Such magnetic fields with a dipole strength of several kG (e.g., Johns-Krull et al. 2009; Johnstone et al. 2014, and references therein) have been observed and the field lines can truncate the disk up to several stellar radii (e.g., Shu et al. 1994; Romanova et al. 2009; Hartmann et al. 2016). Along these field lines the star and the disk can exchange angular momentum and the stellar rotational period can be regulated by the accretion disk (disk-locking model). In more recent studies (Matt & Pudritz 2005b; Zanni & Ferreira 2013), the disk-locking model is shown not to be efficient enough to actually counteract stellar contraction and disk accretion, which in turn would spin up the star. The region, in which the stellar magnetic field lines connect to the disk is too narrow to extract sufficient angular momentum to match the observed stellar rotation periods.

Alternative breaking mechanisms slowing down the star depend on the outflow of material from the star and (or) the disk. First, the so called “X-wind” (Shu et al. 1994) extracts material and angular momentum directly from the disk (close to the co-rotation radius) before it can be accreted onto the star. The amount of angular momentum due to accretion and thus the stellar spin-up is reduced. To launch the “X-wind,” no intrinsic, large scale disk magnetic field is required (e.g., Ferreira et al. 2000) and the field lines are all anchored at the stellar surface. If a large scale magnetic disk field is present, another kind of material outflow can carry away angular momentum from the star. The nature of this outflow depends on the orientation of the stellar magnetic moment compared to the disk’s magnetic moment (Ferreira et al. 2000; Romanova et al. 2009; Zanni & Ferreira 2013). In the case, where the stellar and disk’s magnetic moment are aligned, the so called “ReX-wind” can be launched (Ferreira et al. 2000). At the point inside the disk, the so-called “X-point,” where the stellar poloidal field cancels out the disk field, open field lines can reconnect to the closed stellar magnetic field. Matter parameterized by a specific fraction of the disk’s accretion rate in Ferreira et al. (2000) can be lifted onto and accelerated by stellar rotation along the open field lines at the “X-point” and thus removing angular momentum from the star. Magnetospheric ejections (MEs) can arise in case of a stellar and disk magnetic moment, which are anti-parallel (Romanova et al. 2009; Zanni & Ferreira 2013). A reoccurring process of inflation, opening, and re-connection of field lines leads to an out-bursting behavior, which can transfer angular momentum between the star and the disk. Another model to explain the removal of angular momentum from the star is the accretion powered stellar wind (APSW; Matt & Pudritz 2005a,b, 2008a,b; Matt et al. 2012; Réville et al. 2015; Finley & Matt 2018). The wind is powered by a specific fraction of the disk’s accretion rate and accelerated along the open stellar field lines by stellar rotation removing angular momentum from the star.

An important parameter for the aforementioned theories is the disk’s accretion rate (e.g., Matt & Pudritz 2008b). A precise calculation of the accretion rate including, for example, the influence of stellar magnetic torques and a back-reaction from the influence of stellar rotation on the disk appears rather difficult. Numerical time-step limitations, due to the Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928), in the inner regions of the disk prevent the combined, self-consistent simulation of the stellar rotational evolution and hydrodynamic equations describing the protoplanetary disk. To circumvent this problem, previous studies have parameterized the disk’s accretion rate by a simple power law or exponential decrease with respect to the disk lifetime (e.g., Matt et al. 2010; Gallet et al. 2019). An earlier attempt to combine disk and stellar spin evolution by Armitage & Clarke (1996) has used the diffusion equation to simulate the disk. While the disk’s surface density can be evolved over time, assumptions have to be made, for example, concerning the angular velocity, which has to be made strictly proportional to the Keplerian angular velocity. Thus, a self-consistent calculation, especially of the inner disk parts is not possible as magnetic torques and local pressure gradient can alter the angular disk velocity (e.g., Romanova et al. 2002; Bessolaz et al. 2008; Steiner et al. 2021) and consequently the accretion process itself.

In this study, we want to continue the work of the first paper of this series. Steiner et al. (2021) have shown the importance of stellar magnetic torques for the long-term evolution of protoplanetary disks. We combine hydrodynamic disk and protostellar spin evolution to study the long-term effects of stellar rotation on the protoplanetary disk. The adaptive, implicit hydrodynamic (TAPIR) code for protoplanetary disks (Ragossnig et al. 2020; Steiner et al. 2021), which is based on Dorfi (1998) and Stoekl & Dorfi (2014), is extended with a stellar evolution model. Following Matt et al. (2010), this stellar evolution model includes a simplified description of the protostellar contraction (Kelvin-Helmholtz contraction). Additionally, the stellar spin evolution includes contributions from disk matter accreted onto the star, an accretion powered stellar wind (APSW), as well as magneto-spheric ejections (MEs; e.g., Zanni & Ferreira 2013; Gallet et al. 2019).

By using a hydrodynamic rather than a diffusion approach, a self-consistent calculation of the innermost disk region is possible. This includes, for example, the radial position of the magnetic truncation radius, stellar magnetic torques, the influence of local pressure gradients, as well as episodic accretion outbursts. Those outbursts can be triggered by gravitational fragmentation in the disk outer parts followed by tidal destruction and accretion of inward-migrating clumps (e.g., Vorobyov & Basu 2015), by thermal instabilities (TI, Bell & Lin 1994), magneto-rotational instabilities (MRI; Balbus & Hawley 1991; Armitage et al. 2001) or a joint interaction of gravitational instability and MRI (Martin & Lubow 2011, 2013; Martin et al. 2012; Bae et al. 2014). In this study we focus on the occurrence of MRI triggered outbursts. We refer to the first paper of the series (Steiner et al. 2021) for further details. Additionally, the effect of different stellar rotation periods on the disk’s structure and evolution can be studied.

Due to the implicit nature of the TAPIR code, we can bypass the CFL condition and evolve our model for several mega-year (e.g., Dorfi 1998). Similar to Gallet et al. (2019) and Steiner et al. (2021), we focus on Class II star–disk systems, starting our simulations after the cloud infall has ended and the disk mass corresponds to ≲10% of the stellar mass. Following Gallet et al. (2019), we choose an age of t₀ = 1 Myr as a start time for our simulations.

Finally, we note that “X-winds” and “ReX-winds” are not modeled in this study. A suitable expression for “X-winds” is currently not available. The formulation for “ReX-winds” presented in Ferreira et al. (2000) is very similar to the torque exerted by APSW and we follow Gallet et al. (2019) including APSW and omitting “ReX-winds” in our model.

2 Model description

In this chapter, the model used in our simulations is outlined. The focus of this study is the combination of a full hydrodynamic disk simulation and a stellar spin evolution model. The disk evolution is treated with the implicit hydrodynamics TAPIR code (Ragossnig et al. 2020; Steiner et al. 2021) and the stellar spin model is formulated in detail in Matt et al. (2010, 2012) and Gallet et al. (2019). Here, we want to give a recap of the respective model’s key features. In Sects. 2.1, 2.2 and 2.3 some key features of the TAPIR code, the treatment of the inner disk boundary are shown and the boundary conditions are outlined, respectively. The stellar spin model is depicted in Sect. 2.4.

2.1 Hydrodynamic disk evolution with the TAPIR code

Contrary to long-term global disk simulations based on the diffusion evolution approach for the surface density Σ (e.g., Pringle 1981; Armitage et al. 2001; Zhu et al. 2007, 2010), the TAPIR code allows the treatment of hydrodynamic equations. This includes a deviation from the Keplerian angular velocity u_φ ∝ r^−1/2. Compared to the model presented in Ragossnig et al. (2020), we added the influence of the stellar magnetic field, which is especially important in the inner part of the disk (e.g., Romanova et al. 2002; Bessolaz et al. 2008). The disk is assumed to be axisymmetric in angular direction and to be in hydrostatic equilibrium in vertical direction. Our model is formulated in cylindrical coordinates with the planar components (r, φ), which are then discretized (for a comprehensive description regarding discretization and numerical properties, see Ragossnig et al. 2020). The protoplanetary disk in our simulations is described in Eqs. (1)-(3) as a time-dependent, vertically integrated, viscous accretion disk (e.g., Shakura & Sunyaev 1973; Armitage et al. 2001), $$\frac{\partial }{\partial t}\sum +\nabla \cdot \left(\sum u\right)=0,$$(1) $$\begin{array}{lll}\frac{\partial }{\partial t}\left(\sum u\right)\hfill & +\hfill & \nabla \cdot \left(\sum u:u\right)-\frac{B_{2}B}{2\pi }\hfill \\ \hfill & +\hfill & \nabla P_{\text{gas}}+\nabla \cdot Q+\sum \nabla \psi +H_{\text{p}}\nabla \left(\frac{B_z^2}{4\pi }\right)=0.\hfill \end{array}$$(2) $$\begin{array}{lll}\frac{\partial }{\partial t}\left(\sum e\right)\hfill & +\hfill & \nabla \cdot \left(\sum ue\right)+P_{\text{gas}}\nabla \cdot u\hfill \\ \hfill & +\hfill & Q:\nabla u-4\pi \sum {\kappa }_{\text{R}}\left(J-S\right)+{\dot{E}}_{\text{rad}}=0,\hfill \end{array}$$(3)

where Σ, u, e, P_gas and H_P denote the gas column density, gas velocity in the planar components u = (u_r, u_φ), specific internal energy density, vertically integrated gas pressure and the vertical scale height of the gas disk, respectively. The gradient in planar cylindrical coordinates reads ∇ = (∂/∂r, r−¹∂/∂φ) with ∂/∂φ = 0 for our axisymmetric model. P_gas is utilized by the ideal equation of state P_gas = Σe(γ – 1), with the adiabatic coefficient γ = 5/3 (e.g., Bessolaz et al. 2008) and κ_R denotes the Rosseland-mean opacity (e.g., Mihalas & Mihalas 1984), which is composed of a gaseous component κ_R,gas (based on Caffau et al. 2011) and a dust-dominated component κ_R,dust (based on Pollack et al. 1985); with κ_r = κ_R,gas + f_dustκ_R,dust, where f_dust = 0.01 is the dust-to-gas mass ratio. The gravitational potential of the star-disk system is denoted by ψ. Additionally, the effects of a stellar magnetic field on the disk are included. The vertical component of the magnetic field B_z is assumed to be a dipole and constant within the disk’s vertical extent but can change radially. The planar magnetic field components B = (B_r, B_φ) are taken at the surface of the disk. We note that we neglect the radial component of the magnetic field in this study (B_r = 0, see Rappaport et al. 2004 or Kluzniak & Rappaport 2007). The last term on the left-hand side of Eq. (2) is the vertically integrated magnetic pressure gradient and B_zB denotes the magnetic stress induced by a Lorentz force (as in Vorobyov et al. 2020), which exerts a magnetic torque per unit area on the disk ~rB_zB (e.g., Lubow et al. 1994). The field is assumed to be in rigid rotation with the central star and is able to exchange angular momentum between the star and the disk (Sect. 2.4). Viscous angular momentum transport and viscous heating are represented by the viscous stress tensor Q in their corresponding terms in Eqs. (2) and (3), respectively, which reads as, $$Q=\frac{\mu }{2}\left[\nabla u+{\left(\nabla u\right)}^T-\frac{2}{3}\left(\nabla \cdot u\right)1\right].$$(4)

The factor µ denotes the vertically integrated dynamical viscosity, which yields, $$v=\alpha c_{\text{S}}{H}_{\text{p}},$$(5) $$\mu =v\sum ,$$(6)

where α, c_s, H_p and ν depict the viscosity parameter (Shakura & Sunyaev 1973), the speed of sound, the vertical pressure scale height and the kinematic viscosity, respectively.

The disk is irradiated by the star and cools radiatively, which is described by the net stellar heating and cooling rate per unit surface area Ė_rad. Ė_rad is obtained by balancing stellar irradiation, radiative cooling at the disk surface and irradiation from an ambient medium (see Steiner et al. 2021). The (J – S)-term in Eq. (3) depicts radial radiative transport and is modeled in a radiative diffusion approximation with an Eddington factor of 1/3 (e.g., Ragossnig et al. 2020; Steiner et al. 2021). We note that disk heating due to magnetic coupling to the star is not included in Eq. (3). Work done by the magnetic star-disk interaction is assumed to be comparable to PAV work only close to the inner disk edge as the stellar field strength decreases with r⁻³ (see Sect. 2.2). In this region close to the star, however, the disk temperature is dominated by stellar irradiation. If large scale disk fields are taken into account, disk heating due to magnetic fields can affect the disk’s structure and evolution and is added to Eq. (3) in future studies (Steiner et al. in prep.). Furthermore, resistive dissipation in the disk is assumed to be negligible compared to PAV work (e.g., Miller & Stone 1997). Additionally, we utilized the adaptive grid formulation introduced by Dorfi & Drury (1987) to provide adequate radial grid point resolution and allow the inner radius to move according to the disk evolution (see Sect. 2.2). The grid is solved together with Eqs. (1)-(3) to prevent artificial numerical oscillations (e.g., Dorfi & Drury 1987; Ragossnig et al. 2020) from developing.

2.2 Inner disk boundary

In the innermost disk region, inside the co-rotation radius the stellar magnetic field decelerates the disk material. As the disk’s angular velocity u_φ becomes sub-Keplerian, the disk material accelerates radially toward the star (e.g., Romanova et al. 2002; Steiner et al. 2021). This further leads to a drop in the surface density Σ (e.g., Bessolaz et al. 2008; Steiner et al. 2021) toward the star and strong localized pressure gradients dP_gas/dr. The stellar magnetic torque as well as pressure gradients imply a significant deviation from the commonly used 1D diffusion equation for describing the long-term evolution of Σ (e.g., Bell & Lin 1994; Armitage et al. 2001; Zhu et al. 2009a), therefore it is necessary to fully solve the equations of hydrodynamics. The disk is assumed to eventually be disrupted by a strong stellar magnetic field at the magnetic truncation radius r_trunc Inside r_trunc the gas does not behave as an accretion disk but flows along magnetic funnels toward the star (e.g., Romanova et al. 2002; Romanova & Kurosawa 2014). This magnetically dominated region cannot be described by a 1D model, which justifies the choice of our inner boundary r_in to be at r_in = r_trunc The stellar magnetic field B is approximated as a dipole field being in corotation with the star. The vertical component B_z and angular component B_φ of B are modeled as (e.g., Rappaport et al. 2004; Kluzniak & Rappaport 2007) $$B_z\left(r\right)=B_{*}{\left(\frac{R_{*}}{r}\right)}^3\text{and}$$(7) $$B_{\phi }\left(r\right)\simeq -{\alpha }_{\text{cor}}{B}_z\left(r\right)\left[1-{\left(\frac{{\text{Ω}}_{*}}{\text{Ω}\left(r\right)}\right)}^{{\alpha }_{\text{cor}}}\right],$$(8)

with $${\alpha }_{\text{cor}}\equiv \{\begin{array}{ll}+1\hfill & \text{if\hspace{0.17em}}r<r_{\text{cor}}\hfill \\ -1\hfill & \text{if\hspace{0.17em}}r\ge r_{\text{cor}},\hfill \end{array}$$(9, 10)

as well as the co-rotation radius $$r_{\text{cor}}={\left(\frac{G{M}_{*}{P}_{*}^2}{4{\pi }^2}\right)}^{1/3}.$$(11)

Here, B_*, M_*, Ω_*, P_* and G denote the magnetic field strength at the stellar surface, the stellar mass, the stellar angular velocity, the stellar rotation period and Newton’s gravitational constant, respectively. The factor α_cor switches sign at the corotation radius r_cor and accounts for the braking region inside r_cor (field winds up as a trailing spiral) and the propeling region outside r_cor (field winds up as a leading spiral). Additionally, α_cor ensures that ∣B_φ∣ ≃ B_z for all radii r ≪ r_cor and r ≫ r_cor, which is a physically motivated model, since a tightly wound up field leads to open field lines due to magnetic buoyancy and therefore limits the maximum value of ∣B_φ∣ (e.g., Lynden-Bell & Boily 1994; Rappaport et al. 2004).

The magnetic truncation radius r_trunc denotes the point where the magnetic pressure P_mag equals either the ram pressure P_ram of the in-falling material or the gas pressure P_gas, depending on which pressure contribution dominates (e.g., Koldoba et al. 2002; Romanova et al. 2002; Bessolaz et al. 2008), $$P_{\text{magn}}\left(r_{\text{trunc}}\right)=\max \left[P_{\text{ram}}\left(r_{\text{trunc}}\right),P_{\text{gas}}\left(r_{\text{trunc}}\right)\right].$$(12)

This yields for r_in = r_tmnc the following description for the case P_ram ≥ P_gas (Hartmann et al. 2016), $$\begin{array}{c}{r}_{\text{trunc}}\left(P_{\text{ram}}\ge P_{\text{gas}}\right)\approx 18\xi R_{\odot }{\left(\frac{B_{*}}{{10}^{3}G}\right)}^{4/7}{\left(\frac{R_{*}}{2R_{\odot }}\right)}^{12/7}{\left(\frac{M_{*}}{0.5M_{\odot }}\right)}^{-1/7}\\ {\left(\frac{{\dot{M}}_{*}}{{10}^{-8}{M}_{\odot }\text{yr}}\right)}^{-2/7},\end{array}$$(13)

with the accretion rate onto the star Ṁ_* and a correction factor ξ = 0.7 (following Hartmann et al. 2016) accounting for non-spherical accretion and the oversimplified star-disk interaction. The case P_ram < P_gas follows from equating the magnetic pressure of the vertical field component B_z (Eq. (7)) with the gas pressure, $$r_{\text{trunc}}\left(P_{\text{ram}}<P_{\text{gas}}\right)=\frac{B_{*}{}^{1/3}{R}_{*}}{P_{\text{gas}}{}^{1/6}{\left(8\pi \right)}^{1/6}}.$$(14)

2.3 Numerical boundary conditions

The numerical boundary conditions for the hydrodynamical quantities used in the TAPIR code (cp. Eqs. (1)–(3)) are chosen according to Romanova et al. (2002, 2004) as “free” (Neumann) boundary conditions (e.g., ∂Σ/∂r = 0) at the inner boundary. At the outer boundary, the surface density and the internal energy are also treated with “free” boundary conditions, while the angular velocity u_φ is fixed to the Keplerian value and radial velocity u_r is set to zero representing no inflow of material from an parent cloud or outer envelope, as usual in late Class II systems. The position of the outer boundary is set to the position where Σ = 0.5 g cm⁻² (similar to Vorobyov et al. 2020) and is allowed to move during the simulation according to the disk evolution.

2.4 Stellar spin model

The spin evolution model in this work considers the torques acting on the young star from external contributors as well as the internal stellar evolution. The temporal derivative of angular momentum equals the sum of external torques ${\dot{J}}_{*}={\text{Γ}}_{\text{ext}}$. Following Matt et al. (2010), we assume rigid body rotation with a moment of inertia $I_{\text{*}}=k^2{M}_{\text{*}}{R}_{*}^2$ with the normalized radius of gyration k (k² = 0.2, e.g., Armitage & Clarke 1996; Matt et al. 2010, 2012). The temporal derivative of the stellar angular momentum J_* = 1_*Ω_* reads $${\dot{J}}_{*}=I_{*}{\dot{\text{Ω}}}_{*}+{\dot{I}}_{*}{\dot{\text{Ω}}}_{*}={\text{Γ}}_{\text{ext}}.$$(15)

The evolution of the stellar spin results from rewriting Eq. (15) as follows $${\dot{\text{Ω}}}_{*}=\frac{{\text{Γ}}_{\text{ext}}}{I_{*}}-\frac{{\dot{I}}_{*}}{I_{*}}{\text{Ω}}_{*},$$(16)

with $$\frac{{\dot{I}}_{*}}{I_{*}}\text{=}\frac{{\dot{M}}_{*}}{M_{*}}\text{+2}\frac{{\dot{R}}_{*}}{R_{*}}.$$(17)

Following the model of Gallet et al. (2019), the external torque Γ_ext acting on the star consists of three components; Γ_ext = Γ_acc + Γ_w + Γ_ME. The torques resulting from accreting disk material (gas and dust) onto the star Γ_acc, stellar winds and outflows Γ_w and magnetospheric ejections Γ_ME are described in this section. As we focus on young class II disk systems, we assume a fully convective star. Thus, internal torques resulting from the interaction between a radiative core and a convective envelope (as describe in Gallet et al. 2019) are not considered in this work. The change of the stellar moment of inertia due to contraction of the stellar radius is included in our model.

2.4.1 Stellar radius evolution

We assume a fully convective young star (e.g., Matt et al. 2010; Pantolmos et al. 2020) contracting during the pre-main-sequence until a radiative core develops (Siess et al. 2000). During this phase, the stellar radius evolution can be approximated by a Kelvin-Helmholtz contraction (e.g., Collier Cameron & Campbell 1993; Matt et al. 2010, 2012) resulting in $${\dot{R}}_{*}=2\frac{R_{*}}{M_{*}}{\dot{M}}_{*}-\frac{28\pi \sigma R_{*}^4{T}_{\text{eff}}^4}{3G{M}_{*}^2},$$(18)

with the stellar effective temperature T_eff and the Stefan–Boltzmann constant σ. The first term on the left hand side of Eq. (18) accounts for accreted disk material. The second term represents the contraction on the Kelvin-Helmholtz timescale.

We set the initial stellar mass of the star to M_* = 0.7 M_⊙. A broader range of stellar masses are included into future studies. The respective R_*,init and T_eff,init are taken from the isochrones of Baraffe et al. (2015). The stellar radius then evolves according to Eq. (18) and the effective temperature updates according to the results of Baraffe et al. (2015) with time. We want to note that this stellar evolution model is a simple approximation as several important aspects of the pre-main-sequence evolution are ignored (e.g., Deuterium burning slows down the stellar contraction; Siess et al. 2000; Matt et al. 2010). The maximum difference between Eq. (18) and the results in Baraffe et al. (2015) is <2.5% for a 0.7 M_⊙ star during the first 10 Myr of its lifetime, assuming no accretion.

2.4.2 Magnetospheric ejections

In star-disk systems, with a stellar magnetic field strong enough to connect to the disk material, angular momentum can be transferred between the star and the disk by ME (Zanni & Ferreira 2013), which is nicely summarized in Gallet & Bouvier (2013). Differential rotation between the star and the disk generates a toroidal magnetic field. The resulting magnetic pressure inflates the magnetic field lines, followed by a re-connection and contraction (see Fig. 2 of Zanni & Ferreira 2013). The magnetic torque felt by an annulus (Δr) of the disk (e.g., Armitage & Clarke 1996; Matt & Pudritz 2005a; Gallet et al. 2019; Ireland et al. 2021) is given by $$\text{Δ}{\Gamma }_{\text{ME}}=q{r}^2{B}_{\text{z}}^2\text{Δ}r,$$(19)

with the twist of magnetic field lines $$q\propto \frac{B_{\phi }}{B_z}~\left[1-{\left(\frac{r_{\text{trunc}}}{r_{\text{cor}}}\right)}^{3/2}\right].$$(20)

The additional assumption that the radius of the annulus of the disk connected to the MEs is of the order of (Gallet et al. 2019) yields the total torque felt by the star due to MEs, $${\text{Γ}}_{\text{ME}}=K_{\text{ME}}\frac{B_{*}^2{R}_{*}^6}{r_{\text{trunc}}^3}\left[K_{\text{rot}}-{\left(\frac{r_{\text{trunc}}}{r_{\text{cor}}}\right)}^{3/2}\right],$$(21)

with K_ME = 0.21 and K_rot = 0.7. This parameter choice represents the fact that the MEs rotate at a sub-Keplerian rate (see Gallet et al. 2019). We want to note that Pantolmos et al. (2020) found smaller values for K_ME = 0.13 and K_rot = 0.46. The differences are explained by slightly different boundary definitions for the stellar wind and MEs. Furthermore, Ireland et al. (2021) include multi-dimensional and non-linear effects to the calculation of Γ_ME and also derive lower values for the constant of proportionality compared to Gallet et al. (2019). In this study we use the values given by Gallet et al. (2019), which can be seen as upper limits for Γ_ME.

The effect of the star-disk interaction depends on the position of the inner radius in relation to the corotation radius r_cor. If r_trunc < $K_{\text{rot}}^{\text{2/3}}$ r_cor, the stellar angular momentum increases and the star spins up, otherwise the star spins down. As argued in Zanni & Ferreira (2013), MEs reduce the effect of the accretion torque Γ_acc. In Gallet et al. (2019), a constant K_acc < 1 qualitatively represents this reduction of Γ_acc (see Sect. 2.4.3). K_acc serves as an approximation to quantify small-scale, time-dependent magnetically mediated mass transfer between the star and the disk. Following Gallet et al. (2019), we choose K_acc = 0.4 for all our simulations. As shown in Romanova et al. (2009), the process is very dynamic and thus cannot be adequately described in a 1D model. The following approximations should be seen as time averaged values. We refer the reader to Romanova et al. (2009) and Zanni & Ferreira (2013) for a more elaborated view on the topic.

2.4.3 Accretion

During the accretion of material onto the young star, mass and angular momentum is added to the central object. The accretion torque Γ_acc is given by the product of the specific angular momentum of disk material at the inner edge of the disk and the accretion rate; $${\text{Γ}}_{\text{acc}}=K_{acc}{\dot{M}}_{*}{r}_{\text{trunc}}^2{\text{Ω}}_{\text{disk,in}}$$(22)

with r_trancΩ_disk,in = u_φ,in.

2.4.4 Stellar outflow

In a magnetized stellar outflow, material is transported along magnetic field lines. These field lines are rooted on the star and rotate with Ω_*. As the wind material is moving outwards, it magneto-centrifugally accelerates and thus removes angular momentum from the star.

In this study, we consider the APSW (Matt & Pudritz 2005a), which takes its driving force from the accreting disk material. Consequently, the mass loss due to this wind is a fraction W of the disk accretion rate; Ṁ_W = WṀ_disk. The torque, caused by a stellar outflow (Weber & Davis 1967), acting on the star, can be written as $${\text{Γ}}_{\text{w}}={\dot{M}}_{\text{w}}{r}_{\text{A}}^2{\text{Ω}}_{*},$$(23)

with the Alfvén radius $$r_{\text{A}}=K_1{\left[\frac{B_{*}^2{R}_{*}^2}{{\dot{M}}_{\text{wind}}\sqrt{K_2^2{v}_{\text{esc}}^2+{\text{Ω}}_{*}^2{R}_{*}^2}}\right]}^m{R}_{*},$$(24)

where $v_{\text{esc}}=\sqrt{2\text{G}{M}_{*}/R_{*}}$, m = 0.2177, K₂ = 0.0506 and K₁ = 1.7 arising as best fit values from magneto-hydrodynamic stellar wind simulations (see Matt et al. 2012; Gallet et al. 2019). We want to note that Γ_w acts as an angular momentum “sink” in our model. The exact driving mechanisms of the APSW and the order of W are still under debate (we refer to the discussion in Gallet et al. 2019 for further details).

2.4.5 Numerical integration

To include the spin evolution model in the TAPIR code, we updated the stellar parameters after each timestep At as follows. The current accretion rate resulting from solving the disk equation Eqs. (1)–(3) updates the stellar mass and radius. With the updated values the respective torque contributions acting during this timestep are calculated according to the formulation above. Finally, $\dot{\text{Ω}}$ and the new stellar angular velocity ${\text{Ω}}_{\text{new}}{\text{=Ω}}_{\text{old}}+\dot{\text{Ω}}\text{Δ}t$ is derived and used in the next timestep. Ω_old and Ω_new denote the stellar angular velocity used for the current and the new timestep, respectively.

3 Results

In this section, we want to define the parameter space and a reference model (Sect. 3.1), followed by the results of our simulations: the back-reaction of different stellar rotation rates on the disk (Sect. 3.4), the construction of a fast and a slow rotating protostar within our parameter space (Sect. 3.3), as well as the influence of stellar rotation on episodic accretion outburst events (Sect. 3.5).

3.1 Parameter space

We define a reasonable parameter space for stellar and disk properties covered in our simulations, which is based on observations and previous theoretical models. First we define the stellar parameters appropriate for a 1 Myr old T Tauri star. As described in Sect. 2, the initial stellar mass is M_* = 0.7 M_⊙ and the initial stellar radius, as well as the effective temperature are taken from isochrones provided by Baraffe et al. (2015). The initial stellar rotation period for a 1 Myr old star can vary significantly. We choose a range for the initial period from 2 to 15 days. This choice is motivated by observations of the Orion Nebular Cluster (ONC; Herbst et al. 2002). The stellar dipole magnetic field strength has to be estimated. Measurements by Johnstone et al. (2014) suggest that the range of the dipole component of the magnetic field can range from several hundred Gauss to ~kG. Furthermore, M dwarfs can have very large magnetic field strengths of several thousand Gauss (e.g., Kochukhov & Reiners 2020). Based on these measurements we assume B_* to be in a range between 500 and 2000 G The magnetic field strength remains constant during our simulations. We note that T Tauri stars are usually magnetically saturated (Rossby number ≲0.1, e.g., Briggs et al. 2007; Johnstone et al. 2021) and the stellar magnetic field is roughly independent of the rotational period (e.g., Reiners et al. 2009; Lavail et al. 2019). During later evolutionary stages, this relation is reverted and fast rotating stars show more magnetic activity up to the state of magnetic saturation (e.g., Wright et al. 2011). Finally, the efficiency W of the APSW has to be chosen. We choose three different values for W. Based on the theoretical model by Cranmer (2008) and simulations by Pantolmos et al. (2020), we choose W = 2% as a reference value. W = 0% and W = 5% represent a weak and strong APSW, respectively.

The hydrodynamic disk model requires an initial disk accretion rate as well as description of the viscous α parameter. Based on observations (e.g., Gullbring et al. 1998; Manara et al. 2016) and reproduced by numerical simulations (e.g., Vorobyov et al. 2017b), the accretion rate from the disk onto the star ranges between 10⁻¹⁰ and 10⁻⁷ M_⊙ yr⁻¹ for late-type star-disk systems. We therefore choose the initial accretion rate Ṁ_int = 1 × 10⁻⁸ M_⊙ yr⁻¹. The effects of a variation of Ṁ_int are covered in a further study. The viscous α-parameter ranges between ≲0.1 and ~0.0001 (e.g., Zhu et al. 2007, 2009b; Vorobyov & Basu 2009; Yang et al. 2018). We adopt initially a constant value of α = 0.005 used in Sects. 3.4 and 3.3, unless otherwise stated.

To generate episodic accretion outbursts (Sect. 3.5) due to MRI (Balbus & Hawley 1991), we adapt the layered disk model (Gammie 1996). A detailed description of this model in the TAPIR code is presented in Steiner et al. (2021). In the layered disk model, the disk is divided in two layers. On the one hand, we assume an MRI active layer, which is sufficiently ionized due to stellar irradiation, cosmic rays or collisional ionization of, for example, potassium (Gammie 1996) with an effective viscosity α_MRI. On the other hand, we have an MRI inactive layer, the dead zone (DZ), with α_Dz ≪_MRI, in which the mass transport rate is strongly diminished. This leads to mass piling up in the DZ up to a point where the increasing mass leads to increasing temperature and eventually to onset of the MRI above a critical temperature T_crit. The turbulent viscosity and viscous heating increase sharply, which further causes ionization of the whole inner region and a highly increased mass transport rate Ṁ. This ionized inner region is then accreted in an strong accretion event onto the star and the inner disk is replenished by material from the outer disk, which completes a burst cycle.

Depending on the choice of the previously mentioned parameters, the evolution of the accretion disk and therefore the spin evolution of the star can vary. One example is the viscosity parameter a that influences the viscous timescale of the accretion process; τ_v = r²/v ~ 1/α. A large (small) α-value corresponds to short (long) disk lifetime, respectively, and the accretion rate changes accordingly (e.g., Vorobyov & Basu 2009). In our model, we can include these effects self-consistently. In Fig. 1, the effect of different α-values on the disk evolution is shown. The red, black and blue lines show the temporal evolution of the disk’s accretion rate Ṁ_disk for a constant α-value of 0.01, 0.005 and 0.0025, respectively. The time for Ṁ_disk to drop differs by mega-years, which has an effect on the stellar spin evolution. Furthermore, our model enables the consideration of episodic accretion outbursts (see Steiner et al. 2021), based on the layered disk model (Gammie 1996) with α_MRI = 0.02 and α_DZ = 0.0004 (similar to Steiner et al. 2021). During such an outburst, the accretion rate can rise for orders of magnitude (yellow line in Fig. 1) and the disk’s structure and temperature changes significantly. These effects are also calculated in a self-consistent manner. Previous stellar spin evolution studies tend to use a simplified description of the disk’s accretion rate. Either an exponential decay (e.g., Matt et al. 2010) or a power law description (Gallet et al. 2019) is used, which require free parameters to define for example the disk lifetime. Additionally, accretion outbursts as well as the back-reactions of such events on the disk cannot be included into these models.

Fig. 1 Disk accretion rate Ṁdisk with respect to the simulation time t in years starting from t0 = 1 Myr. The red, black, blue lines correspond to a constant α value of 0.01, 0.005 and 0.0025, respectively. The yellow line corresponds to a layered disk model (Gammie 1996) with αMRI = 0.02 and αDZ = 0.0004 (similar to Steiner et al. 2021).

3.2 Reference model

Based on the stellar and disk parameter ranges, defined in Sect. 3.1, we define a reference model for both a constant α-value and a layered disk model.

3.2.1 Constant α-value

The stellar and disk parameters for the reference model using a constant α-value (Mref) are summarized in Table 1. The stellar spin evolution for this model is shown in Fig. 2. During the simulation time of 5 Myr, the stellar rotation period decreases slightly from the initial value of 6 toward 4 days (see panel a). The individual external contributions to the total torque Γ_tot are shown in panel b (black line). A positive (negative) torque value adds (removes) angular momentum to the star and the star spins up (spins down), respectively. The accretion torque Γ_acc (blue line) spins up the star during the whole simulation time, whereas the APSW (Γ_w, yellow line) spins down the star. If MEs (Γ_ME, red line) spin the star up or down, depends on the position of the inner disk radius r_trunc with respect to the co-rotation radius r_cor (see Eq. (21)), which is shown in panel c. For comparison, the value $K_{\text{rot}}^{\text{2/3}}$ is indicated by the dashed black line.

During the first 2 × 10⁵ yr, r_trunc is close to r_cor and r_trunc/r_cor > $K_{\text{rot}}^{\text{2/3}}$. Γ_ME is negative and the MEs remove angular momentum from the star. The influence of the APSW and the ME, however, are not sufficient to counteract accretion and contraction and the star spins up slightly (Γ_tot > 0 after 10⁴ yr). While the star contracts and the stellar radius as well as the inner disk radius decreases, the co-rotation radius resides at approximately the same positions as the stellar rotational period only changes slightly during the first couple 10⁵ yr. Thus, r_trunc/r_cor decreases until r_trunc/r_com < $K_{\text{rot}}^{\text{2/3}}$ and the MEs also spin up the star. After 2 Myr of accretion, the disk is nearly depleted and the disk accretion rates drops (see Fig. 1). Consequently, r_trunc moves outwards and the MEs remove angular momentum from the star again.

Fig. 2 Stellar spin evolution for Mref with respect to the simulation time t in years starting from t0 = 1 Myr. Panel a shows the stellar rotation period P* in days. The dashed black line indicates the initial period of 6 days. Panel b shows the external torque contributions (see Sect. 2.4) as well as the total torque acting on the star (black line). The blue, red and yellow line represent the accretion torque Γacc, torque due to MEs ΓME and the APSW torque Γw, respectively. The dashed line at Γ = 0.0 indicate the boundary between a spin up and a spin down contribution. Panel c shows the position of the inner disk radius rtrunc with respect to the position of the co-rotation radius rcor The dashed line indicate $K_{\text{rot}}^{2/3}$.

3.2.2 Layered disk model

The reference parameters (MrefB) used in combination with the layered disk model differ from Mref only by the viscous α-value. Instead of a constant α-value trough out the disk, we use α_MRI = 0.02 and α_DZ = 0.0004 (similar to Steiner et al. 2021). The effect of different stellar spin evolutionary tracks based on this reference model are shown in Sect. 3.5. In Fig. 3, a detailed view on accretion bursts and the respective contributions to the total torque acting on the star is shown. In panel a, the disk’s accretion rate is shown. Once the disk temperature reaches T_crit = 1000 K and the MRI is triggered, the disk’s accretion rate Ṁ_diSk increases sharply over an order of magnitude until the material piled up in the DZ is depleted and Ṁ_disk decreases again until the next outburst is triggered. In panel b, the position of the inner radius r_trunc with respect to the co-rotation radius r_cor is shown. During the time of increased Ṁ_disk, r_trunc is pushed toward the star (see Eq. (13)). This again effects the total torque acting on the star (black line in panel c). During the first several 10¹ yr, the total as well as the external torque contributions are similar to the model Mref and the star spins down, which is indicated by Γ_tot < 0. When the accretion outburst starts, Γ_acc (Γ_w) increases (decreases) due to the higher accretion rate, respectively, without changing their sign. The contribution of the MEs to Γ_tot, however, changes significantly. Before the outburst, r_trunc is close to r_cor and Γ_ME < 0 resulting in a stellar spin down effect. During the burst the distance between r_trunc and r_cor increases drastically and r_trunc/r_cor < $K_{\text{rot}}^{\text{2/3}}$. Now, the MEs spin up the star and the APSW alone cannot counteract the combined spin up torque consisting of Γ_acc, Γ_ME and the stellar contraction. As a result, during the accretion outburst, spins up the star. While the peak accretion rate of the outburst in Fig. 3 is similar to EXor outbursts, the duration is reminiscent of FUor outbursts (e.g., Audard et al. 2014). We note that the peak and duration of these outbursts can be influenced by the choice of α_MRI, α_DZ and T_crit (e.g., Steiner et al. 2021) and the effects of different outburst types are investigated in appropriate further studies.

Fig. 3 Stellar spin evolution for MrefB with respect to the simulation time t in years starting from t0 = 1 Myr. Panel a shows the disk’s accretion rate Ṁdisk. Panel b shows the position of the inner disk radius rtrunc with respect to the positions of the co-rotation radius rcor. The dashed line indicate $K_{\text{rot}}^{\text{2/3}}$ Panel c show the total torque acting on the star Γtot as well as the external torque contributions, which are the same as in Fig. 2. Additionally, the dashed line at Γ = 0.0 indicate the boundary between a spin up and a spin down contribution.

Fig. 4 Evolution of the stellar rotational period P* with respect to the simulation time t in years starting from t0 = 1 Myr for the models Mref (solid line), Mfast (dashed line) and Mslow (dash-dotted line). Model MrefW5 (solid blue line) highlights the dependence on B* and the APSW value W (see text). The area between Mfast and Mslow is colored for better visualization. The observational data from four young stellar clusters are shown as “crosses” with the 50th percentile (big circle) and the 95th (5th) percentile (downward and upward pointing triangles). From left to right the open clusters are: ONC (Herbst et al. 2002), NGC 6530 (Henderson & Stassun 2012), NGC 2264 (Affer et al. 2013) and NGC 2362 (Irwin et al. 2008).

3.3 Fast and slow rotating stars

Within the predefined parameter space (Sect. 3.1), we define models for a fast and slow rotating star to test the limits of our approach and compare the results to observations of stellar rotation periods of young clusters. Based on our reference model Mref, we change the stellar magnetic field strength B_*, the APSW value W and the initial stellar period P_init. The fast (slow) rotation model Mfast (Mslow) has a stellar magnetic field strength of B_* = 500 (2000) G, W = 0 (5)% and P_init =2 (15) days, respectively¹.

In Fig. 4, the rotational periods of Mref (solid line), Mfast (dashed line) and Mslow (dash-dotted line) are shown. The stellar rotation periods of four young clusters are also shown to compare the results to the observations. The Orion Nebular Cluster (ONC) with an age of 1 Myr (Herbst et al. 2002) was already used to define the initial parameters used in this study and is located at t = 0 yr. The following clusters from left to right are NGC 6530 (Henderson & Stassun 2012) with an age of 1.65 Myr, NGC 2264 (Affer et al. 2013) with an age of 2 Myr and NGC 2362 (Irwin et al. 2008) with an age of 5 Myr. The solid black circles and the downward (upward) pointing triangles represent the 50th and 95th (5th) percentiles, respectively. The reference model Mref remains close to the 50th percentile value (mean value) of all four clusters. Additionally, the majority of stellar rotation periods in NGC 6530, NGC 2264, as well as NGC 2362 are enclosed between Mfast and Mslow, which is visualized by the blue area between the respective curves. The spin evolution of Mfast and Mslow depends on the choice of B_* and W². A weaker (stronger) magnetic field strength result in a shift toward shorter (longer) periods. In a similar way, a smaller (larger) APSW value W result in shorter (longer) periods. The order of magnitude of this dependence is symbolized by model MrefW5 (solid blue line in Fig. 4). Reducing B_* by a factor of 2 spins up the star by 92% after 5 Myr (compare models Mslow and MrefW5). An increase of W from 2% to 5% spins down the star by 60% after 5 Myr (compare models Mref and MrefW5).

3.4 Back-reaction on the disk

In this section, we show the effects stellar rotation has on the protoplanetary accretion disk. First, the effects of different spin evolutionary tracks are analyzed (Sect. 3.4.1), by comparing a star that spins up with a star that spins down. Second, we show the effects that different initial rotation periods have on the accretion disk (Sect. 3.4.2).

3.4.1 The effect of different stellar spin evolutionary tracks

To study the effects a star that spins up has on the accretion disk, compared with a star that spins down, we change the APSW efficiency value W. While the initial stellar and disk properties are identical the APSW removes more (less) angular momentum from the star that spins down (up), respectively. Thus, the effect of different stellar spin evolutionary tracks on the accretion disk, starting from the same initial model, can be shown.

We choose, based on the reference value of W = 2%, W = 0%, as well as W = 5% representing a star that spins up and down, respectively. In Fig. 5, we show the resulting evolution of the stellar rotational period (Panel a). The model Mref (red line) evolves as already shown in Fig. 2. As expected, a low APSW value W = 0% (black line) does not remove enough angular momentum from the star to counteract accretion and contraction and the star spins up during the whole simulation time. A large APSW value W = 5% (blue line) removes considerably more stellar angular momentum and the star spins down toward 8 days. The resulting effects of these different stellar spin evolutionary tracks on the accretion disk are shown in panels b and c. In panel b, the surface density distributions of the innermost disk part (<0.5 AU) are shown for different time snapshots: the initial model at t = 0.0 Myr (solid line), t = 0.1 Myr (dashed lines) and t = 1.0 Myr (dotted lines). Additionally the co-rotation radii are shown for the respective time snapshots (vertical lines). Starting from the same initial model (the blue and black solid line are identical), the surface density structure of the inner disk part differs clearly when changing the APSW value from 0 to 5%. Responsible for this difference is the position of the co-rotation radius r_cor. For a star that spins down (blue lines), the co-rotation radius moves outwards away from the star and a larger part of the disk in located inside r_cor. Disk material inside r_cor rotates faster compared to the star and is slowed down due to magnetic torques. Thus, the disk’s angular velocity (µ_φ) decreases below the Keplerian velocity and the disk material inside r_cor accelerates in radial direction toward the star (similar to Steiner et al. 2021). For a star that spins up, r_cor moves inward toward the star. The disk material that is now located outside r_cor rotates slower compared to the star and is accelerated by magnetic torques to velocities faster than Keplerian. As a result, the radial inward velocity of the disk outside r_cor is decreased. Material transported from the outer parts toward the inner disk region, is piled up outside r_cor, pushing the inner disk radius r_trunc closer to the star.

With the position of r_cor also the position of the inner disk radius r_trunc changes. With r_trunc moving closer to or away from the star, the inner disk is heated differently due to stellar irradiation. This effect is shown in panel c. The maximum disk midplane temperature T_{disk, max} shows a different behaviour depending on the stellar spin evolutionary track. As expected, T_disk,max decreases for models that show a stellar spin down (blue and red lines). For the model with W = 0% showing a stellar spin up, T_{disk, ax} increases for 1.0 Myr and differs by over 300 K compared to the model with W = 5%. Such a difference can have a considerable effect on the disk’s ionization fraction that is strongly temperature dependent in these regions of T_disk,max (e.g., Fromang et al. 2002; Vorobyov et al. 2020) as well as on the occurrence of episodic accretion outbursts.

Fig. 5 Panel a shows the evolution of the stellar rotational period P* with respect to the simulation time t in years starting from t0 = 1 Myr. The reference model Mref with W = 2%, the models with W = 0%, and W = 5% are given by the red, black and blue lines, respectively. Panel b shows the surface density of the innermost disk part (0.5 AU) for the models with W = 0% (black lines) and W = 5% (blue lines) at three different time snapshots: The initial model at t = 0.0 Myr (solid line), t = 0.1 Myr (dashed lines) and t = 1.0 Myr (dotted lines). The vertical lines represent the position of the co-rotation radius for the given APSW efficiency value and time snapshot. Note that the black and blue solid line are identical as the same initial model (t = 0.0 Myr) is used for this comparison. Panel c shows the evolution of the maximum disk temperature Tdisk,max at the disk’s midplane with respect to the simulation time t in years starting from t0 = 1 Myr. The line colors represent the same models as in panel a.

Fig. 6 Evolution of the stellar rotational period P* (panel a) and the maximum disk midplane temperature Γdisk,max (panel b) with respect to the simulation time t in years starting from t0 = 1 Myr. Different colors represent different initial stellar rotational periods: Pinit = 2 days (red line), Pinit = 4 days (magenta line), Pinit = 6 days (black line), Pinit = 10 days (blue line) and Pinit = 15 days (yellow line).

3.4.2 The effect of different initial stellar periods

When changing the initial stellar rotation period P_init we expect the period to converge toward a common value in the course of ~Myr as seen in Armitage & Clarke (1996) or Matt et al. (2010). In Fig. 6, we compare the stellar spin evolution for initial stellar periods ranging from 2 to 15 days. Panel a shows the rotation period for P_init = 2 days (red line), P_init = 4 days (magenta line), P_init = 6 days (black line), P_init = 10 days (blue line) and P_init = 15 days (yellow line). While the stellar periods approach each other, we still see a spread of 2 days after 2 Myr. The reason for this difference is the back-reaction of different stellar spin evolutionary tracks on the accretion disk, which is explained in Sect. 3.4.1.

The maximum disk midplane temperature T_disk,max for each model is shown in panel b. Similar to Fig. 5, a fast rotating star with an inner disk and co-rotation radius close to the star has a higher disk temperature. T_disk,max varies by a factor of 2 for the models shown. Again, such a difference strongly affects the disk’s ionization as well as influence the occurrence of accretion outbursts.

3.5 Effect on accretion outbursts

Finally, we want to show the effects of different stellar spin evolutionary tracks on the occurrence of accretion outbursts using the layered disk reference model MrefB. In Fig. 7, similar to Sect. 3.4.1, different stellar spin evolutionary tracks are shown. Starting from a common initial rotational period of 6 days, the stellar rotational period (panel a) of the reference model MrefB (red line) does not change significantly over a simulation time of 1 Myr, whereas an APSW value of W = 0% results in a stellar spin up (black line) and W = 5% results in a slow rotating star (blue line). In panel b, the respective disk accretion rates are shown. Similar to the models in Steiner et al. (2021), outbursts occur over a time span of ≲1 Myr. The number of bursts as well as the time span on which outbursts are triggered depend on the stellar spin evolution. Fast rotating stars tend to have more mass located closer to the star over longer periods of time (see Fig. 5). Thus, the disk temperature remains higher within these disks and the conditions that trigger an outburst are also met over longer periods of time. This results in a longer period as well as a larger number of outbursts for stars that spin up and rotate faster. In the model MrefB (red line in Fig. 7) 231 outbursts occur over a time span of 0.96 Myr. In the model with W = 5%, representing a star that spins down, 85 outbursts occur over a time span of 0.53 Myr, whereas in the model with W = 0%, representing a fast rotating star, 835 outbursts occur over a time span of 1.11 Myr. Each outburst influences the temperature structure (e.g., Steiner et al. 2021) and leads to a change in the chemical composition of the disk (e.g., Rab et al. 2017, and references therein). This can change observables such as the snow line (e.g., Cieza et al. 2016), the abundance of HCO⁺ (Jørgensen et al. 2013) or CO observations (Rab et al. 2017).

Fig. 7 Stellar spin evolution for different stellar spin evolutionary tracks of the layered disk reference model with respect to the simulation time t in years starting from t0 = 1 Myr. Panel a shows the stellar rotational period P* for the reference model MrefB with W = 2% (red line) as well as for models with W = 0% (black line) and W = 5% (blue line). Panel b shows the disk accretion rate Ṁdisk for the respective models.

4 Discussion

To explain the rotational behaviour of young stars, still surrounded by a protoplanetary disk, the star-disk interaction and its associated torques have to be considered (see Sect. 2.4). The disk component in previous studies, however, has been modeled in a simplistic way and important disk features and parameters have not been included in the stellar spin evolution (e.g., Matt et al. 2010, 2012; Gallet et al. 2019). To improve the treatment of the disk, we combine hydrodynamic disk and stellar spin evolution to simulate the disk’s influence on the stellar rotation and the stellar back-reaction on the disk. The influence of stellar magnetic torques on the disk, the position of the inner disk radius, disk radial density, temperature profiles and the effects on episodic accretion outbursts can be calculated self-consistently during the whole simulation time.

A self-consistent treatment of both stellar and disk evolution introduces further complexity. As the stellar spin evolution depends on a precise formulation of the processes of the inner disk, this region has to be well resolved during the whole simulation ranging over several mega-year. Additionally, magnetospheric accretion is an intrinsically multi-dimensional problem (e.g., Romanova et al. 2002; Bessolaz et al. 2008; Romanova & Kurosawa 2014). In addition to timestep limitations imposed by the CFL condition (Courant et al. 1928), current explicit, multi-dimensional simulations either exclude the innermost disk regions or treat them in a simplified manner (e.g., Vorobyov et al. 2019) or converge toward a steady-state solution within several rotational periods of the inner disk (e.g., Romanova et al. 2004; Zhu et al. 2020; Pantolmos et al. 2020). These scenarios are not ideal to study the dynamic stellar spin evolution over several mega-year. The TAPIR code is able to spatially resolve the inner disk region sufficiently well over the whole disk lifetime. Solving the hydrodynamic equations (cp. Eqs. (1)–(3)), the accretion disk can be modeled in a more self-consistent manner compared to previous stellar spin evolution studies (e.g., Matt et al. 2010, 2012; Gallet et al. 2019). We have to note that the one dimensional TAPIR code requires a parameterization of the multi-dimensional processes occurring in the inner disk. Consequently, the results presented above are based on simplifications. Nevertheless, the general trends and dependencies shown in this study should be valid.

Young stars are assumed to have strong magnetic fields (e.g., Johns-Krull et al. 2009; Gregory et al. 2012; Johnstone et al. 2014). The complex field topology and temporal variability of the field, however, is difficult to observe or to predict and equally challenging to develop in a numerical model, which describes the field dynamics correctly and simultaneously maintains the ability to carry out long-term simulations. Based on the measurement technique, different values for the stellar magnetic field are observed for the same star. This issue can be explained on the basis of the young M dwarf AU Mic (Kochukhov & Reiners 2020). A mean magnetic field of B_mean = 2.3 kG results from Stokes I profile measurements. With a more detailed Zeeman Doppler imaging technique based on the Stokes V parameter, a global field strength of B_global = 0.09 kG is observed. Finally, based on the Stokes Q and U parameters a dipole field strength of B_dipole = 2.0 kG is obtained. Thus, the observed magnetic field can change based on the used measurement technique by a factor of 20. Finley & Matt (2018), however, have shown in stellar wind calculations that usually the lowest order field component dominates Γ_w. Additionally, stellar magnetic fields can vary on short timescales. Johnstone et al. (2014) have shown that the dipole magnetic field strength of V2129 Oph varies by a factor of almost 4 over a four year period.

Stellar winds, in this study represented by APSW, strongly affect the long-term stellar spin evolution. While an APSW influences not only the stellar spin evolution providing a possible solution to the occurrence of slow spinning stars, it also affects the disk structure as well as its long term evolution (see Figs. 5 and 7). The exact nature and strength of these winds, however, are still debated (see the discussion in Gallet et al. 2019). Previous authors need large APSW values of W ≳ 10% to explain slowly rotating stars (e.g., Matt et al. 2012; Gallet et al. 2019; Ireland et al. 2021), which contradicts the results of Cranmer (2008) and Pantolmos et al. (2020), limiting W ~ 1%. We can reduce this value to W ≲ 5% to reproduce the observed rotational period of most stars in young stellar clusters (see Fig. 4). Apart from APSW, other wind mechanisms can be important for the stellar and disk evolution. Especially magnetic disk winds affect the disk evolution noticeably (e.g., Königl et al. 2011). By removing matter and, in the presence of a large-scale magnetic field threading the disk, angular momentum from the disk before it can be accreted onto the star, the stellar spin evolution is affected. Similar, photo-evaporation (internal or external) can influence the disks lifetime and thus, the stellar spin evolution (e.g., Roquette et al. 2021). These effects are added to our model in further studies.

The stellar spin evolution is strongly influenced by the initial stellar parameters. Besides the magnetic field strength, the initial radius as well as the temperature play an important role. The change of the stellar moment of inertia (İ_*) can vary significantly based on the choice of initial parameters (Eq. (18)).

We take initial stellar values from Baraffe et al. (2015) and evolve the star using a rather simplified model to describe Ṙ_*. As described in Matt et al. (2010) such assumptions reproduce more sophisticated stellar evolution models to a reasonable degree of accuracy. Siess et al. (2000), however, concludes that the variation of, for example, stellar metallicity can vary T_* by ≳200 K. Additionally, a spin up or spin down torque can alter the rotational period and thus stellar evolution (Zeipel 1924; Endal & Sofia 1976; Maeder & Meynet 2000; Tassoul 2000; Paxton et al. 2019; Amard et al. 2019). Moreover, the amount of accreting energy that is radiated away from the star influences the stellar evolution. The deviation in stellar radius and luminosity between models, in which most of the accreted energy is radiated away (cold accretion) and a certain amount of accreted energy is added to the star (hybrid or hot accretion), can reach up to 50% (e.g., Baraffe et al. 2009; Vorobyov et al. 2017b). Especially during an accretion outburst, the stellar radius and effective temperature is influenced (e.g., Vorobyov et al. 2017b), which is not yet included into our model. Thus, the combination of our disk hydrodynamics model with a more elaborate stellar evolution code (e.g., the MESA code; Paxton et al. 2019) could include these effects in a self-consistent way (similar to e.g., Vorobyov et al. 2017a) and would be a further improvement of the stellar spin evolution model.

5 Summary and conclusion

In this study, we have enhanced the hydrodynamic disk evolution TAPIR code (Ragossnig et al. 2020; Steiner et al. 2021) with a stellar spin evolution model (see Sect. 2.4). Several properties of the disk (e.g., the disk’s accretion rate or the disk lifetime) that have been approximated in previous studies (e.g., Matt et al. 2010; Gallet et al. 2019) are now calculated in a self-consistent way. Additionally, the back reaction of stellar rotation on the disk can be included. Different stellar spin evolutionary tracks and the associated movement of the co-rotation radius influences the inner disk structure. The inner disk radius of a star that spins up is located closer to the star (see Fig. 5). More disk material is heated up due to the proximity to the star and the disk temperature increases. Similarly, the disk around a star with large initial rotational period remains relatively cold, compared to a initially fast rotating star (see Fig. 6). From our range of parameters presented in Sect. 3.1, we can produce fast and slow rotating models (see Fig. 4), which are able to enclose most of the stellar rotational periods observed in young stellar clusters (e.g., Henderson & Stassun 2012; Gallet & Bouvier 2013) with smaller APSW values compared to previous studies (e.g., Matt et al. 2012; Gallet et al. 2019; Ireland et al. 2021). Furthermore, for the first time, the effects of episodic outburst events on the stellar rotation and vice versa can be studied over several mega-year. Disks around stars that spin up show a larger number of outbursts as well as a longer outbursting period, compared to slowly rotating stars. While we have shown the effects of different stellar spin evolutionary tracks on the occurrence of outbursts, further studies have to be conducted to analyze a wider range of parameters such as, for example, different stellar and disk masses, outburst strengths as well as including the early evolution of the star-disk system (setting t₀ to ≲ 10⁴ yr). Our model that combines the stellar spin and hydrodynamic disk evolution enables us to include these effects in further studies and work toward a better understanding of the stellar spin evolution.

Acknowledgements

The authors are thankful to the anonymous referee for useful suggestions that helped to improve the manuscript. E.I.V. acknowledges support of Ministry of Science and Higher Education of the Russian Federation under the grant 075-15-2020-780 (N13.1902.21.0039; Sects. 3 and 4).

References

All Tables

All Figures

	Fig. 1 Disk accretion rate Ṁdisk with respect to the simulation time t in years starting from t0 = 1 Myr. The red, black, blue lines correspond to a constant α value of 0.01, 0.005 and 0.0025, respectively. The yellow line corresponds to a layered disk model (Gammie 1996) with αMRI = 0.02 and αDZ = 0.0004 (similar to Steiner et al. 2021).
In the text

Fig. 2 Stellar spin evolution for Mref with respect to the simulation time t in years starting from t0 = 1 Myr. Panel a shows the stellar rotation period P* in days. The dashed black line indicates the initial period of 6 days. Panel b shows the external torque contributions (see Sect. 2.4) as well as the total torque acting on the star (black line). The blue, red and yellow line represent the accretion torque Γacc, torque due to MEs ΓME and the APSW torque Γw, respectively. The dashed line at Γ = 0.0 indicate the boundary between a spin up and a spin down contribution. Panel c shows the position of the inner disk radius rtrunc with respect to the position of the co-rotation radius rcor The dashed line indicate $K_{\text{rot}}^{2/3}$.

In the text

Fig. 3 Stellar spin evolution for MrefB with respect to the simulation time t in years starting from t0 = 1 Myr. Panel a shows the disk’s accretion rate Ṁdisk. Panel b shows the position of the inner disk radius rtrunc with respect to the positions of the co-rotation radius rcor. The dashed line indicate $K_{\text{rot}}^{\text{2/3}}$ Panel c show the total torque acting on the star Γtot as well as the external torque contributions, which are the same as in Fig. 2. Additionally, the dashed line at Γ = 0.0 indicate the boundary between a spin up and a spin down contribution.

In the text

Fig. 4 Evolution of the stellar rotational period P* with respect to the simulation time t in years starting from t0 = 1 Myr for the models Mref (solid line), Mfast (dashed line) and Mslow (dash-dotted line). Model MrefW5 (solid blue line) highlights the dependence on B* and the APSW value W (see text). The area between Mfast and Mslow is colored for better visualization. The observational data from four young stellar clusters are shown as “crosses” with the 50th percentile (big circle) and the 95th (5th) percentile (downward and upward pointing triangles). From left to right the open clusters are: ONC (Herbst et al. 2002), NGC 6530 (Henderson & Stassun 2012), NGC 2264 (Affer et al. 2013) and NGC 2362 (Irwin et al. 2008).

In the text

Fig. 5 Panel a shows the evolution of the stellar rotational period P* with respect to the simulation time t in years starting from t0 = 1 Myr. The reference model Mref with W = 2%, the models with W = 0%, and W = 5% are given by the red, black and blue lines, respectively. Panel b shows the surface density of the innermost disk part (0.5 AU) for the models with W = 0% (black lines) and W = 5% (blue lines) at three different time snapshots: The initial model at t = 0.0 Myr (solid line), t = 0.1 Myr (dashed lines) and t = 1.0 Myr (dotted lines). The vertical lines represent the position of the co-rotation radius for the given APSW efficiency value and time snapshot. Note that the black and blue solid line are identical as the same initial model (t = 0.0 Myr) is used for this comparison. Panel c shows the evolution of the maximum disk temperature Tdisk,max at the disk’s midplane with respect to the simulation time t in years starting from t0 = 1 Myr. The line colors represent the same models as in panel a.

In the text

Fig. 6 Evolution of the stellar rotational period P* (panel a) and the maximum disk midplane temperature Γdisk,max (panel b) with respect to the simulation time t in years starting from t0 = 1 Myr. Different colors represent different initial stellar rotational periods: Pinit = 2 days (red line), Pinit = 4 days (magenta line), Pinit = 6 days (black line), Pinit = 10 days (blue line) and Pinit = 15 days (yellow line).

In the text

Fig. 7 Stellar spin evolution for different stellar spin evolutionary tracks of the layered disk reference model with respect to the simulation time t in years starting from t0 = 1 Myr. Panel a shows the stellar rotational period P* for the reference model MrefB with W = 2% (red line) as well as for models with W = 0% (black line) and W = 5% (blue line). Panel b shows the disk accretion rate Ṁdisk for the respective models.

In the text