# The first two-dimensional stellar structure and evolution models of rotating stars. Calibration to $\beta$ Cephei pulsator HD 192575

J. S. G. Mombarg, M. Rieutord, F. Espinosa Lara

► To cite this version:

J. S. G. Mombarg, M. Rieutord, F. Espinosa Lara. The first two-dimensional stellar structure and evolution models of rotating stars. Calibration to $\beta$ Cephei pulsator HD 192575. Astronomy and Astrophysics - A&A, 2023, 677, 10.1051/0004-6361/202347454. insu-04473160

HAL Id: insu-04473160

https://insu.hal.science/insu-04473160

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER TO THE EDITOR

# The first two-dimensional stellar structure and evolution models of rotating stars

## Calibration to $\beta$ Cephei pulsator HD 192575

J. S. G. Mombarg<fnref n="1" />, M. Rieutord<fnref n="1" />, and F. Espinosa Lara<fnref n="2" />

Received 13 July 2023 / Accepted 14 August 2023

## ABSTRACT

**Contact.** Rotation is a key ingredient in the theory of stellar structure and evolution. Until now, stellar evolution codes operate in a one-dimensional 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 self-consistent 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_{\odot}$ 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 $\beta$ 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 two-dimensional 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

## 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<fnref n="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 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

<footnote n="1">¹ IRAP, Université de Toulouse, CNRS, UPS, CNES, 14 Avenue Édouard Belin, 31400 Toulouse, France
e-mail: jmombarg@irap.omp.eu</footnote>

<footnote n="2">² Space Research Group, University of Alcalá, 28871 Alcalá de Henares, Spain</footnote>

<footnote n="1">¹ https://github.com/ester-project/ester</footnote>
LETTER TO THE EDITOR

# The first two-dimensional stellar structure and evolution models of rotating stars

## Calibration to $\beta$ Cephei pulsator HD 192575

J. S. G. Mombarg<fnref n="1" />, M. Rieutord<fnref n="1" />, and F. Espinosa Lara<fnref n="2" />

Received 13 July 2023 / Accepted 14 August 2023

## ABSTRACT

**Contact.** Rotation is a key ingredient in the theory of stellar structure and evolution. Until now, stellar evolution codes operate in a one-dimensional 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 self-consistent 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_{\odot}$ 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 $\beta$ 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 two-dimensional 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

## 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<fnref n="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 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

<footnote n="1">¹ IRAP, Université de Toulouse, CNRS, UPS, CNES, 14 Avenue Édouard Belin, 31400 Toulouse, France
e-mail: jmombarg@irap.omp.eu</footnote>

<footnote n="2">² Space Research Group, University of Alcalá, 28871 Alcalá de Henares, Spain</footnote>

<footnote n="1">¹ https://github.com/ester-project/ester</footnote>
successfully performed 2D evolution of a slow and a moderate rotator with a mass of 12 $M_\odot$, 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,

$$ \frac{\partial X}{\partial t} + v \cdot \nabla X = \frac{1}{\rho} \nabla(\rho[D]\nabla X) + \dot{X}_{\text{nuc}}, \quad (1) $$

where $\rho$ 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 $\dot{X}_{\text{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,

$$ D_v = \eta \langle N_0^2 \rangle_v^{-1} \langle K r^2 (n \cdot \nabla \Omega)^2 \rangle_\theta. \quad (2) $$

Here, $K$ is the thermal diffusivity, $n$ is the unit vector normal to the isobars, and $\Omega$ 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^2$ that often give rise to numerical difficulties. Furthermore, $\eta$ 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,

<figure><img src="image_1.png" /><figcaption>Fig. 1. HRD showing the evolution tracks of the 2D models for a 12 $M_\odot$ star ($Z = 0.012$, $X_i = 0.71$) for $(\Omega_{eq}/\Omega_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 $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.</figcaption></figure>

$\eta$ needs to be increased to compensate for the smaller shear compared to model M2. The values of $\eta$ are given in Table C.1. Furthermore, we take the angular average of the term between $\langle \cdot \rangle_\theta$ in Eq. (2) and assume a fixed horizontal diffusion coefficient $D_h = 10^5 \text{ cm}^2 \text{ 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_\odot$ star: once for an initial rotation (at the equator) of 15% the critical angular velocity $\Omega_c$ (model M1, $v_{eq} = 112 \text{ km s}^{-1}$) and once for 50% (model M2, $v_{eq} = 358 \text{ km s}^{-1}$)$^2$. 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 $\Omega/\Omega_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

$$ \tau_{\text{baro}} = \left\langle \frac{N^2}{\Omega^2 K} \right\rangle_v \frac{d_{\text{env}}^2}{(\pi^2 + 4)\pi^2}, \quad (3) $$

following the work of Busse (1981). The $\tau_{\text{baro}}$ is usually of the order of the Eddington–Sweet timescale (Rieutord et al. 2006b). In its expression, $d_{\text{env}}$ is the thickness of the radiative envelope in the polar direction. The nuclear evolution timescale,

<footnote>$^2$ The models are available on https://doi.org/10.5281/zenodo.8228904.</footnote>
successfully performed 2D evolution of a slow and a moderate rotator with a mass of 12 $M_\odot$, 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,

$$ \frac{\partial X}{\partial t} + v \cdot \nabla X = \frac{1}{\rho} \nabla(\rho[D]\nabla X) + \dot{X}_{\text{nuc}}, \quad (1) $$

where $\rho$ 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 $\dot{X}_{\text{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,

$$ D_v = \eta \langle N_0^2 \rangle_v^{-1} \langle K r^2 (n \cdot \nabla \Omega)^2 \rangle_\theta. \quad (2) $$

Here, $K$ is the thermal diffusivity, $n$ is the unit vector normal to the isobars, and $\Omega$ is local angular velocity. The term $\langle N_0^2 \rangle_\nu$ 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^2$ that often give rise to numerical difficulties. Furthermore, $\eta$ 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,

<figure><img src="image_2.png" /><figcaption>Fig. 1. HRD showing the evolution tracks of the 2D models for a 12 $M_\odot$ star ($Z = 0.012$, $X_i = 0.71$) for $(\Omega_{\text{eq}}/\Omega_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 $X_c/X_{\text{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.</figcaption></figure>

$\eta$ needs to be increased to compensate for the smaller shear compared to model M2. The values of $\eta$ are given in Table C.1. Furthermore, we take the angular average of the term between $\langle \cdot \rangle_\theta$ in Eq. (2) and assume a fixed horizontal diffusion coefficient $D_h = 10^5 \text{ cm}^2 \text{ 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_\odot$ star: once for an initial rotation (at the equator) of 15% the critical angular velocity $\Omega_c$ (model M1, $v_{\text{eq}} = 112 \text{ km s}^{-1}$) and once for 50% (model M2, $v_{\text{eq}} = 358 \text{ km s}^{-1}$)$^2$. 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 $\Omega/\Omega_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

$$ \tau_{\text{baro}} = \left\langle \frac{N^2}{\Omega^2 K} \right\rangle_{\nu} \frac{d_{\text{env}}^2}{(\pi^2 + 4)\pi^2}, \quad (3) $$

following the work of Busse (1981). The $\tau_{\text{baro}}$ is usually of the order of the Eddington–Sweet timescale (Rieutord et al. 2006b). In its expression, $d_{\text{env}}$ is the thickness of the radiative envelope in the polar direction. The nuclear evolution timescale,

<footnote>$^2$ The models are available on https://doi.org/10.5281/zenodo.8228904.</footnote>
<figure><img src="image_2.png" /><figcaption>**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 $(\Omega_{eq}/\Omega_c)_i = 0.15$.</figcaption></figure>

$\tau_{evol} = X_c/\dot{X}_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 $\tau_{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 $\bar{\Omega}(r)$ decreases from the core towards the surface until roughly 60% of the fractional radius, after which $\bar{\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, $\bar{\Omega}$ 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.

<figure><img src="image_3.png" /><figcaption>**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_v \langle N_0^2 \rangle / (\eta N^2)$ 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 Brunt-Väisälä frequency instead of the volume averaged value.</figcaption></figure>

### 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 $\beta$ 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_\star$), 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

$$ \Omega(r) = \begin{cases} \frac{\Omega_{\text{core}} - \Delta\Omega}{\Omega_{\text{surf}}} & r \le r_{\text{core}}, \\ \frac{r-r_{\text{core}}}{r_{\text{shear}}-r_{\text{core}}} & r_{\text{core}} < r < r_{\text{shear}}, \\ r & r \ge r_{\text{shear}}, \end{cases} \quad (4) $$

where $\Delta\Omega = \Omega_{\text{core}} - \Omega_{\text{surf}}$, $r_{\text{core}}$ is the radius of the convective core, and $r_{\text{shear}}$ is the outer radius of the shear zone where the rotation frequency decreases linearly until it reaches the surface value. The values of $\Omega_{\text{core}}$ and $\Omega_{\text{surf}}$ are then optimised to best reproduce the observations for a given assumption for $r_{\text{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
<figure><img src="image_3.png" /><figcaption>**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 $(\Omega_{eq}/\Omega_c)_i = 0.15$.</figcaption></figure>

$\tau_{evol} = X_c/\dot{X}_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 $\tau_{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 $\bar{\Omega}(r)$ decreases from the core towards the surface until roughly 60% of the fractional radius, after which $\bar{\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, $\bar{\Omega}$ 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.

<figure><img src="image_4.png" /><figcaption>**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_v\langle N_0^2 \rangle_v/(\eta N^2)$ 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 Brunt-Väisälä frequency instead of the volume averaged value.</figcaption></figure>

### 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 $\beta$ 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_\star$), 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

$$ \Omega(r) = \begin{cases} \frac{\Omega_{\text{core}} - \Delta\Omega}{\Omega_{\text{surf}}} & r \le r_{\text{core}}, \\ \frac{r-r_{\text{core}}}{r_{\text{shear}}-r_{\text{core}}} & r_{\text{core}} < r < r_{\text{shear}}, \\ r & r \ge r_{\text{shear}}, \end{cases} \quad (4) $$

where $\Delta\Omega = \Omega_{\text{core}} - \Omega_{\text{surf}}$, $r_{\text{core}}$ is the radius of the convective core, and $r_{\text{shear}}$ is the outer radius of the shear zone where the rotation frequency decreases linearly until it reaches the surface value. The values of $\Omega_{\text{core}}$ and $\Omega_{\text{surf}}$ are then optimised to best reproduce the observations for a given assumption for $r_{\text{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
<table><caption>**Table 1.** Comparison of ESTER with the observations.</caption><thead><tr><th>Parameter</th><th>Burssens et al. (2023)</th><th>2D model</th></tr></thead><tbody><tr><td>M<sub>★</sub></td><td>12 ± 1.5 M<sub>☉</sub></td><td>12 M<sub>☉</sub></td></tr><tr><td>Z</td><td>0.012<sup>+0.004</sup><sub>-0.000</sub></td><td>0.012</td></tr><tr><td>X</td><td>0.71</td><td>0.71</td></tr><tr><td>X<sub>c</sub></td><td>0.176<sup>+0.035</sup><sub>-0.045</sub></td><td>0.168</td></tr><tr><td>f<sub>CBM</sub></td><td>0.030<sup>+0.005</sup><sub>-0.025</sub></td><td>—</td></tr><tr><td>τ</td><td>17<sup>+4.7</sup><sub>-5.4</sub> Myr</td><td>14.75 Myr</td></tr><tr><td>m<sub>cc</sub></td><td>2.9<sup>+0.5</sup><sub>-0.8</sub> M<sub>☉</sub></td><td>2.59 M<sub>☉</sub></td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td>4.30 ± 0.07</td><td>4.36</td></tr><tr><td>T<sub>eff,p</sub></td><td>23 900 ± 900 K</td><td>25 892 K</td></tr><tr><td>T<sub>eff,e</sub></td><td></td><td>25 604 K</td></tr><tr><td>Ω<sub>core</sub>/(2π)</td><td>0.30<sup>+0.09</sup><sub>-0.05</sub> d<sup>-1</sup></td><td>0.40 d<sup>-1</sup></td></tr><tr><td>Ω<sub>surf</sub>/(2π)</td><td>0.20 ± 0.01 d<sup>-1</sup></td><td>0.20 d<sup>-1</sup></td></tr><tr><td>(Ω/Ω<sub>c</sub>)<sub>ini</sub></td><td>—</td><td>0.15</td></tr><tr><td>R<sub>★,p</sub></td><td>9.1<sup>+0.8</sup><sub>-1.7</sub> R<sub>☉</sub></td><td>7.57 R<sub>☉</sub></td></tr><tr><td>R<sub>★,e</sub></td><td></td><td>7.64 R<sub>☉</sub></td></tr></tbody></table>

**Notes.** Stellar parameters of HD 192575 derived by Burssens et al. (2023) and those of the ESTER model for this same star. The parameters above the horizontal line are input for the model, while those below are output.

axis<fnref n="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 quasi-steady state.

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 ($\mu$). 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 $\Omega_{core}$ in the ESTER models is defined as the value of $\Omega(r_{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_{\star}$ (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 $(\Omega_{eq}/\Omega_{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} \text{ cm}^{2} \text{ s}^{-1}$ compared to the value of $10^{7} \text{ cm}^{2} \text{ 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).

## 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 pp-chain and CNO cycle (Kippenhahn & Weigert 1990; Rieutord et al. 2016). To track the abundance of $^{14}$N as a result of hydrogen-burning 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

$$ \rho \frac{dX(^{14}\text{N})}{dt} = \lambda(^{16}\text{O} \rightarrow ^{17}\text{F}) n(^{16}\text{O}) n(\text{H}) m(^{14}\text{N}), \quad (5) $$

where $n(\cdot)$ denotes the number densities, $m$ the atomic mass, and $\lambda(^{16}\text{O} \rightarrow ^{17}\text{F})$ is the Maxwellian-average reaction rate $\langle\sigma V\rangle$ taken from Angulo et al. (1999),

$$ N_A \lambda(^{16}\text{O} \rightarrow ^{17}\text{F}) = 7.37 \times 10^7 e^{-16.696 T_9^{-1/3}} T_9^{-0.82}, \quad (6) $$

where $N_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 \times 10^{-3} \text{ 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: $\Delta[N/H] < 0.05 \text{ 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 stellar configuration, the angular dependency of the meridional circulation is mostly described by a spherical harmonic of

<footnote n="3">³ 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 $\rho v$ is used instead of $v$.</footnote>
<table><caption>**Table 1.** Comparison of ESTER with the observations.</caption><thead><tr><th>Parameter</th><th>Burssens et al. (2023)</th><th>2D model</th></tr></thead><tbody><tr><td>M<sub>★</sub></td><td>12 ± 1.5 M<sub>☉</sub></td><td>12 M<sub>☉</sub></td></tr><tr><td>Z</td><td>0.012<sup>+0.004</sup><sub>-0.000</sub></td><td>0.012</td></tr><tr><td>X</td><td>0.71</td><td>0.71</td></tr><tr><td>X<sub>c</sub></td><td>0.176<sup>+0.035</sup><sub>-0.045</sub></td><td>0.168</td></tr><tr><td>f<sub>CBM</sub></td><td>0.030<sup>+0.005</sup><sub>-0.025</sub></td><td>—</td></tr><tr><td>τ</td><td>17<sup>+4.7</sup><sub>-5.4</sub> Myr</td><td>14.75 Myr</td></tr><tr><td>m<sub>cc</sub></td><td>2.9<sup>+0.5</sup><sub>-0.8</sub> M<sub>☉</sub></td><td>2.59 M<sub>☉</sub></td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td>4.30 ± 0.07</td><td>4.36</td></tr><tr><td>T<sub>eff,p</sub></td><td>23 900 ± 900 K</td><td>25 892 K</td></tr><tr><td>T<sub>eff,e</sub></td><td></td><td>25 604 K</td></tr><tr><td>Ω<sub>core</sub>/(2π)</td><td>0.30<sup>+0.09</sup><sub>-0.05</sub> d<sup>-1</sup></td><td>0.40 d<sup>-1</sup></td></tr><tr><td>Ω<sub>surf</sub>/(2π)</td><td>0.20 ± 0.01 d<sup>-1</sup></td><td>0.20 d<sup>-1</sup></td></tr><tr><td>(Ω/Ω<sub>c</sub>)<sub>ini</sub></td><td>—</td><td>0.15</td></tr><tr><td>R<sub>★,p</sub></td><td>9.1<sup>+0.8</sup><sub>-1.7</sub> R<sub>☉</sub></td><td>7.57 R<sub>☉</sub></td></tr><tr><td>R<sub>★,e</sub></td><td></td><td>7.64 R<sub>☉</sub></td></tr></tbody></table>

**Notes.** Stellar parameters of HD 192575 derived by Burssens et al. (2023) and those of the ESTER model for this same star. The parameters above the horizontal line are input for the model, while those below are output.

axis<fnref n="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 quasi-steady state.

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 ($\mu$). 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 $\Omega_{core}$ in the ESTER models is defined as the value of $\Omega(r_{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_{\star}$ (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 $(\Omega_{eq}/\Omega_{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} \text{ cm}^{2} \text{ s}^{-1}$ compared to the value of $10^{7} \text{ cm}^{2} \text{ 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).

## 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 pp-chain and CNO cycle (Kippenhahn & Weigert 1990; Rieutord et al. 2016). To track the abundance of $^{14}$N as a result of hydrogen-burning 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

$$ \rho \frac{dX(^{14}\text{N})}{dt} = \lambda(^{16}\text{O} \rightarrow ^{17}\text{F}) n(^{16}\text{O}) n(\text{H}) m(^{14}\text{N}), \quad (5) $$

where $n(\cdot)$ denotes the number densities, $m$ the atomic mass, and $\lambda(^{16}\text{O} \rightarrow ^{17}\text{F})$ is the Maxwellian-average reaction rate $\langle \sigma V \rangle$ taken from Angulo et al. (1999),

$$ N_A \lambda(^{16}\text{O} \rightarrow ^{17}\text{F}) = 7.37 \times 10^7 e^{-16.696 T_9^{-1/3}} T_9^{-0.82}, \quad (6) $$

where $N_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 \times 10^{-3} \text{ 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: $\Delta[N/H] < 0.05 \text{ 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 stellar configuration, the angular dependency of the meridional circulation is mostly described by a spherical harmonic of

<footnote n="3">³ 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 $\rho v$ is used instead of $v$.</footnote>
<figure><img src="image_4.png" /></figure>

**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 $\theta$) 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><img src="image_5.png" /></figure>

**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 $\mu$-gradient zone, and $r_{shear} = R_*$, respectively.

degree $\ell = 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_\odot$ 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 $\beta$ 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_\odot$ 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.
<figure><img src="image_5.png" /></figure>

**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 $\theta$) 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><img src="image_6.png" /></figure>

**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 $\mu$-gradient zone, and $r_{shear} = R_*$, respectively.

degree $\ell = 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_\odot$ 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 $\beta$ 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_\odot$ 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.
*Acknowledgements.* 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

Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3

Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55

Bouchaud, K., Domiciano de Souza, A., Rieutord, M., Reese, D. R., & Kervella, P. 2020, A&A, 633, A78

Brott, I., Evans, C. J., Hunter, I., et al. 2011, A&A, 530, A116

Burssens, S., Bowman, D. M., Michielsen, M., et al. 2023, Nat. Astron., 3, 760

Busse, F. 1981, Geophys. Astrophys. Fluid Dyn., 17, 215

Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173

Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19

Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96

Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146

Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35

Gagnier, D., & Rieutord, M. 2020, J. Fluid Mech., 904, A35

Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019a, A&A, 625, A88

Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019b, A&A, 625, A89

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368

Howarth, I. D., Bailey, J., Cotton, D. V., & Kedziora-Chudczer, L. 2023, MNRAS, 520, 1193

Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15

Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution (New York: Springer)

Mathis, S., Palacios, A., & Zahn, J. P. 2004, A&A, 425, 243

Mombarg, J. S. G. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202345956

Mombarg, J. S. G., Dotter, A., Rieutord, M., et al. 2022, ApJ, 925, 154

Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10

Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603

Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3

Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4

Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15

Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34

Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10

Rieutord, M. 2006a, A&A, 451, 1025

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

Rieutord, M. 2015, Fluid Dynamics: An Introduction (New York: Springer), 508

Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277

Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1

Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902

Talon, S., & Charbonnel, C. 1998, A&A, 335, 959

Vanlaer, V., Aerts, C., Bellinger, E. P., & Christensen-Dalsgaard, J. 2023, A&A, 675, A17

Varghese, A., Ratnasingam, R. P., Vanon, R., Edelmann, P. V. F., & Rogers, T. M. 2023, ApJ, 942, 53

Zahn, J. P. 1992, A&A, 265, 115
*Acknowledgements.* 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

Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3

Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55

Bouchaud, K., Domiciano de Souza, A., Rieutord, M., Reese, D. R., & Kervella, P. 2020, A&A, 633, A78

Brott, I., Evans, C. J., Hunter, I., et al. 2011, A&A, 530, A116

Burssens, S., Bowman, D. M., Michielsen, M., et al. 2023, Nat. Astron., 3, 760

Busse, F. 1981, Geophys. Astrophys. Fluid Dyn., 17, 215

Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173

Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19

Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96

Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146

Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35

Gagnier, D., & Rieutord, M. 2020, J. Fluid Mech., 904, A35

Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019a, A&A, 625, A88

Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019b, A&A, 625, A89

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368

Howarth, I. D., Bailey, J., Cotton, D. V., & Kedziora-Chudczer, L. 2023, MNRAS, 520, 1193

Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15

Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution (New York: Springer)

Mathis, S., Palacios, A., & Zahn, J. P. 2004, A&A, 425, 243

Mombarg, J. S. G. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202345956

Mombarg, J. S. G., Dotter, A., Rieutord, M., et al. 2022, ApJ, 925, 154

Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10

Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603

Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3

Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4

Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15

Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34

Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10

Rieutord, M. 2006a, A&A, 451, 1025

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

Rieutord, M. 2015, Fluid Dynamics: An Introduction (New York: Springer), 508

Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277

Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1

Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902

Talon, S., & Charbonnel, C. 1998, A&A, 335, 959

Vanlaer, V., Aerts, C., Bellinger, E. P., & Christensen-Dalsgaard, J. 2023, A&A, 675, A17

Varghese, A., Ratnasingam, R. P., Vanon, R., Edelmann, P. V. F., & Rogers, T. M. 2023, ApJ, 942, 53

Zahn, J. P. 1992, A&A, 265, 115
# Appendix A: Stellar equations

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

- Mass conservation

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (A.1) $$

where $\rho$ is the density and $\mathbf{v}$ the meridional velocity.

- Meridional momentum equation

$$ \frac{1}{\rho} \nabla P + \nabla \phi - s \Omega^2 \hat{\mathbf{e}}_s = \mathbf{F}_{\text{visc}}^{\text{merid}}, \qquad (A.2) $$

where $P$ is the pressure, $\phi$ the gravitational potential, $s$ the radial distance to the rotation axis, $\Omega$ the local angular velocity, and $\mathbf{F}_{\text{visc}}^{\text{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

$$ \frac{\partial s^2 \Omega}{\partial t} + \mathbf{v} \cdot \nabla (s^2 \Omega) = \frac{1}{\rho} \nabla (\rho v s^2 \nabla \Omega), \qquad (A.3) $$

where $\mathbf{v}$ 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 \text{ cm}^2 \text{ s}^{-1}$, 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

$$ \rho T \left( \frac{\partial S}{\partial t} + \mathbf{v} \cdot \nabla S \right) = \nabla \cdot (\chi \nabla T) + \rho \varepsilon, \qquad (A.4) $$

where $S$ is the entropy, $\chi$ the thermal conductivity, and $\varepsilon$ the energy generation rate per unit mass.

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

$$ \varepsilon_*(\rho, T_9, X, Z) = \varepsilon_0(X, Z) \rho T_9^{-2/3} \exp(-A/T_9^{1/3}) (1 + C(T_9)), \qquad (A.5) $$

from which we deduce the hydrogen consumption $\dot{X}_{\text{nuc}}$. Here, $T_9 = T/10^9 \text{ 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 (($\Omega_{\text{eq}}/\Omega_c$)$_i = 0.5$).

<figure><img src="image_6.png" /><figcaption>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).</figcaption></figure>
# Appendix A: Stellar equations

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

- Mass conservation

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (A.1) $$

where $\rho$ is the density and $\mathbf{v}$ the meridional velocity.

- Meridional momentum equation

$$ \frac{1}{\rho} \nabla P + \nabla \phi - s \Omega^2 \hat{\mathbf{e}}_s = \mathbf{F}_{\text{visc}}^{\text{merid}}, \qquad (A.2) $$

where $P$ is the pressure, $\phi$ the gravitational potential, $s$ the radial distance to the rotation axis, $\Omega$ the local angular velocity, and $\mathbf{F}_{\text{visc}}^{\text{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

$$ \frac{\partial s^2 \Omega}{\partial t} + \mathbf{v} \cdot \nabla (s^2 \Omega) = \frac{1}{\rho} \nabla (\rho v s^2 \nabla \Omega), \qquad (A.3) $$

where $\mathbf{v}$ 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 \text{ cm}^2 \text{ s}^{-1}$, 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

$$ \rho T \left( \frac{\partial S}{\partial t} + \mathbf{v} \cdot \nabla S \right) = \nabla \cdot (\chi \nabla T) + \rho \varepsilon, \qquad (A.4) $$

where $S$ is the entropy, $\chi$ the thermal conductivity, and $\varepsilon$ the energy generation rate per unit mass.

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

$$ \varepsilon_*(\rho, T_9, X, Z) = \varepsilon_0(X, Z) \rho T_9^{-2/3} \exp(-A/T_9^{1/3}) (1 + C(T_9)), \qquad (A.5) $$

from which we deduce the hydrogen consumption $\dot{X}_{\text{nuc}}$. Here, $T_9 = T/10^9 \text{ 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 (($\Omega_{\text{eq}}/\Omega_c$)$_i = 0.5$).

<figure><img src="image_7.png" /><figcaption>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).</figcaption></figure>
## 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.

<table><caption>**Table C.1.** Diffusivity factor and characteristic timescales of the models.</caption><thead><tr><th>Model</th><th>(Ω/Ω<sub>c</sub>)<sub>i</sub></th><th>η</th><th>Age</th><th>τ<sub>nuc</sub><br>(Myrs)</th><th>τ<sub>baro</sub><br>(Myrs)</th></tr></thead><tbody><tr><td rowspan="2">M1</td><td rowspan="2">0.15</td><td rowspan="2">9.48</td><td>ZAMS</td><td>2.72</td><td>0.95</td></tr><tr><td>TAMS</td><td>0.12</td><td>7.55</td></tr><tr><td rowspan="2">M2</td><td rowspan="2">0.50</td><td rowspan="2">0.72</td><td>ZAMS</td><td>2.80</td><td>0.09</td></tr><tr><td>TAMS</td><td>0.10</td><td>0.80</td></tr></tbody></table>

## 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 heliostatic 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

$$D_{\text{IGW}}(r) = D_0 \left( \frac{\rho_0}{\rho(r)} \right), \qquad (\text{D.1})$$

where $D_0$ is a free parameter we set to $10^3 \text{ cm}^2 \text{ s}^{-1}$ (same as Burssens et al. (2023)), and $\rho_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<fnref n="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).

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

<footnote n="4">⁴ am_D_mix_factor in MESA.</footnote>
## 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.

<table><caption>**Table C.1.** Diffusivity factor and characteristic timescales of the models.</caption><thead><tr><th>Model</th><th>(Ω/Ω<sub>c</sub>)<sub>i</sub></th><th>η</th><th>Age</th><th>τ<sub>nuc</sub><br>(Myrs)</th><th>τ<sub>baro</sub><br>(Myrs)</th></tr></thead><tbody><tr><td rowspan="2">M1</td><td rowspan="2">0.15</td><td rowspan="2">9.48</td><td>ZAMS</td><td>2.72</td><td>0.95</td></tr><tr><td>TAMS</td><td>0.12</td><td>7.55</td></tr><tr><td rowspan="2">M2</td><td rowspan="2">0.50</td><td rowspan="2">0.72</td><td>ZAMS</td><td>2.80</td><td>0.09</td></tr><tr><td>TAMS</td><td>0.10</td><td>0.80</td></tr></tbody></table>

## 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 heliostatic 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

$$D_{\text{IGW}}(r) = D_0 \left( \frac{\rho_0}{\rho(r)} \right), \qquad (\text{D.1})$$

where $D_0$ is a free parameter we set to $10^3 \text{ cm}^2 \text{ s}^{-1}$ (same as Burssens et al. (2023)), and $\rho_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<fnref n="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).

<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

<footnote n="4">⁴ am_D_mix_factor in MESA.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. E.2. Same as Fig E.1, but for model M2.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. E.2. Same as Fig E.1, but for model M2.</figcaption></figure>