© The Authors 2026

1. Introduction

Modeling turbulent convection inside stellar cores and envelopes is one of the most important problems in astrophysics. While most studies have focused on phenomenological improvements to the mixing-length theory of Böhm-Vitense (1958) (Kupka & Muthsam 2017), in the case of radial stellar pulsations, time-dependent convection models were necessary (Houdek & Dupret 2015) to describe the significant effects of the narrow convective zones in the envelope of these stars. In these objects, envelope convection is confined to partial ionization zones, which also drive radial pulsations via the κ–γ mechanism. Despite the envelope being mostly radiative, these zones provide major damping for the pulsations.

The most straightforward way of modeling stellar convection is to use 3D numerical hydrocodes. Geroux & Deupree (2011, 2013) developed the multidimensional stellar pulsation code SPHERLS, which was the first to successfully model RR Lyrae stars at maximum amplitudes in 2D (Geroux & Deupree 2013, 2014) and 3D (Geroux & Deupree 2015), although this was achieved using relatively simple numerical methods¹. Mundprecht et al. (2013) extended the ANTARES (Muthsam et al. 2010) code to model radial stellar pulsations; they compared the 2D structure to 1D treatments for a Cepheid model (Mundprecht et al. 2015). Vasilyev et al.; Vasilyev et al. (2017; 2018, and references therein) used the COBOLD5 software to model a 2D Cepheid variable and then used it as input for spectral synthesis models. Recently, Stuck et al. (2025) performed a 2D model analysis of the convective region in Cepheid models using the MUSIC code, while Kovács et al. (2025) performed a similar study of RR Lyrae stars using the SPHERLS code. However, both studies investigated these stars in the static limit.

All of the aforementioned studies either focus on a special case or calculate only a few stellar models due to the extreme computational cost of these simulations. Multidimensional simulations of these stars remain challenging, because adaptive grids are required to keep the partial ionization regions (which drive the pulsation) resolved as the star undergoes rapid expansion and contraction. Meanwhile, to reach maximum amplitudes, one must simulate several (thousands) pulsation cycles. This makes improving 1D convection theories an important task (Kupka & Muthsam 2017; Kovács et al. 2023), as these models can be used in large model surveys (e.g., Smolec 2016) and also serve as inputs for multidimensional simulations.

Traditionally, asteroseismology has been studied primarily using linear approximations (Aerts et al. 2010), which have also provided important information on radial pulsators (Cox 1980). While linear theory can provide accurate pulsation frequencies, pulsation amplitudes enter the equations only as a normalization factor. Linear pulsation analysis can be performed in the adiabatic approximation, in which case the theory provides only frequencies. Linear nonadiabatic calculations also provide information about the growth rates of pulsation modes (i.e., whether a mode is pulsationally unstable) and result in somewhat different frequencies. Meanwhile, one can only approximately describe coupling between modes through amplitude equations. Nonlinear calculations, on the other hand, include every amplitude-dependent effect present in the system (and more accurate coupling with convection). This modifies the linear frequencies and provides additional results comparable to observations, such as light curves, amplitudes, and the actual nonlinearly selected pulsational modes. This is particularly important, because we often observe monomode pulsators where multiple modes are linearly excited. Therefore, mode selection studies investigate the nonlinear behavior of pulsation modes to determine which pulsation mode is excited, which survives if multiple modes are present, or whether they can coexist in double-mode pulsation. This, along with other nonlinear phenomena, has been central to radial pulsation modeling over the past two decades (see, e.g., Kolláth et al. 2002; Szabó et al. 2004; Smolec 2016). Nevertheless, both linear and nonlinear studies have shown that including turbulent convection is crucial for describing radial stellar pulsations (Baker 1987).

Building on the pioneering work of Gough (1965), Unno (1967), and Castor (1968), the nonlocal models² of Stellingwerf (1982a), Kuhfuss (1986), and Gehmeyr & Winkler (1992) are used in current state-of-the-art nonlinear radial stellar pulsation hydrocodes developed in the 1990s and early 2000s (Bono & Stellingwerf 1994; Kolláth et al. 2002; Smolec & Moskalik 2008a; Paxton et al. 2019). For a detailed summary, see Kovács et al. (2023).

Meanwhile, these models have their own limitations: they have many free parameters; nonlinear phenomena such as mode selection and double-mode pulsations are model-dependent (Szabó et al. 2004; Smolec & Moskalik 2008b; Kovács et al. 2023, 2024); light and radial velocity curves show discrepancies with observations (Marconi 2017; Kovács et al. 2023) and the models produce unphysical stratification in the stellar envelope (Canuto & Dubovikov 1998; Kupka et al. 2022; Kovács et al. 2023).

Many of these problems stem from the so-called downgradient approximation of the convective flux and the resulting nonphysical stratification of turbulent kinetic energy (TKE) (Canuto & Dubovikov 1998; Kupka et al. 2022; Kovács et al. 2023). In the downgradient approximation, the convective flux is considered proportional to the super-adiabatic gradient. There is a further assumption that in stable regions (i.e., where the gradient is negative), the gradient is taken as zero (Stellingwerf 1982a). This approximation also affects the TKE profile, which in turn influences mode selection through the turbulent pressure and viscous stress terms. This effect determined whether simulations reproduce or fail to exhibit double-mode pulsation (Smolec & Moskalik 2008a,b; Kovács et al. 2023). Neglecting negative convective flux leads to insufficient turbulent damping, causing overshooting into the deep stellar interiors. This provides additional damping that can stabilize double-mode pulsations (Smolec & Moskalik 2008b; Kovács et al. 2023). Meanwhile, maintaining a negative flux prevents TKE leakage to deeper regions but leads to unstable double-mode pulsations, where one mode is eventually damped despite a “physically more correct” picture. This occurs because negative flux terms cause strong turbulent damping in the overshooting region, where convective elements behave as if colliding with a rigid wall (Kupka et al. 2022; Kovács et al. 2023). Both pictures are nonphysical, illustrating the limitations of one-equation models. While the aforementioned problems depend primarily on the coupling between convection and radial pulsation, other observed features – such as the Blazhko effect (Blažko 1907) and additional modes (the so-called f_x modes; see, e.g., Szabó et al. 2015) – can be explained by resonances with non-radial modes (see, e.g., Van Hoolst et al. 1998; Dziembowski & Mizerski 2004; Dziembowski 2016). Although such resonances require multidimensional treatment, improving the 1D description cannot be ruled out as a means to explain these features.

The analysis of multidimensional models by Mundprecht et al. (2013, 2015) and Kovács et al. (2025) reveals that a natural progression is to incorporate nonlocal time-dependent differential equations describing the convective flux and temperature fluctuations. Developed by Xiong (1980, 2021), Xiong et al. (1998b) and Kuhfuß (1987), these so-called three-equation models were incorporated into stellar evolution codes (see Flaskamp 2002, 2003; Weiss & Schlattl 2008) and used in linear non-radial pulsation theory (Xiong et al. 1998a, 2000, 2016, 2018, and references therein). However, further adjustments to these models are necessary, since convection in classical pulsators substantially differ from standard astrophysical convection:

• Highly dynamical regime. The convection timescale is comparable with the cycle of the pulsation, making a time-dependent description important, while the large structural changes during pulsation create more complex coupling between pulsation and convection.

• Anisotropy effects. Convective eddies convert kinetic energy into TKE through shear. This eddy viscosity is crucial in damping convection (Gonczi 1981; Stellingwerf 1982a), and it arises from the anisotropy of convective motions (Gehmeyr & Winkler 1992). Consequently, modeling anisotropies is also important (Mundprecht et al. 2015).

• Sharp ionization regions, especially in the case of RR Lyraes. Most of these stars exhibit a sharp HI ionization region. Small perturbations in ionization levels have a drastic effect on specific heat and opacity, which must be taken into account (Kovács et al. 2025).

We introduce an extended version of the Kuhfuss three-equation model in this paper (Kuhfuß 1987)³, which provides a basis for improving 1D nonlinear stellar pulsation codes.

2. Extended Kuhfuss three-equation model

We provide the basic equations of the extended Kuhfuß (1987) model below, while a derivation consistent with our notation is given in Appendix A. The equations are the following:

$$\begin{array}{cc}\hfill {\partial }_t\rho +\frac{1}{r^2}{\partial }_r(r^2\rho v)& =0,\hfill \end{array}$$(1) $$\begin{array}{cc}\hfill {\mathrm{d}}_{t}v& =-\frac{1}{\rho }{\partial }_r(p+p_t)-U_{\nu }-g,\hfill \end{array}$$(2) $$\begin{array}{cc}\hfill {\mathrm{d}}_{t}e+p\frac{1}{r^2}{\partial }_r(r^{2}v)& =-\frac{1}{\rho r^2}{\partial }_r\left[r^2(F_{\mathrm{r}}+F_{\mathrm{c}})\right]-\frac{2T\mathrm{\Phi }}{c_p{\tau }_{\kappa }}-S+\epsilon ,\hfill \end{array}$$(3) $$\begin{array}{cc}\hfill {\mathrm{d}}_t\tilde{\omega }+p_t\frac{1}{r^2}{\partial }_r(r^{2}v)& =-\frac{1}{\rho r^2}{\partial }_r\left(r^2{\mathcal{F}}_{\tilde{\omega }}\right)+E_{\nu }+S-\epsilon ,\hfill \end{array}$$(4) $$\begin{array}{cc}\hfill {\mathrm{d}}_t\mathrm{\Phi }& =-\frac{1}{\rho r^2}{\partial }_r\left(r^2{\mathcal{F}}_{\mathrm{\Phi }}\right)+{\mathcal{P}}_{\mathrm{\Phi }}-{ϵ}_{\mathrm{\Phi }}-\frac{2\mathrm{\Phi }}{{\tau }_{\kappa }},\hfill \end{array}$$(5) $$\begin{array}{cc}\hfill {\mathrm{d}}_t\mathrm{\Pi }& =-\frac{1}{\rho r^2}{\partial }_r\left(r^2{\mathcal{F}}_{\mathrm{\Pi }}\right)+{\mathcal{P}}_{\mathrm{\Pi }}+S_{\mathrm{\Pi }}-{ϵ}_{\mathrm{\Pi }}-\frac{\mathrm{\Pi }}{{\tau }_{\kappa }}.\hfill \end{array}$$(6)

Here, d_t = ∂_t + v∂_r denotes the Stokes derivative, and ∂_r is the spatial derivative relative to the r, radius coordinate⁴. The thermodynamic quantities follow standard notation: p is pressure, ρ is the density, T is the temperature, e is specific internal energy, and c_p is specific heat at constant pressure. The variable v describes the velocity of a given spherical shell, while g = Gm/r² denotes the gravitational acceleration, where m is the mass inside a sphere of radius r within the star. F_r is the radiative energy flux in the diffusion approximation:

$$\begin{array}{c}\hfill F_r=K_r{\partial }_{r}T=\frac{16{\sigma }_{\mathrm{B}}{T}^3}{3\rho \kappa }{\partial }_{r}T,\end{array}$$

where κ(T, ρ) is the Rosseland-mean opacity and σ_B is the Stefan-Boltzman constant. The remaining quantities relate to turbulent convection and its coupling to the dynamics. The specific TKE is denoted by $\stackrel{\sim }{\omega }$, the entropy fluctuation is denoted by Φ, and Π denotes the velocity entropy covariance.

The convective flux, F_c, is represented by a second-order description:

$$\begin{array}{c}\hfill F_c=T\rho \mathrm{\Pi }+{\alpha }_c\frac{c_{\mathit{pT}}T}{c_p}{\mathcal{F}}_{\mathrm{\Phi }},\end{array}$$(7)

where c_pT = (∂lnc_p/∂lnT)_p and α_c is a dimensionless factor of order of unity.

The term $p_t=(2/3)\rho \stackrel{\sim }{\omega }$ denotes turbulent pressure, while terms E_ν and U_ν denote the turbulent viscous energy and momentum transfer rates (Wuchterl & Feuchtinger 1998):

$$\begin{array}{cc}\hfill U_{\nu }& =\frac{1}{\rho }{\partial }_r\left[(\xi -\frac{1}{3})2\tilde{\omega }\rho \right]+(\xi -\frac{1}{3})\frac{6\tilde{\omega }}{r},\hfill \end{array}$$(8) $$\begin{array}{cc}\hfill E_{\nu }& =2\tilde{\omega }[(\xi -\frac{1}{3}){\partial }_{r}v+2(\frac{1}{3}-\xi )\frac{v}{r}].\hfill \end{array}$$(9)

Here, ξ denotes the flow anisotropy (see Appendix A). In Eqs. (3)–(6), we refer to buoyancy-driven terms as source terms denoted by S, while shear-related terms are referred to as production terms (Pope 2000) denoted by 𝒫:

$$\begin{array}{cccc}\hfill S& =-\frac{\delta }{\rho c_p}\mathrm{\Pi }{\partial }_{r}p,\hfill & \hfill S_{\mathrm{\Pi }}& =-\frac{\delta }{\rho c_p}\mathrm{\Phi }{\partial }_{r}p,\hfill \end{array}$$(10) $$\begin{array}{cccc}\hfill {\mathcal{P}}_{\mathrm{\Phi }}& =-\mathrm{\Pi }{\partial }_{r}s,\hfill & \hfill {\mathcal{P}}_{\mathrm{\Pi }}& =-2\xi \tilde{\omega }{\partial }_{r}s,\hfill \end{array}$$(11)

where ∂_rs is the specific entropy.

The dissipation terms are modeled locally through the dissipation and radiative timescales. First, the dissipation of TKE (ε) is modeled using the convective timescale. This is constructed from a characteristic length (i.e., mixing length) and characteristic velocity (approximated as $\sqrt{2\omega }$). This yields ${\tau }_{\mathrm{conv}}=\mathrm{\Lambda }/{(2\stackrel{\sim }{\omega })}^{-1/2}$, $\epsilon =\stackrel{\sim }{\omega }/{\tau }_{\mathrm{conv}}$, and

$$\begin{array}{c}\hfill \epsilon ={\alpha }_d\frac{{\tilde{\omega }}^{3/2}}{\mathrm{\Lambda }},\end{array}$$(12)

where the factor $\sqrt{2}$ is incorporated into the dimensionless constant α_d of order unity. The dissipation of the entropy fluctuations (ϵ_Φ) and velocity-entropy covariance (ϵ_Π) has two components: a radiative ($\propto {\tau }_{\mathrm{r}}^{-1}$) and a viscous part (∝Λ⁻¹):

$$\begin{array}{cccc}\hfill {ϵ}_{\mathrm{\Phi }}& =\frac{2\mathrm{\Phi }}{{\tau }_r}+{\alpha }_{d,\mathrm{\Phi }}\frac{2\mathrm{\Phi }{\tilde{\omega }}^{1/2}}{\mathrm{\Lambda }},\hfill & \hfill {ϵ}_{\mathrm{\Pi }}& =\frac{\mathrm{\Pi }}{{\tau }_r}+{\alpha }_{d,\mathrm{\Pi }}\frac{\mathrm{\Pi }{\tilde{\omega }}^{1/2}}{\mathrm{\Lambda }}.\hfill \end{array}$$(13)

The radiation timescale can be determined from the Péclet number, since $Pe=\mathrm{\Lambda }\sqrt{2\stackrel{\sim }{\omega }}\rho c_p/K_r={\tau }_r/{\tau }_{\mathrm{conv}}$. Thus,

$$\begin{array}{c}\hfill {\tau }_r=Pe\times {\tau }_{\mathrm{conv}}={\alpha }_r^{-1}\frac{{\rho }^2{\mathrm{\Lambda }}^2{c}_p\kappa }{{\sigma }_{\mathrm{B}}{T}^3},\end{array}$$(14)

where α_r = 3/16, is a dimensionless constant⁵.

The final timescale introduced is the opacity-induced dissipation-production timescale, which provides additional dissipation or production of entropy fluctuations. It is defined by τ_κ = τ_r/(3 − κ_T), where κ_T = (∂lnκ/∂lnT)_p.

Finally, the third-order moments (TOMs) describe the turbulent flux of turbulent quantities and are approximated by the down-gradient approximation⁶ (Castor 1968; Xiong 1980):

$$\begin{array}{cccccc}\hfill {\mathcal{F}}_{\tilde{\omega }}& =-{\alpha }_{\tilde{\omega }}{\mu }_t{\partial }_r\tilde{\omega },\hfill & \hfill {\mathcal{F}}_{\mathrm{\Phi }}& =-{\alpha }_{\mathrm{\Phi }}{\mu }_t{\partial }_r\mathrm{\Phi },\hfill & {\mathcal{F}}_{\mathrm{\Pi }}& =-{\alpha }_{\mathrm{\Pi }}{\mu }_t{\partial }_r\mathrm{\Pi },\end{array}$$(15)

where ${\mu }_t=\mathrm{\Lambda }\rho \sqrt{2\xi \stackrel{\sim }{\omega }}$ is the turbulent viscosity.

3. Modifications

3.1. Enhanced dissipation

The local model of viscous dissipation, ε, is known to be insufficient in the dynamically stable regions (Canuto 1993) when using the standard prescription for the mixing length Λ = Λ₀ = α_ΛH_p, where H_p = pr²/(ρGm) is the pressure scale height and α_Λ is a dimensionless parameter. This problem arises because gravity waves provide extra damping in the stable regions. Kupka et al. (2022) derived a correction to the mixing length, which we adopt here using the formula introduced of Ahlborn et al. (2022):

$$\begin{array}{c}\hfill \mathrm{\Lambda }(r)=\{\begin{array}{cc}{\displaystyle -\frac{r+{\mathrm{\Lambda }}_0}{2{\tau }^{\star }}+\sqrt{{\left[\frac{r+{\mathrm{\Lambda }}_0}{2{\tau }^{\star }}\right]}^2+\frac{{\mathrm{\Lambda }}_{0}r}{{\tau }^{\star }}}}& \text{if}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}})<0,\\ \\ {\displaystyle \frac{{\alpha }_{\beta }r}{{\mathrm{\Lambda }}_0+{\alpha }_{\beta }r}{\mathrm{\Lambda }}_0}& \text{if}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}})\ge 0,\end{array}\end{array}$$(16)

where

$$\begin{array}{c}\hfill {\tau }^{\star }={\alpha }_{\tau }{\tilde{\omega }}^{-1/2}{\mathrm{\Lambda }}_{0}g\sqrt{\rho p^{-1}({\mathrm{\nabla }}_{\mathrm{ad}}-\mathrm{\nabla })}\end{array}$$(17)

is the ratio of the dissipative to buoyancy timescales in the stable regions, and α_β and α_τ are dimensionless parameters. Here, ∇ − ∇_ad = H_pT⁻¹∂_rs is the dimensionless super-adiabatic gradient, with ∇ = (∂lnT/∂lnp) and ∇_ad = (∂lnT/∂lnp)_s.

3.2. Local anisotropy equation

We follow the traditional splitting of the Reynolds-tensor into isotropic turbulent pressure and anisotropic turbulent viscous stress tensor. The latter is usually modeled using the eddy-viscosity hypothesis (Pope 2000), which yields the anisotropy parameter

$$\begin{array}{c}\hfill \xi \approx \frac{1}{3}-\frac{{\alpha }_{\nu }\mathrm{\Lambda }}{2{\tilde{\omega }}^{1/2}}{\partial }_{r}v.\end{array}$$(18)

Kuhfuß (1987) applied the above formulation only to the turbulent viscous stress tensor, setting ξ = 1/3 elsewhere⁷. Recently Kovács et al. (2023, 2024) showed that the α_ν parameter correlates with α_d and significantly affects mode selection (Kovács et al. 2024) and nonlinear behavior (Kolláth et al. 2011). Therefore, instead of Eq. (18), we propose a local prescription derived from the nonlocal model of Canuto (1998) using the method given in Kupka et al. (2022), but retaining the shear terms (see appendix B):

$$\begin{array}{c}\hfill \xi =\frac{1}{3}+\frac{{\tau }_{\xi }}{3\tilde{\omega }}[{\mathcal{D}}_{\omega }-(2{\beta }_5+1)S+\epsilon ]-\frac{2{\tau }_{\xi }}{9}[({\alpha }_1-\frac{1}{5}){\partial }_{r}v-2\frac{v}{r}],\end{array}$$(19)

equivalent to Kupka et al.’s (2022) solution in the static case. Here, ${\mathcal{D}}_{\omega }=\sqrt{2/3}{\alpha }_t{\rho }^{-1}{r}^{-2}{\partial }_r(r^2\mathrm{\Lambda }\rho {\stackrel{\sim }{\omega }}^{1/2}{\partial }_r\stackrel{\sim }{\omega })$ describes diffusion of ξ through isotropic TKE flux, ${\tau }_{\xi }^{-1}=(2/3)[{\alpha }_1(v/r)+(2{\alpha }_1-1){\partial }_{r}v]-5\epsilon /(2\stackrel{\sim }{\omega })$, while α₁ = 1.08 and β₅ = 0.7 are constant parameters given by Canuto et al. (1994).

3.3. Ionization energy transport

The second term in Eq. (7) provides a second-order correction representing the ionization energy flux carried by ions transported via convection. In fact, this reflects a change in chemical potential, even though ionization effects are usually absorbed into c_p and κ (Kippenhahn & Weigert 1990). It becomes a significant component of the total enthalpy flux due to the sharp stratification present in RR Lyrae stars (Kovács et al. 2025). Since its primary effect is on specific heat, we parameterize it through temperature fluctuations (see Appendix C).

3.4. Opacity effects

Sharp stratification in the ionization zones also results in strong opacity bumps; therefore, opacity fluctuations cannot be neglected. Thus, we must track changes in radiative conductivity K_r, approximated by a first-order Taylor series of K_r(T, κ), and κ(T) (Xiong 1980). While the opacity bumps enhance radiative losses, they also generate further fluctuations analogous to the κ-mechanism. This is handled through the opacity-induced dissipation-production timescale, τ_κ.

3.5. Viscous dissipation

The sometimes short dynamical timescales alongside enhanced opacity (i.e., less efficient cooling) make the otherwise insignificant viscous damping of the entropy fluctuations and flux relevant. Assuming scale similarity with kinetic energy damping and the convective timescale ${\tau }_{\mathrm{conv}}=\stackrel{\sim }{\omega }/\epsilon$, the second terms in (13) follow.

4. Discussion

The model proposed here reduces to the standard three-equation model of Kuhfuß (1987) when using (18) in E_ν and U_ν, setting ξ = 1/3 elsewhere, neglecting τ_κ terms, using Λ = α_ΛH_p, and also setting the parameters α_c,α_d, Π,α_d, Φ to zero. Some of our extensions can also be incorporated into one-equation models, such as the local anisotropy formulation.

Most of the free dimensionless parameters can be adopted from the existing literature. Assuming isotropy in the dissipation regime, we can adopt the dissipation parameter of Kuhfuß (1987) calibrated to the mixing-length theory: ${\alpha }_d=8/3\sqrt{2/3}\approx 2.17$.

The diffusion parameters can be obtained from those of Kuhfuß (1987) by dividing them by $\sqrt{2\xi }=\sqrt{2/3}$. Therefore, we obtain α_Π = 6 and α_Φ = 4. In the case of α_t, we can start from the local condition of the one-equation model $\sqrt{3{\alpha }_t/(2{\alpha }_d)}\le (3/8)\sqrt{2}$, which gives α_t ≤ 4/8.

Regarding the remaining dissipation terms, Kuhfuss sets the viscous dissipation timescale of Π and Φ to ε/(c_pT) and then neglects it in the local limit (Kuhfuß 1987). This is equivalent to scaling α_d, Π ≈ α_d, Φ parameters with α_d such that ${\alpha }_{\mathrm{d},\mathrm{\Pi }}/{\alpha }_d\sim \stackrel{\sim }{\omega }/(c_{p}T)$, which is between 10⁻¹ and 10⁻² in the ionization region. Thus, we adopt α_d, Φα_d, Π = 0.0217 ≈ α_d/100.

The mixing-length parameters, α_Λ and α_β, can follow the standard values used in existing 1D codes (Kupka et al. 2022; Kovács et al. 2024): α_Λ = 1.5 and α_β = 1. The value α_τ = 0.2 is given by Kupka et al. (2022). The authors also give α_τ = 5/(32α_d). These two values match only when α_d = 0.8, derived from Kolmogorov-scaling (Canuto & Dubovikov 1998).

We chose α_r = 3/16, consistent with simple Péclet number dependence for the radiative timescale, rather than the mixing-length calibrated value of (Kuhfuß 1987) (equivalent to α_r = 48 in our notation).

Lastly, α_c = 1 effectively captures ionization energy transport effects (see Appendix C), whereas Montgomery & Kupka (2004) used α_c = 0.5 from second-order Taylor series enthalpy fluctuations.

Although our modifications improve the original Kuhfuß (1987) model for nonlinear pulsation codes, we retain a mixing length approach and a local description of turbulent dissipation and anisotropy. A more complete model (e.g., Canuto 1998) is not recommended at this stage, because a gradual improvement in the treatment help identify the key phenomena governing nonlinear pulsation. Similar to its predecessor, our model is based on the assumption of local thermodynamic equilibrium and neglects radiative transfer in optically thin regions. As these are essential for accurately modeling light curves, the next step should implement them into the models⁸.

In conclusion, implementing this model into a new hydrocode should be done gradually. Successive modifications may likely have diminishing effects: the first two have the largest effect, followed by the convection correction and opacity effects. The effect of the additional dissipation terms is likely to have an effect only in mode selection and nonlinear phenomena. We also note that the chosen parameters given here also require observational calibration.

5. Conclusions

In this paper, we presented our extended version of the Kuhfuss three-equation turbulent convection model (Kuhfuß 1987), tailored for use in the field of nonlinear radial stellar pulsations. Our main modifications include:

1. the adoption of the enhanced turbulent dissipation from Kupka et al. (2022),

2. the derivation of a local non-static anisotropy model from Canuto (1998),

3. a second-order correction to the convective flux,

4. a first-order correction accounting for opacity fluctuations. This and the flux correction represent effects of sharp ionization regions and opacity sensitivity (Kovács et al. 2025).

5. the retention of turbulent dissipation for entropy fluctuations and flux, expected to have additional stabilizing effect on nonlinear behavior.

Numerical comparison with state-of-the-art one-equation hydrocodes (as analyzed in Kovács et al. 2023, 2024) and a detailed analysis of the model approximations, based on multidimensional results, will be presented in upcoming papers.

Acknowledgments

We are grateful to the anonymous referee for the careful and thorough reading of our manuscript, and whose suggestions improved the quality of this paper. 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 project has been supported by the ‘SeismoLab’ KKP-137523 Élvonal grant, Lendület program of the Hungarian Academy of Sciences project No. LP2025-14/2025, OTKA projects K-129249, SNN-147362, K-147131 and NN-129075, as well as the NKFIH excellence grant TKP2021-NKTA-64. 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: Deriving the model

We start from the Navier-Stokes equations, of the continuity, momentum, and energy conservation in index notation. As we focus on envelope convection in a pulsating envelope, we neglect any energy production from nucleosynthesis and rotation, while we keep all time derivatives. Thus, we have:

$$\begin{array}{cc}\hfill {\partial }_t\rho +{\partial }_k(v_k\rho )& =0,\hfill \end{array}$$(A.1) $$\begin{array}{cc}\hfill {\mathrm{D}}_t{v}_j& =-\frac{1}{\rho }{\partial }_{j}p+\frac{1}{\rho }{\partial }_k{\sigma }_{\mathit{kj}}+g{\delta }_{j3},\hfill \end{array}$$(A.2) $$\begin{array}{cc}\hfill {\mathrm{D}}_t(e+\frac{1}{2}{v}_j^2)& =-\frac{1}{\rho }{\partial }_j(p{v}_j+F_j-{\sigma }_{\mathit{jk}}{v}_k)+g_j{v}_j.\hfill \end{array}$$(A.3)

Here, ρ denotes the density, v denotes the velocity, p is the thermodynamic pressure, σ is the viscous stress tensor, g = −Gm/r² is the gravitational acceleration, where m is the mass inside the spherical shell with radius r. F denotes radiative flux and e is the internal energy. The symbol ${\mathrm{D}}_t={\partial }_t+{\overline{v}}_k{\partial }_k$. Additionally, we define the heat sources as the radiative flux and viscous heat production, that is, the specific entropy (s) production is

$$\begin{array}{c}\hfill T{\mathrm{D}}_{\mathrm{t}}s=-\frac{1}{\rho }{\partial }_i{F}_i+{\sigma }_{\mathit{ij}}{\partial }_i{v}_j.\end{array}$$(A.4)

We define Reynolds-average as the horizontal average of a quantity, namely

$$\begin{array}{c}\hfill \overline{q}\frac{1}{4\pi }{\int }_0^{2\pi }{\int }_0^{\pi }q(r,\theta ,\varphi )sin\theta \mathrm{d}\theta \mathrm{d}\varphi .\end{array}$$

This operation is commutative with (Eulerian) spatial and time derivatives. Using this operation, we split quantities into mean and fluctuation parts (denoted as ( ⋅ )′), for example:

$$\begin{array}{c}\hfill v=\overline{v}+v^{\prime },p=\overline{p}+p^{\prime },\cdots \end{array}$$

The operator has the following identities:

$$\begin{array}{c}\hfill \overline{\mathit{ab}}=\overline{a}\overline{b}+\overline{a^{\prime }{b}^{\prime }},\overline{a^{\prime }}=0,\overline{\overline{a}}=\overline{a}.\end{array}$$

We denote the Stokes derivative related to the average field as (Kuhfuss 1986)

$$\begin{array}{c}\hfill {\mathrm{d}}_{t}q={\partial }_{t}q+{\overline{v}}_k{\partial }_{k}q={\mathrm{D}}_{t}q-v_k^{\prime }{\partial }_{k}q.\end{array}$$(A.5)

Now we derive equations for the mean and fluctuating quantities. First, we neglect ρ′ except buoyancy, then write Eq. (A.1) as

$$\begin{array}{cc}\hfill \partial \overline{\rho }+{\partial }_k(\overline{\rho }{\overline{v}}_k)& =0,\hfill \end{array}$$(A.6) $$\begin{array}{cc}\hfill {\partial }_k(\overline{\rho }{v}_k^{\prime })& =0.\hfill \end{array}$$(A.7)

Here, (A.7) is called anelastic approximation.

Taking the average of Eq. (A.2), we get the momentum equation of the pulsation (i.e., large-scale motion)

$$\begin{array}{c}\hfill d_t{\overline{v}}_j+\overline{v_k^{\prime }{\partial }_k{v}_j^{\prime }}=-\frac{1}{\overline{\rho }}{\partial }_j\overline{p}+\frac{1}{\overline{\rho }}{\partial }_k{\overline{\sigma }}_{\mathit{kj}}+g{\delta }_{j3}.\end{array}$$(A.8)

Here, the second term can be converted into a tensor divergence, using Eq. (A.7), one can show that $\overline{v{\prime }_k{\partial }_k{v}_j\prime }=(1/\overline{\rho }){\partial }_k\overline{\rho v_k\prime v_j\prime }\equiv -(1/\overline{\rho })\partial {\stackrel{\sim }{\sigma }}_{\mathit{kj}}$, where we define the Reynolds tensor as (Kuhfuss 1986)

$$\begin{array}{c}\hfill {\tilde{\sigma }}_{\mathit{ij}}=-\overline{\rho v_i^{\prime }{v}_j^{\prime }}.\end{array}$$(A.9)

After rearranging, Eq. (A.8).

$$\begin{array}{c}\hfill d_t{\overline{v}}_j=-\frac{1}{\overline{\rho }}{\partial }_j\overline{p}+\frac{1}{\overline{\rho }}{\partial }_k({\overline{\sigma }}_{\mathit{kj}}+{\tilde{\sigma }}_{\mathit{kj}})+g{\delta }_{j3}.\end{array}$$(A.10)

We derive the momentum equation of the fluctuating field by subtracting Eq. (A.8) from Eq. (A.2), thus

$$\begin{array}{c}\hfill d_t{v}_j^{\prime }+v_k^{\prime }{\partial }_k{\overline{v}}_j+{(v_k^{\prime }{\partial }_k{v}_j^{\prime })}^{\prime }=-{\left(\frac{1}{\rho }{\partial }_{j}p\right)}^{\prime }+\frac{1}{\overline{\rho }}{\partial }_k{\sigma }_{\mathit{kj}}^{\prime }.\end{array}$$(A.11)

Here, the first term of the right-hand side can be expanded as

$$\begin{array}{c}\hfill -\frac{{\rho }^{\prime }}{{\overline{\rho }}^2}{\partial }_j\overline{p}+\frac{1}{\overline{\rho }}{\partial }_j{p}^{\prime }\end{array}$$

, where the second term is usually neglected in the small Mach-number approximation (Gough 1977). On the other hand, Canuto (1992) gives an explanation of the importance of pressure fluctuations, while Viallet et al. (2013) showed that pressure transport can be a significant contribution to envelope convection.

Next, we derive the mean total energy equation similarly to Kuhfuss (1986). Applying averaging, the left-hand side of Eq. (A.3) becomes (Kuhfuss 1986; Kuhfuß 1987):

$$\begin{array}{cc}\hfill \overline{{\mathrm{D}}_t((e+\frac{1}{2}{v}_j^2)}& ={\mathrm{d}}_{\mathrm{t}}(\overline{e}+\frac{1}{2}{({\overline{v}}_j)}^2+\frac{1}{2}\overline{v_j^{\prime }{v}_j^{\prime }})\hfill \\ \hfill & +\overline{v_k^{\prime }{\partial }_k{e}^{\prime }}-\frac{1}{\overline{\rho }}({\tilde{\sigma }}_{\mathit{kj}}{\partial }_k{\overline{v}}_j+{\overline{v}}_j{\partial }_k{\tilde{\sigma }}_{\mathit{kj}})+\frac{1}{2}\overline{v_k^{\prime }{\partial }_k{v}_j^{\prime }{v}_j^{\prime }}\hfill \end{array}$$

Here, the Reynolds stress-related terms are derived by using Eq. (A.7) and the definition (A.9). The term $(1/2)\overline{v_j\prime v_j\prime }=\stackrel{\sim }{\omega }\equiv (1/\overline{\rho }){\stackrel{\sim }{\sigma }}_{\mathit{jj}}$ is the kinetic energy stored in the turbulence, and proved to be crucial in modeling radial stellar pulsations (Stellingwerf 1982b). Therefore, it’s necessary to close it through an exact equation. The third order moment (TOM) $(1/2)\overline{v_k\prime {\partial }_k{v}_j\prime v_j\prime }=(1/2\overline{\rho }){\partial }_k\left(\overline{\rho }\overline{v_k\prime v_j\prime v_j\prime }\right)$, where Eq. (A.7) was used, and it describes the turbulent flux of the TKE.

The RHS of the averaged (A.3) is

$$\begin{array}{c}\hfill -\frac{1}{\overline{\rho }}{\partial }_i(\overline{p}{\overline{v}}_j+\overline{p^{\prime }{v}_j^{\prime }}+{\overline{F}}_j-{\overline{\sigma }}_{\mathit{jk}}{\overline{v}}_k-\overline{{\sigma }_{\mathit{jk}}^{\prime }{v}_k^{\prime }})+g_j{\delta }_{j3}{\overline{v}}_j.\end{array}$$

Multiplying (A.10) with ${\overline{v}}_j$, one can eliminate the kinetic energy of the averaged field. The pressure fluctuation term can be put together on the right-hand side after rearranging the $\overline{v_k\prime {\partial }_{k}e\prime }$ term using Eq. (A.7), because $\overline{\rho }h\prime =\overline{\rho }e\prime +p\prime$, where h is the specific enthalpy. Using this, the average of the total energy becomes

$$\begin{array}{cc}& {\mathrm{d}}_t(\overline{e}+\tilde{\omega })=-\frac{1}{\overline{\rho }}\overline{p}{\partial }_j{\overline{v}}_j\hfill \\ \hfill & -\frac{1}{\rho }{\partial }_j(F_j+\overline{\rho }\overline{h^{\prime }{v}_j^{\prime }}+\frac{1}{2}\overline{\rho }\overline{v_j^{\prime }{v}_k^{\prime }{v}_k^{\prime }}-\overline{{\sigma }_{\mathit{jk}}^{\prime }{v}_k^{\prime }})+\frac{1}{\overline{\rho }}({\overline{\sigma }}_{\mathit{jk}}+{\tilde{\sigma }}_{\mathit{jk}}){\partial }_j{\overline{v}}_k\hfill \end{array}$$(A.12)

Neglecting viscocity transport, we need to determine ${\stackrel{\sim }{\sigma }}_{\mathit{ij}}$, $\stackrel{\sim }{\omega }$, $F_{\mathrm{c}}=\overline{\rho }\overline{h\prime v_j\prime }$, and ${\mathcal{F}}_{\stackrel{\sim }{\omega }}=(1/2)\overline{\rho }\overline{v_j\prime v_k\prime v_k\prime }$ to close Eq. (A.10) and (A.12). Since $\stackrel{\sim }{\omega }$ is half of the trace of $\stackrel{\sim }{\sigma }$, we derive the overall equation, and then derive $\stackrel{\sim }{\omega }$. First, we multiply (A.11) with v_i′, then we take the some of this equation and the one with switched indices. Lastly, we take the average. The result after rearranging is:

$$\begin{array}{cc}& d_t\overline{v_i^{\prime }{v}_j^{\prime }}=-\overline{v_k^{\prime }{\partial }_k(v_i^{\prime }{v}_j^{\prime })}-(\overline{v_i^{\prime }{v}_k^{\prime }}{\partial }_k{\overline{v}}_j+\overline{v_j^{\prime }{v}_k^{\prime }}{\partial }_k{\overline{v}}_i)-\frac{1}{\overline{\rho }}(\overline{v_i^{\prime }{\partial }_j{p}^{\prime }}+\overline{v_j^{\prime }{\partial }_i{p}^{\prime }})\hfill \\ \hfill & +\frac{1}{{\overline{\rho }}^2}(\overline{v_i^{\prime }{\rho }^{\prime }}{\partial }_j\overline{p}+\overline{v_j^{\prime }{\rho }^{\prime }}{\partial }_i\overline{p})+\frac{1}{\overline{\rho }}(\overline{v_i^{\prime }{\partial }_k{\sigma }_{\mathit{kj}}^{\prime }}+\overline{v_j^{\prime }{\partial }_k{\sigma }_{\mathit{ki}}^{\prime }}).\hfill \end{array}$$(A.13)

The first term of the right-hand side is the turbulent flux of the tensor component, the second term is the so-called shear production (Pope 2000), which converts kinetic energy into turbulence. The third term can be split into two parts as v_i′∂_jp′=∂_j(v_i′p′) − p′∂_jv_i′. The first term describes the pressure transport, and the second term is the pressure-velocity rate-of-strain tensor. The role of this latter is to restore the isotropy of the turbulence. The last two terms are the bouyancy source term, which converts internal energy into turbulent energy, and the viscous dissipation, which again is composed of a transport term and a sink term, which generates heat (v_i′∂_kσ_kj′=∂_k(v_i′σ_kj′) − σ_kj′∂_kv′_i). Let the tensor, pressure, and viscosity transport term be note by 𝒯, 𝒯^p, 𝒯^σ, the production term by 𝒫, the pressure-velocity rate-of-strain tensor by Π_ij, the source by 𝒮 and the dissipation by ε. Eq. (A.13) is then

$$d_t\overline{v_i^{\prime }{v}_j^{\prime }}=-\frac{1}{\overline{\rho }}{\partial }_k({\mathcal{T}}_{\mathit{kij}}+{\mathcal{T}}_{\mathit{kij}}^p-{\mathcal{T}}_{\mathit{kij}}^{\sigma })+{\mathcal{P}}_{\mathit{ij}}+{\mathrm{\Pi }}_{\mathit{ij}}+{\mathcal{S}}_{\mathit{ij}}-{\epsilon }_{\mathit{ij}}.$$(A.14)

Eq. (A.14) shows, that the eddy viscosity hipothesis, that ${\stackrel{\sim }{\sigma }}_{\mathit{ij}}-(2/3)\overline{\rho }\stackrel{\sim }{\omega }\propto {\partial }_i{\overline{v}}_j$ is a very crude approximation. Therefore, Canuto & Dubovikov (1998) developed a model that takes into account these terms, including rotation, and calibrated it to observational and simulation results.

The half of the trace of Eq. (A.14) gives the equation for the specific TKE, $\stackrel{\sim }{\omega }$, which is

$$\begin{array}{cc}\hfill d_t\tilde{\omega }& =-\frac{1}{\overline{\rho }}{\partial }_k(\frac{1}{2}\overline{\rho }\overline{v_k^{\prime }{v}_j^{\prime }{v}_j^{\prime }}+\overline{v_j^{\prime }{p}^{\prime }}-\overline{v_j^{\prime }{\sigma }_{\mathit{kj}}^{\prime }})\hfill \\ \hfill & +\frac{{\tilde{\sigma }}_{\mathit{kj}}}{\overline{\rho }}{\partial }_k{\overline{v}}_j+\frac{{\partial }_j\overline{p}}{{\overline{\rho }}^2}\overline{v_j^{\prime }{\rho }^{\prime }}-\overline{{\sigma }_{\mathit{kj}}^{\prime }{\partial }_k{v}_j^{\prime }},\hfill \end{array}$$(A.15)

where Π_jj ≈ 0 was used⁹.

To close our system of equations, we need to define the remaining primed quantities. First, we model the enthalpy flux. Using the Taylor series

$$\begin{array}{c}\hfill h(\overline{s}+s^{\prime },\overline{p}+p^{\prime })\approx h(\overline{s},\overline{p})+{\frac{\partial h}{\partial s}|}_p{s}^{\prime }+{\frac{\partial h}{\partial p}|}_s{p}^{\prime }+\mathcal{O}(s{{}^{\prime }}^2,p{{}^{\prime }}^2),\end{array}$$(A.16)

and under the small Mach-number approximation (Kuhfuss 1986), we have $h\prime \approx \overline{T}s\prime$. One can use this argument to produce the source term of Eq. (A.15) as $\overline{\rho \prime v_j\prime }\approx -(\overline{\delta }\overline{\rho }/\overline{c_p})\overline{s\prime v_j\prime }$. In short, we reduce the two terms to a function of $\mathrm{\Pi }=\overline{v_j\prime s\prime }$. To produce an equation for Π, one needs an equation for s′ as well. First, we divide Eq. (A.4) by T and take the average, thus:

$$\begin{array}{cc}\hfill d_t\overline{s}+\overline{v_k^{\prime }{\partial }_k{s}^{\prime }}=\frac{1}{\overline{\rho }\overline{T}}& (-{\partial }_i\overline{F_i}+\overline{{\sigma }_{\mathit{ij}}}{\partial }_i{\overline{v}}_j+\overline{{\sigma }_{\mathit{ij}}^{\prime }{\partial }_i{v}_j^{\prime }})\hfill \\ \hfill & -\frac{1}{\overline{\rho }{\overline{T}}^2}(-\overline{T^{\prime }{\partial }_i{F}_i^{\prime }}+\overline{T^{\prime }{({\sigma }_{i}j{\partial }_i{v}_j)}^{\prime }}),\hfill \end{array}$$(A.17)

where the second term of the right-hand side provides the coupling with the entropy fluctuations. We get the equation of the fluctuating entropy by taking the difference of Eq. (A.4) and (A.17), thus

$$\begin{array}{c}\hfill d_t{s}^{\prime }+v_k^{\prime }{\partial }_k\overline{s}+{(v_k^{\prime }{\partial }_{k}s{}^{\prime })}^{\prime }=-\frac{1}{\overline{\rho }\overline{T}}[{\partial }_i{F}_i^{\prime }-{({\sigma }_{\mathit{ij}}^{\prime }{\partial }_i{v}_j^{\prime })}^{\prime }]+\text{higher order terms}..\end{array}$$(A.18)

One produces the equation of Π by multiplying Eq. (A.11) by s′ and Eq. (A.18) by v_j′, and take the average of the sum. This is

$$\begin{array}{cc}\hfill {\mathrm{d}}_t\mathrm{\Pi }& =-\frac{1}{\overline{\rho }}{\partial }_k\left(\overline{\rho }\overline{v_k^{\prime }{v}_j^{\prime }{s}^{\prime }}\right)-\frac{1}{\overline{\rho }}{\tilde{\sigma }}_{\mathit{jk}}{\partial }_k\overline{s}-\overline{v_k^{\prime }{s}^{\prime }}{\partial }_k{\overline{v}}_j\hfill \\ \hfill & -\overline{{\left(\frac{1}{\rho }{\partial }_{j}p\right)}^{\prime }{s}^{\prime }}+\frac{1}{\overline{\rho }}\overline{s^{\prime }{\partial }_k{\sigma }_{\mathit{kj}}^{\prime }}-\frac{1}{\overline{\rho }\overline{T}}[\overline{v_j^{\prime }{\partial }_i{F}_i^{\prime }}-\overline{v_j^{\prime }{({\sigma }_{\mathit{ik}}^{\prime }{\partial }_i{v}_k^{\prime })}^{\prime }}].\hfill \end{array}$$(A.19)

Here, Eq. (A.7) was used in the first term. Neglecting p′, the fourth term of the RHS reduces to $\overline{s{\prime }^2}\overline{\delta }/(\overline{\rho }\overline{c_p}){\partial }_j\overline{p}$, hence we produce our last equation to $\mathrm{\Phi }=\overline{s{\prime }^2}/2$ by multiplying Eq. (A.18) by s′ and take the average, that is

$$\begin{array}{c}\hfill {\mathrm{d}}_t\mathrm{\Phi }=-\frac{1}{2\overline{\rho }}{\partial }_k\left(\overline{\rho }\overline{v_k^{\prime }s{{}^{\prime }}^2}\right)-\overline{s^{\prime }{v}_k^{\prime }}{\partial }_k\overline{s}-\frac{1}{\overline{\rho }\overline{T}}[\overline{s^{\prime }{\partial }_i{F}_i^{\prime }}-\overline{s^{\prime }{({\sigma }_{\mathit{ij}}^{\prime }{\partial }_i{v}_j{}^{\prime })}^{\prime }}].\end{array}$$(A.20)

The remaining terms need to be closed algebraically because we want to produce a three-equation model.

Firstly, we consider the average and fluctuating radiative flux. The flux can be written as F_i = K_r(T, κ)∂_iT, then its average and first-order fluctuation is:

$$\begin{array}{cc}\hfill \overline{F_i}& =\overline{K_r}{\partial }_i\overline{T}+\overline{K_r^{\prime }{\partial }_i{T}^{\prime }},\hfill \end{array}$$(A.21) $$\begin{array}{cc}\hfill F_i^{\prime }& =K_r^{\prime }{\partial }_i\overline{T}+\overline{K_r}{\partial }_i{T}^{\prime },\hfill \end{array}$$(A.22)

where, introducing κ_T = (∂lnκ/∂lnT)_p, we have

$$\begin{array}{c}\hfill K_r^{\prime }\approx \frac{\overline{K_r}}{\overline{c_p}}(3-{\kappa }_T)s^{\prime }.\end{array}$$(A.23)

To describe the radiation dissipation effects, we have to model ${\partial }_{i}T\prime \approx \overline{T}/\overline{c_p}{\partial }_{i}s\prime$. To do this, following Kuhfuß (1987), we introduce the radiation dissipation length Λ_r, hence ∂_is′≈s′/Λ_r. This means, that the radiation fluctuation terms are modeled as:

$$\begin{array}{cc}& \frac{1}{\overline{\rho }\overline{T}}\overline{s^{\prime }{\partial }_i\overline{K_r}{\partial }_i{T}^{\prime }}\approx \frac{2\mathrm{\Phi }}{{\tau }_r},\frac{1}{\overline{\rho }\overline{T}}\overline{v_j^{\prime }{\partial }_i\overline{K_r}{\partial }_i{T}^{\prime }}\approx \frac{\mathrm{\Pi }}{{\tau }_r},\hfill \\ \hfill & \frac{1}{\overline{\rho }}{\partial }_i\overline{K_r^{\prime }{\partial }_i{T}^{\prime }}\approx \frac{2\mathrm{\Phi }\overline{T}}{\overline{c_p}{\tau }_{\kappa }},\frac{1}{\overline{\rho }\overline{T}}\overline{s^{\prime }{\partial }_i{K}_r^{\prime }{\partial }_i\overline{T}}\approx \frac{2\mathrm{\Phi }}{{\tau }_{\kappa }},\frac{1}{\overline{\rho }\overline{T}}\overline{v_j^{\prime }{\partial }_i{K}_r^{\prime }{\partial }_i\overline{T}}\approx \frac{\mathrm{\Pi }}{{\tau }_{\kappa }},\hfill \end{array}$$

where

$$\begin{array}{c}\hfill {\tau }_r=\frac{{\mathrm{\Lambda }}_r^2\overline{\rho }\overline{c_p}}{\overline{K_r}},{\tau }_{\kappa }=\frac{{\tau }_r}{3-{\kappa }_T}.\end{array}$$

Here we note that although τ_κ has the dimension of time, it is not a simple dissipation timescale, as when κ_T < −3, it provides an additional source of entropy-fluctuation.

The next term we model is the dissipation of the TKE. We use the local ${\tau }_{\mathrm{conv}}=\mathrm{\Lambda }{\stackrel{\sim }{\omega }}^{-1/2}$ approximation leading to

$$\begin{array}{c}\hfill \overline{{\sigma }_{\mathit{ij}}^{\prime }{\partial }_i{v}_j^{\prime }}\approx \frac{\tilde{\omega }}{{\tau }_{\mathrm{conv}}}={\alpha }_d\frac{{\tilde{\omega }}^{3/2}}{\mathrm{\Lambda }}.\end{array}$$

We consider the dynamic dissipation of Π and Φ to happen on similar scale then on $\stackrel{\sim }{\omega }$ thus

$$\begin{array}{cc}\hfill \frac{1}{\overline{\rho }}\overline{s^{\prime }{\partial }_k{\sigma }_{\mathit{kj}}^{\prime }}-\frac{1}{\overline{\rho }\overline{T}}\overline{v_j^{\prime }{({\sigma }_{\mathit{ik}}^{\prime }{\partial }_i{v}_k^{\prime })}^{\prime }}& \approx -{\alpha }_{d,\mathrm{\Pi }}\frac{\mathrm{\Pi }{\tilde{\omega }}^{1/2}}{\mathrm{\Lambda }}\hfill \\ \hfill -\frac{1}{\overline{\rho }\overline{T}}\overline{s^{\prime }{({\sigma }_{\mathit{ij}}^{\prime }{\partial }_i{v}_j^{\prime })}^{\prime }}& \approx -{\alpha }_{\mathrm{d},\mathrm{\Phi }}\frac{\mathrm{\Phi }{\omega }^{1/2}}{\mathrm{\Lambda }}\hfill \end{array}$$

Now, we consider the property of averaging that only keeps radial terms of the vectors. This means that $\overline{v_{j}s\prime }=\overline{v_{3}s\prime }$, and ${\partial }_j\overline{p}={\partial }_3\overline{p}{\delta }_{j3}$.

Furthermore, we define anisotropy as $\xi =\overline{v_3{\prime }^2}/(2\stackrel{\sim }{\omega })$, and assume horizontal isotropy thus, $\overline{v_1{\prime }^2}=\overline{v_2{\prime }^2}=\stackrel{\sim }{\omega }(1-\xi )$. Using this, we split the Reynolds tensor to the isotropic turbulent pressure, and anisotropic eddy viscous pressure (Gehmeyr & Winkler 1992): $-{\stackrel{\sim }{\sigma }}_{\mathit{ij}}=(2/3)\stackrel{\sim }{\omega }\overline{\rho }{\delta }_{\mathit{ij}}+P_{\mathit{ij}}^{\nu }$, so only P_ij^ν is need to be modeled, which equivalent with $\stackrel{\sim }{\omega }{b}_{\mathit{ij}}$ in Canuto & Dubovikov (1998) notation. Contributions of the eddy viscous pressure in Eq. (A.10) and Eq. (A.12) are called turbulent momentum transfer rate (U_ν) and turbulent viscous energy transfer rate (E_ν), respectively (Wuchterl & Feuchtinger 1998). Their definitions ar

$$\begin{array}{cccc}\hfill U_{\nu }& \frac{1}{\overline{\rho }}{\partial }_k{P}_{\mathit{kj}}^{\nu },\hfill & \hfill E_{\nu }& \frac{1}{\overline{\rho }}{P}_{\mathit{ij}}^{\nu }{\partial }_i{\overline{v}}_j.\hfill \end{array}$$(A.24)

Lastly, we describe the TOMs. In Eq. (A.15) we absorb pressure, viscosity transport, and the TOM into a single diffusion flux term using the diffusion approximation, we introduce the turbulent dynamic viscosity as ${\mu }_{\mathrm{t}}:=\overline{\rho }\mathrm{\Lambda }{\stackrel{\sim }{\omega }}^{1/2}$, thus

$$\begin{array}{c}\hfill {\mathcal{F}}_{\tilde{\omega }}\frac{1}{2}\overline{\rho }\overline{v_k^{\prime }{v}_j^{\prime }{v}_j^{\prime }}+\overline{v_k^{\prime }{p}^{\prime }}+\overline{v_j^{\prime }{\sigma }_{\mathit{kj}}^{\prime }}\approx -{\alpha }_t\xi {\mu }_{\mathrm{t}}{\partial }_r\tilde{\omega }.\end{array}$$(A.25)

The remaining TOMs in eqs. (A.19)-(A.20) become

$$\begin{array}{cc}\hfill {\mathcal{F}}_{\mathrm{\Pi }}& \overline{\rho }\overline{v_k^{\prime }{v}_j^{\prime }{s}^{\prime }}\approx -{\alpha }_{\mathrm{\Pi }}\xi {\mu }_{\mathrm{t}}{\partial }_r\mathrm{\Pi }\hfill \end{array}$$(A.26) $$\begin{array}{cc}\hfill {\mathcal{F}}_{\mathrm{\Phi }}& \frac{\overline{\rho }}{2}\overline{v_k^{\prime }s{{}^{\prime }}^2}\approx -{\alpha }_{\mathrm{\Phi }}\xi {\mu }_t{\partial }_r\mathrm{\Phi }\hfill \end{array}$$(A.27)

The model equations presented in Sect. 2. can be derived by using spherical symmetry and substituting all previously described terms into Eqs. (A.6), (A.10), (A.12), (A.15), (A.19), and (A.20).

Appendix B: Preliminaries for the anisotropy expression

From Canuto (1998) (25a) and (27a), one can derive the evolution equation of $\overline{v_3{\prime }^2}=-{\overline{\rho }}^{-1}{\stackrel{\sim }{\sigma }}_{33}$ using the same steps as Kupka et al. (2022). To derive the shear terms, one has to use the definition of the mean velocity, that is ${\overline{v}}_i={\partial }_{i3}\overline{v}(r)$. Therefore

$${\partial }_i{\overline{v}}_j=\left(\begin{array}{ccc}\overline{v}/r& 0& 0\\ 0& \overline{v}/r& 0\\ 0& 0& {\partial }_r{v}_r\end{array}\right)$$

and the flux term is λ_ih_j = δ_i3δ_ijS. Furthermore, we assume isotropy in the horizontal direction, thus $\overline{v_1{\prime }^2}+\overline{v_2{\prime }^2}=2\stackrel{\sim }{\omega }-\overline{v_3{\prime }^2}$.

Appendix C: Derivation of the value of α_c

The modeling of the convective flux is based on the Taylor expansion of the enthalpy as a function of temperature, which is equivalent to eq (A.16) with using: $T\prime \approx (\overline{T}/\overline{c_p})s\prime$. On the other hand, the transport of ions means that the mean molecular weight also changes, that is,

$$\begin{array}{cc}& h(\overline{T}+T^{\prime },\overline{p}+p^{\prime },\overline{\mu }+{\mu }^{\prime })=h(\overline{T},\overline{p})\hfill \\ \hfill & +{\frac{\partial h}{\partial T}|}_p{T}^{\prime }+{\frac{\partial h}{\partial p}|}_T{p}^{\prime }+{\frac{\partial h}{\partial \mu }|}_{T_{,}p}{\mu }^{\prime }+\mathcal{O}(T{{}^{\prime }}^2,p{{}^{\prime }}^2,\mu {{}^{\prime }}^2)\hfill \end{array}$$(C.1)

The absorption of ionization effects into c_p = (∂h/∂T)_p and (∂h/∂p)_T is a usual process (Kippenhahn & Weigert 1990), but in this case the specific heat is

$$\begin{array}{c}\hfill c_p\equiv {\frac{\partial h}{\partial T}|}_p=h_T+h_{\mu }{\mu }_T,\end{array}$$(C.2)

where h_T = (∂h/∂T)_p, μ, h_μ = (∂h/∂μ)_p, t, μ_T = (∂μ/∂T)_p. If we use Reynolds-averaging, we get

$$\begin{array}{cc}\hfill \overline{c_p}& =\overline{h_T}+\overline{h_{\mu }}\overline{{\mu }_T}+\overline{h_{\mu }^{\prime }{\mu }_T^{\prime }}=c_p(\overline{T},\overline{p})+\overline{h_{\mu }^{\prime }{\mu }_T^{\prime }},\hfill \end{array}$$(C.3) $$\begin{array}{cc}\hfill c_p^{\prime }& =h_T^{\prime }+\overline{h_{\mu }}{\mu }_T^{\prime }+h_{\mu }^{\prime }\overline{{\mu }_T}+{(h_{\mu }^{\prime }{\mu }_T^{\prime })}^{\prime }.\hfill \end{array}$$(C.4)

Now, we assume that all fluctuations are due to temperature fluctuation in the first order. Thus c_p′=c_pTT′, h_T′=h_TTT′ $h_{\mu }\prime =h_{\mu T}T\prime$, where T in the subscript means derivative by T, for example, h_TT = (∂²h/∂T²)_p, and all second derivatives are considered to be functions of average values only¹⁰. Substituting these into (C.4) and dividing by T′ we get

$$\begin{array}{c}\hfill c_{\mathit{pT}}=h_{\mathit{TT}}+\overline{h_{\mu }}{\mu }_{\mathit{TT}}+h_{\mu T}\overline{{\mu }_T}+h_{\mu T}{\mu }_{\mathit{TT}}{T}^{\prime }-h_{\mu T}{\mu }_{\mathit{TT}}\overline{T{{}^{\prime }}^2}/T^{\prime }\end{array}$$(C.5)

Let $c_{p0}=c_p(\overline{T},\overline{p})=\overline{h_T}+\overline{h_{\mu }}\overline{{\mu }_T}$ and $c_{pT0}=\overline{h_{\mathit{TT}}}+\overline{h_{\mu T}}\overline{{\mu }_T}+\overline{h_{\mu }}\overline{{\mu }_{\mathit{TT}}}$. These are the specific values that can be exactly determined from the mean values, and as we can see, c_p ≠ c_p0. Now substitute these back into the definition of h′, and we get

$$\begin{array}{c}\hfill h^{\prime }\approx c_{p0}{T}^{\prime }+h_{\mu T}{\mu }_{\mathit{TT}}\overline{T{{}^{\prime }}^2}{T}^{\prime }+c_{pT0}T{{}^{\prime }}^2+h_{\mu T}{\mu }_{\mathit{TT}}T{{}^{\prime }}^3-h_{\mu T}{\mu }_{\mathit{TT}}\overline{T{{}^{\prime }}^2}{T}^{\prime },\end{array}$$(C.6)

therefore setting α_c = 1 is adequate to follow ionization changes in first-order.