#  Impact of radial truncation on global 2D hydrodynamic simulations for a
Sun-like model

D. G. Vlaykov,1 I. Baraffe,1,2 T. Constantino,1 T. Goffrey,3 T. Guillet,1 A.
Le Saux,1,2 A. Morison,1 and J. Pratt4
1 Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4
4QL,UK
2 École Normale Supérieure, Lyon, CRAL (UMR CNRS 5574), Université de Lyon,
France
3 Centre for Fusion, Space and Astrophysics, Department of Physics, University
of Warwick, Coventry, CV4 7AL, UK
4 Department of Physics and Astronomy, Georgia State University, Atlanta, GA
30303, USA<EMAIL_ADDRESS>
(Accepted 2022 April 29. Received 2022 April 29; in original form 2022
February 22)

###### Abstract

Stellar convection is a non-local process responsible for the transport of
heat and chemical species. It can lead to enhanced mixing through convective
overshooting and excitation of internal gravity waves (IGWs) at convective
boundaries. The relationship between these processes is still not well
understood and requires global hydrodynamic simulations to capture the
important large-scale dynamics. The steep stratification in stellar interiors
suggests that the radial extent of such simulations can affect the convection
dynamics, the IGWs in the stably stratified radiative zone, and the depth of
the overshooting layer. We investigate these effects using two-dimensional
global simulations performed with the fully compressible stellar hydrodynamics
code MUSIC. We compare eight different radial truncations of the same solar-
like stellar model evolved over approximately $400$ convective turnover times.
We find that the location of the inner boundary has an insignificant effect on
the convection dynamics, the convective overshooting and the travelling IGWs.
We relate this to the background conditions at the lower convective boundary
which are unaffected by the truncation, as long as a significantly deep
radiative layer is included in the simulation domain. However, we find that
extending the outer boundary by only a few percent of the stellar radius
significantly increases the velocity and temperature perturbations in the
convection zone, the overshooting depth, the power and the spectral slope of
the IGWs. The effect is related to the background conditions at the outer
boundary, which are determined in essence by the hydrostatic stratification
and the given luminosity.

###### keywords:

stars: interiors, stars: solar-type, hydrodynamics, convection, waves

††pubyear: 2022††pagerange:  Impact of radial truncation on global 2D
hydrodynamic simulations for a Sun-like model–LABEL:lastpage

##  1 Introduction

Stellar convection is a non-local process (Rieutord & Zahn, 1995; Spruit,
1997; Brandenburg, 2016), which can affect stellar structure and evolution
(Spruit, 1997; Gough, 1977; Canuto & Dubovikov, 1997). Geometrical parameters,
such as the aspect ratio and vertical extent can play an important dynamical
role in simulations of boundary-driven convection (Hurlburt et al., 1984;
Cossette & Rast, 2016; Pratt et al., 2016; Wagner & Shishkina, 2013). This
work aims to further qualify this role, by studying the effects of varying the
radial extent in global two-dimensional hydrodynamic simulations of a Sun-like
stellar model on the dynamics of the convective and radiative zones and the
boundary between them.

Global stellar hydrodynamic simulations are computationally challenging
because of the extreme dynamical ranges of length, time, and amplitude scales
found in stellar interiors. Typically, the geometrical domain is restricted
radially at some point below the stellar photosphere and/or above the stellar
core. For instance in global solar simulations, Rogers & Glatzmaier (2006) set
the radial range at
$$0.001\text{\,}\mathrm{R}_{\sun}$-$0.93\text{\,}\mathrm{R}_{\sun}$$ and Brun
et al. (2011) use the radial range
$$0.07\text{\,}\mathrm{R}_{\sun}$-$0.97\text{\,}\mathrm{R}_{\sun}$$. When
focusing on the solar convection zone typical radial ranges are
$$0.72\text{\,}\mathrm{R}_{\sun}$-$0.97\text{\,}\mathrm{R}_{\sun}$$ (Augustson
et al., 2015), or when including convective boundary effects,
$$0.71\text{\,}\mathrm{R}_{\sun}$-$0.96\text{\,}\mathrm{R}_{\sun}$$ (Hotta &
Kusano, 2021) or
$$0.61\text{\,}\mathrm{R}_{\sun}$-$0.96\text{\,}\mathrm{R}_{\sun}$$ (Guerrero
et al., 2013). While some authors consider the effect of boundary conditions
(e. g. by including an outer cooling layer Guerrero et al. 2013; Hotta et al.
2019), the isolated effect of the radial truncation itself has rarely been
investigated to our knowledge (see below). However, both the inner and outer
boundary layers can house important processes with both global and long-term
effects.

The temperature and density stratification in the outer convective layers can
be very steep, dropping by an order of magnitude over a percent of a stellar
radius. Thus a small change in radius leads to a large change in the
background thermodynamic conditions. Due to the inherent non-locality of
boundary-driven convection (Gough, 1977; Spruit, 1997; Canuto & Dubovikov,
1997), this is expected to impact the global dynamics of the simulation
domain. In fact, it has been suggested that the outer solar layers are the
dominant driver of the convective zone (Spruit, 1997; Brandenburg, 2016)
through an ‘entropy rain’ of radiatively cooled dense perturbations111 Note
that this study does not attempt to address the viability of the entropy rain
framework, because the considered simulation domain does not extend to regions
where radiative cooling becomes significant. We focus instead on isolating the
effect of radial extension of the region with efficient convection..

The lower convective boundary is also of particular interest. Due to inertia,
flows from the convective zones can penetrate some distance across the
boundary into the stably-stratified radiative zone. This process has been
referred to as convective penetration, overshooting or convective boundary
mixing (Zahn, 1991; Hurlburt et al., 1994; Brummell et al., 2002) in different
contexts. Here we follow the terminology in Brummell et al. (2002) and use the
term _overshooting_ to indicate convective flows traveling through sub-
adiabatic layers. Despite occurring in a small stellar region (by volume and
mass), convective overshooting can have a quantitative impact on the stellar
evolution. In general, the dynamics at the lower convective boundary are
central to a number of open questions of stellar physics relating to mixing
and transport of material, angular momentum, and magnetic fields. For
instance, the enhanced mixing of material from the well-mixed cooler
convective zone across the lower convective boundary can help deplete the
surface lithium abundance in pre-main-sequence and main-sequence stars
(Baraffe et al., 2017; Constantino et al., 2021). Convective motions along and
across the boundary excite internal gravity waves (IGWs) in the stably-
stratified radiative zone (RZ) (Press, 1981; Goldreich & Kumar, 1990; Dintrans
et al., 2005; Pinçon et al., 2016). IGWs in turn are a sensitive probe of
stellar interior structure, observable in massive stars (Bedding et al., 2010;
Rogers et al., 2013; Bowman et al., 2019), and are involved in the
redistribution of angular momentum (Rogers et al., 2013) and chemical mixing
(Rogers & McElwaine, 2017). Finally, there is an active debate regarding the
amplitude of the large-scale convective velocities in solar interior, the so-
called _solar convective conundrum_. Different helioseismological techniques
point to a bigger (Hanasoge et al., 2012) or smaller (Birch et al., 2018)
discrepancy compared to the typical values obtained in numerical simulations.
Sometimes considered as part of the convective conundrum is also the problem
of reproducing the solar-type of differential rotation in the convection zone,
i. e. fast equator and slow poles222Recent advances suggest that
magnetohydrodynamic simulations with sufficiently high resolution can
reproduce the solar type of differential rotation (Hotta & Kusano, 2021)..

The solutions to such questions require stellar hydrodynamic simulations,
which can capture non-local processes in both space and scale. While three-
dimensional (3D) simulations can be more realistic, their large computational
cost makes extensive parameter studies challenging, especially over long time
intervals. The relaxation time scales may be shortened by artificially
boosting the stellar luminosity and thermal diffusivity. However, running such
simulations self-consistently requires that the convection zone be adiabatic,
to avoid changes in the temperature stratification due to the boosting
(Baraffe et al., 2021). It has also been shown that boosting can have a
significant impact on the dynamics throughout the stellar interior (Baraffe et
al., 2021; Le Saux et al., 2021). Therefore, in this study, we do not resort
to boosted simulations and consider instead two-dimensional (2D) models with
realistic stellar luminosity covering a wide selection of radial truncations.

Several studies have examined the influence of the radial extent using
different sets of approximations. Hurlburt et al. (1984) performed 2D fully
compressible simulations on a Cartesian grid with polytropic stratification
and ideal gas law equation of state (EoS) and found that mean kinetic energy
grows with the density contrast in the domain (which is analogous to the depth
of the convection zone). Cossette & Rast (2016) performed 3D anelastic
simulations on a Cartesian grid with polytropic stratification. They reported
that the super-adiabaticity in the boundary layer affects both the amplitude
of the temperature fluctuations in the convective bulk and the horizontal
length-scale of the velocity field. Such a change in super-adiabaticity at the
outer boundary occurs naturally when varying the radial extent of the
simulation domain using realistic stellar stratification. The reason is that
in this case hydrostatic equilibrium dictates a sharp drop in density close to
the stellar surface. This results in a decrease of the total heat capacity of
the outer layers and thus the amount of luminosity that the convective heat
flux can transport333With realistic stellar EoS the specific heat at constant
pressure is not sensitive enough to the temperature to counterbalance the
effect of the decreasing density.. As a result the temperature stratification
must become steeper than adiabatic, so that the radiative flux can transport
the remaining luminosity. Pratt et al. (2016) performed 2D fully compressible
global simulations with a realistic EoS for a low mass star, opacity and
stratification. They studied the effects on convection and overshooting due to
grid resolution and the radial truncation of the simulation domain in a low-
mass pre-main sequence star, a ‘young sun’. They found that the radial extent
of the simulation domain and the coupling at the upper and lower convective
boundaries has a pronounced effect on the dynamics of the convective zone and
the depth of the overshooting layer. Finally, Hotta et al. (2019) considered
both local Cartesian and global spherical simulations in 3D using a reduced
speed of sound technique and a combination of linear and tabulated EoS. They
found little difference in the convective bulk when comparing the effects on
the convection zone from (i) including a photospheric outer layer with
resolved radiative transport and (ii) adding artificial cooling with a
Gaussian profile around the top of the convection zone. They also found that
local Cartesian and global spherical simulation produce similar results, and
indicate that including a photospheric layer is unlikely to help solve the
solar convective conundrum.

Building on these results, we perform global 2D fully-compressible simulations
of a Sun-like star (with solar mass, luminosity and metallicity) on the main
sequence (at $\sim$4.6\text{\,}\mathrm{G}\mathrm{y}\mathrm{r}$$) using a
realistic EoS, opacity and initial stratification. In order to isolate the
impact of the radial extent of the outer convective and inner radiative zones
we do not consider the effects of rotation, magnetic fields, and cooling from
atmospheric layers. Similarly to Pratt et al. (2016), we study the impact on
the convective velocity and the depth of the overshooting layer in 2D
simulations. In addition, we measure the temperature fluctuations as well as
the impact on the generated IGWs and their potential feedback on the
overshooting. We highlight the role of the outer boundary layers and how they
can qualitatively affect the entire convection zone and the overshooting
region. We use a main-sequence Sun-like model, which has a lower luminosity, a
significantly shallower convection zone and much steeper density
stratification in the radiative layers than the Pratt et al. (2016) model. To
our knowledge this is the first study of its type applied to a Sun-like
stellar model on the main sequence.

The paper is organised as follows. In Section 2 we describe the numerical
methods and the simulations. We examine and discuss the impact of the radial
truncation on the convection zone dynamics, (Section 3), the overshooting
(Section 4) and the IGWs in the radiative zone (Section 5). We conclude with a
summary of the results and a discussion of their implications.

##  2 Methods

###  2.1 Numerics: MUSIC

We use the fully compressible time-implicit code MUSIC. The code is described
in detail in Viallet et al. (2011); Viallet et al. (2016); Goffrey et al.
(2017). Here, we provide a brief description of the main features used in this
study. MUSIC solves the inviscid Euler equations in the presence of thermal
diffusion and external gravity,

$\displaystyle\partial_{t}\left(\rho\right)$
$\displaystyle=-\nabla\cdot\left(\rho\boldsymbol{v}\right),$ (1)
$\displaystyle\partial_{t}\left(\rho\boldsymbol{v}\right)$
$\displaystyle=-\nabla\cdot\left(\rho\boldsymbol{v}\boldsymbol{v}\right)-\nabla
p+\rho\boldsymbol{g},$ (2) $\displaystyle\partial_{t}\left(\rho e\right)$
$\displaystyle=-\nabla\cdot\left(\rho
e\boldsymbol{v}\right)-p\nabla\cdot\boldsymbol{v}+\nabla\cdot\left(\chi\nabla
T\right),$ (3)

with density $\rho$, velocity $\boldsymbol{v}$, specific internal energy
density $e$, gas pressure $p$, gravitational acceleration $\boldsymbol{g}$,
and thermal conductivity $\chi$. In the presented simulations, the thermal
conductivity is dominated by radiative transfer:

$\displaystyle\chi=\frac{16\sigma T^{3}}{3\kappa\rho},$ (4)

with $\kappa$ the Rosseland mean opacity, and $\sigma$ the Stefan-Boltzmann
constant. We use realistic opacities and equation of state for stellar
interiors. The 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 solar-like interior structures.

We do not include explicit viscosity in the equations, because its value in
stellar interiors is too low to dominate over numerical viscosity for
computationally tractable grid resolutions. Instead we rely on the numerical
dissipation of the model and interpret the results in the context of the
implicit large-eddy simulation (ILES) paradigm.

###  2.2 Description of the simulations.

Figure 1:  Initial background radial profile in density, temperature, pressure
and entropy. The first three are normalised to their values at $0.4R_{\ast}$
to better compare the stratification in the upper layers. The solid vertical
line shows the boundary between the convective and radiative zone. The dotted
vertical lines on the left show the truncation radii for $r_{\min}$, the
dashed vertical lines on the right show the truncation radii for $r_{\max}$.

The initial background state is defined by the density and temperature
stratifications shown in Fig. 1 with
$R_{\ast}=$7.0775\text{\times}{10}^{10}\text{\,}\mathrm{c}\mathrm{m}$=$1.0169\text{\,}\mathrm{R}_{\sun}$$,
$L_{\ast}=$3.792\text{\times}{10}^{33}\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\,\mathrm{s}^{-1}$=$0.98772\text{\,}\mathrm{L}_{\sun}$$
and
$M_{\ast}=$1.9891\text{\times}{10}^{33}\text{\,}\mathrm{g}$=$1\text{\,}\mathrm{M}_{\sun}$$,
where $\odot$ represents the solar reference values. The initial conditions
are generated with the one-dimensional Lyon stellar evolution code Baraffe &
El Eid (1991); Baraffe et al. (1998), using the same opacities and equation of
state as implemented in MUSIC. Note that the lower boundary of the convection
zone in the model is at $r_{\textrm{Schw}}\sim$0.72$R_{\ast}$ as defined by
the Schwarzschild criterion, which gives the thermal stratification required
for a convective instability. Explicitly, the Schwarzschild criterion defines
the convective boundary via the super-adiabaticity being non-negative, i. e.
$\nabla-\nabla_{\mathrm{ad}}\geq 0$. Here, $\nabla=\mathrm{d}\log
T/\mathrm{d}\log p$ is the temperature gradient and
$\nabla_{\mathrm{ad}}=\mathrm{d}\log T/\mathrm{d}\log p|_{S}$ is the adiabatic
temperature gradient, i. e. at constant entropy $S$. In practice in the
convection zone the super-adiabaticity is very close to zero, as indicated by
the quasi-isentropic profile shown in Fig. 1. Note however that the convection
zone is not precisely adiabatically stratified with super-adiabaticity
reaching $O(10^{-2})$ close to the surface. A close analogue of this initial
model has been used in Baraffe et al. (2021); Le Saux et al. (2021) to study
the effects of luminosity boosting. This requires minimising the super-
adiabaticity in the convection zone (so that stratification is not affected by
the luminosity boosting). A comparison between the two models shows that this
leads to a significant change in the total stellar radius and the lower
convective boundary and a slight increase in luminosity.

Table 1 summarises the key parameters of the eight performed simulations. We
consider two lower boundaries ($r_{\min}/R_{\ast}\in[0.4,0.6]$) and four upper
boundaries ($r_{\max}/R_{\ast}\in[0.9,0.94,0.97,0.99]$). The one-dimensional
profile is extruded symmetrically in the angular direction. The simulation
domain covers the full angular range of $\theta\in[$$,$$]$ with angular
resolution of $\Delta\theta=$$$. We use a fixed radial resolution for all
simulations of $\Delta
r=$6\text{\times}{10}^{-4}$R_{\ast}\approx$420\text{\,}\mathrm{k}\mathrm{m}$$.
The radial cell size is chosen to ensure good representation of the pressure
scale height $H_{p}(r)=-\partial r/\partial\ln(p(r))$ even at the outermost
simulated layers. Thus, at the outer boundaries the radial cell size
corresponds to $\Delta r=0.225H_{p}(0.99R_{\ast})$ and $\Delta
r=0.017H_{p}(0.9R_{\ast})$. The large range of $\Delta r/H_{p}(r_{\max})$ is
due to the strong stratification which leads the pressure scale height to
decrease sharply in the top $\sim 10\%$ of the stellar radius, as can be
inferred from Fig. 1. At the lower convective boundary, the radial cell size
corresponds to $\Delta r=0.0075H_{p}(r_{\textrm{Schw}})$.

At the polar boundaries ($\theta=$$$ and $\theta=$$$), the normal derivative
is set to zero for all fields. At the radial boundaries, the boundary
conditions are reflective for the velocity. The density and internal energy
are linearly extrapolated in the ghost cells used to calculate the spatial
gradients and fluxes. This condition is well-suited to the stratified
background. Due to the use of staggered grids, the radiative flux at the
radial boundary must also be given. Here, it is set to ensure a constant
uniform energy flux corresponding to the respective luminosity of the 1D model
at that location. In practice, this makes for a slightly different luminosity
at the inner boundary between the $r_{\min}=0.4R_{\ast}$ and
$r_{\min}=0.6R_{\ast}$ simulations (see Table 1) but has no measurable
difference at the outer boundaries where
$L_{\mathrm{out}}=$0.98772\text{\,}\mathrm{L}_{\sun}$$ for all simulations.
The slight difference between the inner and outer boundary drives the
background state to evolve, however this happens on the thermal time-scale
which is much longer than the duration of the simulations, as discussed below.

Table 1:  Simulation parameters: inner and outer simulation boundary, number of cells in the radial direction $N_{r}$, radial cell width in units of the pressure scale height at the outer boundary, total simulation duration $\tau_{\mathrm{tot}}$ in units of the convective turnover time, the steady-state convective turnover time $\tau_{\textrm{conv}}$, the inner luminosity in units of the solar luminosity, the average temperature at the outer boundary $T(r_{\max})$, the average density at the outer boundary $\rho(r_{\max})$, hydrodynamic boundary layer thickness $\delta_{v_{r}}$ in units of the stellar radius.  $r_{\min}/R_{\ast}$ | $r_{\max}/R_{\ast}$ | $N_{r}$ | $\Delta_{r}/H_{p}(r_{\max})$ | $\tau_{\mathrm{tot}}/\tau_{\textrm{conv}}$ | $\tau_{\textrm{conv}}/$\mathrm{1}\mathrm{0}^{5}\mathrm{s}$$ | $L_{\mathrm{in}}/L_{\sun}$ | $T(r_{\max})/$\mathrm{1}\mathrm{0}^{5}\mathrm{K}$$ | $\rho(r_{\max})/$\mathrm{1}\mathrm{0}^{-3}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$ | $\delta_{v_{r}}/R_{\ast}$ |
---|---|---|---|---|---|---|---|---|---|---
0.4 | 0.9 | 840 | 0.017 | $346$ | 5.8 | 0.98778 | 6.10 | 23.7 | 0.010 |
0.4 | 0.94 | 912 | 0.030 | $374$ | 5.4 | 0.98778 | 3.44 | 9.9 | 0.009 |
0.4 | 0.97 | 960 | 0.060 | $547$ | 4.7 | 0.98778 | 1.60 | 2.9 | 0.007 |
0.4 | 0.99 | 996 | 0.220 | $615$ | 4.3 | 0.98778 | 0.51 | 0.3 | 0.002 |
0.6 | 0.9 | 504 | 0.017 | $366$ | 5.9 | 0.98775 | 6.10 | 23.7 | 0.010 |
0.6 | 0.94 | 576 | 0.030 | $730$ | 5.4 | 0.98775 | 3.44 | 9.9 | 0.009 |
0.6 | 0.97 | 624 | 0.061 | $551$ | 4.6 | 0.98775 | 1.60 | 2.9 | 0.007 |
0.6 | 0.99 | 660 | 0.220 | $585$ | 4.2 | 0.98775 | 0.51 | 0.3 | 0.002 |

The duration of each simulation is given in Table 1. As we will see in Section
3, the typical velocity and hence global time-scale in the convection zone
depends on the radial truncation. So, in order to perform meaningful
comparisons, we measure the evolution of the simulation in terms of the
convective turnover time $\tau_{\textrm{conv}}$, which we define as

$\displaystyle\tau_{\textrm{conv}}=\int_{r_{\textrm{Schw}}}^{r_{\max}}\frac{\mathrm{d}r}{\left<v_{\textrm{rms}}\right>_{t}},$
(5)

where
$v_{\textrm{rms}}(r,t)=\sqrt{\left<\boldsymbol{v}^{2}(r,\theta,t)\right>_{\theta}}$
and $\left<\right>_{t}$ and $\left<\right>_{\theta}$ denote the time and
volume-weighted angular average, respectively. This estimates the typical time
to cross the convection zone at the r.m.s. velocity. As shown in Fig. 2, the
evolution is characterised by an initial relaxation phase (approx.
$100\,\tau_{\textrm{conv}}$) during which the spherically symmetric one-
dimensional initial conditions relax and convection develops. During the
subsequent evolution the convection is quasi-steady e. g. in terms of the mean
kinetic energy density
$\left<\textrm{E}_{\mathrm{k}}\right>_{r,\theta}=\left<\rho\boldsymbol{v}^{2}/2\right>_{r,\theta}$.
The value of $\tau_{\textrm{conv}}$ and all following results are based on the
statistics of this quasi-steady period of evolution.

Figure 2:  Time evolution of the average kinetic energy density for all
considered simulations. Solid and dashed lines correspond to simulations with
$r_{\min}=0.4R_{\ast}$ and $r_{\max}=0.6R_{\ast}$, respectively.

We collect statistics for a few hundred $\tau_{\textrm{conv}}$. This is needed
to ensure that the statistics in the convection zone are converged, since the
overshooting layer is characterised by intermittent statistics (Brummell et
al., 2002; Pratt et al., 2017; Baraffe et al., 2021). Thus, the long
simulation time allows us to determine that the simulation is in a
statistically stationary steady state.

Achieving exact thermal relaxation is a common challenge for global
hydrodynamic simulations of convection based on realistic stellar interior
structures (Meakin & Arnett, 2007; Horst et al., 2020; Higl et al., 2021). The
global thermal relaxation or Kelvin-Helmholtz timescale of the initial stellar
structure used for our simulations is given by
$\tau_{\mathrm{th}}=GM^{2}/(RL)$, and ranges between
$5\text{\times}{10}^{4}\text{\,}\mathrm{y}\mathrm{r}$ and
$2\text{\times}{10}^{6}\text{\,}\mathrm{y}\mathrm{r}$. This is computationally
unreachable without changing the underlying thermal time-scale, e. g. by
boosting the luminosity to shorten the KH time-scale. However, as demonstrated
by Baraffe et al. (2021); Le Saux et al. (2021) such boosting causes
significant changes in the convection dynamics, the overshooting layer and the
internal gravity waves in the radiative zone. Consequently, we refrain from
luminosity boosting in this study. While the simulations are not completely
thermally relaxed, the evolution is very slow and smooth over the time period
we investigate. Thus, the time-scales of the processes we follow are
significantly shorter than the simulation duration.

##  3 Results: convection zone

It has been proposed (Spruit, 1997; Brandenburg, 2016; Käpylä et al., 2017;
Anders et al., 2019) that stellar convection is primarily driven by the ‘rain’
of radiatively cooled (low entropy) perturbations which originate at the top
convective layers and buoyantly fall towards the lower convective boundary.
The perturbations traverse the bulk of the convection zone quasi-
adiabatically. It takes the boundary layer dynamics near the bottom convective
boundary to break them up and extract their entropy deficit. The flows through
the convective zone of our simulation are qualitatively consistent with this
picture, specifically in regards to the adiabatic transport through the
convective bulk.

Figure 3 shows the radial profiles of the standard deviation of the radial
velocity
$\sigma(v_{r})=\sqrt{\left<(v_{r}-\left<v_{r}\right>_{\theta})^{2}\right>_{t,\theta}}$
and the temperature fluctuations
$\sigma(T)=\sqrt{\left<(T-\left<T\right>_{\theta})^{2}\right>_{t,\theta}}$.

Figure 3:  Radial profiles of the standard deviation $\sigma$ of the radial
velocity (top) and horizontal temperature fluctuations (bottom). The
Schwarzschild boundary is marked by the vertical grey line.

As expected, both profiles are almost constant in the bulk of the convection
zone. Because the background is stratified quasi-adiabatically, the latter
indicates that the fluctuations also evolve quasi-adiabatically in the
convective bulk. This implies that the fluctuations generated near the
boundaries do not have time to thermalise with the background as they travel
through the bulk. The ILES nature of the simulations prevents an exact
quantification of the time scale of this thermalisation because of the
numerical diffusion contribution. However, an upper bound is given by the
thermal diffusion time scale $\tau_{\chi}$. Here, we compute it based on the
pressure scale-height as $\tau_{\chi}=H_{p}^{2}/k_{\mathrm{rad}}$, with
thermal diffusivity $k_{\mathrm{rad}}=\chi/(\rho c_{p})$ and specific heat at
constant pressure $c_{p}$. Using this definition the thermal diffusion time
scale is of the order of
$$2\text{\times}{10}^{12}\text{\,}$-$1\text{\times}{10}^{13}\text{\,}\mathrm{s}$$
in the convection zone increasing with increasing radius, i. e. $7$ to $8$
orders magnitude larger than $\tau_{\mathrm{conv}}$ (see Table 1). Note that
this is also consistent with previous results (Viallet et al., 2015; Pratt et
al., 2016).

The data show that increasing $r_{\max}$ increases both $\sigma(v_{r})$ and
$\sigma(T)$ (note the logarithmic axis in Fig. 3). At the same time shifting
$r_{\min}$ by $20\%$ of the stellar radius (from $0.4R_{\ast}$ to
$0.6R_{\ast}$) has an almost imperceptible effect on both diagnostics. The
same trend is observed for the average kinetic energy density (see Fig. 2)
which is dominated by the convection zone. This is consistent with the results
of Hurlburt et al. (1984) who studied a convection layer under similar
boundary conditions but with an idealised setup (polytropic stratification,
ideal gas EoS, Cartesian geometry, and maximum density contrast of $21$). They
reported a growth of the mean kinetic energy with the number of pressure scale
heights in the simulation domain.

At first glance the dependence on $r_{\max}$ and the insensitivity to
$r_{\min}$ that we report may be surprising because all simulations are driven
with the same luminosity and share the same stratification and boundary
conditions. Indeed, the main energy balance in steady state in the convection
zone is between the kinetic energy and convective heat fluxes, which must add
up to transport the total stellar luminosity. However, due to mass
conservation and the density stratification, the mean kinetic energy flux is
negative. Consequently, the energy balance is not sufficient to determine the
amplitudes of the two fluxes. Thus, the amplitudes of the velocity and
temperature fluctuations are allowed to vary with $r_{\max}$ while maintaining
the same luminosity. The dependence can be explained by the properties of the
background at the top convective boundary, which change with $r_{\max}$. For
instance the background temperature drops from
$T(r_{\max}=0.9R_{\ast})=$6.1\text{\times}{10}^{5}\text{\,}\mathrm{K}$$ to
$T(r_{\max}=0.99R_{\ast})=$0.51\text{\times}{10}^{5}\text{\,}\mathrm{K}$$, and
the density drops from
$\rho(r_{\max}=0.9R_{\ast})=$23.7\text{\times}{10}^{-3}\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$
to
$\rho(r_{\max}=0.99R_{\ast})=$0.3\text{\times}{10}^{-3}\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$,
see Table 1. At the same time, the properties of the background at the lower
convective boundary are largely the same across all simulations (e. g.
$T(r_{\textrm{Schw}})=$2.15\text{\times}{10}^{6}\text{\,}\mathrm{K}$$ and
$\rho(r_{\textrm{Schw}})=$0.16\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$).
This implies that for a fixed realistic stellar background, the properties of
the radial velocity and temperature fluctuations are determined by the
conditions in the boundary layers, at least for the type of convection zone
that the simulations can represent. This leads us to consider the outer
boundary layer in more detail. Specifically, we consider the fluctuations at
$r_{\max}$ and how they evolve across the boundary layers.

As we evolve the specific internal energy density as a primitive variable,
$\sigma(T)(r_{\max})$ is implicitly set by the stratification and the imposed
luminosity. The larger value of $\sigma(T)(r_{\max})$ for increasing
$r_{\max}$ can be primarily attributed to the corresponding decrease of the
mean heat capacity $\left<\rho c_{p}\right>_{t,\theta}$. The luminosity is
carried predominantly444Because the radial velocity and the temperature are
computed on staggered grids in MUSIC, the radial velocity and hence the
convective heat flux interpolated to the centre of the outermost $r_{\max}$
cell can be non-zero. At this location, the radiative diffusivity is
negligible and hence the radiative flux is insignificant. by the convective
heat flux in this region. Hence, as the mean heat capacity at $r_{\max}$
decreases across simulations, $\sigma(T)(r_{\max})$ must increase to maintain
the same luminosity.

To describe the radial structure of the temperature and radial velocity
fluctuations at the outer boundary, we use the linear rescaling

$\displaystyle\hat{\sigma}=\frac{\sigma(r)-\sigma_{\mathrm{bulk}}}{|\sigma(r_{\max})-\sigma_{\mathrm{bulk}}|},$
(6)

where $\sigma$ stands in for $\sigma(T)$ and $\sigma(v_{r})$ and
$\sigma_{\mathrm{bulk}}$ is the volume-weighted mean value of the respective
$\sigma(r)$ in the convective bulk. Here, we consider the convective ‘bulk’ to
consist of the shell $[r_{\max}-0.08R_{\ast},r_{\max}-0.02R_{\ast}]$.

Figure 4 shows $\hat{\sigma}(T)$ and $\hat{\sigma}(v_{r})$ along with the
associated boundary layer thickness defined similarly to Featherstone &
Hindman (2016) as

$\displaystyle\delta=\int_{r_{\max}-0.08R_{\ast}}^{r_{\max}}\hat{\sigma}(r)dr.$
(7)

Here we set the lower integration boundary well above the lower convective
boundary, so that $\delta$ is not influenced by the boundary layer near
$r_{\textrm{Schw}}$. As the figure shows, the hydrodynamic boundary layer
becomes significantly thinner with increasing $r_{\max}$, as indicated by
$\delta_{v_{r}}$, see also Table 1. In contrast, the rescaled thermal
fluctuations $\hat{\sigma}(T)$ have the same radial dependence across all
simulations with a thermal boundary layer thickness
$\delta_{T}=0.002R_{\ast}$. As a result, the ratio $\delta_{v_{r}}/\delta_{T}$
decreases with increasing $r_{\max}$ reaching unity for the
$r_{\max}=0.99R_{\ast}$ simulations. This can explain the $r_{\max}$
dependence of $\sigma(v_{r})$ and $\sigma(T)$ in the convective bulk, as
follows.

By definition, the hydrodynamic boundary layer is characterised by the growth
of the radial velocity from zero to its bulk value, due to e. g. buoyancy.
Similarly, the thermal boundary layer is characterised by the decay of
temperature fluctuations through e. g. diffusion and horizontal shear. Thus,
as $\delta_{v_{r}}/\delta_{T}$ increases the hydrodynamic boundary layer is
dominated by smaller temperature (and hence density) perturbations and thus
weaker buoyancy. At the same time, the thermal boundary layer is characterised
by weaker radial advection and hence a longer time scale over which the
temperature fluctuations can decay. Hence, as $\delta_{v_{r}}/\delta_{T}$
increases the saturation values of both the temperature and velocity
fluctuations, i. e. their values in the convective bulk, drop. The precise
values of $\delta_{v_{r}}$ are $\delta_{T}$ are a result of a detailed balance
between radial advection and buoyancy, on one hand, and diffusion, horizontal
shear and dissipation on the other. As is typical for global stellar
hydrodynamic simulations, the ILES nature of the simulations precludes an in-
depth discussion of this balance. However, we note that the boundary layers
thicknesses are not set by the local stratification, e. g. they do no scale
with $H_{p}(r_{\max})$. Moreover, the measured value of $\delta_{T}$ is
limited by the grid resolution in all simulations (corresponding to
approximately 3 grid cells), while $\delta_{v_{r}}$ is only resolution-limited
in the largest $r_{\max}$ simulations. Thus, the ratio
$\delta_{v_{r}}/\delta_{T}$ should increase with increasing radial resolution.
We establish this by repeating the $0.4R_{\ast}-0.94R_{\ast}$ and
$0.6R_{\ast}-0.99R_{\ast}$ simulations with $\sim$50\text{\,}\mathrm{\char
37\relax}$$ larger $N_{r}$, while keeping all other parameters the same. The
simulations start from the same initial conditions and evolve for
$1.8\text{\times}{10}^{8}\text{\,}\mathrm{s}$. The thermal boundary layer
shrinks in physical units with the increased resolution, so that $\delta_{T}$
remains approximately the same number of radial grid cells. However, the
$\delta_{v_{r}}$ remains approximately the same thickness (in physical units)
for the $0.4R_{\ast}-0.94R_{\ast}$ truncation and it grows by approximately
$20\text{\,}\mathrm{\char 37\relax}$ for the $0.6R_{\ast}-0.99R_{\ast}$
truncation. Importantly, the radial profiles of $\sigma(T)$ and
$\sigma(v_{r})$ in the convective bulk are not affected by the small increase
in resolution.

Figure 4:  Radial profiles of temperature fluctuations (top) and radial
velocity (bottom) in the upper boundary layer linearly rescaled as indicated
by Eq. 7. The vertical lines indicate the corresponding boundary layer
thickness, see Eq. 7.

In summary, the thermal perturbations originating in the upper boundary layer
are effectively frozen while they travel through the convective bulk, because
of the long diffusion time scale. As a result the speed and temperature
contrast of the radial flows near the lower boundary layer is effectively set
by the conditions at the upper boundary layers and not in the convective bulk.
The key parameters are the density stratification and the ratio between the
hydrodynamic and thermal boundary layer thicknesses. For this set of
simulations, the density stratification is steep enough to lead to decreasing
heat capacity with height, which requires larger temperature fluctuations to
maintain the same luminosity. In addition, the ratio of
$\delta_{v_{r}}/\delta_{T}$ decreases towards unity with increasing
$r_{\max}$, meaning that an increasingly bigger part of the radial velocity is
generated by not-fully-decayed thermal fluctuations with larger buoyancy. This
leads to $\sigma(T)_{\mathrm{bulk}}$ which scales with $\sigma(T)(r_{\max})$
and $\sigma(v_{r})_{\mathrm{bulk}}$ which increases correspondingly. This
qualitative picture illustrates how the boundary layer dynamics can determine
the properties of the convective bulk and illustrates the non-locality of
boundary-driven convection. It is also consistent with the results of Cossette
& Rast (2016), who find that for boundary-driven convection the initial
density and entropy contrast of the adiabatically descending perturbations
determine the dynamical properties of the convective bulk.

##  4 Overshooting layer

###  4.1 Characteristic overshooting depth

As mentioned, the background properties near the lower convective boundary are
the same among all simulations. However, as we show in this section, the non-
locality of convection means that the $r_{\max}$-dependence of the outer
convective layers makes a qualitative difference to the lower convective
boundary layer dynamics and convective overshooting in particular.

We recall that we use the Schwarzschild criterion to locate the lower
convection boundary. The overshooting layer is then the region below this
boundary where convection flows are able to penetrate. Specifically, we
characterise the overshooting in the same way as Baraffe et al. (2021).
Namely, we consider the radial kinetic energy $F^{\mathrm{ke}}_{r}$ and
convective heat flux $F^{\mathrm{conv}}_{r}$, defined by

$\displaystyle\boldsymbol{F}^{\mathrm{conv}}_{r}(r,\theta,t)$
$\displaystyle=h^{\prime}(\rho v_{r})^{\prime},$ (8)
$\displaystyle\boldsymbol{F}^{\mathrm{ke}}_{r}(r,\theta,t)$
$\displaystyle=v_{r}\textrm{E}_{\mathrm{k}}=v_{r}\frac{1}{2}\rho\boldsymbol{v}^{2},$
(9)

where $h=e+p/\rho$ is the specific enthalpy and the primes denote fluctuations
from the angular average, e. g. $h^{\prime}=h-\left<h\right>_{\theta}$. Here,
the convective heat flux contains only the enthalpy fluctuations. From the
hydrostatic energy balance, i. e. hydrostatic equilibrium, it can be seen that
the mean enthalpy balances the gravitational potential (see Appendix A). We
consider $F^{\mathrm{conv}}$ because our simulations solve the internal energy
equations where the enthalpy appears naturally. A quantity that is also often
used in the literature is $\left<F_{\Delta T}\right>_{t,\theta}=\left<\rho
c_{p}\Delta Tv_{r}\right>_{t,\theta}$ (Hurlburt et al., 1984; Zahn, 1991;
Baraffe et al., 2021), where $\Delta T$ is the temperature fluctuations with
respect to a reference model. It is closely related to
$\left<F^{\mathrm{conv}}_{r}\right>_{t,\theta}$ (as can be seen analytically
e.g. for an ideal gas), and in the considered simulations has the same
quantitative behaviour.

To identify an instantaneous local overshooting event, we track the locations
of the first zero-crossing of the two fluxes $F^{\mathrm{conv}}_{r}$ and
$F^{\mathrm{ke}}_{r}$ below the Schwarzschild boundary,
$l^{\mathrm{ke}}(\theta,t)$ and $l^{\mathrm{conv}}(\theta,t)$ respectively.
Specifically, we consider the instantaneous position of the boundary
$r_{\textrm{Schw}}(t)=\min\\{r:\left<\nabla-\nabla_{\mathrm{ad}}\right>_{\theta}\geq
0\\}$. Retaining the time-dependence of $r_{\textrm{Schw}}(t)$ here allows us
to extract the effects of the convection motion in the stably-stratified
region and exclude the motion of the convective boundary itself. Brummell et
al. (2002) distinguishes the two as ‘convective overshooting’ and ‘convective
penetration’, respectively. Note, however that the time variation of the
boundary $r_{\textrm{Schw}}(t)$ is very slow (significantly slower than the
$\tau_{\mathrm{conv}}$ and much slower than the time scale of individual
overshooting flows) and as such, it would be perceived as stationary by
individual penetrating flows. We discuss the modification of the background
and possible convective penetration in more detail in Section 4.2.

The distribution of both $l^{\mathrm{ke}}$ and $l^{\mathrm{conv}}$ can be
associated with significant frequency of extreme events both in time and space
(Pratt et al., 2017). We confirm this in the presented datasets, so we use two
summary statistics to describe the overshooting

$\displaystyle l_{\mathrm{bulk}}(t)$
$\displaystyle=\left<l(\theta,t)\right>_{\theta},$ (10) $\displaystyle
l_{\max}(t)$ $\displaystyle=\max_{\theta}\\{l(\theta,t)\\}.$ (11)

Here $l_{\mathrm{bulk}}$ describes the angular average and is interpreted as
the _typical_ depth of overshooting, while $l_{\max}$ is the angular maximum,
i. e. the instantaneous maximum depth of overshooting and $l(\theta,t)$ can
refer to either $l^{\mathrm{conv}}(\theta,t)$ or $l^{\mathrm{ke}}(\theta,t)$.

Figure 5:  Time distribution of the instantaneous overshooting depth in units
of the pressure scale height at the convective boundary (left axis) and radial
cell width (right axis): mean-in-$\theta$ $l_{\mathrm{bulk}}$ (top row) and
max-in-$\theta$ $l_{\max}$ (bottom row). The left hand side panels use the
convective heat flux and the right column – the kinetic energy flux, i. e. the
radial velocity. The simulations are ordered by decreasing $r_{\max}$ and
$r_{\min}$.

Figure 5 shows the distribution of these quantities across time in the
different simulations in units of the pressure scale height at the
Schwarzschild boundary $H_{p}(r_{\textrm{Schw}})$. We exclude a $None$ angle
around the polar axes, to reduce the impact of the polar boundaries – the
reflective boundary conditions drive strong radial flows which result in
locally deeper overshooting. The figure shows that there is a clear dependence
of the overshooting depth on $r_{\max}$ – larger $r_{\max}$ is associated with
a thicker overshooting layer. There is also such a dependence for both
$l_{\mathrm{bulk}}$ and $l_{\max}$, based on either of the fluxes, independent
of whether we consider the time-average or either of the tails of the temporal
distribution. At the same time there does not seem to be a clear dependence on
$r_{\min}$ shared by all overshooting depth diagnostics, for the two
$r_{\min}$ we examine.

These observations correlate well with the dependence of the convective
intensity (as measured by e. g. $\sigma(v_{r})$ and $\sigma(T)$, i. e. the
radial velocity and temperature fluctuations) on $r_{\max}$ and $r_{\min}$,
see Fig. 3. This suggests that the larger overshooting depth is caused by more
vigorous convection, i. e. faster and colder downflows can reach deeper into
the stably-stratified region and such flows are generated in simulations with
larger $r_{\max}$. For instance, in the $r_{\min}=0.4R_{\ast}$ simulations
$\sigma(v_{r})(r_{\textrm{Schw}})$ takes the values
$[$5\text{\times}{10}^{2}\text{\,}$,$1\text{\times}{10}^{3}\text{\,}$,$3\text{\times}{10}^{3}\text{\,}$,$8\text{\times}{10}^{3}\text{\,}$]$
$\mathrm{c}\mathrm{m}\,\mathrm{s}^{-1}$ for
$r_{\max}/R_{\ast}=[0.9,0.94,0.97,0.99]$, respectively. The lack of dependence
on $r_{\min}$ indicates that either the feedback from the radiative zone is a
negligible influence on the overshooting, or that the feedback itself does not
depend very strongly on the spectral properties of internal gravity waves, at
least up to the variations of the IGWs across the presented simulations. This
will be discussed in more detail in Section 5, where we show that the spectrum
of the travelling waves does not depend qualitatively on $r_{\min}$.

It follows trivially from the definitions Equations 8 and 9 that the loci of
$l^{\mathrm{conv}}(\theta,t)$ and $l^{\mathrm{ke}}(\theta,t)$ can be
identified with the first zero-crossings of the vertical momentum and enthalpy
fluctuations. In other words, they track the bottom end of the downflows as
defined by their velocity and thermal properties. Furthermore, by construction
$l^{\mathrm{ke}}\geq l^{\mathrm{conv}}$, because $l^{\mathrm{ke}}$ depends
only on the zeroes of the radial velocity, while $l^{\mathrm{conv}}$ is also
triggered by the zeroes of the enthalpy fluctuations. Physically, this is
because the radial buoyant acceleration vanishes at
$l^{\mathrm{conv}}(\theta,t)$, while the radial velocity vanishes at
$l^{\mathrm{ke}}(\theta,t)$. The inequality is not saturated in the presented
simulations despite the long run times. Even instantaneously, the difference
between the bulk values
$l^{\mathrm{ke}}_{\mathrm{bulk}}(t)-l^{\mathrm{conv}}_{\mathrm{bulk}}(t)$
remains steady at about $8\text{\,}\mathrm{\char 37\relax}$ to
$10\text{\,}\mathrm{\char 37\relax}$ of their mean value
$(l^{\mathrm{ke}}_{\mathrm{bulk}}(t)+l^{\mathrm{conv}}_{\mathrm{bulk}}(t))/2$.
The difference between the extreme values
$l^{\mathrm{ke}}_{\max}(t)-l^{\mathrm{conv}}_{\max}(t)$ is approximately
$20\text{\,}\mathrm{\char 37\relax}$ to $30\text{\,}\mathrm{\char 37\relax}$
of the corresponding mean.

This indicates that typically the temperature fluctuations dissipate before
the downflows stop descending (whereby we consider the statistical variation
to be negligible). While buoyancy is sufficient to explain this, it is not the
only the process acting on the overshooting flows. More generally for the
inequality $l^{\mathrm{ke}}>l^{\mathrm{conv}}$ to be strict, the mechanisms
removing the heat deficit of downflows must operate on shorter times-scales
than those dissipating the excess radial momentum. Broadly speaking, the heat
deficit can be removed by adiabatic compression and heat exchange with the
background and the IGWs, e. g. enthalpy mixing/entrainment and radiative flux,
wave excitation and breaking. The radial momentum excess is removed by
buoyancy, turbulent mixing, viscous dissipation, and IGW excitation.

###  4.2 Modification of the thermal background

As the downflows travel below their local, instantaneous
$r_{\textrm{Schw}}(\theta,t)$ depth, and the background stratification becomes
sub-adiabatic, they become hotter than the background. Some of this excess
heat may be deposited in the overshooting layer through irreversible processes
(e.g. IGW excitation/breaking, mixing, etc.). Indeed, as Fig. 6 shows, there
is a good correlation between the overshooting depth and a layer which is
slightly hotter than the background at the beginning of the quasi-steady
state.

Figure 6:  (top) Relative difference $\Delta T$ in the mean temperature
background in the quasi-steady state $\left<T\right>_{\theta,t}$ from the
background at the beginning of the quasi-steady state
$\left<T_{0}\right>_{\theta}$. The locations of the overshooting depth
measured by the convective heat flux are marked for reference. (bottom) Radial
profile of the squared rate-of-strain averaged over $\theta$ and time. The
location of the peak $\Delta T$ is marked with $\bullet$ as a visual guide.
The vertical grey line in both panels corresponds to the Schwarzschild
boundary.

Indications for a similar process (a layer with a negative enthalpy flux, an
excess radiative flux, or slightly larger mean temperature) has been
previously reported by e. g. Rogers & Glatzmaier (2006); Brun et al. (2011);
Brun et al. (2017); Korre et al. (2019); Baraffe et al. (2021) and Higl et al.
(2021) in a variety of parameter regimes (see Baraffe et al. 2021, section 3.2
for more details).

Baraffe et al. (2021) note that there is a good correlation between the
location of the radial profile of the heating layer and the trace of the
square rate-of-strain
$\mathrm{Tr}S^{2}=\mathrm{Tr}(\nabla\boldsymbol{v}+\nabla\boldsymbol{v}^{T})^{2}/4$,
where $\nabla\boldsymbol{v}$ denotes the velocity gradient. We can confirm
this correlation as well, see Fig. 6. Note that $\mathrm{Tr}S^{2}$ encodes two
types of material deformation – dilatation and shear. However, because of the
low Mach number $\la 10^{-3}$, $\mathrm{Tr}S^{2}$ is dominated by the
incompressible shear component (not shown). This suggests that the primary
component of the heat deposition is shear-induced mixing of previously
adiabatically compressed hot material.

This local heating causes a dip in the sub-adiabaticity
$\nabla-\nabla_{\mathrm{ad}}$ just below the convective boundary and a
corresponding increase below the dip, around the heat bump. If left unstopped
this process may lead to the so-called convective penetration in the language
of Brummell et al. (2002), i. e. the descent of the convective boundary and
growth of the convection zone.

Energetically, in a sub-adiabatically stratified medium adiabatic compression
of downflows and adiabatic expansion of upflows both lead to downward/negative
enthalpy flux (Hurlburt et al., 1994; Muthsam et al., 1995; Brun et al., 2011;
Pratt et al., 2017; Korre et al., 2019; Käpylä et al., 2019). So, as discussed
in Baraffe et al. (2021), the reason for the local change in the background
stratification must be that the heat deposited by the enthalpy flux cannot be
efficiently evacuated. Considering the internal energy equation, Eq. 3, the
excess heat deposited by the enthalpy flux can only be evacuated through
thermal diffusion, i. e. radiative flux (considering that the local Mach
number is $\la 10^{-3}$ the $p\nabla\cdot\boldsymbol{v}$ term does not play a
significant role). Indeed, as noted by Baraffe et al. (2021), the local
heating must lead to increased thermal diffusivity and hence increased
radiative flux. The feedback should eventually grow sufficiently to balance
the negative enthalpy flux. Indeed, we note a slow growth in the heating rate
throughout the simulated time interval555 The available data is consistent
with a constant convective heat flux and a linearly growing radiative flux in
the overshooting layer. Extrapolating this in a zero-order approximation, the
two fluxes would match and the excess heating would saturate on a time scale
of $O($\mathrm{1}\mathrm{0}^{9}\mathrm{s}$)$. . However, such a regime cannot
be reached with the presented simulations.

##  5 Radiative zone

One of the effects of the convection zone and the overshooting layer is to
excite internal gravity waves (IGW) in the radiative zone. Internal gravity
waves play a crucial role in the internal structure and evolution of stars, as
they redistribute angular momentum, energy and chemical species, deposited in
the overshooting layer by convective flows. It is interesting then to consider
how the changes in the overshooting layer and the convection zone due to the
different radial truncations are reflected in the IGWs, as well as the
strength of any potential feedback onto the convective overshooting.

We focus the analysis on the power spectra of radial velocity. The spectral
coefficients $\widehat{v_{r}}(\ell,\omega)$ are obtained using the conventions
in Le Saux et al. (2021). In short, we perform a Blackman-windowed Fourier
transform with respect to time and a 2D spherical harmonic transform ($m=0$)
with respect to the angular direction $\theta$. The resulting 2D power
spectrum is then given by

$\displaystyle\mathcal{P}[v_{r}](\ell,\omega)=\begin{cases}\left|\widehat{v_{r}}(\ell,\omega)\right|^{2}&\textrm{for
}\omega=0\\\ 2\left|\widehat{v_{r}}(\ell,\omega)\right|^{2}&\textrm{for
}\omega>0,\end{cases}$ (12)

where $\ell$ is the angular degree and $\omega$ is the frequency.

Figure 7 shows the radial velocity power spectra for a representative
selection of simulations at three locations. The spectra are computed from a
subset of the data spanning $2.25\text{\times}{10}^{7}\text{\,}\mathrm{s}$ in
the quasi-steady state. The deepest location is chosen in the bulk of the
radiative zone at $r=0.64R_{\ast}=r_{\textrm{Schw}}-H_{p}(r_{\textrm{Schw}})$
far away from any overshooting activity and we expect it to be dominated by
the IGW signal. The intermediate location is at the convective boundary at
$r_{\textrm{Schw}}$ and as such the spectra there are expected to exhibit
signs of the overshooting dynamics superposed with the IGW signal. The highest
location is in the bulk of the convective zone at $r=0.84R_{\ast}$ and is
shown as a reference to illustrate the length- and time-scales of the
convective flow which is responsible for exciting the IGW signal.

Figure 7:  Power spectra of the radial velocity deep in the radiative zone –
$H_{p}(r_{\textrm{Schw}})$ below the Schwarzschild boundary (left), at the
bottom end of the mixing zone at $l_{\max}$ (middle), and in the bulk of the
convection zone $0.841R_{\ast}$ (right). Each row corresponds to spectra from
the same simulation. All spectra are produced from data spanning the last
$2.25\text{\times}{10}^{7}\text{\,}\mathrm{s}$ of each simulation. The lower
five percentile of the data are rendered in black to improve the colour
contrast and aid visual comparison. The frequency axis extends up to $3/4$ of
the Nyquist frequency $\mathrm{Ny}$ to avoid aliasing artefacts.
$\mathrm{Ny}=(2\delta
t)^{-1}=$100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$ based
on separation between snapshots of $\delta
t\sim$5\text{\times}{10}^{3}\text{\,}\mathrm{s}$$. Note that the power deficit
at the low-$\omega$ high-$\ell$ region of the $r=0.64R_{\ast}$ panels is
associated with resolution constraints related to the finite duration of the
analysed data and the cell size.

The classical pattern of IGWs (Alvan et al., 2014; Horst et al., 2020) is
clearly visible in the radiative zone for all simulations, most prominently in
the low-$\ell$ range. This pattern can be inferred from the dispersion
relation ${\omega^{2}}/{N^{2}}={k_{h}^{2}}/{\boldsymbol{k}^{2}}$ (see e.g.
Press (1981)), where $\boldsymbol{k}$ is the wave vector and
$k_{h}^{2}=\ell\left(\ell+1\right)r^{-2}$ its horizontal component. After re-
arranging, the dispersion relation can be seen as

$\displaystyle\left(\left(\frac{2\ell+1}{2rk_{r}}\right)^{2}+1-\frac{1}{4r^{2}k_{r}^{2}}\right)\left(1-\frac{\omega^{2}}{N^{2}}\right)=1,$
(13)

which describes a family of hyperbolae in $(2\ell+1)^{2}$ and $\omega^{2}$ for
a fixed $r$ (parametrised by $k_{r}$).

To confirm quantitatively that the computed spectra in the radiative zone
contain IGWs, we compute the frequencies of the $g$-modes (i. e. the standing
waves) associated with the radiative zones of the simulation for a few
distinct values of $\ell$. The computation is performed using the stellar
oscillations code GYRE666Version 6.0, see https://gyre.readthedocs.io
(Townsend & Teitler, 2013; Townsend et al., 2018). Since the modes depend only
on the stratification and the geometry of the resonant cavity (i. e. the
radiative zone), we obtain two distinct sets of modes for the
$r_{\min}=0.4R_{\ast}$ and $r_{\min}=0.6R_{\ast}$ sets of simulations. As Fig.
8 shows for the $\ell=2$ modes, the frequencies obtained with GYRE match
almost perfectly the peaks in the spectra for both values of $r_{\min}$ in the
radiative zone and at the convective boundary locations.

Figure 8:  Frequency distribution of the radial velocity power spectrum
$|\widehat{v_{r}}|^{2}(l=2,\omega)$ deep in the radiative zone (left) and at
the convective boundary $r_{\textrm{Schw}}$ (right). Top row shows a
comparison of simulations with $r_{\min}=0.4R_{\ast}$ and the bottom row –
with $r_{\max}=0.97R_{\ast}$. The dashed and dotted vertical lines show the
frequencies of the $g$-modes obtained with GYRE for the $r_{\min}=0.4R_{\ast}$
and $r_{\min}=0.6R_{\ast}$ respectively. The solid coloured vertical lines
show the convective frequency $\omega_{\mathrm{conv}}=1/\tau_{\mathrm{conv}}$
for the respective simulation.

###  5.1 Dependence on $r_{\max}$ and $r_{\min}$

The value of $r_{\min}$ changes the depth of the resonant cavity where IGWs
can propagate – for $r_{\min}=0.6R_{\ast}$ it is very shallow, spanning only
about $0.12R_{\ast}\sim 1.4H_{p}(r_{\textrm{Schw}})$, while for
$r_{\min}=0.4R_{\ast}$ it spans approximately $0.32R_{\ast}\sim
3.9H_{p}(r_{\textrm{Schw}})$. As a result the frequency distribution of the
$g$-modes changes significantly. For instance, for the $r_{\min}=0.6R_{\ast}$
simulations there are fewer $g$ modes at high frequencies than in the
$r_{\min}=0.4R_{\ast}$ simulations, as shown in Fig. 8. Apart from this
however, the power spectrum is quite similar between the two sets of
simulations, particularly for
$\omega\ga$5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{H}\mathrm{z}$$. This
is consistent with the earlier observations that the value of $r_{\min}$ does
not seem to have any significant effect on the depth of the overshooting layer
and the convection zone, for the presented diagnostics at least (see Section
4.1). This suggest that the detailed frequency distribution of the $g$-modes
has no direct impact on the overshooting dynamics.

On the other hand, the value of $r_{\max}$ has no impact on the frequency
distribution of the $g$-modes. However, the amplitude of both travelling and
standing waves in the radiative bulk grows significantly with increasing
$r_{\max}$, as can be seen in both Figs. 7 and 8. The qualitative relation can
be intuitively expected because of the more vigorous convection in the higher
$r_{\max}$ simulations which drives faster overshooting flows and a deeper
overshooting layer. Figure 7 shows that the peak of the power spectrum in the
radiative zone grows by approximately 3 orders of magnitude from the
$r_{\max}=0.9R_{\ast}$ to the $r_{\max}=0.99R_{\ast}$ simulation, from
$2.4\text{\,}$ to
$4.4\text{\times}{10}^{3}\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{s}^{-2}$
respectively. At the same time, in the bulk of the convection zone the peak of
the power spectrum grows by less than an order of magnitude (from
$5.8\text{\times}{10}^{6}\text{\,}$ to
$3.0\text{\times}{10}^{7}\text{\,}\mathrm{c}\mathrm{m}^{2}\,\mathrm{s}^{-2}$)
and so does $\sigma(v_{r})$ (see Fig. 3). We note that $\sigma(v_{r})$ when
restricted only to the downflows in the convective bulk also does not have as
strong $r_{\max}$ dependence as the IGWs amplitude (not shown).

With increasing $r_{\max}$ the energy of the IGWs is shifted consistently to
higher frequencies (cf. Fig. 7). This is consistent with theoretical
expectations. Lecoanet & Quataert (2013) link the frequency of the IGWs
excited by Reynolds stress at the convective boundary with the characteristic
convective frequency $\omega_{\mathrm{conv}}=1/\tau_{\mathrm{conv}}$, which
grows with $r_{\max}$ (see Table 1). Similarly, Pinçon et al. (2016) link the
frequency of the IGWs excited by penetrating plumes with the plume lifetime in
the penetration zone. While it is difficult to obtain a direct estimate of the
lifetime of penetration events with Eulerian statistics, Pinçon et al. (2016)
suggests the convective time scale $\tau_{\mathrm{conv}}$ as a good first
approximation. Under this approximation, the plume lifetime would decrease
with $r_{\max}$. 777 Alternatively, a simple ballistic approximation can be
constructed from the average downward velocity at the Schwarzschild boundary,
$\left<v_{\downarrow}(r_{\textrm{Schw}})\right>_{\theta,t}$ and either of the
penetration depths $l_{\mathrm{bulk}}^{\mathrm{ke}}$ and
$l_{\max}^{\mathrm{ke}}$. For the presented simulations, both
$l_{\mathrm{bulk}}^{\mathrm{ke}}/\left<v_{\downarrow}(r_{\textrm{Schw}})\right>_{\theta,t}$
and
$l_{\max}^{\mathrm{ke}}/\left<v_{\downarrow}(r_{\textrm{Schw}})\right>_{\theta,t}$
decrease as $r_{\max}$ increases, supporting the notion that the plume
lifetime decreases as $r_{\max}$ increases. Thus, both excitation mechanisms
indicate that with increasing $r_{\max}$ more energy should be deposited in
higher frequency IGWs.

In terms of length scales, at a given frequency, with increasing $r_{\max}$,
the energy is shifted to larger length scales (smaller $\ell$). This is also
consistent with theory, since increasing $r_{\max}$ leads to convective eddies
with larger length scale (smaller $\ell_{\textrm{eddy}}$). In that case, they
are expected to excite waves with $\ell\leq\ell_{\textrm{eddy}}$ (Lecoanet &
Quataert, 2013).

Overall, the dependence of the IGW power spectrum in the radiative zone on
$r_{\max}$ are qualitatively analogous to the effect of luminosity boosting
discussed by Le Saux et al. (2021), who also note an increase in power and
shift to higher frequencies and larger length scales with increasing
luminosity.

##  6 Summary and discussion

In this study, we investigate the effects of the radial extent on the dynamics
of hydrodynamic 2D simulations of a solar-like stellar model. We consider the
dynamics in the convective, overshooting and radiative zones.

The location of the outer boundary plays a crucial role in determining the
convective intensity as measured by the radial velocity and temperature
fluctuations, i. e. $\sigma(v_{r})$ and $\sigma(T)$. As $r_{\max}$ is
increased, the steep hydrostatic density stratification leads to a decrease of
the heat capacity of the outer layers and a corresponding increase in the
temperature fluctuations. The latter is required, so that the convective heat
flux can transport the given fixed luminosity. As we do not consider
$r_{\max}$ high enough to include layers with inefficient convection, the
radiative flux is negligible throughout the simulated convection zone. The
larger temperature fluctuations in the outer boundary layers drive faster
radial flows through buoyancy. In turn the faster radial flows advect larger
temperature fluctuations out of the diffusive thermal boundary layer. These
fluctuations then travel quasi-adiabatically in the convective bulk before
dissipating at the lower convective boundary layer and driving convective
boundary mixing.

The more intense convection leads to a deeper overshooting layer below the
lower convective boundary. We measure this by the radial kinetic and
convective heat fluxes. To characterise both the typical and extreme
overshooting events we consider both the mean and the extreme overshooting
depth, following Baraffe et al. (2021). Both diagnostics grow rapidly with
increasing $r_{\max}$.

Over several hundred $\tau_{\textrm{conv}}$ the overshooting leads to a small
but measurable excess heating of the background in the overshooting layer. The
heating is well-correlated with strong mean background shear, indicating that
the process responsible for it is likely mechanical mixing of adiabatically
compressed hot convective material. For a thermally relaxed state to be
reached the excess heating must be balanced by a feedback mechanism, e. g.
naturally with the increased background temperature comes an increase of the
radiative flux. Despite the long simulation times (between $\sim 340$ and
$730\,\tau_{\textrm{conv}}$) a thermally relaxed state could not be reached in
the simulations, indicating that the relaxation time of the excess heating is
significantly longer (at least $O(10^{3}\,\tau_{\textrm{conv}})$).

The location of the inner boundary of the simulations $r_{\min}$ is expected
_a priori_ to change the properties of the internal gravity waves in the
radiative zone. In particular, it affects the frequency of the standing waves
(the $g$-modes) and some wave reflection off the bottom boundary is observed
for shallower radiative zones. However, this appears to have no strong or
consistent impact on the power spectrum of the travelling waves, the
convection zone and overshooting layers. This insensitivity is expected to
hold while the feedback of the IGWs on the overshooting layer is weak, e. g.
the radiative zone is deep enough to allow for significant wave damping and no
direct interaction between the overshooting and bottom boundary layers.

The position of the outer convective boundary impacts significantly the
internal gravity waves in the radiative zone. Even though the frequency
dependence of the $g$-modes remains unaffected, the energy of all IGWs is
increased because of larger overshooting velocities (which also imply a larger
overshooting depth). In particular, larger-frequency waves are more strongly
excited in agreement with theoretical expectations (Lecoanet & Quataert,
2013).

###  6.1 Luminosity boosting analogy

It is interesting to note that the dependence of the overshooting layer on
$r_{\max}$ is qualitatively similar to that on luminosity boosting discussed
by Baraffe et al. (2021). Luminosity boosting is an often-used technique in
stellar hydrodynamics, which shortens the thermal time-scale of the model and
brings it closer to the convective one. This leads to a substantial reduction
in computational costs and makes reaching thermal equilibrium feasible given a
large enough boost factor. However, it has to be done carefully in order to
maintain the original background stratification of the un-boosted model, e. g.
the thermal diffusivity has to be increased proportionally and stratification
in the convection zone has to be adjusted to as close to adiabatic as
numerically feasible. Even then, as discussed in Baraffe et al. (2021), it
comes at the cost of affecting the overshooting dynamics by increasing the
overshooting depth, the local heating in the overshooting layer and the shape
of the IGW spectrum (Lecoanet et al., 2019).

This effect is qualitatively similar to the dependence we find on $r_{\max}$,
especially considering that the stellar models used by Baraffe et al. (2021)
and in this study have significantly different background stratification and
luminosity. Recalling that all the presented simulations have the same
luminosity, this highlights that the underlying dependency is not on the
energy flux _per se_ , but on the convective intensity, in the form of e.g.
the temperature contrast and the radial velocity. Ultimately, because of the
adiabatic evolution in the convective bulk, the dependence is on the
stratification at the outer boundary, which determines the fluctuations
required to transport the luminosity there.

It is difficult to make a direct quantitive comparison between the effects of
increasing $r_{\max}$ in the simulations presented here and the effects of
boosting the luminosity discussed by Baraffe et al. (2021), because the
background stellar models are significantly different. (For instance the
boosted simulations have convection zone with much smaller super-adiabaticity
$\sim 10^{-8}$, deeper convective envelope $r_{\textrm{Schw}}\sim
0.68R_{\ast}$, smaller stellar radius
$R_{\ast}\sim$0.8\text{\,}\mathrm{R}_{\sun}$$ and larger stellar luminosity).
However, as an indication we note that the difference in the overshooting
depths between the $r_{\max}=0.9R_{\ast}$ and $r_{\max}=0.99R_{\ast}$
simulations are comparable to the effect of boosting the luminosity by a
factor of 100 in the boosted simulations. Note that Baraffe et al. (2021)
caution against using significantly larger boost factors in any case, as they
may lead to evolution of the background model and spurious changes of e. g.
the IGW spectrum. Similarly, extrapolating the presented results to values of
$r_{\max}$ significantly above the $0.99R_{\ast}$ level is challenging,
because this would include surface layers where the opacity quickly decreases
and additional physics would need to be considered (e. g. the increasing
impact of radiative cooling requires a detailed treatment of radiative
transport).

###  6.2 Outlook

Because of the steep stratification it is computationally challenging to
perform reasonably well-resolved global stellar hydrodynamic simulations with
$r_{\max}=1R_{\ast}$ with realistic radiation physics and self-consistent
surface boundary layers. Hence quantitative predictions of convective
velocities, IGW spectra and overshooting depths should always be considered
with great care. This is a big part of the challenge of solving the so-called
solar convective conundrum (Hanasoge et al., 2012; Schumacher & Sreenivasan,
2020; Hotta & Kusano, 2021; Vasil et al., 2021) and is the reason to refrain
from a more quantitative analysis of the data. As $r_{\max}$ is increased
towards unity the time and length scales of the boundary layers decrease while
the cooling efficiency increase. Eventually, an asymptotic limit must be
reached for the amplitude of the radial velocity and temperature fluctuations
in the convection zone (and consequently for the overshooting depth and IGW
power spectrum). However, it would be computationally unfeasible to seek such
a limit with the uniform radial grids used in this study. Encouragingly, Hotta
et al. (2019) find similar profiles of the root-mean-square velocities in
local convection zone simulations truncated at $r_{\max}=0.992R_{\ast}$ and
$r_{\max}=1R_{\ast}$, when the former is augmented with an artificial surface
cooling layer. They attribute the similarity of the results to efficient
mixing in the near-surface region, which reduces the non-locality of the
convection. This indicates that it may be possible to model the influence of
the truncated, and difficult to resolve, surface boundary layers. Further
studies are needed to examine the influence that such modelling may have on
the overshooting depth and IGWs. It also remains to be established how the
presented non-local dependencies translate to 3D convection and how they
interact with other important processes like rotation and magnetism, for which
a follow-up study is currently under way. However, the qualitative trends this
study has highlighted (the significant dependence on $r_{\max}$ and the weaker
dependence on $r_{\min}$) should be robust results that can be linked to
physical mechanisms which operate in the stellar hydrodynamics context.

## Acknowledgements

This work is supported by the ERC grant No. 787361-COBOM and the consolidated
STFC grant ST/R000395/1. The authors would like to acknowledge the use of the
University of Exeter High-Performance 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 e-Infrastructure.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Alvan et al. (2014) Alvan L., Brun A. S., Mathis S., 2014, Astronomy & Astrophysics, 565, A42
  * Anders et al. (2019) Anders E. H., Lecoanet D., Brown B. P., 2019, arXiv:1906.02342 [astro-ph, physics:physics]
  * Augustson et al. (2015) Augustson K., Brun A. S., Miesch M., Toomre J., 2015, The Astrophysical Journal, 809, 149
  * Baraffe & El Eid (1991) Baraffe I., El Eid M. F., 1991, Astronomy & Astrophysics, 245, 548
  * Baraffe et al. (1998) Baraffe I., Chabrier G., Allard F., Hauschildt P., 1998, arXiv:astro-ph/9805009
  * Baraffe et al. (2017) Baraffe I., Pratt J., Goffrey T., Constantino T., Folini D., Popov M. V., Walder R., Viallet M., 2017, The Astrophysical Journal, 845, L6
  * Baraffe et al. (2021) Baraffe I., Pratt J., Vlaykov D. G., Guillet T., Goffrey T., Saux A. L., Constantino T., 2021, arXiv:2109.00943 [astro-ph]
  * Bedding et al. (2010) Bedding T. R., et al., 2010, The Astrophysical Journal, 713, L176
  * Birch et al. (2018) Birch A., Duvall T., Gizon L., Hanasoge S., Hindman B. W., Nagashima K., Sreenivasan K. R., 2018, in Catalyzing Solar Connections. Ghent, Belgium, p. 42
  * Bowman et al. (2019) Bowman D. M., et al., 2019, Nature Astronomy, 3, 760
  * Brandenburg (2016) Brandenburg A., 2016, The Astrophysical Journal, 832, 6
  * Brummell et al. (2002) Brummell N. H., Clune T. L., Toomre J., 2002, The Astrophysical Journal, 570, 825
  * Brun et al. (2011) Brun A. S., Miesch M. S., Toomre J., 2011, The Astrophysical Journal, 742, 79
  * Brun et al. (2017) Brun A. S., et al., 2017, The Astrophysical Journal, 836, 192
  * Canuto & Dubovikov (1997) Canuto V. M., Dubovikov M., 1997, The Astrophysical Journal, 484, L161
  * Constantino et al. (2021) Constantino T., Baraffe I., Goffrey T., Pratt J., Guillet T., Vlaykov D., Amard L., 2021, arXiv:2108.10361 [astro-ph]
  * Cossette & Rast (2016) Cossette J.-F., Rast M. P., 2016, The Astrophysical Journal, 829, L17
  * Dintrans et al. (2005) Dintrans B., Brandenburg A., Nordlund A., Stein R. F., 2005, Astronomy & Astrophysics, 438, 365
  * Featherstone & Hindman (2016) Featherstone N. A., Hindman B. W., 2016, The Astrophysical Journal, 818, 32
  * Goffrey et al. (2017) Goffrey T., et al., 2017, Astronomy & Astrophysics, 600, A7
  * Goldreich & Kumar (1990) Goldreich P., Kumar P., 1990, The Astrophysical Journal, 363, 694
  * Gough (1977) Gough D. O., 1977, in Spiegel E. A., Zahn J.-P., eds, , Vol. 71, Problems of Stellar Convection. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 349–363, doi:10.1007/3-540-08532-7_56
  * Guerrero et al. (2013) Guerrero G., Smolarkiewicz P. K., Kosovichev A. G., Mansour N. N., 2013, The Astrophysical Journal, 779, 176
  * Hanasoge et al. (2012) Hanasoge S. M., Duvall T. L., Sreenivasan K. R., 2012, Proceedings of the National Academy of Sciences, 109, 11928
  * Higl et al. (2021) Higl J., Müller E., Weiss A., 2021, Astronomy & Astrophysics, 646, A133
  * Horst et al. (2020) Horst L., Edelmann P. V. F., Andrássy R., Röpke F. K., Bowman D. M., Aerts C., Ratnasingam R. P., 2020, Astronomy & Astrophysics, 641, A18
  * Hotta & Kusano (2021) Hotta H., Kusano K., 2021, Nature Astronomy
  * Hotta et al. (2019) Hotta H., Iijima H., Kusano K., 2019, Science Advances, 5, eaau2307
  * Hurlburt et al. (1984) Hurlburt N. E., Toomre J., Massaguer J. M., 1984, The Astrophysical Journal, 282, 557
  * Hurlburt et al. (1994) Hurlburt N. E., Toomre J., Massaguer J. M., Zahn J.-P., 1994, The Astrophysical Journal, 421, 245
  * Iglesias & Rogers (1996) Iglesias C. A., Rogers F. J., 1996, The Astrophysical Journal, 464, 943
  * Käpylä et al. (2017) Käpylä P. J., Rheinhardt M., Brandenburg A., Arlt R., Käpylä M. J., Lagg A., Olspert N., Warnecke J., 2017, The Astrophysical Journal, 845, L23
  * Käpylä et al. (2019) Käpylä P. J., Viviani M., Käpylä M. J., Brandenburg A., Spada F., 2019, Geophysical & Astrophysical Fluid Dynamics, 0, 1
  * Korre et al. (2019) Korre L., Garaud P., Brummell N. H., 2019, Monthly Notices of the Royal Astronomical Society, 484, 1220
  * Le Saux et al. (2021) Le Saux A., et al., 2021, Astronomy & Astrophysics, submitted
  * Lecoanet & Quataert (2013) Lecoanet D., Quataert E., 2013, Monthly Notices of the Royal Astronomical Society, 430, 2363
  * Lecoanet et al. (2019) Lecoanet D., et al., 2019, The Astrophysical Journal, 886, L15
  * Meakin & Arnett (2007) Meakin C. A., Arnett D., 2007, The Astrophysical Journal, 667, 448
  * Muthsam et al. (1995) Muthsam H. J., Goeb W., Kupka F., Liebich W., Zoechling J., 1995, Astronomy & Astrophysics, 293, 127
  * Pinçon et al. (2016) Pinçon C., Belkacem K., Goupil M. J., 2016, Astronomy & Astrophysics, 588, A122
  * Pratt et al. (2016) Pratt J., et al., 2016, Astronomy & Astrophysics, 593, A121
  * Pratt et al. (2017) Pratt J., Baraffe I., Goffrey T., Constantino T., Viallet M., Popov M. V., Walder R., Folini D., 2017, Astronomy & Astrophysics, 604, A125
  * Press (1981) Press W. H., 1981, The Astrophysical Journal, 245, 286
  * Rieutord & Zahn (1995) Rieutord M., Zahn J.-P., 1995, Astronomy & Astrophysics, 296, 127
  * Rogers & Glatzmaier (2006) Rogers T. M., Glatzmaier G. A., 2006, The Astrophysical Journal, 653, 765
  * Rogers & McElwaine (2017) Rogers T. M., McElwaine J. N., 2017, The Astrophysical Journal, 848, L1
  * Rogers & Nayfonov (2002) Rogers F. J., Nayfonov A., 2002, The Astrophysical Journal, 576, 1064
  * Rogers et al. (2013) Rogers T. M., Lin D. N. C., McElwaine J. N., Lau H. H. B., 2013, The Astrophysical Journal, 772, 21
  * Schumacher & Sreenivasan (2020) Schumacher J., Sreenivasan K. R., 2020, Reviews of Modern Physics, 92, 041001
  * Spruit (1997) Spruit H. C., 1997, Memorie della Societa Astronomica Italiana, 68, 397
  * Townsend & Teitler (2013) Townsend R. H. D., Teitler S. A., 2013, Monthly Notices of the Royal Astronomical Society, 435, 3406
  * Townsend et al. (2018) Townsend R. H. D., Goldstein J., Zweibel E. G., 2018, Monthly Notices of the Royal Astronomical Society, 475, 879
  * Vasil et al. (2021) Vasil G. M., Julien K., Featherstone N. A., 2021, Proceedings of the National Academy of Sciences, 118, e2022518118
  * Viallet et al. (2011) Viallet M., Baraffe I., Walder R., 2011, Astronomy & Astrophysics, 531, A86
  * Viallet et al. (2015) Viallet M., Meakin C., Prat V., Arnett D., 2015, Astronomy & Astrophysics, 580, A61
  * Viallet et al. (2016) Viallet M., Goffrey T., Baraffe I., Folini D., Geroux C., Popov M. V., Pratt J., Walder R., 2016, Astronomy & Astrophysics, 586, A153
  * Wagner & Shishkina (2013) Wagner S., Shishkina O., 2013, Physics of Fluids, 25, 085110
  * Zahn (1991) Zahn J.-P., 1991, Astronomy & Astrophysics, pp 179–188

##  Appendix A Convective heat flux derivation

In this work we use Eq. 8 as a definition of the convective heat flux. This is
derived from the angular average of the total energy density budget

$\displaystyle\partial_{t}\left<\textrm{E}_{\mathrm{t}}\right>_{\theta}$
$\displaystyle=-\partial_{r}\left<\frac{1}{2}\rho
v_{r}\boldsymbol{v}^{2}-\chi\nabla T+\rho v_{r}h+\rho
v_{r}\Phi\right>_{\theta}$ (14)
$\displaystyle=-\partial_{r}\left<\boldsymbol{F}^{\mathrm{ke}}_{r}+\boldsymbol{F}^{\chi}_{r}+\boldsymbol{F}^{\mathrm{conv}}_{r}+\boldsymbol{F}^{\textrm{bkgr}}_{r}\right>_{\theta},$
(15)

where we introduce $\textrm{E}_{\mathrm{t}}$ as the sum of the kinetic,
gravitational and internal energy densities,
$\boldsymbol{F}^{\chi}_{r}=-\chi\nabla T$ as the radiative flux and
$\boldsymbol{F}^{\textrm{bkgr}}_{r}=\left<\rho
v_{r}\right>_{\theta}\left<\Phi+h\right>_{\theta}$ as the flux due to the mean
stratification. We recall that
$\boldsymbol{F}^{\mathrm{conv}}_{r}=(\rho\boldsymbol{v})^{\prime}h^{\prime}$
and ′ designates deviations away from the instantaneous angular average. In
principle, given that in the convection zone $\left<h\right>\gg h^{\prime}$,
the contribution of the mean specific enthalpy to the energy transport may be
non-negligible in the not-perfectly-relaxed case even with a small mean radial
momentum. However, this is largely compensated by the hydrostatic
stratification through the mean potential energy $\left<\Phi\right>$. Using
the _fundamental thermodynamic relation_
$\mathrm{d}h=T\mathrm{d}s+\rho^{-1}\mathrm{d}p$ (with $s$ the specific
entropy) we can express $h^{\prime\prime}=\left<h+\Phi\right>_{\theta}$, the
mean enthalpy excess away from hydrostatic equilibrium (HSE), as

$\displaystyle\partial_{r}\left<h^{\prime\prime}\right>_{\theta}$
$\displaystyle=\left<T\nabla s+\rho^{-1}\nabla p-\rho^{-1}\nabla
p_{HSE}\right>_{\theta}$ (16) $\displaystyle=\left<T\nabla s+\rho^{-1}\nabla
p^{\prime\prime}\right>_{\theta},$ (17)

where we define $\nabla p_{HSE}=-\rho\nabla\Phi$ as the pressure gradient in
HSE and $p^{\prime\prime}=p-p_{HSE}$ as the pressure perturbations away from
it. This illustrates that mean enthalpy flux contains pressure perturbations
with respect to HSE and entropy perturbations. Given a dynamical time scale of
$10^{3}\,$\mathrm{s}$-10^{4}\,$\mathrm{s}$$ in the convection zone (depending
on radius), any deviations of the mean stratification from HSE are quickly
corrected, so the $p^{\prime\prime}$ term can be neglected over the time
scales considered in this study. The mean entropy perturbations are at the
order of $10^{-5}$ and are therefore similarly negligible. This justifies
neglecting the background flux $\boldsymbol{F}^{bkgr}_{r}$ and considering
$\boldsymbol{F}^{\mathrm{conv}}_{r}$ as the convective heat flux.