© ESO 2025

1. Introduction

Since the discovery of the fundamental relationship between their pulsation period and their luminosity (the Leavitt Law, Leavitt & Pickering 1912), classical Cepheids have had a large influence on astronomy and astrophysics. Bright Cepheids are easily detected both inside and outside the Galaxy because the period of their radial pulsations typically ranges from a few days to a few weeks (e.g., Bono et al. 2005; Anderson et al. 2014). Measurements and models for these variable stars provide a foundation for studies with subjects as diverse as the spiral structure of the Milky Way (e.g., Majaess et al. 2009; Dambis et al. 2015), extra-galactic distances (see reviews Feast & Walker 1987; Feast 1999), and estimation of the Hubble constant (Ngeow & Kanbur 2006; Van Leeuwen et al. 2007; Tammann et al. 2008; Freedman & Madore 2010; Freedman 2021). Because of the period–luminosity (PL) relationship, the study of Cepheids has been an active research area for several decades (Madore & Freedman 1991; Turner 2010; Bono et al. 2010), including work on the dependence of their properties on metallicity (Alibert et al. 1999; Persson et al. 2004).

Despite their importance to astronomers, a complete theoretical picture of the interior of Cepheids is lacking (see the recent review by Guzik et al. 2023). One clear indication of the limitations of current models is the mass discrepancy problem. A typical disagreement of 10% to 20% is found when comparing the masses of Cepheids obtained with stellar evolution models and the ones calculated based on stellar pulsation measurements (Bono et al. 2006; Keller 2008; Pietrzyński et al. 2010; Cassisi & Salaris 2011). Dynamical masses of Cepheids in binary systems (Evans et al. 1997, 2024; Gallenne et al. 2018) are lower than the masses provided by evolution models, and more consistent with pulsation calculations (e.g., Stobie 1969; Neilson et al. 2011), highlighting that improvement of the current stellar evolution models of Cepheids is a high priority.

In the stellar interior, convection underpins the fundamental processes of heat transport, fluid mixing, and shear. One-dimensional (1D) parameterizations of convection are thus central to stellar evolution models. The most widely used method to model convection is stellar mixing length theory (MLT; Vitense 1953; Böhm-Vitense 1958; Salaris & Cassisi 2008, 2017). Because MLT relies on a single parameter, the mixing length, it cannot accurately capture the multi-scale, non-local, multi-dimensional, and intermittent nature of turbulent convection. MLT is used to model evolutionary tracks for Cepheids before they enter the instability strip. However, more sophisticated convection models are generally used to evolve stars while they experience radial pulsation. This is necessary because radial pulsations can happen on a time-scale comparable to the convection (e.g., Deupree 1984, 1986; Glasner & Livne 1995), implying that these two effects can have a two-way interaction. Observations have recently verified that such an interaction can occur. A recent study of R Car (Rosales-Guzmán et al. 2024), an M-type Mira variable star with a pulsation period of 314 days that is evolving along the asymptotic giant branch (AGB), has reported that convective structures grow when the star expands and decrease when the star contracts. Because classical Cepheids have a relatively short pulsation period compared to other evolved variable stars, they are likely to be affected more strongly by this interaction. Time-dependent convection (TDC) models have been developed to model convection during radial pulsations (e.g., Marconi 2017). The radial stellar pulsation (RSP) functionality of the Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2019) models nonlinear pulsations of various pulsating stars by including the TDC model described by Kuhfuß (1986), Wuchterl & Feuchtinger (1998), and Smolec & Moskalik (2008). The hydrodynamic simulations presented in this work have the potential to contribute to improvements in both the MLT and TDC models as they are applied to the evolution of Cepheids.

In addition to convection, convective overshooting could play a role in the mass discrepancy problem (Chiosi et al. 1992; Bono et al. 1999). The effects of convective boundary mixing (e.g., Paxton et al. 2010; Freytag et al. 1996; Pratt et al. 2017) have been modeled as an add-on to stellar MLT in two different ways. One is with a simple overshooting length, expressed as a fraction of the pressure scale height (Schwarzschild 1975) evaluated at the convective boundary; this overshooting model based on a step function defines a region in the stable regions (RZ, radiative zone) where the evolution model assumes full convective mixing. A more precise approach is based on the formulation of a one-dimensional diffusion coefficient that describes the enhanced mixing in the overshooting layer due to convective motions as a smooth function of the star’s internal radius (see, e.g., Freytag et al. 2010; Noels et al. 2010; Zhang 2013; Pratt et al. 2017). As with the step-function model, models based on diffusion coefficients require the input of a length scale, and possibly additional parameters, that must be calibrated. Multidimensional hydrodynamic simulations of convection can provide insights that improve one-dimensional stellar evolution models, and calibrate the relevant parameters (e.g., Freytag et al. 1996; Herwig 2000; Rogers et al. 2006; Arnett et al. 2019; Käpylä et al. 2017; Pratt et al. 2017; Baraffe et al. 2017; Pratt et al. 2020; Käpylä 2024).

Work involving multidimensional hydrodynamic simulations of classical Cepheids is limited (see, e.g., Mundprecht et al. 2013, 2015); in addition, these works study Cepheid models that include only a thin outer envelope that is convectively unstable. Many Cepheid models involve an internal convective shell that occurs in the Z-bump region, as well as a convective envelope. In the setting of a thin convective shell, plumes can overshoot the convection zone both above and below the shell. Convective shells, which are common in stars on the AGB, have received less attention than convective envelopes. Stars with convective shells have been investigated in the context of nuclear burning (e.g., Meakin & Arnett 2006, 2007; Mocak et al. 2011; Rizzuti et al. 2022, 2023, 2024; Georgy et al. 2024). Cristini et al. (2017, 2019) and Rizzuti et al. (2022) performed 3D hydrodynamic “box-in-star” simulations of a massive star and found that convective overshooting is weaker at the bottom of the convective shell than at the top with a factor of approximately one fifth. However, the depth and thickness of the convective shells in previous hydrodynamic studies have not been clearly documented. Weaker overshooting beneath a convective shell makes intuitive sense because the fluid there is denser than the material above the shell, which could result in shallower penetration of the convective plumes. The convective shells in AGB stars are often thicker and located deeper within the star than in the Cepheids we examine. Using a custom model in MESA, Li (2012, 2017) and Guo & Li (2021) have found that for thin convective shells in A-type stars, the extent of the overshooting layer above the shell is approximately twice the extent of the one below. However the difference that these works find between overshooting depths above and below the convective shell raises several important questions about (1) whether the overshooting lengths above and below a convective shell can be related to one another, (2) what other stellar properties affect this relationship, and (3) how models of convection and convective overshooting can be improved for Cepheids.

A high radial resolution is necessary to accurately resolve overshooting; a long time sequence of data is also important for studying overshooting so that well-resolved statistics of this intermittent process can be built. For these reasons, we chose to study two-dimensional (2D) simulations in this work. Three-dimensional (3D) simulations are desirable, and will be explored in future work. 2D simulations of stellar convection have resulted in higher velocities than three-dimensional (3D) simulations (Muthsam et al. 1995; Meakin & Arnett 2007; Arnett & Meakin 2011); however, the difference between these velocities is dependent on the stellar structure. Despite this general difference, at the convective boundaries radial velocities are similar in 2D and 3D simulations. At the radial resolutions, and for the more limited time sequence of data, that can be produced in 3D, the measured overshooting lengths are also similar to the 2D results (Pratt et al. 2020; Dethero et al. 2024).

This work is structured as follows. In Sect. 2, we discuss the realistic global fluid simulations that we have performed with the MUlti-dimensional Stellar Implicit Code (MUSIC) for a selection of these stellar structures. In Sect. 3, we examine the shape of the convection in the inner convective shell and its overshooting layers. We propose a new statistical analysis that allows for the direct comparison of overshooting above and below a convective shell. In Sect. 4, we discuss the implications of our results.

2. Hydrodynamic simulations

An interior convective shell is a common feature for Cepheids (see Appendix A for the detailed analysis of a grid of 48 MESA models). This structure is illustrated in Fig. 1. Evolutionary models of stars that may evolve into Cepheids are sensitive to the choice of parameters for convection and overshooting; therefore the presence of an additional convective layer is an important consideration for the improvement of stellar structure and evolution modeling of these stars. Hydrodynamic simulations of Cepheids, or variable stars more broadly (Geroux & Deupree 2011, 2013, 2014, 2015; Mundprecht et al. 2013, 2015) have not yet examined models that include an interior convective shell; the present work is a first effort to characterize it.

Fig. 1. Schematic of the typical structure of the Cepheids we study in this work with MUSIC. Convective regions are colored, including: a tiny convective core, an inner convective shell surrounded by radiative zones, and a thin outer convective envelope. The simulation domain is a spherical shell, and is outlined with black; it encapsulates the interior convection zone along with the surrounding radiative zones, and is truncated just below the outer convective envelope.

We have produced six models of stars of intermediate mass to examine in hydrodynamic simulations; these stars range from 5.75 to 9 M_⊙, and all lie in a blue loop and inside the instability strip. The MESA inlists for these six stars are based on inlists published by Gilkis et al. (2019) and Espinoza-Arancibia et al. (2022). The parameters of the MESA models used in our hydrodynamic simulations are summarized in Table 1. They have a metallicity of z = 0.020, which ensures the existence of the deep inner convective shell for all the stellar masses that we consider, and is representative of Cepheids found in the Milky Way. For the 7 M_⊙ star, we studied two different models, denoted (a) and (b), to allow for a comparison between otherwise similarly structured stars on the second and third crossings of the instability strip. The position of these stars on their evolutionary tracks and the placement of their internal convection zones are provided in Fig. 2. Above the internal convective shell, the structure of the Cepheids that we study involves the presence of a thin outer convective envelope. The bottom of this outer convective envelope approximately coincides with the helium-II ionization point, around a temperature of 5 × 10⁴ K (Seaton 1993b). In this outer region, the temperature and density gradients vary more rapidly, requiring significantly higher resolution; including this zone in our simulations with any accuracy would require us to increase the grid sizes by an order of magnitude. We instead truncate the present simulation volumes at a point just below the outer convective envelope, which allows us to study overshooting above the convective shell without the influence of overshooting below the outer convective envelope.

Fig. 2. (a) Hertzprung-Russel diagram and (b) radial profile of the Schwarzschild discriminant, for the six stars we simulated with MUSIC. Symbols indicate the Cepheids we simulated, and the pink shaded area indicates the instability strip for the evolutionary tracks. Shaded regions also highlight the internal convection zones around which our hydrodynamic simulations are centered.

We have performed 2D Implicit Large Eddy Simulations (ILES; Grinstein et al. 2007; Ritos et al. 2018) of these stars using the MUSIC code. Our simulations in this work only take convection into account; the possibility of studying additional physical effects such as rotation, a tachocline, chemical mixing, and magnetic fields is omitted from the current study. Of particular relevance is the assumption that each star that we simulate has a homogeneous chemical composition. This assumption is consistent with the MESA stellar structures used for these simulations. Because the helium-II ionization point is not included in the simulation volume, radial pulsations should not occur in the present simulations, and we confirm that they do not. The complex interaction between a star’s radial pulsations and convection will be studied in future work.

The MUSIC code solves the inviscid compressible hydrodynamic equations for density ρ, momentum ρu, and internal energy ρe:

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\rho & =-\mathrm{\nabla }·(\rho \mathit{u}),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\rho \mathit{u}& =-\mathrm{\nabla }·(\rho \mathit{u}\mathit{u})-\mathrm{\nabla }p+\rho \mathit{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\rho e& =-\mathrm{\nabla }·(\rho e\mathit{u})-p\mathrm{\nabla }·\mathit{u}+\mathrm{\nabla }·(\chi \mathrm{\nabla }T).\hfill \end{array}$$(3)

The code uses a finite volume method, a MUSCL method (Thornber et al. 2008) of interpolation, and a van Leer flux limiter (as described by Van Leer 1974; Roe 1986; LeVeque et al. 2006). Time integration in the MUSIC code is fully implicit and uses a Jacobian free Newton-Krylov (JFNK) solver (Knoll & Keyes 2004) with physics-based preconditioning (Viallet et al. 2016; Goffrey et al. 2017). The MUSIC code uses an adaptive time step, which is identically constrained for all simulations in this work. The MESA code uses the OPAL Type II opacities and corresponding equation of state to produce these models; our MUSIC simulations use tables generated from MESA. In Eq. (2), g is the gravitational acceleration, a spherically symmetric vector, consistent with that used in the stellar evolution calculation, and not evolved during the present simulations.

We studied the MUSIC simulations described in Table 2. For all of these simulations, the compressible hydrodynamic equations (1)–(3) are solved in a spherical shell using spherical coordinates (r, θ). Each simulation is two dimensional, and assumes azimuthal symmetry. In the table, the inner and outer radii of the spherical shell are noted, as well as the upper and lower boundaries of the convection zone. The simulation domain includes the interior convection zone and is truncated just below the outer convective envelope (see Fig. 1).

In Table 2, the pressure scale height is compared to the grid spacing at a convective boundary to measure how well the physics around that boundary is resolved. This results in two resolution criteria for a convective shell: the upper resolution $H_{\mathsf{p},\mathsf{CB}}^{\mathsf{U}}/\mathrm{\Delta }r$ and the lower resolution $H_{\mathsf{p},\mathsf{CB}}^{\mathsf{L}}/\mathrm{\Delta }r$. These two resolutions are similar. In each case, there are more than 100 grid spaces per pressure scale height, ensuring a high resolution of the overshooting layer. Each of our simulations has sufficient radial resolution to produce a characteristic radial profile for velocity in 2D, and convergence toward a velocity profile is observed as the grid is increased. Simulation ceph8 has a grid of size n_r × n_θ = 1728 × 1024.

Our simulations have boundary conditions that maintain the original radial profiles of density and temperature of the 1D stellar evolution model. We held the energy flux and luminosity constant on the outer radial boundary at values established by the stellar structure. For an examination of how boundary conditions affect the dynamics in stars, we refer to Pratt et al. (2016), Vlaykov et al. (2022). Aside from this surface boundary condition, in velocity we imposed non-penetrative and stress-free boundary conditions in the radial directions. The energy flux and luminosity were held constant at the inner radial boundary, at the value of the energy flux at that radius in the one-dimensional stellar evolution calculation. On both the inner and outer radial boundaries of the spherical shell, we imposed a boundary condition on the density that maintained hydrostatic equilibrium by assuming constant internal energy and constant radial acceleration due to gravity in the boundary cells (Grimm-Strele et al. 2015). We imposed periodicity on all physical quantities at the boundaries in θ.

We define the convective turnover time following Pratt et al. (2016):

$$\begin{array}{c}\hfill {\tau }_{\mathsf{conv}}(t)={\displaystyle {\int }_{\mathsf{CZ}}}\mathsf{dV}\frac{H_{\mathsf{p}}}{\left|u\right|}{\displaystyle {\int }_{\mathsf{CZ}}}\mathsf{dV},\end{array}$$(4)

where H_p is the pressure scale height and ${\displaystyle \left|u\right|=\sqrt{u_r^2+u_{\theta }^2}}$ the magnitude of the velocity computed in 2D from the radial and angular velocities, u_r and u_θ respectively. The integration covers the convection zone and is volume-weighted using the volume element V. We use the convective turnover time to verify that convection has reached a steady-state, a period where the time-averaged value of the convective turnover time is well-defined and not evolving in time. For each simulation, the time-averaged value of τ_conv in the steady-state regime is included in Table 2 alongside its standard deviation. In general, when we describe time-averaged quantities, this refers to a time average taken over the entire time in steady state covered by the simulation. This time span is more than 25τ_conv in every case, ensuring enough data to produce well-converged statistics. The convective boundaries have not evolved during the time span of our hydrodynamic simulations: the temperature at the Schwarzschild boundaries varies by less than 0.8% for all models. This is expected because the Kelvin-Helmholtz timescale of the stars we have examined is at least 4 orders of magnitude longer than the time covered by our simulations.

In our simulations, the maximum amplitude of the convective velocity in the convective shell is approximately proportional to the cube-root of luminosity, ${\displaystyle u_{\mathsf{rms}}\propto L_{\mathsf{star}}^{1/3}}$ (see Fig. 4), a result also found by previous works investigating convective shells (e.g., Jones et al. 2016; Andrássy et al. 2020), convective cores (e.g., Higl et al. 2021; Horst et al. 2020; Baraffe et al. 2023) and convective envelopes (e.g., Baraffe et al. 2021). The root-mean-squared (RMS) velocities just above and below the convective boundaries (indicated by vertical lines in Fig. 4) are caused by convective overshooting that reaches into the neighboring stable regions. Far from the convection zone, the non-zero velocities are due to the propagation of internal gravity waves (IGWs), excited by the convective motions and overshooting plumes. Such waves, as well as the internal convective shell, can be observed in visualizations of vorticity and radial and angular velocities (see Fig. 3). The present paper has focused on convection and convective boundary mixing, and the propagation of IGWs has not been investigated further.

Fig. 3. Visualizations of (from left to right) vorticity magnitude, radial velocity, and angular velocity of the 8 M⊙ simulation after 30τconv of steady-state convection. The zero point for vorticity is in dark red in the left panel. In the center and right panels, outwards flows are in red while inwards flows are in yellow; the zero point in velocity is black. The radial velocity ranges from −9.25 × 105 cm/s to 9.25 × 105 cm/s at this instant in time, while the angular velocity ranges from −1.01 × 106 cm/s to 1.01 × 106 cm/s. Dashed white lines indicate the position of the Schwarzschild boundaries.

Fig. 4. Radial profile of the time-averaged RMS velocity scaled by (Lstar/1035)1/3 for a range of stellar masses. Vertical lines indicate the Schwarzschild convective boundaries for each star.

3. Results

To describe convection in stellar interiors, the concept of a filling factor has been applied widely (e.g., Zahn 1991; Canuto & Dubovikov 1998; Dethero et al. 2024). A filling factor is a crude measure that reduces the complex flow of stellar convection to a one-dimensional measure of the amount of fluid moving inward (or outward) at a given radius. A large part of the value of a filling factor is as a measure of the asymmetry between inflows and outflows (e.g., Zhang et al. 1997; Brown & Ahlers 2007; Dethero et al. 2024).

Several definitions of the filling factor have been proposed in the literature (Dethero et al. 2024). We have considered the filling factor based on a simple volume-percentage of inflows at a given radius (Schmitt et al. 1984; Brummell et al. 2002; Andrássy 2015; Käpylä 2024). The volume-percentage filling factor is defined as the fraction of volume occupied by the inflows:

$$\begin{array}{c}\hfill {\sigma }_{\mathsf{vp},\mathsf{in}}=\frac{V^{\mathsf{in}}(r,t)}{V^{\mathsf{in}}(r,t)+V^{\mathsf{out}}(r,t)}.\end{array}$$(5)

Here Vⁱⁿ(r, t) indicates the total volume of grid cells at a given radius, angular position, and instant in time that has an inward velocity. Similarly, V^out(r, t) indicates the total volume of grid cells with an outward velocity. A value of σ_vp, in = 1/2 indicates symmetry between inflows and outflows, and is expected in the stable radiative regions. In convection zones and overshooting layers, it is possible to obtain highly asymmetric convection patterns, depending on the density stratification of the star in those layers. The radial profiles of the time-averaged volume-percentage filling factor for inflows remain close to 0.5 (see Fig. 5); this indicates only a slight asymmetry between inflows and outflows in the convective shells considered in this work. A filling factor less than one half corresponds to inflows that are thinner, faster, and more concentrated than outflows. For the convective shells that we have studied, σ_vp, in has a distinctive sinusoidal shape within the convectively unstable region; it is highest at the outer boundary of the convective shell, and drops to its lowest value at the lower convective boundary. This trend is clearly identifiable for all six of the stars, and it is the opposite of what is generally found for outer convective envelopes in other star types (Dethero et al. 2024). Observations indicate filling factors less than one half are found at the solar surface; MUSIC simulations have previously reproduced this observational result, and showed that the filling factor rises to one half by the convective boundary (Dethero et al. 2024). Thus, although the convective shells that we consider are in the outer regions of the star, where convective envelopes are commonly studied, the radial trend of asymmetry is more similar to convective cores. Although statistically significant, this trend is slight, and the convective shells are therefore most similar to laboratory experiments of Boussinesq convection.

Fig. 5. Radial profile of the time-averaged volume percentage filling factor of the inward moving plumes σvp, in versus the star’s internal radius, in units of the total stellar radius R for the 6 M⊙ simulation. The shaded area indicates one standard deviation above and below this averaged line and thin vertical lines indicate the Schwarzschild boundaries.

A motivation for the concept of the filling factor is to establish a universal diagnostic to quantify the symmetry or the asymmetry of inflows and outflows. However, as pointed out recently by Dethero et al. (2024), the filling factor is a number that sums up the width of inflows, and is thus related to a low-order statistic. As an attempt to pursue high-order statistics, Dethero et al. (2024) proposed to examine the width of inflows directly. The width of inflows (outflows) is defined as the width of a continuous set of cells in the angular direction, at a given radius, which all have a negative (positive) radial velocity. The multi-scale nature of stellar convection can be understood through the radial distribution of plume widths. For large convection zones, an extreme density stratification can lead to a wide range of convective scales.

In our simulations of relatively thin convective shells, the width of inflows in the convection zones varies only mildly except near the boundaries (see Fig. 6). For example, the time-averaged width of inflows increases from only 1% of the volume at that radius between the two boundary layers of the 5.75 M_⊙ Cepheid, and 3% for the 9 M_⊙ simulation, ultimately leading to convective rolls of almost the same size in the heart of the convection zone. The overshooting layers both above and below the convection zones are readily identified by minima in the widths of convective structures. These minima appear to indicate a radial position where plumes defined by radial velocity dissipate to their thinnest extent; the angular resolution of the present simulations would be able to resolve plumes ten times smaller if they occurred. Our diagnostic continues to calculate widths of radial flows measured further from the convection zone; these are contaminated by the small radial velocities of IGWs, which tend to have larger angular extents. Any detailed study of such wave dynamics will be undertaken in future work. The width of inflows in the top overshooting layer is systematically larger than those from the bottom overshooting layers, as indicated by the relative minimum values. Different scales of the convection rolls exist in the two distinct overshooting layers. This is in agreement with the results obtained with the volume-percentage filling factor, i.e inflows are fatter at the top while outflows are fatter at the bottom of the convective shell.

Fig. 6. Width of inflowing plumes Win for the 5.75 and the 9 M⊙ simulations as a function of the star’s radius. Shaded areas indicate one standard deviation above and below the time-averaged line. Thin vertical lines indicate the radial position of the Schwarzschild boundaries.

Both the volume-percentage filling factor (Fig. 5) and the width of the inflowing plumes (Fig. 6) demonstrate that inside these thin convective shells, apart from the boundary layers, convective flows of approximately the same size dominate. We did not observe a multi-scale convective flow to the extent described in larger convective envelopes. This result can be understood most clearly by looking at the visualization of vorticity in Fig. 3. It shows the presence of convection rolls slightly larger than the width of the convection zone, which overshoot both convective boundaries. Large convective rolls dominate in the convectively unstable region, while thinner and faster inflowing plumes are found close to the Schwarzschild boundaries, consistent with a picture of boundary layer flows.

3.1. Determination of the extent of the overshooting layer

Convective overshooting is a fundamental phenomenon of stellar physics. Its influence on one-dimensional stellar evolution models is typically accounted for through an overshooting length. Determining the extent of convective overshooting requires a method for determining what part of the fluid is inside a convective flow structure and what part is outside. This is potentially a complex problem because convection drives shear flows and turbulence, and because convective flow structures can be defined both by their large-scale radial motions and by their ability to transport heat. Several definitions have been used to measure the depth of the overshooting layer: the thermal flux, the kinetic energy flux, and the vertical kinetic flux. In this work, we have used the vertical kinetic energy flux to define the extent of convective overshooting (e.g., Hurlburt et al. 1986, 1994; Ziegler & Rüdiger 2003; Rogers et al. 2006; Tian et al. 2009; Chan et al. 2010; Pratt et al. 2017). The vertical kinetic energy flux F_k is defined:

$$\begin{array}{c}\hfill F_{\mathsf{k}}=\frac{1}{2}{u}_r\rho {\left|u\right|}^2.\end{array}$$(6)

For convective shells where overshooting can happen both above and below the convective layer, we have examined the depth of the upper and lower overshooting layers using the first zero of the vertical kinetic energy flux in the stable radiative zone. This is a method that has been used by Pratt et al. (2017, 2020), Baraffe et al. (2021), Vlaykov et al. (2022), Baraffe et al. (2023), Dethero et al. (2024).

In our hydrodynamic simulations with no rotation, the overshooting length is statistically independent of the angular position θ. However, convective overshooting varies in θ at any given instant in time (see Fig. 7a). The vertical kinetic flux in the overshooting layer exhibits round shapes over several degrees. These shapes represent large convective rolls protruding from each side of the convection zone that contribute to large overshooting layers (see Fig. 7b). Instead of dissipating right after the convective boundaries, convective motions that are parts of these convection rolls return into the CZ. A straightforward approach to evaluate a general overshooting length is to calculate a simple time and space average. For example, the dashed line in Fig. 7a indicates the average in θ of the overshooting length at the particular instant in time displayed. Despite its simplicity, this straightforward approach hides physically important features of convective overshooting and masks contributions of the convective motions that penetrate further than the average overshooting plume.

Fig. 7. (a) Angular structure of the overshooting layer after 114τconv of steady-state convection in the simulation of the 6 M⊙ Cepheid. The overshooting length in this illustration is determined by zeros of the vertical kinetic flux. The Schwarzschild convective boundaries between the convection zone and the neighboring stable radiative zones are indicated by a solid black line. The stellar radius (r/R) is indicated on the left axis, and the overshooting length is given in units of the pressure scale height Hp, CB at the convective boundaries on the right axis. A dashed line indicates the average overshooting length at this time. (b) The same, but plotted over a visualization of vorticity. White dashed lines indicate the Schwarzschild boundaries.

A method used historically (Browning et al. 2004; Brun et al. 2011) to calculate the radial extent of the overshooting layer uses the time-averaged radial profile of the convective flux. We have defined the convective flux following Freytag et al. (1996) and Pratt et al. (2016):

$$\begin{array}{c}\hfill F_{\mathsf{conv}}=\langle H\rho u_r\rangle -\langle H\rangle \langle \rho u_r\rangle ,\end{array}$$(7)

where ⟨▫⟩ denotes a time average and the enthalpy is H = e + p/ρ. From this flux, the overshooting layer is defined using the negative peak in the convective flux that occurs around the convective boundary; the radial point in the radiative zone where this peak becomes negligible defines the extent of the overshooting layer. This method of calculating the overshooting depth has two disadvantages. Firstly, the volume-averaged profile of the convective flux requires a long-time average to converge to a smooth profile (see Fig. 8). The convective flux thus provides a useful approximation of the width of the overshooting layer, but cannot be regarded as a precise measurement. Secondly, in the overshooting layer and radiative zone, the time-averaged profile of the convective flux is contaminated by small-scale oscillations at the point where it would become negligible; these may be the signature of internal gravity waves. As a consequence, some judgment must be exercized to estimate this point. We have used the convective flux estimate of the overshooting layer as a reasonable check for the more accurate statistical methods that we describe next. Overall, the convective flux shows that the size of the overshooting layers is of the same order of magnitude as the size of the convection zone. Approximate values of the overshooting depth for the upper and lower overshooting layers, calculated on the basis of the convective flux profile are included in Table 3, and are denoted $\overline{{\displaystyle {\ell }_{\mathsf{ov}}}}$. A superscript U or L differentiates the upper and lower overshooting layers.

Fig. 8. Radial profiles of the convective flux for the 6 M⊙ Cepheid. The shaded area indicates one standard deviation around the mean value. Thin vertical lines indicate the convective boundaries determined by the Schwarzschild criterion. Black dashed lines indicates one pressure scale height above and below the convection zone. The inset shows a zoomed in region centered on the red box.

A more precise method for defining the extent of the overshooting layer is to use instantaneous values of the vertical kinetic energy flux. We have defined the position r_o as the position of the first zero of the vertical kinetic energy flux. The probability density functions (PDF) of the overshooting extent of all plumes for the 6 and the 7 M_⊙(b) are shown in Fig. 9; for the other simulations, the PDF(r_o) have a similar shape. As observed by Pratt et al. (2017) in the case of large convective envelopes, these PDFs have a single tail that is heavier than a Gaussian distribution, and are asymmetrical because plumes do not stop inside the convection zone. These heavy tails suggest that strong plumes reach deeply in the radiative zone. A straightforward average of r_o − R_CB is therefore not sufficient to determine an overshooting length to describe how convective overshooting affects the stable radiative layer. If a distribution is Gaussian, it can be fully characterized by a mean and a standard deviation. However, for a non-Gaussian distribution such as the ones produced here, using the average eliminates the most significant part of this intermittent fluid process. In this section, we will first seek to apply the strategy of using extreme value theory (EVT) to evaluate an overshooting length both above and below the convective shell.

Fig. 9. Probability density functions of ro, the extent of all overshooting plumes, for the 6 and the 7 M⊙(b) Cepheid. Thin vertical lines indicate the convective boundaries determined by the Schwarzschild criterion.

3.2. Application of the EVT

To describe the influence of the most energetic overshooting plumes which penetrate further in the RZ and are responsible for the heavy tail in the PDF(r_o), we have first followed the approach developed by Pratt et al. (2017) using the generalized extreme value distribution (GEVD; Castillo et al. 2005). We have defined the maximal overshooting depth, r_max, to be the lowest (highest) position in the radiative zone below (above) the convective shell that is reached at a given time. Below the convective shell, this reads

$$\begin{array}{c}\hfill r_{\mathsf{max}}(t)={\min }_{\theta }(r_{\mathsf{o}}(\theta ,t)),\end{array}$$(8)

and above

$$\begin{array}{c}\hfill r_{\mathsf{max}}(t)={\max }_{\theta }(r_{\mathsf{o}}(\theta ,t)).\end{array}$$(9)

The maximal overshooting length Δr_max is then the difference between r_max and the convective boundary R_CB:

$$\begin{array}{c}\hfill \mathrm{\Delta }{r}_{\mathsf{max}}(t)=|r_{\mathsf{max}}(t)-R_{\mathsf{CB}}|.\end{array}$$(10)

In the following, Δr_max will be made dimensionless with the pressure scale height at the convective boundary, H_p, CB.

3.2.1. Below the convection zone

Applying EVT, the cumulative distribution function (CDF) of the maximal overshooting length should follow a Weibull distribution:

$$\begin{array}{c}\hfill F(\mathrm{\Delta }{r}_{\mathsf{max}})=exp(-{(1+\alpha \left(\frac{\mathrm{\Delta }{r}_{\mathsf{max}}-\mu }{\lambda }\right))}^{-1/\alpha }),\end{array}$$(11)

where α < 0 is the shape parameter, λ is the scale parameter and μ is the location parameter. This distribution provides a good fit to the maximum overshooting data generated by simulations of convective envelopes. Pratt et al. (2017) demonstrated that the association of the overshooting length with the location parameter of the fitted PDF accounts for the effects of the strongest plumes in the stable radiative zones. For overshooting below a convective shell, the CDFs of the maximal penetration length below the convective shell are also fit reasonably well by a Weibull distribution (see Fig. 10). The downward curvature of the CDFs shows that they also follow a Weibull distribution for the current simulations. For each simulation, our EVT analysis provides a maximal width, denoted ℓ_ov^max, and derived from the location parameter. These values are summarized in Table 3. The overshooting lengths we have extracted from the GEVD fit are significantly larger than what is obtained from a straightforward averaged as done in Fig. 7a, in agreement with the conclusions of Pratt et al. (2017).

Fig. 10. Natural logarithm of the negative natural logarithm of the cumulative distribution functions of maximal overshooting length, Δrmax, for each simulation below the convective shell. Black lines indicate the best fit with the Weibull distribution, defined in Eq. (11).

3.2.2. Above the convection zone: Convective boundary mixing throughout the radiative layer

As explained in Sect. 2, the simulation domain has been truncated at a point just below the outer convective envelope. In our simulations, this presents a difficulty because plumes overshooting above the convective shell sometimes reach the outer boundary of the simulation; this is rarely the case below the convective shell. Because sometimes plumes reach the outer simulation boundary, maximal-in-time plumes often reach the simulation boundary. Therefore, the EVT approach produces an overshooting length equal to the depth of the radiative layer for all of the stars we studied except for ceph5. The simulation ceph5 represents the thinnest convective shell; the mass and luminosity of this star are lower than for the other stars we studied, resulting in lower velocities at the convective boundary. Therefore, the overshooting plumes stop closer to the CZ, and the RZ is larger than the full overshooting layer.

However, the majority of the stars we examined have higher mass and luminosity; the convective motions reach further in the RZ. For M > 6 M_⊙, the EVT analysis shows that the extreme events often hit the outer boundary of the spherical shell. The shape of the PDF(Δr_max) approaches a δ distribution; a δ distribution would be achieved if, at any time, there is at least one plume that hits the upper boundary of the domain. A radiative zone that is fully dominated by overshooting plumes has been described in recent work by Guo & Li (2021), who studied A-type stars. The EVT analysis confirms overshooting throughout the RZ. However, it is unable to provide us with a detailed picture of the upper overshooting layer that can be compared with the lower overshooting layer.

3.3. Description of the PDF(r_o) using a mixture distribution of gamma functions

The EVT analysis of overshooting was developed by Pratt et al. (2017) to describe the intermittency of overshooting and the influence of the most energetic overshooting events. Indeed, the intermittency of this process is underestimated when a simple average is used because of the non-Gaussian nature of the PDF(r_o). This approach provides us with a maximal width for the overshooting layers. However, we need a more detailed picture of the upper overshooting layer to compare with overshooting below the convective shell. Here, we propose a new way to characterize the intermittency of overshooting events. Strongly overshooting plumes are responsible for the heavy tail in the PDF(r_o). To model the non-Gaussian nature of the distribution of all plumes, we propose to separate the contribution of the strongest plumes from that of the plumes that simply make it just beyond the boundary of the convective instability. These weakest overshooting plumes play a negligible role in the enhanced diffusion throughout the radiative zone. The PDF(r_o) reveals that most overshooting plumes stop immediately after the convective boundary (see Fig. 9). Taking the average value of these non-Gaussian distributions therefore overestimates the influence of the weakest plumes and underestimates how the strongest overshooting events enhance the mixing mechanisms in the RZ.

3.3.1. Description of the mixture model

In an effort to separate the contributions of the weak overshooting from the strong overshooting, we describe these distributions with a mixture model of two gamma distributions. First, we consider the strictly positive random variable Δr_o, defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{r}_{\mathsf{o}}=\left|\frac{r_{\mathsf{o}}-R_{\mathsf{CB}}}{H_{\mathsf{p},\mathsf{C}\mathsf{B}}}\right|.\end{array}$$(12)

The choice of the gamma distribution was motivated by the fact that it is the maximum entropy probability distribution for a strictly positive random variable, given its mean value and the mean value of its logarithm (Park & Bera 2009). Essentially, the gamma distribution is the most general distribution that describes a strictly positive random variable with a positive skewness.

The probability density function of a gamma distribution has the form

$$\begin{array}{c}\hfill f(x;\alpha ,\lambda ,\mu )=\frac{1}{\mathrm{\Gamma }(\alpha ){\lambda }^{\alpha }}{(x-\mu )}^{\alpha -1}exp\left(\frac{-(x-\mu )}{\lambda }\right).\end{array}$$(13)

Here x is a coordinate that will be substituted for the variable Δr_o in our calculations. The cumulative distribution function of the gamma distribution is

$$\begin{array}{c}\hfill F(x;\alpha ,\lambda ,\mu )=\frac{1}{\mathrm{\Gamma }(\alpha )}\gamma (\alpha ,\frac{x-\mu }{\lambda }),\end{array}$$(14)

where α, λ and μ are parameters respectively known as the shape, the scale, and the location parameter. The lower incomplete gamma function γ is defined for every complex number z with a strictly positive real part

$$\begin{array}{c}\hfill \gamma (z,x)={\displaystyle {\int }_0^x}{t}^{z-1}{e}^{-t}\mathrm{d}t.\end{array}$$(15)

Here Γ is the gamma function, defined as the limit of γ when x tends to infinity. The mixture model that we have applied is composed of two gamma distributions. The first describes the peak of the PDF(Δr_o), and will be referred to as the weak overshooting (WO) distribution (f_wo). The second describes the contribution of the heavy tail of the PDF(Δr_o), which corresponds to the most energetic events. This second distribution will be called the strong overshooting (SO) distribution (f_so). The PDF of the mixture model then reads

$$\begin{array}{c}\hfill f_{\mathsf{mix}}(\mathrm{\Delta }{r}_{\mathsf{o}})=w_{\mathsf{wo}}{f}_{\mathsf{wo}}(\mathrm{\Delta }{r}_{\mathsf{o}})+w_{\mathsf{so}}{f}_{\mathsf{so}}(\mathrm{\Delta }{r}_{\mathsf{o}}).\end{array}$$(16)

Here w_wo and w_so are weights of the WO and SO distributions, respectively, such that w_wo, w_so ≥ 0 and w_wo + w_so = 1.

We have written a fitting program which first calculates the CDF(Δr_o) of the data set and defines the mixture model given in Eq. (13) and Eq. (16) using w_so = 1 − w_wo. Our mixture model has seven free parameters: the location, shape and scale parameters for each of the two gamma distributions and the weight. We then employed a Levenberg–Marquardt algorithm (Levenberg 1944; Marquardt 1963) to find the best fit for the data. This is achieved using a non-linear least squares method to minimize the residual and find the optimal values for the 7 parameters of the mixture distribution. The Levenberg-Marquardt algorithm is an iterative procedure which requires an initial guess. We tested the robustness of this algorithm for the current problem by checking that a range of different initial guesses converge to the same solution. Application of this mixture distribution to the PDF(Δr_o) above the CZ of the 6 M_⊙ star is shown in Fig. 11. The mixture model produces a useful splitting of the original PDF into two gamma distributions. In each case, the WO distribution corresponds to approximately 85% to 95% of the overall PDF, confirming that the weakest overshooting events are more frequent than the more energetic ones. Because the WO distribution includes the measured overshooting of convective rolls that protrude from each side of the CZ, it makes physical sense that it contains a much larger portion of the probability space than the SO distribution. The mixture model therefore allowed us to filter the contribution of the weak overshooting plumes and model exclusively the events that most enhance the mixing mechanisms in the RZ.

Fig. 11. Probability density function (Δro) of a 6 M⊙ Cepheid above the convective shell and best fit with the mixture distribution.

3.3.2. Estimation of the overshooting length with the mixture model and comparison with the EVT

The mixture model uses the heavy tail of the PDF of all plumes to produce a full distribution for the strong overshooting events. The mean value that we extract from this distribution represents an average depth reached by strong overshooting plumes inside the RZ; we denote this depth by ℓ_so. It represents a characteristic length scale associated with enhanced diffusion within the RZ. We found a relatively large standard deviation around the expected value. These values are summarized in Table 3. Using both the mixture model and the EVT provides a detailed picture of the strong overshooting layer. Indeed, the EVT provides a maximal width of the overshooting layer.

Both the EVT and the mixture model yield length scales for overshooting that are larger than what is usually found for either convective envelopes or convective cores. These large length scales are consistent with the approximate width of the overshooting layer given by the averaged profile of the convective flux, as well as visualizations. For example, the visualization of vorticity in Fig. 3 clearly shows that convective rolls can be larger than the width of the CZ itself and can protrude on either side. Convective plumes are projected ballistically away from the convective boundaries with high energy due to the centrifugal force of the rotating vortices. These high energies appear to be increased by the opposite rotation of side-by-side convection rolls. This regime of convection differs significantly from large convective cores or envelopes.

In addition to the mean value, Table 3 includes the quantity ℓ_so^L, 90% (ℓ_so^U, 90%). This length scale is defined as the 90th percentile of the SO gamma distribution, the distance from the convection boundary in the lower (upper) overshooting layers, where 90% of the strongly overshooting plumes have already stopped. The 90% length scale is comparable to the maximal width of the overshooting layer provided by the EVT, and thus provides a link between these two statistical analysis methods. The 90% SO length scale reinforces the idea that the upper overshooting layer fills its radiative zone; it also provides a picture of how many plumes reach the outer convective envelope. 10% of the strongly overshooting plumes move beyond this length, and the strongly overshooting plumes are no more than 15% of the full distribution of overshooting plumes; thus an upper estimate is that 1.5% of overshooting plumes can reach and be mixed directly into the outer convective envelope of any of these Cepheids. By providing the 90% length scale alongside the mean value of the SO distribution, we have indicated a rough rate for which the number of overshooting plumes are able to mix each radial point within the radiative zone.

The ratio of the overshooting lengths for convective shells is uncertain, and is a measurement that can directly benefit stellar evolution calculations. Using 3D hydrodynamic simulations of massive stars, Cristini et al. (2017, 2019) and Rizzuti et al. (2022) found the ratio of lower overshooting length to upper overshooting length to be 0.2. For thin and shallow shells, the k − ω model of Guo & Li (2021) predicted a ratio of 0.5. These two numbers mark out a wide range of values.

To understand where the convective shells in Cepheids fall, we described the amount of convective boundary mixing both above and below the shell using our mixture model. To describe this, we examined the ratio of the lower overshooting length to the upper overshooting length, ℓ_so^L/ℓ_so^U. This ratio is close to 1 for shallow convective shells. We used the position of the lower convective boundary to determine how deep the convective shell lies within the star. The ratio of overshooting lengths decreases as the shell is located closer to the core of the star (see Fig. 12a).

Fig. 12. Ratio of the extent of the lower overshooting layer (ℓsoL) over the extent of the upper overshooting layer (ℓsoU) given by the mixture model, (a) versus the depth of the convective shell, measured as RCBL/R, and (b) versus the width of the convective shell measured as (RCBU−RCBL)/Hp, CBL. Error bars indicate one standard deviation around the mean value. A pink line indicates a linear regression fit.

The depth of the shell within the star also correlates with its width, measured as (R_CB^U−R_CB^L)/H_p, CB^L. In our study, the deepest convective shell is more than four times wider than the smallest, by this measure. The ratio of overshooting lengths decreases as the shell becomes thicker (see Fig. 12b). This makes sense because the density difference above and below a thin, internal convection shell is small; the difference between the RZ above and below the shell is also small. Moreover, for such thin shells, the CZ is filled with a line of convection rolls, producing a flow that is dominated by a single length scale. We thus expect the velocities of the convective rolls to be approximately the same at the bottom and the top of the CZ. When the shell is thicker, the overshooting ratio is smaller, due to the density difference between the lower and upper RZs.

The 8 M_⊙ Cepheid has a deeper shell than the 7 M_⊙ (b), but it is also thinner, demonstrating the combination of these effects. However, the small number of MESA models we selected for our hydrodynamic simulations differ more than just the mass and luminosity of the star. Because this work has included only high resolution simulations of 6 stars, it is unclear to what extent the depth or width of the shell influences the overshooting ratio independently. A more exhaustive parametric study would be required to fully characterize these thin and shallow convection zones and confirm these trends.

3.3.3. Determination of a diffusion coefficient for convection

Based on these results, we propose a new form of the diffusion coefficient used to model convective overshooting, designed for the complex situation of a radiative zone between two convective layers. Pratt et al. (2017) demonstrated that the diffusion coefficient can be expressed using probability distributions of overshooting plumes with the simple formula

$$\begin{array}{c}\hfill D_{\mathsf{EVT}}(r)=D_{\mathsf{MLT}}(1-F(r)).\end{array}$$(17)

In this expression, F(r) is the CDF as a function of the star’s internal radius; in Pratt et al. (2017) the form for the CDF was provided by EVT. The diffusion coefficient given by MLT is evaluated at the convective boundary, and defined D_MLT = 1/3 L_MLTv_MLT ; the MLT length scale L_MLT is proportional to the pressure scale height. Using the CDF of the mixture model, we propose the following form for mixing enhanced by convective overshooting

$$\begin{array}{cc}\hfill D_{\mathsf{shell}}& =D_{\mathsf{MLT}}[1-\frac{1}{\mathrm{\Gamma }\left({\alpha }_{\mathsf{s}}^{\mathsf{L}}o\right)}\gamma ({\alpha }_{\mathsf{s}}^{\mathsf{L}}o,\frac{(R_{\mathsf{C}}^{\mathsf{L}}B-r)/R-{\mu }_{\mathsf{s}}^{\mathsf{L}}o}{{\lambda }_{\mathsf{s}}^{\mathsf{L}}o})\hfill \\ \hfill & -\frac{1}{\mathrm{\Gamma }\left({\alpha }_{\mathsf{s}}^{\mathsf{U}}o\right)}\gamma ({\alpha }_{\mathsf{s}}^{\mathsf{U}}o,\frac{(r-R_{\mathsf{C}}^{\mathsf{U}}B)/R-{\mu }_{\mathsf{s}}^{\mathsf{U}}o}{{\lambda }_{\mathsf{s}}^{\mathsf{U}}o})].\hfill \end{array}$$(18)

The location parameter μ_so^U(L), the shape parameter α_so^U(L), and the scale parameter λ_so^U(L) of the SO distribution of the upper (lower) overshooting layer can be estimated from the simulations we performed in this work, and they give a range (in units of the pressure scale height at the convective boundary)

$$\begin{array}{ccc}& {\mu }_{\mathsf{s}}^{\mathsf{U}}o\approx 0.5-0.9H_{\mathsf{p},\mathsf{C}\mathsf{B}}^{\mathsf{U}},\hfill & \hfill {\mu }_{\mathsf{s}}^{\mathsf{L}}o\approx 0.7-1.0H_{\mathsf{p},\mathsf{C}\mathsf{B}}^{\mathsf{L}},\end{array}$$(19)

$$\begin{array}{ccc}& {\alpha }_{\mathsf{s}}^{\mathsf{U}}o\approx 2.0-5.0,\hfill & \hfill {\alpha }_{\mathsf{s}}^{\mathsf{L}}o\approx 1.5-3.0,\end{array}$$(20)

$$\begin{array}{ccc}& {\lambda }_{\mathsf{s}}^{\mathsf{U}}o\approx 0.2-0.4H_{\mathsf{p},\mathsf{C}\mathsf{B}}^{\mathsf{U}},\hfill & \hfill {\lambda }_{\mathsf{s}}^{\mathsf{L}}o\approx 0.15-0.25H_{\mathsf{p},\mathsf{C}\mathsf{B}}^{\mathsf{L}}.\end{array}$$(21)

Unlike the EVT analysis, neither of the location parameters used in the mixture model are directly equivalent to an overshooting length. This model for the diffusion coefficient predicts a smaller overshooting layer than the EVT, by definition (see Fig. 13). However, rather than filtering out all plumes except the one that penetrates the deepest, the new mixture model method produces a full distribution of overshooting plumes that matter, by filtering out only the weak overshooting plumes that do not. It assumes full mixing in this weak overshooting region. This allows for the amount of overshooting within a radiative layer between two convection zones to be modeled in a more detailed way.

Fig. 13. Diffusion coefficient for the 7 M⊙ (a) Cepheid. Heavy black lines indicate the convective boundaries of both the convective shell and the outer convection zone. The thin black dotted line is the diffusion coefficient of the convective shell based on the mixture model proposed in Eq. (18), and the red dashed line is the diffusion coefficient based on the EVT analysis (Eq. (17)). A gray line extrapolates our model to show a possible diffusion coefficient for the outer convection zone. The solid yellow line is a model for mixing in the radiative zone that experiences overshooting from both above and below, based on Eq. (20). In this example, DMLTenvelope = 0.8 DMLT.

Our simulations did not include the outer convective envelope. However, our overshooting diffusion coefficient based on the mixture model allows us to consider the combined effects of the convective shell and the convective envelope on the radiative layer between them. The EVT analysis predicts full mixing between these two convective layers, creating an outer envelope that is proportionally at least as large as the current sun for these Cepheids. Our new mixture model in Eq. (18) provides a detailed picture of how such overshooting layers could overlap. We have assumed that, for the outer convective envelope, the mixture model would produce a diffusion coefficient of the same form as for the lower overshooting layer of the convective shell, taking into account the different pressure scale height. We propose the following form of the diffusion coefficient for the full star, denoted D_MM:

$$\begin{array}{c}\hfill D_{\mathsf{MM}}=D_{\mathsf{shell}}+D_{\mathsf{envelope}}^{\mathsf{L}},\end{array}$$(22)

where the additional diffusion coefficient is defined as

$$\begin{array}{c}\hfill D_{\mathsf{envelope}}^{\mathsf{L}}=D_{\mathsf{ML}}{T}^{\mathsf{envelope}}[1-\frac{1}{\mathrm{\Gamma }\left({\alpha }_{\mathsf{s}}^{\mathsf{E}}o\right)}\gamma ({\alpha }_{\mathsf{s}}^{\mathsf{E}}o,\frac{(R_{\mathsf{C}}^{\mathsf{E}}B-r)/R-{\mu }_{\mathsf{s}}^{\mathsf{E}}o}{{\lambda }_{\mathsf{s}}^{\mathsf{E}}o})].\end{array}$$(23)

For the purposes of exploring mixing throughout this radiative zone, we assumed μ_so^E = μ_so^L, α_so^E = α_so^L and λ_so^E = λ_so^L and that D_MLT^envelope = 0.8 D_MLT.

With these parameters, the mixing in the radiative zone between two convective layers can be calculated (an example is shown in Fig. 13). In the figure, the two overshooting layers overlap. Close to the upper convective boundary, we observe that D_MM/D_MLT can be greater than 1.0 in some regions. This finding introduces the idea that a “super-mixing layer” could exist between two convection zones that are relatively close to each other; such super-mixing layers could have more intense mixing than in the neighboring convection zones. This idea is consistent with the concept of diffusion that is enhanced by the local properties of convection. Hydrodynamic simulations with high radial resolution, which contain both convective layers would be required to confirm the existence of these “super-mixing layers” and more fully explore the statistics necessary for this overshooting model.

4. Summary and conclusions

We have studied the internal structure of intermediate-mass stars that lie in the instability strip, and are therefore Cepheid variable stars. We produced a range of stellar models using MESA and documented the presence of a thin, shallow interior convective shell as the star crosses the instability strip. Such inner convective shells have been studied only rarely using multidimensional hydrodynamic simulations. We have performed hydrodynamic simulations in spherical shells that cover 40%–60% of the stellar radius, for stellar structures corresponding to masses in the range of 5.75–9 M_⊙. We have studied stellar convection and convective boundary mixing in these simulations using statistical methods. The statistical results that we produce are based on data sets that cover long periods of stellar convection for each star.

The volume-percentage filling factor of inflowing plumes reveals only a slight asymmetry between inflows and outflows in these convective shells. The filling factor has a sinusoidal shape for each simulation, indicating larger inflows in the upper half of the convective shell while outflows are bigger in the lower half. This is opposite to the general picture seen in convective envelopes, where inflows are thicker at the bottom of the convection zone. Although this trend is significant, the filling factor remains relatively close to 0.5, and the sinusoidal shape has an amplitude that makes the filling factor deviate away from perfectly symmetric inflows and outflows by approximately ∼10%. These results are in agreement with the analysis of the average width of inflows. In the unstable region, large convection rolls that fill the width of the layer dominate the flow, and we observe a nearly single-scale convective structure. However, boundary layer flows are strong, and closer to the Schwarzschild boundaries we observe thinner inflows. In both of the overshooting layers, the widths of plumes reach a distinct minimum.

To define the extent of the overshooting layers, we have used the vertical kinetic energy flux. Historically, a simple time and spatial average has often been used to determine the depth of the penetration layer. However, as Pratt et al. (2017) observed for pre-main sequence stars with large convective envelopes, in the stars we study, the PDFs of all penetrating plumes are non-Gaussian. Two distinct layers can be identified in the radiative zone: a shallow layer in which the weakest convective plumes overshoot frequently, and a deeper layer where strong overshooting plumes overshoot more rarely. The weak overshooting layer is characterized by the presence of convection rolls that exceed the convective boundaries before returning back in the CZ, an unusual regime of stellar convection that differs from both convective cores and envelopes. The strong overshooting plumes are responsible for the enhanced diffusion mechanisms in the radiative zone. We therefore applied the extreme value theory model to determine the maximal extent of the overshooting layer. However, although the convective shells in Cepheids are thin, overshooting plumes impact the full extent of the radiative layer above the convective shell and below the outer convective envelope. The crudest picture of these stars is thus one where there is a large outer convective envelope that is fully mixed by convection; this picture, however, conceals interesting and relevant details about how much the convection enhances mixing at a given radius in these stars. The EVT statistical analysis does not provide sufficient details to allow an unprejudiced comparison between the upper and the lower overshooting layers for the convective shell.

This has motivated us to propose an alternative statistical analysis in which weakly and strongly overshooting plumes are separated. We then modeled the PDF of all overshooting plumes using a mixture of gamma distributions. This has allowed us to filter the contribution of the weak overshooting layer from that of the strong overshooting layer. We find that the intermittent, strong overshooting plumes correspond to 5–15% of all plumes and we characterized the strong overshooting distribution by its mean value and standard deviation. We defined a second characteristic length scale as the 90th percentile of the SO distribution; this length scale is comparable with the length scale obtained with the EVT, reinforcing the idea that the upper overshooting layer fills the entire radiative zone. We demonstrated that for the present simulations, the overshooting lengths thus obtained are consistent with the EVT analysis below the convection zone. This approach to measuring the extent of convective boundary mixing is new and further investigation would be required to demonstrate its utility for describing convective overshooting in simulations that include additional physical effects such as rotation, magnetic fields, and pulsations.

We have used the mixture model to determine a ratio of the overshooting lengths for the overshooting layers on each side of the convective shell. For our thin and shallow convective shells, in which one length scale of convective structures dominates, the difference of overshooting above and below the shell is less than what has previously been reported in the literature. As we examine convective shells that are thicker, and deeper inside the star, the ratio of lower overshooting length to upper overshooting length drops. This ratio is always less than one for the simulations we have studied; however the relatively large standard deviation associated with the mixture model allows for the possibility that at some particular instant, overshooting could be greater in the lower overshooting layer than in the upper overshooting layer.

Based on the mixture model, we have proposed a diffusion coefficient for overshooting from a convective shell, as well as a form of a diffusion coefficient for the full star that provides a more detailed picture of the enhanced mixing mechanisms in the RZ between two convective layers. Although further investigation would be necessary to fully characterize the complex mixing in overlapping overshooting layers, both the EVT and our model point out the possibility of a “super-mixing layer” in which more efficient mixing could result in the merging of two adjacent convective zones.

The stellar structures explored in the hydrodynamic simulations in this work were not produced to systematically explore the parameter space of mass or luminosity. Furthermore, we have found no clear differences other than those related to the depth and width of the convective shell between the 7 M_⊙ (a) and (b) Cepheids, which are located respectively in the second and third crossing of the instability strip. In future work, a more systematic study would be required to assess a trend in the decrease of the ratio of overshooting with the depth and width of the convection zone. The present simulations are 2D, and 3D simulations would be desirable, as long as those simulations would have a radial resolution equivalent to those studied here. A 3D simulation of shell convection might allow for the study of instabilities of convection rolls (e.g., Busse & Clever 1979; Clever & Busse 1989) in the regime relevant to stellar interiors; such instabilities could contribute to a more complete picture of convective mixing. A better understanding of the nonlinear interactions between convection and radial pulsations in Cepheid variable stars remains necessary to deepen our theoretical knowledge of these physically complex stars; 3D simulations of the whole star, as well as simulations that include rotation and magnetic fields are necessary to expand on these results. Those directions are ones we intend to pursue.

Acknowledgments

The research leading to these results is partly supported by the ERC grants 320478-TOFU and 787361-COBOM and by the STFC Consolidated Grant ST/Y002164/1. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (http://www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. This work also used the University of Exeter local supercomputer ISCA. Computing support for this work came from the Lawrence Livermore National Laboratory (LLNL) Institutional Computing Grand Challenge program. Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-JRNL-2003963

References

Appendix A: The structure of classical Cepheid variable stars

We used the MESA code (version 24.03.1) to understand typical structures for the interior of classical Cepheids. Our study is not exhaustive, and we limit the range of free parameters in our MESA models. We generated a grid of 48 MESA models with six different metallicities, four different masses, and two different overshooting parameters. This study is targeted toward establishing common structures for Cepheids with parameters documented in the literature.

Classical Cepheids are intermediate mass stars that have evolved beyond the main sequence and are grouped together in the so-called instability strip of the Hertzprung–Russell (HR) diagram (indicated by the shaded areas in Fig. A.1; see, e.g., Bono et al. 2005; Anderson et al. 2014). The evolutionary tracks of classical Cepheids first rapidly cross this instability strip while the hydrogen shell is still burning. When they begin to burn helium in their core, they typically execute a blue loop, crossing the instability strip a second and third time. Stellar structure and evolution models of intermediate mass stars that may evolve to Cepheids are sensitive to both the input parameters (including overshooting, convection, metallicity, and nuclear reaction rates) and the numerical scheme (see, e.g., the recent study of Ziółkowska et al. 2024). The occurrence and duration of blue loops produced by stellar evolution models can be impacted by the choice of input parameters. However the duration of the blue loop tends to be smaller for lower mass stars (see, e.g., Walmswell et al. 2015). Blue loops for stars M < 4 solar masses (denoted M_⊙ thereafter) that are produced by current versions of the MESA code do not cross the instability strip (Anderson et al. 2016; De Somma et al. 2022; Espinoza-Arancibia et al. 2024). For our MESA study, we therefore have focused on stars in the mass range from 6 to 9M_⊙.

Fig. A.1. Hertzprung–Russell diagrams of 6, 7, 8, and 9 M⊙ stars of six different metallicities. The symbols indicate structures of individual Cepheid variable stars that lie in the instability strip (pink shaded area) and that we examine further in Appendix A. The blue and red edges of the instability strips were estimated from the work of Anderson et al. (2016) for z = 0.002, z = 0.006 and z = 0.014, from the work of Alibert et al. (1999) for z = 0.004, and from the theoretical relationships of Bono et al. (2000) for z = 0.008 and z = 0.020. L is the luminosity of the star, expressed in units of the luminosity of the sun L⊙ (we use L⊙ = 3.828 × 1033erg.s−1), and Teff denotes its effective temperature. Evolutionary tracks with different masses are distinguished by color.

Fig. A.2. Schwarschild discriminant as a function of stellar radius for the 6, 7, 8, and 9M⊙ stars presented in Fig. A.1. The radius r is given in units of the total radius of each star, R. The arrow points in the direction of increasing metallicity. Profiles with different metallicities are distinguished by color.

Based on inlists published by Espinoza-Arancibia et al. (2022), we have produced evolutionary tracks in this mass range for a range of stellar metallicities between z = 0.002 and z = 0.020. This range of metallicities has been documented in various Cepheid populations (see, e.g., Baraffe et al. 1998; Anderson et al. 2016; Chiosi et al. 1992; Alibert et al. 1999; Hocdé et al. 2024). For each mass and metallicity, the helium mass fraction y = 0.256 and we have used a mixing length parameter of α_mlt = 1.88. Exponential core and envelope overshooting is used at the convective boundaries. For half of the evolutionary tracks produced, we set the overshooting parameter to f_ov = 0.009; for the other half we set f_ov = 0.019 so that we could examine the effect of overshooting in a binary way. For these parameters, for each mass and metallicity, we always find at least one evolutionary track that executes a blue loop large enough to cross the instability strip twice. The evolutionary tracks of 24 of these stars are plotted in Fig. A.1 with symbols indicating selected models of Cepheids that we examine in detail below. The blue and red edges of the instability strips were estimated from the work of Anderson et al. (2016) for z = 0.002, z = 0.006 and z = 0.014, from the work of Alibert et al. (1999) for z = 0.004, and from the theoretical relationships of Bono et al. (2000) for z = 0.008 and z = 0.020.

Fig. A.3. Radial profiles of opacity (κ) for the 6, 7, 8, and 9M⊙ stars presented in Fig. A.1. The arrow points in the direction of increasing metallicity. Profiles with different metallicities are distinguished by color.

Fig. A.4. Cepheids for an 8M⊙ star that has metallicity z = 0.020 (also shown in Fig. A.1): (a) the H-R diagram showing that this star crosses the instability strip three times. Symbols denote three stellar models corresponding to the Schwarzschild discriminants in (b).

Fig. A.5. Location and width of the convective shell for each of the three crossings of the instability, for the 6, 7, 8, and 9M⊙ stars at z = 0.020.

Because our interest is stellar convection, we have focused on the presence of convectively unstable zones in the interior of Cepheids, and their dependence on mass and metallicity. To characterize whether a layer of the stellar interior is stable against convection, we used the Schwarzschild criterion (Schwarzschild & Kiess 1958). For convectively unstable layers, the Schwarzschild discriminant, defined as the difference between the adiabatic and the radiative temperature gradients, ∇_ad − ∇_rad is negative. We first examined several stars that are crossing the instability strip for the second time. On the evolutionary tracks in Fig. A.1, several points are indicated; for the stellar structures at these points, we have plotted the Schwarzschild discriminant against the star’s internal radius (see Fig. A.2). For each mass and metallicity, these Cepheids exhibit a convective core that is so small that the curve visibly merges with the left axis, and a small outer convective envelope just below the photosphere. For the 8M_⊙ Cepheid with z = 0.020, we found that the convective boundary of the convective core is located at 0.65% of the stellar radius, but represents 9% of the stellar mass. In contrast, the outer convection zone corresponds to the outer 5% of the star’s radius.

The Schwarzschild discriminant also exhibits a dip in the interior of the star, in the range 0.4 < r/R < 0.7 (see Fig. A.2). For sufficiently high metallicity, this dip corresponds to the appearance of a convectively unstable layer. The position of this dip corresponds to the so-called Z-bump (Iglesias et al. 1992; Rogers & Iglesias 1992; Seaton 1993a; Seaton et al. 1994; Daszynska-Daszkiewicz et al. 2019) in the opacity profiles (see Fig. A.3). The Z-bump is a local maximum in the Rosseland mean opacity that occurs around a temperature of 2 × 10⁵K; this phenomenon is caused by a large number of bound-bound but also bound-free transitions in iron-group elements. The metallicity of the star influences the width and depth of this convection zone; as the metallicity decreases, this convective shell becomes thinner until it completely disappears at the lowest metallicities that we consider (see Fig. A.2). Based on the grid of MESA models that we have produced, we found that this convective shell is wider and deeper for more massive stars. Above a minimum mass, the convective shell also exists for lower metallicities. Thus, for the 6 and 7M_⊙ stars, this convectively unstable region exists close to solar abundance and vanishes at lower z; the Schwarzschild discriminant is positive and the opacity Z-bump is less pronounced in these cases. The mass of the star affects the location of the convective shell, although the metallicity appears to have only a limited influence on it. As the mass of the Cepheid increases, the convective shell is shifted deeper in the star. Although further investigation would be required to fully characterize this inner convective shell, the combination of a tiny convective core, an inner convective shell, and an outer convective envelope is a common structure for Cepheids on their second crossing of the instability strip, provided that they have sufficiently high mass and metallicity.

For the first and third crossings of the instability strip (see Fig. A.4) we find similar results. During all three crossings, a convective shell appears in the Z-bump region, (see Fig. A.4 for a comparison of structures for the three crossings of a 8M_⊙ star). For stars that are on the first and third crossings of the instability strip, the same dependence on mass and metallicity exists; this convective shell becomes thinner as the mass decreases, and disappears at low mass and low metallicity (see Fig. A.5). This comparison of the width and location of the convective shell for the 6, 7, 8 and 9M_⊙ stars at z = 0.020, for each of the three crossings of the instability strip, exhibits a clear trend. No convective shell exists for the 6M_⊙ star during the first crossing, even at the highest metallicity studied, but a convective shell appears during this star’s second and third crossings of the instability strip. The location of the unstable region slightly changes as one moves along the evolutionary track in the instability strip, and it is shallower and smaller at the first crossing compared with the second and third crossings. As a consequence, the convective shell tends to disappear at higher metallicities than for the last two crossings. As illustrated in Figs. A.4 and A.5, few differences are observed between the second and third crossings. As we consider stars of higher mass, the convective shell is both wider and deeper, for all three crossings of the instability strip.

All Tables

All Figures

Fig. 1. Schematic of the typical structure of the Cepheids we study in this work with MUSIC. Convective regions are colored, including: a tiny convective core, an inner convective shell surrounded by radiative zones, and a thin outer convective envelope. The simulation domain is a spherical shell, and is outlined with black; it encapsulates the interior convection zone along with the surrounding radiative zones, and is truncated just below the outer convective envelope.

In the text

	Fig. 2. (a) Hertzprung-Russel diagram and (b) radial profile of the Schwarzschild discriminant, for the six stars we simulated with MUSIC. Symbols indicate the Cepheids we simulated, and the pink shaded area indicates the instability strip for the evolutionary tracks. Shaded regions also highlight the internal convection zones around which our hydrodynamic simulations are centered.
In the text

Fig. 3. Visualizations of (from left to right) vorticity magnitude, radial velocity, and angular velocity of the 8 M⊙ simulation after 30τconv of steady-state convection. The zero point for vorticity is in dark red in the left panel. In the center and right panels, outwards flows are in red while inwards flows are in yellow; the zero point in velocity is black. The radial velocity ranges from −9.25 × 105 cm/s to 9.25 × 105 cm/s at this instant in time, while the angular velocity ranges from −1.01 × 106 cm/s to 1.01 × 106 cm/s. Dashed white lines indicate the position of the Schwarzschild boundaries.

In the text

	Fig. 4. Radial profile of the time-averaged RMS velocity scaled by (Lstar/1035)1/3 for a range of stellar masses. Vertical lines indicate the Schwarzschild convective boundaries for each star.
In the text

	Fig. 5. Radial profile of the time-averaged volume percentage filling factor of the inward moving plumes σvp, in versus the star’s internal radius, in units of the total stellar radius R for the 6 M⊙ simulation. The shaded area indicates one standard deviation above and below this averaged line and thin vertical lines indicate the Schwarzschild boundaries.
In the text

	Fig. 6. Width of inflowing plumes Win for the 5.75 and the 9 M⊙ simulations as a function of the star’s radius. Shaded areas indicate one standard deviation above and below the time-averaged line. Thin vertical lines indicate the radial position of the Schwarzschild boundaries.
In the text

Fig. 7. (a) Angular structure of the overshooting layer after 114τconv of steady-state convection in the simulation of the 6 M⊙ Cepheid. The overshooting length in this illustration is determined by zeros of the vertical kinetic flux. The Schwarzschild convective boundaries between the convection zone and the neighboring stable radiative zones are indicated by a solid black line. The stellar radius (r/R) is indicated on the left axis, and the overshooting length is given in units of the pressure scale height Hp, CB at the convective boundaries on the right axis. A dashed line indicates the average overshooting length at this time. (b) The same, but plotted over a visualization of vorticity. White dashed lines indicate the Schwarzschild boundaries.

In the text

	Fig. 8. Radial profiles of the convective flux for the 6 M⊙ Cepheid. The shaded area indicates one standard deviation around the mean value. Thin vertical lines indicate the convective boundaries determined by the Schwarzschild criterion. Black dashed lines indicates one pressure scale height above and below the convection zone. The inset shows a zoomed in region centered on the red box.
In the text

	Fig. 9. Probability density functions of ro, the extent of all overshooting plumes, for the 6 and the 7 M⊙(b) Cepheid. Thin vertical lines indicate the convective boundaries determined by the Schwarzschild criterion.
In the text

	Fig. 10. Natural logarithm of the negative natural logarithm of the cumulative distribution functions of maximal overshooting length, Δrmax, for each simulation below the convective shell. Black lines indicate the best fit with the Weibull distribution, defined in Eq. (11).
In the text

	Fig. 11. Probability density function (Δro) of a 6 M⊙ Cepheid above the convective shell and best fit with the mixture distribution.
In the text

Fig. 12. Ratio of the extent of the lower overshooting layer (ℓsoL) over the extent of the upper overshooting layer (ℓsoU) given by the mixture model, (a) versus the depth of the convective shell, measured as RCBL/R, and (b) versus the width of the convective shell measured as (RCBU−RCBL)/Hp, CBL. Error bars indicate one standard deviation around the mean value. A pink line indicates a linear regression fit.

In the text

Fig. 13. Diffusion coefficient for the 7 M⊙ (a) Cepheid. Heavy black lines indicate the convective boundaries of both the convective shell and the outer convection zone. The thin black dotted line is the diffusion coefficient of the convective shell based on the mixture model proposed in Eq. (18), and the red dashed line is the diffusion coefficient based on the EVT analysis (Eq. (17)). A gray line extrapolates our model to show a possible diffusion coefficient for the outer convection zone. The solid yellow line is a model for mixing in the radiative zone that experiences overshooting from both above and below, based on Eq. (20). In this example, DMLTenvelope = 0.8 DMLT.

In the text

Fig. A.1. Hertzprung–Russell diagrams of 6, 7, 8, and 9 M⊙ stars of six different metallicities. The symbols indicate structures of individual Cepheid variable stars that lie in the instability strip (pink shaded area) and that we examine further in Appendix A. The blue and red edges of the instability strips were estimated from the work of Anderson et al. (2016) for z = 0.002, z = 0.006 and z = 0.014, from the work of Alibert et al. (1999) for z = 0.004, and from the theoretical relationships of Bono et al. (2000) for z = 0.008 and z = 0.020. L is the luminosity of the star, expressed in units of the luminosity of the sun L⊙ (we use L⊙ = 3.828 × 1033erg.s−1), and Teff denotes its effective temperature. Evolutionary tracks with different masses are distinguished by color.

In the text

	Fig. A.2. Schwarschild discriminant as a function of stellar radius for the 6, 7, 8, and 9M⊙ stars presented in Fig. A.1. The radius r is given in units of the total radius of each star, R. The arrow points in the direction of increasing metallicity. Profiles with different metallicities are distinguished by color.
In the text

	Fig. A.3. Radial profiles of opacity (κ) for the 6, 7, 8, and 9M⊙ stars presented in Fig. A.1. The arrow points in the direction of increasing metallicity. Profiles with different metallicities are distinguished by color.
In the text

	Fig. A.4. Cepheids for an 8M⊙ star that has metallicity z = 0.020 (also shown in Fig. A.1): (a) the H-R diagram showing that this star crosses the instability strip three times. Symbols denote three stellar models corresponding to the Schwarzschild discriminants in (b).
In the text

	Fig. A.5. Location and width of the convective shell for each of the three crossings of the instability, for the 6, 7, 8, and 9M⊙ stars at z = 0.020.
In the text