Issue 
A&A
Volume 654, October 2021



Article Number  A126  
Number of page(s)  13  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202140441  
Published online  22 October 2021 
Twodimensional simulations of solarlike models with artificially enhanced luminosity
I. Impact on convective penetration
^{1}
University of Exeter, Physics and Astronomy, EX4 4QL Exeter, UK
email: i.baraffe@ex.ac.uk
^{2}
École Normale Supérieure, Lyon, CRAL (UMR CNRS 5574), Université de Lyon, Lyon, France
^{3}
Department of Physics and Astronomy, Georgia State University, Atlanta, GA 30303, USA
^{4}
Centre for Fusion, Space and Astrophysics, Department of Physics, University of Warwick, Coventry CV4 7AL, UK
Received:
28
January
2021
Accepted:
20
August
2021
We performed twodimensional, fully compressible, timeimplicit simulations of convection in a solarlike model with the MUSIC code. Our main motivation is to explore the impact of a common tactic adopted in numerical simulations of convection that use realistic stellar conditions. This tactic is to artificially increase the luminosity and to modify the thermal diffusivity of the reference stellar model. This work focuses on the impact of these modifications on convective penetration (or overshooting) at the base of the convective envelope of a solarlike model. We explore a range of enhancement factors for the energy input (or stellar luminosity) and confirm the increase in the characteristic overshooting depth with the increase in the energy input, as suggested by analytical models and by previous numerical simulations. We performed highorder moments analysis of the temperature fluctuations for moderate enhancement factors and find similar flow structure in the convective envelope and the penetration region, independently of the enhancement factor. As a major finding, our results highlight the importance of the impact of penetrative downflows on the thermal background below the convective boundary. This is a result of compression and shear which induce local heating and thermal mixing. The artificial increase in the energy flux intensifies the heating process by increasing the velocities in the convective zone and at the convective boundary, revealing a subtle connection between the local heating of the thermal background and the plume dynamics. This heating also increases the efficiency of heat transport by radiation which may counterbalance further heating and helps to establish a steady state. We suggest that the modification of the thermal background by penetrative plumes impacts the width of the overshooting layer. Additionally, our results suggest that an artificial modification of the radiative diffusivity in the overshooting layer, rather than only accelerating the thermal relaxation, could also alter the dynamics of the penetrating plumes and thus the width of the overshooting layer. Results from simulations with an artificial modification of the energy flux and of the thermal diffusivity should thus be regarded with caution if used to determine an overshooting distance.
Key words: convection / hydrodynamics / instabilities / stars: evolution
© ESO 2021
1. Introduction
One of the major uncertainties in stellar physics is the treatment of mixing taking place at convective boundaries. Convective motions do not abruptly stop at the classical Schwarzschild boundary, but extend beyond it. The complex dynamics resulting from convective penetration in stable layers is a major process in stellar interiors that drives the transport of chemical species and heat, strongly affecting the structure and the evolution of many types of stars. This process is usually called overshooting, convective penetration or convective boundary mixing. The same complex dynamics can also drive the transport of angular momentum, affecting the rotational evolution of stars, the generation of a magnetic field in their interior and their magnetic activity. Overshooting is one of the oldest unsolved problems of stellar structure and evolution theory (e.g., Shaviv & Salpeter 1973) and affects all stars that develop a convective envelope or core, meaning all stars with a mass above ∼0.4 M_{⊙}.
Analytical and semianalytical approaches have been developed to describe this process and estimate the width of the overshooting layer (e.g., Schmitt et al. 1984; Zahn 1991; Rempel 2004). With the improvement of computational methods and resources, an increasing number of studies have been devoted to numerical simulations of convection and overshooting using realistic stellar conditions (geometry, luminosity, thermal diffusivity, equation of state, opacities, etc.). A commonly used tactic to increase the efficiency and improve the stability of these simulations is to artificially increase the luminosity (or nuclear energy for convective cores or burning shells) and to modify the thermal diffusivity of the reference stellar model. This tactic is common and has been used, for example, in Rogers et al. (2006, 2013), Meakin & Arnett (2007), Tian et al. (2009), Brun et al. (2011, 2017), Hotta (2017), Cristini et al. (2017), Edelmann et al. (2019), Horst et al. (2020). This approach is used to increase the Mach number of the convective flow, reducing the disparity between advective and acoustic timescales, improving the efficiency of timeexplicit codes limited by the Courant–Friedrich–Levy constraint. It is also used to provide numerical stability or to accelerate the thermal relaxation. But no examination of its potentially farreaching impact has been conducted. Rempel (2004) pointed out that enhanced energy flux in numerical simulations could lead to unrealistically vigourous convection, which could impact the properties of the overshooting layer and could explain some of the discrepancies between analytical models and numerical simulations. Numerical simulations also suggest an increase in the overshooting depth with increasing flux (Hotta 2017; Käpylä 2019). Determining scaling laws of the overshooting depth as a function of the energy input could thus allow for an extrapolation of the results to more realistic stellar conditions and help to estimate the overshooting depth in real stars, as suggested by Hotta (2017), for example, for the Sun. But Käpylä (2019) also shows that an artificial modification of the heat conductivity in the radiative and overshooting regions could impact the overshooting process. In some computational studies, both the luminosity and the thermal diffusivity are enhanced by the same factor to ensure that the thermal structure is unchanged compared to the reference stellar structure and with the expectation that the larger thermal diffusivity counterbalances the larger energy flux. This procedure has been proposed as a way to provide a good representation of the true dynamics of the system (e.g., Rogers et al. 2006, 2013; Tian et al. 2009). But this expectation has never been demonstrated. Another expectation concerns internal waves, excited by convective motions and by flows penetrating the convective boundary. Simulations with artificially modified luminosity and thermal diffusivity are also used to perform the analysis of internal waves, either for convective envelopes (e.g., Rogers et al. 2006; Brun et al. 2011; Alvan et al. 2014) or for convective cores (e.g., Rogers & McElwaine 2017; Edelmann et al. 2019; Horst et al. 2020). None of these works have examined whether the wave spectrum of a realistic star is accurately predicted by such simulations.
The modification of the energy flux and thermal diffusivity commonly performed in stellar hydrodynamics simulations of convection and overshooting thus raises several questions. Firstly, to which extent can scaling laws of the overshooting depth with luminosity be extrapolated down to realistic stellar luminosities? Secondly, does the enhancement of the luminosity and thermal diffusivity only impact the timescale of the simulation, that is accelerates it, or does it have other effects on the dynamics of the region of interest? Thirdly, do such boosted models predict a realistic spectrum of internal waves generated at the convective boundary and propagating in the stable region?
We attempt to address these questions by performing a suite of twodimensional simulations for a solarlike star and analysing the impacts of artificially enhanced luminosity and thermal diffusivity. These experiments are restricted for now to twodimension models, which allow for longer simulations in order to calculate accurate statistics, and the examination of multiple simulations to explore a wider parameter range. We will discuss in Sect. 6 the impact of the dimensionality on our main results. In this work we examine the properties of convection and convective penetration at the base of the convective envelope. Two followup studies are underway. The first one is devoted to the properties of internal waves generated at the convective boundary (Le Saux et al., in prep.). The second one is devoted to the analysis of the overshooting depth based on a combination of the simulations presented in this work and Lagrangian tracer particles (Guillet et al., in prep.).
2. Numerical simulations
In this work we use the fully compressible timeimplicit code MUSIC. A full description of MUSIC and of the timeimplicit integration can be found in Viallet et al. (2011, 2016), Goffrey et al. (2017). Here, we provide a brief description of its main characteristics. MUSIC solves the inviscid Euler equations in the presence of external gravity and thermal diffusion:
where ρ is the density, e the specific internal energy, v the velocity, p the gas pressure, T the temperature, g the gravitational acceleration, and χ the thermal conductivity. The symbol ⊗ is the outer product. All hydrodynamical simulations presented in this work are performed assuming spherically symmetric gravitational acceleration that does not evolve with time, meaning that g is calculated at time t = 0 using the initial density profile and remains constant with time. For the stellar model that we consider in this work, radiative transfer is the major heat transport that contributes to the thermal conductivity, which is given for photons by,
where κ is the Rosseland mean opacity, and σ the Stefan–Boltzmann constant. In the following, thermal diffusivity and radiative diffusivity are used interchangeably. Realistic stellar opacities and equation of states appropriate for the description of stellar interiors are implemented in MUSIC. Opacities are interpolated from the OPAL tables (Iglesias & Rogers 1996) for solar metallicity and the equation of state is based on the OPAL EOS tables of Rogers & Nayfonov (2002), which are appropriate for the description of solarlike interior structures.
2.1. Initial stellar models
The multiD simulations require as initial input a radial profile of density and internal energy; for this work, those initial profiles are provided by the onedimensional Lyon stellar evolution code (Baraffe & El Eid 1991; Baraffe et al. 1998), using the same opacities and equation of state as implemented in MUSIC. We have chosen stellar interior structures close to the Sun’s structure, that is a solar mass star on the Main Sequence with a convective envelope covering ∼30% of the stellar radius. Our motivation is to use initial structures as close as possible to realistic stellar interior structures, as done in our previous studies (see Pratt et al. 2016, 2017). But some caution is needed to construct an initial model using a stellar evolution code if the goal is to test the impact of artificially enhanced luminosity and radiative diffusivity. In a typical solar model calculated with the Mixing Length Theory and a mixing length l_{mix} = 1.9H_{P}, the superadiabaticity (see definition below) in the bulk of the convective zone is very small, typically ≪10^{−4}, but in the outer convective zones the superadiabaticity is high, typically > 10^{−2}. We recall that the superadiabaticity is defined as (∇ − ∇_{ad}) with the temperature gradient and the adiabatic gradient. The Schwarzschild boundary is defined as the transition layer between convective instability (∇ > ∇_{ad}) and stability (∇ < ∇_{ad}). The outer structure is thus very sensitive to any change of the opacity (and thus of the thermal diffusivity) and of the luminosity, as such changes will modify the superadiabaticity and thus the temperature stratification. Therefore, in order to avoid a readjustment of the model structure when starting a hydrodynamical simulation using MUSIC with enhanced luminosity and thermal diffusivity, the profile of the stellar structure model must be adiabatic. We have thus constructed an artificial solarlike model with our stellar evolution code enforcing a very small superadiabaticity (< 10^{−8}) throughout the convective zone. In this case, an increase in the luminosity and of the radiative diffusivity (or a decrease of the opacity by the same factor) has no impact on the model structure (in terms of density and temperature radial profiles). This yields a reference initial model slightly more compact and hotter than a model for the current Sun (that is a 1 M_{⊙} model with the luminosity and radius of the current Sun) calculated with an initial helium abundance of 0.28, metallicity Z = 0.02 and l_{mix} = 1.9H_{P}, as shown in Fig. 1.
Fig. 1.
Internal structure of the reference model (solid lines) used as initial model for the 2D simulations compared to a model for the current Sun (dashed lines). The figure shows the density ρ (in g cm^{−3}), temperature T (in K) and radiative diffusivity κ_{rad} = χ/(ρc_{P}) (in cm^{2} s^{−1}), with c_{P} the specific heat at constant pressure. The Schwarzschild boundary is indicated by the vertical solid line for the reference model and the vertical dashed line for the Sunlike model. 
The properties of the reference model are displayed in Table 1. We stress that starting numerical hydrodynamic simulations from realistic stellar structures with artificially enhanced luminosities and thermal diffusivities requires the superadiabaticity in the convective zone of the initial model to be extremely small. If not, hydrodynamical models will relax towards a structure that departs from the initial 1D structure, and the larger the luminosity enhancement factor, the higher the departure because of the shorter thermal relaxation timescale. Consequently, the conclusions one may derive by comparing simulations with different luminosity enhancement factors will be affected by the additional impact of having different structures. This would be equivalent to comparing different stellar models with different luminosities.
Properties of the initial reference model: mass, luminosity, radius, depth of the convective envelope and pressure scale height at the Schwarzschild boundary.
2.2. Sphericalshell geometry and boundary conditions
We perform twodimensional simulations in a spherical shell using spherical coordinates, namely r the radius and θ the polar angle, and assuming azimuthal symmetry in the ϕdirection. The inner radius r_{in} is defined at 0.4 R_{star} and the outer radius at r_{out} = 0.9 R_{star}. The angular extent ranges from θ = 0° to θ = 180°, including the full hemisphere. We use a uniform grid resolution of r × θ = 512 × 512 cells. This provides a good resolution of the pressure scale height at the Schwarzschild boundary H_{P, CB}/Δr ∼ 92, with Δr = 5.5 10^{7} cm (550 km) the radial grid spacing. In the θ direction, the typical size of a cell is 2300 km. The choice of the resolution in the θ direction is set by the requirement to preserve a good aspect ratio of the grid cells on the whole domain on a spherical grid. We note that increasing the resolution by a factor two (1024 × 1024) does not change our results and conclusions (see discussion in Sect. 6). In terms of boundary conditions, we impose a constant radial derivative on the density on the inner and outer radial boundaries, as discussed in Pratt et al. (2016). For the reference model ref, the energy flux at the inner and outer radial boundaries are set to the value of the energy flux at that radius in the onedimensional stellar evolution model. For the artificially boosted simulations, the energy flux, and equivalently the luminosity, at the boundaries is multiplied by an enhancement factor, and the Rosseland mean opacities κ in MUSIC are decreased by the same factor. In this work we analyse the impact of enhancing the luminosity and thermal diffusivity by factors 10, 10^{2} and 10^{4}. Larger enhancement factors (up to 10^{6}) can be found in previous works (e.g., Rogers et al. 2006; Hotta 2017), but our selected range and values are appropriate for our general purpose, namely exploring the impact of artificially increasing the luminosity on the properties of overshooting and of internal waves. In velocity, we impose on the radial boundaries nonpenetrative condition for the radial velocity and stressfree boundary condition for the angular velocity. At the boundaries in θ, because of the extension of the angular domain to the “poles”, we use reflective boundary conditions for the density and energy, meaning that they are mirrored at the boundary. We adopt a stress free boundary condition for the radial velocity and a reflecting boundary condition for the angular velocity, to ensure it is equal to zero at the boundary.
3. Results: Average dynamics and fluxes
3.1. General properties
Our four simulations are summarized in Table 2. Figure 2 shows the time evolution of the total kinetic energy for the four models over the entire spherical volume with ρ the density, v the total velocity and dV the volume element. After a relaxation phase characterised by the propagation of strong acoustic waves and the onset of convection, E_{kin} reaches a plateau which characterises the steady state for the convection. The time t_{steady} to reach this state is indicated in Table 2 for each model. Our numerical simulations are not thermally relaxed. Achieving thermal relaxation is a wellknown challenge for global hydrodynamical simulations of convection based on realistic stellar interior structures (e.g., Meakin & Arnett 2007; Horst et al. 2020; Higl et al. 2021). The global thermal relaxation or KelvinHelmholtz timescale of the initial stellar structure used for our simulations is given by τ_{th} ∼ GM^{2}/(RL) ∼ 4 × 10^{6} yr. Our total simulation times (see Table 2) remain much shorter than τ_{th}, even for the most boosted model analysed here and for which τ_{th} ∼ 4 × 10^{2} yr. By specifically choosing as reference model a stellar model with no enhancement factor for the luminosity, this is unavoidable. As a consequence all these simulations are expected to maintain a secular drift. For the ref, boost1d1 and boost1d2 models, the drift is so slow that calculating statistical data and thermal properties during this very slowly changing transitional state is sensible. However Fig. 2 shows that the most boosted model is so far out of balance that it is continuously evolving. As will be discussed in Sects. 4 and 5, time averages for this model have limited meaning.
Fig. 2.
Evolution of the total kinetic energy (in erg) as a function of time (in s) of the simulations described in Table 2. The dotted lines correspond to the value of the total kinetic energy at the beginning of the steady state for convection. 
Summary of the 2D simulations.
In the following, time averages are denoted by ⟨⟩_{t} and calculated between t_{steady} and t_{sim} (see values in Table 2), i.e., for any quantity X we define:
We estimate a global convective turnover time τ_{conv} based on the rms velocity v_{rms}(r, t) at radius r and time t, which characterises a bulk convective velocity:
where the rms velocity is given by,
with , v_{r} and v_{θ} being the radial and angular velocities, respectively. ⟨⟩_{θ} denotes a volumeweighted average in the angular (θ) direction and is defined for any quantity X as,
The values of τ_{conv} are provided in Table 2 and the corresponding rms velocity profiles are shown in the top panel of Fig. 3, for reference. The turnover time follows a scaling τ_{conv} ∝ L^{−1/3}, as expected from the scaling of v_{rms} with luminosity v_{rms} ∝ L^{1/3} based on theoretical arguments from the mixinglength theory (Biermann 1932) and confirmed by many hydrodynamical simulations (e.g., Porter & Woodward 2000; Viallet et al. 2013; Jones et al. 2017; Edelmann et al. 2019; Andrassy et al. 2020). Our simulations reproduce this scaling as illustrated in the bottom panel of Fig. 3 which displays the rms velocity and rms radial velocity for the four models. The rms velocities observed in the stably stratified region are due to the propagation of internal waves excited by the convective motions and penetrative plumes at the convective boundary. Interestingly, the amplitude of these waves also follow a scaling with the luminosity. We find that the rms velocities follow very closely a scaling v_{rms} ∝ L^{0.61} for all models and for the radial velocity, the scaling is close to v_{r, rms} ∝ L^{0.86}. Analysis of the internal waves is beyond the scope of this work and will be presented in a followup study. Our focus in the present study is on the process of overshooting. Figure 3 clearly shows deeper penetration of the convective motions beyond the convective boundary with increasing energy input. This is predicted by analytic models of overshooting (e.g., Zahn 1991; Rempel 2004) and found by numerical simulations in the literature (e.g., Hotta 2017; Käpylä 2019). In the following section (see Sect. 4), we will quantitatively estimate overshooting depths based on our approach developed in Pratt et al. (2017). This approach relies on high order statistics of energy fluxes, such as the kinetic energy and heat fluxes, in contrast to the previous estimates which use only average quantities (e.g., Hurlburt et al. 1994; Brummell et al. 2002; Rogers et al. 2006; Käpylä et al. 2017).
Fig. 3.
Top panel: radial profile of the time averaged rms velocity (solid lines) and rms radial velocity (dashed lines). Velocities are in cm/s. The curves from bottom to top correspond to the four models: ref (black), boost1d1 (blue), boost1d2 (magenta) and boost1d4 (red), respectively. Bottom panel: velocities scaled by the enhancement factor (L/L_{star})^{1/3}. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 
An inspection of averaged fluxes is still interesting to further illustrate the dependence of penetration on luminosity. At each radius r, we define a mean (time average and volumeweighted average in the θ direction) vertical kinetic energy flux F_{k}, enthalpy flux F_{h} or equivalently convective heat flux F_{δT} by:
where H = e + P/ρ is the specific enthalpy and the second term in the r.h.s. of Eq. (10) subtracts any nonzero mean vertical mass flux (Freytag et al. 1996). The temperature fluctuation δT for each gridcell and time t is defined as,
The luminosities L = 4πr^{2}F corresponding to these fluxes are displayed in Fig. 4 for the four models. Note that in the literature, the enthalpy flux is calculated using either Eqs. (10) or (11). Both expressions give similar results as shown by the excellent agreement between these two quantities in Fig. 4. We show in addition the mean radiative flux F_{rad} given by:
Fig. 4.
Radial profile of luminosities for the four models as indicated in each panel, corresponding to various fluxes: vertical enthalpy flux (L_{h}, dashdot, blue), radiative flux (L_{rad}, long dashdot, black), vertical kinetic energy flux (L_{k}, dot, red), convective heat flux (L_{δT}, dash, magenta). The total luminosity L_{tot} = L_{h} + L_{rad} + L_{k} (solid, green) is also displayed. Luminosities are divided by the luminosity of the corresponding model. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 
with . The layer below the convective boundary where convective penetration proceeds can be characterised by the negative peak of the enthalpy flux (Hurlburt et al. 1986; Muthsam et al. 1995; Brun et al. 2011; Pratt et al. 2017; Korre et al. 2019; Käpylä 2019). Convective downflows transport low entropy (cool material) from the stellar surface down to the bottom of the convective zone. Inspection of the temperature fluctuations indeed indicates that downflows are characterised by negative δT and upflows by positive δT in the convective zone (see Fig. 5 and Sect. 5). But when the downflows cross the convective boundary and penetrate the stably stratified medium, they are adiabatically compressed and therefore get hotter (positive δT) and less dense (negative density fluctuation δρ) than the subadiabatically stratified environment. Upward flows have the reverse pattern so both flow types contribute to a negative enthalpy flux in the layer where the bulk of the convective plumes penetrate.
Fig. 5.
Visualisation of the radial velocity and the temperature fluctuations close to the convective boundary (horizontal black line) for the ref model at a given time. The xaxis represents the colatitude (in terms of cos θ). 
The properties of the flow and of the temperature fluctuations close to the convective boundary are shown in Fig. 5 at a given time for the ref model, illustrating the patterns above described which are common to all models. The figure shows that the present simulations can capture the front of plumes penetrating across the convective boundary and the waves excited and propagating below it, even if our resolution does not allow to fully resolve their detailed structure. The behaviour of the enthalpy flux in the penetration layer above described is a known feature and has already been discussed in detail in the literature (e.g., Hurlburt et al. 1986; Muthsam et al. 1995). However, less attention has been given to the impact on the thermal profile and more specifically to the feedback between the radiative flux and the temperature profile in this layer. This is the focus of the next section.
3.2. Modification of the thermal background
The position of the negative peak of the enthalpy flux moves deeper inward as the luminosity of the model increases, confirming the larger extension of the overshooting depth with input flux already illustrated with the analysis of the rms velocities (see Fig. 3). We can clearly see for the boosted models that the position of the negative peak of the enthalpy flux corresponds to the position of a peak of the radiative flux. A peak of the radiative flux also exists for the ref model, although much less pronounced. This feature has also been observed in previous numerical simulations (Rogers et al. 2006; Brun et al. 2011, 2017; Käpylä 2019; Cai 2020). In Brun et al. (2011, 2017), the authors explain that ideally the system would compensate for the negative heat flux, considered as a flux deficit, and adjust to a new equilibrium by modifying the background thermal stratification in a thermal timescale. In order to achieve satisfactory flux balance, Brun et al. (2011, 2017) arbitrarily increase the radiative diffusivity near the base of the convective zone, in order to accelerate the thermal relaxation process. This thermal adjustment yields a peak of the radiative flux in order to compensate for the flux deficit near the base of the convective zone.
In the simulations of Rogers et al. (2006) devoted to a solar model, similar results are found with the temperature gradient in the penetration region becoming slightly steeper, a feature that is interpreted in the same way as Brun et al. (2011, 2017), namely that the upward diffusive heat flux increases to compensate for the negative convective heat flux in that region. Rogers et al. (2006) also find that the region heats up because convective motions are continually pumping heat into this region and the timescale for this heat transfer is much shorter than the diffusive timescale. This local heating causes a steeper temperature profile and a less subadiabatic temperature gradient, i.e., ∇−∇_{ad} decreases, immediately below the convective boundary. Immediately below this region, the temperature profile is slightly flatter (hotter) and the temperature gradient becomes more subadiabatic, i.e., ∇−∇_{ad} increases (see Fig. 3 in Rogers et al. 2006). Note that Korre et al. (2019) also find that the fluid motions affect the thermal stratification at the base of the convective zone in the penetration region. But since Korre et al. (2019) use the Boussinesq approximation, it is not clear how the modification of the thermal background is accounted for in such framework. Finally, in 2D simulations for convective core overshooting, Higl et al. (2021) note also a modification of the temperature stratification in the penetration region compared to the initial 1D model. Despite the fact that their simulations are not thermally relaxed, they suggest that penetration will very likely impact the final thermal background in the overshooting region.
We find similar features (see Fig. 6): a slight heating of the layers below the position of the negative peak of the heat (enthalpy) flux. The peak and the nonmonotonic profile of the radiative flux is a consequence of the slight change of the averaged thermal profile, which impacts the radiative conductivity χ and the radial temperature gradient .
Fig. 6.
Radial profile of time and space averages of various quantities close to the convective boundary. From top to bottom: the temperature T (in units of 6 × 10^{6} K); the radiative conductivity χ divided by the luminosity enhancement factor of the model and by the constant 4 × 10^{15}; the radiative flux divided by the luminosity enhancement factor of the model and by the constant 2 × 10^{12}; the subadiabaticity ∇ − ∇_{ad}; the radial temperature gradient ∂T/∂r multiplied by the constant 2.8 × 10^{3}. The coloured curves correspond to the four models: ref (black), boost1d1 (blue), boost1d2 (magenta) and boost1d4 (red). The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 
In order to understand the origin of the local heating, we examine the rateofstrain tensor s which has the following components in 2D spherical coordinates:
Figure 7 shows the trace of s^{2}, which combines the contributions from compression and from shear, in the region below the convective boundary. Both contributions are due to the penetrative downflows that are compressed and braked by the stratification, inducing vertical and horizontal shears in this region. The region with the largest departure of the mean temperature profile compared to the initial profile coincides with the region with the largest trace of s^{2}, suggesting that compression and shear induce local heating and thermal mixing (through mixing of hot material). The coincidence is not particularly obvious for the ref model, but there is an overlap between the two regions. We note that for this model the variation of temperature is very small, of the order of 0.01%. The clear presence of these features which are intensified in the three boosted models helps identifying them in the ref model.
Fig. 7.
Radial profile of the time and space averages of the trace of s^{2}, (solid blue lines) in the overshooting layer. The red dashed lines indicate the profile of the relative difference between the time/space average of the temperature ⟨⟨T⟩_{θ}⟩_{t} and the initial profile T_{0}. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 
The additional heat injected by the downflows cannot be efficiently evacuated by radiative transport. Even in the model boost1d4, despite the enhancement of the radiative diffusivity, radiative transport remains insufficient to transport the flux excess on the timescales considered in the present simulation. At the base of the convective boundary, the radiative diffusivity κ_{rad} of the reference model is ∼3 × 10^{6} cm^{2} s^{−1} (see Fig. 1). Over a characteristic length scale H_{P}, the radiative timescale for the model boost1d4 is ∼8 × 10^{8} s at the convective boundary, which is ∼40 times greater than the simulation time (see Table 2). Heat thus accumulates and the structure is adapting in order to evacuate this heat excess. The local modification of the thermal profile for the models boost1d2 and boost1d4, and the subsequent impact on the temperature gradient and the radiative flux, are clearly seen in Fig. 6. For the models ref and boost1d1 the features are much less pronounced, but are present.
The local heating and the modification of the thermal background should proceed until the increase in radiative flux counterbalances the negative enthalpy flux and an equilibrium can be reached. The modification of the temperature gradient impacts the dynamics of the plumes. Larger subadiabaticity ∇−∇_{ad} will provide a stronger resistance to penetrating motions, since the adiabatically compressed penetrating plume will have a larger temperature excess w.r.t. its surrounding and will thus be more efficiently braked by buoyancy. This behaviour has been illustrated by studies analysing the sensitivity of the penetration process as a function of the socalled stability (or stiffness) parameter S, which measures the ratio of subadiabaticity to superadiabaticity (Hurlburt et al. 1994; Brummell et al. 2002). Conversely, smaller subadiabaticity reduces the braking effect and allows plumes to go deeper, broadening the overshooting layer and propagating the heating in deeper layers. The layers where ∇−∇_{ad} increase w.r.t. to the initial profile (larger subadiabaticity) is a region which will be able to efficiently brake the most vigorous plumes. An equilibrium is reached when the modification of the thermal profile is sufficient to increase the radiative flux and to stop the inward progression of the overshooting layer and of the local heating. There is thus a complex connection between the strength of the plumes, their thermal properties and the thermal background. This was noted by Viallet et al. (2013) when comparing the properties of convective penetration in the simulations of a red giant convective envelope and of an oxygenburning shell. They found a quasisteady state in the red giant model in which nonadiabatic processes due to radiative transport can counterbalance the effects of turbulent entrainment and stop the growth in size of the convective region.
In contrast, our results show that for the most boosted model, the local heating of the thermal profile is continuously increasing with time, with the temperature gradient getting steeper and closer to the adiabatic gradient. The thermal profile continues heating without reaching an equilibrium state within the simulated time scales, and the temperature gradient may at some point become locally unstable, as can been seen in Fig. 6 for the model boost1d4. We consider this as an extreme and artificial situation, due to combined effects, and interpret this behaviour as follows. The enhanced luminosity yields stronger penetrative plumes, which enhance the heating by compression and shears. The modification of the thermal background is more pronounced and proceeds faster due to the enhanced thermal diffusivity (and thus smaller thermal relaxation time). It increases locally the radiative flux and the heating could stop if the radiative flux could evacuate the heat excess. But it also flattens the subadiabaticity, reducing the braking of penetrative plumes and the strongest plumes progress deeper, broadening the penetration region. This is a kind of runaway situation that may stop with a significant extension of the convectively unstable region to deeper levels, which seems unrealistic for a “real” stellar model. Note that even larger enhancement factors for the luminosity, with factors > 10^{6}, will produce even more unrealistic results with convective velocities in the outer part of the domain that would become close to the speed of sound.
An artificial modification of the heat conductivity in the penetration region, as done in Brun et al. (2011, 2017), will obviously have a significant impact on the thermal profile, the temperature gradient and thus the dynamics of the penetrative plumes. Increasing the radiative transport efficiency in the overshooting region will not only accelerate the thermal relaxation process, but may also change the subtle balance mentioned above and consequently the characteristic penetration depths. There is thus a subtle equilibrium between thermal background modification and plume dynamics, which will eventually determine the maximal depths of penetrative plumes (see Sect. 4).
Since our simulations are not thermally relaxed, it is an open question whether the effects that we describe above and the modification of the thermal background that is linked to the bump of the radiative flux are transient and vanish once thermal relaxation is reached. But the bump of the radiative flux is also observed in thermally relaxed simulations (e.g., Rogers et al. 2006; Kupka et al. 2018; Käpylä 2019). Our experiments thus link the radiative flux bump to the modification of the local thermal background due to the penetrating downflows. They point the way towards a deeper understanding of the underlying physical processes responsible for this bump also observed in thermally relaxed state. One has, however, to take with caution the direct application of the overshooting depths predicted by these simulations to “real” stars since the final relaxed state for these simulations may have different properties from present non thermally relaxed states. This does not however preclude analysing the effect of luminosity and radiative diffusivity enhancement during the slowly evolving transient phase during which convection is considered to be in steady state, at least for the ref, boost1d1 and boost1d2 models.
4. Results: Statistics of convective penetration depths
In order to determine the extent of convective penetration, we adopt the approach developed in Pratt et al. (2017) based on a statistical analysis of the depth of all convective plumes that penetrate below the convective boundary. We use two criteria to determine the depth of a penetrative plume at a given angle θ and time t. The first one is given by the first zero below the convective boundary of the vertical kinetic energy flux and the second one considers instead the vertical heat flux f_{δT}(r, θ, t) = ρc_{P}(δT)v_{r}.
Pratt et al. (2017) show that the statistical analysis of penetration depths calculated for each angular grid cell at each time step reveals a nonGaussian probability distribution and the tail of the distribution, due to extreme events of plume penetration, plays a key role in defining the extent of the penetration layer. The width of the penetration layer, where mixing proceeds, should thus be better characterised by the maximal depths of plume penetration than by the average depth. Additionally, while there is significant discrepancy between the penetration depth defined from the averaged vertical kinetic energy flux (e.g., based on Eq. (9)) and the averaged vertical heat flux (e.g., based on Eq. (11)), Pratt et al. (2017) finds that the PDFs of convective penetration depths calculated from the local quantities f_{k} and f_{δT} are characteristically similar.
Pratt et al. (2017) use extreme value theory to derive the cumulative distribution function of the maximal penetration depth obtained from numerical simulations. This method allows to infer a diffusion coefficient describing the mixing driven by the convective plumes in the penetration layer. This diffusion coefficient can be used in stellar evolution codes (see e.g., Baraffe et al. 2017). The fit to the generalised extreme value distribution, which is used to model the probability of maximal events, depends in particular on a location parameter μ which is approximately equivalent to a maximal penetration depth and provides a simple overshooting length (Pratt et al. 2017, 2020). Our goal presently is not to determine diffusion coefficients but to analyse the impact of enhanced luminosity and radiative diffusivity on penetration depths. We thus focus on the comparison and characterisation of the distribution of all penetration depths and of maximal depths. We use the same criteria as Pratt et al. (2017) based on the two fluxes f_{k} and f_{δT}, respectively, to define the penetration depth of a plume. At each time t, we calculate at each angle θ the penetration depth r_{0}(θ, t) of a plume (i.e., the radial position of the first zero of f_{k} and f_{δT} below r_{conv}, the position of the convective boundary given by the position of the Schwarzschild boundary from the 1D initial model) and define the corresponding penetration length l_{0} with respect to r_{conv} by,
Note that with this definition the length l_{0} is always positive. The maximal penetration length is defined by the maximum over all angles θ:
The simulation data are sampled at a fixed time interval, typically at a rate of ∼τ_{conv}/10^{2}. Larger time intervals would not allow to capture the fastest moving plumes entering the penetration layer.
Following Pratt et al. (2017), we consider two different layers. We define a first layer of characteristic length l_{bulk} where convective plumes frequently penetrate and with l_{bulk} given by the average of all lengths l_{0}. The second, deeper, layer is characterised by convective plumes that penetrate intermittently and its length l_{max} is given by the average of all , providing an effective width of the overshooting layer. In this approach, the overshooting length l_{max}, thus defined, is expected to be close to the location parameter μ inferred from a more complete extreme value statistical analysis. Table 3 displays l_{bulk} and l_{max} in units of the total stellar radius and of the pressure scale height at the convective boundary for the four models.
Characteristic lengths l_{bulk} and l_{max}, in units of the total stellar radius and of the pressure scale height at the convective boundary, for the four models considered in this study.
Figures 8–11 show the distributions of maximal lengths based on the criterion for f_{k} and f_{δT}, respectively, for our four models. For the models boost1d1 and boost1d2, the histograms of derived from the vertical kinetic energy flux and the vertical heat flux, respectively, are similar, confirming the findings of Pratt et al. (2017). For these two models, we obtain ∼12% difference between the mean value for based on f_{k}, denoted l_{max}(f_{k}), and the one based on f_{δT}, denoted l_{max}(f_{δT}) (see Table 3). Note that when the same calculation is performed using only the first half of both simulations, i.e., N_{conv} ∼ 200 for boost1d1 and N_{conv} ∼ 220 for boost1d2, the difference in means between the two methods is ∼18% in both cases. Simulations over several hundred convective turnover times (≳400 × τ_{conv}) are thus necessary for the statistics to converge. We also calculate the location parameter μ of the generalised extreme value distribution for penetration depths based on f_{k}. We find ∼1% difference between l_{max} and μ for the model boost1d1 and ∼5% for the model boost1d2. This confirms that our approach to estimate an overshooting length from the mean of all maximal lengths is a good first approximation to the location parameter derived from the more complete statistical analysis based on extreme value theory. Note that this is valid if the goal is to perform a simple comparison between different models, but the mean of all will not reveal the wealth of information on extreme plume events contained in the generalised extreme value distribution.
Fig. 8.
Histograms of maximal penetration lengths , in units of the stellar radius R_{star}, for the reference model ref, based on the vertical kinetic energy flux f_{k} (left panel) and the vertical heat flux f_{δT} (right panel). 
Fig. 11.
Same as Fig. 8, but for the model boost1d4. Note that the bins have the same width as in Figs. 8–10. 
For the reference model, the shape of the histograms differ (see left and right panels of Fig. 8) and the difference between l_{max}(f_{k}) and l_{max}(f_{δT}) is ∼20 %. But the two distributions become more similar as the simulation time increases. In particular the distribution based on f_{k} is getting narrower with time and closer to the distribution based on f_{δT}. We stopped this simulation after ∼565 convective turnover times, considering that a difference of ≤20% between l_{max}(f_{k}) and l_{max}(f_{δT}) is satisfactory for our present purpose.
As expected from the results in Sect. 3.2, the model boost1d4 shows a different behaviour compared to the less boosted models. This model does not reach a steady state since the thermal background in the overshooting region continues to evolve and the two penetration lengths l_{bulk} and l_{max} continue increasing, indicating that the penetration layer moves deeper inwards with time. The plume depth analysis as performed here has to be taken with caution and is essentially illustrative. The histograms and overshooting lengths based on f_{k} and f_{δT}, respectively, significantly differ, and the discrepancy does not decrease with time. The overshooting length predicted by the vertical kinetic energy flux is significantly larger than the one predicted by the vertical heat flux. This suggests that the velocity magnitude v^{2} and the temperature fluctuation δT react differently in this highly boosted model compared to the less boosted ones. Strong internal waves are generated at the convective boundary by the convective motions and by the penetrating plumes, as illustrated by the rms velocity of the boost1d4 model in the radiative zone (see Fig. 3). The large discrepancy between l_{max}(f_{k}) and l_{max}(f_{δT}) may be indicative of the interference of waves with the plume dynamics, which may undermine the criterion presently used to define the depth of a plume. We will analyse in more details these effects in two followup studies devoted to internal waves (Le Saux et al., in prep.) and to the analysis of mixing in the overshooting layer based on tracer particles (Guillet et al., in prep.).
Our data are not extensive enough to robustly determine a power law scaling of the overshooting depth with the luminosity of the model, as done in e.g., Hotta (2017) and Käpylä (2019). It is however still interesting to examine the scaling suggested by our results. This is illustrated in Fig. 12 which shows the results for both l_{bulk} and l_{max}. The results for models boost1d4 are also shown for illustration in Fig. 12, but they should not be used to infer a scaling since the penetration layer progresses deeper with time. Scaling of the overshooting depth l_{ov} reported in the literature from numerical simulations can vary from l_{ov} ∝ L^{0.08} to l_{ov} ∝ L^{0.31} and strongly depends on the heat conductivity profile at the convective boundary, as shown by Käpylä (2019). The scaling also depends on the method used to define l_{max}, as illustrated in Fig. 12. The results for l_{max} do not show a simple scaling law with L. The results for l_{bulk} show an approximate scaling ∝L^{0.3}. This suggests that the dynamics of the bulk of the plumes is primarily linked to the convective velocity v_{rms} ∝ L^{1/3}, whereas the determination of l_{max} could be impacted by other factors, which will be the focus of a followup study based on tracer particles. In the next section we analyse in more details the statistics of key quantities in the overshooting layer to gain more insight about the properties of this region.
Fig. 12.
Maximal penetration length l_{max} (red) and length of the layer where the bulk of the plumes penetrate l_{bulk} (blue curves), in units of the pressure scale height at the convective boundary, as a function of the luminosity of the model. The solid curves and circles are the values derived from the use of f_{k} and the dashed curves and triangles corresponds to the values derived from f_{δT}. The dotted curves show power laws L^{0.23} and L^{0.3} as guide for the eyes. 
The results for overshooting lengths, listed in Table 3, correspond to an initial model that has differences from a solar model, as explained in Sect. 2.1. These lengths therefore should not be used more generally to characterise the overshooting in other stars.
5. Results: Statistical characteristics of the overshooting layer
A more indepth analysis of the temperature fluctuations and associated second, third and fourthorder moments can provide additional information about the flow structure in the convection and penetration regions. It also provides additional elements to understand our results and their interpretation. Because of the unsteady evolution of model boost1d4, this model is excluded from a highorder moment analysis. For the three remaining models, we calculate the standard deviation, rms, skewness and kurtosis of the temperature fluctuation, defined at time t and at each radius r by:
The temperature fluctuation δT is defined by Eq. (12). Figure 13 shows the time average of these quantities. The skewness reflects the asymmetry of the temperature fluctuation distribution around the mean value. The kurtosis provides a measure of the significance of extreme temperature fluctuations in the temperature distribution. In the convective zone, the skewness is negative due to the descending cold plumes. At the bottom of the convective zone the kurtosis increases in amplitude, indicating stronger temperature deviations. This behaviour is also due to the cold downflows which are immersed in increasingly hotter background as they descend towards the convective boundary. This is confirmed by the negative bump of the skewness where the kurtosis has its maximum. We have indicated in Fig. 13 the position of l_{bulk} and l_{max}. The position of the negative peak of the enthalpy flux coincides with the level l_{bulk} where the plumes penetrate frequently. This confirms the interpretation of this negative peak (see Sect. 3.1). It is interesting to note that the general behaviour of these highorder moments is very similar across the range of luminosity enhancement factors, both within the convective zone and at the convective boundary.
Fig. 13.
Time average of standard deviation (solid red), rms (dashdot magenta), skewness (long dash blue) and kurtosis (long dashdot black) of the temperature fluctuations for the models indicated in the panels. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. The two vertical dotted lines in each panel indicate the extension of the penetration layer l_{max} based on f_{k} and f_{δT}, respectively. Similarly, the two vertical dashed lines correspond to l_{bulk}. 
For the three models, the strongest departure between σ_{δT} and rms_{δT} lies between the location where the bulk of the plumes penetrates and the layer reached by the extreme plumes. The difference between σ_{δT} and rms_{δT} is due to the time evolution of the mean temperature profile. Indeed if ⟨T⟩_{θ} = ⟨⟨T⟩_{θ}⟩_{t}, then the standard deviation and the rms are equivalent. This strong departure thus indicates the modification of the temperature background due to the slight heating mentioned in Sect. 3.2. The depth reached by the most vigorous plumes l_{max} coincides with the position where σ_{δT} and rms_{δT} get closer again, i.e., where the background is not significantly evolving. A closer look at the penetration layers shows that the position of l_{max} is located in the region where the subadiabaticity ∇−∇_{ad} increases compared to the initial profile, implying more efficient braking of descending downflows. This is illustrated in Fig. 14 which shows the profile of ⟨⟨(∇ − ∇_{ad})⟩_{θ}⟩_{t} − (∇−∇_{ad})_{init} for the three models. This is consistent with the picture described in Sect. 3.2 and highlights the importance of the thermal profile modified by the penetrating plumes and affecting in return their dynamics and the effective overshooting depth.
Fig. 14.
Solid lines (red): Departure of the subadiabaticity (time and horizontal average) from its initial value δ = ⟨⟨(∇ − ∇_{ad})⟩_{θ}⟩_{t} − (∇−∇_{ad})_{init} as a function of radius in the overshooting region for the three models, as indicated in the panels. Dashed lines (blue): radial profile of ⟨cor⟩_{t}, the time average of the radial velocity autocorrelation defined by Eq. (23). Since these quantities are very small, they are multiplied by different constants, as indicated in each panel, in order to see the details. The vertical lines are the same as in Fig. 13. 
The connections found between the position of l_{max} and the behaviour of σ_{δT} and rms_{δT} and between l_{max} and the subadiabaticity are of particular interest and will be the focus of future studies as these quantities may provide other criteria to determine the depth of the overshooting layer. The statistical approach performed in Sect. 4 requires simulation times covering several hundreds of convective turnover timescales (≳400 × τ_{conv}) to reach statistical convergence. Criteria based on the behaviour of σ_{δT} and rms_{δT} and/or on the profile of (∇ − ∇_{ad}) would not require such long simulation times and would provide an advantageous alternative approach, in terms of computation time.
Korre et al. (2019) and Cai (2020) suggested to use the radial velocity correlation with the velocity at the convective boundary as an alternative measure of the extent of the overshooting distance. They argue that this quantity should favour the strongest downflows and capture extreme events. Cai (2020) suggests that the secondzero of this velocity autocorrelation is a good indicator for the measurement of upward overshooting. We examine this suggestion and, as Cai (2020), we define the velocity autocorrelation at a given time and at each radius r by:
with v_{r, rms} the rms of the radial velocity following the definition given in Eq. (7) and r_{conv} the radius at the convective boundary. We show the radial profile of the time average ⟨cor⟩_{t} in Fig. 14 and compare it to the position of l_{max} and the profile of the subadiabaticity. We find that the velocity autocorrelation has several zeros between l_{bulk} and l_{max}. The second zero lies between these two lengths, but its position shows no systematic link with the position of l_{max} and with the behaviour of the subadiabaticity. A more thorough analysis of the most relevant criterion to characterise the mixing in the overshooting layer and its length is beyond the scope of this study and will be the focus of a followup study.
6. Discussion and conclusion
The experiments with artificially enhanced luminosity and thermal diffusivity performed in this study have several implications. Our results suggest a link between the properties of penetrating plumes and the thermal structure of the overshooting layer; this confirms the findings of earlier works (e.g., Hurlburt et al. 1986; Muthsam et al. 1995; Rempel 2004; Käpylä 2019). They stress the importance of the impact of the penetrative downflows on the local thermal background, as a result of compression and shear which induce local heating and thermal mixing in the overshooting layer. This effect is not pronounced in the ref model, but is present. The artificial increase in the energy flux intensifies the effect by increasing the velocities in the convective zone and at the convective boundary. This artefact reveals the subtle connection between the local heating of the thermal background and the plume dynamics as the heating alters the subadiabaticity in the penetration region. This local heating increases also the efficiency of the heat transport by radiation, which may counterbalance further heating and help establishing a steady state. We suggest that the modification of the thermal background by the penetrative plumes and its impact on the radiative flux and on the plume dynamics play a role on the efficiency of the mixing driven by the plumes and eventually on the width of the overshooting layer.
This finding suggests that an artificial modification of the radiative diffusivity in the overshooting layer, rather than only accelerating the thermal relaxation, may also alter this balance and most likely affect the efficiency of the overshooting mixing and thus the effective width of the overshooting layer. In the most boosted model presently studied, namely boost1d4, this balance seems difficult to establish because of the strong plumes driving efficient heating and progression of the penetration layer deeper inwards. In the context of the study of overshooting, the argument that a simulation with enhanced energy flux and radiative diffusivity compared to a reference model is an accelerated version of this model, or in other words that the boost1d4 simulation at a given time provides a good representation of the ref simulation at much later times, is difficult to justify in light of present results.
Admittedly, stellar hydrodynamic simulations based on realistic energy flux and radiative diffusivity still remain very far away from realistic stellar conditions. For models with no or moderate enhancement of the luminosity, the numerical diffusion due to truncation errors likely dominates the physical radiative diffusion. Consequently, the effective Péclet number, which measures the relative importance of advection versus diffusion, is set by the numerical diffusion and thus underestimated. For large enhancement of the luminosity and of the radiative diffusivity, as in the boost1d4 model, the physical radiative diffusion can start to dominate over the numerical diffusion. In these models, the heat diffusion is hence better described. But unfortunately, these models drift away strongly from their initial state and their relaxed state will not describe the initial one. Artificial enhancement of the energy flux or modification of the thermal diffusivity profile therefore pushes the simulated conditions away from the original target star. Based on present results, we do not advise using luminosity enhancement factors ≳10^{4} as they induce a significant drift from the initial stellar structure. Factors even larger (> 10^{6}) could produce unrealistic supersonic velocities close to the stellar surface.
The modification of the local thermal background in the overshooting region has been reported in several works (e.g., Rogers et al. 2006; Viallet et al. 2013; Korre et al. 2019; Higl et al. 2021), but has not been the subject of deeper scrutiny to understand its origin. This lack of deepening is most likely due to the fact that many of these simulations, including our work, have not achieved thermal relaxation, raising a concern about the physical relevance of this feature. Paradoxically, many numerical simulations acknowledge the existence of a bump of the radiative flux in the penetration region (e.g., Brun et al. 2011, 2017; Hotta 2017; Käpylä 2019; Cai 2020), needed to counterbalance the negative enthalpy flux. But no details are given regarding what causes this bump, except that it is there to ensure thermal equilibrium. There must be physical processes which yield such a peak in the radiative flux. Our simulations also feature this bump (see Fig. 4). Several thermal relaxation timescales are required in order for the radiative flux to compensate exactly for the negative enthalpy flux. Our reference and moderately boosted models are not able to reach this state, whereas as discussed above, the most boosted simulation is evolving towards an extreme state. Because our experiments link the radiative flux bump to the modification of the local thermal background due to the penetrating downflows, they could help understanding the underlying physical processes responsible for this bump.
The sensitivity of this feature to some of the input physics included in the models or to the dimensionality of the numerical simulations should be examined. In order to test the robustness of the local heating process against the equation of state, we constructed an artificial 1D stellar model close to the stage of evolution of the ref model but with an ideal gas equation of state with constant adiabatic index γ = 5/3. Such a simple equation of state is implemented in the MUSIC code and is straightforward to implement in a 1D stellar evolution code. This provides an unrealistic stellar model because γ varies in the stellar interior due to recombination or ionisation processes in the convective envelope and in the overshooting region. Indeed, the radiative–convective transition at the bottom of the convective envelope corresponds to the last bump in opacity due to partial ionisation of metals (Rogers & Iglesias 1992). However, such a model provides an interesting case to test the sensitivity of the heating process to the equation of state and, for the most boosted model, the fast inward progression of the penetration layer. We ran two cases with the MUSIC code starting from this artificial stellar model with luminosity enhanced by a factor 10 and 10^{4}, respectively, and find the same modification of the thermal profile and same heating in the overshooting region, and for the most boosted model, the same inward progression of the penetration layer with time as found for models with a realistic stellar equation of state.
We performed an additional test to check the robustness of the results against the treatment of gravity. Given the nonnegligible modification of the background profile found in the boosted models, we have verified for the boost1d1 and boost1d4 models that updating the calculation of the spherically symmetric gravitational acceleration g during the simulation has no impact on the dynamics of the plumes and on the heating process. Since there is no need to recalculate g at each numerical timestep, it was recalculated every time interval Δt, with Δt fixed during a whole test simulation. The typical dynamical timescale of the stellar model , with ρ_{mean} the mean density and G the gravitational constant, is of the order of 10^{3} s. We have thus performed a number of test simulations with Δt varied between 10^{2} and 10^{5} s and found no impact on our main results.
One may also question the impact of the resolution and of the dimensionality on the local modification of the thermal background. We performed two test simulations for boost1d1 and boost1d4 with double resolution (1024 × 1024) over a few convective turnover timescales and find the same local heating and modification of the thermal background over similar timescales compared to the 512 × 512 counterpart. Even though such resolution tests will not bring the simulations close to the true parameter regime, they remain useful as they provide a test for the sensitivity of the results to the numerical diffusion.
We are currently performing threedimensional simulations of a solarlike model but they are not yet exploitable to analyse the time evolution of the thermal background in the overshooting layer. In the meantime, we can only speculate on the impact of dimensionality on these results. Brummell et al. (2002) find different structures of the penetrative convection in 3D compared to 2D simulations. They still find significant overshoot of the convective motions in highly turbulent 3D simulations, but these motions do not establish an adiabatic penetration region as previously found in 2D simulations (e.g., Hurlburt et al. 1994). Even if the structure and the geometry of the penetrating downflows are modified in 3D, these motions should still produce compression and shear in the overshooting region. We should thus still expect them to induce local heating and thermal mixing, and consequently a modification of the thermal background. We note that the alteration of the thermal background has been reported in other numerical simulations based on very different numerical approaches, namely a 2D anelastic method in Rogers et al. (2006) and a 3D approach within the Boussinesq approximation in Korre et al. (2019). The resolution and dimensionality may however impact the establishment of the balance between the local heating, the increase in the radiative flux and the inward progression of the penetration layer. We suspect that this balance is subtle in numerical simulations and could be impacted by various factors and uncertainties. We are not able to be more specific based on the present study and this is the subject of future experiments we are performing for convective envelopes and convective cores.
Simulations with an artificial increase in the energy flux may be of limited relevance for the quantitative determination of an overshooting distance. Because of additional artefacts introduced for large enhancement factors, scaling laws based on such simulations need to be applied with caution. But such simulations are still valuable experiments to perform for moderate enhancement factors. They may be used for the general analysis of convective flow structure based on highorder statistics of temperature fluctuations, given the similarities found across a range of modest luminosity enhancement factors. These experiments may also be used for the predictions of internal wave spectrum and properties, modulo an appropriate rescaling. This has never been proven or disproven and is the purpose of a followup study (Le Saux et al., in prep.). They also provide excellent laboratories to compare a method based on fluxes in an Eulerian approach and a method based on Lagrangian tracer particles to define an overshooting depth characterising material mixing.
Lastly, as suggested by e.g., Rempel (2004) and Zhang et al. (2012), the modification of the thermal background at the convective boundary due to convective penetration could improve the sound speed discrepancy between helioseismology data and solar models. Our results suggest a complex behaviour of the temperature gradient at the convective boundary that deserves further attention, not only in the context of helioseismology but also for the general understanding of convective penetration.
Acknowledgments
This work is supported by the ERC grant No. 787361COBOM and the consolidated STFC grant ST/R000395/1. IB thanks the Max Planck Institut für Astrophysics (Garching) for warm hospitality during completion of part of this work. The authors would like to acknowledge the use of the University of Exeter HighPerformance Computing (HPC) facility ISCA and of the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility. The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National eInfrastructure.
References
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2020, MNRAS, 491, 972 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., & El Eid, M. F. 1991, A&A, 245, 548 [NASA ADS] [Google Scholar]
 Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403 [Google Scholar]
 Baraffe, I., Pratt, J., Goffrey, T., et al. 2017, ApJ, 845, L6 [Google Scholar]
 Biermann, L. 1932, ZAp, 5, 117 [NASA ADS] [Google Scholar]
 Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [Google Scholar]
 Cai, T. 2020, ApJ, 888, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Ludwig, H. G., & Steffen, M. 1996, A&A, 313, 497 [NASA ADS] [Google Scholar]
 Goffrey, T., Pratt, J., Viallet, M., et al. 2017, A&A, 600, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Higl, J., Müller, E., & Weiss, A. 2021, A&A, 646, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&A, 641, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
 Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J.P. 1994, ApJ, 421, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, S., Andrassy, R., Sandalski, S., et al. 2017, MNRAS, 465, 2991 [Google Scholar]
 Käpylä, P. J. 2019, A&A, 631, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Rheinhardt, M., & Brandenburg, A. 2017, ApJ, 845, L23 [CrossRef] [Google Scholar]
 Korre, L., Garaud, P., & Brummell, N. H. 2019, MNRAS, 484, 1220 [Google Scholar]
 Kupka, F., Zaussinger, F., & Montgomery, M. H. 2018, MNRAS, 474, 4660 [NASA ADS] [CrossRef] [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [Google Scholar]
 Muthsam, H. J., Goeb, W., Kupka, F., Liebich, W., & Zoechling, J. 1995, A&A, 293, 127 [Google Scholar]
 Porter, D. H., & Woodward, P. R. 2000, ApJS, 127, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Pratt, J., Baraffe, I., Goffrey, T., et al. 2016, A&A, 593, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pratt, J., Baraffe, I., Goffrey, T., et al. 2017, A&A, 604, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pratt, J., Baraffe, I., Goffrey, T., et al. 2020, A&A, 638, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rempel, M. 2004, ApJ, 607, 1046 [Google Scholar]
 Rogers, F. J., & Iglesias, C. A. 1992, ApJ, 401, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
 Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1 [CrossRef] [Google Scholar]
 Rogers, T. M., Glatzmaier, G. A., & Jones, C. A. 2006, ApJ, 653, 765 [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Schmitt, J. H. M. M., Rosner, R., & Bohn, H. U. 1984, ApJ, 282, 316 [NASA ADS] [CrossRef] [Google Scholar]
 Shaviv, G., & Salpeter, E. E. 1973, ApJ, 184, 191 [Google Scholar]
 Tian, C.L., Deng, L.C., & Chan, K.L. 2009, MNRAS, 398, 1011 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Baraffe, I., & Walder, R. 2011, A&A, 531, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Goffrey, T., Baraffe, I., et al. 2016, A&A, 586, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zahn, J. P. 1991, A&A, 252, 179 [Google Scholar]
 Zhang, C., Deng, L., Xiong, D., & ChristensenDalsgaard, J. 2012, ApJ, 759, L14 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Properties of the initial reference model: mass, luminosity, radius, depth of the convective envelope and pressure scale height at the Schwarzschild boundary.
Characteristic lengths l_{bulk} and l_{max}, in units of the total stellar radius and of the pressure scale height at the convective boundary, for the four models considered in this study.
All Figures
Fig. 1.
Internal structure of the reference model (solid lines) used as initial model for the 2D simulations compared to a model for the current Sun (dashed lines). The figure shows the density ρ (in g cm^{−3}), temperature T (in K) and radiative diffusivity κ_{rad} = χ/(ρc_{P}) (in cm^{2} s^{−1}), with c_{P} the specific heat at constant pressure. The Schwarzschild boundary is indicated by the vertical solid line for the reference model and the vertical dashed line for the Sunlike model. 

In the text 
Fig. 2.
Evolution of the total kinetic energy (in erg) as a function of time (in s) of the simulations described in Table 2. The dotted lines correspond to the value of the total kinetic energy at the beginning of the steady state for convection. 

In the text 
Fig. 3.
Top panel: radial profile of the time averaged rms velocity (solid lines) and rms radial velocity (dashed lines). Velocities are in cm/s. The curves from bottom to top correspond to the four models: ref (black), boost1d1 (blue), boost1d2 (magenta) and boost1d4 (red), respectively. Bottom panel: velocities scaled by the enhancement factor (L/L_{star})^{1/3}. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 

In the text 
Fig. 4.
Radial profile of luminosities for the four models as indicated in each panel, corresponding to various fluxes: vertical enthalpy flux (L_{h}, dashdot, blue), radiative flux (L_{rad}, long dashdot, black), vertical kinetic energy flux (L_{k}, dot, red), convective heat flux (L_{δT}, dash, magenta). The total luminosity L_{tot} = L_{h} + L_{rad} + L_{k} (solid, green) is also displayed. Luminosities are divided by the luminosity of the corresponding model. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 

In the text 
Fig. 5.
Visualisation of the radial velocity and the temperature fluctuations close to the convective boundary (horizontal black line) for the ref model at a given time. The xaxis represents the colatitude (in terms of cos θ). 

In the text 
Fig. 6.
Radial profile of time and space averages of various quantities close to the convective boundary. From top to bottom: the temperature T (in units of 6 × 10^{6} K); the radiative conductivity χ divided by the luminosity enhancement factor of the model and by the constant 4 × 10^{15}; the radiative flux divided by the luminosity enhancement factor of the model and by the constant 2 × 10^{12}; the subadiabaticity ∇ − ∇_{ad}; the radial temperature gradient ∂T/∂r multiplied by the constant 2.8 × 10^{3}. The coloured curves correspond to the four models: ref (black), boost1d1 (blue), boost1d2 (magenta) and boost1d4 (red). The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 

In the text 
Fig. 7.
Radial profile of the time and space averages of the trace of s^{2}, (solid blue lines) in the overshooting layer. The red dashed lines indicate the profile of the relative difference between the time/space average of the temperature ⟨⟨T⟩_{θ}⟩_{t} and the initial profile T_{0}. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. 

In the text 
Fig. 8.
Histograms of maximal penetration lengths , in units of the stellar radius R_{star}, for the reference model ref, based on the vertical kinetic energy flux f_{k} (left panel) and the vertical heat flux f_{δT} (right panel). 

In the text 
Fig. 9.
Same as Fig. 8, but for the model boost1d1. 

In the text 
Fig. 10.
Same as Fig. 8, but for the model boost1d2. 

In the text 
Fig. 11.
Same as Fig. 8, but for the model boost1d4. Note that the bins have the same width as in Figs. 8–10. 

In the text 
Fig. 12.
Maximal penetration length l_{max} (red) and length of the layer where the bulk of the plumes penetrate l_{bulk} (blue curves), in units of the pressure scale height at the convective boundary, as a function of the luminosity of the model. The solid curves and circles are the values derived from the use of f_{k} and the dashed curves and triangles corresponds to the values derived from f_{δT}. The dotted curves show power laws L^{0.23} and L^{0.3} as guide for the eyes. 

In the text 
Fig. 13.
Time average of standard deviation (solid red), rms (dashdot magenta), skewness (long dash blue) and kurtosis (long dashdot black) of the temperature fluctuations for the models indicated in the panels. The convective boundary corresponding to the Schwarzschild boundary from the 1D initial model is indicated by the vertical solid line. The two vertical dotted lines in each panel indicate the extension of the penetration layer l_{max} based on f_{k} and f_{δT}, respectively. Similarly, the two vertical dashed lines correspond to l_{bulk}. 

In the text 
Fig. 14.
Solid lines (red): Departure of the subadiabaticity (time and horizontal average) from its initial value δ = ⟨⟨(∇ − ∇_{ad})⟩_{θ}⟩_{t} − (∇−∇_{ad})_{init} as a function of radius in the overshooting region for the three models, as indicated in the panels. Dashed lines (blue): radial profile of ⟨cor⟩_{t}, the time average of the radial velocity autocorrelation defined by Eq. (23). Since these quantities are very small, they are multiplied by different constants, as indicated in each panel, in order to see the details. The vertical lines are the same as in Fig. 13. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.