SPHERLS Code: Multidimensional Pulsation

• SPHERLS Code is an open-source multidimensional hydrocode that couples radial pulsation and turbulent convection using a hybrid Lagrangian–Eulerian approach.
• It employs finite-difference discretization with a Smagorinsky-type SGS model to resolve convective eddy scales and energy fluxes in RR Lyrae envelope simulations.
• The study challenges 1D pulsation-convection closures by revealing complex flux decompositions and convective morphologies, advocating for nonlocal convection treatments.

Searching arXiv for SPHERLS and related RR Lyrae pulsation papers to ground the article in published work. SPHERLS is an open-source multidimensional numerical hydrocode, released under the MIT license and developed for modeling stellar pulsation using LES principles, with a particular emphasis on radial pulsation coupled to multidimensional turbulent convection in stellar envelopes. In the 2025 RR Lyrae study, it is used in a 2D, quasi-static configuration to construct static envelope models, vary horizontal wedge size and resolution, quantify convective eddy scales and energy transport, and optimize the code’s operating regime for later comparisons with 1D pulsation-convection formalisms such as MESA-RSP and the Budapest–Florida codes (Kovács et al., 15 Jul 2025).

1. Definition and computational role

SPHERLS is specifically designed to model radial stellar pulsation together with multidimensional convection. Its central architectural choice is a hybrid Lagrangian–Eulerian moving-mesh formulation: the radial grid moves so that the mass in each shell remains constant, while the multidimensional flow is evolved on that moving background. In the RR Lyrae envelope application, this design is used to preserve resolution in the partial ionization zones while accommodating the large-amplitude radial shell motion characteristic of pulsating stars (Kovács et al., 15 Jul 2025).

The code is formulated as a multidimensional extension of classical radial pulsation hydrocodes rather than as a high-resolution convection code of the type used in modern solar granulation studies. The 2025 analysis explicitly characterizes SPHERLS as “an extended version of [1D pulsation codes], but without the Mixing Length Theory and other 1D approximations.” Earlier SPHERLS work, identified there as SPHERLS I–IV and associated with Geroux & Deupree between 2011 and 2015, had already demonstrated 2D and 3D convection plus radial pulsation, as well as full-amplitude RR Lyrae models and RR Lyrae light curves consistent with observations. The repository cited for the code is https://github.com/cgeroux/SPHERLS (Kovács et al., 15 Jul 2025).

In that sense, SPHERLS occupies a specific methodological niche. It is intended to retain the envelope-physics priorities of stellar pulsation modeling while replacing 1D convective closures with explicitly evolved multidimensional flow plus SGS modeling. This suggests that its main scientific value lies in diagnosing where local 1D closures fail, rather than in providing a DNS-level treatment of stellar turbulence.

2. Governing equations and physical closures

SPHERLS solves the continuity, momentum, and internal-energy equations in spherical geometry. In the notation given for the code, these are

$${\rm D}_t \rho + \rho\,{\rm div}\,\boldsymbol{U} = 0,$$

$${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$

and

$${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$

The derivative operator is the Stokes derivative,

$${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$

where $\boldsymbol{U}_0$ is the grid velocity. This term encodes the moving radial mesh and is therefore central to the code’s hybrid Lagrangian–Eulerian character (Kovács et al., 15 Jul 2025).

Gravity is imposed as spherically symmetric, through the term $-GM_r/r^2\,\boldsymbol{\hat r}$. The equation of state and opacity treatment follow the SPHERLS II setup, using OPAL tables for high and intermediate temperatures and Alexander & Ferguson low-temperature opacities. Radiative transfer is treated in the diffusion approximation,

$$\boldsymbol{F}_r = - \frac{16\sigma T^3}{3 \rho \kappa_R}\,{\rm grad}\,T.$$

The paper notes that this approximation fails in optically thin surface layers and yields atmospheres that are too hot, but also states that SPHERLS is not intended to produce accurate atmospheres for synthetic spectra and that classical pulsation is driven deeper in the envelope by the $\kappa$ mechanism (Kovács et al., 15 Jul 2025).

Unresolved turbulence is represented with a Smagorinsky-type SGS model adapted to compressible stellar flow. The SGS viscosity is written in terms of the grid-scale filter length $l=\sqrt{\Delta r\,r\,\Delta\theta}$, with parameters $L=3.75\,l$, $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$0, $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$1, $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$2, and $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$3. The paper also states that the SGS closure assumes incompressible turbulence at subgrid scales. Its stated regime of validity is therefore one in which turbulent velocities at SGS scales remain subsonic even when the resolved large-scale radial motion becomes supersonic (Kovács et al., 15 Jul 2025).

A key physical diagnostic is the convective flux, defined through a Reynolds/Favre decomposition as

$${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$4

The code analysis then decomposes the enthalpy fluctuation into temperature, pressure, and ionization/composition contributions. In the partial hydrogen ionization region, the study finds that approximately half of the convective flux is supplied by ionization energy transport, implying that the common identification $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$5 is not valid there (Kovács et al., 15 Jul 2025).

3. Geometry, discretization, and boundary conditions

For the RR Lyrae envelope calculations, SPHERLS is used in spherical coordinates $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$6 but restricted to a 2D wedge in $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$7 with no $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$8 dependence. The code therefore implements a box-in-a-star approach: it models a slice of an angular shell rather than a global 3D convective envelope. In this study, the horizontal extent is parameterized by an angular opening $${\rm D}_t \boldsymbol{U} = - \frac{1}{\rho}\,{\rm grad}\,p - \frac{G M_r}{r^2}\,\boldsymbol{\hat r} + {\rm div}\,\boldsymbol{\underline{\tau}}^{\rm SGS},$$9, with explored domain sizes from $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$0 to $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$1 and horizontal resolutions from $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$2 to $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$3 (Kovács et al., 15 Jul 2025).

The reference setup inherited from SPHERLS II is

$${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$4

The radial domain spans the full convective envelope of the RR Lyrae models, including the H and He partial ionization zones, and extends into a deeper radiative region to permit overshooting. The top boundary is placed at the stellar surface with a constant mass closing shell and pressure, while the bottom boundary is located near $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$5 and uses an incompressible boundary condition. Periodic boundary conditions are imposed in the angular direction (Kovács et al., 15 Jul 2025).

Spatial discretization is finite-difference on a staggered mesh, with velocities located at cell interfaces and scalar quantities such as $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$6, $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$7, $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$8, and $${\rm D}_t e = -\frac{p}{\rho}\,{\rm div}\,\boldsymbol{U} - \frac{1}{\rho}\,{\rm div}(\boldsymbol{F}_r + \boldsymbol{F}_{\rm SGS}) + D_{\rm SGS}.$$9 stored at cell centers. Simple derivatives use second-order central differencing, whereas advection uses a donor-cell upwind scheme with a blending parameter $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$0 determined from the local Mach number. At high Mach number the method becomes first-order upwind; at low Mach number it blends toward central differencing. The paper emphasizes that the truncation errors associated with this advection operator behave as an additional nonlinear viscosity and diffusivity, and explicitly derives numerical viscosity and diffusivity terms in Appendix B (Kovács et al., 15 Jul 2025).

Time integration is formally explicit with respect to the advective terms and is constrained by the CFL condition. In the near-surface layers, the study reports time steps as small as $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$1 because of the small scale heights, high sound speeds, and fine horizontal resolution. Initial conditions are obtained by mapping a 1D static radiative envelope into the 2D wedge, setting all velocities to zero, and allowing numerical perturbations to seed convection. The models are then evolved for 40 days, corresponding to roughly $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$2–$${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$3 convective turnover times, followed by a 9-day statistical averaging window (Kovács et al., 15 Jul 2025).

4. Turbulence diagnostics and convective-zone morphology

The 2025 analysis is notable for turning SPHERLS into a diagnostic framework for stellar-envelope convection. Horizontal Reynolds and Favre averages are defined, temporal averaging is introduced under an ergodic assumption, and the resulting mean and fluctuating fields are used to compute RMS velocities, anisotropy, integral length scales, spectra, and flux decompositions. The RMS fluctuations are

$${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$4

with anisotropy parameter

$${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$5

The study finds $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$6 in the main convective regions, consistent with Rayleigh–Bénard-like cells dominated by vertical motion, and $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$7 in counter-gradient and overshooting layers, where the turbulence is nearly isotropic (Kovács et al., 15 Jul 2025).

Convective eddy sizes are characterized with horizontal autocorrelation functions and integral length scales $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$8, using the first zero or first minimum of the autocorrelation to avoid wrap-around effects from periodic boundaries. The corresponding eddy diameters are identified as $${\rm D}_t = \frac{\partial}{\partial t} + (\boldsymbol{U} - \boldsymbol{U}_0)\cdot{\rm grad},$$9 and $\boldsymbol{U}_0$0. This procedure shows that the characteristic eddy sizes in the He II zone are sensitive to the domain width until the wedge becomes large enough to contain the largest structures without artificial truncation (Kovács et al., 15 Jul 2025).

The paper also computes the horizontal turbulent kinetic-energy spectrum from the Fourier transform of the autocorrelation function. In adequately resolved models, the spectrum displays a low-wavenumber driving range, a peak associated with the largest eddies, and a decaying range that is approximately compatible with $\boldsymbol{U}_0$1 or $\boldsymbol{U}_0$2 behavior before entering a steep dissipation range dominated by numerical viscosity. In poorly resolved cases, especially for $\boldsymbol{U}_0$3, the spectrum is truncated at the Nyquist frequency before any inertial-like decay can develop (Kovács et al., 15 Jul 2025).

By combining convective flux, superadiabatic gradient, and turbulent kinetic energy, the study defines a structural classification of the envelope. The convective region is identified by $\boldsymbol{U}_0$4 and $\boldsymbol{U}_0$5, the counter-gradient or Deardorff region by $\boldsymbol{U}_0$6 and $\boldsymbol{U}_0$7, and the overshooting region by $\boldsymbol{U}_0$8 and $\boldsymbol{U}_0$9. On that basis, the quasi-static RR Lyrae envelopes are found to contain two distinct dynamically unstable regions, associated with He II and H I partial ionization, that are loosely connected by a counter-gradient layer. A substantial overshoot layer extends downward from the He II convection zone (Kovács et al., 15 Jul 2025).

5. Numerical regime, effective dissipation, and optimized usage

A central conclusion of the SPHERLS analysis is that the horizontal extent of the computational wedge matters strongly below a threshold of about $-GM_r/r^2\,\boldsymbol{\hat r}$0. For $-GM_r/r^2\,\boldsymbol{\hat r}$1, quantities such as overshoot depth, non-zero turbulence extent, convective velocities, and Nusselt numbers depend strongly on domain size. For $-GM_r/r^2\,\boldsymbol{\hat r}$2, those quantities become more stable, indicating that the largest He II eddies are no longer being artificially constrained by the boundaries (Kovács et al., 15 Jul 2025).

The study quantifies the balance between modeled SGS transport and truncation-error transport by defining effective Reynolds, Péclet, Prandtl, and Nusselt numbers together with ratios comparing numerical and SGS contributions. The reported estimates are

$-GM_r/r^2\,\boldsymbol{\hat r}$3

These values show that numerical diffusion dominates SGS dissipation in the reported runs. The paper therefore interprets SPHERLS, in this parameter regime, as effectively functioning as an iLES in which the numerical scheme itself provides the principal filtering and damping (Kovács et al., 15 Jul 2025).

Resolution requirements are then inferred directly from the convection diagnostics. The analysis states that $-GM_r/r^2\,\boldsymbol{\hat r}$4 is needed to capture granulation-like eddies and avoid excessive numerical damping, whereas $-GM_r/r^2\,\boldsymbol{\hat r}$5 suppresses granulation characteristics. It further suggests that $-GM_r/r^2\,\boldsymbol{\hat r}$6 would yield a more realistic onset of convection. By contrast, the vertical resolution typical of 1D pulsation modeling, approximately 150–200 radial shells, is described as sufficient for resolving the ionization fronts and convective layers relevant to the study. The practical optimum proposed for future RR Lyrae work is either $-GM_r/r^2\,\boldsymbol{\hat r}$7 with $-GM_r/r^2\,\boldsymbol{\hat r}$8 or, as a compromise, $-GM_r/r^2\,\boldsymbol{\hat r}$9 with $$\boldsymbol{F}_r = - \frac{16\sigma T^3}{3 \rho \kappa_R}\,{\rm grad}\,T.$$0 (Kovács et al., 15 Jul 2025).

A common misconception would be to read the presence of an SGS closure as evidence that SPHERLS is operating in a cleanly resolved LES regime. The numerical analysis in the paper argues otherwise for the present setup: the SGS model is formally present, but numerical damping is the dominant effective dissipation channel.

6. Relation to 1D pulsation theory, limitations, and scientific significance

The SPHERLS results are explicitly framed as preparatory for later 2D–1D comparisons, but the implications for 1D pulsation-convection closures are already substantial. First, the envelope structure recovered in SPHERLS is more complex than a local MLT picture. Instead of a single simply parameterized convective layer, the quasi-static RR Lyrae models display two unstable zones, associated with He II and H I ionization, connected by a Deardorff region in which the convective flux remains positive while the superadiabatic gradient is negative. This violates local down-gradient assumptions commonly embedded in MLT and in many time-dependent convection formalisms (Kovács et al., 15 Jul 2025).

Second, the flux decomposition shows that in the H I partial ionization region roughly half of the convective transport is supplied by the ionization or latent-heat term $$\boldsymbol{F}_r = - \frac{16\sigma T^3}{3 \rho \kappa_R}\,{\rm grad}\,T.$$1. That result directly challenges closures that identify the convective flux primarily with a temperature covariance. Pressure-flux terms are also found to matter, depending on whether the enthalpy fluctuation is written in entropy-based or temperature-based form. The paper therefore concludes that the structure of the convective zone urges reconsideration of recent approaches to convective-flux modeling in radial stellar pulsation codes (Kovács et al., 15 Jul 2025).

Third, the treatment has important limits. The calculations are 2D rather than 3D, and the paper notes that 2D turbulence has a different cascade and longer-lived structures than 3D turbulence. The models are static or quasi-static rather than fully nonlinear pulsation calculations, because the runs are terminated before pulsation amplitudes become large. Radiation is handled with diffusion rather than a detailed non-gray treatment in optically thin layers. The wedge geometry excludes non-radial global modes and fully global convective organization. These caveats do not negate the results, but they restrict how directly the diagnosed convective morphology can be mapped onto fully developed stellar turbulence (Kovács et al., 15 Jul 2025).

Within those limits, SPHERLS functions as a physically consistent multidimensional extension of 1D pulsation hydrocodes. Its strongest contribution is not the realization of fully resolved stellar turbulence, but the extraction of geometric, statistical, and energetic information—eddy sizes, anisotropy, overshoot depth, counter-gradient transport, and flux decomposition—that exposes where standard 1D assumptions break down and where future radial pulsation codes may need nonlocal or non-diffusive convection formalisms.