Issue 
A&A
Volume 677, September 2023



Article Number  L5  
Number of page(s)  9  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/202347454  
Published online  31 August 2023 
Letter to the Editor
The first twodimensional stellar structure and evolution models of rotating stars
Calibration to β Cephei pulsator HD 192575
^{1}
IRAP, Université de Toulouse, CNRS, UPS, CNES, 14 Avenue Édouard Belin, 31400 Toulouse, France
email: jmombarg@irap.omp.eu
^{2}
Space Research Group, University of Alcalá, 28871 Alcalá de Henares, Spain
Received:
13
July
2023
Accepted:
14
August
2023
Contact. Rotation is a key ingredient in the theory of stellar structure and evolution. Until now, stellar evolution codes operate in a onedimensional framework for which the validity domain in regards to the rotation rate is not well understood.
Aims. In this Letter, we present the first results of selfconsistent stellar models in two spatial dimensions, which compute the time evolution of a star and its rotation rate along the main sequence (MS). We also present a comparison to observations.
Methods. We make use of an extended version of the ESTER code, which solves the stellar structure of a rotating star in two dimensions with time evolution, including chemical evolution, and an implementation of rotational mixing. We computed evolution tracks for a 12 M_{⊙} model, once for an initial rotation rate equal to 15% of the critical frequency, and once for 50%.
Results. We first show that our model initially rotating at 15% of the critical frequency is able to reproduce all the observations of the β Cephei star HD 192575, which was recently studied with asteroseismology. Beyond the classical surface parameters, such as effective temperature or luminosity, our model also reproduces the core mass along with the rotation rate of the core and envelope at the estimated age of the star. This particular model also shows that the meridional circulation has a negligible influence on the transport of chemical elements such as nitrogen, for which the abundance may be increased at the stellar surface. Furthermore, it shows that in the late MS, nuclear evolution is faster than the relaxation time needed to reach a steady state of stellar angular momentum distribution.
Conclusions. We demonstrate that we have successfully taken a new step towards twodimensional evolutionary modelling of rotating stars. This opens new perspectives on the understanding of the dynamics of fast rotating stars and on the way rotation impacts stellar evolution.
Key words: asteroseismology / stars: interiors / stars: massive / stars: rotation / stars: evolution
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Over the next decades, stellar structure and evolution (SSE) theory will begin to advance from onedimensional 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 sheardriven 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 selfconsistent manner in two dimensions, that is, with the ESTER^{1} 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 earlytype 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 twodimensional steadystate models by solving, in addition, the equation of chemical evolution. The overall aim of the present work is to demonstrate the capabilities of the stateoftheart 2D evolution ESTER models of a rotating massive star, and present predictions for the rotation and chemical profiles. This includes fully solving the largescale 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 zeroage main sequence (ZAMS) up to near the terminalage 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 largescale flows, namely the differential rotation and meridional circulations. The first version of the ESTER code, which computes the steady states of earlytype 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,
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 steadystate 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 firstorder 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,
Here, K is the thermal diffusivity, n is the unit vector normal to the isobars, and Ω is local angular velocity. The term is the volumeaveraged squared Brunt–Väisälä frequency of the initial steadystate model, where the average is taken over the volume where . We therefore smooth out the rapid radial variations of N^{2} that often give rise to numerical difficulties. Furthermore, η is a free nondimensional 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^{5} cm^{2} s^{−1}. 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^{−1}) and once for 50% (model M2, v_{eq} = 358 km s^{−1})^{2}. The evolution models start from a steadystate model at the ZAMS with a uniform chemical composition with X = 0.71 and Z = 0.012 (no assumption of solidbody rotation). The evolution in the HertzsprungRussell 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 solarlike, and gradually evolve towards a more shellular rotation with a slightly antisolar surface rotation. The same behaviour is also seen in M2 (see Appendix B). Contrary to an evolution modelled with a series of steadystate 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
Fig. 1. HRD showing the evolution tracks of the 2D models for a 12 M_{⊙} star (Z = 0.012, X_{i} = 0.71) for (Ω_{eq}/Ω_{c})_{i} = 0.15 (dasheddotted), 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 X_{c}/X_{ini} = 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 decreases from the core towards the surface until roughly 60% of the fractional radius, after which 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, keeps decreasing all the way towards the surface. Figure 3 also shows the difference between taking an average value for N^{2} 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 at a few ages (from left to right: X_{c}/X_{ini} = 0.14, 0.30, 0.58, 0.82, 0.98), that is, if we were to take the local value of the BruntVä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 nonrotating 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 hydrogenmass fraction (X) fixed to the values found by B23. We then picked the model for which, after evolution, the hydrogenmass 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
Comparison of ESTER with the observations.
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 selfconsistently. 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 solarlike to antisolar. This change may be interpreted as follows. The solarlike 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^{3}, 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 quasisteady 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 bestmatching 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 piecewise profile given in Eq. (4) when r_{shear} is taken equal to the outer radius of the region with a nonzero 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 , 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 coreboundary 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^{9} cm^{2} s^{−1} compared to the value of 10^{7} cm^{2} s^{−1} 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 10^{7} cm^{2} s^{−1}, while the dashed line is for 10^{9} cm^{2} 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 r_{shear} equal to the outer boundary of the μgradient zone, and r_{shear} = R_{⋆}, respectively. 
4. Nitrogen enhancement at the surface
The observed abundance of ^{14}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 ppchain and CNO cycle (Kippenhahn & Weigert 1990; Rieutord et al. 2016). To track the abundance of ^{14}N as a result of hydrogenburning via the CNO cycle, we make the following assumptions. First, at the start of the evolution, all ^{12}C present in the core is converted to ^{14}N. The reaction that controls the rate at which ^{14}N is generated is therefore the proton capture of ^{16}O to create ^{17}F, which then rapidly decays into ^{17}O, itself reacting on protons to yield ^{14}N and ^{4}He. Therefore, the evolution of the mass fraction of ^{14}N due to nuclear reactions is described by
where n(⋅) denotes the number densities, m the atomic mass, and λ(^{16}O→^{17}F) is the Maxwellianaverage reaction rate ⟨σV⟩ taken from Angulo et al. (1999),
where 𝒩_{A} is Avogadro’s number. With this simple modelling, we can reproduce the evolution of ^{14}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}^{2}/(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^{−4} to 3 × 10^{−3} cm s^{−1}, 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 timedependent ones take into account the centrifugal flattening of the star as well as the largescale 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 steadystate models to monitor the rotational evolution of earlytype 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 quasisteady 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 timedependent ESTER models will be presented in forthcoming articles.
The models are available on https://doi.org/10.5281/zenodo.8228904.
This is exactly true for a steady solution of an incompressible inviscid rotating fluid when the Coriolis term dominates all other terms (Rieutord 2015), but it can be extended to fluids of varying density if momentum ρv is used instead of v.
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 (ANR21CE31001802). The authors thank Siemen Burssens for providing the MESA setup. Computations of ESTER 2Dmodels have been possible thanks to HPC resources from CALMIP supercomputing center (Grant 2023P0107).
References
 Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [Google Scholar]
 Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55 [Google Scholar]
 Bouchaud, K., Domiciano de Souza, A., Rieutord, M., Reese, D. R., & Kervella, P. 2020, A&A, 633, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brott, I., Evans, C. J., Hunter, I., et al. 2011, A&A, 530, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Burssens, S., Bowman, D. M., Michielsen, M., et al. 2023, Nat. Astron., 3, 760 [NASA ADS] [Google Scholar]
 Busse, F. 1981, Geophys. Astrophys. Fluid Dyn., 17, 215 [Google Scholar]
 Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173 [NASA ADS] [Google Scholar]
 Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19 [Google Scholar]
 Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gagnier, D., & Rieutord, M. 2020, J. Fluid Mech., 904, A35 [NASA ADS] [CrossRef] [Google Scholar]
 Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019a, A&A, 625, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019b, A&A, 625, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Howarth, I. D., Bailey, J., Cotton, D. V., & KedzioraChudczer, L. 2023, MNRAS, 520, 1193 [NASA ADS] [CrossRef] [Google Scholar]
 Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution (New York: Springer) [Google Scholar]
 Mathis, S., Palacios, A., & Zahn, J. P. 2004, A&A, 425, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mombarg, J. S. G. 2023, A&A, in press, https://doi.org/10.1051/00046361/202345956 [Google Scholar]
 Mombarg, J. S. G., Dotter, A., Rieutord, M., et al. 2022, ApJ, 925, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603 [CrossRef] [EDP Sciences] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Rieutord, M. 2006a, A&A, 451, 1025 [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2006b, in Stellar Fluid Dynamics and Numerical Simulations: From the Sun to Neutron Stars, eds. M. Rieutord, & B. Dubrulle, EAS Publ., 21, 275 [NASA ADS] [Google Scholar]
 Rieutord, M. 2015, Fluid Dynamics: An Introduction (New York: Springer), 508 [Google Scholar]
 Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1 [CrossRef] [Google Scholar]
 Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902 [Google Scholar]
 Talon, S., & Charbonnel, C. 1998, A&A, 335, 959 [NASA ADS] [Google Scholar]
 Vanlaer, V., Aerts, C., Bellinger, E. P., & ChristensenDalsgaard, J. 2023, A&A, 675, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Varghese, A., Ratnasingam, R. P., Vanon, R., Edelmann, P. V. F., & Rogers, T. M. 2023, ApJ, 942, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J. P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
Appendix A: Stellar equations
In addition to Eq. (1), the ESTER code solves the following set of fundamental stellar equations.

Mass conservation
where ρ is the density and v the meridional velocity.

Meridional momentum equation
where P is the pressure, ϕ the gravitational potential, s the radial distance to the rotation axis, Ω the local angular velocity, and 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
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^{7} cm^{2} s^{−1}, as this was found to be the order of magnitude needed to explain the rotation profiles of Ftype stars by Mombarg (2023).

Entropy equation
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
from which we deduce the hydrogen consumption Ẋ_{nuc}. Here, T_{9} = T/10^{9} K, A = 15.228, and C(T_{9}) 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, dasheddotted 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.
Diffusivity factor and characteristic timescales of the models.
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
where D_{0} is a free parameter we set to 10^{3} cm^{2} s^{−1} (same as Burssens et al. (2023)), and ρ_{0} 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^{4} (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 counterrotating 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. 
All Tables
All Figures
Fig. 1. HRD showing the evolution tracks of the 2D models for a 12 M_{⊙} star (Z = 0.012, X_{i} = 0.71) for (Ω_{eq}/Ω_{c})_{i} = 0.15 (dasheddotted), 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 X_{c}/X_{ini} = 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 at a few ages (from left to right: X_{c}/X_{ini} = 0.14, 0.30, 0.58, 0.82, 0.98), that is, if we were to take the local value of the BruntVä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 bestmatching 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 10^{7} cm^{2} s^{−1}, while the dashed line is for 10^{9} cm^{2} 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 r_{shear} equal to the outer boundary of the μgradient zone, and r_{shear} = 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, dasheddotted 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 counterrotating 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.