{'fraction_non_alphanumeric': 0.04322114203488576, 'fraction_numerical': 0.025304460744744323, 'mean_word_length': 4.329374357056018, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-dependent, long-term hydrodynamic simulations of the inner protoplanetary disk II: The importance of stellar rotation

The spin evolution of young protostars, surrounded by an accretion disk, still poses problems for observations and theoretical models. In recent studies, the importance of the magnetic star-disk interaction for stellar spin evolution has been elaborated. The accretion disk in these studies, however, is only represented by a simplified model and important features are not considered. We combined the implicit hydrodynamic TAPIR disk code with a stellar spin evolution model. The influence of stellar magnetic fields on the disk dynamics, the radial position of the inner disk radius, as well as the influence of stellar rotation on the disk were calculated self-consistently. Within a defined parameter space, we can reproduce the majority of fast and slow rotating stars observed in young stellar clusters. Additionally, the back reaction of different stellar spin evolutionary tracks on the disk can be analyzed. Disks around fast rotating stars are located closer to the star. Consequently, the disk midplane temperature in the innermost disk region increases significantly compared to slow rotating stars. We can show the effects of stellar rotation on episodic accretion outbursts. The higher temperatures of disks around fast rotating stars result in more outbursts and a longer outbursting period over the disk lifetime. The combination of a long-term hydrodynamic disk and a stellar spin evolution model allows the inclusion of previously unconsidered effects such as the back-reaction of stellar rotation on the long-term disk evolution and the occurrence of accretion outbursts. However, a wider parameter range has to be studied to further investigate these effects. Additionally, a possible interaction between our model and a more realistic stellar evolution code (e.g., the MESA code) can improve our understanding of the stellar spin evolution and its effects on the pre-main sequence star.


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 Article number, page 1 of 12 arXiv:2208.08852v1 [astro-ph.SR] 18 Aug 2022 A&A proofs: manuscript no. main 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, 2008aRé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 magnetospheric 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 possi-ble. 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;Martin et al. 2012;Martin & Lubow 2013;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 0 = 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.

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);  and Gallet et al. (2019). Here, we want to give a recap of the respective model's key features. In Sec. 2.1, Sec. 2.2 and Sec. 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 Sec. 2.4.

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. 2007Zhu et al. , 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, L. Gehrig et al.: The importance of stellar rotation for protoplanetary disk evolution viscous accretion disk (e.g., Armitage et al. 2001), ∂ ∂t (Σ e) + ∇ · (Σ u e) + P gas ∇ · u 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 −1 ∂/∂ϕ) 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. 2004or Kluźniak & Rappaport 2007). The last term on the left-hand side of Eq. 2 is the vertically integrated magnetic pressure gradient and B z B 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 z B (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 (Sec. 2.4). Viscous angular momentum transport and viscous heating are represented by the viscous stress tensor Q in their corresponding terms in Eq. 2 and Eq. 3, respectively, which reads as, The factor µ denotes the vertically integrated dynamical viscosity, which yields, where α, c s , H p and ν depict the viscosity parameter , 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 P∆V work only close to the inner disk edge as the stellar field strength decreases with r −3 (see Sec. 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 P∆V 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 Sec. 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.

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;Kluźniak & Rappaport 2007) B z (r) = B R r 3 and (7) with as well as the co-rotation radius 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 magn (r trunc ) = max P ram (r trunc ), P gas (r trunc ) .
This yields for r in = r trunc the following description for the case P ram ≥ P gas (Hartmann et al. 2016), with the accretion rate onto the starṀ and a correction factor ξ = 0.7 (following Hartmann et al. 2016) accounting for nonspherical 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 trunc (P ram < P gas ) = B 1/3 R P gas 1/6 (8π) 1/6 .

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. (2002Romanova et al. ( , 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 2 (similar to Vorobyov et al. 2020) and is allowed to move during the simulation according to the disk evolution.

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 torquesJ = Γ ext . Following Matt et al. (2010), we assume rigid body rotation with a moment of inertia I = k 2 M R 2 with the normalized radius of gyration k (k 2 = 0.2, e.g., Armitage & Clarke 1996;Matt et al. 2010;. The temporal derivative of the stellar angular momentum J = I Ω readṡ The evolution of the stellar spin results from rewriting Eq. 15 as followṡ witḣ 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.

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

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 magnetospheric ejection (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 with the twist of magnetic field lines The additional assumption that the radius of the annulus of the disk connected to the MEs is of the order of r trunc (Gallet et al. 2019) yields the total torque felt by the star due to MEs, 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 r trunc in relation to the corotation radius r cor . If r trunc < K 2/3 rot 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 Sec. 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.

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; with r trunc Ω disk,in = u ϕ,in .

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 with the Alfvén radius where v esc = √ 2GM /R , m = 0.2177, K 2 = 0.0506 and K 1 = 1.7 arising as best fit values from magneto-hydrodynamic stellar wind simulations (see 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).

Numerical integration
To include the spin evolution model in the TAPIR code, we updated the stellar parameters after each timestep ∆t 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,Ω and the new stellar angular velocity Ω new = Ω old +Ω∆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.

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

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 Sec. 2, the initial stellar mass is M = 0.7M 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 −10 and 10 −7 M yr −1 for late-type star-disk systems. We therefore choose the initial accretion ratė M init = 1 · 10 −8 M /yr. The effects of a variation ofṀ init are covered in a further study. The viscous α-parameter ranges between 0.1 and ∼ 0.0001 (e.g., Zhu et al. 2007Zhu et al. , 2009bVorobyov & Basu 2009;Yang et al. 2018). We adopt initially a constant value of α = 0.005 used in Sec. 3.4 and Sec. 3.3, unless otherwise stated.
To generate episodic accretion outbursts (Sec. 3.5) due to magnetorotational instabilities (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 α that influences the viscous timescale of the accretion process; τ ν = r 2 /ν ∼ 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 selfconsistent 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.

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

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 P 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 2/3 rot is indicated by the dashed black line. During the first 2×10 5 yr, r trunc is close to r cor and r trunc /r cor > K 2/3 rot . Γ 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 4 years). 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 5 years. Thus, r trunc /r cor decreases until r trunc /r cor < K 2/3 rot 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.

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 Sec. 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 ratė M 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 1 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 2/3 rot . 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.

Fast and slow rotating stars
Within the predefined parameter space (Sec. 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 1 .
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 rot . 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. and W 2 . 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).

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 (Sec. 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 (Sec. 3.4.2). 2 The influence of the initial stellar period is studied in Sec. 3.4.2.

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 P (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 (u ϕ ) 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, max 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.

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

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 Sec. 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 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 6. Evolution of the stellar rotational period P (Panel (a)) and the maximum disk midplane temperature T disk, max (Panel (b)) with respect to the simulation time t in years starting from t 0 = 1 Myr. Different colors represent different initial stellar rotational periods: 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).

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 Sec. 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;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, multidimensional 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. 2019;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 disk can be modeled in a more self-consistent manner compared to previous stellar spin evolution studies (e.g., Matt et al. 2010;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 Fig. 5 and Fig. 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., 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.

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 Sec. 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 tempera-ture 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 Sec. 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., 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 0 to 10 4 years). 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.