© The Authors 2025

1. Introduction

Convection is a buoyancy-driven transport process (Glatzmaier 2014), formed when the temperature gradient is greater than the adiabatic temperature gradient in stars¹. Hence, small fluctuations in density cause mass elements to rise and sink nearly adiabatically, creating a 3D dimensional flow that is random, but has quasi-stationary structures called Rayleigh-Bénard cells. The flow inside these cells carries energy from the hotter bottom of the cell to the cooler top where the extra energy is released to the surroundings through heat diffusion and viscosity. This way, it transports energy, decreasing the temperature gradient of the system. In a closed system, this can sufficiently decrease this gradient to a level where the convection can stop, and the kinetic energy of the cells would dissipate. However, in an open system, for example a stellar envelope, the temperature difference between the two sides of the convectively unstable region does not disappear, so buoyancy maintains the convective cells. This phenomenological picture was formulated in many different versions of mixing length theories (MLT).

However, in stellar envelopes (also in stellar cores and other astrophysical fluids), the material viscosity of the fluid is very low. This causes the structure of Rayleigh-Bénard cells to decompose into smaller eddies, which also decompose, and so on, until they reach small-scale motions that are effectively braked by viscous forces. This cascade of eddies shows an inherently unpredictable chaotic flow structure that is called turbulent flow (Pope 2000), and its appearance can be characterized by a dimensionless number, the Reynolds number, determined by the fluid physical parameters². This number, by definition, is the ratio of the inertial to viscous forces, Re = uL/ν, with a characteristic velocity, u, the characteristic length scale, L, and the kinetic viscosity, ν, of the flow. If Re ≫ 1, a flow becomesturbulent (e.g., its value is Re ∼ 10¹² (Kupka & Muthsam 2017) in the solar envelope). Hence, convection in stars is always turbulent, but we emphasize that the two phenomena are different, and one can occur without the other, outside the scope of stellar astrophysics.

These are genuinely 3D phenomena, and their description in 1D pulsation models can only be achieved through approximations. First, one needs to use averaging on the 3D Navier-Stokes equations to derive horizontal mean quantities that are changing with the pulsations (such a derivation is presented by Kuhfuss 1986). These are called ensemble averages (Pope 2000 and for usage in other astrophysics problems, see Kupka & Muthsam 2017). Unfortunately, the method always includes unspecified terms; therefore, the 1D equations are not closed. This closure problem can be solved with simplifications, and the level of complexity depends on the number and quality of the used approximations. In pulsation theory, these include the Boussinesq approximation meaning neglecting density fluctuations except for buoyancy (first used for pulsations by Gough 1969), the neglect of acoustic modes, the eddy viscosity hypothesis (Kuhfuss 1986) simplifying pressure terms, and the gradient-diffusion (or down-gradient) hypothesis (Stellingwerf 1982), which relates the unsolved transport terms to the gradient of the mean-field. We refer to Kupka & Muthsam (2017) and Houdek & Dupret (2015) for background material.

The other solution for the problem would be 3D modeling, directly solving the Navier-Stokes equations. The direct numerical simulation (DNS) is still too computationally expensive on modern computers as one would need to resolve scales of the flow from scales of several million kilometers to the scales of meters (Kupka & Muthsam 2017). Instead, one can use the so-called large eddy simulation (LES) method, pioneered by Smagorinsky (1963), Lilly (1969), and Deardorff (1970). The main idea is to resolve the large-scale, mainly anisotropic (in our case, buoyancy-driven convection) motions and model only those on this scale. These modeling results are similar to measurements for a given resolution, or more likely, a filtered view of good resolution data. Regardless of the details, the filtering is similar to the ensemble averaging of the 1D models, although it has different properties (Pope 2000). The method itself has been widely used in engineering and meteorology since the sixties (Pope 2000), while it has had many successes in the simulation of solar granulation and solar physics (Nordlund et al. 2009), it is used for modeling A type stellar envelopes (Fabbian et al. 2023), studies of stellar cores (see, e.g., Rizzuti et al. 2023, Georgy et al. 2024), AGB stars (e.g., Freytag et al. 2012), and stellar convection in general (Meakin & Arnett 2007; Arnett et al. 2009; Viallet et al. 2013)³. Kupka & Muthsam (2017) gives an overview of the history of LES in astrophysics⁴. However, we emphasize that expectations that the 3D simulations will solve all the problems via brute force are misleading.

There are two approaches that are widely used in the aforementioned studies: the so-called box-in-a-star and the star-in-a-box method. In the former, one defines a simulation box(horizontally with periodic boundary conditions) inside the star (see, e.g., Stein & Nordlund 2001; Muthsam et al. 2010) because the length scales of the convective motions are not comparable to the size of the star (this is the case in main-sequence stars and classical pulsators). The latter case is used for red and asymptotic giant branch stars, where the convective length scales are comparable to the stellar radius (e.g., Freytag et al. 2012).

In the case of radial stellar pulsations of classical variables, a further problem occurs, namely the large amplitude radial oscillation of the outer layers, where the convection zone resides. This means that any codes need to follow the movement of the stellar layers; meanwhile, convection needs to remain resolved. If the former condition is not met, one cannot reach full-amplitude pulsation with the model, as the important partially ionization zones moves out from the grid (Deupree 1977a). The first attempt to reproduce 2D pulsation and convection together was done by Deupree (1977a), successfully reproducing the red edge of the instability strip (Deupree 1977b, 1980, 1985). Recently Mundprecht et al. (2013) implemented an adaptive method into the ANTARES code (Muthsam et al. 2010) which was used to simulate Cepheids (Mundprecht et al. 2015). Vasilyev et al. (2018) used the COBOLD5 software (Vasilyev et al. 2017, and references therein) to simulate a 2D Cepheid atmosphere as an input for spectral synthesis models. Geroux & Deupree (2011), Geroux & Deupree (2013) developed the SPHERLS code based on the ideas of Deupree (1977a) and successfully reached full-amplitude RR Lyrae models (Geroux & Deupree 2014, 2015), and we are using this code in our current study.

In this paper in our series, we optimize the SPHERLS code input parameters (resolution, model size) and determine the structure of the convective zone through a series of 2D quasi-static simulations to establish a theoretical ground for the comparison with 1D. For this, we use statistical methods and dimensionless numbers largely used in turbulence theory (Pope 2000), and define some new to follow the effects of the numerics. For example, we are using nondimensional numbers, such as the Prandtl number, that characterize the fluid and can be compared to numbers of real stellar flows derived from material science (Kupka & Muthsam 2017), while we are using numbers that characterize the flow, for example the Reynolds number, which is also dependent on the simulation setup.

In addition to the statistical quantities, we define other quantities by describing the size and structure of the convective zone, as these are comparable directly to 1D models. Only the consistency of these and statistical quantities can ensure that a comparison with 1D calculations, and the conclusions derived from it, would be meaningful.

In our previous papers (Kovács et al. 2023, 2024, and references therein), we confronted two 1D models (calculated with the Budapest-Florida and MESA-RSP codes) of RR Lyrae stars in M3 globular cluster with observed radial velocities, calibrated their free parameters, and analyzed the processes occurring in them during pulsation. This work focuses on technical questions of 2D modeling of these stars, and the following papers will discuss in detail the results provided by the comparison of 1D codes and the multidimensional SPHERLS in 2D for static and pulsating states.

The structure of this paper is the following. First, we introduce the main concepts behind the SPHERLS code alongside with the stellar parameters of the input models in Section 2. Then we define our analytical methods in Section 3: The modeling procedure (Sect. 3.1), the statistical methods that are describing turbulent flows (Sect. 3.2), and the various dimensionless quantities related to turbulence, numerics, and 1D treatment (Sect. 3.3). This is followed by the description of our results. First, we present the numerical costs of our calculations in Sect. 4.1, then we describe the vertical structure of the convective zone in the studied quasi-static RR Lyrae models in Section 4.2. We studied the different components of the convective flux, which is presented in Section 4.3. The effects of resolution and horizontal model size on the structure of the velocity field, on the size of the eddies and on the behavior of the energy cascade is given in Sections 4.4, 4.6, and 4.7, respectively. Our results and their consequences are discussed in Section 5, and our conclusions and a summary are given in Section 6.

2. The multidimensional SPHERLS model

SPHERLS is an open-source multidimensional numerical hydrocode⁵ that was developed by Geroux & Deupree (2011, 2013) for modeling stellar pulsation using LES principles. This code is based on a hybrid Lagrangian-Eulerian scheme, which makes it possible to follow the long-term and short-term changes in radially pulsating variable stars while maintaining a good enough resolution of the ionization regions. We repeat the governing equations in Appendix A to match our notations.

In the field of classical stellar pulsations, the main focus is on the optically thick stellar envelope, in which the partial ionization zones that drive the pulsations reside. This means that the atmosphere itself is used as a boundary condition and radiation flux is described through the diffusion approximation. SPHERLS has limited capabilities regarding horizontal resolution (see Sect. 4.1), even with simple numeric methods (see Appendices B and C). With horizontal resolution being magnitudes lower than current hydro-codes used in solar convection modeling, and following the same principles as 1D pulsation codes, this modeling method can be seen as an extended version of the latter, but without the Mixing Length Theory and other 1D approximations. So this code is suitable to study stellar convection in pulsating stars, keeping in mind the aforementioned limitations we have.

The three primary numerical concepts crucial to our study are the box-in-a-star nature of SPHERLS, the utilization of the subgrid scale viscosity model and the effects of numericaldamping.

The box-in-a-star approach of SPHERLS means that we are modeling a slice of an angular shell of the star. Geroux & Deupree (2015) have shown that 2D calculations are adequate to describe RR Lyrae variables. Therefore, we calculated 2D wedges, keeping in mind that in 2D the average lifetime of a convective structure can be longer (Kupka 2020) and the cascade behavior is different than in the 3D case (Kupka 2009). On the horizontal boundaries, we use periodic boundary conditions; the top boundary of the star is determined at the surface with a constant mass closing shell and pressure (Geroux & Deupree 2013). The bottom boundary is taken at 2 million K temperature with an incompressible boundary. These vertical boundaries are the same which are used in radial pulsation models in general (see, e.g., Paxton et al. 2019), and the stellar parameters define the vertical extent of the model. In contrast, the horizontal extent is always chosen to an arbitrary angular width that we call θ_width. This horizontal extent is defined by the number of horizontal cells, N_θ, and the horizontal size of a single cell, Δθ, as θ_width = N_θΔθ. The choice of these values is not straightforward, as one has to choose them to have a θ_width large enough to capture the large-scale convective processes but with small enough Δθ to maintain the necessary resolution.

Another essential concept is the subfilter or subgrid scale (SGS) model. Since the simulation only solves equations on the so-called resolved scales, the effects on smaller scales than the grid scales have to be modeled somehow. Such an effect is, for instance, the turbulent cascade itself: the kinetic energy of the resolved scales is transferred to these SGS scales, and without proper treatment, the simulations fail. The SGS model usually defines the filter function of the LES equations, as well. The length scale of this function Δ_filt is coupled to the resolution where the LES method gives information (Magnient et al. 2007). The length scale regarding the dissipation of energy (Δ_diss) is related to Δ_filt by the applied SGS model. This, in practice, means that an LES calculation shows the dissipation of kinetic energy at a scale of a different length than the actual flow. In the case of the SPHERLS, chose Δ_filt = l where l is the length scale of the grid (see Appendix A).

SPHERLS uses the Smagorinsky-model (Smagorinsky 1963), which connects the viscous effects of SGS turbulence to the velocity stress tensor, similar to the eddy viscosity models in 1D; this also assumes incompressible turbulence. Although we perform compressible hydrodynamics, one can assume that on the subgrid scales, the fluid is incompressible. This assumption is valid if the related turbulent velocities are subsonic (Garnier et al. 2009) even in the case when the background radial pulsation has velocities above the speed of sound.

Lastly, numerical damping is a serious issue, especially in low-resolution models. This damping is caused by the numerical scheme of every hydrocode, and its main source is the truncation error of the spatial derivatives. This effect also can be stronger than the effects of the SGS models. We discuss this in detail in Appendix B. If this damping is too strong, shockwaves cannot be modeled. This could be an issue in the modeling of the ionization front, for which case we performed additional numerical tests next to the original ones of Geroux & Deupree (2011), which we present in Appendix C.

Similar to Geroux & Deupree (2013), we use tabulated opacity and equation of state values by OPAL (Iglesias & Rogers 1996) and the low temperature Alexander & Ferguson (1994) tables.

We ran a series of 2D models with two different sets of stellar parameters. These represent the RRab stars v036 and v046 in the M3 globular cluster, which were used for 1D model calibration by Kovács et al. (2023). We present the input parameters in Table 1.

Geroux & Deupree (2013) used originally an angular size of θ_width = 6° and horizontal cell number of N_θ = 20, claiming that it is enough for one upward and downward moving convective cells to form. We refer to this model as the reference model throughout our paper. We chose a number of different angular widths for the models, covering the range from 4° to 45°. The number of horizontal cells chosen is between 20 and 60; this way we ran bad-resolution and small-sized models to test the effects of these arbitrary parameters.

The vertical resolution of our models is based on the 1D resolutions of pulsating stars (Paxton et al. 2019): we have a total N_r number of radial Lagrangian shells from which N_u is uniformly dispatched above the so-called anchor temperature, which is chosen to be T_anchor = 11 000 K in the partial hydrogen ionization region. Below this temperature, the shell sizes are increased according to a geometrical series. We ran 3 different vertical resolutions: N_r = 150, N_u = 50 is our baseline resolution (N_r = 152 in the case of v046 models), we have a few models calculated with N_r = 122, N_u = 25, and a high-resolution model with N_r = 200, N_u = 75, T_anchor = 50 000 K. Meanwhile, due to the very small cell sizes and the otherwise damped horizontal motions, SPHERLS uses 1D cells at the bottom of the model. We have 25 1D cells in the bottom of our baseline of low vertical resolution models, and 50 cells in the larger vertical resolution models. The different model dimensions are presented in Table 2.

Due to this setup, the aspect ratio of the grid cells are changing throughout the stellar envelope: since the horizontal size is ∼r, we have a finer horizontal grid deeper in the star, while in the region above T_anchor the vertical resolution is finer, with aspect ratios Δr: rΔθ ∼ 1: 40 or larger. This directly affects the numerical damping and SGS viscosity: in some regions we have stronger horizontal damping then in other regions. We present aspect ratios of our models in Table D.1, in Appendix D.

3. Methods

We restricted our investigation to quasi-static envelopes before pulsation is triggered in them. In these envelopes, convection, and turbulence already develop spontaneously, which can thus be studied purely in themselves and their characterization is our aim in this study.

3.1. The modeling procedure

The discussed models are quasi-static which means that we started the hydrodynamic calculations from the zero velocity state in contrast with the method used in pulsation physics (see, e.g., Smolec & Moskalik 2008a; Geroux & Deupree 2013; Kovács et al. 2023). As these are pulsationally unstable envelopes, which means that a small perturbation (caused by either the convective motions or the numerical errors) can excite them, the models could not be run for too long; some models (e.g., v036_6_150x40) showed small amplitude but detectable pulsations at the end of the simulations. As the main focus was on the LES methodology of the SPHERLS code and not on the pulsation for the current study, we adjusted the simulation times accordingly.

All models were started from a fully radiative static initial structure (following the procedure of Geroux & Deupree 2013), and convection is ignited automatically by the numerical truncation errors in quantities of the models. Since the convective regions in RR Lyrae variables are relatively shallow, the convective turnover timescale (τ_conv) is 2–5 days. To reach a steady state of the convection, we ran the models for 40 days of simulation time, which sets the relaxation timescale (Kupka & Muthsam 2017) τ_rel ∼ 8 − 10 τ_conv. The statistical analysis of the models started after 40.5 days, and the statistical timescale was chosen to be 9 days (τ_stat ∼ 2 − 3 τ_conv); therefore, the statistical data were collected over the interval of ∼(40 − 50) days. The pulsation amplitude starts to increase only after that point, and was where we stopped the simulations.

In Fig. 1. we present the total turbulent kinetic energy (hereafter TKE) evolution of six different realizations of the v036 model. The differences are in the angular width (left panel) and angular resolution (right panel). The start of the statistical analysis is denoted by the vertical dashed line. The total turbulent kinetic energy is the kinetic energy of the shells projected onto the entirety of the star,

Fig. 1. Time evolution of the total turbulent kinetic energy (TKE). Left panel: Comparison of the TKE evolution for different opening angles with resolution of Δθ = 0.3°. The models are v036_6_150x20 (blue solid line), v036_9_150x30 (yellow dashed line), and v036_12_150x40 (green dotted line). Right panel: Same as left panel, but with the same opening angle (θwith = 6°) and different resolutions: v036_6_150x20 (blue solid line), v036_6_150x30 (red dash-dotted line), and v036_6_150x40 (purple dotted line). The vertical dashed line denotes the end of τrelax and the start of τstat.

$$\begin{array}{c}\hfill E_{\mathrm{TKE}}={\mathcal{S}}^{\star }\frac{1}{2}{\displaystyle \underset{R_0}{\overset{R}{\int }}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_0+{\theta }_{\mathrm{width}}}{\int }}}\rho (r,\theta ){|\mathit{U}(r,\theta )-{\mathit{U}}_0(r)|}^2{r}^{2}sin\theta \mathrm{d}\theta \mathrm{d}r,\end{array}$$(1)

where ρ is the density, U is the velocity field, U₀ is the grid velocity (i.e., the velocity which is used to maintain the shell masses, see Appendix A) and S^⋆ = 2π/(θ_width∫_θ0^{θ₀ + θ_width}sinθdθ) is a scaling factor used to upscale the energy from the shell to the entire star. Here θ_width is the angular size of the model, R₀ is the inner radius of the envelope, R is the radius of the star, θ₀ = π/2 − θ_width/2 is the starting horizontal (θ) coordinate of the model.

The evolution of this quantity over time gives us a picture of the state of the model. First, it shows a somewhat steady rising phase, which comes from the initial small perturbation, and after that a relaxation phase in which it reaches a quasi-steady phase in which it fluctuates but has an almost constant average in time. Essentially, we can say that the first two phases are somewhat dependent on the resolution of the model, which is related to the numerical viscosity of the given model. This can be seen in Fig. 1 right panel: The larger numerical viscosity coming from the larger horizontal cell sizes of the model v036_6_150x20 (N_θ = 20, Δθ = 0.3°) delay the rising phase. Nevertheless all models reached the quasi-steady state after τ_rel.

3.2. Basic statistical quantities

There are several methods to analyze turbulent convection. In this section, we present some of the statistical methods that are broadly used in the field of turbulence.

In the literature of numerical hydrodynamical modeling of pulsation, the meaning of the terms of “ergodicity” and “ensemble average” are slightly different than in other branches of physics, for example in statistical physics (Kupka & Muthsam 2017). In this type of numerical solution, the computation time has to cover many periods or characteristic time intervals to provide good enough statistics, which requires a long computation time. Another possibility is that one runs several copies of the model, and instead of taking the average of one long sample, one takes the average of these copies. The ergodic hypothesis of numerical turbulence states that these two methods give identical results.

Hence, to collect statistical data about the turbulent convection, we assume ergodicity; so, ensemble averages are taken as time averages of τ_stat = 9-day interval (i.e., instead of averaging a static q quantity of N number of simulations, we take the time average of a single simulation). The ensemble averages are denoted by

$$\begin{array}{c}\hfill ⟨q⟩=\frac{1}{N}{\displaystyle \sum _i^N}{q}_i\approx \frac{1}{{\tau }_{\mathrm{stat}}}{\displaystyle {\int }_{t_0}^{t_0+{\tau }_{\mathrm{stat}}}}q(t)\mathrm{d}t.\end{array}$$(2)

In different cases, we use horizontal Reynolds and Favre averages. The former is denoted by the overbar,

$$\begin{array}{c}\hfill \overline{q}\equiv \frac{1}{\mathrm{\Delta }\mathrm{\Omega }}{\displaystyle \stackrel{\mathrm{\Delta }\mathrm{\Omega }}{\int }}q(r,\theta ,\varphi ,t)\mathrm{d}\mathrm{\Omega }=\frac{{\displaystyle \stackrel{{\theta }_{\mathrm{width}}}{\int }}q(r,\theta ,t)sin\theta \mathrm{d}\theta }{{\displaystyle \stackrel{{\theta }_{\mathrm{width}}}{\int }}sin\theta \mathrm{d}\theta },\end{array}$$(3)

and the Favre average is defined by the density-weighted Reynolds average and denoted by a tilde:

$$\begin{array}{c}\hfill \stackrel{\sim }{q}=\frac{\overline{\rho q}}{\overline{\rho }}·\end{array}$$(4)

We note that the ensemble and horizontal averaging are commutable. To mitigate the fluctuations in the statistics caused by numerical uncertainties, we use ensemble averages of the horizontally averaged quantities as well, a procedure is widely used in the field (Kupka & Muthsam 2017). Using the above averaging methods, we can split quantities into mean and fluctuating parts:

$$\begin{array}{c}\hfill q=\overline{q}+q\prime =\stackrel{\sim }{q}+q″.\end{array}$$(5)

We note that using fluctuating quantities allows the computation of other statistical moments: for instance, the horizontal standard deviation of a quantity q is defined by $\overline{{\left(q\prime \right)}^2}$.

The Favre and Reynolds averaging and their fluctuations are related to each other as

$$\begin{array}{cc}\hfill \stackrel{\sim }{q}& =\overline{q}+\frac{\overline{\rho \prime q\prime }}{\overline{\rho }};\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill q″& =q\prime -\frac{\overline{\rho \prime q\prime }}{\overline{\rho }}·\hfill \end{array}$$(7)

One can see that in the case of ρ′/ρ ≪ 1, the additive terms on the righthand side disappear, thus the two averaging and corresponding fluctuations are equivalent. These averages are usually used to analyze the multidimensional models in 1D and also to construct 1D models (Viallet et al. 2013).

To tackle the turbulent nature of stellar convection, we use statistical quantities usually used to describe turbulence. The structure of the velocity field can be described by the auto-correlation function⁶. The auto-correlation function at position x is the ensemble average of the velocity fluctuations and the fluctuation at x + r_d, where r_d is a displacement vector (Pope 2000):

$$\begin{array}{c}\hfill R_{\mathit{ij}}({\mathit{r}}_d)=⟨u_i(\mathit{x}+{\mathit{r}}_d)u_j(\mathit{x})⟩.\end{array}$$(8)

Here u_i = U_i − ⟨U_i⟩ is the ensemble fluctuation of the velocity in i direction. Since most Reynolds-averaged models used in the pulsation-literature (see Baker 1987) assume horizontal symmetry, and we mostly consider the horizontal (θ_d) displacements, r_d = rθ_de_θ. On the other hand, numerical fluctuations break the strict homogeneity; thus we use the Reynolds averages as well: $\overline{R_{\mathit{ij}}}$.

We show vorticity (ω = rot(U − U₀)) maps of different models in Fig. 2. All models are different realizations of v036. The left panel shows realizations with the same resolution and different widths, while the right panel shows different resolutions with the same angular width. We note that we included an animated version of this figure as supplementary material, in which one can follow the change of the vorticity field in τ_stat time. In this video, one can see that the different blobs move upward and downward, and there are larger ones that also rotate. To quantify the characteristic length scales of these eddies, we use the integral length scale:

Fig. 2. Snapshot of five different simulations, taken at the same simulation time. The 2D simulations are in polar coordinates; the radial coordinate is in 1011 cm units; and the angular coordinates are in degrees. The color-coding denotes the vorticity of the fluid flow The simulation domains are rotated to fit on the figure. Left panel: Models v036_6_150x20 (top), v036_9_150x30 (middle), and v036_12_150x40 (bottom). These models have different angle widths, but the same resolution. Right panel: v036_6_150x20 (top), v036_6_150x30 (middle), v036_6_150x40 (bottom). These models have the same angle width, but different horizontal resolutions. In the left panel we can see that with increasing opening angle, the number of structures is also increasing in the HI region (at around 3.4 × 1011 cm radius.). In the HeII zone (at around 3.2 × 1011 cm), the size of the vortices increase. This means that the initial model with of 6° degrees used by Geroux & Deupree (2013) is inadequate to properly capture the properties of the flow structure in this area. Meanwhile, we can see on the right-hand side that the number of structures in the HI region is similar with higher resolutions (although more details appear). An animated version of this figure is available online.

$$\begin{array}{c}\hfill L_{\mathit{ij}}={\displaystyle \underset{0}{\overset{r_{max}}{\int }}}{R}_{\mathit{ij}}(r_d)/R_{\mathit{ij}}(0)\mathrm{d}{r}_d.\end{array}$$(9)

In engineering cartesian scenarios, the upper boundary r_max is set to infinity (Pope 2000). This cannot be maintained in stellar astrophysical scenarios, since in the radial direction it makes no sense to continue the integration into the interstellar medium, while in the horizontal domain the star has a finite extent of 2πR, which is a periodic boundary resulting in a self interaction. In our case the horizontal extent is much smaller being θ_widthR. The periodic boundary also modifies the shape of the autocorrelation function, thus we select r_max based on the following conditions (O’Neill et al. 2004):

• we integrate until the auto-correlation function reaches zero,

• or until it reaches a minimum value.

These statistics give us a quantity that can be dependent on the horizontal width and resolution. If the angular size of the model is too small, it prevents the creation of larger eddies. The nonlinear interaction between convection and pulsation causes cycle-to-cycle variations in the pulsation (alongside with granulation noise), which can be severely modified by unphsysical self-interaction caused by the too small horizontal size. Having the current limitations regarding θ_width, we use L_ij to minimize this effect, as this length scale does not increase further after reaching a minimum optimal angular size (O’Neill et al. 2004).

In addition to the length scale of the eddies, the overall horizontal fluctuation of the velocity field (or root-mean-square difference velocity) is an important quantity, as it is used directly in the 1D pulsation codes, so it can be compared to 1D results. Hence, it is crucial to map this dependence on the horizontal model size, θ_width and horizontal resolution, Δθ. We define these velocities by:

$$\begin{array}{cc}\hfill u_{r,\mathrm{rms}}& \equiv \sqrt{\langle \overline{{\left({U\prime }_r\right)}^2}\rangle }=\sqrt{2e_{22}},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill u_{\theta ,\mathrm{rms}}& \equiv \sqrt{\langle \overline{{\left({U\prime }_{\theta }\right)}^2}\rangle }=\sqrt{2e_{11}}.\hfill \end{array}$$(11)

We note that by definition, these are equivalent to the square root of twice the horizontal and radial specific turbulent kinetic energy densities (e₁₁ and e₂₂ , respectively)⁷. The total specific kinetic energy is defined by

$$\begin{array}{c}\hfill e_{t,0}=e_{11}+e_{22}=\frac{1}{2}\langle \overline{{({U\prime }_r)}^2+{({U\prime }_{\theta })}^2}\rangle .\end{array}$$(12)

Since the Rayleigh-Bénard cells of the convection are genuinely anisotropic, and these anisotropies have their effects in the pressure terms of the turbulence (Pope 2000), the large scale anisotropy of the flow is also interesting. This is defined by the ratio of the rms velocities:

$$\begin{array}{c}\hfill {\sigma }_{r\theta }\equiv \frac{u_{r,\mathrm{rms}}}{u_{\theta ,\mathrm{rms}}}·\end{array}$$(13)

Convection in the stellar matter is always turbulent, as the viscous forces are negligibly small (Kupka & Muthsam 2017). Meanwhile, the radiative diffusion is large and the Prandtl numbers are in the 10⁻⁷ − 10⁻¹⁰ range. Turbulent convection with these low Prandtl numbers follows Kolmogorov-scaling (Verma et al. 2011), meaning that kinetic energy of the large scale convective motions are transferred through a cascade of eddies to scales where viscosity is relevant (Pope 2000). In this case, the kinetic energy in the cascade follows a power law of −5/3. Thus, the spatial spectra of the turbulent kinetic energy can be divided into three spatial ranges: the energy-containing range, which contains the bulk (> 80% Pope 2000) of the kinetic energy; the inertial subrange follows the Kolmogorov-law and the dissipation range, where the kinetic energy decreases exponentially, and is transformed into heat. In LES simulations, the main focus is to resolve the energy-containing range, which is done by a filtering function; however, this filter affects damping energies. To minimize the damping effect, one has to maintain a filter resolution in which Δ_filt ≤ L₁₁. In our case ${\mathrm{\Delta }}_{\mathrm{filter}}=l=\sqrt{r\mathrm{\Delta }r\mathrm{\Delta }\theta }$ (Magnient et al. 2007).

The other source of damping is numerical diffusion or numerical viscosity, which is caused by the truncation errors of the numerical scheme of the code. This actually can be more severe than the damping effect of the filter (see Appendices B and C). Therefore, it is worth to study the energy cascade of these models and compare it to the literature.

To quantitatively describe the energy cascade, we define the energy spectrum tensor by the Fourier transform of the auto-correlation function (Pope 2000),

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathit{ij}}(\kappa )=\frac{1}{\mathcal{L}}{\displaystyle \underset{0}{\overset{\mathcal{L}}{\int }}}{R}_{\mathit{ij}}(r_d,t)e^{-i\kappa r}\mathrm{d}{r}_d,\end{array}$$(14)

where $\mathcal{L}$ is the horizontal width of the model ℒ = rθ_width With having N horizontal cells, the grid spacing horizontally is θ_width/N_θ = Δθ, which translates to a cut-off scale of rΔθ ≡ Δ_c. In practice, the largest horizontal frequency (Nyquist-frequency) is at wave number κ_c = π/Δ_c, and one can obtain N/2 number of spectrum points. Following this concept, the spatial spectrum of the turbulent kinetic energy is

$$\begin{array}{c}\hfill e_t(\kappa )=\frac{1}{2}({\mathrm{\Phi }}_{11}(\kappa )+{\mathrm{\Phi }}_{22}(\kappa )).\end{array}$$(15)

We note that by using horizontal means, all previously described quantities will be dependent on the radius coordinate inside the star.

3.3. Dimensionless quantities

Although we do not have direct measurements of RR Lyrae subsurface convection, one can derive dimensionless numbers from material physics to characterize flow in these layers (Kupka & Muthsam 2017)⁸. The other and commonly used approach is comparing modeling results to observations. However, we need to follow many pulsation periods in nonlinear radial pulsation modeling, which has an unfeasible numerical cost in the case of two or three dimensional computations. One would need multiple observed stars, and for each, a dozen models followed through at least 400–500 cycles, which calculation for a single model can take months ( for the computational cost in SPHERLS, see Section 4.1 and also the original SPHERLS papers of Geroux & Deupree 2011, 2013, 2014; Geroux & Deupree 2015). Hence, we use dimensionless numbers that characterize the flow and underlying numerics.

The dimensionless numbers characterizing the flow are the ratio of effective kinematic to radiative energy diffusion (effective Prandtl number, Pr_eff), the ratio of inertial to viscous forces (Reynolds number, Re), the ratio of total energy flux to radiative energy transport (Nusselt number, Nu) and the ratio of convective velocities to radiative diffusion (Péclet number, Pe). We also define the horizontal Reynolds number (Re_hor) in which we only consider horizontal velocities and horizontal viscosity.

Their definitions:

$$\begin{array}{cc}\hfill P{r}_{\mathrm{eff}}& =\frac{{\nu }_{\mathrm{num}}+{\nu }_{\mathrm{SGS}}}{{\chi }_{\mathrm{tot}}},Re=\frac{u_{\mathrm{rms}}2L_{11}}{{\nu }_{\mathrm{num}}+{\nu }_{\mathrm{SGS}}},\hfill \\ \hfill R{e}_{\mathrm{hor}}& =\frac{u_{\theta ,\mathrm{rms}}2L_{11}}{{\nu }_{\mathrm{num},\theta }+{\nu }_{\mathrm{SGS}}},Nu=\frac{F_r+F_c}{F_r},Pe=\frac{u_{\mathrm{rms}}2L_{11}}{{\chi }_{\mathrm{tot}}},\hfill \end{array}$$

where ν_num is the kinetic diffusivity of the numerics, ${\chi }_{\mathrm{tot}}=\chi +{\chi }_{\mathrm{SGS}}+{\chi }_{\mathrm{E}}^{\mathrm{num}}+{\chi }_{\mathrm{T}}^{\mathrm{num}}$, with χ = 16σT³/(3κc_pρ²) which is the thermal diffusivity. The quantities ${\chi }_{\mathrm{E}}^{\mathrm{num}}$ and ${\chi }_{\mathrm{T}}^{\mathrm{num}}$ are the internal energy and heat diffusivity of the numerics, while χ_SGS = ν_SGS/(Pr_t) which is the subgrid scale model diffusivity with ν_SGS as the SGS viscosity and Pr_t is a constant (see Appendix A). The characteristic length scale of the flow is chosen to be the size of the eddies, L₁₁. The numerical diffusivities are discussed in detail in Appendix B.

The last undefined term in the above definitions is F_c, the convective flux. This can be described as the covariance of the velocity and enthalpy (h) fields:

$$\begin{array}{c}\hfill F_c=⟨\overline{\rho }\stackrel{\sim }{h″{U″}_r}⟩.\end{array}$$(16)

To measure the effects of the numerical methods and the resolutions, we define the following quantities: the ratio of resolved to SGS kinetic energies (q_SGS), the ratio of the numerical to SGS viscosities (q_num) and the ratio of the numerical energy diffusivity to SGS energy diffusivity (q_diff):

$$\begin{array}{c}\hfill q_{\mathrm{SGS}}=\frac{e_{\mathrm{SGS}}}{e_{t,0}},q_{\mathrm{num}}=\frac{{\nu }_{\mathrm{num}}}{{\nu }_{\mathrm{SGS}}},q_{\mathrm{diff}}=\frac{{\chi }_{\mathrm{E}}^{\mathrm{num}}+{\chi }_{\mathrm{T}}^{\mathrm{num}}}{{\chi }_{\mathrm{SGS}}},\end{array}$$

where χ_num is the numerical energy diffusivity (see Appendix B).

A significant problem of stellar convection modeling is that we cannot directly simulate Pr ∼ 10⁻⁷ flows, as the numerical dissipation would act as an additional viscosity and increase the effective Prandtl number. In SPHERLS, the numerical viscosity strongly depends on the model resolution and the current flow velocities (see Appendix B); therefore, it depends on location, and we cannot determine it before calculations. To circumvent this problem and describe expected flow characteristics before simulation, we define a characteristic Prandtl and Reynolds number that is used to describe expected Prandtl and Reynolds number differences via different zoning and resolution. For this, we use the minimum effective numerical viscosity needed for stable solution of hyperbolic partial differential equations (Strikwerda 2004), in the form of

$$\begin{array}{c}\hfill {\nu }_{\mathrm{eff}}=\frac{\mathrm{\Delta }{U}_{\mathrm{ch}}{l}_{\mathrm{ch}}}{2},\end{array}$$(17)

where ΔU_ch is the characteristic velocity differences on the grid scale (1.5 km/s expected at this resolution with U_ch ≈ 10 km/s characteristic velocities for RR Lyrae stars based on 1D calculations), and l_ch is the characteristic cell size of the model taken in the HI ionization zone. Since the radial grid is finer in the HI zone, we set l_ch = RΔθ. The characteristic length scales of the eddies, L_ch is chosen to be 1.5H_p, based on the 1D mixing length theory, and usual parametrization of 1D codes (Paxton et al. 2019) The radiative thermal diffusivity χ_ch, and pressure scale height H_p is calculated from 1D models at 10 000 K. So, our characteristic Prandtl number is defined as

$$\begin{array}{c}\hfill P{r}^{\mathrm{ch}}=\frac{{\nu }_{\mathrm{eff}}}{{\chi }_{\mathrm{ch}}}=\frac{\mathrm{\Delta }{U}_{\mathrm{ch}}R\mathrm{\Delta }\theta }{2{\chi }_{\mathrm{ch}}},\end{array}$$(18)

and the Reynolds number as

$$\begin{array}{c}\hfill R{e}_{\mathrm{ch}}=\frac{U_{\mathrm{ch}}{L}_c}{{\nu }_{\mathrm{eff}}}=\frac{3U_{\mathrm{ch}}{H}_p}{\mathrm{\Delta }{U}_{\mathrm{ch}}R\mathrm{\Delta }\theta }·\end{array}$$(19)

We use these quantities to get a predominant view about the expected behavior of the simulations, as in this form they only depend on Δr and Δθ which are defined by our choice of cell numbers (N_r and N_θ respectively) and the horizontal size of the domain (θ_width) in the latter.

The final goal of our study is the comparison of the 1D and 2D pulsation models; hence, it is crucial to determine how the 1D properties of the convective regions depend on the domain size and resolution. We define five length scale parameters characterizing the size of the different structural components of the convective zone in RR Lyraes; see Tables 3 and 4.

These length scale parameters are given in pressure scale height units (H_p = pr²/(ρGM_r), where p is the pressure, M_r is the mass inside the sphere with radius r and G is the gravitational constant). We use the H_p at the bottom of the HeII partial ionization zone.

The length scales are defined by the relation of the convective flux F_c and the nondimensional super-adiabatic gradient Y = −(H_p/c_p)(∂s/∂r) similarly to Käpylä (2019), where c_p is the specific heat at constant pressure and s is the entropy. The length scales Λ_HI, and Λ_HeII are defined by the buoyancy-driven convectively unstable region, that is, where F_c > 0 and Y > 0. We define the so-called Deardorff’s or counter-gradient⁹ regions by the condition of F_c > 0 and Y < 0. The overshooting zone, Λ_OV, is described by the negative convective flux, that is F_c < 0 and Y < 0. We note that since F_c → 0 asymptotically, we have an extra condition as |F_c|/|maxF_c|> 0.001. The length scale of the size of the convective region is Λ_conv = L_c/H_p, where L_c is the size of the convective region, which is defined from the bottom of the HeII ionization zone (R_bottom) to the top of the HI ionization zone (R_top): L_c = R_top − R_bottom and $u_{\mathrm{rms}}=\sqrt{u_{\theta ,\mathrm{rms}}^2+u_{r,\mathrm{rms}}^2}$. The effectively turbulent region, Λ_t , is defined by the condition e_t, 0 ≥ 10⁶ erg. This is an important length scale to compare with the 1D calculations, as the 1D turbulent convection models of radial stellar pulsations are explicitly solving equations for e_t, 0, and many other terms of these models depend on this quantity¹⁰ (Kolláth et al. 2002; Paxton et al. 2019; Kovács et al. 2023).

Due to the choice of the H_p parameter and the fact that the two buoyancy-driven convective regions are separated, the length scales have this expected relation:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_t>{\mathrm{\Lambda }}_{\mathrm{conv}}>{\mathrm{\Lambda }}_{\mathrm{OV}}\gtrsim {\mathrm{\Lambda }}_{\mathrm{HeII}}>{\mathrm{\Lambda }}_{\mathrm{HI}}.\end{array}$$

Alongside with the length scales, we also define the mean convective velocity, $u_{\mathrm{rms}}^{\mathrm{tot}}$, and the convective turnover timescale τ_conv. These values give us an overall estimate about the timescale of convective processes, and the total kinetic energy stored in the convective region, and were used for a red giant model by Viallet et al. (2013). The mean convective velocity, $u_{\mathrm{rms}}^{\mathrm{tot}}$, is defined similarly to Arnett et al. (2009),

$$\begin{array}{c}\hfill \frac{1}{2}{M}_{\mathrm{conv}}{(u_{\mathrm{rms}}^{\mathrm{tot}})}^2=E_{\mathrm{t},\mathrm{conv}},\end{array}$$(20)

where E_t, conv is the total turbulent kinetic energy in the convective shells, and M_conv is the total mass of the convective shells:

$$\begin{array}{cc}\hfill E_{\mathrm{t},\mathrm{conv}}& =2\pi {\displaystyle \underset{R_{\mathrm{bottom}}}{\overset{R_{\mathrm{top}}}{\int }}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_0+{\theta }_{\mathrm{width}}}{\int }}}\rho (r,\theta ){|\mathit{U}(r,\theta )-{\mathit{U}}_0(r)|}^2{r}^{2}sin\theta \mathrm{d}\theta \mathrm{d}r;\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill M_{\mathrm{conv}}& =2\pi {\displaystyle \underset{R_{\mathrm{bottom}}}{\overset{R_{\mathrm{top}}}{\int }}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_0+{\theta }_{\mathrm{width}}}{\int }}}\rho (r,\theta )r^{2}sin\theta \mathrm{d}\theta \mathrm{d}r.\hfill \end{array}$$(22)

This is used to define the convective turnover timescale: ${\tau }_{\mathrm{conv}}=L_c/u_{\mathrm{rms}}^{\mathrm{tot}}$.

In summary, we use different statistical and horizontal averaged quantities that describe the size of the turbulent eddies, the behavior of the velocity field, the turbulent cascade, general flow properties, and the size and behavior of the convective region. We are interested in determining whether these quantities are affected by the sizes and the resolutions of the 2D models, and if they are, then which model parameters give consistent results, and what is the best compromise between resolution and numerical performance. In the next section, we first show the fundamental computational limitations of the calculations, and then the different results of the aforementioned quantities.

4. Results

4.1. Calculation times

In this section, we present some benchmark information that are relevant to our study. The main reason behind resolution tests is the calculation times. We aim at finding a model resolution that:

• (a)

  can resolve horizontal scales that are relevant in radial pulsation physics

• (b)

  can ensure reaching full amplitude pulsations within reasonable computation times.

Therefore, here we study the evolution of time steps in the different scales and the scaling of computation times on different numbers of cores and different resolutions on two different machines.

In Fig. 3. we show the calculation time on different number of cores in two different environments. The filled squares represent calculations on the HUN-REN Cloud (Héder et al. 2022), with AMD EPYC-Rome CPU, the filled circles are points calculated on the HUN-REN CSFK Malor computer with a CPU of Intel Xeon Gold 5220R CPU @ 2.20 GHz. One can see that in most cases, after reaching twelve cores, the computational time increases, which is a phenomenon called API overhead. This problem is caused here by the OpenMPI environment, in which parallelization is solved by having the parallel threads in different processes. These processes have their own memory and do not see into each other’s, hence, a massage interface is used to communicate between the different parallel processes. This allows efficient usage of large CPU clusters. However, in this case, due to the continuous synchronizing between cores in each time step, after (in this case) twelve parallel processes, the communication between cores takes more time than the calculation.

Fig. 3. Calculation times for 10 000 time steps, on HUN-REN Cloud (dashed line, filled squares) and the HUN-REN Malor computer (solid line, filled circles). The API overhead becomes a problem after 12 cores, while in the low vertical resolution case, this happens only for 18 cores.

Since SPHERLS uses an explicit method to update velocities (Geroux & Deupree 2011), the time steps must be kept below the Courant-criterion. This causes two problems: first, as we reach the surface, the pressure scale height drops, decreasing the sound speed and thus the time step; secondly, the better horizontal resolution also decreases the time step for a given velocity. These effects can decrease time steps to the level of Δt ∼ 10⁻³ seconds, which means that in the case of 10 000 time steps, one propagates the model by 10 seconds only, but to do this the calculation time amounts to 8 minutes.

4.2. Structure of the convective zone

The convective zone in the observed quasi-static RR Lyrae stellar models can be divided into four main regions based on the conditions listed in Table 3. This structural division is presented together with the components of the convective flux in Fig. 4. The zones can be followed through the relation of the super-adiabatic (or dimensionless entropy) gradient, Y (presented by a green dotted line, values corresponding to the right vertical axis), and the convective flux, F_c (presented by the blue curve, values belonging to the left vertical axis). The identified regions are the convectively unstable zone at the second helium ionization zone (we denote it with HeII), the unstable zone in the overlapping hydrogen and first helium partial ionization zones (HI), a counter-gradient layer (CG) which connects the two unstable zones, and a large overshooting region (OV) penetrating into the lower layers of the star. Above HI, one can also see a thin counter-gradient and overshooting region, which are very thin compared to the regions below.

Fig. 4. Structure of the convective zone and the flux components. The left-hand side is the vicinity of the HeII zone, and the right-hand side is the vicinity of the HI zone. The left-hand side vertical axis is the flux, and the horizontal axis is the radius. The green dotted curve with the right-hand axis denotes the dimensionless entropy gradient. Comparing the sign of this and the convective flux (blue curve) gives the borders of the different regions, which are separated by the dashed gray vertical lines. The condition of the regions is given in Section 3, and also in Table 3. The other flux components are the temperature flux (orange curve), the pressure flux (red curve), and the ionization energy flux, Fμ (yellow dash-dotted curve). The sum of these three components gives the black dashed curve. The kinetic energy flux is shown as a purple line.

In the animated version of Fig. 2, available online, one can also see plumes in the CG region, originating from the HI region and penetrating into the HeII region. The HeII region has larger rotating structures, and blobs from these can break off, escaping in either direction. Those in the CG usually connect into the HI region, while in the OV they are dissolved.

The sizes of the different layers for every model are given in Table 5. We plot these and $u_{\mathrm{rms}}^{\mathrm{tot}}$ against the θ_width in Fig. 5. The color code denotes the characteristic Prandtl numbers of the different models. Models of v046 are denoted by plus signs, and models with radial cell numbers of 122 are denoted by crosses; v0366200x20 is presented by the gray triangle. For $u_{\mathrm{rms}}^{\mathrm{tot}}$, Λ_conv, Λ_HI and Λ_HeII, there is a decreasing trend with the horizontal width, while Λ_OV and Λ_t increases when we increase the θ_width. While in the θ_width < 10° region, the resolution has larger effects than in the larger θ_width regime.

Fig. 5. Different properties of convection vs. the horizontal size of the model. These data are also presented in Table 5. The color-coding denotes the resolution through the characteristic Prandtl number. The standard vertical resolution v036 models are denoted by circles, and the reduced vertical resolution (Nr = 122) models by crosses. Pluses denote the v046 runs, and the v036_6_200x20 model is represented by the gray triangle. Top left panel: Average convective velocity; top right panel: Size of the turbulent region (Λt) in pressure scale height (Hp) units. Middle left: Size of the total convective region in Hp units. Middle right: Size of the hydrogen ionization zone in Hp units. Bottom left: Size of the overshooting region in Hp units. Bottom right: Size of the second helium ionization zone in Hp units.

The consistency and amplitude of differences are also different for the different Λs. The largest change happens with Λ_OV, changing by a factor of 3. This means that θ_width > 6° is a crucial condition for having meaningful simulations in the HeII region. The models with low resolution (Pr_ch > 0.3) show very low velocities, but otherwise, they fit in the trend in the case of the other parameters. v046 models show the same trends as v036 models, while the vertical resolution has a small effect.

Interestingly, we did not find a counter-gradient region between the overshooting and HeII zones, not even in the higher resolution models. In low-resolution models, the CG zone can be “squeezed out”, meanwhile convection is already insufficient in the HeII region, causing an insufficient transport of flux downward. On the bottom-line, this feature needs further investigation in quasi-static RR Lyrae models in the future.

4.3. Components of the convective flux

As our final goal is the comparison with 1D models, we study the different decompositions of the convective flux, that are used to derive 1D models. The energy transported by the convective motions is typically the fluctuating enthalpy flux F_c (see Eq. (16) and also Kupka & Muthsam 2017). Regularly, the two main methods using thermodynamic equities to split this term into its components are (Viallet et al. 2013)

$$\begin{array}{cc}\hfill h\prime & =\overline{T}s\prime +\overline{\left(\frac{1}{\rho }\right)}p\prime +\overline{{\frac{\partial h}{\partial \mu }|}_{s,p}}\mu \prime ,\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill h\prime & =\overline{c_p}T\prime +\overline{\left(\frac{1-\delta }{\rho }\right)}p\prime +\overline{{\frac{\partial h}{\partial \mu }|}_{T,p}}\mu \prime ,\hfill \end{array}$$(24)

where δ = (ρ/T)(∂lnρ/∂lnT)_p is the isobaric compressibility, and μ is the mean molecular weight. In most cases, the pressure terms are neglected based on filtering out acoustic waves (this is called low Mach number approximation Gough 1977; Stellingwerf 1982; Kuhfuss 1986). Eq. (23) is usually used in the descriptions of the current state-of-the-art pulsation codes (Bono & Stellingwerf 1994; Yecko et al. 1998; Paxton et al. 2019). Viallet et al. (2013) showed that in the case of Eq. (23) the pressure term is, in fact, significant in the stellar envelope; meanwhile, in Eq. (24), the compressibility mostly compensates the pressure term (δ ≈ 1). This latter equation is used many times for analysis (Geroux & Deupree 2013), and for modeling, as well (Canuto & Dubovikov 1998).

The last term of Eqs. (23) and (24) are also neglected in the stellar envelopes, usually because the mean molecular weight changes only because of the ionization, which can be incorporated into c_p and δ or s. On the other hand, in the case of the sharp stratification of the ionized zones in RR Lyrae stars, the ionized matter that moves upward by the convective motions carries a significant amount of latent heat: their ionization energy. We found this in our study to be a significant term, giving up to 50% of the convective flux in the v036 models. This is mostly significant in the HI region.

We present this effect in Fig. 4. The two panels show the vicinity of the two dynamically unstable regions: the left-hand side is the HeII zone, and the right-hand side is the HI zone. The left vertical axis is the flux, while on the right vertical axis, we present the nondimensional super-adiabatic gradient Y with a green dotted curve, and the horizontal axis is the radial coordinate inside the star. The real enthalpy flux, F_c, is denoted by the blue line. We separate the different regions described in the previous section by vertical dashed black lines. We also show the different terms of the convective energy flux: the temperature component $F_T=\overline{\rho }\overline{c_p}\stackrel{\sim }{T″{U″}_r}$ is presented by an orange curve, the pressure term $F_p=\overline{\rho }\overline{(1-\delta )/\rho }\stackrel{\sim }{p″{U″}_r}$ is given by the red solid line, and the molecular weight term, $F_{\mu }=\overline{\rho }\overline{{(\partial h/\partial \mu )}_{T,p}}\stackrel{\sim }{\mu ″{U″}_r}$ is denoted by the dash-dotted yellow line. This last term is calculated only approximately in our study, and the sum of the aforementioned terms (black dashed curve) is not fully equivalent with F_c. Lastly, the purple curve denotes the turbulent kinetic energy flux defined by (Viallet et al. 2013) $F_{\mathit{kin}}=\overline{\rho }\stackrel{\sim }{{e″}_{t,0}{U″}_r}$. This term does not take part in the enthalpy transport, and it is sometimes neglected (see, e.g., parameter set A and B in Paxton et al. 2019), but as we can see in Fig. 4, it is more significant than F_p. This flux term appears near the edges of the total convective zone (i.e., not in the inter-zone counter-gradient layer). Interestingly in the top of the HI zone, F_kin becomes positive, which phenomenon is not present in solar and red giant convection, but was found earlier in A-type stars, which have thin convection zones similarly to the models presented here (Kupka et al. 2009).

4.4. Structure and anisotropy of the velocity field

Now, we turn our attention to the horizontal fluctuations of the velocity fields. In Fig. 6. we present the velocity fluctuations ($u_{r,\mathrm{rms}}$ with solid and $u_{\theta ,\mathrm{rms}}$ with dashed line) versus the radius of four different models of v036: v0366150x20 (top left), v03612150x20 (top right), v0366150x50 (bottom left) and v03612150x40 (bottom right). These four models cover the observed trends. The rms velocities show two peaks, a smaller one in the HeII ionization zone at R ∼ 3.2 × 10¹¹ cm, and a larger one in the overlapping HeI and HI ionization zone at R ∼ 3.4 × 10¹¹ cm. The slope of the radial velocity shows a hump at roughly the middle of the HeI region. The figure shows the presence of the Rayleigh-Bénard cells: we can see that in the middle of the convection zones, one has higher radial than horizontal velocities.

Fig. 6. Velocity fluctuations in different models. The solid line denotes the radial direction, and the dashed lines denote the horizontal direction. Top left: v036_6_150x20; top right: v036_12_150x20; bottom left: v036_6_150x50; bottom right: v036_12_150x40.

This structure also mirrors the structure seen in Fig. 4. The two dynamically unstable regions generate their own convective shells. Still, in the inter-zonal region, these layers overshoots toward each other, establishing a loose connection between the two zones.

The anisotropy of the flow follows the same pattern in every model. Generally, the convective zones are anisotropic with σ_rθ ∼ 2, while in the counter-gradient and overshoot layers, it is almost isotropic. Mostly, we can say that the counter-gradient and overshooting layers show isotropic decaying of the granulation structure, as shown in Fig. 7. In the case of bad resolution Δθ > 0.2°, we experience a drop in the horizontal velocities and an increase in the radial velocities. In the case of resolutions with Δθ > 1°, numerical viscosity is so large that the overall convective velocities drop. Compared to the reference case of θ_width = 6°, N_θ = 20, increasing the size of the computational domain and increasing the resolution has a similar effect on $u_{\theta ,\mathrm{rms}}$: it is increased, especially near the surface.

Fig. 7. Anisotropy factors, σrθ, of four different v036 models vs. the radial coordinate. The models are the reference model, v036_6_150_20 (blue solid line); v036_12_150x20 (orange dashed line); v036_6_150x50 (green dotted line); v036_12_150x40 (red dash-dotted line).

4.5. General properties of the flow structure

The dimensionless numbers describing the properties of the flow are presented against the Pr_ch and θ_width in Fig. 8, and are listed in Table 2. We also present a list of grid aspect ratios, horizontal Reynolds numbers, and effective Prandtl and Péclet numbers at different radial positions in Appendix D. From the characteristic Prandtl and Reynolds number, one can conclude that the resolution presented here is too low to show actual turbulence in the simulations. The effective Prandtl numbers in the middle of the HI zone (Pr_HI) are in the same range for every model: Pr_HI ∈ [0.02; 0.06] (top left panel). This is systematically lower than Pr_ch by a factor of 2–10, which is a combined effect of our arbitrarily chosen χ_ch, the fact that viscosity ν ∝ U, and that $|\mathit{U}|<U_{\mathrm{ch}}$¹¹. This can be seen especially in the case of the very low-resolution models, with θ_width = 45°. The ratio of numerical to SGS diffusion is in the same range, weakly anti-correlating with the opening angle (top right panel), and similar behavior can be observed in the case of the ratio of numerical to SGS viscosities (middle left panel). Both quantities are extremely large (q_diff ∼ 𝒪(1000); q_num ∼ 𝒪(100)) which also indicates very large numerical damping, which easily outweighs the SGS model. Therefore, the ratio of SGS energy to total kinetic energy is low: q_SGS ≈ 10⁻³ (middle right panel), and it is slightly increased in the bad resolution models, which is expected.

Fig. 8. Dimensionless properties of the simulations vs. the angle width. The color-coding denotes the characteristic Prandtl numbers. The symbols are the same as in Figure 5. The corresponding data can be found in Table 2. Upper left: Prandtl number in the middle of the HI zone; Upper right: Ratio of numerical to SGS diffusivity; Middle left: Ratio of numerical viscosity to SGS viscosity; Middle right: Ratio of the SGS kinetic energy to the turbulent kinetic energy of the resolved scales; Lower left: Maximum Nusselt number; Lower right: Péclet number in the middle of the HI zone.

In the v036 models, the convective flux is inefficient and strongly depends on the opening angle (bottom left panel); meanwhile, Péclet numbers are similar in the different models: Pe ≈ 0.75, except for the low-resolution model. The v046 models have larger convective fluxes (Nu > 1.1), but still smaller than one would expect from 1D models (Kovács et al. 2023).

It can be said altogether that the numerical properties and Nusselt number (similarly to $u_{\mathrm{rms}}^{\mathrm{tot}}$, see Fig. 5) anti-correlate with the opening angle of the models. At the same time, the other dimensionless quantities do not depend significantly on them. The resolution also has a much smaller effect, and it decreases Pr_HI (which is expected) while increasing $u_{\mathrm{rms}}^{\mathrm{tot}}$. Meanwhile q_diff ≫ 1 and q_num ≫ 1 meaning that even though q_SGS ≪ 1, this means in this case that the velocity is so much damped by the numerics on the grid scale, that the SGS model is just a second-order problem next to it. This means that q_SGS cannot be used to determine the quality of the resolution.

Comparing these results with the calculations of v046, we can see that as it has a lower effective temperature, it also shows the same trends with the opening angle. Meanwhile, having increased (v0366200x20) or decreased vertical resolution does not produce significantly different results.

4.6. Length scales of the turbulent convection

As shown in Fig. 2 and the supplementary animation (available online), the horizontally larger models do not contain more vorticity structures in the HeII region; instead, the size of the structures is increased. To further delve into this problem, we calculated the longitudinal and transversal integral length scales of the velocity field, as they give us information about the size of the largest eddies. In this section, we refer to the diameters of the eddies as 2L₁₁ and 2L₂₂. We note that we choose the reference models to Δθ = 0.3° and θ_witdh = 6°, that was used in Geroux & Deupree (2013).

We present the results of seven chosen models in Fig. 9 to differentiate between the effects of different opening angles and resolutions. This figure shows 2L₁₁ and 2L₂₂ values in the top and bottom panels, respectively. We chose the models for the presentation to cover the problem of resolution and opening angles. Hence, we compare the reference models to different θ_width runs in the left-hand side panels, while the effects of resolution change compared to the reference model are presented in the right-hand side panels.

Fig. 9. Longitudinal (top panels) and transversal (bottom panels) integral length scales of different realizations as a function of the radial coordinate. Left panels: Comparison of different θwidth models; Right panels: Comparison of θwidth = 6° models with different resolutions. The vertical dashed lines separate the different regions of the convection zones: the dynamically stable zone, the overshooting region, the convection zone corresponding to the HeII ionization layer, the counter gradient zone that connects the HeI ionization and the HeII ionization layers, and the convective zone belonging to the overlapping HeI and HI ionization. The thin dotted lines denote the horizontal length of a cell (bottom) and the total horizontal domain in the 6° wide model (top).

In Figure 9 top left panel we compare the 2L₁₁ values of the model v0364150x20 (blue solid line), v0366150x20 (yellow dashed line), v0369150x30 (green dotted line) and v03612150x40 (red dash-dotted line). The horizontal axis is the radial coordinate, and we labeled the different subregions of the convective zone. We plotted the pressure scale height with a thinner purple dashed line for comparison. The resolution of the Δθ = 0.3° models is shown with a black dotted line at the bottom of the figure, while the total domain size of the θ_width = 6° model in the different depths is denoted by the top black dotted line. We can see that the horizontal velocities, given a large enough domain, create eddies at the bottom and top of the HeII zone around 2/3 of the reference (θ_width = 6°) models. In other words, the HeII zone’s boundary regions produce eddies with sizes correlating with the computational domain size.

In the meantime, there is no significant change in the longitudinal integral length scales when we increase the resolution, as it can be seen in the top right panel of Fig. 9. Here, the reference model (v0366150x20) is the blue solid line, v0366150x30 is denoted by the yellow dashed line, a green dotted line is the v0366150x40 and red dash-dotted line presents the higher resolution v0366150x50. As one can see, there are no significant differences between the longitudinal integral length scales of these models.

The angle width effect is present, albeit with much less efficiency, when we compare the transversal lengths scale in the bottom left panel of Fig. 9; the models are the same as those in the top left panel. We note that in homogenous isotropic turbulence (Pope 2000), L₂₂ = L₁₁/2. Comparing the left-hand side panels to each other, and also with the anisotropies in Fig. 7, one can deduce that the flow is indeed almost isotropic in the counter-gradient and overshooting regions. Corresponding to the larger anisotropy in the dynamically unstable regions, the 2L₂₂ values are also larger.

The increased resolution did not affect the 2L₂₂ values, similarly to the case of 2L₁₁. On the other hand, we detected a strange anomaly in the case of v0366150x30 and v0366150x40, in which the transversal integral length scales are enhanced at the top boundary of the HeII zone. We show this phenomenon in the bottom right panel of Fig. 9. This behavior is puzzling and will need further investigation.

Another interesting feature of the length scales is the return to the horizontal size of the model in the stable region. As $|\mathit{U}|\to 0$, the R_rr(r_d)/R_rr(0)→1, the transversal integral length scale quickly becomes the total horizontal width. Meanwhile, it seems that there is some remaining horizontal velocity in the deeper regions, which still shows a decaying auto-correlation. In the deeper 2D regions, 2L₂₂(r) = rθ_width while 2L₁₁(r)≈rθ_width/3.

In summary, we can conclude that the diameter of the larger eddies is greatly dependent on the angular size of the modeling box until reaching θ_width > 9°. After that point, the effect diminishes. The larger resolution than the reference Δθ = 0.3° does not affect the larger eddies, and there is some remainder horizontal velocity in the dynamically stable region.

4.7. The behavior of the energy cascade

An LES simulation aims to resolve the large, possibly anisotropic scales of motion, that is, the energy-containing ranges (Pope 2000) of the turbulent cascade. The filtering function that appears in the derivation of such a method (Garnier et al. 2009) damps the energy cascade at the filtering scales (in our case, it is l, the grid scale length). The energy of these scales and (those that are below l) are modeled through the e_SGS models. Having reached ∼80% of the kinetic energy in the resolved scales implies that the largest structures are resolved (Pope 2000), but these will be severely damped by the numerical (and SGS) viscosities; therefore, a real turbulent energy cascade cannot be formed. For a good resolution, one would need at least 95–99% in the resolved TKE, which also means that q_num ≪ 1 and q_SGS ≪ 1. As we have seen in Section 4.5, our simulations are far from this.

Nevertheless, numerical simulations of this scale can show an energy cascade that eventually have a power law subrange. This does not mean that we reached automatically the inertial subrange. In fact, the granulation pattern of stars (slow uprising and strong confined downflows) show a similar energy spectrum (Nordlund et al. 1997).

To study the presence of granulation, we present the normalized horizontal kinetic energy spectra (Φ₁₁/(2e₁₁)) and normalized turbulent kinetic energy spectra (e_t(κ)/e_t, 0) of fourrepresentative models in Fig. 10 (top and bottom panels, respectively). We chose two different depths for the comparison: the middle of the HeII zone and the middle of the HI zone; the orange curves represent the former, while the blue curves denote the latter. The dashing type of the curves denotes the models: results of v0366150x20 is denoted by solid curves, v0366150x50 is presented by the dotted curve, dashed curve presents v03615150x40 and lastly the dash-dotted curve is the v03645150x20.

Fig. 10. Logarithm of the horizontal kinetic energy spectrum normalized by 2e11 (top panel) and the total normalized turbulent kinetic energy spectrum normalized by et, 0 (bottom panel) vs. the logarithm of the horizontal wave-number in the HI zone (blue) and the HeII zone (orange) for different models: v036_6_150x20 (solid line); v036_6_150x40 (dotted line); v036_12_150x60 (dashed line); v036_45_150x20 (dash-dotted line). We show the Kolmogorov −5/3 slope as a dashed black line for comparison. The low-resolution model does not show the granulation energy cascade as it is damped already at the energy maximum. Simultaneously, the amount of total kinetic energy is lower in the HeII region, and L11 and L22 are larger. Thus, the spectrum has κmax at a lower wavenumber.

Taking the strong damping into account, we accept a resolution (and horizontal size) if the structure of the energy spectrum is the following:

• There is a driving range from the largest scales (smallest κ), in which the spectrum increases, reaching its maximum.

• This maximum (κ_max) corresponds to the largest scale eddies, and the energy is determined by the number of up and downflows that can be created.

• At κ > κ_max, e(κ) starts decreasing, which range contains a part where e(κ)∝κ^−5/3 or even κ⁻² (Nordlund et al. 1997).

• This is followed by a steeper decay, where numerical/SGS viscosity takes over.

In the top panel of Fig. 10, we can see that in the bad-resolution models, the large-scale part of the spectrum is present, but it is sharply cut at the Nyquist-frequency, without reaching the decaying range. In addition, the spectrum fluctuates, and HI and HeII spectra are the same. The other models show that the real granulation spectra starts at wave numbers κ_max ∼ 3 × 10⁻³ km⁻¹. The driving range start at lower wave numbers and are not present in the 6° wide models (lack of rising part). Otherwise, the larger resolution model resembles more closely the spectra presented by Nordlund et al. (1997).

If we look at the total turbulent kinetic energy spectrum in the bottom panel of Fig. 10, we can see that the vertical kinetic energy has different spectra in the two different regions. The resolution features are similarly present, though we can see that the dissipation range of the HI zone is not that strong with models Δθ = 0.3. Meanwhile, in the HeII zone, κ_max ∼ 1.7 × 10⁻³ km⁻¹, that is, motions here are on a larger scale, which can be explained by the larger pressure scale height in this region.

5. Discussion

We organize our discussion around the three main results: (1) the convective zone structure, (2) the domain size and resolution effects, and (3) the optimal size of RR Lyrae models based on our results.

5.1. The consequences of the convective zone structure

As we have seen in Sections 4.2 and 4.3, the convective zone structure differs from the 1D picture, such as provided by the hydro-codes of Yecko et al. (1998) or Smolec & Moskalik (2008a). Recently, these two codes were compared by Kovács et al. (2023). Primarily, the two convectively unstable regions at the first hydrogen ionization zone and the second helium ionization zone are not completely separated in the static case (as in the code of Smolec & Moskalik 2008a), but they are connected through decaying eddies, in a counter-gradient region.

The presence of counter-gradient layers is commonly known in astrophysics (Kupka & Muthsam 2017), but their actual presence in a flow can depend on various parameters. For instance, Käpylä (2019) found a correlation between the extent of this zone and the effective Prandtl number in solar-like convection simulations. Meanwhile, we do not find a CG region at the bottom of the HeII zone, not even in higher resolution models, suggesting that this region is very narrow, and we couldn’t resolve it. In the RGB model of Viallet et al. (2013), the transition between buoyancy driven and overshooting region is also very sharp. We conclude that this region would need an infeasibly large vertical resolution.

Nevertheless, the presence of CG layers means that the general mixing length view currently employed in pulsation theory is not a good approximation in these regions. This approach assumes that F_c ∝ Y, that is, the convective flux is aligned with the super-adiabatic temperature gradient. The usage of this type of gradient-diffusion hypothesis creates various problems: overestimated convective velocities (Wuchterl & Feuchtinger 1998), unphysical stratification (Kupka et al. 2022; Kovács et al. 2023), missing turbulent kinetic energy in the overshooting regions (Kovács et al. 2023). This theoretical problem has effects on the nonlinear behavior of pulsating stars (see, e.g., Smolec & Moskalik 2008b), and probably is one of the main reasons for the reproduction difficulties of secondary light curve features (Marconi 2017).

Our other main result is that due to the sharp stratification of the hydrogen partial ionization zone, the temperature flux cannot reproduce the whole convective flux: $\overline{\rho }\overline{c_P}\stackrel{\sim }{U″T″}\not\approx \overline{\rho }\stackrel{\sim }{U″h″}=F_c$ (see Fig. 4). Instead, we have latent heat transport, which originates from the ionization energy of the transported HII and HeII ions. This energy transport is significant, giving about half of the convective flux. In the meantime, Viallet et al. (2013) showed that the pressure term is nonnegligible if we use the entropy flux only (see the derivation of Kuhfuss 1986). This means that convection models that describe convective flux through solving temperature fluctuation (e.g., Canuto & Dubovikov 1998) have to deal with this problem, as well. One way of doing this is to incorporate the fluctuations in the specific heat (in which ionization is usually incorporated) into the definition of F_c. This method was investigated by Montgomery & Kupka (2004), who found it to be a correction of 15% in A-type stars as well as in DA and DB white dwarfs. We note that this is a specific feature of sharp ionization fronts, and this problem probably does not occur in solar or red giant models. On the other hand, nonlocality of the convective-flux is a more serious problem in 1D models.

Another feature we found is the sign inversion of kinetic energy flux at the top of the HI zone. This phenomenon is not present in solar or red giant convection, instead it was reported in A-type star convection by Kupka et al. (2009). Kupka et al. (2009) gives an explanation by pointing out that in their model the convective region is very shallow, causing heating from the bottom and cooling from the top acting on the same distance around one pressure scale height. Although they have presented this in a model without helium, our model has lower density, and an even sharper temperature gradient.

The structure of the velocity field is similar to the red giant and other models, similar to those in Viallet et al. (2013), but the two convective zones can be separated as two maxima in the rms velocity profile (Fig. 6). In general, the large-scale convective motions are anisotropic (see Fig. 7), they become almost isotropic in the inter-zone counter-gradient layer and the overshooting region as well, which is supported by the behavior of integral length scales of this region (L₁₁ ≈ 2L₂₂). The question of how this anisotropy influences the mode selection through pressure-like terms needs further study.

Altogether, the overall structure of the convective zones in static RR Lyrae models is composed of the two distinct dynamically unstable regions (also the engines of the pulsation),connected through a counter-gradient layer. This creates a very complicated coupling between these zones, which needs further studies for both the static and the dynamical models. The first step in unfolding this problem is determining the optimal model size and resolution to allow affordable computations with enough details. In the next section, we summarize and discuss our experiences with these models’ size and resolution constraints.

5.2. Effects of the domain size and resolution

Nearly all the important quantities (size of the different structures in the convective zones, velocities, effectiveness of the convection, sizes of the eddies, anisotropy) showed some dependence on the horizontal size and resolution of the models.

In the case of the structure of the convective zone, the resolution effects become less important at larger models (see Sect. 4.2 and Fig. 5). The horizontal width of the models had larger effects, and universally, we see that the mean convective velocity, $u_{\mathrm{rms}}^{\mathrm{tot}}$ decreases with the horizontal width of the models (θ_width) while the length scale of the overshooting zone, Λ_OV increases. This, together with the shape of the energy spectrum (i.e., κ_max ∼ 1.7 × 10⁻³ km⁻¹) and the fact that the eddy sizes (Fig. 9) became unaffected by the angular sizes, show us that in the lower θ_width models, the “box” of our simulation is too small, the largest two eddies interact with themselves through the periodic boundary condition. We note that with periodic boundary conditions, one always has at least two counter-rotating eddies (Glatzmaier 2014). As we presented in Sect. 4.6, the horizontal size severely affects the length scale of the overshooting (Λ_OV) and nonzero turbulence (Λ_t) regions. The size of these regions are connected to the eddies formed in the HeII region. If θ_width is too small, the corresponding eddies will be small as well, resulting in a lower penetration depth.

The change in the other parameters can also be connected to the eddy sizes. The smaller eddies must increase the velocities to cover up the kinetic energy budget missing from the larger scales. This enhances convective flux in these models.

Meanwhile, we have a stronger resolution problem in the horizontal direction, which increases numerical viscosity, outweighing the SGS model. This enhanced convection (especially in the HeII zone) will increase the Nusselt number. Due to the technical details of SPHERLS (Geroux & Deupree 2011), the numerical viscosity is also dependent on the velocity, which will also be enhanced in this case (see Appendix B). Interestingly, we experience that the background pulsation is more effectively excited with smaller θ_width, which is probably caused by the stronger instabilities of the flow. The supplementary animation also showcases the more violent nature (see also Fig. 2) , which can be responsible for the increase of total turbulent kinetic energy of the v066150x40 model seen in Fig. 1.

The v046 models follow the same trends as the v036 models too, indicating that this might be a general feature of the RR Lyrae LES models. On the other hand, it is important to note that despite being too small, the reference models that were also used in Geroux & Deupree (2013, 2014), Δθ = 0.3° ; θ_width = 6° models reproduced M3 RR Lyrae light curves reasonably well, leading us to the conclusion that the effects are probably small compared to the classical pulsation instabilities themselves. This is due to the fact that classical pulsation itself is governed by the kappa-mechanism which needs good radial resolution, while 2D simulation even with low resolution can produce more accurate convective effects then 1D mixing length theory. These box sizes and resolutions are nevertheless inadequate to model strongly nonradial phenomena, such as p-mode pulsation (see, e.g., the resolution of the models of Stein & Nordlund 2001). Meanwhile, studying the nonlinear, chaotic effects and the comparison with 1D models need the full view of which effects are numerical and which are physical; our further studies will need larger θ_width models.

The larger L₂₂ integral length scales of the models v0366150x30 and v0366150x40 in the bottom right panel of Fig. 9 remain a puzzle. The asymmetry of the resolution is an unlikely explanation as the v0366150x50 model has a similar L₂₂ profile, than v0366150x20. We hypothesize that some quasi-stationary structures are developed in these models. This can be a feature of the 2D modeling as well, since, 2D convection models are mostly less ergodic than the 3D ones (Kupka 2020). Studying stochastic and quasi-static structures would be an interesting future project, as these can cause modulations, for example, but are beyond of the scope of the present work.

The behavior of horizontal structures that lower L₁₁ in the stable zone is probably a numerical artifact inherent to the periodic boundary condition. The original numerical fluctuations can normally cause a net drifting of the Rayleigh-Bénard cells (Glatzmaier 2014). Meanwhile, we note that the net velocities in the stable region are negligible ($u_{\theta ,\mathrm{rms}}\ll 1$ km/s).

In the case of the velocity profiles, we see that increasing horizontal resolution and increasing the horizontal size of the models both decrease the anisotropy of the models. The effect in the bad resolution models is interesting: at around Δθ ≈ 0.6°, we had enlarged radial velocities and decreased horizontal velocities at the convective boundary. In contrast, above Δθ > 1°, all kinetic energy terms are decreased. To understand this feature, we should look at the kinetic energy spectra of the models, as we can see in Fig. 10. The Δθ > 1° model cannot catch the granulation, while the θ_width = 12° , Δθ = 0.3° model shows similarities with the granulation spectra presented by Nordlund et al. (1997).

Another problem with the large models with bad resolution is the numerical damping. With q_num ≫ 1 and q_diff ≫ 1, numerical damping acts on all lengthscales. Comparing our energy spectra to the one of Chan & Sofia (1996) or Porter & Woodward (2000), we cannot reproduce the Kolmogorov cascade, and our energy spectra resembles the granulation pattern (Nordlund et al. 1997). Even with codes producing less damping, one needs at least a decade of wavenumbers, that are affected by it (Kritsuk et al. 2011). In our models wavenumbers barely exceed one decade, and therefore all scales are affected by the damping. In these circumstances, the effects of the SGS model cannot be studied.

The vertical resolution has smaller effects, as it has a constraint originating from the pulsation calculation itself, and therefore it is larger than the horizontal resolution. The increased vertical resolution model (v0366200x20 with gray triangle in Figs. 5 and 8) does not show significant differences from the standard 150 radial zone models. We ran some decreased vertical resolution models, with 122 radial zones (denoted by crosses in Figs. 5 and 8). In these cases, having more horizontal cells had a larger effect than having less vertical cells. This mainly means that a model grid with resolution adequate for pulsation modeling (see, for example, some description in Paxton et al. 2019) is adequate to resolve the convective eddies, as well.

5.3. Optimal size of 2D RR Lyrae models

The largest obstacle is the calculation time of the models. As we previously presented in Fig. 3, the SPHERLS code reaches API overhead above twelve cores. This number mostly depends on the vertical resolution, as horizontal parallelization is not implemented in the code due to calculating grid velocities (Geroux & Deupree 2013). This means that roughly above 40 horizontal cells, we reach real-time computation speed. For a classical pulsator to reach maximum amplitude, it is needed to calculate hundreds or thousands of pulsation cycles (Kolláth et al. 2002). To study mode selection, we need more (Szabó et al. 2004; Paxton et al. 2019). Therefore, a real-time calculation with an RR Lyrae, for 1000 pulsation cycles could take almost one and a half years, which is not feasible.

Our other obvious constraints are the resolution and minimum horizontal sizes needed for a meaningful simulation. We found that the minimum angular size that is needed for the largest eddies in the second helium ionization zone is at least 9°–10°, supported by the behavior of the longitudinal length scale and the appearance of the second up and downflow structures. The longitudinal integral length scale can be also used for a resolution constraint (Pope 2000), which is 2L₁₁ ≈ 1 × 10¹⁰ cm which is around Δθ ≈ 0.3°. Meanwhile, the delay in the initial convective velocity growth (see Fig. 1) indicates, that Δθ < 0.2° is needed for realistic driving. This way, we can conclude that Geroux & Deupree (2013) had an adequate resolution to observe surface eddies, but for a meaningful comparison with 1D models, we need at least 1.5–2 times larger models with slightly better resolution. The N_θ = 40 θ_width = 12° models are currently the maximum affordable models for our available hardware regarding domain size, while reaching stricter Δθ constraints one would need N_θ = 45 θ_width = 9°. An optimistic choice of N_θ = 30 θ_width = 9° models can be seen as a compromise solution for quicker calculations.

6. Summary and conclusions

We performed a series of quasi-static envelope simulations with different angular widths and resolutions of two M3 RR Lyrae stars to determine the main features of the convection zone and to acquire an optimal setup for comparison with 1D model codes. The input parameters were previously derived by Kovács et al. (2023).

To achieve our goals, we analyzed the resulting models through a series of statistical methods used in turbulence theory (Pope 2000). Under the assumption of the ergodic hypothesis, we calculated the fluctuating velocity (or root-mean-square difference) profiles $\sqrt{\stackrel{\sim }{{(U″)}^2}}$, the convective flux, and its components; we determined the size of the larger convective eddies by using longitudinal and transversal integral length scales. We also calculated the kinetic energy spectrum of the models in the two main convective zones. Our results are the following:

1. The convective zone in the studied RR Lyrae stars is complex. In the rms velocity profile, one can determine the two unstable zones at the overlapping first partial ionization regions of hydrogen and helium and the second ionization region of helium. We call these dynamically unstable regions HI and HeII, respectively. However, these two regions are connected in the velocity profiles by a counter-gradient layer.

2. This inter-zonal counter-gradient layer and the bottom overshooting region are almost isotropic. The size of the overshooting region correlates with the angular size of the model, with a difference in the region’s size with a full pressure scale length.

3. We found positive kinetic energy flux in the top of the HI zone, which indicates very shallow convection, similar toA-type stars (Kupka et al. 2009).

4. Due to the sharp stratification of the HI region, the traditional assumption that the convective flux can be described by temperature fluctuation alone does not hold. Approximately half of the energy in this zone is transported by the ionized matter as latent heat. This is because the stratification is so sharp and the ionization energies are so large that a small amount of displacement causes a large amount of energy difference.

5. The presence of a counter-gradient layer suggests that the relaxation of local diffusion based description of convection in 1D models should be a top priority in new radial pulsation codes.

6. The smaller angular size models had larger convective velocities, more efficient convective transport, and more unstable flow structure. They had larger contamination of background pulsation, as well.

7. This was caused by the enhanced convection that was an effect of the small box size, where eddies are not large enough, and the kinetic energy missing from the larger scales is compensated at the scales available in the model. The enhancements were negated at larger angular sizes. This process could also be followed by the change in the longitudinal and transversal integral length scales.

8. Based on the correlation behavior (O’Neill et al. 2004) between the integral length scales (L₁₁ and L₂₂) and horizontal model sizes (θ_width), as well as the appearance of the second up- and down-flow structures in the HeII region, we estimate a minimal necessary angular size of θ_width ≈ 9°.

9. Based on the condition of Δ_filt < =L₁₁ (Garnier et al. 2009; Pope 2000), the minimal necessary angular resolution was determined as Δθ = 0.3°, while adequate driving of the convection needs at least Δθ ≲ 0.2°, based on the initial velocity growth time. In addition, at these resolutions the energy spectra resembles the formation of granules (Nordlund et al. 1997).

10. The size difference in the eddies of the different convective zones can also be followed in the horizontal and full kinetic energy spectra. The wavenumber of maximum energy (κ_max is smaller in the HeII region, explained by the larger pressure scale heights in this region.

The models presented here show a severe numerical damping on all scales, which can limit their usefulness. Based on the previous results of Geroux & Deupree (2014, 2015), we can conclude that even with these problems, the optimal size and resolution presented here is enough to make conclusions about the large-scale convection and its representation in 1D models of classical variables.

Our current results of quasi-static RR Lyrae envelope convection pave the way toward a full-scale comparison with the 1D static envelope models, calculated with the Budapest-Florida Code (Yecko et al. 1998; Kolláth et al. 2002) and MESA RSP (Smolec & Moskalik 2008a; Paxton et al. 2019). The SPHERLS code already showed promising 2D results (Geroux & Deupree 2013, 2014, 2015). In the following papers in our series, we provide a full analysis and comparison with the 1D models in addition to the static model structure differences covering the fully developed pulsation states.

This project helps make our current knowledge more robust regarding the interaction between turbulent convection and pulsation. It also provides a valuable asset to improve our models in the future in their physical content and numerical realizations, exploiting the recently available hardware and software tools.

Data availability

Movie associated to Fig. 2 is available at https://www.aanda.org

Acknowledgments

We are grateful to the referee, Dr. Friedrich Kupka for the careful and thorough reading of our manuscript, and whose detailed suggestions greatly improved the quality of this paper. We thank Robert Deupree and Chris Geroux for their generous help with the SPHERLS code. This project has been supported by the ‘SeismoLab’ KKP-137523 Élvonal grant, OTKA projects K-129249, SNN-147362, K-147131 and NN-129075, as well as the NKFIH excellence grant TKP2021-NKTA-64. The research was supported by the EKÖP-24 University Excellence Scholarship Program of the Ministry for Culture and Innovation from the source of the National Research, Development and Innovation Fund, through the EKÖP-24-4-I-ELTE-363 grant. This research was also supported by the International Space Science Institute (ISSI) in Bern/Beijing through ISSI/ISSI-BJ International Team project ID #24-603 – “EXPANDING Universe” (EXploiting Precision AstroNomical Distance INdicators in the Gaia Universe). On behalf of Project ‘Hydrodynamical modeling of classical pulsating variables with SPHERLS’ we are grateful for the usage of HUN-REN Cloud (see Héder et al. 2022, https://science-cloud.hu/) which helped us achieve the results published in this paper.

References

Appendix A: Governing equations of SPHERLS

Here we present the governing equations of the SPHERLS code using a notation consistent with that one used in the paper. Further details are presented in the original instrument papers of Geroux & Deupree (2011, 2013). The main three governing equations are the conservation of mass, momentum, and energy:

$$\begin{array}{cc}\hfill {\mathrm{D}}_t\rho +\rho \mathrm{div}\mathit{U}& =0\text{,}\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill {\mathrm{D}}_t\mathit{U}& =-\frac{1}{\rho }\mathrm{grad}p-\frac{G{M}_r}{r^2}\widehat{\mathit{r}}+\mathrm{div}{{\displaystyle \underline{\mathit{\tau }}}}^{\mathrm{SGS}}\text{,}\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill {\mathrm{D}}_{t}e& =-\frac{p}{\rho }\mathrm{div}\mathit{U}-\frac{1}{\rho }\mathrm{div}({\mathit{F}}_r+{\mathit{F}}_{\mathrm{SGS}})+D_{\mathrm{SGS}}\text{.}\hfill \end{array}$$(A.3)

Here ${\mathrm{D}}_t=\partial /\partial t+(\mathit{U}-{\mathit{U}}_0)·\mathrm{grad}$ is the Stokes derivative, U is the velocity, U₀ is the grid velocity, p is the thermodynamic pressure, ρ is the density, G is the gravitational constant, M_r is the total mass inside radius r, $\widehat{\mathit{r}}$ is the unity vector in radial direction, ${{\displaystyle \underline{\mathit{\tau }}}}^{\mathrm{SGS}}$ is the subgrid scale turbulent stress tensor. In eq. (A.3) e denotes thermodynamic internal energy, F_SGS the subgrid scale energy flux, and D_SGS describes energy transfer (dissipation) from subgrid scales to resolved scales. The radiative energy flux, F_r, is described by the diffusion approximation,

$$\begin{array}{c}\hfill {\mathit{F}}_r=-\frac{16\sigma T^3}{3\rho {\kappa }_R}\mathrm{grad}T\text{,}\end{array}$$(A.4)

where T is the temperature, σ is the Stefan-Boltzmann constant and κ_R denotes the Rosseland mean opacity. This approximation breaks down in the near-surface layers, where the atmosphere becomes optically thin, resulting in too hot layers, but classical pulsation is driven by the kappa mechanism deeper in the envelope. Since the kinematics and their coupling with convection matters more, 1D pulsation codes (and SPHERLS) have no aim to produce accurate atmospheres suitable for applications such as synthetic spectra calculations. Instead, all of these code produce bolometric light-curves which are then transformed to photometric bands via quasi-static atmosphere models (see Bono & Stellingwerf 1994).

The grid velocity is used to maintain radial resolution, especially in pulsating models, where due to the structural differences between the initial and full-amplitude models, a simple Eulerian grid wouldn’t been able to maintain the resolution of the partial ionization zones, which drive the pulsation. This grid velocity is used to update the shell radii, while Lagrangian (mass) derivatives (i.e., $\partial /\partial r=4\pi r^2\overline{\rho }\partial /\partial M_r$) are used during the calculation in the radial direction. The grid velocity is derived from the net inflow and outflow of the total spherical shell (see for details Geroux & Deupree 2011):

$$\begin{array}{c}\hfill 0={\displaystyle \int \frac{\partial \rho }{\partial t}}\mathrm{d}V+{\displaystyle \oint \mathit{U}\rho }\mathrm{d}A-{\displaystyle \oint {\mathit{U}}_0}\rho \mathrm{d}A\text{.}\end{array}$$(A.5)

The horizontal boundary condition is periodic, so U₀ only has a radial component, and the grid radius coordinates are updated by

$$\begin{array}{c}\hfill \frac{\mathit{dr}}{\mathit{dt}}={\mathit{U}}_0\widehat{\mathit{r}}\text{.}\end{array}$$(A.6)

The subgrid scale effects are described by the Smagorinsky-model (Smagorinsky 1963), and their formulas are the following (Geroux & Deupree 2013, and references therein):

$$\begin{array}{cc}& {\nu }_{\mathrm{SGS}}=\frac{{\mathrm{C}}^2{l}^2}{\sqrt{2}}{(\mathrm{grad}\mathit{U}:[\mathrm{grad}\mathit{U}+{(\mathrm{grad}\mathit{U})}^T])}^{\frac{1}{2}}\text{,}\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill & {{\displaystyle \underline{\mathit{\tau }}}}^{\mathrm{SGS}}=\rho {\nu }_{\mathrm{SGS}}[\mathrm{grad}\mathit{U}+{(\mathrm{grad}\mathit{U})}^T-\frac{2}{3}{\displaystyle \underline{\mathit{I}}}\mathrm{div}\mathit{U}]\text{,}\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill & {\mathit{F}}_{\mathrm{SGS}}=-\frac{{\nu }_{\mathrm{SGS}}}{P{r}_{\mathrm{t}}}\mathrm{grad}e\text{,}\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill & D_{\mathrm{SGS}}=\frac{{\mathrm{d}}_{\mathrm{t}}}{L}{\left(\frac{{\nu }_{\mathrm{SGS}}}{{\mathrm{A}}_{\mathrm{t}}L}\right)}^3\text{.}\hfill \end{array}$$(A.10)

Here ${\displaystyle \underline{\mathit{I}}}$ denotes identity matrix, l is the cell scale length¹² ($l=\sqrt{\mathrm{\Delta }rr\mathrm{\Delta }\theta }$, with vertical grid spacing of Δr and horizontal grid spacing of Δθ), L = 3.75 × l, A_t = 0.117, d_t = 1.4,Pr_t = 0.7, C = 0.17 are constants. The subgrid scale turbulent kinetic energy is related to the subgrid scale kinematic viscosity by the relation: e_SGS = [ν_t/(A_tL)]².

Appendix B: Numerical viscosity and artificial diffusion

Since there is no modeled molecular viscosity in the SPHERLS models (since it affects only scales that are much smaller than the model grid), to actually make our models comparable to real stars and other models, we need to derive some other measures instead.

In numerical calculations, the other source of dissipation next to molecular (and SGS) viscosity is numeric viscosity or numerical damping, which is always unique to the numerical method. Its source is the truncation errors caused by finite differences.

These truncation errors can be calculated using Taylor series (Hirsch 2007). Since we deal with polar coordinates, one has to take into account the transform of the unit vector. Hence, we write the general Taylor series as:

$$\begin{array}{cc}\hfill f(\mathit{x}+& \mathbf{\Delta }\mathit{x})=f(\mathit{x})+\mathbf{\Delta }\mathit{x}·\mathrm{grad}f+\frac{1}{2}{\mathbf{\Delta }\mathit{x}}^T·(\mathrm{grad}\mathrm{grad}f)·\mathbf{\Delta }\mathit{x}+\hfill \\ \hfill & +\frac{1}{6}\mathbf{\Delta }\mathit{x}·[{\mathbf{\Delta }\mathit{x}}^T·(\mathrm{grad}\mathrm{grad}\mathrm{grad}f)·\mathbf{\Delta }\mathit{x}]+{\mathcal{O}(|\mathbf{\Delta }\mathit{x}|}^4)\text{.}\hfill \end{array}$$(B.1)

The displacement vectors are $\mathbf{\Delta }\mathit{x}=\mathrm{\Delta }r\widehat{\mathit{r}}$ in the radial and $\mathbf{\Delta }\mathit{x}=r\mathrm{\Delta }\theta \widehat{\mathit{\theta }}$ in the horizontal direction (with $\widehat{\mathit{\theta }}$ as the covariant unit vector in the horizontal direction).

SPHERLS uses different discretization techniques in different situations. In general, it uses a staggered mesh, with fluxes defined at cell interfaces, internal variables defined at cell centers (Deupree 1977a; Geroux & Deupree 2011). Simple derivatives are calculated by the second-order accurate central method, while for the advection terms, the donor-cell algorithm is used with weighted upwind and central schemes. The upwind scheme is the CIR scheme (Courant et al. 1952). The scheme is in 1D (with positive velocity direction),

$$\begin{array}{c}\hfill \frac{\mathrm{d}f}{\mathrm{d}x}\equiv (1-\vartheta )\frac{f_{i+1}-f_{i-1}}{2\mathrm{\Delta }x}+\vartheta \frac{f_i-f_{i-1}}{\mathrm{\Delta }x}-\frac{\vartheta \mathrm{\Delta }x}{2}\frac{d^{2}f}{d{x}^2}+\mathcal{O}(h^3)\text{,}\end{array}$$(B.2)

where x denotes one of the coordinates, and the last term on the right-hand side is the second-order truncation error. The parameter ϑ is a weighting factor; its value is dependent on the velocities on the grid:¹³

$$\begin{array}{c}\hfill \vartheta =\{\begin{array}{ccc}1& \mathrm{if}& |\mathit{U}|/c_{\mathrm{s}}\ge 1,\\ |\mathit{U}|/c_{\mathrm{s}}& \mathrm{if}& 1>|\mathit{U}|/c_{\mathrm{s}}\ge 0.2,\\ 0.2& \mathrm{if}& |\mathit{U}|/c_{\mathrm{s}}\le 0.2.\end{array}\end{array}$$(B.3)

The second-order truncation error is equivalent with an additional nonlinear diffusion term (Pope 2000), which is in 2D Cartesian (assuming ϑ = 1):

$$\begin{array}{cc}\hfill {\mathrm{D}}_{t}f& ={\partial }_{t}f+U_x{\partial }_{x}f+U_y{\partial }_{y}f\hfill \\ \hfill & \approx {\partial }_{t}f+U_x{\delta }_{x}f+U_y{\delta }_{y}f-\frac{U_x\mathrm{\Delta }x}{2}{\partial }_x^{2}f-\frac{U_y\mathrm{\Delta }y}{2}{\partial }_y^{2}f\hfill \\ \hfill & ={\partial }_{t}f+U_x{\delta }_{x}f+U_y{\delta }_{y}f-\mathrm{div}{\mathit{F}}_{\mathbf{num}}.\hfill \end{array}$$(B.4)

Here δ denotes the finite difference operator describing the algorithm. The diffusivity is linear and χ_f^num = UΔx/2 in 1D (Pope 2000). However, in more than one dimension it becomes nonlinear and direction dependent. This nonlinearity is also enhanced by the unit vector transformations in polar coordinates.

Based on this, we first derive the second-order truncation errors in each direction for the velocity and energy advection terms. The velocity damping is called numerical viscosity. The damping terms can be derived from Eq. (B.1), and have the following forms:

$$\begin{array}{cc}\hfill {\mathcal{D}}_r^U& =\frac{(U_r-U_0)\mathrm{\Delta }r}{2}\frac{{\partial }^2{U}_r}{\partial r^2}+\frac{U_{\theta }\mathrm{\Delta }\theta }{2}\frac{{\partial }^2{U}_{\theta }}{\partial r\partial \theta }\hfill \\ \hfill & -\frac{U_{\theta }\mathrm{\Delta }\theta }{2r}\frac{\partial U_{\theta }}{\partial \theta }+\frac{U_{\theta }\mathrm{\Delta }\theta }{2}\frac{\partial U_r}{\partial r}-\frac{U_{\theta }{U}_r\mathrm{\Delta }\theta }{2r}\text{,}\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\theta }^U& =\frac{(U_r-U_0)\mathrm{\Delta }r}{2r}\frac{{\partial }^2{U}_r}{\partial r\partial \theta }+\frac{U_{\theta }\mathrm{\Delta }\theta }{2r}\frac{{\partial }^2{U}_{\theta }}{\partial {\theta }^2}+\frac{(U_r-U_0)\mathrm{\Delta }r}{2r^2}\frac{\partial U_r}{\partial \theta }\hfill \\ \hfill & +\frac{U_{\theta }\mathrm{\Delta }\theta }{2}\frac{\partial U_{\theta }}{\partial r}-\frac{(U_r-U_0)\mathrm{\Delta }r}{2r}\frac{\partial U_{\theta }}{\partial r}+\frac{U_{\theta }^2\mathrm{\Delta }\theta }{2r}\text{,}\hfill \end{array}$$(B.6)

$$\begin{array}{c}\hfill {\mathcal{D}}_e=\frac{(U_r-U_0)\mathrm{\Delta }r}{2}\frac{{\partial }^{2}e}{\partial r^2}+\frac{U_{\theta }\mathrm{\Delta }\theta }{2r}(\frac{{\partial }^{2}e}{\partial {\theta }^2}+\frac{\partial e}{\partial r})\text{.}\end{array}$$(B.7)

Here 𝒟_r^U, 𝒟_θ^U, and 𝒟_e are the damping of U_r, U_θ, and e, respectively. Assuming Δr/r ≪ 1 and rΔθ/r ≪ 1 we can neglect the last term in (B.5) and the last two terms in (B.6), arriving at

$$\begin{array}{c}\hfill D_t\mathit{U}\approx \partial \mathit{U}+\mathit{U}\mathit{\delta }\mathit{U}-({\mathit{\nu }}_{\text{num}}\mathrm{grad}){(\mathrm{grad}\mathit{U})}^T,\end{array}$$(B.8)

where the direction dependent numerical viscosity is

$$\begin{array}{c}\hfill {\mathit{\nu }}_{\text{num}}=(\frac{U_r\mathrm{\Delta }r}{2},\frac{U_{\theta }r\mathrm{\Delta }\theta }{2}).\end{array}$$(B.9)

Then we take numerical viscosity as ${\nu }_{\text{num}}=|{\mathit{\nu }}_{\text{num}}|$. Similarly one can derive from (B.7), that 𝒟_e = (ν_numgrad )grad e ≈ χ_e^numdivgrad e with the numerical energy diffusivity χ_e^num = ν_num.

Since the radiative heat conductivity is very large in our models, it is worth to study the numerical dissipation of the heat diffusion term divF_r. This is calculated by the second-order accurate scheme (Hirsch 2007),

$$\begin{array}{cc}& \mathrm{div}{(K\mathrm{grad}T)}_{i,j}=\hfill \\ \hfill & \frac{1}{\mathrm{\Delta }{x}^2}[K_{i+1/2,j}(T_{i+1,j}-T_{i,j})-K_{i-1/2,j}(T_{i,j}-T_{i-1,j})]\hfill \\ \hfill & +\frac{1}{\mathrm{\Delta }{y}^2}[K_{i,j+1/2}(T_{i,j+1}-T_{i,j})+K_{i,j-1/2}(T_{i,j}-T_{i,j-1})]\hfill \\ \hfill & +\mathcal{O}(\mathrm{\Delta }{x}^2,\mathrm{\Delta }{y}^2),\hfill \end{array}$$(B.10)

where K = χρc_p, and the radiative conductivity, and half-index values are K_{i + 1/2, j} = (K_i + 1 + K_i)/2. We note that equation (B.10) is in Cartesian coordinates for simplicity. Truncation error here involves the third derivatives of the flux (i.e., fourth derivatives of temperature). Since this is a higher-order error than in the case of advection, we simplify the problem by considering F_r components as scalars. Therefore, the truncation error becomes:

$$\begin{array}{cc}\hfill {\mathcal{D}}_T& =\frac{\mathrm{\Delta }{r}^2}{6}\frac{{\partial }^{3}K{(\mathrm{grad}T)}_r}{\partial r^3}+\frac{\mathrm{\Delta }{\theta }^2}{6r}\frac{{\partial }^{3}K{(\mathrm{grad}T)}_{\theta }}{\partial {\theta }^3}\hfill \\ \hfill & +\frac{\mathrm{\Delta }{\theta }^2}{2}\frac{{\partial }^{2}K{(\mathrm{grad}T)}_{\theta }}{\partial r\partial \theta }-\frac{\mathrm{\Delta }{\theta }^2}{6r}\frac{\partial K{(\mathrm{grad}T)}_{\theta }}{\partial \theta },\hfill \end{array}$$(B.11)

from this we approximate numerical heat diffusivity as

$$\begin{array}{c}\hfill {\chi }_T^{\mathrm{num}}\approx \frac{{\mathcal{D}}_T}{\mathrm{div}\mathrm{grad}T}\text{.}\end{array}$$(B.12)

This way, we provided all necessary numerical values. We note that these are usually only good approximations, as truncation errors are highly nonlinear.

On the other hand, to provide stability during shocks, SPHERLS also includes an artificial viscosity term (Geroux & Deupree 2011, and references therein) as an additional pressure term in equations (A.2) and (A.3). However, we found this term is negligible in smooth regions.

Appendix C: Numerical testing

The presence of the ionization front of RR Lyrae and other classical variables may question the usage of simple numerical methods that are used in SPHERLS. To investigate this problem, we perform a numerical test of the donor-cell method against the Sod shock tube problem (Sod 1978) on an 1D Eulerian grid, solving Euler equations. The equation is given in the conservative form:

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\mathit{q}+\frac{\partial }{\partial x}\mathit{F}=0\text{,}\end{array}$$(C.1)

where:

$$\begin{array}{c}\hfill \mathit{q}=\left[\begin{array}{c}\rho \\ \rho v\\ \frac{p}{\gamma -1}+\frac{1}{2}\rho v^2\end{array}\right],\mathit{F}=\left[\begin{array}{c}\rho v\\ \rho v^2+p\\ (\frac{\gamma }{\gamma -1}p+\frac{1}{2}\rho v^2)v\end{array}\right].\end{array}$$(C.2)

Here ρ is the density, v is the velocity, p is the pressure, and γ = 1.4 is the adiabatic constant. The equation in the test is solved explicitly by this scheme:¹⁴

$$\begin{array}{c}\hfill q_i^{\ast ,n}=q_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}(F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}}),\end{array}$$(C.3)

$$\begin{array}{cc}& q_i^{n+1}=q_i^{\ast ,n}\hfill \\ \hfill & +\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\xi [|q_{i+1}^{\ast ,n}-q_i^{\ast ,n}|(q_{i+1}^{\ast ,n}-q_i^{\ast ,n})-|q_i^{\ast ,n}-q_{i-1}^{\ast ,n}|(q_i^{\ast ,n}-q_{i-1}^{\ast ,n})].\hfill \end{array}$$(C.4)

The second equation is the artificial viscosity in diffusion form (Sod 1978). The half-centered flux terms are given by B.2.

We present our result of a test¹⁵ on the Sod shock tube problem in Figure C.1. The Sod shock tube is characterized with a discontinuity with

$$\begin{array}{c}\hfill {\rho }_l=1,p_l=1,v_l=0;{\rho }_r=1/8,p_r=0.1,v_r=0.1\text{.}\end{array}$$(C.5)

After the simulation is started, a shockwave starts to move to the right, which is followed by a slower contact discontinuity, while a rarefaction fan propagates to the left. This problem is a common test for numerical hydrocodes.

We performed one test with good resolution (Δx = 0.02), and three tests on a sparser grid (Δx = 0.12). All four tests were run to t = 1.4 s, the time steps were determined by the Courant–Friedrichs–Lewy condition, with different CFL numbers (0.1 for the dense grid, and 0.5, 0.25, and 0.1 in the three sparse cases, respectively), such that

$$\begin{array}{c}\hfill CFL=\frac{max(|v_i|+c_i)\mathrm{\Delta }t}{\mathrm{\Delta }x}\text{.}\end{array}$$(C.6)

In Fig. C.1, we present the density, velocity, pressure, and entropy of the system (top to bottom). In the case of the dense grid, the scheme is highly unstable, as we needed to use a strong artificial viscosity with ξ = 2.5 to achieve not disastrous result, but the oscillations are still visible, especially for the velocity and entropy. The sparse grids are more stable with ξ = 1, we only experienced oscillations with CFL = 0.5, while over and undershooting was present. The rarefaction fan is only slightly diffused, while the discontinuities show diffusion as well. Overall, this scheme shows similar features than the Ricthmyer scheme, but with stronger overshooting and oscillations.

Fig. C.1. Results of the Sod shock-tube test for four scenarios: a higher resolution (Δx = 0.02, first column), and three sparser grids (Δx = 0.12) with three different numbers: 0.5, 0.25, and 0.1. The rows represent the density (ρ), velocity (v), pressure (p), and specific entropy (s) (from top to bottom).

Other numerical tests on an adaptive grid were made by Buchler et al. (1997) against the Noh (Noh 1987) and Sedov (Sedov 1959) problems, while in SPHERLS, the 2D Sedov-problem was investigated by Geroux & Deupree (2011).

It is worth to mention that in 1D pulsation codes (e.g., Stellingwerf 1975; Bono & Stellingwerf 1994; Smolec & Moskalik 2008a) traditionally Lagrangian grids are used with simple upwind based flux terms. The sole exception is the adaptive code of the Budapest-Florida code (Buchler et al. 1997; Yecko et al. 1998; Kolláth et al. 2002), where the van Leer (1977) scheme is used. Their numerical tests showed that the simple CIR scheme is inadequate to properly follow shock fronts (Buchler et al. 1997). This is an interesting result, considering the successes of 1D pulsation modeling (see Kovács et al. 2023, Sec. 2. and references therein). We propose a few answers to resolve this contradiction:

• The ionization front is only investigated in the stellar envelope, while shock fronts are created above the photosphere, where the speed of sound drops. Therefore, the ionization front is mathematically not a discontinuity, and it can be resolved.

• All of these codes (including SPHERLS) use Lagrangian grids in which they resolve the ionization front, and its steep temperature gradient.

• The energy equation is solved implicitly, mostly using the Stellingwerf (1975) method, or based on the method of Fraley (1968).

At a bottom line, SPHERLS can be considered as a pulsation hydrocode that is expanded into the horizontal direction. We can expect it to give better results for the convective processes, which happens naturally here, even if (especially horizontal) numerical dissipation is much stronger than in state-of-the-art LES/VLES/iLES codes.

Appendix D: Flow and grid characteristics

In this appendix, we provide the aspect ratios and Reynolds numbers of the models (Table D.1.), as well as the Prandtl and Péclet numbers (Table D.2.) in the two convective regions and in the overshooting region. We note that our models have very inefficient convection, which means that radiative diffusion largely exceeds the numerical (and SGS) diffusion caused by calculation of the advective fluxes. This feature represents a flow very different from those at the bottom of the solar convection zone or above a convective core of intermediate to high-mass stars, resulting in small Prandtl and Péclet numbers.

All Tables

All Figures

Fig. 1. Time evolution of the total turbulent kinetic energy (TKE). Left panel: Comparison of the TKE evolution for different opening angles with resolution of Δθ = 0.3°. The models are v036_6_150x20 (blue solid line), v036_9_150x30 (yellow dashed line), and v036_12_150x40 (green dotted line). Right panel: Same as left panel, but with the same opening angle (θwith = 6°) and different resolutions: v036_6_150x20 (blue solid line), v036_6_150x30 (red dash-dotted line), and v036_6_150x40 (purple dotted line). The vertical dashed line denotes the end of τrelax and the start of τstat.

In the text

Fig. 2. Snapshot of five different simulations, taken at the same simulation time. The 2D simulations are in polar coordinates; the radial coordinate is in 1011 cm units; and the angular coordinates are in degrees. The color-coding denotes the vorticity of the fluid flow The simulation domains are rotated to fit on the figure. Left panel: Models v036_6_150x20 (top), v036_9_150x30 (middle), and v036_12_150x40 (bottom). These models have different angle widths, but the same resolution. Right panel: v036_6_150x20 (top), v036_6_150x30 (middle), v036_6_150x40 (bottom). These models have the same angle width, but different horizontal resolutions. In the left panel we can see that with increasing opening angle, the number of structures is also increasing in the HI region (at around 3.4 × 1011 cm radius.). In the HeII zone (at around 3.2 × 1011 cm), the size of the vortices increase. This means that the initial model with of 6° degrees used by Geroux & Deupree (2013) is inadequate to properly capture the properties of the flow structure in this area. Meanwhile, we can see on the right-hand side that the number of structures in the HI region is similar with higher resolutions (although more details appear). An animated version of this figure is available online.

In the text

	Fig. 3. Calculation times for 10 000 time steps, on HUN-REN Cloud (dashed line, filled squares) and the HUN-REN Malor computer (solid line, filled circles). The API overhead becomes a problem after 12 cores, while in the low vertical resolution case, this happens only for 18 cores.
In the text

Fig. 4. Structure of the convective zone and the flux components. The left-hand side is the vicinity of the HeII zone, and the right-hand side is the vicinity of the HI zone. The left-hand side vertical axis is the flux, and the horizontal axis is the radius. The green dotted curve with the right-hand axis denotes the dimensionless entropy gradient. Comparing the sign of this and the convective flux (blue curve) gives the borders of the different regions, which are separated by the dashed gray vertical lines. The condition of the regions is given in Section 3, and also in Table 3. The other flux components are the temperature flux (orange curve), the pressure flux (red curve), and the ionization energy flux, Fμ (yellow dash-dotted curve). The sum of these three components gives the black dashed curve. The kinetic energy flux is shown as a purple line.

In the text

Fig. 5. Different properties of convection vs. the horizontal size of the model. These data are also presented in Table 5. The color-coding denotes the resolution through the characteristic Prandtl number. The standard vertical resolution v036 models are denoted by circles, and the reduced vertical resolution (Nr = 122) models by crosses. Pluses denote the v046 runs, and the v036_6_200x20 model is represented by the gray triangle. Top left panel: Average convective velocity; top right panel: Size of the turbulent region (Λt) in pressure scale height (Hp) units. Middle left: Size of the total convective region in Hp units. Middle right: Size of the hydrogen ionization zone in Hp units. Bottom left: Size of the overshooting region in Hp units. Bottom right: Size of the second helium ionization zone in Hp units.

In the text

	Fig. 6. Velocity fluctuations in different models. The solid line denotes the radial direction, and the dashed lines denote the horizontal direction. Top left: v036_6_150x20; top right: v036_12_150x20; bottom left: v036_6_150x50; bottom right: v036_12_150x40.
In the text

	Fig. 7. Anisotropy factors, σrθ, of four different v036 models vs. the radial coordinate. The models are the reference model, v036_6_150_20 (blue solid line); v036_12_150x20 (orange dashed line); v036_6_150x50 (green dotted line); v036_12_150x40 (red dash-dotted line).
In the text

Fig. 8. Dimensionless properties of the simulations vs. the angle width. The color-coding denotes the characteristic Prandtl numbers. The symbols are the same as in Figure 5. The corresponding data can be found in Table 2. Upper left: Prandtl number in the middle of the HI zone; Upper right: Ratio of numerical to SGS diffusivity; Middle left: Ratio of numerical viscosity to SGS viscosity; Middle right: Ratio of the SGS kinetic energy to the turbulent kinetic energy of the resolved scales; Lower left: Maximum Nusselt number; Lower right: Péclet number in the middle of the HI zone.

In the text

Fig. 9. Longitudinal (top panels) and transversal (bottom panels) integral length scales of different realizations as a function of the radial coordinate. Left panels: Comparison of different θwidth models; Right panels: Comparison of θwidth = 6° models with different resolutions. The vertical dashed lines separate the different regions of the convection zones: the dynamically stable zone, the overshooting region, the convection zone corresponding to the HeII ionization layer, the counter gradient zone that connects the HeI ionization and the HeII ionization layers, and the convective zone belonging to the overlapping HeI and HI ionization. The thin dotted lines denote the horizontal length of a cell (bottom) and the total horizontal domain in the 6° wide model (top).

In the text

Fig. 10. Logarithm of the horizontal kinetic energy spectrum normalized by 2e11 (top panel) and the total normalized turbulent kinetic energy spectrum normalized by et, 0 (bottom panel) vs. the logarithm of the horizontal wave-number in the HI zone (blue) and the HeII zone (orange) for different models: v036_6_150x20 (solid line); v036_6_150x40 (dotted line); v036_12_150x60 (dashed line); v036_45_150x20 (dash-dotted line). We show the Kolmogorov −5/3 slope as a dashed black line for comparison. The low-resolution model does not show the granulation energy cascade as it is damped already at the energy maximum. Simultaneously, the amount of total kinetic energy is lower in the HeII region, and L11 and L22 are larger. Thus, the spectrum has κmax at a lower wavenumber.

In the text

	Fig. C.1. Results of the Sod shock-tube test for four scenarios: a higher resolution (Δx = 0.02, first column), and three sparser grids (Δx = 0.12) with three different numbers: 0.5, 0.25, and 0.1. The rows represent the density (ρ), velocity (v), pressure (p), and specific entropy (s) (from top to bottom).
In the text