3D hydrodynamics simulations of a 3 M ⊙ core-helium burning star

The inner structure of core-helium burning (CHeB) stars remains uncertain due to the yet unknown nature of mixing at the boundary of their cores. Large convective cores beyond a bare Schwarzschild model are favoured both from theoretical arguments and from asteroseismological constraints. However, the exact nature of this extra mixing, and in particular the possible presence of semiconvective layers, is still debated. In this work, we approach this problem through a new avenue by performing the first full-sphere 3D hydrodynamics simulations of the interiors of CHeB stars. We use the PPMstar explicit gas dynamics code to simulate the inner 0.45 𝑀 ⊙ of a 3 𝑀 ⊙ CHeB star. Simulations are performed using different Cartesian grid resolutions (768 3 , 1152 3 and 1728 3 ) and heating rates. We use two different initial states, one based on MESA ’s predictive mixing scheme (which yields a large overshoot region) and one based on the convective premixing approach (which exhibits a semiconvective interface). The general behaviour of the flow in the convective core and in the stable envelope (where internal gravity waves are observed) is consistent with our recent simulations of core convection in massive main-sequence stars, and so are the various scaling relations. The semiconvective layers are dominated by strong internal gravity waves that do not produce measurable species mixing, but overshooting motions from the convective core gradually homogenize the semiconvective interface. This process can possibly completely erase the semiconvective layers, which would imply that CHeB stars do not harbour a semiconvection zone.


INTRODUCTION
Core helium burning (CHeB) stars are characterized by a central convective He-burning core surrounded by a convectively stable Herich envelope.Observationally, CHeB stars are known as red clump stars, secondary clump stars, RR-Lyrae, horizontal branch stars, or subdwarfs B. These various classes of CHeB stars descend from different evolutionary pathways, but they all have a He-fusing core where the triple- reaction produces carbon and 12 C(, ) 16 O makes oxygen.
In low-and intermediate-mass stars, the treatment of convective boundary mixing (CBM) at the edge of the He-fusing core is particularly challenging for 1D stellar evolution (e.g., see the review by Salaris & Cassisi 2017).The generation of C and O inside the convective core enhances its opacity , and the radiative gradient, increases with time in the convection zone.In Equation (1),  is the radiation constant,  is the speed of light,  is the gravitational constant,  is the luminosity,  is the pressure,  is the mass enclosed within the radius at which ∇ rad is calculated, and  is the temperature.This ∇ rad increase leads to the formation of a discontinuity, such that ∇ rad > ∇ ad just inside the core and ∇ rad = ∇ ad just outside (blue ★ E-mail: sblouin@uvic.casolid line in Figure 1), where ∇ ad is the adiabatic temperature gradient.By definition, the convective boundary is located at the radius where radiation carries out all the energy.Within the mixing length theory (MLT) framework used in many 1D stellar evolution calculations, this condition is satisfied where the Schwarzschild criterion ∇ rad = ∇ ad is met (Biermann 1932).It is of course well established that convective penetration can significantly shift the location of the convective boundary with respect to the Schwarzschild boundary (for recent multi-dimensional simulations demonstrating this point, see Käpylä 2019;Anders et al. 2022a;Andrassy et al. 2023; Baraffe et al. 2023;Blouin et al. 2023b;Mao et al. 2023), but it is instructive in what follows to consider the consequences of applying the standard Schwarzschild criterion to the He-burning core.
Within this framework, ∇ rad = ∇ ad must be satisfied on the convective side of the boundary, since otherwise the convective flux cannot be zero (Gabriel et al. 2014).In CHeB stars, we encounter in 1D stellar evolution a situation where ∇ rad = ∇ ad on the radiative side of the ∇ rad discontinuity, but where ∇ rad > ∇ ad on the convective side.Therefore, the convective flux is not zero at that location, and this cannot be the convective boundary.This conclusion can also be reached by realizing that the Schwarzschild boundary is in an unstable equilibrium (Schwarzschild & Härm 1969;Castellani et al. 1971a).If the layer just above the convection zone is mixed with the core, its opacity increases due to the inflow of opaque C/O-rich material, ∇ rad surpasses the adiabatic gradient, and the convective core grows.Mixing of the layer immediately above the core is inevitable: any amount of convective overshooting can accomplish that, and even atomic diffusion alone could be sufficient (Michaud et al. 2007).
To find a stable boundary, the convective core must be extended until ∇ rad = ∇ ad on the convective side of the Schwarzschild boundary, but this leads to further complications.Given the behaviour of the opacity and the thermodynamic structure of CHeB stars, extending the core leads to the formation of a local minimum in ∇ rad within the convective region.Eventually, this local minimum in ∇ rad reaches ∇ ad , effectively splitting the convective region in two.What does that mean for mixing in the star?Full mixing of the gap now separating the convective core from the convective shell above is problematic as it would imply the formation of a region with ∇ rad < ∇ ad at the local minimum, in contradiction with the assumption of full mixing.
Two families of solutions have been proposed to solve this problem within the 1D MLT framework.The first consists of the formation of a partially mixed (or semiconvective) region between the ∇ rad local minimum and the radiative envelope (Schwarzschild & Härm 1969;Paczyński 1970;Castellani et al. 1971b).In this scenario, an extended semiconvective region where ∇ rad = ∇ ad = ∇ T (∇ T is the actual temperature gradient) separates the convective core from the stable layers (grey dash-dotted line in Figure 1).This CBM prescription has been implemented in many stellar evolution codes (e.g., Dorman & Rood 1993;Cassisi et al. 2003;Constantino et al. 2015), including in MESA through the convective premixing scheme (CPM, Paxton et al. 2019;Ostrowski et al. 2021).Another solution consists of artificially halting the growth of the convective core before the splitting occurs.This is known as maximal overshoot (Constantino et al. 2015;Denissenkov et al. 2017;Li et al. 2018, in the case where the growth is stopped just before the point where splitting would take place) or predictive mixing (PM, Paxton et al. 2018;Ostrowski et al. 2021).1In this scenario, there is no partially mixed zone and, problematically, ∇ rad remains greater than ∇ ad on the inner side of the convective boundary (orange dashed line in Figure 1).
It is unclear which of these two approaches should be preferred in 1D stellar evolution calculations, but some insights can be gained from observational constraints.Asteroseismological studies of various flavours of CHeB stars have clearly established that a bare Schwarzschild core (i.e., no extra CBM, as in the blue solid line of Figure 1) can be ruled out.A more extended convective core is definitely required to reproduce the observed pulsation periods (Van Grootel et al. 2010a,b;Charpinet et al. 2011Charpinet et al. , 2019;;Montalbán et al. 2013;Bossini et al. 2015;Constantino et al. 2015;Uzundag et al. 2021).A similar conclusion is reached via cluster star counts, which are used to infer the lifetimes of horizontal branch stars (Constantino et al. 2016).But beyond this finding, CHeB asteroseismological constraints do not yet offer a clear answer regarding which CBM scheme should be adopted.Constantino et al. (2015Constantino et al. ( , 2016) ) conclude that both the maximal overshoot and semiconvective prescriptions can be compatible with observations of horizontal branch stars.Similarly, for subdwarfs B, Uzundag et al. (2021) find that both the PM (or maximal overshoot) and CPM (or semiconvective) schemes implemented in MESA produce core masses that are in better agreement with the observations than bare Schwarzschild models (although even then the core sizes remain below the seismically derived values).Another promising observational window is white dwarf asteroseismology.
The composition profiles of white dwarf C-O cores bear the imprint of the CHeB phase (Straniero et al. 2003;Salaris et al. 2010;Chidester et al. 2023), and empirically derived white dwarf internal stratifications provide valuable constraints on CHeB CBM.Giammichele et al. (2018Giammichele et al. ( , 2022) ) have mapped the core compositions of a handful of white dwarfs and found large O-rich central regions that require extra CBM during the CHeB phase compared to standard evolution models (De Gerónimo et al. 2019).
Better understanding CBM in CHeB stars would have important ramifications beyond CHeB stars.In particular, extra mixing during the CHeB phase ultimately leads to the formation of white dwarfs with more O-rich cores, and the exact O abundance profile of white dwarfs is a key determinant of their cooling rates.Not only does it determine the total thermal energy content of the star (Fontaine et al. 2001), but it also controls how much energy is released by fractionation processes during core crystallization (Segretain et al. 1994;Montgomery et al. 1999;Salaris et al. 1997Salaris et al. , 2000;;Althaus et al. 2012;Blouin et al. 2020Blouin et al. , 2021)).Due to current uncertainties related to CBM during the CHeB phase, the O abundance profile of white dwarfs remains poorly constrained (Salaris et al. 2010(Salaris et al. , 2022)), and this injects systematic errors in the white dwarf cooling models (Renedo et al. 2010;Bédard et al. 2020;Salaris et al. 2022;Bauer 2023) used in diverse age-dating applications (Hansen et al. 2013;Fantin et al. 2019;Boylan-Kolchin & Weisz 2021;Kaiser et al. 2021;Cimatti & Moresco 2023).For the oldest white dwarfs in the Milky Way, CHeB CBM uncertainties result in errors of the order of a Gyr on inferred ages.
In this work, we use full-sphere 3D stellar hydrodynamics simulations of the interiors of CHeB stars to confront the CHeB CBM problem from a yet unexplored angle.Our simulations, performed with the PPMstar explicit gas dynamics code, follow for dozens of convective turnover timescales the 3D hydrodynamics response of the gas to a given CHeB structure.We study the behaviour of the flow for the two commonly used 1D CBM prescriptions, MESA's PM and CPM schemes.We describe these MESA models in Section 2, where we also explain how our PPMstar simulations are initialized.The general properties of the 3D simulations (flow morphology, convergence with respect to grid resolution, scaling relations) are explored in Section 3. In Section 4, we study mixing and entrainment in our simulations and discuss the implications of our findings in the context of the CHeB CBM problem.Finally, we conclude in Section 5.

METHODS
In this section, we first describe the calculations we have performed with the MESA code to initialize our PPMstar simulations.In Section 2.2 we then explain how the MESA models were mapped into our 3D simulations.Finally, we detail the PPMstar simulations themselves in Section 2.3.

MESA→PPMstar mapping
As in our recent PPMstar simulations of various types of stars (Blouin et al. 2023a,b;Herwig et al. 2023;Mao et al. 2023), the initial 3D state is reconstructed from the MESA entropy () and mean molecular weight () profiles.The PPMstar base state is calculated using these profiles and enforcing hydrostatic equilibrium.This integration is performed using PPMstar's equation of state, which guarantees that our initial state is precisely in hydrostatic equilibrium.Note that Age (Myr) small-scale noise was filtered out from the raw MESA  and  profiles to avoid injecting spurious small-scale structures into the 3D base state.In addition, we smoothed the  and  profiles at the outer edge of the semiconvection zone (in the CPM case) and of the convective core (in the PM case) as the CPM and PM schemes each produce unphysical discrete jumps.These transition regions were flattened so that they span at least 10 grid cells in our PPMstar simulations.This procedure is shown in the top panels of Figures 3 and 4.
The equation of state currently implemented in PPMstar includes the pressure contributions from the ideal gas and from the radiation field, where  is the ideal gas constant.In the central layers of CHeB stars, additional contributions from electron degeneracy pressure and ionion nonideal interactions come into play.As a result there is a ≃ 10% mismatch between the central pressure in our MESA models and that in our PPMstar setups (see third panels of Figures 3 and 4).A difference of that magnitude is not expected to impact the dynamics of the simulations in any meaningful way.Indeed, in our previous work on massive main-sequence stars, we have found the neglect of the radiation pressure term (which accounts for ≃ 20% of the pressure in the 25  ⊙ star we have studied) to be an excellent approximation (Herwig et al. 2023;Mao et al. 2023).However, we will see in Section 4 that this mismatch complicates the interpretation of the observed migration of the convective/semiconvective boundary over the course of our simulations.
Our PPMstar simulations include radiative diffusion via a radiative flux term in the energy equation (Mao et al. 2023), To model the opacity , we use the OPAL tables (Iglesias & Rogers 1996).Since table look-ups would be too inefficient for a highly optimized code like PPMstar, we have built a polynomial expression that reliably approximates the OPAL tables in the limited composition-temperature-density space explored in our simulations (Appendix A).As shown in the bottom panels of Figures 3 and 4, this simple prescription reliably captures the dependence of  on the composition and physical conditions.The green dotted lines obtained with this opacity model closely match the opacities calculated by MESA (blue solid lines).The opacity profile actually used in our simulations (orange dashed lines) departs more significantly from the MESA profiles.This is caused by the differences between the ther-  modynamic structures of the star in MESA and in PPMstar due to the incomplete equation of state.

PPMstar simulations
PPMstar is an explicit gas dynamics code where the conservation equations are solved on a 3D Cartesian grid (Woodward et al. 2015(Woodward et al. , 2018(Woodward et al. , 2019;;Jones et al. 2017;Andrassy et al. 2020;Herwig et al. 2023).In our simulations, the nuclear energy source from He burning in the core is modelled by a constant volume heating following a Gaussian radial profile that matches the shape of the MESA heating profile (our use of boost factors, described below, implies that we do not match its magnitude).Two fluids are included in the calculations: one having the mean molecular weight of the C/O-rich core ( core = 1.5845) and one having the mean molecular weight of the almost pure-He envelope ( env = 1.3359).
All our simulations are performed with heating luminosities that exceed the nominal He burning luminosity  ★ of the star.MLT predicts a convective Mach number smaller than 10 −4 in the Heburning cores of our 3  ⊙ stars.As PPMstar is an explicit gas dynamics code, accurately resolving such slow flows would demand prohibitively small simulation grid cells.To circumvent this problem, we apply a boost factor to  ★ .We will present heating series (i.e., series of simulations that are identical except for their heating boost factors) in Section 3.3 that can be used to extrapolate our results to nominal luminosity.Another benefit of calculating heating series is that deviations from established scaling laws at low luminosities can be used to identify numerical resolution issues (e.g., when the flow becomes too slow for our explicit gas dynamics code; Herwig et al. 2023).Note that the radiative conductivity, is always multiplied by the same boost factor to ensure energy conservation in the star.If more heat is generated in the central layers, then it must be transported more efficiently by radiation.We also perform simulations for three different grid resolutions to assess the numerical convergence of our calculations (Section 3.2).All simulations discussed in this work are listed in Table 1.The analysis of our simulations relies on three types of outputs (Andrassy et al. 2020;Stephens et al. 2021;Herwig et al. 2023).Every 1000-3000 time steps (depending on the grid resolution), a detailed output ("dump") is written to disk.Each of these dumps contains spherically averaged profiles, high-precision 3D briquette data (on a grid that is four times smaller in each direction than the simulation grid), and full-resolution byte-sized data cubes that we use to generate qualitative visualizations of the flow.
As shown in Table 1, each simulation is run for 25-60 h of star time.As we will see in the next section, the rms velocity in the convective core for the log / ★ = 5 simulations is of the order of  = 4 km s −1 .This implies a convective turnover timescale of ∼ 2 × 25 Mm/4 km s −1 ∼ 3 h.Our simulations therefore span ∼ 10 − 50 convective turnover timescales, depending on the heating rate3 and total simulation length.This is sufficient to robustly measure the mean properties of the flow.We show in Figure 5 the time evolution of the spherically averaged convective velocity at a radius located 2   (15 Mm) inside the convective boundary.We can see that for almost all simulations, a state that appears to be stationary on the convective timescale is reached after just a few turnover timescales (∼ 10 h).The initial transient before that time, when the convective flow is still building up, is discarded from our analysis in the rest of this work.
To quantify the absence of drifts in the velocity time series, we used linear regression to compute the slopes and estimate the percentage changes over the duration of each simulation.If we discard the first 10 h of each run from the linear regression, we find that only four simulations (W12, W17, W21 and W31) exhibit a change of more than 15% over their respective time durations.Most simulations are therefore unaffected by significant velocity drifts, and reliable time averages can be calculated.Note that two of the drifting simulations (W17 and W31) are low-heating runs at log / ★ = 4.5 that take more time to reach a stationary state due to the slower dynamics.As we will see, these low-heating simulations also suffer from numerical resolution issues (Section 3.3).The results of these runs have to be interpreted with caution, and we will ignore them in most of our analysis.and the much smaller velocities that characterize the stable envelope (which are barely visible in these renderings).In contrast, in the CPM simulations (Figure 6), the semiconvection zone imprints a region of moderate velocities between the high velocities of the convective core and low velocities of the stable envelope.It is evident from the ring-like structure of the flow in this region that the semiconvection zone is not dominated by turbulent convective motions.This will be investigated in more detail in Section 3.4.

Convergence with respect to the grid resolution
In Figure 8 we compare the spherically averaged velocity profiles of our log / ★ = 5 simulations performed using three different grid resolutions.For both setups, our results indicate that the properties of the flow in the convective core and in the boundary region are already well converged with respect to the grid resolution at 1152 3 since there is little difference between the 1152 3 and 1728 3 cases.The situation is less favourable in the envelope.With the CPM setup, there is no sign of convergence of the envelope velocities with respect to the simulation grid.The situation is better in the case of the PM setup, where the decreasing velocity difference between successive grid refinements suggests that the velocities are approaching convergence.The slower convergence of the velocities in the envelope compared to the core (or lack of convergence in the CPM case) is at least partially due to the smaller velocities in that region of the star.At the log / ★ = 5 heating rate shown in Figure 8, the Mach number in the envelope is only of the order of 10 −4 , a challenging regime for an explicit gas dynamics code.This also implies that a better convergence under grid refinement is expected for our simulations with higher heating luminosities, as the flow in the envelope will be more vigorous (Section 3.3).In any case, this issue is not a major concern in what follows, as the central goal of this work is to study the boundary region where for both setups the velocities converge much better under grid refinement.
We have seen that our PM simulations develop a fully convective core separated from the stable He envelope by a well-defined boundary (Figures 7-8).It is unsurprising to find a large fully mixed core with the PM prescription, but it is still interesting to note that we cannot identify any impact on the gas dynamics resulting from the presence of a region in the core where the Schwarzschild criterion is barely satisfied.Indeed, at  = 21.3Mm, ∇ rad − ∇ ad reaches its minimum in the convective core of just 0.002, and yet there is no obvious slowdown of the convective motions in that region (right panel of Figure 8).The global dipolar circulation pattern is completely oblivious to the presence of this minimum.At least at the heating rates of our simulations, this result confirms that it is consistent to assume full mixing in regions that have ∇ rad − ∇ ad > 0 by an arbitrarily small margin, as assumed in the maximum overshoot prescription and with MESA's PM scheme.

Heating series
As explained in Section 2.3, all our simulations are performed using heating luminosities that far exceed the nominal luminosity of the star.It is therefore important to understand how the properties of the flow scale with respect to the heating rate in order to properly extrapolate to the nominal case.To do so, we show in Figure 9 the rms velocity || in the convective core at  = 15 Mm as a function of / ★ .For the CPM setup, we precisely recover a  1/3 scaling law, as expected based on previous results of 3D hydrodynamics simulations of convection in stars (e.g., Porter & Woodward 2000;Müller et al. 2016;Jones et al. 2017;Baraffe et al. 2021;Herwig et al. 2023).Note also the excellent convergence with respect to the grid resolution: there is virtually no difference between results obtained using a 768 3 , 1152 3 or 1728 3 grid.However, the situation is different for the PM setup.First, departures from the  1/3 power law appear at log / ★ ≤ 5. Second, while there is excellent convergence with respect to the grid resolution at log / ★ > 5 and good convergence at log / ★ = 5 (consistent with Figure 8), there is a more significant difference between the 768 3 and 1152 3 runs at log / ★ = 4.5.For reference, the dotted line shows a  1/3 power law.All simulations listed in Table 1 were used to generate this figure.The rms velocities were averaged over the last 15 h of each simulation.
Given that the convective velocities are similar for both setups, it is surprising to see this departure from the expected scaling law in the PM setup while the expected behaviour is recovered down to at least log / ★ = 4.5 in the CPM case.Nevertheless, the fact that the  1/3 power law still applies for log / ★ ≥ 5 supports the use of this scaling relation to extrapolate to nominal luminosity using the log / ★ ≥ 5 simulations.
Figure 10 repeats the same exercise in the stable envelope at  = 33 Mm.Based on the results of Herwig et al. (2023), we expect a  2/3 scaling relation in this region where internal gravity wave (IGW) motions excited by the convective core characterize the flow (we will demonstrate the IGW-dominated nature of the flow in Section 3.4).This is indeed what we recover for both setups at high luminosities, and as expected we observe a departure from this scaling relation at lower luminosities.In addition, there is a systematic offset between the 768 3 and 1152 3 scaling laws.This resolution dependence is a reflection of the fact that the velocities in the envelope are not converged with respect to the grid resolution (Section 3.2).

Power spectra
We now go beyond spherical averages of the velocity field and investigate how much power the flow contains across different length scales.Figure 11 shows the power spectra of the tangential velocity magnitude |  | ([|  |], top row) and radial velocity   ([  ], bottom row), where [] is the power spectrum operator.Power spectra are displayed for different radii in the star and for different Cartesian grid resolutions.The power is binned as a function of the angular degree ℓ and is calculated using SHTools (Wieczorek & Meschede 2018).The maximum ℓ value shown in Figure 11 depends on the radius at which the power spectrum is calculated.A larger radius or grid resolution allows us to capture smaller scale features, since the angular resolution of the projection of the Cartesian grid on a sphere increases.Note that the spectra are calculated using the filtered briquette data, which is down-sampled by a factor of 4 in each direction.
Both for |  | and   , the power spectra in the convective core (blue lines,  = 15 Mm in Figure 11) is close to the expected Kolmogorov ℓ −5/3 cascade from large-scale modes (ℓ = 1 − 2) down to very small scales (ℓ ∼ 100 for high-resolution runs) where dissipation takes place.The convective power spectra are consistently extended to higher ℓ values when the grid resolution is increased.This reflects the fact that the turbulent cascade extends down to the smallest spatial scales resolved in the simulation.In the stable envelope (green lines,  = 35 Mm) and in the semiconvective region (orange lines,  = 25 Mm), the   power spectra display a very different shape compared to the convective core, with a shallow increasing slope up to ℓ ∼ 60 followed by a rapid drop at higher ℓ.This shape is reminiscent of that found in the stable layers of our recent mainsequence and red giant branch PPMstar simulations (Herwig et al. 2023;Blouin et al. 2023b) and hints at the IGW-dominated nature of the flow at those radii.Contrarily to the convective power spectra, the power at all ℓ ≲ 200 continuously increases upon grid refinement.This behaviour is consistent with the previously noted absence of convergence for the velocities in the stable envelope of the CPM runs (Section 3.2).We speculate that this increase in IGW power could be the result of the extension of the convective power spectra to larger ℓ values with increasing resolution.IGWs are possibly mainly excited by small-scale convective motions close to the core boundary, which would imply stronger IGW motions at high resolutions due to the enhanced power at high ℓ in the convective spectra.
For both setups, the power spectra at  = 35 Mm are approaching convergence for ℓ ≳ 100: the separation between the 1152 3 and 1728 3 spectra is smaller than that between the 768 3 and 1152 3 spectra.This is to be contrasted with the behaviour observed in the convective layers, where the power spectrum always extends to larger ℓ values upon grid refinement.This is most likely the result of radiative damping.Wave velocities are expected to be damped when radiative diffusion is taken into account, as the temperature of oscillating fluid parcels is equilibrated with their surroundings.This damping is stronger at large ℓ (Zahn et al. 1997), since small fluid parcels can lose their heat more quickly than large ones.This is consistent with the saturation of the IGW power spectrum observed in Figure 11.This result implies that radiative damping in our simulations only affects very small spatial scales (ℓ ≳ 100) that require large Cartesian grids to be properly resolved (above 1152 3 ).
We have so far been discussing various properties of the IGWs in the stable envelope, but without explicitly demonstrating that IGWs indeed dominate fluid motions in that region of the star.We remediate this issue in Figure 12.We show   power spectra taken at three different radii in run W20 (CPM setup, 1728 3 grid).Here, the power is binned both as a function of ℓ and the temporal frequency (the methodology used to calculate these spectra is detailed in Thompson et al. 2023).The power spectra calculated in the convective core shows power spread out over a wide range of spatial and temporal frequencies with no clear structure, as expected for turbulent motions.In contrast, the power spectra calculated in the stable envelope (bottom panel) displays distinct ridges constituted of discrete modes.These discrete modes quickly disappear for frequencies that surpass the local Brunt-Väisälä frequency (represented here by a thin dotted line).This is precisely what is expected for IGWs, which are damped when the Brunt-Väisälä frequency is exceeded.Figure 13 further demonstrates the IGW nature of these modes.Here, we zoom in on the low-ℓ portion of the bottom panel of Figure 12 and overlay eigenmodes predictions from the stellar oscillation code GYRE (Townsend & Teitler 2013;Townsend et al. 2018).This analysis is performed as in Thompson et al. (2023) and is based on spherically averaged radial profiles from run W20.Two GYRE calculations were performed: one using the structure corresponding to the beginning of the series of dumps used in calculating the power spectrum from the simula-tion data and one using the last dump of that series.The frequencies shown in Figure 13 are the averages of these two calculations.The excellent agreement between the GYRE predictions definitely confirms that these modes are IGWs, and it even allows the identification of individual radial orders.
Going back to Figure 12, we now focus on the particularly interesting power spectrum of the semiconvection zone (middle panel).We see the same discrete ridges as in the stable envelope, signalling the presence of IGWs and indicating that the semiconvective layers are dominated by IGW motions.The power distribution is admittedly more smeared out than in the stable envelope where the ridges are well separated, but this is presumably due to shorter mode lifetimes.The lack of power at high ℓ and low frequencies is further indication of the absence of convective motions in this region.We will investigate the potential implications of the IGW-dominated nature of the semiconvection zone on CBM in Section 4.2.

ON THE NATURE OF THE CONVECTIVE BOUNDARY
Having described the main properties of the flow, established the numerical convergence properties of our simulations and obtained scaling relations, we now shift our focus to the CHeB CBM problem introduced in Section 1.As a reminder, we have performed two series of simulations using two different initial setups with different CBM prescriptions because the accurate boundary mixing scheme is unknown.Can we use these simulations to infer the best CBM prescription to employ?In Section 4.1, we explore the most direct approach to tackle this question, namely investigating the evolution of the stratification in the boundary region.

Long-term evolution of the convective boundary
Figure 14 shows the time evolution of the  profile in the boundary region for our log / ★ = 6 simulations performed on a 1152 3 grid.For the CPM setup, we see that both the convective-semiconvective (left inset in the left panel) and semiconvective-radiative (right inset in the left panel) boundaries appear to be slowly migrating outward with time.The boundary is also changing for the PM setup, where the  profile is steepening with the outer portion of the boundary moving inward (right inset in the right panel).These changes are admittedly small compared to the 0.07 Mm grid cell size of these simulations, but they nevertheless constitute a robust result that we recover at other grid resolutions and heating rates.
Note that these migrations do not show sign of slowing down and approaching equilibrium after 40 h of simulation time.At the boundary, the radiative diffusivity is  rad ∼ 10 11 cm 2 s −1 in our log / ★ = 6 simulations, which implies that heat diffusion had time to operate over a length scale of only √︁  rad • 40 h ∼ 1 Mm.Hence, there is not enough time to reinstate thermal equilibrium following a local perturbation in the thermal structure of the stable layers due to a displacement of the convective boundary.It is therefore unsurprising that the boundary is still evolving at the end of our simulations, and reaching a stable boundary would demand much longer simulations (Anders et al. 2022a;Herwig et al. 2023;Mao et al. 2023).Nevertheless, it is tempting to extrapolate these trends to find the final equilibrated state towards which the simulations are evolving.For example, the apparent outward migration of the CPM setup may imply a larger mixed core than the one present in the initial MESA stratification.
However, such extrapolations are hazardous as it is not clear that these trends also apply to the real star.We have seen in Section 2.2 that the mapping between the MESA stratifications and the PPMstar base states is not perfect.The main problem is the omission of electron degeneracy pressure and ion-ion nonideal interactions in the equation of state currently employed by PPMstar, which leads to a 10% offset of the pressure profile.This inaccuracy affects the whole stratification, including the temperature and opacity profiles.This in turn perturbs the thermal balance of the star in a way that can be expected to shift its convective boundary.Another way to see this is to realize that a hypothetical MESA calculation with a different equation of state and/or opacity tables would yield a different convective boundary, since the Schwarzschild criterion would not be satisfied at the same radius.Hence, the small convective boundary reconfigurations observed in Figure 14 can at least partially be attributed to the response of the simulations to a change in the equation of state that has induced a state of thermal imbalance in the star.There is no obvious way to separate this behaviour from the (more interesting) 3D hydrodynamical response to the initial MESA base state, which is thermally balanced within the MLT framework.
All things considered, it is not currently possible to directly extrapolate the long-term evolution of our simulations to infer the "correct" stratification at the convective boundary of CHeB stars.Nevertheless, we will see in Sections 4.2 and 4.3 that we can gain further insights on the evolution of the convective boundary in our CPM simulations by isolating individual mixing processes through measurements of the diffusivity profile.

IGW mixing
We have described in Section 3.4 how the semiconvection zone of our CPM simulations is dominated by strong IGW motions.Our recent simulations of core convection in massive main-sequence stars have revealed species mixing in stable layers with strong IGW motions, presumably because of IGW-induced mixing (Herwig et al. 2023).A priori, the same phenomenon could conceivably occur inside the semiconvection zones of CHeB stars.If it is the case, then this would provide a mechanism to homogenize and destroy the semiconvective interface.
We have measured the diffusion coefficient in our CPM simulations using the method described in Jones et al. (2017) to invert the observed evolution of the FV profile, where FV is the fraction that represents the contribution of the envelope fluid to the two-fluid mixture. 5The resulting diffusivity profile for a log / ★ = 6 run is shown in black in Figure 15.Deep inside the convective core, the measured diffusivity is consistent with a simple MLT prescription (green dashed line). becomes much smaller than the MLT value closer to the convective-semiconvective boundary, a well-known limitation of this simple prescription (Eggleton 1972;Jones et al. 2017;Herwig et al. 2023;Blouin et al. 2023a,b).Once we reach the semiconvection zone,  drops precipitously and we cannot measure any mixing above  = 23.7 Mm.The mixing measured near the convectivesemiconvective boundary is due to the overshooting motions that we will explore further in Section 4.3.Beyond this boundary region, we measure no mixing in the semiconvection zone.This implies that despite their vigour, IGWs do not produce measurable mixing: we can establish an upper limit of  IGW ∼ 10 9 cm 2 s −1 in the semiconvection of our log / ★ = 6 simulations.Given the high heating rates of our simulations, we have every reason to believe that we are overestimating IGW mixing as we are overestimating IGW velocities by orders of magnitude.The upper limit on  IGW implied by our simulations is therefore much smaller at nominal luminosity.In the absence of measurable IGW mixing in our simulations, it is of course impossible to establish a scaling relation for  IGW that could be used to extrapolate to nominal luminosity.As an alternative, we can assume that the  IGW ∝  4/3 relation found in the IGW-dominated convective boundary layers of the H-core burning simulations without radiation diffusion of Herwig et al. (2023) still holds here.With this assumption, the upper limit on IGW mixing at nominal luminosity would be of only ∼ 10 cm 2 s −1 .To evaluate whether this level of mixing could have an impact on the star, we estimate how long it would take to homogenize the semiconvective layers that are located above the radius where we can still measure  in Figure 15.It would take more than (2 Mm) 2 / IGW ≃ 130 Myr to homogenize this 2 Mm-thick shell.This is a timescale comparable to the entire CHeB lifetime: IGW mixing in the semiconvection zone is negligible assuming the  IGW ∝  4/3 scaling relation is applicable to these simulations.
This result is consistent with the classical picture of semiconvection, where the diffusion coefficient in the semiconvective layers is taken to be (Langer et al. 1985) where  sc is a free parameter that controls the mixing timescale,   is the heat capacity at constant pressure, and ∇ L is the Ledoux gradient.
Once an adiabatic stratification is reached (∇ − ∇ ad = 0, as is the case in the semiconvective region of our CPM setup), no further mixing takes place and the composition profile remains constant if secular changes to the structure are neglected.In this picture, the amplitudes of IGWs saturate before they become strong enough to overturn (e.g., Merryfield 1995) and trigger nonlinear turbulent mixing.is convenient as our PPMstar calculations follow the evolution of two fluids (a C/O-rich core fluid and an almost pure-He envelope fluid), but it will inevitably underestimate the true rate at which the convective core is eroding the semiconvective region in our simulations.Indeed, the semiconvective region is closer in composition to the convective core than to the stable envelope (FV < 0.5), meaning that most of the entrained mass is not captured by our definition.As we will see, this issue has no impact on the main conclusion of the mass entrainment rate analysis.

Convective
The top panel of Figure 19 demonstrates how the mass entrainment rate is measured for a given simulation, and the bottom panel shows the mass entrainment scaling relation thus found.To stay clear from the initial transient, we have evaluated the mass entrainment rates using only  > 20 h, and for more robust measurements we only employed simulations longer than 40 h.As in previous works (Jones et al. 2017;Andrassy et al. 2020;Herwig et al. 2023;Mao et al. 2023), we find that the entrainment rate scales linearly with heating (if we exclude the lowest-heating run, consistent with our discussion in Section 3.3).Extrapolating this relation to nominal luminosity  yields an entrainment rate of 4.5 × 10 −8  ⊙ yr −1 .The 4 Mm-thick semiconvection zone contains 0.08  ⊙ : assuming entrainment at this constant rate, it would take less than 2 Myr to completely erase it.This is much shorter than any relevant evolutionary timescale.In particular, it is significantly faster than the rate at which the convective core grows in CHeB stars (a growth rate of ≃ 3 × 10 −9  ⊙ yr −1 can be inferred from Figure 2).According to this analysis, we should therefore expect the semiconvection zone to be completely erased, implying that CHeB stars cannot have a semiconvection zone.In other words, a semiconvection zone cannot form in the first place, and the initial CPM setup is incorrect.Note that our restrictive definition of the entrained mass, which necessarily underestimates the total entrained mass in our simulations, can only strengthen this conclusion.
However, there is a major caveat with this analysis.As explained in Section 4.1, the convective boundary may also be migrating for The log / ★ = 4.5 run was omitted from the fit; its departure from the linear relation is likely a manifestation of the numerical issues that appear at low-Mach numbers with our explicit gas dynamics code (Section 3.2).
other spurious reasons related to the current PPMstar equation of state.When we measure an entrainment rate, we capture the sum of all the processes at play.We may therefore be overestimating the rate at which the semiconvection zone is eroded due to real mixing processes as opposed to the boundary migration due to the equation of state).That being said, independent results also indicate that the semiconvection zone should be quickly erased.Anders et al. (2022b) have recently performed hydrodynamics simulations of an idealized plane-parallel setup containing a semiconvection zone between a convective and a stable region.They describe a similar behaviour to that detailed above, with overshooting convective motions gradually entraining low- material into the convective region and homogenizing the semiconvective layers.Ultimately, after thousands of convective turnover timescales, they report the disappearance of the semiconvection zone.

CONCLUSION
We have presented the first full-sphere 3D hydrodynamics simulations of the interior of a CHeB star.For dozens of convective turnover timescales, we followed the hydrodynamical response of our highresolution simulations to two different initial stratifications (one with a semiconvective region and one without) for a 3  ⊙ CHeB star.
We have recovered many of the key findings of our recent simulations of core convection in massive main-sequence stars, including the presence of a large dipole circulation pattern in the convective core, the excitation of a rich spectrum of IGWs in the stable envelope consistent with the eigenfrequencies predicted by GYRE for the same stratification, and the  1/3 and  2/3 scaling of the convective and IGW velocities, respectively.We found that the extended core of the PM prescription remains fully convective even if ∇ rad − ∇ ad becomes very close to 0 (while remaining positive).We have described how an hypothetical semiconvection zone would be dominated by IGW motions.Despite the high amplitudes of these waves, we see no evidence for IGW mixing, in contrast with our recent finding for the core boundary region of massive main-sequence stars (Herwig et al. 2023).
For simulations initialized with a semiconvective interface, we have also observed the incursion of convective motions inside the semiconvective zone, a phenomenon that gradually erases this region.While the efficiency of this mixing process at nominal luminosity remains unclear, it could be sufficient to homogenize the semiconvection zone of a CHeB star much more rapidly than any relevant evolutionary timescale.This would imply, as recently suggested by Anders et al. (2022b), that CHeB stars cannot harbour a semiconvective interface between their C/O-rich cores and He envelopes.
Future research should investigate other CBM prescriptions beyond the MESA CPM and PM schemes to model the CHeB phase.A growing body of work (Anders et al. 2022a;Andrassy et al. 2023;Baraffe et al. 2023;Blouin et al. 2023b;Mao et al. 2023) is suggesting that the stratifications implied by these two CBM schemes is not appropriate.In particular, multi-dimensional hydrodynamics simulations consistently point to the formation of a convective penetration zone where the temperature gradient smoothly transitions from ∇ ad to ∇ rad over a fraction of a pressure scale height (∼ 0.1 − 1   ).This is in clear tension with the discontinous transitions implied by MESA's CPM and PM schemes.

Figure 1 .
Figure1.Qualitative representation of the radiative gradient profiles described in Section 1.The actual temperature gradient ∇ T is equal to ∇ ad in the regions where ∇ rad /∇ ad > 1 and equal to ∇ rad otherwise.

Figure 2 .
Figure2.Kippenhahn diagram of the CHeB phase as calculated in our PM MESA run.The dark blue regions show where nuclear burning takes place (He burning at the center and H burning further out), and the grey region shows the extent of the convective core.The dashed vertical line marks the model used for our 3D PPMstar calculations; its vertical extent corresponds to the total mass included in the simulations.

Figure 3 .
Figure 3. Radial profiles for our CPM setup.The MESA profile (blue solid line) and the base state of our PPMstar simulations (dashed orange line) are shown for each quantity.In the last panel, the dotted green line shows the opacity profile obtained using the MESA stratification and the PPMstar opacity module.The semiconvection zone is the region between  ≃ 22 and 26 Mm characterized by a shallow  gradient.

Figure 4 .
Figure 4. Same as Figure but for our PM setup.

Figure 5 .Figure 7 .
Figure 5.Time evolution of the spherically averaged rms convective velocity | | at a radius located 2   (15 Mm) inside the convective boundary.The left panel shows the runs using the CPM setup and the right panel the PM setup.Most simulations reach a stationary state after a few convective turnover timescales.

Figure 8 .Figure 9 .
Figure 8. Spherically averaged tangential (orange) and radial (blue) velocity profiles for the CPM (left) and PM (right) setups for different grid resolutions.Simulations W12, W13, W20, W21, W22 and W23 were used to generate this figure (log / ★ = 5 heating).The velocity profiles were obtained by averaging over the last 400 dumps of each simulation.

Figure 10 .Figure 11 .
Figure 10.Same as Figure 9, but this time looking at the velocities at  = 33 Mm in the stable envelope.

Figure 12 .
Figure 12.Left: Centre-plane slice rendering of the vorticity magnitude | ∇ ×  | of run W20 (CPM setup, 1728 3 grid) at dump 840 ( = 20.3h).Dark blue, turquoise, yellow, red, and dark red represent a sequence of increasing vorticity magnitudes.Unlike in Figure 6 and 7, the full simulation domain is shown.The white circles indicate the radii at which power spectra were calculated.Right: Power spectra of the radial velocity components at  = 20, 25, 30 Mm.The spectra were calculated over dumps 480 to 880 ( = 11.6 − 21.3h).The local value of the Brunt-Väisälä frequency is shown by the dotted horizontal line.

Figure 13 .
Figure 13.Power spectra of the radial velocity component at  = 30 Mm for run W20 calculated over dumps 480 to 880 ( = 11.6 − 21.3h).The local value of the Brunt-Väisälä frequency is shown by the dotted horizontal line.The overlaid solid white lines show the eigenmodes predicted by GYRE for the stellar structure corresponding to that of run W20.The GYRE analysis was carried out using the spherically averaged structures of dumps 480 and 880: the frequencies shown here are averages of these two calculations.The labels next to the white lines indicate the radial orders  of the modes (their negative signs indicate that they are  modes).

Figure 14 .Figure 15 .
Figure 14.Evolution of the mean molecular weight profile in the boundary region for a CPM (left panel, run W28) and a PM (right panel, run W29) simulation.These runs employed a 1152 3 grid and a log / ★ = 6 heating rate.The legends above the panels indicate the simulation times at which the profiles are shown.

Figure 16 .
Figure16.Center-plane slice rendering of the vorticity magnitude in run W32 (CPM setup, 768 3 grid, log / ★ = 7) at dump 1650 ( = 39.9 h).In the southwest quadrant, note the intrusion of the dipole circulation pattern of the convective core into the semiconvective layers.This is to be contrasted with the log / ★ = 5 simulation rendered in Figures6 and 12, where the turbulent motions are confined to a spherical region circumscribed by a clear boundary at  ≃ 22 Mm.White arrows indicate physical distances from the centre of the star.

Figure 17 .
Figure17.FV profile at  = 20 h for our CPM runs performed on a 768 3 Cartesian grid and with different heating rates (see legend).The black solid line corresponds to the initial setup.The semiconvection zone in the initial setup is shaded in grey.Note how the FV profile in the semiconvective layers is homogenized in the highest-heating run (W32).

Figure 18 .
Figure 18.Measured diffusion coefficient in our CPM simulations (see legend for grid resolution and heating rate).For each run,  was calculated by inverting the evolution of the FV profile between  = 20.0 − 22.5 h and  = 22.5 − 25.0 h.

Figure 19 .
Figure19.Top: Time evolution of the entrained mass in run W10 (blue line, log / ★ = 6 and 768 3 grid) and best-fit linear entrainment rate (black dotted line).Bottom: Heating series of the measured mass entrainment rate (blue and orange symbols) and best-fit linear scaling law (black dotted line).The log / ★ = 4.5 run was omitted from the fit; its departure from the linear relation is likely a manifestation of the numerical issues that appear at low-Mach numbers with our explicit gas dynamics code (Section 3.2).

Table 1 .
Summary of simulations used in this work.

Table A1 .
Opacity model fit parameters