© The Authors 2026

1. Introduction

The origin of the rotation rates seen for massive stars remains unclear. Surveys of rotation rates of O and B type stars typically show a bimodal distribution of rotation rates (Holgado et al. 2022; Britavskiy et al. 2023; Huang & Gies 2008; Huang et al. 2010), where fast rotators are assumed to have been formed by binary interaction (de Mink et al. 2013). Furthermore, Holgado et al. (2022) and Huang et al. (2010) estimated the initial distribution of rotation rates of O and B type stars at the beginning of the main sequence at around 10% to 20% of the critical rotation rate (i.e., the rotation rate at which gravito-centrifugal equilibrium is achieved). However, this critical rotation rate is reached very quickly when modeling protostars accreting through a disk (Haemmerlé et al. 2017) and during protostellar binary mergers (Oliva & Kuiper 2020). Thus, we must conclude that angular momentum removal mechanisms are at work over the course of pre-main sequence evolution.

Modeling the origin and evolution of stellar spins requires a coordinated description of the (proto)stellar interior, its surrounding material, and the ways in which angular momentum is transferred from the latter to the former. Matt & Pudritz (2005, 2008a,b) proposed accretion-powered winds as a solution to the stellar angular momentum problem in the case of low-mass stars. More recently, Gallet et al. (2019), Amard & Matt (2023) and Gehrig & Vorobyov (2023) linked star-disk interaction simulations by Ireland et al. (2022) to stellar evolution models; however, they did not include a detailed treatment of angular momentum transport within the star. Instead, they assumed a solid-body instantaneous redistribution. Moreover, those models require the presence of large-scale stellar magnetic fields. Only about 10% of massive stars exhibit such magnetic fields during their main sequence evolution (see, e.g., Donati & Landstreet 2009). Furthermore, an external convective envelope capable of sustaining a solar-like dynamo is not necessarily present across the entire pre-main sequence evolution of a massive star due to strong changes in the stellar structure (see Bernasconi & Maeder 1996; Hosokawa & Omukai 2009; Haemmerlé et al. 2017). Therefore, the principal mechanism responsible for angular momentum removal from the protostar while in the embedded phase must work independently of any coupling of the stellar magnetic field (if it exists) and the accretion disk. This is consistent with the findings of Holgado et al. (2022), where no strong magnetic braking effects have been observed in massive stars throughout their main sequence evolution (i.e., once the natal cloud has been dispersed).

In the context of low-mass star formation, simulations of star-disk interaction including disk- and stellar winds such as Romanova et al. (2005), Zanni & Ferreira (2013), Réville et al. (2015), Ireland et al. (2022) typically start from an initial model for the accretion disk and protostellar outflows. In recent years, however, it has been possible to model the formation of the disk and protostellar outflows starting directly from the gravitational collapse of a molecular cloud (Oliva & Kuiper 2023b; Mignon-Risse et al. 2021), including comparison against observations (Moscadelli et al. 2024; Sanna et al. 2025). However, an accurate description of the innermost ∼3 au in the cloud has been lacking until now.

In this Letter, we study the evolution of the spin of massive (proto)stars by performing and connecting two different computer models: a disk-jet ab initio simulation (Sect. 2) and a stellar evolution simulation (Sect. 3). To this end, we first calculated the mass and angular momentum transfer exerted by the disk and jet onto the protostar and then we introduced our results in the stellar evolution model.

Our main goal in this work is to test whether magnetocentrifugally driven protostellar jets are capable of keeping the protostar rotating well below the critical velocity, without any additional assumptions on the presence of large-scale stellar magnetic fields or of other stellar winds. A discussion of our results and conclusions are offered in Sect. 4.

2. Angular momentum removal by a protostellar jet

We performed two simulations of the formation of the accretion disk and the magnetically driven jet starting from the collapse of a slowly rotating molecular cloud (which we call the disk-jet ab-initio simulations). The circumstellar material was modeled using magnetohydrodynamics, including diffusive radiation transport, Ohmic dissipation, and self-gravity.

The simulations were initiated from identical initial conditions, so that we could obtain an accretion rate of ≈5 ⋅ 10⁻⁴ M_⊙ yr⁻¹. In an earlier work, we used one of the simulations (referred to as run A) to trace detailed streamlines of the magnetocentrifugal jet originating in the inner disk, leading to strong agreement with observations of water masers in the star-forming region IRAS 21078+5211 (Moscadelli et al. 2022). Run A has an inner boundary of radius 3 au and we assume that matter that enters this region is accreted by the protostar. To check the convergence of our results and understand the mass and angular momentum transfer from the inner disk onto the surface of the star, we reduced the inner boundary in run B down to a radius of r = 150 R_⊙ ≈ 0.7 au¹. This new reduced value is within the expected order of magnitude for the radius of the bloated phase of the protostar (see, e.g., Bernasconi & Maeder 1996; Hosokawa & Omukai 2009). More details on the setup of the simulations can be found in Appendix A, along with the convergence study given in Appendix B.

2.1. Mass and angular momentum flow along streamlines

We computed the angular momentum carried away by the jet by following several of its streamlines within a specific time window. To illustrate the process, we discuss how one of these streamlines ends up removing angular momentum from one point in the inner disk below. We note that we repeated the process for several launching points surrounding the protostar.

We started by focusing on the streamline shown in Fig. 1a, computed in the frame of reference co-rotating with the launching point, that is, the vector field, v − ω_lpRe_ϕ, where v is the velocity, ω_lp is the angular velocity of the launching point, R is the cylindrical radius, and e_ϕ is the unit vector in the azimuthal direction. The material was launched from the surface of the inner disk at a distance of R_lp = 3 au from the protostar and at a time t = 11.50 kyr after the start of the simulation². In this region, the magnetic field is strong enough and the density low enough for the flow to become sub-Alfvénic (corresponding to the interior of the yellow surface in Fig. 1a; see a more detailed discussion in Moscadelli et al. 2023, 2024; Oliva & Kuiper 2023b). As a consequence, the ejecta is forced to follow a magnetic field line while inside of the Alfvén surface and when viewed from the co-rotating frame of reference.

Fig. 1. Overview of the angular momentum transport done by the streamlines of the magneto-centrifugal jet, using the data from run A. (a) 3D view of one streamline launched from the inner disk at Rlp = 3 au, traced using the velocity field that co-rotates with the launching point. (b) Specific angular momentum transported by the streamline depicted in panel a. (c) Lever arm Λ of jet streamlines launched at r = 3 au from different launching angles within a time window. (d) Ratio ξ of the mass ejected through the Alfvén surface to the accreted mass. The calculation of the ejected mass for the light-blue line includes all streamlines of the jet and the gray line excludes those originated at the surface of the disk. The thick purple line indicates the median value considering the results of both calculations.

Along a streamline, the value expressed via

$$\begin{array}{c}\hfill L=r{v}_{\varphi }-\frac{B_r}{4\pi \rho v_r}r{B}_{\varphi }\end{array}$$(1)

is a constant of motion when the magnetic field only has radial and azimuthal components (see Appendix C for a discussion on how to arrive at this expression). We checked the simulations to ensure that in the jet, B_r dominates over B_θ and B_ϕ (see Oliva & Kuiper 2023b). To check that this quantity is a constant of motion also in our simulations, we calculated L, with the results shown in Fig. 1b. The specific angular momentum, rv_ϕ, increases as the ejecta are accelerated along the field line while the magnetic torque decreases, keeping L approximately constant along the streamline.

When the flow accelerates and becomes super-Alfvénic, it becomes capable of dragging the magnetic field lines with it and, thus, the angular momentum can be transferred onto large scales originating from the launching point. The torque, dJ/dt, enacted by the magnetic outflow can then computed by (see, e.g., Maeder 2009, Sect. 13.2) via

$$\begin{array}{c}\hfill \frac{\mathrm{d}J}{\mathrm{d}t}=\frac{2}{3}\frac{\mathrm{d}M}{\mathrm{d}t}{R}_{\mathrm{lp}}^2{\mathrm{\Omega }}_{\mathrm{lp}}{\left(\frac{R_A}{R_{\mathrm{lp}}}\right)}^2,\end{array}$$(2)

where dM/dt ≡ M_⊙e is the mass flow that crosses the Alfvén surface at a cylindrical radius, R_A (see Fig. 1a), R_lp is the cylindrical radius of the launching point, and Ω_lp is its angular velocity. The ratio R_A/R_lp ≡ Λ is known as the magnetic lever arm and expresses the amplification of the angular momentum outflow due to the magnetic torque. For the streamline shown in Fig. 1, we found that the Alfvén point is located at R_A = 26.8 au, yielding a lever arm of 8.93.

We repeated this process for many launch points situated at a distance of 3 au away from the protostar (limited by the inner boundary of our run A), using several launch angles, within the time window t ∈ [7, 11] kyr. We then computed the lever arm, with the results shown in Fig. 1c. The time average of the lever arm was between 10 and 20 for the different launch angles. In producing this result, however, we made the following simplification: since we do not need the azimuthal position of the Alfvén point, we traced only the poloidal magnetic field lines and computed their intersection with the Alfvén surface.

In our model of the torque, we express the ejected mass flow as a fraction ξ of the accretion rate Ṁ_a, that is, Ṁ_e ≡ ξ Ṁ_a. We integrated the mass flux that exits through the Alfvén surface (∫ρv ⋅ dA, where dA points outwards with respect to the Alfvén surface) for the same time window, resulting in an average mass outflow rate of ≈17% the mass inflow rate passing through the disk (see Fig. 1d). This calculation sums all the mass lifted by the jet, including the mass launched from the disk surface and the surroundings of the protostar. We also repeated the calculation while excluding mass from streamlines originating directly on the disk surface, yielding a time average of 0.3%. This value has to be considered a lower bound for ξ because our inner boundary condition, which prevents mass from leaving the central sink cell where the protostar forms. The median value, taking into account the results of both calculations, is ξ = 4.5%. We took this value as a reference for our models.

2.2. Discussion: extension of our calculations to the stellar surface

Strictly speaking, the surface of the protostar is not included within the computational domain of the disk-jet simulations. Next, we discuss how applicable these ideas are to the angular momentum transfer onto and from the surface of the forming protostar.

In low-mass stars, the torque resulting from the interaction of the (poloidal) stellar magnetic field and the magnetic field of the disk must be considered. Due to the fact that in massive protostars, the conditions for producing a stable solar-like dynamo are not fulfilled over their entire pre-main sequence evolution, we simply did not include in our models any poloidal stellar magnetic field. Instead, we examined the question of whether the magnetic fields of the jet and the cloud alone are already able to exert a sufficiently strong torque to keep the protostar below the critical surface rotation rate at all times. The scenario we investigated is similar to the accretion-powered wind described in Matt & Pudritz (2005); however, here we propose that the cloud magnetic field determines the Alfvén radius, R_A, instead of the stellar magnetic field and, thus, the lever arm can be ultimately linked to the conditions of the natal environment of the star (i.e., the properties of molecular clouds). This is because the gravitational collapse of the cloud naturally produces the kind of strong radial magnetic fields around the protostar that are at the core of how the (2) was derived. This happens a similar way to how a dipolar stellar magnetic field behaves in a magnetized stellar wind.

3. Angular momentum evolution of the protostar

The angular momentum accretion and its redistribution along the structure of forming massive star was investigated using the Geneva stellar evolution code (Genec, see Eggenberger et al. 2008). We modified the code to include the angular momentum extraction by the jet according to (2) as well as the angular momentum delivered by the Keplerian-like accretion disk. We started the calculations from the pre-main sequence evolution of a protostar and sustained the accretion of mass and angular momentum even after the onset of hydrogen burning, as expected during the formation of massive stars (see, e.g., Haemmerlé et al. 2017; Kuiper et al. 2010).

Motivated by the results of our ab initio disk-jet simulations, we set a constant accretion rate of Ṁ_a = 5 · 10⁻⁴ M_⊙/yr and chose two parametrizations of the jet, both with a lever arm of Λ = 20: a “strong” jet of ejection-to-accretion rate of ξ = 0.05 and a “weak” jet with ξ = 0.01 (cf. Fig. 1). Other values that are consistent with our disk-jet simulations will be discussed in future work; we chose only these two cases to explain the main idea for the purposes of this work. Once it is inside of the star, the accreted angular momentum is transported diffusively by convection, horizontal turbulence, and the secular shear instability.

The evolution of the surface angular velocity of the (proto)star is presented in Fig. 2. We defined the rotation criticality as Ω_surf/Ω_c. The critical angular velocity is computed as in Ekström et al. (2012, taking the deformation of the star as it spins into account), but for context, it is roughly ${\mathrm{\Omega }}_c=\sqrt{G{M}_{★}/R_{★}^3}$. Our results show that rotation criticality is determined by the strength of the jet, and remains constant in time for a given jet strength. This does not mean that the surface angular velocity is constant over time (Fig. 2b) because Ω_c also changes in time.

Fig. 2. (a) Evolution of the rotation criticality as a function of the instantaneous stellar mass under two different jet parametrizations. (b) Evolution of the surface and central angular velocity of the protostar as the structure changes from fully convective to fully radiative and then a convective core and radiative envelope. The opaque lines are running averages.

The rotation criticality is found to be independent of the stellar radius and the stellar structure. Even during the strong changes in stellar radius occurring during the bloating phase, the rotation criticality does not change. We interpret this result to be a consequence of the disk and jet torques reaching a rotational equilibrium very quickly, coupled with the fact that we used constant values of ξ and Λ. In Appendix D, we discuss this situation more carefully using a simple analytical argument and present further details on the setup of the stellar evolution simulations.

Our models reveal that even though a spin-up and a spin-down happen during changes in the moment of inertia of the star (e.g., during the bloating phase), the surface angular momentum is dominated by the accretion of the incoming material. The difference between the surface and stellar angular velocity (Fig. 2b) shows how the accreted angular momentum is then redistributed very efficiently in the case of a fully convective structure (resulting in solid body rotation) and less efficiently with the presence of a radiative zone (resulting in differential rotation).

4. Discussion and conclusions

The idea that an important part of the angular momentum removal process is carried out by accretion-powered winds is not a novel one (Gallet et al. 2019; Amard & Matt 2023; Matt & Pudritz 2008a,b; Gehrig & Vorobyov 2023, etc.). However, previous works have been focused on low-mass star formation and, thus, poloidal stellar magnetic fields were considered, which determine the magnetic lever arm (see, e.g., Ireland et al. 2022). From the results obtained in this work, we argue that the magnetic field of the cloud (i.e., the magnetic field of the jet) has an adequate structure and capacity to independently produce a sufficient level of magnetic braking to keep the protostar below the critical rotation rate at all times as long as there is accretion (or that the same magnetic field structure remains attached to the stellar surface). Moreover, Fig. 2 shows that in principle, with a sufficiently strong jet, we can obtain any surface rotation rate at the end of the star formation process. The jet strength is a direct result of the initial conditions of the gravitational collapse of the cloud (Oliva & Kuiper 2023b); thus, we can link the rotation rate of a massive star just after the formation process directly to its natal environment. According to our results, during the accretion epoch and just after it ends, the rotation rate of massive stars is not, in principle, correlated to stellar properties such as mass.

However, in reaching this conclusion, we did not consider the possible additional effects of spin-up or spin-down due to binary interactions, nor the additional spin-down effects expected from the interaction of any existing (proto)stellar magnetic field or the disk. Additionally, our 2D setup ends up neglecting 3D effects in angular momentum transport (i.e., breaking of axial symmetry due to disk fragmentation or turbulence).

In future works, we will address the dependence of the angular momentum removal carried out by the jet with the initial conditions for star formation, how this determines the rotation rate of the formed massive star, and what the roles of stellar feedback and metallicities are.

Acknowledgments

AO has received support from ESKAS No 2023.0405 and the University of Costa Rica project No 829-C6-133. LS and SE have received support from the SNF project No 212143. PE have received support from the SNF project No 219745.

References

Appendix A: Disk-jet simulations

The disk-jet simulations start from the collapse of a gravitational cloud initially containing 100 M_⊙ enclosed within a radius of 0.1 pc. This initial configuration produces an accretion rate of ≈5 ⋅ 10⁻⁴ M_⊙ yr⁻¹ including the effects of the outflows (see Fig. A.1 for more details). The motion of the gas is modeled with the methods of non-ideal magnetohydrodynamics using the Pluto code (Mignone et al. 2007) with the module Haumea for self-gravity (Kuiper et al. 2010) and the Machida et al. (2007) prescription for Ohmic resistivity as a function of density only (with a fixed temperature of 10 K). Radiation transport is taken into account using the gray flux-limited diffusion approximation using the module Makemake as described in Kuiper et al. (2020). More details on the methods used in the disk-jet simulations can be found in Oliva & Kuiper (2023a) and Oliva & Kuiper (2023b).

Fig. A.1. (a) Mass accretion rate and (b) stellar mass as a function of time using the data from run A.

Initially, the density in the cloud is distributed according to the power law ρ ∝ r^−1.5. It is slowly rotating as a solid body with its total rotational energy being 4% of its total gravitational energy. There is an initially uniform magnetic field of 68.31 μG oriented parallel to the rotation axis, which corresponds to a mass-to-flux ratio of 20 times the critical (collapse-preventing) value.

We used a time-independent, 2D, axisymmetric grid in spherical coordinates: the radial coordinate increases logarithmically with distance from the center, and the polar coordinate increases linearly from the rotation axis until the midplane, where we assume equatorial symmetry. Two different grids were used in this study. The first case (run A) corresponds to the fiducial run with the grid labelled as x8 in Oliva & Kuiper (2023b)³, that is, a radial coordinate that extends from r = 3 au to r = 0.1 pc divided in 448 cells, and a polar angle divided in 80 cells. In the second case (run B), the grid extends from r = 0.698 au ≈ 150 R_⊙ to 2500 au, while maintaining the same initial density distribution than in the first case. In run B, there are 160 cells in the radial coordinate and 32 cells in the polar coordinate.

The goal of this study is to understand the angular momentum removal done by the magnetically driven outflows, so no stellar feedback is taken into account. This approximation is valid for the early stages of the formation of a massive star. We will address the effects of stellar feedback in a future study, but we anticipate that they will increase the angular momentum removal done by the total protostellar outflow and thus help further reducing the rotation rate of the (proto)star. At any point in time, the disk-jet simulation computes the instantaneous stellar mass as the mass that has entered the sink cell, and interpolates the stellar evolution tracks of Hosokawa & Omukai (2009) to compute the stellar radius. This is the way we obtained the stellar radii quoted in Sect. 2, but not in Sect. 3, which were computed directly with the Geneva stellar evolution code.

Appendix B: Convergence of our results with sink cell size

Since the inner boundary of run A is located at r = 3 au, far away still from the stellar surface, we produced run B to study the convergence of the results of Sect. 2.1 when the launching points get closer to the star. A reduced sink cell size dramatically increases computing time due to the decrease in the orbital timescale in the inner disk. This meant that we were able to compute the evolution of the disk-jet system only up to the initial stages of jet launch, that is, until 7.5 kyr after the onset of the gravitational collapse. At this stage, the Alfvén surface has not yet grown to its full size, and thus the Alfvén radii are smaller than in the time window considered in the rest of the work.

We repeated the calculations of the lever arm for several launch angles around the protostar, for both run A and run B, in the time window t ∈ [5.7,  7.5] kyr (the results are shown in Fig. B.1). The reduction in sink cell size made for run B does not significantly affect the average value of the lever arm. In run A, the jet is narrower than in run B (therefore, there are only launching points for θ_launch ≲ 30° in run A). In run B in turn, the smaller sink cell means that the inner disk is better resolved in that time window, and thus the jet is wider in comparison to run A as the disk also launches streamlines into the jet. We conclude from this exercise that the lever arm values deduced from run A with a launch radius of 3 au hold for smaller launch radii, and therefore it is reasonable to use them to compute the wind torque down to the stellar surface.

Fig. B.1. Convergence of the lever arm with sink cell size.

Appendix C: Angular momentum transported by a magnetized wind

Consider an axisymmetric, fully ionized plasma with density ρ, velocity v = v_re_r + v_ϕe_ϕ and magnetic field B = B_re_r + B_ϕe_ϕ surrounding a star (subindex ★). In the steady state, the momentum equation in the azimuthal direction is

$$\begin{array}{cc}& \rho {[(\mathit{v}·\mathrm{\nabla })\mathit{v}]}_{\varphi }=\frac{1}{4\pi }{[(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}]}_{\varphi }\hfill \end{array}$$(C.1)

$$\begin{array}{cc}& \Rightarrow \frac{d}{\mathit{dr}}(r{v}_{\varphi })=\frac{B_r}{4\pi \rho v_r}\frac{d}{\mathit{dr}}(r{B}_{\varphi }).\hfill \end{array}$$(C.2)

The solenoidality of the magnetic field (∇ ⋅ B = 0) implies that

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\partial }{\partial r}(r^2{B}_r)=0\Rightarrow r^2{B}_r=\mathrm{const},\end{array}$$(C.3)

and the conservation of mass implies that

$$\begin{array}{c}\hfill r^2\rho v_r=\mathrm{const}.\end{array}$$(C.4)

We can use (C.3) and (C.4) to integrate (C.2) to obtain a constant of motion,

$$\begin{array}{c}\hfill r{v}_{\varphi }-\frac{B_r}{4\pi \rho v_r}r{B}_{\varphi }=\mathrm{const}\equiv L.\end{array}$$(C.5)

The first term in this equation is the specific angular momentum and the second term is the torque exerted by magnetic stresses. Using the frame which co-rotates with the launching point, and introducing the definition of the Alfvén speed, we find L to be (Weber & Davis 1967)

$$\begin{array}{c}\hfill L={\mathrm{\Omega }}_{\mathrm{surf}}{R}_A^2,\end{array}$$(C.6)

where Ω_surf is the angular velocity of the surface of the star (which is the same angular velocity as the launching point) and R_A is the radius where the flow switches from sub-Alfvénic to super-Alfvénic. From this, we conclude that the torque done onto the star by the flow should have the form

$$\begin{array}{c}\hfill \frac{d{J}_{★}}{\mathit{dt}}\sim {\mathrm{\Omega }}_{\mathrm{surf}}{R}_A^2\frac{\mathit{dM}}{\mathit{dt}}\end{array}$$(C.7)

as proposed in, for example, Weber & Davis (1967) and Maeder (2009, Sect. 13.2). The full derivation of this expression (see Mestel 1984) contains an integral of the streamlines originating in the surface of the star and which cross the Alfvén surface (to account for the different latitudes of origin, we multiply an extra geometrical factor of 2/3 like in Weber & Davis 1967; Matt & Pudritz 2005). Even though some streamlines of the jet originate close to the poles, they can still carry a significant angular momentum flow, as long as the Alfvén radius reached is sufficiently large.

Appendix D: Stellar evolution models

We modified the Geneva stellar evolution code to include the (spin-down) torque done by the protostellar jets using (2). Introducing the definitions of ξ and Λ, and taking the launching point as an average across the surface of the star (see Gallet et al. 2019), we obtain the prescription

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{jet}}=\frac{2}{3}\xi {\dot{M}}_a{\mathrm{\Lambda }}^2{R}_{★}^2{\mathrm{\Omega }}_{\mathrm{surf}},\end{array}$$(D.1)

where we take $R_{★}=\sqrt{{\mathcal{L}}_{★}/(4\pi \sigma )}/T_{\mathrm{eff}}^2$ (ℒ_★ being the stellar luminosity, σ the Stefan-Boltzmann constant and T_eff the effective temperature) and Ω_surf is the angular velocity of the surface of the star just before it is entrained by the jet⁴. The (spin-up) torque delivered by the disk onto the protostar is included as (Haemmerlé et al. 2017)

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{disk}}={\dot{M}}_a\sqrt{G{M}_{★}{R}_{★}}.\end{array}$$(D.2)

The net torque onto the stellar surface is then $\dot{J}={\dot{J}}_{\mathrm{disk}}-{\dot{J}}_{\mathrm{jet}}$. As noted by Matt & Pudritz (2005), a consequence of writing the ejection rate as proportional to the accretion rate, and the Alfvén radius as proportional to the stellar radius is the following.

If we assume that the mass and radius of the star change slowly, then the accretion torque (D.2) is roughly constant. However, the jet torque (D.1) depends on Ω_surf, which means that there is a value of Ω_surf for which ${\dot{J}}_{\mathrm{disk}}={\dot{J}}_{\mathrm{jet}}$, that is, rotational equilibrium is produced. The equality of the torques leads very easily to the expression

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{surf}}=\frac{3}{2}\frac{1}{\xi \mathrm{\Lambda }}{\mathrm{\Omega }}_c,\end{array}$$(D.3)

which means that a constant ξ and Λ produce a constant rotation criticality. In our disk-jet and stellar evolution models, we wish to test whether we reach a rotational equilibrium where Ω_surf < Ω_c.

Rotation inside of the star is treated with the shellular rotation hypothesis, that is, that the angular velocity is constant on isobars. The angular momentum is transported within the star with a diffusive scheme, neglecting advection by meridional currents. We include the diffusion of angular momentum by horizontal turbulence with the diffusion coefficient of Zahn (1992, which depends on the velocity of the meridional current). Additionally, we consider the effects of secular shear turbulence with the diffusion coefficient by Maeder (1997). A large diffusion coefficient models the efficient angular momentum transport in convective regions (Eggenberger et al. 2008).

We start the calculation with the adiabatic contraction of a protostellar seed of 1.5 M_⊙ in very slow rotation forming from material with solar metallicity (i.e., star formation in the present-day universe). After some tests we will present in future work, we concluded that the mass and angular velocity of the seed do not significantly change our results; this is in agreement with previous findings by Gehrig & Vorobyov (2023).

As a test of the stellar angular momentum problem, we ran the model with a weak jet and switched it off at an arbitrary time, in this case just after the protostar has reached a mass of 2 M_⊙. The results of this test are shown in Fig. D.1. The seed for the Genec model is built from a polytrope, therefore it does not include the effects of accretion or rotation. Because of this, when we switch on accretion and rotation at the beginning of the calculation, we let the model redistribute instantly the accreted angular momentum as a solid body for a short time, while thermal adjustment is achieved. After this, the angular momentum is transported diffusively in the stellar interior and the disk and jet set the rotation rate of the protostar as discussed in Sect. 3. When the protostar reaches a mass of 2 M_⊙, we switched off the jet and allowed only angular momentum accretion transferred from the Keplerian-like disk. The result is that the protostellar surface reaches very quickly the critical speed: at the assumed accretion rate of 5 ⋅ 10⁻⁴ M_⊙ yr⁻¹, only ∼0.1 M_⊙ could be accreted before the critical limit. This result supports our conclusion that the angular momentum of the stellar surface is completely dominated by the angular momentum of the accreted material.

Fig. D.1. Example of the stellar angular momentum problem using Genec: after switching off the jet, the protostar reaches very quickly the critical speed.

The general features of the structure and evolution of the (proto)star are examined in detail in Fig. D.2. The Kippenhahn diagram (Fig. D.2a) reveals the locations and times at which the transfer of energy (and angular momentum) is radiative and convective. The strong convective motions present in the entirety of the protostar at the beginning of its evolution permit a very efficient redistribution of the angular momentum and produce the solid-body rotation profile seen in Fig. D.2c for M_★ = 2.87 M_⊙. After ≈11.8 kyr of evolution, shell deuterium burning begins and with it, the bloating phase. The star becomes radiative and the angular momentum is not efficiently distributed within the structure, giving rise to the differential rotation profile of Fig. D.2c for M_★ = 12 M_⊙.

Fig. D.2. Overview of the evolution of the protostar. (a) Kippenhahn diagram showing deuterium and hydrogen burning, with the red line being the total instantaneous stellar mass and the shadowed regions indicate convection. An element is considered to be burning when the energy generation rate is at least 102 erg g s−1. The dotted lines indicate the limits of the region where the corresponding element is burned, and the solid lines indicate the center of the burning region for each element (colors). (b) Evolution of the stellar radius. (c) Internal distribution of the stellar angular velocity for different times (expressed as instant stellar masses). (d) Hertzsprung–Russell diagram of the accreting protostar (the bottom right corresponds to the initial time).

During the bloating phase, the star spins down (Fig. 2 lower panel), but because the stellar radius increases, the critical angular velocity also decreases. After the Kelvin-Helmholtz contraction at ≈38 kyr, the star starts burning hydrogen and develops a convective interior. During this time, the star spins up (although not substantially), but the stellar radius decreases, which means that the critical angular velocity also increases. A more massive star has a higher critical angular velocity than a lower mass star, so for a fixed criticality, the rotation rate of the star when it reaches hydrogen burning is larger than in the pre-main sequence phase. At this stage (M_★ = 40.3 M_⊙ in Fig. D.2c), the angular momentum structure inside of the star is also distributed depending on whether the regions are convective (leading to flat regions in the curve, that is, solid body rotation) or radiative (where differential rotation can develop). Our models reveal that despite the changes in radius and mass of the star, a sub-critical rotational equilibrium is found and maintained, and that the rotation criticality is controlled by ξ and Λ. This is a consequence of the surface angular momentum being dominated by the angular momentum of the accreting material and that any spin-up of the star also makes ${\dot{J}}_{\mathrm{jet}}$ more efficient.

All Figures

Fig. 1. Overview of the angular momentum transport done by the streamlines of the magneto-centrifugal jet, using the data from run A. (a) 3D view of one streamline launched from the inner disk at Rlp = 3 au, traced using the velocity field that co-rotates with the launching point. (b) Specific angular momentum transported by the streamline depicted in panel a. (c) Lever arm Λ of jet streamlines launched at r = 3 au from different launching angles within a time window. (d) Ratio ξ of the mass ejected through the Alfvén surface to the accreted mass. The calculation of the ejected mass for the light-blue line includes all streamlines of the jet and the gray line excludes those originated at the surface of the disk. The thick purple line indicates the median value considering the results of both calculations.

In the text

	Fig. 2. (a) Evolution of the rotation criticality as a function of the instantaneous stellar mass under two different jet parametrizations. (b) Evolution of the surface and central angular velocity of the protostar as the structure changes from fully convective to fully radiative and then a convective core and radiative envelope. The opaque lines are running averages.
In the text

	Fig. A.1. (a) Mass accretion rate and (b) stellar mass as a function of time using the data from run A.
In the text

	Fig. B.1. Convergence of the lever arm with sink cell size.
In the text

	Fig. D.1. Example of the stellar angular momentum problem using Genec: after switching off the jet, the protostar reaches very quickly the critical speed.
In the text

Fig. D.2. Overview of the evolution of the protostar. (a) Kippenhahn diagram showing deuterium and hydrogen burning, with the red line being the total instantaneous stellar mass and the shadowed regions indicate convection. An element is considered to be burning when the energy generation rate is at least 102 erg g s−1. The dotted lines indicate the limits of the region where the corresponding element is burned, and the solid lines indicate the center of the burning region for each element (colors). (b) Evolution of the stellar radius. (c) Internal distribution of the stellar angular velocity for different times (expressed as instant stellar masses). (d) Hertzsprung–Russell diagram of the accreting protostar (the bottom right corresponds to the initial time).

In the text