# Stellar convective penetration: parameterized theory and dynamical
simulations

Evan H. Anders CIERA, Northwestern University, Evanston IL 60201, USA  Adam S.
Jermyn Center for Computational Astrophysics, Flatiron Institute, New York, NY
10010, USA  Daniel Lecoanet CIERA, Northwestern University, Evanston IL 60201,
USA  Department of Engineering Sciences and Applied Mathematics, Northwestern
University, Evanston IL 60208, USA  Benjamin P. Brown Department Astrophysical
and Planetary Sciences & LASP, University of Colorado, Boulder, CO 80309, USA
Evan H. Anders<EMAIL_ADDRESS>
(Received July 28, 2021; Revised October 19, 2021)

###### Abstract

Most stars host convection zones in which heat is transported directly by
fluid motion. Parameterizations like mixing length theory adequately describe
convective flows in the bulk of these regions, but the behavior of convective
boundaries is not well understood. Here we present 3D numerical simulations
which exhibit penetration zones: regions where the entire luminosity _could_
be carried by radiation, but where the temperature gradient is approximately
adiabatic and convection is present. To parameterize this effect, we define
the “penetration parameter” $\mathcal{P}$ which compares how far the radiative
gradient deviates from the adiabatic gradient on either side of the
Schwarzschild convective boundary. Following Roxburgh (1989) and Zahn (1991),
we construct an energy-based theoretical model in which the extent of
penetration is controlled by $\mathcal{P}$. We test this theory using 3D
numerical simulations which employ a simplified Boussinesq model of stellar
convection. We find significant convective penetration in all simulations. Our
simple theory describes the simulations well. Penetration zones can take
thousands of overturn times to develop, so long simulations or accelerated
evolutionary techniques are required. In stellar contexts, we expect
$\mathcal{P}\approx 1$ and in this regime our results suggest that convection
zones may extend beyond the Schwarzschild boundary by up to $\sim$20-30% of a
mixing length. We present a MESA stellar model of the Sun which employs our
parameterization of convective penetration as a proof of concept. We discuss
prospects for extending these results to more realistic stellar contexts.

UAT keywords

††journal: ApJ

##  1 Introduction

###  1.1 Context

Convection is a crucial mechanism for transporting heat in stars (Woosley et
al., 2002; Hansen et al., 2004; Christensen-Dalsgaard, 2021), and convective
dynamics influence many poorly-understood stellar phenomena. For example,
convection drives the magnetic dynamo of the Sun, leading to a whole host of
emergent phenomena collectively known as solar activity (Brun & Browning,
2017). Convection also mixes chemical elements in stars, which can modify
observed surface abundances or inject additional fuel into their cores,
thereby extending stellar lifetimes (Salaris & Cassisi, 2017). Furthermore,
convective motions excite waves, which can be observed and used to constrain
the thermodynamic structure of stars (Aerts et al., 2010; Basu, 2016). A
complete and nuanced understanding of convection is therefore crucial for
understanding stellar structure and evolution, and for connecting this
understand to observations.

Despite decades of study, robust parameterizations for the mechanisms broadly
referred to as “convective overshoot” remain elusive, and improved
parameterizations could resolve many discrepancies between observations and
structure models. In the stellar structure literature, “convective overshoot”
refers to any convectively-driven mixing which occurs beyond the boundaries of
the Ledoux-unstable zone. This mixing can influence, for example, observed
surface lithium abundances in the Sun and solar-type stars, which align poorly
with theoretical predictions (Pinsonneault, 1997; Carlos et al., 2019; Dumont
et al., 2021). Furthermore, modern spectroscopic observations suggest a lower
solar metallicity than previously thought, and models computed with modern
metallicity estimates and opacity tables have shallower convection zones than
helioseismic observations suggest (Basu & Antia, 2004; Bahcall et al., 2005;
Bergemann & Serenelli, 2014; Vinyoles et al., 2017; Asplund et al., 2021);
modeling and observational discrepancies can be reduced with additional mixing
below the convective boundary (Christensen-Dalsgaard et al., 2011).

Beyond the Sun, overshooting in massive stars with convective cores must be
finely tuned as a function of stellar mass, again pointing to missing physics
in our current parameterizations (Claret & Torres, 2018; Jermyn et al., 2018;
Viani & Basu, 2020; Martinet et al., 2021; Pedersen et al., 2021). Since core
convective overshoot increases the reservoir of fuel available for nuclear
fusion at each stage in stellar evolution, improved models of core convective
boundary mixing could have profound impacts on the post-main sequence
evolution and remnant formation of massive stars (Farmer et al., 2019; Higgins
& Vink, 2020).

In order to ensure that models can be evolved on fast (human) timescales, 1D
stellar evolution codes rely on simple parameterizations of convection (e.g.,
mixing length theory, Böhm-Vitense, 1958) and convective overshoot (Shaviv &
Salpeter, 1973; Maeder, 1975; Herwig, 2000; Paxton et al., 2011, 2013, 2018,
2019). While some preliminary work has been done to couple 3D dynamical
convective simulations with 1D stellar evolution codes (Jørgensen & Weiss,
2019), these calculations are prohibitively expensive to perform at every
timestep in a stellar evolution simulation. To resolve discrepancies between
stellar evolution models and observations, a more complete and
_parameterizeable_ understanding of convective overshoot is required.

The broad category of “convective overshoot” in the stellar literature is an
umbrella term for a few hydrodynamical processes (Zahn, 1991; Brummell et al.,
2002; Korre et al., 2019). Motions which extend beyond the convective boundary
but do not adjust the thermodynamic profiles belong to a process called “
_convective overshoot_ ” in the fluid dynamics literature. Convection zones
can expand through a second process called “ _entrainment_ ,” through which
motions erode composition gradients or modify the radiative gradient (Meakin &
Arnett, 2007; Viallet et al., 2013; Cristini et al., 2017; Fuentes & Cumming,
2020; Horst et al., 2021). The primary focus of this work is a third process
called “ _convective penetration_ ”. Convective penetration occurs when
motions mix the entropy gradient towards the adiabatic in a region that is
stable by the Schwarzschild criterion.

Figure 1:  Vertical slice through a simulation with $\mathcal{R}=6.4\times
10^{3}$, $\mathcal{P}_{D}=4$ and $\mathcal{S}=10^{3}$ (see Sec. 21). The
dashed horizontal line denotes the Schwarzschild convective boundary where
$\nabla_{\rm{ad}}=\nabla_{\rm{rad}}$. The top of the penetrative zone
($\delta_{0.1}$, see Sec. 21) is shown by a solid horizontal line. (Left)
Vertical velocity is shown; orange convective upflows extend far past the
Schwarzschild boundary of the convection zone but stop abruptly at the top of
the penetration zone where $\nabla$ departs from $\nabla_{\rm{ad}}$. (Right)
Temperature fluctuations, normalized by their average magnitude at each height
to clearly display all dynamical features.

Convective overshoot, entrainment, and penetration have been studied in the
laboratory and through numerical simulations for decades, and the state of the
field has been regularly reviewed (e.g., Marcus et al., 1983; Zahn, 1991;
Browning et al., 2004; Rogers et al., 2006; Viallet et al., 2015; Korre et
al., 2019). Experiments exhibiting extensive expansion of convection zones via
entrainment have a long history (e.g., Musman, 1968; Deardorff et al., 1969;
Moore & Weiss, 1973, and this process is often confusingly called
“penetration”). Modern numerical experiments often examine the importance of
the “stiffness” $\mathcal{S}$ of a radiative-convective interface.
$\mathcal{S}$ compares the relative stability of a radiative zone and an
adjacent convection zone according to some measure like a dynamical frequency
or characteristic entropy gradient. Some recent studies in simplified
Boussinesq setups exhibit stiffness-dependent convection zone expansion via
entrainment (Couston et al., 2017; Toppaladoddi & Wettlaufer, 2018); others
find stiffness-dependent pure overshoot (Korre et al., 2019). A link between
$\mathcal{S}$ and the processes of entrainment and overshoot has seemingly
emerged, but a mechanism for penetration remains elusive.

Many studies in both Cartesian and spherical geometries have exhibited hints
of penetrative convection. Some authors report clear mixing of the entropy
gradient beyond the nominal convecting region (Hurlburt et al., 1994; Saikia
et al., 2000; Brummell et al., 2002; Rogers & Glatzmaier, 2005; Rogers et al.,
2006; Kitiashvili et al., 2016), but it is often unclear how much mixing is
due to changes in the location of the Schwarzschild boundary (entrainment) and
how much is pure penetration. Other authors present simulations with dynamical
or flux-based hints of penetration such as a negative convective flux or a
radiative flux which exceeds the total system flux, but do not clearly report
the value of the entropy gradient (Hurlburt et al., 1986; Singh et al., 1995;
Browning et al., 2004; Brun et al., 2017; Pratt et al., 2017). Still other
simulations show negligible penetration (e.g., Higl et al., 2021). Even
detailed studies which sought a relationship between penetration depth and
stiffness $\mathcal{S}$ have presented contradictory results. Early work by
e.g., Hurlburt et al. (1994) and Singh et al. (1995) hinted at a link between
$\mathcal{S}$ and penetration length, at least for low values of
$\mathcal{S}$. Subsequent simulations by Brummell et al. (2002) exhibit a weak
scaling of penetration depth with $\mathcal{S}$; the authors interpret this
scaling as a sign of pure overshoot and claim their simulations do not achieve
adiabatic convective penetration. Still later simulations by Rogers &
Glatzmaier (2005) demonstrate a negligible scaling of the penetration depth
against $\mathcal{S}$ at moderate values of $\mathcal{S}$. Prior simulations
thus consistently show hints of penetration at low $\mathcal{S}$ (where
results may not be relevant for stars, Couston et al., 2017), but present
confusing and contradictory results at moderate-to-high $\mathcal{S}$.

There are hints in the literature that convective penetration may depend on
energy fluxes. Roxburgh (1978, 1989, 1992, 1998) derived an “integral
constraint” from the energy equation and found that a spatial integral of the
flux puts an upper limit on the size of a theoretical penetrative region. Zahn
(1991) theorized that convective penetration should depend only on how steeply
the radiative temperature gradient varies at the convective boundary.
Following Zahn (1991)’s work, Rempel (2004) derived a semianalytic model and
suggested that inconsistencies seen in simulations of penetrative dynamics can
be explained by the magnitude of the fluxes or luminosities driving the
simulations. Indeed, some simulations have tested this idea, and found that
penetration lengths depend strongly on the input flux (Singh et al., 1998;
Käpylä et al., 2007; Tian et al., 2009; Hotta, 2017; Käpylä, 2019).
Furthermore, in the limit of low stiffness, the simulations of Hurlburt et al.
(1994) and Rogers & Glatzmaier (2005) may agree with Zahn’s theory (although
at high stiffness they disagree). In light of these results, and the possible
importance of energy fluxes, Roxburgh’s integral constraint and Zahn’s theory
deserve to be revisited.

###  1.2 Convective penetration & this study’s findings

Convective penetration is the process by which convective motions extend
beyond the Schwarzschild-stable boundary and mix the entropy gradient to be
nearly adiabatic.

> _In this paper, we present simulations which exhibit convective
> penetration._

This process is phenomenologically described in Sec. 2. In this work, the
convection zone lies beneath an adjacent stable layer and convection
penetrates upwards; our results equally apply to the reversed problem.

In order to understand this phenomenon, we derive theoretical predictions for
the size of the penetrative zone based on the ideas of Roxburgh (1989) and
Zahn (1991).

> _We find that the extent of convective penetration depends strongly on the
> shape and magnitude of the radiative gradient near the convective boundary._

Thus, the penetration length can be calculated using the radiative
conductivity (or opacity) _profile_ near the convective boundary.

We present these findings as follows. In Sec. 2, we present the central
finding of this work: penetration zones in nonlinear convective simulations.
In Sec. 3, we describe the equations used and derive a parameterized theory of
convective penetration. In Sec. 21, we describe our simulation setup and
parameters. In Sec. 5, we present the results of these simulations, with a
particular focus on the height of the penetrative regions. In Sec. 6, we
create and discuss a stellar model in MESA which has convective penetration.
Finally, we discuss pathways for future work in Sec. 7.

##  2 Central Result: Convective Penetration

In Fig. 1, we display a snapshot of dynamics in an evolved simulation which
exhibits convective penetration. The simulation domain is a 3D Cartesian box,
and this figure shows a vertical slice through the center of the domain. In
the left panel, we display the vertical velocity. We see that convective
motions extend beyond the Schwarzschild boundary of the convection zone, which
is denoted by a horizontal dashed grey line. These motions stop at the top of
a penetration zone, denoted by a solid horizontal line, where the temperature
gradient departs from adiabatic towards the radiative gradient. In the right
panel, we display temperature perturbations away from the time-evolving mean
temperature profile. We see that warm upwellings in the Schwarzschild-unstable
convection zone (below the dashed line) become cold upwellings in the
penetration zone (above the dashed line), and these motions excite gravity
waves in the stable radiative zone (above the solid line).

We further explore the simulation from Fig. 1 in Fig. 2 by displaying time-
and horizontally-averaged 1D profiles of the temperature gradient $\nabla$
(defined in Sec. 3). The adiabatic gradient $\nabla_{\rm{ad}}$ (purple) has a
constant value in the simulation. Also shown is the radiative gradient
$\nabla_{\rm{rad}}$ (orange). The domain exhibits a classical Schwarzshild-
unstable convection zone (CZ) for $z\lesssim 1.04$ where
$\nabla_{\rm{rad}}>\nabla_{\rm{ad}}$; the upper boundary of this region is
denoted by a dashed vertical line. Above this point,
$\nabla_{\rm{rad}}<\nabla_{\rm{ad}}$ and the domain would be considered stable
by the Schwarzschild criterion. However, the evolved convective dynamics in
Fig. 1 have raised $\nabla\rightarrow\nabla_{\rm{ad}}$ in an extended
penetration zone (PZ) which extends from $1.04\lesssim z\lesssim 1.3$. Above
$z\gtrsim 1.4$, $\nabla\approx\nabla_{\rm{rad}}$ in a classical stable
radiative zone (RZ). Between $1.3\lesssim z\lesssim 1.4$, there is a PZ-RZ
boundary layer (referred to as the “thermal adjustment layer” in some prior
studies) where convective motions give way to conductive transport and
$\nabla$ adjusts from $\nabla_{\rm{ad}}$ to $\nabla_{\rm{rad}}$.

Our goals in this paper are to understand how these PZs form and to
parameterize this effect so that it can be included in 1D stellar evolution
calculations.

Figure 2:  Horizontally- and temporally-averaged profiles of the thermodynamic
gradients from the simulation in Fig. 1. We plot $\nabla$ (green) compared to
$\nabla_{\rm{ad}}$ (purple, a constant) and $\nabla_{\rm{rad}}$ (orange); note
the extended penetration zone (PZ) where
$\nabla\approx\nabla_{\rm{ad}}>\nabla_{\rm{rad}}$. The dashed vertical line
denotes the Schwarzschild boundary of the convection zone (CZ), the solid
vertical line denotes the bottom of the radiative zone (RZ), and the greyed
region denotes the PZ-RZ boundary layer.

##  3 Theory

In this section we derive a theoretical model of convective penetration by
examining the energetics and energy fluxes in the Schwarzschild-unstable
convection zone (CZ) and penetration zone (PZ). In Sec. 3.1, we describe our
equations and problem setup and define the heat fluxes. In Sec. 3.2, we build
a parameterized theory based on the kinetic energy (KE) equation. We find that
imbalances in KE source terms within the CZ determine the extent of the PZ. By
balancing the excess KE generation in the CZ with buoyant deceleration and
dissipation work terms in the PZ, we are able to derive the size of the PZ. We
find that a description of the size of a theoretical PZ does not depend on the
often-considered stiffness, which measures the relative stability between the
convection zone and an adjacent radiative zone.

###  3.1 Equations & flux definitions

Throughout this work, we will utilize a modified version of the incompressible
Boussinesq equations,

$\displaystyle\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0$ (1)
$\displaystyle\partial_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}=-\frac{1}{\rho_{0}}\boldsymbol{\nabla}p+\frac{\rho_{1}}{\rho_{0}}\boldsymbol{g}+\nu\boldsymbol{\nabla}^{2}\boldsymbol{u}$
(2)
$\displaystyle\partial_{t}T+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}T+w\nabla_{\rm{ad}}+\boldsymbol{\nabla}\boldsymbol{\cdot}[-k\boldsymbol{\nabla}\overline{T}]=\chi\boldsymbol{\nabla}^{2}T^{\prime}+Q$
(3) $\displaystyle\frac{\rho_{1}}{\rho_{0}}=-|\alpha|T.$ (4)

Here, the density is decomposed into a uniform, constant background $\rho_{0}$
with fluctuations $\rho_{1}$ which appear only in the buoyancy force and
depend on the temperature $T$ and the coefficient of thermal expansion
$\alpha=\partial\ln\rho/\partial T$. We define the velocity vector
$\boldsymbol{u}$, the pressure $p$, the viscous diffusivity $\nu$, the thermal
diffusivity $\chi$, the bulk internal heating $Q$, the adiabatic gradient
$\nabla_{\rm{ad}}$, and a height-dependent thermal conductivity $k$. We will
consider Cartesian coordinates $(x,y,z)$ with a constant vertical gravity
$\boldsymbol{g}=-g\hat{z}$. Throughout this work, we will represent horizontal
averages with bars ($\overline{\,\cdot\,}$) and fluctuations away from those
averages with primes (′). Thus, in Eqn. 3, $\overline{T}$ is the horizontally
averaged temperature and $T^{\prime}$ are fluctuations away from that; both of
these fields evolve in time according to Eqn. 3.

Assuming convection reaches a time-stationary state, the heat fluxes are found
by horizontally-averaging then vertically integrating Eqn. 3 to find

$\overline{F_{\rm{tot}}}=\overline{F_{\rm{rad}}}+\overline{F_{\rm{conv}}}=\int
Qdz+F_{\rm{bot}},$ (5)

where $F_{\rm{bot}}$ is the flux carried at the bottom of the domain, and
$\overline{F_{\rm{tot}}}$ is the total flux, which can vary in height due to
the heating $Q$. The mean temperature profile $\overline{T}$ carries the
radiative flux $\overline{F_{\rm{rad}}}=-k\boldsymbol{\nabla}\overline{T}$. We
note that $k$ and $-\partial_{z}\overline{T}$ fully specify
$\overline{F_{\rm{rad}}}$ and in turn the convective flux,
$\overline{F_{\rm{conv}}}=\overline{F_{\rm{tot}}}-\overline{F_{\rm{rad}}}$. We
define the temperature gradient and radiative temperature gradient

$\nabla\equiv-\partial_{z}\overline{T}\qquad\nabla_{\rm{rad}}\equiv\frac{\overline{F_{\rm{tot}}}}{k}.$
(6)

We have defined the $\nabla$’s as positive quantities to align with stellar
structure conventions and intuition. Marginal stability is achieved when
$\nabla=\nabla_{\rm{ad}}$, which we take to be a constant. We note that the
classical Schwarzschild boundary of the convection zone is the height
$z=L_{s}$ at which $\nabla_{\rm{rad}}=\nabla_{\rm{ad}}$ and
$\overline{F_{\rm{conv}}}=0$.

The addition of a nonzero $\nabla_{\rm{ad}}$ to Eqn. 3 was derived by Spiegel
& Veronis (1960) and utilized by e.g., Korre et al. (2019). In this work, we
have decomposed the radiative diffusivity into a background portion
($\boldsymbol{\nabla}\boldsymbol{\cdot}\overline{F_{\rm{rad}}}$) and a
fluctuating portion ($\chi\boldsymbol{\nabla}^{2}T^{\prime}$); by doing so, we
have introduced a height-dependent $\nabla_{\rm{rad}}$ to the equation set
while preserving the diffusive behavior on fluctuations felt by classical
Rayleigh-Bénard convection. Here, we will assume a model in which an unstable
convection zone ($\nabla_{\rm{rad}}>\nabla_{\rm{ad}}$) sits below a stable
radiative zone ($\nabla_{\rm{rad}}<\nabla_{\rm{ad}}$), but in this
incompressible model where there is no density stratification to break the
symmetry of upflows and downflows, precisely the same arguments can be applied
to the inverted problem.

###  3.2 Kinetic energy & the dissipation-flux link

Figure 3:  Temporally- and horizontally-averaged profiles from Eqn. 8 in the
simulation in Fig. 1. The vertical dashed line denotes the Schwarzschild CZ
boundary, and the vertical solid line corresponds to the top of the PZ.
(upper) Kinetic energy fluxes $\overline{\mathcal{F}_{z}}$, which go to zero
at the bottom boundary and the top of the PZ. (bottom) Source terms from Eqn.
8 normalized by the maximum of $\overline{\mathcal{B}}$
($\overline{\mathcal{F}_{z}}$ in the upper panel is similarly normalized). The
buoyancy source $\overline{\mathcal{B}}$ changes sign at the Schwarzschild
boundary, and $\overline{\Phi}$ is positive-definite.

Taking a dot product of the velocity and Eqn. 2 reveals the kinetic energy
equation,

$\frac{\partial\mathcal{K}}{\partial
t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathcal{F}}=\mathcal{B}-\Phi,$
(7)

where we define the kinetic energy $\mathcal{K}\equiv|\boldsymbol{u}|^{2}/2$,
the fluxes of kinetic energy
$\boldsymbol{\mathcal{F}}\equiv\left[\boldsymbol{u}(\mathcal{K}+p/\rho_{0})-\nu\boldsymbol{u}\boldsymbol{\times}\boldsymbol{\omega}\right]$,
the buoyant energy generation rate $\mathcal{B}\equiv|\alpha|gwT^{\prime}$,
and the viscous dissipation rate $\Phi\equiv\nu|\boldsymbol{\omega}|^{2}$
where
$\boldsymbol{\omega}=\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{u}$ is
the vorticity and
$|\boldsymbol{u}|^{2}=\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}$ &
$|\boldsymbol{\omega}|^{2}=\boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\omega}$.
We next take a horizontal- and time-average of Eqn. 7 (we absorb the time-
average into the horizontal-average $\overline{\,\cdot\,}$ notation for
simplicity). Assuming that $\overline{\mathcal{K}}$ reaches a statistically
stationary state, convective motions satisfy

$\frac{d\overline{\mathcal{F}_{z}}}{dz}=\overline{\mathcal{B}}-\overline{\Phi},$
(8)

where $\mathcal{F}_{z}$ is the z-component of $\boldsymbol{\mathcal{F}}$. Each
profile in Eqn. 8 is shown in Fig. 3 for the simulation whose dynamics are
displayed in Fig. 1. As in Fig. 2, the Schwarzschild CZ boundary is plotted as
a dashed line, and the top of the PZ is plotted as a solid vertical line. In
the top panel, we display $\overline{\mathcal{F}_{z}}$, neglecting the viscous
flux term which is only nonzero in a small region above the bottom boundary.
We see that $\overline{\mathcal{F}_{z}}$ is zero at the bottom boundary (left
edge of plot) and at the top of the PZ. In the bottom panel, we plot
$\overline{\mathcal{B}}$ and $\overline{\Phi}$; we see that
$\overline{\mathcal{B}}$ changes sign at the Schwarzschild CZ boundary, and
that $\overline{\Phi}$ is positive-definite.

At the boundaries of the convecting region, $\overline{\mathcal{F}_{z}}$ is
zero (Fig. 3, upper panel). We integrate Eqn. 8 vertically between these zeros
to find

$\int\overline{\mathcal{B}}\,dz=\int\overline{\Phi}\,dz.$ (9)

Integral constraints of this form are the basis for a broad range of analyses
in Boussinesq convection (see e.g., Ahlers et al., 2009; Goluskin, 2016) and
were considered in the context of penetrative stellar convection by Roxburgh
(1989). Eqn. 9 is the straightforward statement that work by buoyancy on large
scales must be balanced by viscous dissipation on small scales.

We break up the convecting region into a Schwarzschild-unstable “convection
zone” (CZ) and an extended “penetration zone” (PZ); we assume that convective
motions efficiently mix $\nabla\rightarrow\nabla_{\rm{ad}}$ in both the CZ and
PZ. The buoyant energy generation is proportional to the convective flux,
$\overline{\mathcal{B}}=|\alpha|g\overline{wT^{\prime}}=|\alpha|g\overline{F_{\rm{conv}}}$,
and is positive in the CZ and negative in the PZ (see Fig. 3, bottom panel).
Breaking up Eqn. 9, we see that

$\int_{\rm{CZ}}\overline{\mathcal{B}}\,dz\,\,=\,\,\int_{\rm{CZ}}\overline{\Phi}\,dz+\int_{\rm{PZ}}\overline{\Phi}\,dz+\int_{\rm{PZ}}(-\overline{\mathcal{B}})\,dz.$
(10)

Eqn. 10 is arranged so that the (positive) buoyant engine of convection is on
the left-hand side, and the (positive) sinks of work are on the RHS. If
viscous dissipation in the CZ does not balance the buoyant generation of
energy in the CZ, the kinetic energy of the convective flows grows, resulting
in a penetrative region. This region grows with time until Eqn. 10 is
satisfied. We see that the viscous dissipation and buoyant deceleration felt
by flows in the PZ determine its size. We now define

$f\equiv\frac{\int_{\rm{CZ}}\overline{\Phi}\,dz}{\int_{\rm{CZ}}\overline{\mathcal{B}}\,dz},$
(11)

the measurable fraction of the buoyant engine consumed by CZ dissipation. Eqn.
10 can then be rewritten as

$\frac{\int_{\rm{PZ}}(-\overline{\mathcal{B}})\,dz}{\int_{\rm{CZ}}\overline{\mathcal{B}}\,dz}+\frac{\int_{\rm{PZ}}\overline{\Phi}\,dz}{\int_{\rm{CZ}}\overline{\mathcal{B}}\,dz}=(1-f).$
(12)

We will measure and report the values of $f$ achieved in our simulations in
this work. Eqn. 12 provides two limits on a hypothetical PZ:

  1. 1.

In the limit that $f\rightarrow 0$, viscous dissipation is inefficient.
Reasonably if we also assume that
$\int_{\rm{PZ}}\overline{\Phi}\,dz\rightarrow 0$, Eqn. 12 states that the PZ
must be so large that its negative buoyant work is equal in magnitude to the
positive buoyant work of the CZ. This is the integral constraint on the
maximum size of the PZ that Roxburgh (1989) derived.

  2. 2.

In the limit that $f\rightarrow 1$, viscous dissipation efficiently
counteracts the buoyancy work in the CZ. Per Eqn. 12, the positive-definite PZ
terms must approach zero and no PZ develops in this limit. This is
mathematically equivalent to standard boundary-driven convection experiments.

In general, we anticipate from the results of e.g., Currie & Browning (2017)
that $f$ is closer to 1 than 0, but its precise value must be measured from
simulations. Indeed, we find that $f\gg 0$ but $f<1$ in our simulations (see
e.g., Fig. 3, bottom panel111the bulk dynamics suggest by eye $f\sim 0.5$, but
due to e.g., the height dependence of $\overline{\mathcal{B}}$ in our
simulations we measure $f\approx 0.74$.). Our simulations produce typical
values of $f\sim 0.7$.

Assuming that a PZ of height $\delta_{\rm{p}}$ develops above a CZ of depth
$L_{\rm{CZ}}$, we model the PZ dissipation as

$\int_{\rm{PZ}}\overline{\Phi}\,dz=\xi\frac{\delta_{\rm{p}}}{L_{\rm{CZ}}}\int_{\rm{CZ}}\overline{\Phi}\,dz=\xi\delta_{\rm{p}}\Phi_{\rm{CZ}}.$
(13)

Here $\Phi_{\rm{CZ}}$ is the volume-averaged dissipation rate in the CZ and
$\xi$ is a measurable parameter in ${[0,1]}$ that describes the shape of the
dissipation profile as a function of height in the PZ. In words, we assume
that $\overline{\Phi}(z=L_{s})\approx\Phi_{\rm{CZ}}$ at the CZ-PZ boundary and
that $\overline{\Phi}$ decreases with height in the PZ. The shape of
$\overline{\Phi}$ determines $\xi$; a linear falloff gives $\xi=1/2$, a
quadratic falloff gives $\xi=2/3$, and $\xi=1$ assumes no falloff. With this
parameterization, and $\overline{\mathcal{B}}\propto\overline{F_{\rm{conv}}}$,
we rewrite Eqn. 12,

$-\frac{\int_{\rm{PZ}}\overline{F_{\rm{conv}}}\,dz}{\int_{\rm{CZ}}\overline{F_{\rm{conv}}}\,dz}+f\xi\frac{\delta_{\rm{p}}}{L_{\rm{CZ}}}=(1-f).$
(14)

The fundamental result of this theory is Eqn. 14, which is a parameterized and
generalized form of Roxburgh (1989)’s integral constraint. This equation is
also reminiscent of Zahn (1991)’s theory, and says that the size of a PZ is
set by the profile of $\nabla_{\rm{rad}}$ near the convective boundary. A
parameterization like Eqn. 14 can be implemented in stellar structure codes
and used to find the extent of penetration zones under the specification of
$f$ and $\xi$. We note that an implementation of Eqn. 14 likely requires an
_iterative_ solve, as the penetration zone depth ($\delta_{\rm{p}}$) and thus
the PZ integral of the flux, are not known _a-priori_. The parameters $f$ and
$\xi$ are measurables which can be constrained by direct numerical
simulations, and we will measure their values in this work. In general, we
expect that $f$ and $\xi$ should not change too drastically with other
simulation parameters.

In order to derive a specific prediction for the PZ height, one must specify
the vertical shape of $\overline{F_{\rm{conv}}}$. We will study two cases in
this work, laid out below. In both of these cases, we define a nondimensional
“Penetration Parameter” whose magnitude is set by the ratio of the convective
flux slightly above and below the Schwarzschild convective boundary $L_{s}$
(assuming $\nabla=\nabla_{\rm{ad}}$ in the CZ and PZ),

$\mathcal{P}\equiv-\frac{\overline{F_{\rm{conv}}}_{\rm{CZ}}}{\overline{F_{\rm{conv}}}_{\rm{PZ}}}.$
(15)

Since $F_{\rm{conv}}<0$ in the PZ, the sign of $\mathcal{P}$ is positive.
Intuitively, $\mathcal{P}$ describes which terms are important in Eqn. 12.
When $\mathcal{P}\ll 1$, the buoyancy term dominates in the PZ and dissipation
can be neglected there. When $\mathcal{P}\gg 1$, buoyancy is negligible and
dissipation constrains the size of the PZ. When $\mathcal{P}\sim 1$, both
terms matter. In this work, we have assumed that $\mathcal{P}$ and $\xi$ are
fully independent parameters. We make this choice because $\mathcal{P}$ can be
determined directly from a known conductivity profile or stratification,
whereas $\xi$ is a measurable of evolved nonlinear convective dynamics.
However, it is possible that there is an implicit relationship between these
parameters (as $\mathcal{P}$ increases, so too does the extent of the PZ,
which likely in turn modifies the value of $\xi$).

####  3.2.1 Case I: Discontinuous flux

We first consider a model which satisfies

$\overline{F_{\rm{conv}}}(z)=F_{\rm{cz}}\begin{cases}1&z\leq L_{s},\\\
-\mathcal{P}_{D}^{-1}&z>L_{s}\end{cases}.$ (16)

Here, $F_{\rm{cz}}$ is a constant value of flux carried in the convection zone
and $\mathcal{P}_{D}$ is the penetration parameter (subscript D for
discontinuous case). Plugging this functional form of the flux into Eqn. 14,
and integrating the CZ over a depth $L_{\rm{CZ}}$ below $L_{s}$ and the PZ
over a height $\delta_{\rm{p}}$ above $L_{s}$, we predict

$\frac{\delta_{\rm{p}}}{L_{\rm{CZ}}}=\mathcal{P}_{D}\frac{1-f}{1+\xi
f\mathcal{P}_{D}}.$ (17)

Assuming that $f$ and $\xi$ are weak functions of $\mathcal{P}_{D}$, we see
that, for small $\mathcal{P}_{D}$, the size of the penetration region is
linearly proportional to $\mathcal{P}_{D}$, but saturates as
$\mathcal{P}_{D}\rightarrow\infty$ due to dissipation. Intuitively, this
result makes sense: as $\mathcal{P}_{D}$ grows, the magnitude of
$\overline{F_{\rm{conv}}}$ and the deceleration caused by buoyancy in the PZ
shrink, resulting in larger penetrative regions (but this growth cannot extend
indefinitely).

####  3.2.2 Case II: Piecewise linear flux

We next examine a model where the derivative of $\overline{F_{\rm{conv}}}(z)$
may be discontinuous at the CZ-PZ boundary,

$\overline{F_{\rm{conv}}}(z)=\frac{\partial F_{\rm{rad}}}{\partial
z}\bigg{|}_{\rm{CZ}}\begin{cases}(L_{s}-z)&z\leq L_{s}\\\
-\mathcal{P}_{L}^{-1}(z-L_{s})&z>L_{s}\end{cases},$ (18)

where $(\partial F_{\rm{rad}}/\partial z)|_{\rm{CZ}}$ is a constant and
$\mathcal{P}_{L}$ is the penetration parameter (subscript L for linear case).
When $\mathcal{P}_{L}=1$, $\overline{F_{\rm{conv}}}$ is a linear profile that
crosses through zero at $z=L_{s}$. Solving Eqn. 14 with Eqn. 18 and
integrating over $L_{\rm{CZ}}$ in the CZ and $\delta_{\rm{p}}$ in the PZ, we
retrieve a quadratic equation. This equation has two solution branches, only
one of which corresponds to a positive value of $\delta_{\rm{p}}$. On that
branch, we find

$\frac{\delta_{\rm{p}}}{L_{\rm{CZ}}}=\sqrt{\mathcal{P}_{L}(1-f)}\,\,(\sqrt{\zeta^{2}+1}-\zeta),$
(19)

where $\zeta\equiv(\xi f/2)\sqrt{\mathcal{P}_{L}/(1-f)}$. We expect the
penetration height to be proportional to $\sqrt{\mathcal{P}_{L}}$ for small
values of $\mathcal{P}_{L}$, and to again saturate at large values of
$\mathcal{P}_{L}$ (as $\mathcal{P}_{L}\rightarrow\infty$, so too
$\zeta\rightarrow\infty$, and $(\sqrt{\zeta^{2}+1}-\zeta)\rightarrow 0$).

In this work, we will test Eqn. 14 through the predictions of Eqns. 17 and 19.
Our goals are to see if the predicted scalings with the penetration parameter
$\mathcal{P}$ are realized in simulations, and to measure the values of $f$
and $\xi$.

##  4 Simulation Details

We will now describe a set of simulations that test the predictions in Sec. 3.
While many simulations of convection interacting with radiative zones have
been performed by previous authors, ours differ in two crucial ways. First, we
construct our experiments so that $\mathcal{P}$ and $\mathcal{S}$ can be
varied separately by driving convection with internal heating, thus avoiding
strongly superadiabatic boundary layers where
$\nabla\rightarrow\nabla_{\rm{rad}}$. $\mathcal{P}$ is the “Penetration
Parameter,” defined in Eqn. 15, which compares the magnitude of the convective
flux in the CZ and PZ; $\mathcal{S}$ is the “stiffness,” defined in Eqn. 27,
and compares the buoyancy frequency in the stable radiative zone to the
convective frequency. We suspect that some past experiments have implicitly
set $\mathcal{P}\approx\mathcal{S}^{-1}$, which would result in negligible
penetration for high stiffness (see discussion following Eqn. 27). Second, as
we will show in Sec. 5, the development of penetrative zones is a slow process
and many prior studies did not evolve simulations for long enough to see these
regions grow and saturate.

Appealing to the Buckingham $\pi$ theorem (buckingham_1914), we count nine
fundamental input parameters in Eqns. 1-4: $\rho_{0}$, $\alpha g$, $L_{s}$,
$\nu$, $\chi$, $Q$, $\nabla_{\rm{ad}}$, $k_{\rm{CZ}}$, and $k_{\rm{RZ}}$.
There are four fundamental dimensions (mass, length, time, and temperature),
and so we are left with five independent prognostic parameters in setting up
our system. For two of these parameters, we will choose the freefall Reynolds
number and the Prandtl number, which are analagous to the Rayleigh and Prandtl
numbers in Rayleigh-Bénard convection. The remaining three parameters are
$\mathcal{S}$, $\mathcal{P}$, and an additional parameter $\mu$, which we will
hold constant and which sets the ratio between $\nabla_{\rm{rad}}$ and
$\nabla_{\rm{ad}}$ in the convection zone.

We nondimensionalize Eqns. 1-4 on the length scale of the Schwarzschild-
unstable convection zone $L_{s}$, the timescale of freefall across that
convection zone

$\tau_{\rm{ff}}=\left(\frac{L_{s}}{|\alpha|gQ_{0}}\right)^{1/3},$ (20)

and the temperature scale of the internal heating over that freefall time
$\Delta T$; mass is nondimensionalized so that the freefall ram pressure
$\rho_{0}(L_{s}/\tau_{\rm{ff}})^{2}=1$,

$\begin{split}&T^{*}=(\Delta T)T=Q_{0}\tau_{\rm{ff}}T,\qquad Q^{*}=Q_{0}Q,\\\
&\partial_{t^{*}}=\tau_{\rm{ff}}^{-1}\partial_{t},\qquad\qquad\qquad\,\,\,\boldsymbol{\nabla}^{*}=L_{s}^{-1}\boldsymbol{\nabla},\,\,\,\,\,\\\
&\boldsymbol{u}^{*}=u_{\rm{ff}}\boldsymbol{u}=\frac{L_{s}}{\tau_{\rm{ff}}}\boldsymbol{u},\qquad\qquad
p^{*}=\rho_{0}u_{\rm{ff}}^{2}\varpi,\\\
&k^{*}=(L_{s}^{2}\tau_{\rm{ff}}^{-1})k,\qquad\mathcal{R}=\frac{u_{\rm{ff}}L_{s}}{\nu},\qquad\mathrm{Pr}=\frac{\nu}{\chi}.\end{split}$
(21)

For convenience, here we define quantities with $*$ (e.g., $T^{*}$) as being
the “dimensionful” quantities of Eqns. 1-4. Henceforth, quantities without $*$
(e.g., $T$) are dimensionless. The dimensionless equations of motion are

$\displaystyle\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0$ (22)
$\displaystyle\partial_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}=-\boldsymbol{\nabla}\varpi+T\hat{z}+\mathcal{R}^{-1}\boldsymbol{\nabla}^{2}\boldsymbol{u}$
(23)
$\displaystyle\begin{split}\partial_{t}T+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}T+w\boldsymbol{\nabla}_{\rm{ad}}+\boldsymbol{\nabla}\boldsymbol{\cdot}[-k\boldsymbol{\nabla}\overline{T}]\qquad\qquad\\\
\qquad\qquad\qquad\qquad\qquad=(\mathrm{Pr}\mathcal{R})^{-1}\boldsymbol{\nabla}^{2}T^{\prime}+Q.\end{split}$
(24)

We construct a domain in the range $z\in[0,L_{z}]$ and choose $L_{z}\geq 2$ so
that the domain is at least twice as deep as the Schwarzschild-unstable
convection zone. We decompose the temperature field into a time-stationary
initial background profile and fluctuations,
$T(x,y,z,t)=T_{0}(z)+T_{1}(x,y,z,t)$. $T_{0}$ is constructed with
$\nabla=\nabla_{\rm{ad}}$ for $z\leq L_{s}$, and $\nabla=\nabla_{\rm{rad}}$
above $z>L_{s}$. We impose a fixed-flux boundary at the bottom of the box
($\partial_{z}T_{1}=0$ at $z=0$) and a fixed temperature boundary at the top
of the domain ($T_{1}=0$ at $z=L_{z}$). We generally impose impenetrable, no-
slip boundary conditions at the top and bottom of the box so that
$\boldsymbol{u}=0$ at $z=[0,L_{z}]$. For a select few simulations, we impose
stress-free instead of no-slip boundary conditions ($w=0$ and
$\partial_{z}u=\partial_{z}v=0$ at $z=[0,L_{z}]$).

We impose a constant internal heating which spans only part of the convection
zone,

$Q=\begin{cases}0&z<0.1\,\,\rm{or}\,\,z\geq 0.1+\Delta_{H},\\\
Q_{\rm{mag}}&0.1\leq z\leq 0.1+\Delta_{H}\end{cases}.$ (25)

The integrated flux through the system from heating is
$F_{H}=\int_{0}^{L_{z}}Q_{\rm{mag}}dz=Q_{\rm{mag}}\Delta_{H}$. Throughout this
work we choose $Q_{\rm{mag}}=1$ and $\Delta_{H}=0.2$ so $F_{H}=0.2$. We offset
this heating from the bottom boundary to $z=0.1$ to avoid heating within the
bottom impenetrable boundary layer where velocities go to zero and $k$ is
small; this prevents strong temperature gradients from establishing there.
Furthermore, since the conductivity is not zero at the bottom boundary, the
adiabatic temperature gradient there carries some flux. We specify the flux
using

$\mu\equiv\frac{F_{\rm{bot}}}{F_{H}}$ (26)

and we choose $\mu=10^{-3}$ so that most of the flux in the convection zone is
carried by the convection.

Throughout this paper, we assume that the convection zone is roughly
adiabatically stratified. We therefore define a dynamical measure of the
stiffness, rather than one based on e.g., the superadiabaticity of
$\nabla_{\rm{rad}}$ in the convection zone. The average convective velocity
depends on the magnitude of the convective flux,
$\langle\boldsymbol{|u|}\rangle\approx
F_{H}^{1/3}=(Q_{\rm{mag}}\Delta_{H})^{1/3}$. The characteristic convective
frequency is ${f_{\rm{conv}}=\langle\boldsymbol{|u|}\rangle/L_{s}}$.
Empirically we find that for our choice of parameters,
$\langle\boldsymbol{|u|}\rangle\approx 1$, so going forward we define
$f_{\rm{conv}}=1$. The stiffness is defined,

$\mathcal{S}\equiv\frac{N^{2}}{f_{\rm{conv}}^{2}}=N^{2},$ (27)

where $N^{2}$ is the Brunt-Väisälä frequency in the radiative zone. In our
nondimensionalization, $N^{2}=\nabla_{\rm{ad}}-\nabla_{\rm{rad}}$ in the
radiative zone. We use $\mathcal{S}$ as a control parameter.

In many prior studies, the stiffness has been set by the ratio of the
subadiabaticity of $\nabla_{\rm{rad}}$ in the RZ to the superadiabaticity of
$\nabla_{\rm{rad}}$ in the CZ,

$\tilde{\mathcal{S}}=\frac{|\nabla_{\rm{rad}}-\nabla_{\rm{ad}}|_{\rm{RZ}}}{|\nabla_{\rm{rad}}-\nabla_{\rm{ad}}|_{\rm{CZ}}}=\frac{N^{2}}{|\nabla_{\rm{rad}}-\nabla_{\rm{ad}}|_{\rm{CZ}}}.$
(28)

In those studies, $\tilde{\mathcal{S}}$ primarily describes the stratification
of the initial state, but it also describes the stratification in
superadiabatic boundary layers which drive convection. In this work, we
maintain a nearly adiabatic convection zone without strongly superadiabatic
regions by driving convection with an internal heating function which is
offset from the lower boundary.

Previous work has not defined $\mathcal{P}$, but its definition in our current
study should apply to previous studies,

$\mathcal{P}=-\frac{k_{\rm{CZ}}(\nabla_{\rm{rad}}-\nabla_{\rm{ad}})_{\rm{CZ}}}{k_{\rm{RZ}}(\nabla_{\rm{rad}}-\nabla_{\rm{ad}})_{\rm{RZ}}}.$
(29)

We note that $\mathcal{P}$ can be related to $\mathcal{S}$ and
$\tilde{\mathcal{S}}$,
$\mathcal{P}=(k_{\rm{CZ}}/k_{\rm{RZ}})\tilde{\mathcal{S}}^{-1}=(k_{\rm{CZ}}/k_{\rm{RZ}})(\nabla_{\rm{rad}}-\nabla_{\rm{ad}})_{\rm{CZ}}\mathcal{S}^{-1}$.
Our use of internal heating to decouple convective perturbations from
$\nabla_{\rm{rad}}$ in the CZ allows us to separately specify these
nondimensional parameters. The distinction between $\mathcal{S}$ and
$\mathcal{P}$ is perhaps clearer in the language of stellar evolution, where
$\mathcal{S}$ is roughly the inverse square Mach number of the convection
while $\mathcal{P}$ is set by the ratio of $\nabla_{\rm{rad}}$ and
$\nabla_{\rm{ad}}$.

Aside from $\mathcal{S}$, $\mathcal{P}$, and $\mu$, the two remaining control
parameters $\mathcal{R}$ and $\mathrm{Pr}$ determine the properties of the
turbulence. The value of $\mathcal{R}$ corresponds to the value of the
Reynolds number $\rm{Re}=\mathcal{R}|\boldsymbol{u}|$, and we will vary
$\mathcal{R}$. Astrophysical convection exists in the limit of $\mathrm{Pr}\ll
1$ (Garaud, 2021); in this work we choose a modest value of $\mathrm{Pr}=0.5$
which slightly separates the thermal and viscous scales while still allowing
us to achieve convection with large Reynolds and Péclet numbers.

We now describe the two types of simulations conducted in this work (Case I
and Case II). We provide Fig. 4 to visualize the portion of the parameter
space that we have studied. We denote two “landmark cases” using a purple box
(Case I landmark) and an orange box (Case II landmark). These landmark cases
will be mentioned throughout this work.

###  4.1 Case I: Discontinuous flux

Most of the simulations in this paper have a discontinuous convective flux at
the Schwarzschild convective boundary. We achieve this by constructing a
discontinuous radiative conductivity,

$k=\begin{cases}k_{\rm{CZ}}&z<1\\\ k_{\rm{RZ}}&z\geq 1\end{cases},$ (30)

where CZ refers to the convection zone and RZ refers to the radiative zone
(some of which will be occupied by the penetrative zone PZ). Using
$\mathcal{S}$ and $\mathcal{P}_{D}$ as inputs and specifying the radiative
flux at the bottom boundary and in the RZ defines this system222 We solve the
system of equations
$\mathcal{S}=(\nabla_{\rm{ad}}-\nabla_{\rm{rad}})/f_{\rm{conv}}^{2}$,
$\mathcal{P}_{D}=F_{H}/(k_{\rm{RZ}}[\nabla_{\rm{ad}}-\nabla_{\rm{rad}}])$,
$F_{\rm{bot}}=k_{\rm{CZ}}\nabla_{\rm{ad}}$, and
${F_{\rm{bot}}+F_{H}=k_{\rm{RZ}}\nabla_{\rm{rad}}}$ to arrive at Eqns. 31. ,

$\begin{split}&k_{\rm{RZ}}=\frac{F_{H}}{f_{\rm{conv}}^{2}\mathcal{S}\mathcal{P}_{D}},\\\
&k_{\rm{CZ}}=k_{\rm{RZ}}\frac{\mu}{1+\mu+\mathcal{P}_{D}^{-1}},\\\
&\nabla_{\rm{ad}}=f_{\rm{conv}}^{2}\mathcal{S}\mathcal{P}_{D}(1+\mu+\mathcal{P}_{D}^{-1}),\\\
&\nabla_{\rm{rad}}=\nabla_{\rm{ad}}-f_{\rm{conv}}^{2}\mathcal{S}.\end{split}$
(31)

Figure 4:  Each marker denotes a simulation conducted in this work in the
$\mathcal{R}-\mathcal{P}$ parameter space at $\mathcal{S}=10^{3}$. Purple
circles represent Case I (Sec. 4.1) simulations and orange diamonds represent
Case II (Sec. 4.2) simulations; empty circular markers have stress-free (SF)
boundary conditions and all other simulations have no-slip boundaries. The
green “x” at $\mathcal{P}=4$ and $\mathcal{R}=400$ denotes the location in
$\mathcal{R}-\mathcal{P}$ parameter space where we vary $\mathcal{S}$ in
select Case I simulations. Boxes denote the two “landmark” simulations. The
landmark Case I simulation has $\mathcal{R}=400$ and $\mathcal{P}_{D}=4$. The
landmark Case II simulation has $\mathcal{R}=800$ and $\mathcal{P}_{L}=4$.
Both landmark simulations have $\mathcal{S}=10^{3}$ and no-slip boundary
conditions.

We study a sweep through each of the ($\mathcal{P}_{D}$, $\mathcal{S}$,
$\mathcal{R}$) parameter spaces while holding all other parameters constant
(see Fig. 4). We study an additional sweep through $\mathcal{R}$ parameter
space using stress-free boundaries to compare to our no-slip cases. According
to Eqn. 17, we expect $\delta_{\rm{p}}\propto\mathcal{P}_{D}$.

###  4.2 Case II: Piecewise linear flux

We also study simulations where the flux’s gradient may be discontinuous at
the Schwarzschild convective boundary. We achieve this by constructing a
radiative conductivity with a piecewise discontinuous gradient,

$\partial_{z}k=\partial_{z}k_{0}\begin{cases}1&z<1\\\
\mathcal{P}_{L}^{-1}&z\geq 1\end{cases}$ (32)

Since $k$ varies with height, formally the values of $\mathcal{S}$ and
$\mathcal{P}$ also vary with height; we specify their values at $z=2$. By this
choice333 We solve the system of equations where
${F_{\rm{bot}}=k_{\rm{bot}}\nabla_{\rm{ad}}}$,
${F_{\rm{bot}}+F_{H}=k_{\rm{ad}}\nabla_{\rm{ad}}}$,
${k_{\rm{ad}}=k_{\rm{bot}}+\partial_{z}k_{0}L_{s}}$,
${\mathcal{S}=(\nabla_{\rm{ad}}-\nabla_{\mathrm{rad},z=2L_{s}})/f_{\rm{conv}}^{2}}$,
and ${\nabla_{\rm{rad}}=F_{\rm{tot}}/k(z)}$. , we require

$\partial_{z}k_{0}=\frac{F_{H}}{f_{\rm{conv}}^{2}L_{s}\mathcal{S}\psi},\,\,k_{b}=\frac{F_{H}\mu}{f_{\rm{conv}}^{2}\mathcal{S}\psi},\,\,\nabla_{\rm{ad}}=f_{\rm{conv}}^{2}\mathcal{S}\psi,$
(33)

where $\psi\equiv 1+\mathcal{P}_{L}(1+\mu)$. We will study one sweep through
$\mathcal{P}_{L}$ space at fixed $\mathcal{R}$ and $\mathcal{S}$ (see Fig. 4).
According to Eqn. 19, we expect $\delta_{p}\propto\mathcal{P}_{L}^{1/2}$.

###  4.3 Numerics

We time-evolve equations 22-24 using the Dedalus pseudospectral solver (Burns
et al., 2020)444we use commit efb13bd; the closest stable release to this
commit is v2.2006. using timestepper SBDF2 (Wang & Ruuth, 2008) and safety
factor 0.35. All fields are represented as spectral expansions of $n_{z}$
Chebyshev coefficients in the vertical ($z$) direction and as
($n_{x}$,$n_{y}$) Fourier coefficients in the horizontal ($x$,$y$) directions;
our domains are therefore horizontally periodic. We use a domain aspect ratio
of two so that $x\in[0,L_{x}]$ and $y\in[0,L_{y}]$ with $L_{x}=L_{y}=2L_{z}$.
To avoid aliasing errors, we use the 3/2-dealiasing rule in all directions. To
start our simulations, we add random noise temperature perturbations with a
magnitude of $10^{-3}$ to a background temperature profile $\overline{T}$; we
discuss the choice of $\overline{T}$ in appendix A. In some simulations we
start with $\overline{T}=T_{0}$, described above, and in others we impose an
established penetrative zone in the initial state $\overline{T}$ according to
Eqn. A1.

Spectral methods with finite coefficient expansions cannot capture true
discontinuities. In order to approximate discontinuous functions such as Eqns.
25, 30, and 32, we must use smooth transitions. We therefore define a smooth
Heaviside step function,

$H(z;z_{0},d_{w})=\frac{1}{2}\left(1+\mathrm{erf}\left[\frac{z-z_{0}}{d_{w}}\right]\right).$
(34)

where erf is the error function. In the limit that $d_{w}\rightarrow 0$, this
function behaves identically to the classical Heaviside function centered at
$z_{0}$. For Eqn. 25 and Eqn. 32, we use $d_{w}=0.02$; while for Eqn. 30 we
use $d_{w}=0.075$. In all other cases, we use $d_{w}=0.05$.

A table describing all of the simulations presented in this work can be found
in Appendix C. We produce the figures in this paper using matplotlib (Hunter,
2007; Caswell et al., 2021). All of the Python scripts used to run the
simulations in this paper and to create the figures in this paper are publicly
available in a git
repository555https://github.com/evanhanders/convective_penetration_paper, and
in a Zenodo repository (Anders et al., 2021).

###  4.4 Penetration height measurements

In our evolved simulations, the penetrative region has a nearly adiabatic
stratification $\nabla\approx\nabla_{\rm{ad}}$. To characterize the height of
the penetrative region, we measure how drastically $\nabla$ has departed from
$\nabla_{\rm{ad}}$. We define the difference between the adiabatic and
radiative gradient,

$\Delta\equiv\nabla_{\rm{ad}}-\nabla_{\rm{rad}}(z).$ (35)

We measure penetration heights in terms of “departure points,” or heights at
which the realized temperature gradient $\nabla$ has evolved away from the
adiabatic $\nabla_{\rm{ad}}$ by some fraction $h<1$ of $\Delta$. Specifically,

$L_{s}+\delta_{h}=\mathrm{max}(z)\,\,\mid\,\,\nabla>(\nabla_{\rm{ad}}-h\,\Delta).$
(36)

In this work, we measure the 10% ($\delta_{0.1}$, $h=0.1$), 50%
($\delta_{0.5}$, $h=0.5$), and 90% ($\delta_{0.9}$, $h=0.9$) departure points.
Using Zahn (1991)’s terminology, $\delta_{0.5}$ is the mean value of the top
of the PZ while $\delta_{0.9}-\delta_{0.1}$ represents the width of the PZ-RZ
boundary layer. We find that these measurements based on the (slowly-evolving)
thermodynamic profile provide a robust and straightforward measurement of
penetration height (for a discussion of alternate measurement choices, see
Pratt et al., 2017).

Figure 5:  Temperature anomalies in vertical slices through the landmark
simulations. (top) Case I landmark ($\mathcal{R}=400$, $\mathcal{P}_{D}=4$,
$\mathcal{S}=10^{3}$) and (bottom) Case II landmark ($\mathcal{R}=800$,
$\mathcal{P}_{L}=4$, $\mathcal{S}=10^{3}$). The temporally- and volume-
averaged Reynolds number in the CZ is $\rm{Re}\sim 250$ in the top panel and
$\rm{Re}\sim 350$ in the bottom panel. A dashed horizontal line denotes the
Schwarzschild convective boundary. A solid line denotes the boundary between
the penetrative and radiative zones. The Case II simulation has an additional
Schwarzschild boundary near the bottom of the domain due to the conductivity
linearly increasing below the internal heating layer. As in Fig. 1,
temperature anomalies have different signs in the bulk CZ and PZ.  Figure 6:
Time evolution of the landmark Case I simulation ($\mathcal{R}=400$,
$\mathcal{P}_{D}=4$, $\mathcal{S}=10^{3}$). In the top left panel, we plot the
PZ height $\delta_{0.5}$ vs. time. Also shown are the time evolution of $f$
(middle left panel, defined in Eqn. 11) and $\xi$ (bottom left panel, defined
in Eqn. 13). Thin horizontal lines denote the equilibrium values of each
trace. Vertical lines denote when each trace first converges to within 1% of
its equilibrium value. (right panel) The vertical profile of
$\nabla/\nabla_{\rm{ad}}$ is plotted against height at regular time intervals.
The line color denotes the time, following the time traces in the left panels.
A horizontal dashed grey line denotes the constant value of
$\nabla_{\rm{ad}}$. The solid grey curve denotes the profile of
$\nabla_{\rm{rad}}$. The location of the Schwarzschild convective boundary is
displayed as a vertical dashed black line. The top-of-PZ departure points
(Eqn. 36) are plotted over the profile evolution ($\delta_{0.1}$ and
$\delta_{0.9}$ as red lines, $\delta_{0.5}$ as a black line).

##  5 Results

We now describe the results of the 3D dynamical simulations described in the
previous section. Fig. 1 displays the dynamics in one of these simulations.
While we will briefly examine dynamics here, our primary goal in this section
is to quantitatively compare our simulations to the theory of Sec. 3 using
temporally averaged measures.

###  5.1 Dynamics

In Fig. 5 we display snapshots of the temperature anomalies in the two
“landmark” simulations denoted by boxes in Fig. 4. We display the temperature
anomaly in the top panel of the Case I simulation with $\mathcal{R}=400$,
$\mathcal{P}_{D}=4$, and $\mathcal{S}=10^{3}$; this simulation is included in
all three of our parameter space sweeps and represents the point where our
($\mathcal{R},\mathcal{P},\mathcal{S}$) cuts converge in Fig. 4. We display
the temperature anomaly in the bottom panel of the Case II simulation with
$\mathcal{R}=800$, $\mathcal{P}_{L}=4$, and $\mathcal{S}=10^{3}$. The bulk
Reynolds number in the convection zones of these simulations are (top)
$\rm{Re}\sim 250$ and (bottom) $\rm{Re}\sim 350$. Thus, these simulations are
less turbulent than the simulation in Fig. 1 (bulk Re $\sim$ 5000). Aside from
the degree of turbulence, the dynamics are very similar in Figs. 1 & 5. In
particular, we observe that relatively hot plumes in the CZ turn into
relatively cold plumes in the PZ (as they cross the dashed horizontal lines),
and relatively hot regions in the PZ lie above relatively cold regions in the
CZ. Convective plumes extend through the penetrative region and impact the
stable radiative zone (above the solid horizontal line). The convective
motions excite waves at a shallow angle above the stiff radiative-convective
boundary. We note that the Case II simulation has an additional temperature
inversion at the base of the simulation. Case II simulations have a linearly
increasing conductivity $k$ in the convection zone, so there is formally a
small penetrative region where
$\nabla\approx\nabla_{\rm{ad}}>\nabla_{\rm{rad}}$ at the base of the domain
below the internal heating layer (lower dotted line in bottom panel of Fig.
5).

While the landmark simulations in Fig. 5 are not as turbulent as the dynamics
in Fig. 1, they are sufficiently nonlinear to be interesting. Importantly,
these simulations develop large penetration zones, and can be evolved for tens
of thousands of convective overturn times. As we will demonstrate in the next
section, the formation timescale of penetrative zones can take tens of
thousands of convective overturn times.

Figure 7:  Simulation results vs. $\mathcal{P}$ for both Case I (left panels;
solid purple circles in Fig. 4) and Case II (right panels; solid orange
diamonds in Fig. 4). Boxed data points denote landmark simulations from Fig.
4. The top panels show the penetration height according to Eqn. 36. The Case I
penetration heights (upper left) vary linearly with $\mathcal{P}$, in line
with the prediction of Eqn. 17. The Case II penetration heights (upper right)
vary like $\sqrt{\mathcal{P}}$, in line with the prediction of Eqn. 19. In the
middle panels, we measure $f$ according to Eqn. 11. We find values of
$f\in[0.6,0.9]$, and changes in $f$ are secondary to changes in $\mathcal{P}$
for determining penetration heights. In the bottom panels, we measure $\xi$
according to Eqn. 13. We find characteristic values of $\xi\in[0.5,0.75]$,
suggesting that the falloff of the $\overline{\Phi}$ in the PZ is well
described by a linear function (at high $\mathcal{P}$ when $\xi\approx 1/2$),
or by a cubic function (at low $\mathcal{P}$ when $\xi\approx 3/4$).

###  5.2 Qualitative description of simulation evolution

In Fig. 6, we show the time evolution of the landmark Case I simulation
($\mathcal{R}=400$, $\mathcal{S}=10^{3}$, and $\mathcal{P}_{D}=4$) whose
initial temperature profile sets $\nabla=\nabla_{\rm{ad}}$ in the convection
zone ($z\lesssim 1$) and $\nabla=\nabla_{\rm{rad}}$ in the radiative zone
($z\gtrsim 1$). In the top left panel, we display the height of the
penetrative region $\delta_{\rm{0.5}}$ vs. time. This region initially grows
quickly over hundreds of freefall times, but this evolution slows down;
reaching the final equilibrium takes tens of thousands of freefall times. The
evolution of the other parameters in our theory ($f$, $\xi$) are shown in the
middle and bottom left panels of Fig. 6. We plot the rolling mean, averaged
over 200 freefall time units. We see that the values of $f$ and $\xi$ reach
their final values ($f\approx 0.67$, $\xi\approx 0.58$) faster than the
penetration zone evolves to its full height. We quantify this fast evolution
by plotting vertical lines in each of the left three panels corresponding to
the first time at which the rolling average converges to within 1% of its
equilibrated value. The equilibrated value is averaged over the final 1000
freefall times of the simulation and plotted as a grey horizontal line. The
evolved value of $f$ indicates that roughly 2/3 of the buoyancy driving is
dissipated in the bulk CZ, so that 1/3 is available for PZ dissipation and
negative buoyancy work. The evolved value of $\xi$ indicates that the shape of
dissipation in the PZ is slightly steeper than linear.

Figure 8:  Case I simulations results vs. $\mathcal{S}$ at $\mathcal{R}=400$,
$\mathcal{P}=4$. Boxed data points denote the landmark simulation from Fig. 4.
(Left panel) Penetration heights vs. $\mathcal{S}$. While $\delta_{0.1}$ and
$\delta_{0.9}$ show some variation, the mean penetration height
($\delta_{0.5}$) is roughly constant. (Right panel) The width of the thermal
transition layer ($\delta_{0.9}-\delta_{0.1}$) vs. $\mathcal{S}$. We roughly
observe a $\mathcal{S}^{-1/2}$ scaling.

In the right panel of Fig. 6, we plot the profile of $\nabla/\nabla_{\rm{ad}}$
in our simulation at regular time intervals, where the color of the profile
corresponds to time, as in the left panels. $\nabla_{\rm{ad}}$ is plotted as a
dashed horizontal line while $\nabla_{\rm{rad}}$ is plotted as a grey solid
line which decreases with height around $z\approx 1$ and satures to a constant
above $z\gtrsim 1.1$. The location of the Schwarzschild boundary, $L_{s}$, is
overplotted as a black vertical dashed line. We note that the Schwarzschild
boundary does not move over the course of our simulation, so the extention of
the convection zone past this point is true penetration and not the result of
entrainment-induced changes in the Schwarzschild (or Ledoux) convective
boundaries. The traces of $\delta_{0.1}$ and $\delta_{0.9}$ are overplotted as
red lines while that of $\delta_{0.5}$ is plotted as a black line. We see that
the fast initial evolution establishes a sizeable PZ (denoted by purple
$\nabla$ profiles), but its final equilibration takes much longer (indicated
by the separation between the purple, green, and yellow profiles decreasing
over time).

This long evolution is computationally expensive; for this modest simulation
(256x642 coefficients), this evolution takes roughly 24 days on 1024 cores for
a total of $\sim$600,000 cpu-hours. It is not feasible to perform simulations
of this length for a full parameter space study, and so we accelerate the
evolution of most of the simulations in this work. To do so, we take advantage
of the nearly monotonic nature of the evolution of $\delta_{\rm{p}}$ vs. time
displayed in Fig. 6. We measure the instantaneous values of
$(\delta_{0.1},\delta_{0.5},\delta_{0.9})$, as well as their instantaneous
time derivatives. Using these values, we take a large “time step” forward to
evolve $\delta_{\rm{p}}$. While doing so, we preserve the width of the
transition from the PZ to the RZ, and we also adjust the solution so that
$\nabla=\nabla_{\rm{rad}}$ in the RZ, effectively equilibrating the RZ
instantaneously. In other words, we reinitialize the simulation’s temperature
profile with a better guess at its evolved state based on its current
dynamical evolution. For details on how this procedure is carried out, see
Appendix A.

###  5.3 Dependence on $\mathcal{P}$

We find that the height of the penetration zone is strongly dependent on
$\mathcal{P}$. In the upper two panels of Fig. 7, we plot the penetration
height ($\delta_{0.1}$, $\delta_{0.5}$, $\delta_{0.9}$ from Eq. 36) from Case
I simulations (discontinuous $k$, upper left) and Case II simulations
(discontinuous $\partial_{z}k$, upper right). The fixed values of
$\mathcal{R}$ and $\mathcal{S}$ are shown above these panels. We find that the
leading-order $\mathcal{P}$ scaling predictions of Eqns. 17 & 19 describe the
data well at intermediate values of $\mathcal{P}$ (orange lines). At small
values of $\mathcal{P}$ we see somewhat weaker scalings than these
predictions, because the profiles of $k$ and $\partial_{z}k$ are not truly
discontinuous but jump from one value in the CZ to another in the RZ over a
finite width (see e.g., the $\nabla_{\rm{rad}}$ profile in Figs. 2 & 6 and
Sec. 4.3). At large values of $\mathcal{P}$, the penetration height falls off
of these predicted scaling laws. In this regime, dissipation dominates over
buoyancy in the PZ, so the PZ height saturates.

The middle and bottom panels of Fig. 7 demonstrate that that $f$ and $\xi$ are
to leading order constant with $\mathcal{P}$. However, we find that $f$ has
slightly smaller values in the Case I simulations (left) than in the Case II
simulations (right). We measure characteristic values of $f\in[0.6,0.9]$,
signifying that 60-90% of the buoyant work is balanced by dissipation in the
convection zone, depending on the simulation. We note a weak trend where $f$
decreases as $\mathcal{P}$ increases. As $\mathcal{P}$ increases, we find that
CZ velocities decrease, leading to a decrease in the dissipation rate. When
$\mathcal{P}$ is small, the PZ-RZ boundary (which acts like a wall, left panel
of Fig. 1) efficiently deflects convective velocities sideways resulting in
increased bulk-CZ velocities. As $\mathcal{P}$ grows, the velocities have
access to an extended PZ in which to buoyantly decelerate before deflection,
resulting in slightly lower bulk velocities. A similar trend of $\xi$
decreasing as $\mathcal{P}$ increases can be seen. Recall that smaller values
of $\xi$ indicate the dissipative dynamics are rather different in the PZ and
CZ. As the size of the PZ grows, the dynamical structures of the PZ shift from
what is found in the CZ, and so $\xi$ shrinks.

###  5.4 Dependence on $\mathcal{S}$

We find that the height of the penetration zone is weakly dependent on
$\mathcal{S}$. In the left panel of Fig. 8, we plot the penetration height of
a few Case I simulations with $\mathcal{P}_{D}=4$ and $\mathcal{R}=400$ but
with different values of $\mathcal{S}$. The mean penetration height
$\delta_{0.5}$ varies only weakly with changing $\mathcal{S}$, but that the
values of $\delta_{0.1}$ and $\delta_{0.9}$ vary more strongly. The PZ-RZ
boundary layer in which $\nabla$ changes from $\nabla_{\rm{ad}}$ to
$\nabla_{\rm{rad}}$ becomes narrower as $\mathcal{S}$ increases. To quantify
this effect, we plot $\delta_{0.9}-\delta_{0.1}$ in the righthand panel of
Fig. 8. We find that the width of this region varies roughly according to a
$\mathcal{S}^{-1/2}$ scaling law, reminiscent of the pure-overshoot law
described by Korre et al. (2019).

Note that if the enstrophy, $\omega^{2}$ in the convection zone exceeds the
value of the square buoyancy frequency $N^{2}$ in the radiative zone, the
gravity waves in the RZ become nonlinear. We therefore restrict the
simulations in this study to relatively large666These values are large for
nonlinear simulations, but modest compared to astrophysical values. While
there is observational uncertainty about the magnitude of deep convective
velocities in the Sun, in the MESA model presented in Sct. 6,
$f_{\rm{conv}}\approx 10^{-6}$ s-1 and $N\approx 10^{-3}$ s-1, so
$\mathcal{S}\approx 10^{6}$. values of $10^{2}\leq\mathcal{S}<10^{4}$ in order
to ensure $N^{2}>\omega^{2}$ even in our highest enstrophy simulations.

###  5.5 Dependence on $\mathcal{R}$

We find that the height of the penetration zone is weakly dependent on
$\mathcal{R}$. In the upper left panel of Fig. 9, we find a logarithmic
decrease in the penetration height with the Reynolds number. In order to
understand how this could happen at fixed $\mathcal{P}$, we also plot the
output values of $f$ (upper middle) and $\xi$ (upper right). We find that $f$
increases with increasing $\mathcal{R}$, but is perhaps leveling off as
$\mathcal{R}$ becomes large. We find that $\xi$ does not increase strongly
with $\mathcal{R}$ except for in the case of laminar simulations with
$\mathcal{R}<200$. Eqn. 17 predicts that $\delta_{\rm{p}}$ should change at
fixed $\mathcal{P}$ and $\xi$ if $f$ is changing. In the bottom left panel, we
show that the change in $\delta_{\rm{p}}$ is due to this change in $f$. We
find that this is true both for simulations with stress-free dynamical
boundary conditions (open symbols, SF) and for no-slip conditions (closed
symbols, NS).

Figure 9:  (Upper left panel) Penetration heights vs. $\mathcal{R}$ for Case I
simulations (vertical cuts in Fig. 4). Empty markers denote stress-free
boundaries (SF) and filled markers denote no-slip boundaries (NS). In both
cases, we see a roughly logarithmic decrease of $\delta_{\rm{p}}$ vs.
$\mathcal{R}$. (Upper middle panel) $f$ increases with $\mathcal{R}$. (Upper
right panel) $\xi$ does not change appreciably with $\mathcal{R}$ for
turbulent simulations with $\mathcal{R}\geq 200$. (Lower left panel) There is
a strong correlation between $\delta_{0.5}$ and $f$, agreeing with our
theoretical model of Eqn. 17. (Lower middle panel) Changes in $f$ are roughly
linearly proportional to the depth of the viscous boundary layer,
$\ell_{\nu}$, at the bottom of the domain. (Lower right panel) $\ell_{\nu}$
follows a well-known convective scaling law, so $\delta_{0.5}$ and $f$ should
saturate as $\mathcal{R}\rightarrow\infty$ and $\ell_{\nu}\rightarrow 0$.
Boxed data points denote the landmark simulation from Fig. 4.

We now examine why $f$ increases as $\mathcal{R}$ increases. In the SF
simulations, within the CZ, we can reasonably approximate $\overline{\Phi}$ as
a constant $\Phi_{\rm{CZ}}$ in the bulk and zero within the viscous boundary
layer,

$\overline{\Phi}(z)=\begin{cases}\Phi_{\rm{CZ}}&z>\ell_{\nu}\\\
0&z\leq\ell_{\nu}\end{cases},$ (37)

where $\ell_{\nu}$ is the viscous boundary layer depth. We have visualized a
NS dissipation profile in the bottom panel of Fig. 3; SF simulations look
similar in the bulk, but drop towards zero at the bottom boundary rather than
reaching a maximum. Then, we have

$\int_{\rm{CZ}}\overline{\Phi}\,dz\approx\Phi_{\rm{CZ}}\left(L_{s}-\ell_{\nu}\right),$
(38)

and so per Eqn. 11,

$f=f_{\infty}\left(1-\frac{\ell_{\nu}}{L_{s}}\right),$ (39)

where $f_{\infty}$ is the expected value of $f$ at $\mathcal{R}=\infty$ when
$\ell_{\nu}=0$. So we see that the CZ dissipation and therefore $f$ vary
linearly with $\ell_{\nu}$.

In the bottom middle panel of Fig. 9, we find that Eqn. 39 with
$f_{\infty}=0.755$ captures the high-$\mathcal{R}$ behavior. To measure
$\ell_{\nu}$, we first measure the height of the extremum of the viscous
portion of the kinetic energy flux $\overline{\mathcal{F}}$ near the boundary,
and take $\ell_{\nu}$ to be the twice that height. We find that Eqn. 39 is a
slightly better description for the SF simulations than the NS simulations; NS
simulations have maximized dissipation in the boundary layer, and therefore
Eqn. 37 is a poor model for $z\leq\ell_{\nu}$. In the bottom right panel of
Fig. 9, we demonstrate that the depth of the viscous boundary layer follows
classical scaling laws from Rayleigh-Bénard convection777 If you assume the
Nusselt Number dependence on the Rayleigh number is throttled by the
boundaries, Nu $\propto$ Ra1/3 (as is frequently measured), and the Reynolds
number is Re $\propto$ Ra1/2, you retrieve Nu $\propto$ Re2/3. The Nusselt
number generally varies like the inverse of the boundary layer depth, Nu
$\propto$ $\ell^{-1}$, and so we expect $\ell_{\nu}\propto\mathcal{R}^{-2/3}$.
(Ahlers et al., 2009; Goluskin, 2016). Combining these trends, we expect

$f=f_{\infty}(1-C\mathcal{R}^{-2/3})$ (40)

for a constant $C$. Thus as $\mathcal{R}\rightarrow\infty$, $f\rightarrow
f_{\infty}$.

We use the fitted function of $f$ from the bottom middle panel, along with
Eqn. 17, to estimate $\delta_{0.5}$ in the bottom left panel. We need to
multiply this equation by a factor of 0.9, which accounts for some differences
between the simulations and the idealized “discontinuous flux” theoretical
model. First, due to internal heating and the finite width of the conductivity
transition around the Schwarzschild boundary, the convective flux is not truly
constant through the full depth of the CZ. Thus, we expect $L_{\rm{CZ}}$ in
Eqn. 17 to be smaller than 1. Furthermore, the theory is derived in the limit
of an instantaneous transition from $\nabla_{\rm{ad}}$ to $\nabla_{\rm{rad}}$
where $\delta_{0.1}=\delta_{0.5}=\delta_{0.9}$; our simulations have a finite
transition width. Despite these subtle differences, we find good agreement.

Using $f_{\infty}=0.755$ we estimate that $\delta_{0.5}\approx 0.31$ for
$\mathcal{R}\rightarrow\infty$ and plot this as a horizontal orange line on
the upper left panel of Fig. 9. This value is coincidentally very near the
value of $\delta_{0.5}$ achieved in our highest-$\mathcal{R}$ simulations.
Unfortunately, we cannot probe more turbulent simulations. We can only run the
$\mathcal{R}=6.4\times 10^{3}$ simulation for a few hundred freefall times.
Our accuracy in measuring results from this simulation is limited by the long
evolutionary timescales of the simulation (see Fig. 6 for similar evolution in
a less turbulent, $\mathcal{R}=400$ case). Even accounting for our accelerated
evolutionary procedure, we can only be confident that the PZ heights of this
simulation are converged to within a few percent. Future work should aim to
better understand the trend of PZ height with turbulence. However, the
displayed relationships between $\delta_{\rm{p}}$ and $f$, $f$ and
$\ell_{\nu}$, and $\ell_{\nu}$ and $\mathcal{R}$ — all of which are effects we
largely understand — suggest that PZ heights should saturate at high
$\mathcal{R}$.

In summary, we find that $\delta_{\rm{p}}$ decreases as $\mathcal{R}$
increases. We find that these changes are caused by increases in $f$. In our
simulations, $f$ seems to have a linear relationship with the size of the
viscous boundary layer $\ell_{\nu}$. By measuring $f$ and $\ell_{\nu}$ in a
simulation, the value of $f_{\infty}$ can be found from Eqn. 39. Stellar
convection zones are not adjacent to hard walls888 Core convection zones have
no lower boundary due to geometry; flows pass through the singular point at
$r=0$. Convective shells in should be bounded both above and below by
penetrative regions. , so $f_{\infty}$ and the limit $\ell_{\nu}\rightarrow 0$
applies to stellar convection.

While we have examined a Case I simulation with $\mathcal{P}=4$ here, we
expect the simulation with $\mathcal{P}_{L}=1$ (a linear radiative
conductivity profile) to be the most representative of conditions near a
stellar convective boundary. In this simulation, we measure $\xi\approx 0.6$,
$f\approx 0.785$, $\ell_{\nu}\approx 0.08$, and $L_{s}=1$. Using Eqn. 39, we
estimate that

$f_{\infty}=0.86\qquad\rm{and}\qquad\xi=0.6$ (41)

are good first estimates for $f$ and $\xi$ when applying our theory of
penetrative convection to stellar models.

##  6 Testing our parameterization in a simple stellar model of the Sun

Figure 10:  (top) Profiles of $\nabla$ (green), $\nabla_{\rm{ad}}$ (purple),
and $\nabla_{\rm{rad}}$ (orange) in a 1 $M_{\odot}$ MESA stellar model with a
penetration zone. (bottom) Sound speed differences between the model shown in
the top panel and a standard (std) model run at identical parameters but
without a PZ. The addition of a PZ creates an acoustic glitch, raising the
sound speed by $\mathcal{O}$(2%) below the convection zone.

Our simulation results present a strong case for a flux- and dissipation-based
model of convective penetration, similar to those considered by Zahn (1991)
and Roxburgh (1989). In this section, we discuss a simple stellar model of the
Sun which we have created by implementing our parameterization into MESA (see
Appendix B). We of course note that the theory and 3D simulations in this work
do not include many of the complications of stellar convection like density
stratification, sphericity, rotation, magnetism, etc. We present this model as
a proof of concept and to inspire further work.

In order to implement our theory into MESA, we need to extend Eqn. 14 to
spherical geometry. To do so, we replace horizontal averages in Eqn. 9 with
integrals over latitude and longitude, and find that the relevant integral
constraint contains the convective luminosity,

$\int|\alpha|gL_{\rm{conv}}\,dr=\int_{V}\rho_{0}\Phi\,dV,$ (42)

where $L_{\rm{conv}}=4\pi\rho_{0}r^{2}\overline{F_{\rm{conv}}}$, $r$ is the
radial coordinate, and we write the RHS as a volume integral. We next define
$f$ in the same way as in Eqn. 11 and define $\xi$ similarly to Eqn. 13,

$\int_{\rm{PZ}}\rho_{0}\Phi\,dV=\xi\frac{V_{\rm{PZ}}}{V_{\rm{CZ}}}\int_{\rm{CZ}}\rho_{0}\Phi\,dV,$
(43)

where $V_{\rm{PZ}}$ and $V_{\rm{CZ}}$ are the volumes of the PZ and CZ
respectively. Eqn. 43 generalizes Eqn. 13 outside of the assumption of a
plane-parallel atmosphere. Thus Eqn. 14 in spherical geometry is

$-\frac{\int_{\rm{PZ}}L_{\rm{conv}}\,dr}{\int_{\rm{CZ}}L_{\rm{conv}}\,dr}+f\xi\frac{V_{\rm{PZ}}}{V_{\rm{CZ}}}=(1-f),$
(44)

We implemented Eqn. 44 in MESA (see Appendix B for details) and evolved a
$1M_{\odot}$ model to an age of 4.56 Gyr with $f=0.86$ and $\xi=0.6$ (Eqn. 41)
to qualitatively understand how our penetration parameterization modifies a
stellar model. In the top panel of Fig. 10 we display $\nabla\equiv d\ln
T/d\ln P$ from the model which includes convective penetration. Note that
$\nabla$ (green) remains close to $\nabla_{\rm{ad}}$ (purple) below the
Schwarzschild convective boundary ($\nabla_{\rm{ad}}=\nabla_{\rm{rad}}$) in a
penetration zone. After some depth $\nabla\rightarrow\nabla_{\rm{rad}}$
(orange) in the star’s interior. We additionally evolved a standard 1
$M_{\odot}$ MESA model to a 4.56 Gyr age without the inclusion of a PZ. We
compare the sound speed $c$ profiles of the PZ and standard (std) model in the
bottom panel of Fig. 10. When a PZ is present beneath a CZ, $\nabla$
experiences a sharp jump from $\nabla_{\rm{ad}}$ to $\nabla_{\rm{rad}}$ (Fig.
10, top panel), resulting in an acoustic “glitch” in the sound speed profile.

In the model shown in Fig. 10, we find $H_{p}\approx 0.082R_{\odot}$ at the
Schwarzschild CZ boundary, and the depth of the penetration zone in Fig. 10 is
$0.042R_{\odot}\sim 0.5H_{p}$. The inclusion of this PZ leads to an
$\mathcal{O}$(2%) increase in $c$ near the base of the solar convection zone.
Helioseismic observations suggest a similar increase below the base of the
solar convection zone (e.g., Christensen-Dalsgaard et al., 2011, their Fig.
17). The difference $\Delta c=c_{\rm{PZ}}-c_{\rm{std}}$ that we see in this
stellar model of the Sun (Fig. 10) has the same sign and roughly the same
shape. However, the magnitude of the change in $c$ is larger than is observed;
literature values include $\Delta c/c\approx\mathcal{O}(1\%)$ (Bergemann &
Serenelli, 2014) and $\Delta c^{2}/c^{2}\approx\mathcal{O}(0.4\%)$
(Christensen-Dalsgaard et al., 2011), and our sound speed bump is located at a
different radius than the observed bump. Other helioseismic studies have
argued that that the solar PZ depth cannot be larger than $\mathcal{O}$(0.05
$H_{p}$), because larger PZs would result in larger glitches than are detected
(see Sct. 7.2.1 of Basu, 2016, for a nice review). It is interesting, however,
that the width of the PZ in Fig. 10 is strikingly similar to the inferred
width of the tachocline $(0.039\pm 0.013)R_{\odot}$ that is reported by
Charbonneau et al. (1999).

It is unsurprising that our Boussinesq-based model only qualitatively matches
observational constraints for the solar CZ. The solar convection zone is
highly stratified ($\sim$14 density scale heights), and we neglected density
stratification in this work. Furthermore, the solar model used here is
essentially a “stock” MESA model and has obvious disagreements with the solar
model S (see Fig. 1 in Christensen-Dalsgaard et al., 2011, where the
Schwarzschild base of the CZ is $r/R_{\odot}\approx 0.712$, whereas the one in
Fig. 10 is at $r/R_{\odot}\approx 0.75$). Despite the limitations of this
minimal proof of concept, Fig. 10 shows that our parameterization can produce
penetration zones in 1D models with measurable acoustic glitches. In a future
paper, we will produce more realistic models by building upon our
parameterization to include the crucial effects of density stratification. We
note briefly that the theory in e.g., Eqn. 44 only knows about integral
quantities of the convection and does not therefore know about quantities like
the filling factor of upflows and downflows which stratification would modify.
We suspect that dynamical differences that arise from including stratification
would manifest as changes in $f$ and $\xi$, but a detailed exploration is
beyond the scope of this work.

##  7 Discussion

In this work, we presented dynamical simulations of convective penetration, in
which convection mixes ${\nabla\rightarrow\nabla_{\rm{ad}}}$ beyond the
Schwarzschild boundary. To understand these simulations, we used an integral
constraint (reminiscent of Roxburgh, 1989) and flux-based arguments (similar
to Zahn, 1991) to derive a parameterization of convective penetration
according to the convective flux and viscous dissipation. In doing so, we have
laid down the first steps (Eqns. 14 & 44) towards incorporating convective
penetration into stellar structure codes. We parameterized the viscous
dissipation into a bulk-CZ portion ($f$) and a portion in the extended
penetrative region ($\xi$), and derived predictions for how the height of a
penetrative region $\delta_{\rm{p}}$ should scale with these measurable
parameters and a new flux-based “penetration parameter” $\mathcal{P}$. We
designed and analyzed two sets of simulations which showed good agreement with
these theoretical predictions. These simulations differ from past studies
because we separately specify $\mathcal{P}$ and the stiffness $\mathcal{S}$,
and we allow the simulations to evolve for a very long time or use numerical
techniques for rapid evolution. We briefly examined what the impliciations of
this theory could be for a simple stellar model.

Our simulation results suggest that stellar convection zones could be bounded
by sizeable penetration zones. In extreme simulations, we observe penetration
zones which are as large as the convection zones they accompany; however, for
realistic stellar values ($\mathcal{P}\approx 1$), we find that they may be as
large as 20-30% of the convective zone length scale ($\sim$the mixing length).

The simulations we presented in this work use a simplified setup to test the
basic tenets of our theory. In particular, they demonstrate that the shape of
the flux near the convective boundary and the viscous dissipation together
determine the height of the penetration zone. The precise values of the
parameters $f$ and $\xi$ achieved in natural, turbulent, fully compressible,
spherical stellar convection may be different from those presented in e.g.,
Fig. 7 and Eqn. 41 here. Future work should aim to understand how these
parameters and the theory presented in e.g., Eqn. 44 change when more
realistic effects are taken into account.

Stellar opacities and thus stellar radiative conductivities are functions of
thermodynamic variables rather than radial location. The formation of a
penetration zone will therefore affect the conductivity profile and
$\nabla_{\rm{rad}}$, which will in turn affect the location of the
Schwarzschild boundary and the estimate of how deep the penetration zone
should be. In other words, convective penetration and entrainment _both_ occur
in realistic settings, and their combined effects should be studied. Future
work should follow e.g., Käpylä et al. (2017) and implement realistic opacity
profiles which evolve self-consistently with the thermodynamic state in order
to understand how these effects feedback into one another.

Our work here assumes a uniform composition through the convective and
radiative region. Convective boundaries often coincide with discontinuities in
composition profiles (Salaris & Cassisi, 2017). Future work should determine
if stabilizing composition gradients can prevent the formation of the
penetration zones seen here.

Furthermore, stellar fluid dynamics exist in the regime of Pr$\,\ll 1$
(Garaud, 2021). Dynamics in this regime may be different from those in the
regime of Pr $\lesssim 1$ that we studied here, which in theory could affect
$f$ and $\xi$. Recently, Käpylä (2021) found that convective flows exhibited
more penetration at low Pr than high Pr. Future work should aim to understand
whether $f$ and/or $\xi$ depend strongly on $\mathrm{Pr}$ in the turbulent
regime.

Two other interesting complications in stellar contexts are rotation and
magnetism. In the rapidly rotating limit, rotation creates quasi-two-
dimensional flows, which could affect the length scales on which dissipation
acts and thus modify $f$. Furthermore, magnetism adds an additional ohmic
dissipation term, which could in theory drastically change our hydrodynamical
measurement of $f$.

In summary, we have unified Roxburgh (1989)’s integral constraint with Zahn
(1991)’s theory of flux-dependent penetration into a parameterized theory of
convective penetration. We tested this theory with simulations and found good
agreement between the theory and our simulations. In future work, we will use
simulations to test some of the complicating factors we discussed here and aim
to more robustly implement convective penetration into MESA.

We thank Keaton Burns, Matt Browning, Matteo Cantiello, Geoff Vasil, and Kyle
Augustson for useful discussions and/or questions which improved the content
of this manuscript. Ben Brown thanks Jeffrey Oishi for many years of
discussions about overshooting convection. EHA is funded as a CIERA
Postdoctoral fellow and would like to thank CIERA and Northwestern University.
We acknowledge the hospitality of Nordita during the program “The Shifting
Paradigm of Stellar Convection: From Mixing Length Concepts to Realistic
Turbulence Modelling,” where the groundwork for this paper was set. This work
was supported by NASA HTMS grant 80NSSC20K1280 and NASA SSW grant
80NSSC19K0026. Computations were conducted with support from the NASA High End
Computing (HEC) Program through the NASA Advanced Supercomputing (NAS)
Division at Ames Research Center on Pleiades with allocation GID s2276. The
Flatiron Institute is supported by the Simons Foundation.

##  Appendix A Accelerated Evolution

As demonstrated in Fig. 6, the time evolution of simulations which start from
a state based on the Schwarzschild criterion can be prohibitively long. In
Anders et al. (2018), we explored the long time evolution of simple convective
simulations and found that fast-forwarding the evolution of a convective
simulation’s internal energy and thermal structure can be done accurately.
This can be done because the convective dynamics converge rapidly even if the
thermal profile converges slowly. This same separation of scales is observed
in the penetrative dynamics in this work, and so similar techniques should be
applicable.

To more quickly determine the final size of the evolved penetration zones we
use the following algorithm.

  1. 1.

Once a simulation has a volume-averaged Reynolds number greater than 1, we
wait 10 freefall times to allow dynamical transients to pass.

  2. 2.

We measure the departure points ($\delta_{0.1}$, $\delta_{0.5}$,
$\delta_{0.9}$) every freefall time, and store this information for 30
freefall times.

  3. 3.

We linearly fit each of the departure points’ evolution against time using
NumPy’s polyfit function. We assume that convective motions influence
$\delta_{0.1}$ and $\delta_{0.5}$ more strongly than $\delta_{0.9}$. We
measure the time-evolution of the convective front
$\frac{d\delta_{\rm{p}}}{dt}$ by averaging the slope of the linear fits for
$\delta_{0.1}$ and $\delta_{0.5}$.

  4. 4.

We take a large “time step” of size $\tau_{\rm{AE}}$ forward. We calculate
$\Delta\delta_{p}=\tau_{\rm{AE}}\frac{d\delta_{\rm{p}}}{dt}$.

     * •

If $\Delta\delta_{p}<0.005$, we erase the first 15 time units worth of
departure point measures and return to step 2 for 15 time units.

     * •

If $\Delta\delta_{p}$ is large, we adjust the top of the PZ by setting
$\delta_{0.5,\rm{new}}=\langle\delta_{0.5}\rangle_{t}+\Delta\delta_{p}$
(angles represent a time average). If $|\Delta\delta_{p}|>0.05$, we limit its
value to 0.05. We calculate the width of the PZ-RZ boundary layer $d_{w}$ as
the minimum of $\langle\delta_{0.9}-\delta_{0.5}\rangle_{t}$ and
$\langle\delta_{0.5}-\delta_{0.1}\rangle_{t}$. We adjust the mean temperature
gradient to

$\nabla=\nabla_{\rm{ad}}+H(z;\delta_{\rm{0.5,\rm{new}}},d_{w})\Delta\nabla,$
(A1)

where $H$ is defined in Eqn. 34 and
$\Delta\nabla=\nabla_{\rm{rad}}-\nabla_{\rm{ad}}$. We also multiply the
temperature perturbations and full convective velocity field by
$(1-H(z;1,0.05))$. This sets all fluctuations above the nominal Schwarzschild
convection zone to zero, thereby avoiding any strange dynamical transients
caused by the old dynamics at the radiative-convective boundary (which has
moved as a result of this process).

  5. 5.

Return to step 1.

In general, the initial profile of $\overline{T}$ that we use when we start
our simulations is given by Eqn. A1 with a value $\delta_{0.5,\rm{new}}\geq
0$. We then evolve $\overline{T}$ towards a statistically stationary state
using the above algorithm and standard timestepping. If a simulation returns
to step 2 from step 4 ten times over the course of its evolution, we assume
that it has converged near its answer, stop this iterative loop, and allow the
simulation to timestep normally. Additionally, in some simulations, we ensure
that this process occurs no more than 25 times. This process effectively
removes the long diffusive thermal evolution on display in the upper left
panel of Fig. 6 by immediately setting the mean temperature profile to the
radiative profile above the PZ.

In Fig. 11, we plot in black the time evolution of $\delta_{\rm{p}}$ and $f$
in Case I simulations with $\mathcal{S}=10^{3}$, $\mathcal{R}=400$, and
$\mathcal{P}_{D}=[1,2,4]$. We overplot the evolution of simulations which use
this accelerated evolution (AE) procedure using orange and green lines. Time
units on the x-axis are normalized in terms of the total simulation run time
in order to more thoroughly demonstrate the evolutionary differences between
standard timestepping and AE. However, the AE simulations are much shorter:
the vertical green-and-yellow lines demonstrate how long the AE simulation ran
compared to the standard timestepping simulation (so for $\mathcal{P}_{D}=1$,
the AE simulations only took $\sim 1/4$ as long; for $\mathcal{P}_{D}=2$, they
took $\sim 1/10$ as long; for $\mathcal{P}_{D}=4$, they took $\sim 1/20$ as
long). AE simulations with orange lines start with PZ heights which are much
larger than the final height, while green line solution start with initial PZ
heights which are smaller than the expected height. Regardless of our choice
of initial condition, we find that this AE procedure quickly evolves our
simulations to within a few percent of the final value. After converging to
within a few percent of the proper penetration zone height, this AE procedure
continues to iteratively “jitter” around the right answer until the
convergence criterion we described above are met. These jitters can be seen in
the top panels of Fig. 11, where the solution jumps away from the proper
answer in one AE iteration before jumping back towards it in the next
iteration. If the PZ height continues to noticeably vary on timescales of a
few hundred freefall times, we continue to timestep the simulations until the
changes of $\delta_{\rm{p}}$ have diminished.

Figure 11:  (top row) Time traces of $\delta_{0.5}$ for simulations using
standard timestepping (black lines), accelerated simulations with large
initial values of $\delta_{\rm{p}}$ (orange lines), and accelerated
simulations with small initial values of $\delta_{\rm{p}}$ (green lines). Thin
horizontal lines denote the equilibrated value of $\delta_{0.5}$. Accelerated
evolution timesteps can be seen as jumps in the $\delta_{\rm{p}}$ trace. After
converging to within a few percent, the accelerated evolution procedure
“jitters” around the equilibrated value. Time units are normalized by the
total run time of the simulation. Accelerated simulations were run for
$t_{\rm{sim}}=3000$ freefall times. The standard timestepping (black line)
simulations were run for $t_{\rm{sim}}=1.2\times 10^{4}$
($\mathcal{P}_{D}=1$), $t_{\rm{sim}}=3.2\times 10^{4}$ ($\mathcal{P}_{D}=2$),
and $t_{\rm{sim}}=6.7\times 10^{4}$ ($\mathcal{P}_{D}=4$) freefall times. The
vertical green-and-yellow lines show the total simulation time of the
accelerated simulation in terms of the direct simulation time; i.e., the
accelerated simulation converged in only $\sim$ 5% of the simulation time of
the direct simulation for $\mathcal{P}=4$. (Bottom row) Rolling average of $f$
over 200 freefall times, plotted in the same way as $\delta_{0.5}$.

##  Appendix B MESA implementation

Our 1D stellar evolution calculations were performed using the Modules for
Experiments in Stellar Astrophysics software instrument (Paxton et al., 2011,
2013, 2015, 2018, 2019, MESA).

###  B.1 Input Physics

The MESA EOS is a blend of the OPAL (Rogers & Nayfonov, 2002), SCVH (Saumon et
al., 1995), FreeEOS (Irwin, 2004), HELM (Timmes & Swesty, 2000), and PC
(Potekhin & Chabrier, 2010) EOSes.

Radiative opacities are primarily from OPAL (Iglesias & Rogers, 1993, 1996),
with low-temperature data from Ferguson et al. (2005) and the high-
temperature, Compton-scattering dominated regime by Buchler & Yueh (1976).
Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al., 2010) plus
additional tabulated weak reaction rates Fuller et al. (1985); Oda et al.
(1994); Langanke & Martínez-Pinedo (2000). (For MESA versions before 11701):
Screening is included via the prescriptions of Salpeter (1954); Dewitt et al.
(1973); Alastuey & Jancovici (1978); Itoh et al. (1979).

###  B.2 Penetration Implementation

Here we describe a first implementation of Eqn. 44 in MESA. We note that this
impelementation is likely not universal or robust enough to be used in most
complex stellar models, but it is robust enough to time-step stably and
produce the results displayed in Sct. 6. Future work should improve upon this
model.

To find the extent of the penetrative region we write Eqn. (44) as

$\displaystyle(1-f)\int_{\rm CZ}L_{\rm conv}dr=\int_{\rm PZ}\left(\xi fL_{\rm
conv,avg,CZ}+L_{\rm conv}\right)dr,$ (B1)

where $L_{\rm conv,avg,CZ}$ is the average of $L_{\rm conv}$ in the convection
zone and $L_{\rm conv}$ in the penetrative region is given by

$\displaystyle L_{\rm conv}=\frac{L_{\rm
rad}}{\nabla_{r}}(\nabla_{a}-\nabla_{r}),$ (B2)

which is the excess luminosity carried if the temperature gradient in the
radiative zone is adiabatic.

We first integrate the left-hand side of Eqn. (B1) over the convection zone
and further use that to evaluate $L_{\rm conv,avg,CZ}$. Next we integrate the
right-hand side of the same away from the convective boundary into the
radiative zone until the equation is satisfied. The point where this
integration stops is the edge of the penetrative region.

We then implement convective penetration in stellar evolution with two
modifications. First, we add an extra chemical mixing term in the penetration
zone with a scale of $D\approx H_{p}(L/4\pi r^{2}\rho)^{1/3}$, which is
roughly the scale of the convective diffusivity. The precise choice of
diffusivity here does not matter, as any plausible scale will be enough to
eliminate any composition gradient on evolutionary time-scales. Secondly, we
override the default routine in MESA for determining $\nabla$ and instead have
the solver reduce $\nabla_{a}-\nabla$ by 90 per cent in the penetrative zone.

Using this procedure with $f=0.86$ and $\xi=0.6$, and timestepping a solar
model to the age of the current Sun ($\sim$ 4.5 Gyr), we find the profile
displayed in Sec. 6.

###  B.3 Models

Models were constructed to reasonably reproduce the present-day Sun and based
on the 2019 MESA summer school lab by Pinsonneault (2019). Inlists and the
run_star_extras source code are available in a Zenodo repository (Anders et
al., 2021).

##  Appendix C Table of simulation parameters

Input parameters and summary statistics of the simulations presented in this
work are shown in Table 1.

Table 1: Table of simulation information.

Type | $\mathcal{P}$ | $\mathcal{S}$ | $\mathcal{R}$ | $nx\times ny\times nz$ | $t_{\rm{sim}}$ | $(\delta_{0.1},\delta_{0.5},\delta_{0.9})$ | $f$ | $\xi$ | $\langle u\rangle$
---|---|---|---|---|---|---|---|---|---
“Standard timestepping” simulations |  |  |  |  |  |  |
D | $1.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $12347$ | (0.078, 0.112, 0.136) | 0.810 | 0.682 | 0.618
D | $2.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $32057$ | (0.200, 0.230, 0.254) | 0.749 | 0.601 | 0.639
D | $4.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $66557$ | (0.445, 0.472, 0.496) | 0.668 | 0.562 | 0.619
“Accelerated Evolution” simulations |  |  |  |  |  |  |
D | $4.0$ | $10^{2}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $5000$ | (0.377, 0.505, 0.581) | 0.654 | 0.526 | 0.617
D | $4.0$ | $3.0\cdot 10^{2}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $5000$ | (0.420, 0.477, 0.514) | 0.663 | 0.551 | 0.618
D | $10^{-1}$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $4561$ | (0.017, 0.042, 0.069) | 0.831 | 0.769 | 0.588
D | $3.0\cdot 10^{-1}$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $4681$ | (0.030, 0.064, 0.092) | 0.814 | 0.804 | 0.620
D | $1.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $3000$ | (0.082, 0.116, 0.140) | 0.804 | 0.690 | 0.624
D | $2.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $5000$ | (0.199, 0.228, 0.252) | 0.750 | 0.597 | 0.638
D | $4.0$ | $10^{3}$ | $2.5\cdot 10^{1}$ | 16x16x256 | $3000$ | (0.321, 0.379, 0.437) | 0.772 | 0.274 | 0.343
D | $4.0$ | $10^{3}$ | $5.0\cdot 10^{1}$ | 32x32x256 | $3000$ | (0.398, 0.442, 0.487) | 0.732 | 0.358 | 0.423
D | $4.0$ | $10^{3}$ | $10^{2}$ | 32x32x256 | $3000$ | (0.469, 0.503, 0.534) | 0.672 | 0.464 | 0.484
D | $4.0$ | $10^{3}$ | $2.0\cdot 10^{2}$ | 64x64x256 | $3000$ | (0.485, 0.515, 0.542) | 0.648 | 0.546 | 0.548
D | $4.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $5000$ | (0.452, 0.480, 0.505) | 0.667 | 0.553 | 0.617
D | $4.0$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $3000$ | (0.407, 0.434, 0.455) | 0.689 | 0.566 | 0.678
D | $4.0$ | $10^{3}$ | $1.6\cdot 10^{3}$ | 128x128x256 | $3000$ | (0.366, 0.397, 0.419) | 0.709 | 0.574 | 0.720
D | $4.0$ | $10^{3}$ | $3.2\cdot 10^{3}$ | 256x256x256 | $3235$ | (0.321, 0.358, 0.381) | 0.723 | 0.605 | 0.746
D | $4.0$ | $10^{3}$ | $6.4\cdot 10^{3}$ | 384x384x384 | $414$ | (0.277, 0.315, 0.335) | 0.744 | 0.605 | 0.757
D | $6.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $6000$ | (0.620, 0.647, 0.667) | 0.635 | 0.532 | 0.597
D | $8.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 128x128x512 | $4357$ | (0.732, 0.759, 0.779) | 0.640 | 0.481 | 0.592
D | $10^{1}$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 128x128x512 | $4226$ | (0.858, 0.885, 0.904) | 0.630 | 0.453 | 0.587
D | $4.0$ | $3.0\cdot 10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x512 | $1170$ | (0.437, 0.454, 0.469) | 0.672 | 0.581 | 0.619
D/SF | $4.0$ | $10^{3}$ | $5.0\cdot 10^{1}$ | 32x32x256 | $5000$ | (0.435, 0.477, 0.516) | 0.680 | 0.418 | 0.505
D/SF | $4.0$ | $10^{3}$ | $10^{2}$ | 32x32x256 | $5000$ | (0.482, 0.516, 0.547) | 0.638 | 0.543 | 0.573
D/SF | $4.0$ | $10^{3}$ | $2.0\cdot 10^{2}$ | 64x64x256 | $5000$ | (0.490, 0.520, 0.547) | 0.634 | 0.589 | 0.640
D/SF | $4.0$ | $10^{3}$ | $4.0\cdot 10^{2}$ | 64x64x256 | $8000$ | (0.474, 0.502, 0.531) | 0.651 | 0.588 | 0.693
D/SF | $4.0$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $5000$ | (0.410, 0.437, 0.461) | 0.683 | 0.587 | 0.732
D/SF | $4.0$ | $10^{3}$ | $1.6\cdot 10^{3}$ | 128x128x256 | $5710$ | (0.368, 0.400, 0.426) | 0.703 | 0.590 | 0.758
D/SF | $4.0$ | $10^{3}$ | $3.2\cdot 10^{3}$ | 256x256x256 | $3917$ | (0.320, 0.357, 0.388) | 0.725 | 0.595 | 0.772
L | $10^{-2}$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $1139$ | (0.017, 0.030, 0.051) | 0.873 | 0.783 | 0.445
L | $3.0\cdot 10^{-2}$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $929$ | (0.020, 0.044, 0.070) | 0.863 | 0.782 | 0.448
L | $10^{-1}$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $1142$ | (0.081, 0.076, 0.102) | 0.848 | 0.725 | 0.450
L | $3.0\cdot 10^{-1}$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $1109$ | (0.076, 0.129, 0.157) | 0.825 | 0.655 | 0.451
L | $1.0$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $3000$ | (0.182, 0.225, 0.251) | 0.787 | 0.599 | 0.442
L | $2.0$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $3000$ | (0.278, 0.315, 0.340) | 0.759 | 0.570 | 0.436
L | $4.0$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $10000$ | (0.399, 0.431, 0.455) | 0.737 | 0.518 | 0.428
L | $8.0$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $5000$ | (0.519, 0.545, 0.562) | 0.718 | 0.484 | 0.421
L | $1.6\cdot 10^{1}$ | $10^{3}$ | $8.0\cdot 10^{2}$ | 128x128x256 | $8000$ | (0.687, 0.709, 0.723) | 0.700 | 0.442 | 0.417

Note. — Simulation type is specified as “D” for discontinuous/Case I or “L”
for linear/Case II. “D/SF” simulations have stress-free boundary conditions.
Input control parameters are listed for each simulation: the penetration
parameter $\mathcal{P}$, stiffness $\mathcal{S}$, and freefall Reynolds number
$\mathcal{R}$. We also note the coefficient resolution (Chebyshev coefficients
$nz$ and Fourier coeficients $nx$, $ny$). We report the number of freefall
time units each simulation was run for $t_{\rm{sim}}$. Time-averaged values of
the departure heights ($\delta_{0.1}$, $\delta_{0.5}$, $\delta_{0.9}$), the
dissipation fraction $f$, and the dissipation fall-off $\xi$, as well as the
average convection zone velocity $\langle u\rangle$ are reported. We take
these time averages over the final 1000 freefall times or half of the
simulation, whichever is shorter.

## References

  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Springer)
  * Ahlers et al. (2009) Ahlers, G., Grossmann, S., & Lohse, D. 2009, Reviews of Modern Physics, 81, 503, doi: 10.1103/RevModPhys.81.503
  * Alastuey & Jancovici (1978) Alastuey, A., & Jancovici, B. 1978, ApJ, 226, 1034, doi: 10.1086/156681
  * Anders et al. (2018) Anders, E. H., Brown, B. P., & Oishi, J. S. 2018, Physical Review Fluids, 3, 083502, doi: 10.1103/PhysRevFluids.3.083502
  * Anders et al. (2021) Anders, E. H., Jermyn, A. S., Lecoanet, D., & Brown, B. P. 2021, Supplemental Materials for ”Convective penetration: It exists and we found it”, 1.0.0, Zenodo, doi: 10.5281/zenodo.5131118
  * Asplund et al. (2021) Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, arXiv e-prints, arXiv:2105.01661.  https://arxiv.org/abs/2105.01661
  * Bahcall et al. (2005) Bahcall, J. N., Serenelli, A. M., & Basu, S. 2005, ApJ, 621, L85, doi: 10.1086/428929
  * Basu (2016) Basu, S. 2016, Living Reviews in Solar Physics, 13, 2, doi: 10.1007/s41116-016-0003-4
  * Basu & Antia (2004) Basu, S., & Antia, H. M. 2004, ApJ, 606, L85, doi: 10.1086/421110
  * Bergemann & Serenelli (2014) Bergemann, M., & Serenelli, A. 2014, Solar Abundance Problem, 245–258, doi: 10.1007/978-3-319-06956-2_21
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Browning et al. (2004) Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512, doi: 10.1086/380198
  * Brummell et al. (2002) Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825, doi: 10.1086/339626
  * Brun & Browning (2017) Brun, A. S., & Browning, M. K. 2017, Living Reviews in Solar Physics, 14, 4, doi: 10.1007/s41116-017-0007-8
  * Brun et al. (2017) Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192, doi: 10.3847/1538-4357/aa5c40
  * Buchler & Yueh (1976) Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440, doi: 10.1086/154847
  * Burns et al. (2020) Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., & Brown, B. P. 2020, Physical Review Research, 2, 023068, doi: 10.1103/PhysRevResearch.2.023068
  * Carlos et al. (2019) Carlos, M., Meléndez, J., Spina, L., et al. 2019, MNRAS, 485, 4052, doi: 10.1093/mnras/stz681
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Caswell et al. (2021) Caswell, T. A., Droettboom, M., Lee, A., et al. 2021, matplotlib/matplotlib: REL: v3.3.4, v3.3.4, Zenodo, doi: 10.5281/zenodo.4475376
  * Charbonneau et al. (1999) Charbonneau, P., Christensen-Dalsgaard, J., Henning, R., et al. 1999, ApJ, 527, 445, doi: 10.1086/308050
  * Christensen-Dalsgaard (2021) Christensen-Dalsgaard, J. 2021, Living Reviews in Solar Physics, 18, 2, doi: 10.1007/s41116-020-00028-3
  * Christensen-Dalsgaard et al. (2011) Christensen-Dalsgaard, J., Monteiro, M. J. P. F. G., Rempel, M., & Thompson, M. J. 2011, MNRAS, 414, 1158, doi: 10.1111/j.1365-2966.2011.18460.x
  * Claret & Torres (2018) Claret, A., & Torres, G. 2018, ApJ, 859, 100, doi: 10.3847/1538-4357/aabd35
  * Couston et al. (2017) Couston, L. A., Lecoanet, D., Favier, B., & Le Bars, M. 2017, Physical Review Fluids, 2, 094804, doi: 10.1103/PhysRevFluids.2.094804
  * Cristini et al. (2017) Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279, doi: 10.1093/mnras/stx1535
  * Currie & Browning (2017) Currie, L. K., & Browning, M. K. 2017, ApJ, 845, L17, doi: 10.3847/2041-8213/aa8301
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240, doi: 10.1088/0067-0049/189/1/240
  * Deardorff et al. (1969) Deardorff, J. W., Willis, G. E., & Lilly, D. K. 1969, Journal of Fluid Mechanics, 35, 7, doi: 10.1017/S0022112069000942
  * Dewitt et al. (1973) Dewitt, H. E., Graboske, H. C., & Cooper, M. S. 1973, ApJ, 181, 439, doi: 10.1086/152061
  * Dumont et al. (2021) Dumont, T., Palacios, A., Charbonnel, C., et al. 2021, A&A, 646, A48, doi: 10.1051/0004-6361/202039515
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53, doi: 10.3847/1538-4357/ab518b
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fuentes & Cumming (2020) Fuentes, J. R., & Cumming, A. 2020, Physical Review Fluids, 5, 124501, doi: 10.1103/PhysRevFluids.5.124501
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1, doi: 10.1086/163208
  * Garaud (2021) Garaud, P. 2021, Physical Review Fluids, 6, 030501, doi: 10.1103/PhysRevFluids.6.030501
  * Goluskin (2016) Goluskin, D. 2016, Internally heated convection and Rayleigh-Bénard convection (Springer International Publishing), arXiv:1506.01656, doi: 10.1007/978-3-319-23941-5
  * Hansen et al. (2004) Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar interiors : physical principles, structure, and evolution (Springer)
  * Herwig (2000) Herwig, F. 2000, A&A, 360, 952.  https://arxiv.org/abs/astro-ph/0007139
  * Higgins & Vink (2020) Higgins, E. R., & Vink, J. S. 2020, A&A, 635, A175, doi: 10.1051/0004-6361/201937374
  * Higl et al. (2021) Higl, J., Müller, E., & Weiss, A. 2021, A&A, 646, A133, doi: 10.1051/0004-6361/202039532
  * Horst et al. (2021) Horst, L., Hirschi, R., Edelmann, P. V. F., Andrassy, R., & Roepke, F. K. 2021, arXiv e-prints, arXiv:2107.02199.  https://arxiv.org/abs/2107.02199
  * Hotta (2017) Hotta, H. 2017, ApJ, 843, 52, doi: 10.3847/1538-4357/aa784b
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hurlburt et al. (1986) Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563, doi: 10.1086/164796
  * Hurlburt et al. (1994) Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J.-P. 1994, ApJ, 421, 245, doi: 10.1086/173642
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752, doi: 10.1086/172958
  * Iglesias & Rogers (1996) —. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Irwin (2004) Irwin, A. W. 2004, The FreeEOS Code for Calculating the Equation of State for Stellar Interiors.  http://freeeos.sourceforge.net/
  * Itoh et al. (1979) Itoh, N., Totsuji, H., Ichimaru, S., & Dewitt, H. E. 1979, ApJ, 234, 1079, doi: 10.1086/157590
  * Jermyn et al. (2018) Jermyn, A. S., Tout, C. A., & Chitre, S. M. 2018, MNRAS, 480, 5427, doi: 10.1093/mnras/sty1831
  * Jørgensen & Weiss (2019) Jørgensen, A. C. S., & Weiss, A. 2019, MNRAS, 488, 3463, doi: 10.1093/mnras/stz1980
  * Käpylä (2019) Käpylä, P. J. 2019, A&A, 631, A122, doi: 10.1051/0004-6361/201834921
  * Käpylä (2021) —. 2021, arXiv e-prints, arXiv:2105.08453.  https://arxiv.org/abs/2105.08453
  * Käpylä et al. (2007) Käpylä, P. J., Korpi, M. J., Stix, M., & Tuominen, I. 2007, in Convection in Astrophysics, ed. F. Kupka, I. Roxburgh, & K. L. Chan, Vol. 239, 437–442, doi: 10.1017/S1743921307000865
  * Käpylä et al. (2017) Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23, doi: 10.3847/2041-8213/aa83ab
  * Kitiashvili et al. (2016) Kitiashvili, I. N., Kosovichev, A. G., Mansour, N. N., & Wray, A. A. 2016, ApJ, 821, L17, doi: 10.3847/2041-8205/821/1/L17
  * Korre et al. (2019) Korre, L., Garaud, P., & Brummell, N. H. 2019, MNRAS, 484, 1220, doi: 10.1093/mnras/stz047
  * Langanke & Martínez-Pinedo (2000) Langanke, K., & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481, doi: 10.1016/S0375-9474(00)00131-7
  * Maeder (1975) Maeder, A. 1975, A&A, 40, 303
  * Marcus et al. (1983) Marcus, P. S., Press, W. H., & Teukolsky, S. A. 1983, ApJ, 267, 795, doi: 10.1086/160915
  * Martinet et al. (2021) Martinet, S., Meynet, G., Ekström, S., et al. 2021, A&A, 648, A126, doi: 10.1051/0004-6361/202039426
  * Meakin & Arnett (2007) Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448, doi: 10.1086/520318
  * Moore & Weiss (1973) Moore, D. R., & Weiss, N. O. 1973, Journal of Fluid Mechanics, 61, 553, doi: 10.1017/S0022112073000868
  * Musman (1968) Musman, S. 1968, Journal of Fluid Mechanics, 31, 343, doi: 10.1017/S0022112068000194
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231, doi: 10.1006/adnd.1994.1007
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pedersen et al. (2021) Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, arXiv e-prints, arXiv:2105.04533.  https://arxiv.org/abs/2105.04533
  * Pinsonneault (1997) Pinsonneault, M. 1997, ARA&A, 35, 557, doi: 10.1146/annurev.astro.35.1.557
  * Pinsonneault (2019) Pinsonneault, M. 2019, Evolution of low mass stars, Zenodo, doi: 10.5281/zenodo.3374959
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82, doi: 10.1002/ctpp.201010017
  * Pratt et al. (2017) Pratt, J., Baraffe, I., Goffrey, T., et al. 2017, A&A, 604, A125, doi: 10.1051/0004-6361/201630362
  * Rempel (2004) Rempel, M. 2004, ApJ, 607, 1046, doi: 10.1086/383605
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064, doi: 10.1086/341894
  * Rogers & Glatzmaier (2005) Rogers, T. M., & Glatzmaier, G. A. 2005, ApJ, 620, 432, doi: 10.1086/423415
  * Rogers et al. (2006) Rogers, T. M., Glatzmaier, G. A., & Jones, C. A. 2006, ApJ, 653, 765, doi: 10.1086/508482
  * Roxburgh (1978) Roxburgh, I. W. 1978, A&A, 65, 281
  * Roxburgh (1989) —. 1989, A&A, 211, 361
  * Roxburgh (1992) —. 1992, A&A, 266, 291
  * Roxburgh (1998) —. 1998, Ap&SS, 261, 43, doi: 10.1023/A:1002029105802
  * Saikia et al. (2000) Saikia, E., Singh, H. P., Chan, K. L., Roxburgh, I. W., & Srivastava, M. P. 2000, ApJ, 529, 402, doi: 10.1086/308249
  * Salaris & Cassisi (2017) Salaris, M., & Cassisi, S. 2017, Royal Society Open Science, 4, 170192, doi: 10.1098/rsos.170192
  * Salpeter (1954) Salpeter, E. E. 1954, Australian Journal of Physics, 7, 373, doi: 10.1071/PH540373
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Shaviv & Salpeter (1973) Shaviv, G., & Salpeter, E. E. 1973, ApJ, 184, 191, doi: 10.1086/152318
  * Singh et al. (1995) Singh, H. P., Roxburgh, I. W., & Chan, K. L. 1995, A&A, 295, 703
  * Singh et al. (1998) —. 1998, A&A, 340, 178
  * Spiegel & Veronis (1960) Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442, doi: 10.1086/146849
  * Tian et al. (2009) Tian, C.-L., Deng, L.-C., & Chan, K.-L. 2009, MNRAS, 398, 1011, doi: 10.1111/j.1365-2966.2009.15178.x
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Toppaladoddi & Wettlaufer (2018) Toppaladoddi, S., & Wettlaufer, J. S. 2018, Physical Review Fluids, 3, 043501, doi: 10.1103/PhysRevFluids.3.043501
  * Viallet et al. (2013) Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1, doi: 10.1088/0004-637X/769/1/1
  * Viallet et al. (2015) Viallet, M., Meakin, C., Prat, V., & Arnett, D. 2015, A&A, 580, A61, doi: 10.1051/0004-6361/201526294
  * Viani & Basu (2020) Viani, L. S., & Basu, S. 2020, ApJ, 904, 22, doi: 10.3847/1538-4357/abba17
  * Vinyoles et al. (2017) Vinyoles, N., Serenelli, A. M., Villante, F. L., et al. 2017, ApJ, 835, 202, doi: 10.3847/1538-4357/835/2/202
  * Wang & Ruuth (2008) Wang, D., & Ruuth, S. J. 2008, Journal of Computational Mathematics, 26, 838.  http://www.jstor.org/stable/43693484
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015, doi: 10.1103/RevModPhys.74.1015
  * Zahn (1991) Zahn, J. P. 1991, A&A, 252, 179