ESTER: 2D Stellar Modeling for Fast Rotators

• ESTER is a two-dimensional stellar-structure framework designed for modeling rapidly rotating stars, accounting for centrifugal distortion, gravity darkening, and differential rotation.
• It employs a multidomain spectral method to solve coupled partial differential equations that self-consistently capture internal dynamics and observable surface quantities.
• Recent extensions enable time-dependent evolution, chemical mixing, and asteroseismic applications, bridging models to interferometric and spectroscopic observations.

Evolution STEllaire en Rotation (ESTER) is a two-dimensional stellar-structure framework for rotating stars that computes the self-consistent structure and large-scale flows of centrifugally deformed, baroclinic stellar models. In its original form, ESTER is a steady-state, axisymmetric, non-magnetic solver for early-type stars with convective cores and radiative envelopes; later developments extended it to time-dependent two-dimensional stellar evolution with chemical evolution and rotational mixing. Its scientific role is to replace spherical or perturbative descriptions when rapid rotation makes oblateness, gravity darkening, meridional circulation, and differential rotation observationally and structurally important (Rieutord et al., 2016, Rieutord et al., 2013, Mombarg et al., 2023).

1. Origin, scope, and scientific motivation

Fast rotation breaks spherical symmetry through centrifugal distortion, producing an equatorial radius larger than the polar radius and a latitude-dependent effective gravity. In radiative envelopes this leads to gravity darkening, with hotter, brighter poles and cooler, dimmer equatorial regions. Rotation also renders the stellar interior baroclinic, so isobars and isotherms are misaligned; the resulting baroclinic torques drive meridional circulation and differential rotation even without magnetic fields. These are intrinsically multidimensional phenomena, which is why one-dimensional rotating models and perturbative treatments lose validity for sufficiently rapid rotators (Bouchaud et al., 2019, Rieutord et al., 2013, Rieutord et al., 2013).

ESTER was developed specifically for isolated, rapidly rotating early-type stars with radiative envelopes. In the steady-state formulation it targets main-sequence objects and treats the convective core as isentropic. This domain makes it directly relevant to A- and B-type stars observed with optical or infrared interferometry, spectroscopy, photometry, polarimetry, and asteroseismology. Several papers describe ESTER as the first code able to compute self-consistent two-dimensional stellar models of fast rotators including their large-scale flows, and later work emphasizes that it is open source and freely available to the community (Rieutord et al., 2016, White et al., 24 Sep 2025).

The framework emerged in response to concrete observational demands. Interferometry resolves oblateness and gravity darkening; spectroscopic line profiles depend on the latitude dependence of temperature and gravity; pulsation spectra depend on the full deformed structure and rotation profile; and, in some cases, precise polarimetry is sensitive to the asymmetric surface radiation field of near-critical rotators. ESTER was designed to provide a common physical model linking these observables to stellar structure and internal dynamics (Bouchaud et al., 2019, Howarth et al., 2023).

2. Governing equations and physical content

In its steady-state form, ESTER solves four coupled partial differential equations in two spatial dimensions for a self-gravitating rotating star. These are Poisson’s equation for the gravitational potential,

$\nabla^2 \Phi = 4\pi G\,\rho,$

the momentum equation,

$$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$

the entropy or energy equation,

$$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$

and mass conservation,

$\nabla\cdot(\rho\,\mathbf{v}) = 0.$

In radiative regions the heat flux is represented by radiative diffusion, for example

$$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$

while the convective core is assumed perfectly isentropic. In barotropic limits or interpretive reductions, hydrostatic support can be written with an effective gravity or effective potential, but ESTER is fundamentally built for baroclinic configurations in which the scalar reduction is not exact (Rieutord et al., 2016, Bouchaud et al., 2019).

A central concept is the thermal-wind relation linking baroclinicity to differential rotation. In the small-Ekman-number limit, the meridional Euler balance and the associated angular-momentum transport equation determine the two-dimensional rotation profile $\Omega(r,\theta)$ together with the meridional circulation. Viscosity is required to lift the inviscid degeneracy of the thermal-wind balance and to determine a unique differential rotation profile. ESTER therefore does not impose solid-body rotation except in special cases; it derives the internal rotation law self-consistently from structure, thermal stratification, and transport (Rieutord et al., 2016, Rieutord et al., 2013).

The microphysics used in the published implementations include OPAL equations of state and opacities, together with nuclear energy generation from pp-chain and CNO reactions. Earlier descriptions emphasize NACRE reaction rates; later formulations also use analytic pp and CNO prescriptions or simplified networks depending on the application. The outer boundary is an isobar chosen to approximate the photosphere. Boundary conditions impose vacuum gravity at large radius, stress-free flow at the stellar surface, and a radiative boundary condition equivalent to blackbody emission. In the educational guide, the visible photosphere is obtained by extrapolating from the bounding surface with an $n=3$ polytropic atmosphere, yielding latitude-dependent effective temperatures and fluxes without imposing a single gravity-darkening exponent (Rieutord et al., 2016, Rieutord et al., 2013, White et al., 24 Sep 2025).

Gravity darkening is therefore a consequence of the two-dimensional radiative structure rather than a free power law. For comparison with older literature, one may write

$$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$

but ESTER-based studies repeatedly stress that the latitudinal flux distribution follows from the full two-dimensional solution. In practical applications, the framework is often compared with the $\omega$-model gravity-darkening prescription of Espinosa Lara and Rieutord, which assumes radiative flux nearly antiparallel to the effective gravity and is especially useful in Roche-geometry forward models (Bouchaud et al., 2019, Howarth et al., 2023, Rieutord et al., 2013).

3. Numerical formulation and solver architecture

ESTER uses a multidomain spectral method adapted to a centrifugally distorted stellar geometry. The star is decomposed into spheroidal shells bounded by isobars, and a nonlinear coordinate mapping aligns the computational domains with those surfaces while reducing to spherical coordinates near the center. Radially, the solution is represented with Chebyshev polynomials and Gauss–Lobatto collocation points; latitudinally, axisymmetric fields are expanded in spherical harmonics. An external compactified domain is added for the gravitational potential so that the vacuum condition at infinity can be imposed consistently (Rieutord et al., 2016).

This geometry-aware mapping is one of the defining numerical features of ESTER. It allows the free surface of the oblate star to be solved as part of the problem rather than prescribed in advance. Domain interfaces are arranged so that strong radial stratification is distributed across shells with similar pressure ratios, which improves convergence and keeps spectral accuracy acceptable even when realistic tabulated microphysics reduces smoothness (Rieutord et al., 2016).

The nonlinear solution strategy depends on the physics. Picard iteration works efficiently for simple rotating polytropes and can achieve virial accuracies of order $\epsilon_2 \sim 10^{-9}$, but realistic stellar models with tabulated opacities and equations of state require Newton’s method. In that regime ESTER solves a large Jacobian system of the form

$$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$0

including variations of both the physical fields and the spheroidal mapping. The linear algebra combines LU factorization of the block-tridiagonal domain structure with Conjugate Gradient Squared iterations used as a preconditioned Krylov step. Reported accuracies for realistic models include virial errors of at most a few $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$1 and global energy-balance errors around $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$2–$$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$3, with the opacity tables being a principal limiter of spectral convergence. With increased angular resolution and domain count, the code reaches rotations around $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$4–$$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$5, where the equatorial cusp becomes the dominant numerical difficulty (Rieutord et al., 2016).

Later documentation of version 1.1.0rc2 describes a practical steady-state workflow in which one first computes a one-dimensional non-rotating or slowly rotating model, then uses that solution to initialize a two-dimensional rotating calculation, and finally walks through parameter space by small changes in rotation or composition until convergence is achieved. That guide gives representative resolutions such as 15 domains, 60 radial points per domain, 24 latitudinal points, and 30 points in the external domain, and notes that Newton–Raphson convergence can fail when parameter jumps are too large (White et al., 24 Sep 2025).

4. Coupling to atmospheres and observables

A major strength of ESTER is that it delivers surface quantities as latitude-dependent functions: radius $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$6, effective temperature $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$7, surface gravity $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$8 or $$\rho\,\mathbf{v}\cdot\nabla\mathbf{v} = -\nabla P - \rho\,\nabla\Phi + \mathbf{F}_v,$$9, and rotation $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$0. These outputs can then be coupled to atmosphere models to compute interferometric images, spectroscopic line profiles, photometric fluxes, and, in suitable cases, polarimetric signals. In the Altair analysis, this coupling was implemented in the ESTERIAS pipeline, which constructs a surface grid, interpolates PHOENIX specific intensities and spectra over $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$1, $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$2, and metallicity, applies local projected areas and Doppler shifts, and computes observables such as squared visibilities, closure phases, high-resolution spectra, and pulsation visibilities (Bouchaud et al., 2019).

The Altair study is a benchmark for ESTER’s multi-technique use. An initial $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$3-model fit to interferometry constrained the geometric parameters, after which full ESTER models were explored over mass, metallicity, hydrogen abundance, and core hydrogen mass fraction. The interferometric and spectroscopic parameter search used the MCMC sampler Emcee, while the expensive ESTER grid was interpolated by Delaunay triangulation inside ESTERIAS. Asteroseismology then broke the degeneracy among mass, metallicity, hydrogen abundance, and evolutionary state. The resulting model yielded $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$4, $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$5, $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$6, $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$7, and a flattening $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$8. With $$\rho\,T\,\mathbf{v}\cdot\nabla S = -\nabla\cdot\mathbf{F} + \varepsilon_*,$$9, the inferred equatorial velocity is about $\nabla\cdot(\rho\,\mathbf{v}) = 0.$0, corresponding to a rotation period of about $\nabla\cdot(\rho\,\mathbf{v}) = 0.$1. The same analysis constrained the mass to $\nabla\cdot(\rho\,\mathbf{v}) = 0.$2, with $\nabla\cdot(\rho\,\mathbf{v}) = 0.$3 and $\nabla\cdot(\rho\,\mathbf{v}) = 0.$4, implying an age of about $\nabla\cdot(\rho\,\mathbf{v}) = 0.$5 Myr. ESTER also predicted that Altair’s core rotates approximately $\nabla\cdot(\rho\,\mathbf{v}) = 0.$6 faster than the envelope, while the surface differential rotation does not exceed $\nabla\cdot(\rho\,\mathbf{v}) = 0.$7 (Bouchaud et al., 2019).

The $\nabla\cdot(\rho\,\mathbf{v}) = 0.$8 Aql study illustrates a different role for ESTER. There the primary forward model for polarization and flux was Roche geometry with $\nabla\cdot(\rho\,\mathbf{v}) = 0.$9-model gravity darkening, but rigid-rotation fits produced rotation periods only marginally consistent with the TESS photometric period. ESTER was then used to compute self-consistent two-dimensional structures and surface velocity fields. The predicted surface rotation departs from solid-body rotation at only the $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$0 level, with modest super-rotation at temperate latitudes, and this small shear is sufficient to bring the equatorial rotation period into agreement with the observed $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$1 hr without fine tuning. In the preferred ESTER+Gaia solution, the star has $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$2, $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$3, $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$4, $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$5, $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$6, $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$7, and subsolar photospheric abundances with $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$8 (Howarth et al., 2023).

These applications also clarify a recurrent misconception: ESTER is not merely a more detailed Roche model. Roche geometry with a prescribed gravity-darkening law remains useful for exploratory fitting, but ESTER adds the self-consistent baroclinic structure, flow field, and latitude-dependent surface properties that can shift derived flattenings, modify line-profile interpretations, and predict differential rotation rather than assume it (Bouchaud et al., 2019, Howarth et al., 2023).

5. Evolutionary extensions and oscillation calculations

For many years ESTER was a steady-state code whose evolutionary state was mimicked by varying the core hydrogen abundance $$\mathbf{F} = -\frac{16\,\sigma\,T^3}{3\,\kappa\,\rho}\,\nabla T,$$9 or the contrast between core and envelope composition. This quasi-static strategy was already sufficient to study how rotation changes along the main sequence. In one early application to Be stars, sequences at constant total angular momentum showed that the criticality $\Omega(r,\theta)$0 increases with evolution, suggesting that initially fast-rotating B stars can approach critical rotation during the main sequence if mass loss is weak enough (Rieutord et al., 2013).

That quasi-steady approach was later extended to include anisotropic mass and angular-momentum loss. For models with $\Omega(r,\theta)$1 and masses between $\Omega(r,\theta)$2 and $\Omega(r,\theta)$3, ESTER predicted that stars reach criticality during the main sequence provided their initial angular velocity exceeds $\Omega(r,\theta)$4 of the Keplerian value. For more massive stars, radiation-driven winds and the associated angular-momentum extraction can prevent critical rotation, especially in the two-wind regime associated with a local bi-stability jump of the latitude-dependent mass-flux density. In this framework, the bi-stability jump does not appear as a global discontinuity in $\Omega(r,\theta)$5, but as a latitudinal discontinuity in $\Omega(r,\theta)$6 that migrates with evolution and can strongly enhance low-latitude angular-momentum loss (Gagnier et al., 2019).

A major milestone was the extension of ESTER to genuine time-dependent two-dimensional stellar evolution. The 2023 implementation solves the time evolution of the two-dimensional structure, rotation profile, and composition, including meridional advection, viscous angular-momentum transport, and anisotropic chemical diffusion. For a $\Omega(r,\theta)$7 star, models initialized at $\Omega(r,\theta)$8 and $\Omega(r,\theta)$9 of the critical rotation rate were evolved along the main sequence. The slowly rotating track reproduced the observed properties of the $n=3$0 Cephei star HD 192575, including luminosity, effective temperature, convective-core mass, and the core and surface rotation rates inferred by Burssens et al. The same work found that meridional circulation has a negligible influence on the surface nitrogen enrichment in that case, and that near the terminal-age main sequence the nuclear timescale becomes shorter than the baroclinic relaxation timescale, invalidating any treatment based purely on a succession of steady states (Mombarg et al., 2023).

The intermediate-mass extension sharpened the connection to Be-star phenomenology. Two-dimensional evolutionary tracks for $n=3$1–$n=3$2 stars with initial rotation rates of $n=3$3, $n=3$4, and $n=3$5 of critical were computed at Galactic and SMC metallicity. These models predict that a minimum initial rotation threshold is required to reach near-critical rotation during the main sequence, and that the threshold increases with mass. At Galactic metallicity, $n=3$6–$n=3$7 models can reach critical rotation with $n=3$8, $n=3$9 requires $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$0, and stars of at least $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$1 do not reach near-critical rotation even for $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$2. At SMC metallicity, $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$3–$$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$4 models can reach critical rotation with $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$5, but stars of at least $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$6 still move away from critical rotation. This suggests that the larger Be-star fraction at low metallicity can be reproduced by self-consistent two-dimensional evolution without invoking metallicity-dependent winds in this mass range (Mombarg et al., 2024).

ESTER models have also become equilibrium backgrounds for non-perturbative oscillation calculations. In the adiabatic TOP calculations, continuous and discontinuous ESTER models of rapidly rotating $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$7 stars were used to study island modes, generalized rotational splittings, and acoustic glitches. The resulting spectra show that the island-mode asymptotic pattern remains sufficiently preserved even in the presence of sharp gradients or discontinuities, while the generalized splitting of prograde and retrograde counterparts is well represented by weighted integrals of the ESTER rotation profile except during avoided crossings. This supports the use of ESTER-based rotation kernels for probing $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$8 in rapidly rotating pulsators (Reese et al., 2020).

6. Validation, limitations, and scientific significance

ESTER has been tested against several classes of observations. Early interferometric validation showed that the framework reproduces the fundamental parameters of the two main components of $$T_{\rm eff}(\theta)\propto g_{\rm eff}(\theta)^\beta,$$9 Velorum, while the inability to reproduce Achernar under main-sequence assumptions led to the conclusion that Achernar had likely left the main sequence and was crossing the Hertzsprung gap. Other summaries cite successful applications to Altair, Alderamin, Vega, Rasalhague, Regulus, and the $\omega$0 Vel Aa/Ab pair, establishing ESTER as a practical tool for stars whose oblateness and gravity darkening are directly measurable (Rieutord et al., 2013, Rieutord et al., 2013).

Its limitations are equally explicit in the literature. The standard code is axisymmetric, non-magnetic, and steady-state; convective envelopes are not modeled; the convective core is treated as isentropic; and post-main-sequence structures lie outside the original domain of applicability. The 2016 algorithm paper notes that the Stewartson shear layer on the tangent cylinder is not captured because viscous terms are removed from the meridional Euler balance to avoid Ekman layers. Near-critical rotation remains numerically difficult because the equatorial cusp demands very high angular resolution. Later guides add that the outer atmosphere is simplified, for example through an $\omega$1 polytropic extrapolation to the photosphere, and that stars near or cooler than the Kraft break are presently out of scope because latitude-dependent outer convection is not yet included (Rieutord et al., 2016, White et al., 24 Sep 2025).

Application papers expose further modeling tensions rather than hiding them. In Altair, the atmospheric metallicity required by the Mg II 4481 Å line is higher than the bulk metallicity inferred from the interior model, and the analysis therefore treats atmospheric $\omega$2 and bulk $\omega$3 as independent parameters. In $\omega$4 Aql, the interior ESTER calculations are done at solar abundance while the atmosphere modeling requires $\omega$5; the interpretation is that the star’s subsolar photospheric abundances are likely surface-layer effects and not representative of the bulk interior composition. A plausible implication is that ESTER applications often have to separate interior structure physics from detailed photospheric abundance diagnostics when rapid rotation and peculiar surface chemistry coexist (Bouchaud et al., 2019, Howarth et al., 2023).

Scientifically, ESTER occupies a specific niche. It provides a physically consistent route from baroclinic stellar interiors to observable surface maps, and from those maps to interferometric, spectroscopic, polarimetric, and seismic diagnostics. The steady-state framework established that rapid rotators require true two-dimensional modeling; the evolutionary extensions show that time dependence, chemical evolution, and angular-momentum redistribution can now be treated in the same geometric setting. Across these developments, ESTER has become a reference framework for analyzing fast-rotating stars, especially when the target problem depends simultaneously on shape, flux anisotropy, internal differential rotation, and the coupling of stellar structure to modern high-precision data (Rieutord et al., 2016, Mombarg et al., 2023, Reese et al., 2020).