© The Authors 2023

1. Introduction

Over the next decades, stellar structure and evolution (SSE) theory will begin to advance from one-dimensional models towards solving the stellar structure equations in three spatial dimensions. So far, all SSE codes rely on the approximation of spherical symmetry as this greatly simplifies the numerics. But asteroseismic and interferometric studies of stars on the main sequence (MS) have revealed that a significant fraction of stars rotate with velocities where the approximation of spherical symmetry is not justified. Many SSE codes that account for rotational effects are based on the pioneering works by Zahn (1992) and Chaboyer & Zahn (1992), and sequential works (e.g., Talon & Charbonnel 1998; Palacios et al. 2003; Mathis et al. 2004), which provide expressions for the efficiency of chemical mixing due to both shear-driven turbulence and advective transport (cf. Ekström et al. 2012). As this paradigm is constructed for one dimension, several assumptions have been made regarding the relative strengths of the horizontal and vertical diffusion, and the validity domain of these assumptions is currently unknown, except that they should apply in the limit of slow rotation.

Attempts at solving the stellar structure of rotating stars in more than one dimension began in the sixties, but it only recently became possible to compute the stellar structure in a self-consistent manner in two dimensions, that is, with the ESTER¹ code (Espinosa Lara & Rieutord 2013; Rieutord et al. 2016). These models have been restricted to a steady state, where the hydrogen mass fraction profile is modelled as a step function, namely a constant value in the convective core and a constant value in the radiative envelope (see e.g., Gagnier et al. 2019a,b; Bouchaud et al. 2020; Howarth et al. 2023). Furthermore, ESTER is designed for early-type stars and the applicability is currently limited to stars with a convective core and radiative envelope.

In this Letter, we present the results of a novel approach whereby the temporal dimension is added to the two-dimensional steady-state models by solving, in addition, the equation of chemical evolution. The overall aim of the present work is to demonstrate the capabilities of the state-of-the-art 2D evolution ESTER models of a rotating massive star, and present predictions for the rotation and chemical profiles. This includes fully solving the large-scale velocity field inside the star. We successfully performed 2D evolution of a slow and a moderate rotator with a mass of 12 M_⊙, from the zero-age main sequence (ZAMS) up to near the terminal-age main sequence (TAMS). For the predicted rotation profile, we also aim to make a comparison with observations. In Sect. 2, we describe the fundamentals of the ESTER code and the improvements made. In Sect. 3, we compare the predictions of the 2D models, especially the predicted rotation profile, with asteroseismic measurements of a massive star. We show that the predictions of the 2D model are consistent with the observations. In Sect. 4, we discuss the contribution of meridional circulation to nitrogen enrichment as predicted by 2D models. Finally, we conclude in Sect. 5.

2. ESTER models

In this work, we make use of new ESTER models where time evolution has been implemented. We recall that ESTER models are 2D models of rotating stars that include the centrifugal distortion of the star and the associated baroclinic flows. These models therefore predict the structure (pressure, density, and temperature distributions) and the associated large-scale flows, namely the differential rotation and meridional circulations. The first version of the ESTER code, which computes the steady states of early-type star models (Espinosa Lara & Rieutord 2013), has now been completed with a time evolutive version that includes the unsteady terms needed to follow the thermal and nuclear evolution of a star. The code solves the set of equations given in Appendix A, still assuming the axisymmetry of the star. Additionally, in this new version, the evolution of the mass fraction of hydrogen is solved,

$$\begin{array}{c}\hfill \frac{\partial X}{\partial t}+\mathit{v}·\mathrm{\nabla }X=\frac{1}{\rho }\mathrm{\nabla }(\rho [D]\mathrm{\nabla }X)+{\dot{X}}_{\mathrm{nuc}},\end{array}$$(1)

where ρ is the density, v is the meridional velocity, and [D] is a tensor constructed from the horizontal and vertical chemical diffusion coefficients.

As for the steady-state version, we are using OPAL for opacities and equation of state (Rogers et al. 1996). Nuclear energy production is modelled via a simple law for the CNO cycle (see appendix) from which we derive Ẋ_nuc.

The spatial discretization is based on a spectral element method where spectral elements are spheroidal shells bounded by isobars Rieutord et al. (2016). Typically, a model uses 12 spheroidal shells with 30 points in a Gauss–Lobatto grid radially and 24 points in a Gauss–Legendre grid in latitude. Time evolution is insured by a first-order backward Euler method where the time step is manually set to 0.5 Myr, and is decreased when the model fails to converge.

Following the theoretical work of Zahn (1992) on the predicted diffusion constant resulting from rotationally induced chemical mixing, and the 1D implementation by Mombarg et al. (2022), the vertical chemical diffusivity is taken as,

$$\begin{array}{c}\hfill D_{\mathrm{v}}=\eta {⟨N_0^2⟩}_V^{-1}{⟨K{r}^2{(\mathit{n}·\mathrm{\nabla }\mathrm{\Omega })}^2⟩}_{\theta }.\end{array}$$(2)

Here, K is the thermal diffusivity, n is the unit vector normal to the isobars, and Ω is local angular velocity. The term ${\langle N_0^2\rangle }_V$ is the volume-averaged squared Brunt–Väisälä frequency of the initial steady-state model, where the average is taken over the volume where $N_0^2>0$. We therefore smooth out the rapid radial variations of N² that often give rise to numerical difficulties. Furthermore, η is a free non-dimensional parameter, which we adjust to avoid numerical instabilities arising if chemical diffusion is too low, while keeping it close to unity. For model M1, η needs to be increased to compensate for the smaller shear compared to model M2. The values of η are given in Table C.1. Furthermore, we take the angular average of the term between ⟨ ⋅ ⟩_θ in Eq. (2) and assume a fixed horizontal diffusion coefficient D_h = 10⁵ cm² s⁻¹. In the ESTER evolution models, there is no ad hoc enhanced mixing at the core boundary (overshooting) as is typically applied in 1D evolution models to account for the discrepancy between predicted and observationally inferred core masses.

We successfully ran 2D evolution models for a 12 M_⊙ star: once for an initial rotation (at the equator) of 15% the critical angular velocity Ω_c (model M1, v_eq = 112 km s⁻¹) and once for 50% (model M2, v_eq = 358 km s⁻¹)². The evolution models start from a steady-state model at the ZAMS with a uniform chemical composition with X = 0.71 and Z = 0.012 (no assumption of solid-body rotation). The evolution in the Hertzsprung-Russell diagram (HRD) of the models computed with ESTER is shown in Fig. 1 for the effective temperature at the pole and at the equator. Figure 2 shows the angular velocity distribution at 1 Myr and 15.75 Myr for model M1. This stellar model starts with equatorial regions rotating more rapidly than polar ones, namely solar-like, and gradually evolve towards a more shellular rotation with a slightly anti-solar surface rotation. The same behaviour is also seen in M2 (see Appendix B). Contrary to an evolution modelled with a series of steady-state models, where one decreases the hydrogen mass fraction in the core (Gagnier et al. 2019b), our model evolution shows that Ω/Ω_c decreases as the star evolves. We understand this behaviour to be a consequence of the slow redistribution of AM through baroclinic modes, which are damped on a timescale similar to the nuclear one. The timescale on which baroclinic modes are damped is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{baro}}={\langle \frac{N^2}{{\mathrm{\Omega }}^{2}K}\rangle }_V\frac{d_{\mathrm{env}}^2}{({\pi }^2+4){\pi }^2},\end{array}$$(3)

Fig. 1. HRD showing the evolution tracks of the 2D models for a 12 M⊙ star (Z = 0.012, Xi = 0.71) for (Ωeq/Ωc)i = 0.15 (dashed-dotted), and 0.5 (dashed). The black lines represent the properties at the pole, the grey lines the properties at the equator. The tracks are terminated slightly before the TAMS when the solver can no longer converge, which occurs at Xc/Xini = 0.141 and 0.131 for the respective aforementioned rotation rates. The dots mark the model of HD 192575 that is discussed in Sect. 3.

Fig. 2. Map of the angular velocity as a fraction of the critical angular velocity for a model of 1 Myr (top panel) and 15.75 Myr (bottom panel). The plots show cuts in the meridian plane. These plots are for a 12 M⊙ star with (Ωeq/Ωc)i = 0.15.

following the work of Busse (1981). The τ_baro is usually of the order of the Eddington–Sweet timescale (Rieutord et al. 2006b). In its expression, d_env is the thickness of the radiative envelope in the polar direction. The nuclear evolution timescale, τ_evol = X_c/Ẋ_c, is found to be roughly 3 and 30 times larger at the start of the MS for M1 and M2, respectively. Hence, and especially for M2, baroclinic modes that may be excited by initial conditions can be damped during nuclear evolution. We can therefore assume that initial conditions are of little importance and are forgotten during the first part of the MS. Yet, at the end of the MS, the ratio of timescales is reversed: the damping of baroclinic modes happens on a timescale that is 100 (M1) and 10 (M2) times longer than the nuclear one. This implies that the dynamical evolution, and especially the rotation rates, cannot be computed as a succession of stationary states as in Gagnier et al. (2019b). Table C.1 lists the values of these time scales for our two models. We note that the growth of the stellar radius is the main contributor to the growth of τ_baro with MS evolution.

Figure 3 shows the evolution of the vertical diffusion coefficient as described by Eq. (2) for model M1. At the start of the MS, the latitudinally average angular velocity $\overline{\mathrm{\Omega }}(r)$ decreases from the core towards the surface until roughly 60% of the fractional radius, after which $\overline{\mathrm{\Omega }}$ increases towards the surface (see right panel of Fig. 4), creating a layer where the shear is weak and thus D_v reaches a local minimum. Further along the MS evolution, $\overline{\mathrm{\Omega }}$ keeps decreasing all the way towards the surface. Figure 3 also shows the difference between taking an average value for N² as we do and using the full profile. This latter case shows the rapid variations of D that occur near the core and give rise to numerical problems as mentioned above.

Fig. 3. Profile of the vertical diffusion coefficient according to Eq. (2) for a 12 M⊙ model with an initial rotation frequency of 15% the critical angular velocity. The grey dotted lines show $D_{\text{v}}{\langle N_0^2\rangle }_V/(\eta N^2)$ at a few ages (from left to right: Xc/Xini = 0.14, 0.30, 0.58, 0.82, 0.98), that is, if we were to take the local value of the Brunt-Väisälä frequency instead of the volume averaged value.

3. Calibration to HD 192575

To calibrate our rotating models, we use the results of Burssens et al. (2023, hereafter B23), who performed an asteroseismic modelling of the β Cephei pulsator HD 192575 based on 1D non-rotating stellar evolution models. This star is a unique case with which to test the theory of angular momentum transport as it has a relatively precise age estimate, and an inferred rotation profile from the rotational splittings of the observed multiplets. The stellar parameters of HD 192575 derived by B23 are summarised in the middle column of Table 1. We computed ESTER evolution tracks with a mass (M_⋆), metallicity (Z), and initial hydrogen-mass fraction (X) fixed to the values found by B23. We then picked the model for which, after evolution, the hydrogen-mass fraction in the convective core is closest to that inferred by these authors. As the shape of the internal rotation profiles in massive stars is currently unknown, these authors modelled the rotational splittings by assuming a rotation profile that is described as

$$\begin{array}{c}\hfill \mathrm{\Omega }(r)=\{\begin{array}{cc}{\mathrm{\Omega }}_{\mathrm{core}}\hfill & r\le r_{\mathrm{core}},\hfill \\ {\mathrm{\Omega }}_{\mathrm{core}}-\mathrm{\Delta }\mathrm{\Omega }\frac{r-r_{\mathrm{core}}}{r_{\mathrm{shear}}-r_{\mathrm{core}}}\hfill & r_{\mathrm{core}}<r<r_{\mathrm{shear}},\hfill \\ {\mathrm{\Omega }}_{\mathrm{surf}}\hfill & r\ge r_{\mathrm{shear}},\hfill \end{array}\end{array}$$(4)

where ΔΩ = Ω_core − Ω_surf, r_core is the radius of the convective core, and r_shear is the outer radius of the shear zone where the rotation frequency decreases linearly until it reaches the surface value. The values of Ω_core and Ω_surf are then optimised to best reproduce the observations for a given assumption for r_shear. However, in the ESTER models, the rotation profile is computed self-consistently. Figure 4 shows the predicted evolution of the rotation profile for model M1. As can be seen from this figure (left), the latitudinal differential rotation changes from solar-like to anti-solar. This change may be interpreted as follows. The solar-like latitudinal profile is indeed the relaxed steady baroclinic state of a rotating radiative envelope (Espinosa Lara & Rieutord 2013). As time evolution proceeds, the core shrinks and spins up due to angular momentum conservation. Thanks to the Taylor–Proudman theorem, which states that the velocity field cannot vary in the direction of the rotation axis³, we understand that polar regions tend to follow the core and therefore rotate more rapidly than equatorial regions. This is a consequence of the rapid nuclear evolution, which prevents the star from relaxing to a quasi-steady state.

Fig. 4. Evolution of the rotation profile of model M1. Left panel: predicted surface rotation as a function of colatitude throughout the evolution (indicated by colour) normalised by the rotation frequency at the pole. Right panel: profile of the angular velocity (averaged over θ) throughout the evolution. The thick coloured line indicates the location of the convective core boundary. The black lines in both panels correspond to the best-matching model for HD 192575 (last column in Table 1).

Figure 4 (right panel) also shows the ‘radial’ rotation profile averaged over the latitude. This latter is similar in appearance to the linear piece-wise profile given in Eq. (4) when r_shear is taken equal to the outer radius of the region with a non-zero gradient in the mean molecular weight (μ). This behaviour appears to be independent of the initial rotation, because model M2 shows the same profile. Therefore, we compare our results with the core and surface rotation frequencies derived with this assumption. The value of Ω_core in the ESTER models is defined as the value of $\overline{\mathrm{\Omega }}(r_{\mathrm{core}})$, where an average is taken over all points in latitude. The last column in Table 1 shows the parameters of the ESTER M1 – model. The predicted core- and surface rotation of this model are consistent with the observationally derived values for HD 192575. Moreover, this model is able to reproduce, at the inferred age, the core mass, the asteroseismic radius from the 1D modelling, the astrometric luminosity from Gaia (Gaia Collaboration 2023), and the spectroscopically derived effective temperature within the uncertainties quoted in B23. It should be noted that with the implemented rotational mixing, the age is consistent with the age predicted from 1D (MESA) models with core-boundary mixing, although B23 were not able to precisely constrain its efficiency. In summary, the 2D evolution model M1 is able to explain the measured core and surface rotation of HD 192575 at the asteroseismically inferred age.

Additionally, we computed 1D SSE models with MESA (r22.11.1; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) using the same physics as B23, except that we also include rotation in our models (more details in Appendix D). The rotation profiles predicted by the 1D MESA models are similar to taking r_shear = R_⋆ (Eq. (4)). Therefore, we compare the MESA profiles with the core and surface rotation frequencies derived by B23 with this assumption. When the same uniform viscosity is assumed as the one used in the ESTER models, we find that the 1D MESA model for (Ω_eq/Ω_c)_i = 0.15 does not match the core and surface rotation rates. To reproduce the measured core and surface rotation frequencies at the current age with the 1D MESA model, the viscosity has to be increased to 10⁹ cm² s⁻¹ compared to the value of 10⁷ cm² s⁻¹ used in the ESTER models (vertical and horizontal components). While both the 1D MESA and 2D ESTER models can reproduce the core and surface frequencies of HD 192575, the shape of the rotation profile is significantly different between the two (see Fig. 5). Future studies on asteroseismic rotation inversion might be able to rule out one of these profiles (Vanlaer et al. 2023).

Fig. 5. Predicted rotation profile of HD 192575 with ESTER (red solid line) and MESA (black lines). The profiles indicated by a solid line and dotted line are for a constant viscosity of 107 cm2 s−1, while the dashed line is for 109 cm2 s−1. The dots correspond to the location of the core boundary. The red and grey shaded areas indicate the measured core (top one) and surface (bottom one) rotation frequencies by Burssens et al. (2023), assuming rshear equal to the outer boundary of the μ-gradient zone, and rshear = R⋆, respectively.

4. Nitrogen enhancement at the surface

The observed abundance of ¹⁴N at the surface is often used to probe the efficiency of (rotational) chemical mixing in massive stars (e.g., Brott et al. 2011). The nuclear reactions in the ESTER models are described by analytical formulae for the energy generation rates of the pp-chain and CNO cycle (Kippenhahn & Weigert 1990; Rieutord et al. 2016). To track the abundance of ¹⁴N as a result of hydrogen-burning via the CNO cycle, we make the following assumptions. First, at the start of the evolution, all ¹²C present in the core is converted to ¹⁴N. The reaction that controls the rate at which ¹⁴N is generated is therefore the proton capture of ¹⁶O to create ¹⁷F, which then rapidly decays into ¹⁷O, itself reacting on protons to yield ¹⁴N and ⁴He. Therefore, the evolution of the mass fraction of ¹⁴N due to nuclear reactions is described by

$$\begin{array}{c}\hfill \rho \frac{\mathrm{d}X{(}^{14}\mathrm{N})}{\mathrm{d}t}=\lambda {(}^{16}\mathrm{O}{\to }^{17}\mathrm{F})n{(}^{16}\mathrm{O})n(\mathrm{H})m({}^{14}\mathrm{N}),\end{array}$$(5)

where n(⋅) denotes the number densities, m the atomic mass, and λ(¹⁶O→¹⁷F) is the Maxwellian-average reaction rate ⟨σV⟩ taken from Angulo et al. (1999),

$$\begin{array}{c}\hfill {\mathcal{N}}_{\mathrm{A}}\lambda {(}^{16}\mathrm{O}{\to }^{17}\mathrm{F})=7.37\times {10}^7{e}^{-16.696T_9^{-1/3}}{T}_9^{-0.82},\end{array}$$(6)

where 𝒩_A is Avogadro’s number. With this simple modelling, we can reproduce the evolution of ¹⁴N core abundance as calculated through a more realistic network of nuclear reactions, such as in the MESA code.

In 1D stellar evolution codes, the effect of chemical transport via meridional circulation is typically treated in a diffusive way by adding an extra term to the vertical diffusion coefficient. This additional term takes on the form D_eff = |ru_v|²/(30D_h) (Chaboyer & Zahn 1992), where u_v is the vertical component of the meridional flow, and D_h the horizontal diffusion coefficient (see Eq. (18) in Palacios et al. 2003).

In the present work, we tested the contribution of the meridional circulation to the transport of chemicals using 2D models. Typical values for the maximum of u_v obtained from the 2D models range from 10⁻⁴ to 3 × 10⁻³ cm s⁻¹, which turn out to be too small to have any noticeable influence on element transport. As a result, surface nitrogen abundance for the ESTER model of HD 192575 is weak: Δ[N/H]< 0.05 dex compared to the initial value.

In Appendix E, we show the stream function of the meridional velocity field for a star at the start of the MS, and for a star near the end. As the rotation profile gradually evolves to a more shellular configuration, the angular dependency of the meridional circulation is mostly described by a spherical harmonic of degree ℓ = 2. Therefore, the 2D models show that the assumption that higher order spherical harmonics can be neglected in the expansion of the meridional velocity field is not true in stars at the beginning of the MS, even at the low rotation rate of model M1.

5. Conclusions

In this Letter, we present the first results of a 2D modelling of rotating stars including time evolution. These new models are an update of the 2D ESTER models designed by Espinosa Lara & Rieutord (2013), which compute steady 2D models of fast rotating stars. As in the steady models, the time-dependent ones take into account the centrifugal flattening of the star as well as the large-scale flows (differential rotation and meridional circulation) driven by the baroclinicity of the star.

We ran two 12 M_⊙ models of a massive star with initial angular velocities of 15 and 50% of the critical one, respectively. The first model has a rather mild rotation rate but can be compared to the recent observation of a massive star, while the second model allowed us to test the performance of the code and revealed some new features of the internal dynamics of a massive star (see below).

The first model was actually designed to reproduce the observations of the β Cephei pulsator HD 192575 as derived by Burssens et al. (2023). We found that our model gives a luminosity, effective temperature, and core mass that are consistent with the observationally derived values. Moreover, the rotation profile derived from the ESTER model is also in accordance with the measurements of the core and envelope rotation rates of Burssens et al. (2023). We note that these authors used a simplified rotation profile (e.g., Eq. (4)) to derive the rotation rates. Fortunately, this profile turns out to be similar to the actual one predicted by the models, allowing a consistent comparison.

Another result provided by the above 2D models is the weakness of the meridional circulation. In the dynamics of baroclinic flows, this is controlled by viscosity to ensure the balance of angular momentum flux. The transport of chemicals by this flow seems to be negligible, but this needs to be confirmed by a more detailed analysis, because our models do not include the jump in viscosity expected at the core–envelope boundary, which is expected to drive a Stewartson layer along the tangential cylinder of the core (Rieutord 2006a; Gagnier & Rieutord 2020).

Our 2D models raise new questions as to the dynamics of rotating stars. In particular, the possibility of using a succession of steady-state models to monitor the rotational evolution of early-type stars, as done in Gagnier et al. (2019b), is now questionable and requires fresh investigation. Our results indeed show that for the 12 M_⊙ model we computed, nuclear evolution is only slow enough to relax the star to a quasi-steady state at the beginning of the MS. When the star nears the end of the MS, the nuclear evolution becomes faster than the damping time of baroclinic modes (Busse 1981). The question then arises as to when a succession of steady models is liable to represent the evolution of a rotating star. This is a complex question related to the unsolved general problem of angular momentum transport in stars known since the first measurements of differential rotation in red giant stars (e.g., Beck et al. 2012; Deheuvels et al. 2012, 2015; Mosser et al. 2012). New investigations with the present 2D time-dependent ESTER models will be presented in forthcoming articles.

Acknowledgments

The authors are grateful to the referee prof. Georges Meynet for his comments and suggestions. The research leading to these results has received funding from the French Agence Nationale de la Recherche (ANR), under grant MASSIF (ANR-21-CE31-0018-02). The authors thank Siemen Burssens for providing the MESA setup. Computations of ESTER 2D-models have been possible thanks to HPC resources from CALMIP supercomputing center (Grant 2023-P0107).

References

Appendix A: Stellar equations

In addition to Eq. (1), the ESTER code solves the following set of fundamental stellar equations.

• Mass conservation

  $$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathbf{v})=0,\end{array}$$(A.1)

  where ρ is the density and v the meridional velocity.

• Meridional momentum equation

  $$\begin{array}{c}\hfill \frac{1}{\rho }\mathrm{\nabla }P+\mathrm{\nabla }\varphi -s{\mathrm{\Omega }}^2{\widehat{\mathbf{e}}}_s={\mathbf{F}}_{\mathrm{visc}}^{\mathrm{merid}},\end{array}$$(A.2)

  where P is the pressure, ϕ the gravitational potential, s the radial distance to the rotation axis, Ω the local angular velocity, and ${\mathbf{F}}_{\mathrm{visc}}^{\mathrm{merid}}$ the meridional components of the viscous force. In this equation, the time derivative of the meridional circulation has been neglected because of its extremely small value compared to other terms.

• Angular momentum equation

  $$\begin{array}{c}\hfill \frac{\partial s^2\mathrm{\Omega }}{\partial t}+\mathbf{v}·\mathrm{\nabla }(s^2\mathrm{\Omega })=\frac{1}{\rho }\mathrm{\nabla }(\rho \nu s^2\mathrm{\nabla }\mathrm{\Omega }),\end{array}$$(A.3)

  where ν is the kinematic viscosity. The kinematic viscosities in both the horizontal and vertical directions are free parameters in ESTER. For both, we take 10⁷ cm² s⁻¹, as this was found to be the order of magnitude needed to explain the rotation profiles of F-type stars by Mombarg (2023).

• Entropy equation

  $$\begin{array}{c}\hfill \rho T(\frac{\partial S}{\partial t}+\mathbf{v}·\mathrm{\nabla }S)=\mathrm{\nabla }·(\chi \mathrm{\nabla }T)+\rho \epsilon ,\end{array}$$(A.4)

  where S is the entropy, χ the thermal conductivity, and ε the energy generation rate per unit mass.

• Nuclear energy generation is computed with a simple law given by Kippenhahn & Weigert (1990), namely

  $$\begin{array}{c}\hfill {\epsilon }_{\ast }(\rho ,T_9,X,Z)={\epsilon }_0(X,Z)\rho T_9^{-2/3}exp(-A/T_9^{1/3})(1+C(T_9)),\end{array}$$(A.5)

  from which we deduce the hydrogen consumption Ẋ_nuc. Here, T₉ = T/10⁹ K, A = 15.228, and C(T₉) is a correction term.

Appendix B: Rotation profiles

In this Appendix, we also show the rotation profiles of the M2 model ((Ω_eq/Ω_c)_i = 0.5).

Fig. B.1. Profiles of the angular velocity as a function of the radial coordinate in the polar direction (top panel), equatorial direction (middle panel), and averaged over latitude (bottom panel). The solid lines correspond to model M1, the dashed lines to model M2. The locations of the outer edge of the convective core are shown as vertical lines (dotted lines for model M1, dashed-dotted for model M2).

Appendix C: Timescales

In this Appendix, we show the nuclear and baroclinic timescales (see Sect. 2), once close to the ZAMS, and once close to the TAMS.

Appendix D: MESA setup

In this Appendix, we provide a short summary of the physics used for the 1D MESA model of HD192575 discussed in Section 3. The transport of AM in MESA is treated as a diffusive process (Eq. (B4) in Paxton et al. (2013)) and shellular rotation is imposed. As our aim here is to compare the 2D ESTER models with the 1D physics used by Burssens et al. (2023), we use the same description for the chemical mixing as these authors. This description is based on predictions of simulations of internal gravity waves (e.g. Rogers & McElwaine 2017; Varghese et al. 2023), where the chemical diffusion coefficient takes the form of

$$\begin{array}{c}\hfill D_{\mathrm{IGW}}(r)=D_0\left(\frac{{\rho }_0}{\rho (r)}\right),\end{array}$$(D.1)

where D₀ is a free parameter we set to 10³ cm² s⁻¹ (same as Burssens et al. (2023)), and ρ₀ the density at the core boundary. This means that we set the factor f_C —which accounts for the different efficiencies between the transport of AM and chemical elements⁴ (Heger et al. 2000)— equal to zero.

Appendix E: Meridional circulation

In this Appendix, we show the stream lines of the meridional flow for models M1 ((Ω_eq/Ω_c)_i = 0.15) and model M2 ((Ω_eq/Ω_c)_i = 0.5).

Fig. E.1. Isocontours of the stream function for model M1 of 1 Myr (top panel) and 15.75 Myr (bottom panel). The different colours indicate counter-rotating cells, where the red cells rotate clockwise in the first quadrant of the plot. These models correspond to those shown in Fig. 2. The isobars at the edge of the convective core and at the surface are shown in grey.

Fig. E.2. Same as Fig E.1, but for model M2.

All Tables

All Figures

Fig. 1. HRD showing the evolution tracks of the 2D models for a 12 M⊙ star (Z = 0.012, Xi = 0.71) for (Ωeq/Ωc)i = 0.15 (dashed-dotted), and 0.5 (dashed). The black lines represent the properties at the pole, the grey lines the properties at the equator. The tracks are terminated slightly before the TAMS when the solver can no longer converge, which occurs at Xc/Xini = 0.141 and 0.131 for the respective aforementioned rotation rates. The dots mark the model of HD 192575 that is discussed in Sect. 3.

In the text

	Fig. 2. Map of the angular velocity as a fraction of the critical angular velocity for a model of 1 Myr (top panel) and 15.75 Myr (bottom panel). The plots show cuts in the meridian plane. These plots are for a 12 M⊙ star with (Ωeq/Ωc)i = 0.15.
In the text

Fig. 3. Profile of the vertical diffusion coefficient according to Eq. (2) for a 12 M⊙ model with an initial rotation frequency of 15% the critical angular velocity. The grey dotted lines show $D_{\text{v}}{\langle N_0^2\rangle }_V/(\eta N^2)$ at a few ages (from left to right: Xc/Xini = 0.14, 0.30, 0.58, 0.82, 0.98), that is, if we were to take the local value of the Brunt-Väisälä frequency instead of the volume averaged value.

In the text

Fig. 4. Evolution of the rotation profile of model M1. Left panel: predicted surface rotation as a function of colatitude throughout the evolution (indicated by colour) normalised by the rotation frequency at the pole. Right panel: profile of the angular velocity (averaged over θ) throughout the evolution. The thick coloured line indicates the location of the convective core boundary. The black lines in both panels correspond to the best-matching model for HD 192575 (last column in Table 1).

In the text

Fig. 5. Predicted rotation profile of HD 192575 with ESTER (red solid line) and MESA (black lines). The profiles indicated by a solid line and dotted line are for a constant viscosity of 107 cm2 s−1, while the dashed line is for 109 cm2 s−1. The dots correspond to the location of the core boundary. The red and grey shaded areas indicate the measured core (top one) and surface (bottom one) rotation frequencies by Burssens et al. (2023), assuming rshear equal to the outer boundary of the μ-gradient zone, and rshear = R⋆, respectively.

In the text

Fig. B.1. Profiles of the angular velocity as a function of the radial coordinate in the polar direction (top panel), equatorial direction (middle panel), and averaged over latitude (bottom panel). The solid lines correspond to model M1, the dashed lines to model M2. The locations of the outer edge of the convective core are shown as vertical lines (dotted lines for model M1, dashed-dotted for model M2).

In the text

	Fig. E.1. Isocontours of the stream function for model M1 of 1 Myr (top panel) and 15.75 Myr (bottom panel). The different colours indicate counter-rotating cells, where the red cells rotate clockwise in the first quadrant of the plot. These models correspond to those shown in Fig. 2. The isobars at the edge of the convective core and at the surface are shown in grey.
In the text

	Fig. E.2. Same as Fig E.1, but for model M2.
In the text