\section{Introduction}
\label{sec:introduction}

\subsection{Context}
Convection is a crucial mechanism for transporting heat in stars~\citep{woosley_etal_2002, hansen_etal_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 \citep{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 \citep{salaris_cassisi_2017}.
Furthermore, convective motions excite waves, which can be observed and used to constrain the thermodynamic structure of stars \citep{aerts2010, basu2016}.
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 \citep{pinsonneault1997, carlos_etal_2019, dumont_etal_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 \citep{basu_antia_2004, bahcall_etal_2005, bergemann_serenelli_2014, vinyoles_etal_2017, asplund_etal_2021}; modeling and observational discrepancies can be reduced with additional mixing below the convective boundary \citep{christensen-dalsgaard_etal_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 \citep{claret_torres_2018, jermyn_etal_2018, viani_basu_2020, martinet_etal_2021, pedersen_etal_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 \citep{farmer_etal_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 \citep[e.g., mixing length theory,][]{bohm-vitense1958} and convective overshoot \citep{shaviv_salpeter_1973, maeder1975, herwig2000, paxton_etal_2011, paxton_etal_2013, paxton_etal_2018, paxton_etal_2019}.
While some preliminary work has been done to couple 3D dynamical convective simulations with 1D stellar evolution codes \citep{jorgensen_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 \emph{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 \citep{zahn1991, brummell_etal_2002, korre_etal_2019}.
Motions which extend beyond the convective boundary but do not adjust the thermodynamic profiles belong to a process called ``\emph{convective overshoot}'' in the fluid dynamics literature.
Convection zones can expand through a second process called ``\emph{entrainment},'' through which motions erode composition gradients or modify the radiative gradient \citep[][]{meakin_arnett_2007, viallet_etal_2013, cristini_etal_2017, fuentes_cumming_2020, horst_etal_2021}.
The primary focus of this work is a third process called ``\emph{convective penetration}''.
Convective penetration occurs when motions mix the entropy gradient towards the adiabatic in a region that is stable by the Schwarzschild criterion.

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{vertical_dynamics_panels.pdf}
\caption{
Vertical slice through a simulation with $\ensuremath{\mathcal{R}} = 6.4 \times 10^3$, $\ensuremath{\mathcal{P}}_D = 4$ and $\ensuremath{\mathcal{S}} = 10^3$ (see Sec.~\ref{sec:simulation_details}).
The dashed horizontal line denotes the Schwarzschild convective boundary where $\ensuremath{\nabla_{\rm{ad}}} = \ensuremath{\nabla_{\rm{rad}}}$.
The top of the penetrative zone ($\delta_{0.1}$, see Sec.~\ref{sec:simulation_details}) 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 $\ensuremath{\nabla}$ departs from $\ensuremath{\nabla_{\rm{ad}}}$.
(Right) Temperature fluctuations, normalized by their average magnitude at each height to clearly display all dynamical features.
\label{fig:vertical_dynamics_panels}
}
\end{figure*}

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 \citep[e.g.,][]{marcus_etal_1983, zahn1991, browning_etal_2004, rogers_etal_2006, viallet_etal_2015, korre_etal_2019}.
Experiments exhibiting extensive expansion of convection zones via entrainment have a long history \citep[e.g.,][and this process is often confusingly called ``penetration'']{musman1968, deardorff_etal_1969, moore_weiss_1973}.
Modern numerical experiments often examine the importance of the ``stiffness'' $\ensuremath{\mathcal{S}}$ of a radiative-convective interface.
$\ensuremath{\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 \citep{couston_etal_2017, toppaladoddi_wettlaufer_2018}; others find stiffness-dependent pure overshoot \citep{korre_etal_2019}.
A link between $\ensuremath{\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 \citep{hurlburt_etal_1994, saikia_etal_2000, brummell_etal_2002, rogers_glatzmaier_2005, rogers_etal_2006, kitiashvili_etal_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 \citep{hurlburt_etal_1986, singh_etal_1995, browning_etal_2004, brun_etal_2017, pratt_etal_2017}.
Still other simulations show negligible penetration \citep[e.g.,][]{higl_etal_2021}.
Even detailed studies which sought a relationship between penetration depth and stiffness $\ensuremath{\mathcal{S}}$ have presented contradictory results.
Early work by e.g., \citet{hurlburt_etal_1994} and \citet{singh_etal_1995} hinted at a link between $\ensuremath{\mathcal{S}}$ and penetration length, at least for low values of $\ensuremath{\mathcal{S}}$.
Subsequent simulations by \citet{brummell_etal_2002} exhibit a weak scaling of penetration depth with $\ensuremath{\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 \citet{rogers_glatzmaier_2005} demonstrate a negligible scaling of the penetration depth against $\ensuremath{\mathcal{S}}$ at moderate values of $\ensuremath{\mathcal{S}}$.
Prior simulations thus consistently show hints of penetration at low $\ensuremath{\mathcal{S}}$ \citep[where results may not be relevant for stars,][]{couston_etal_2017}, but present confusing and contradictory results at moderate-to-high $\ensuremath{\mathcal{S}}$.

There are hints in the literature that convective penetration may depend on energy fluxes.
\citet{roxburgh1978, roxburgh1989, roxburgh1992, roxburgh1998} 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.
\citet{zahn1991} theorized that convective penetration should depend only on how steeply the radiative temperature gradient varies at the convective boundary.
Following \citet{zahn1991}'s work, \citet{rempel2004} 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 \citep{singh_etal_1998, kapyla_etal_2007, tian_etal_2009, hotta2017, kapyla2019}.
Furthermore, in the limit of low stiffness, the simulations of \citet{hurlburt_etal_1994} and \citet{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.


$\,$
\subsection{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.
\begin{quote}
\emph{
    In this paper, we present simulations which exhibit convective penetration.
}
\end{quote}
This process is phenomenologically described in Sec.~\ref{sec:central_results}.
\editone{
    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 \citet{roxburgh1989} and \citet{zahn1991}.
\begin{quote}
\emph{
We find that the extent of convective penetration depends strongly on the shape and magnitude of the radiative gradient near the convective boundary.
}
\end{quote}
Thus, the penetration length can be calculated using the radiative conductivity (or opacity) \emph{profile} near the convective boundary.

We present these findings as follows.
In Sec.~\ref{sec:central_results}, we present the central finding of this work: penetration zones in nonlinear convective simulations.
In Sec.~\ref{sec:theory}, we describe the equations used and derive a parameterized theory of convective penetration.
In Sec.~\ref{sec:simulation_details}, we describe our simulation setup and parameters.
In Sec.~\ref{sec:results}, we present the results of these simulations, with a particular focus on the height of the penetrative regions.
In Sec.~\ref{sec:solar_model}, we create and discuss a stellar model in MESA which has convective penetration.
Finally, we discuss pathways for future work in Sec.~\ref{sec:discussion}.

\section{Central Result: Convective Penetration}
\label{sec:central_results}

In Fig.~\ref{fig:vertical_dynamics_panels}, 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.~\ref{fig:vertical_dynamics_panels} in Fig.~\ref{fig:grad_profiles} by displaying time- and horizontally-averaged 1D profiles of the temperature gradient $\ensuremath{\nabla}$ (defined in Sec.~\ref{sec:theory}).
The adiabatic gradient $\ensuremath{\nabla_{\rm{ad}}}$ (purple) has a constant value in the simulation.
Also shown is the radiative gradient $\ensuremath{\nabla_{\rm{rad}}}$ (orange).
The domain exhibits a classical Schwarzshild-unstable convection zone (CZ) for $z \lesssim 1.04$ where $\ensuremath{\nabla_{\rm{rad}}} > \ensuremath{\nabla_{\rm{ad}}}$; the upper boundary of this region is denoted by a dashed vertical line.
Above this point, $\ensuremath{\nabla_{\rm{rad}}} < \ensuremath{\nabla_{\rm{ad}}}$ and the domain would be considered stable by the Schwarzschild criterion.
However, the evolved convective dynamics in Fig.~\ref{fig:vertical_dynamics_panels} have raised $\ensuremath{\nabla} \rightarrow \ensuremath{\nabla_{\rm{ad}}}$ in an extended penetration zone (PZ) which extends from $1.04 \lesssim z \lesssim 1.3$.
\editone{
    Above $z \gtrsim 1.4$, $\ensuremath{\nabla} \approx \ensuremath{\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 $\ensuremath{\nabla}$ adjusts from $\ensuremath{\nabla_{\rm{ad}}}$ to $\ensuremath{\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.

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{grad_profiles.pdf}
\caption{
Horizontally- and temporally-averaged profiles of the thermodynamic gradients from the simulation in Fig.~\ref{fig:vertical_dynamics_panels}.
We plot $\ensuremath{\nabla}$ (green) compared to $\ensuremath{\nabla_{\rm{ad}}}$ (purple, a constant) and $\ensuremath{\nabla_{\rm{rad}}}$ (orange); note the extended penetration zone (PZ) where $\ensuremath{\nabla} \approx \ensuremath{\nabla_{\rm{ad}}} > \ensuremath{\nabla_{\rm{rad}}}$.
    \editone{
        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.
    }
\label{fig:grad_profiles}
}
\end{figure}



\section{Theory}
\label{sec: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.~\ref{sec:theory_equations}, we describe our equations and problem setup and define the heat fluxes.
In Sec.~\ref{sec:theory_energy}, 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 \editone{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.

\subsection{Equations \& flux definitions}
\label{sec:theory_equations}
Throughout this work, we will utilize a modified version of the incompressible Boussinesq equations,
\begin{align}
&\vec{\nabla}\dot\vec{u} = 0 
\label{eqn:incompressible} \\
&\partial_t \vec{u} + \vec{u}\dot\vec{\nabla}\vec{u} = -\frac{1}{\rho_0}\vec{\nabla} p + \frac{\rho_1}{\rho_0}\vec{g} + \nu\vec{\nabla}^2 \vec{u} 
\label{eqn:momentum} \\
&\partial_t T + \vec{u}\dot\vec{\nabla} T + w \ensuremath{\nabla_{\rm{ad}}} + \vec{\nabla}\dot[-k \vec{\nabla} \overline{T}] = \chi\vec{\nabla}^2 T' + Q
\label{eqn:temperature} \\
&\frac{\rho_1}{\rho_0} = -|\alpha| T.
\label{eqn:boussinesq}
\end{align}
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 $\vec{u}$, \editone{the pressure $p$,} the viscous diffusivity $\nu$, the thermal diffusivity $\chi$, the bulk internal heating $Q$, the adiabatic gradient $\ensuremath{\nabla_{\rm{ad}}}$, and a height-dependent thermal conductivity $k$.
We will consider Cartesian coordinates $(x, y, z)$ with a constant vertical gravity $\vec{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.~\ref{eqn:temperature}, $\bar{T}$ is the horizontally averaged temperature and $T'$ are fluctuations away from that; both of these fields evolve in time according to Eqn.~\ref{eqn:temperature}.

Assuming convection reaches a time-stationary state, the heat fluxes are found by horizontally-averaging then vertically integrating Eqn.~\ref{eqn:temperature} to find
\begin{equation}
\overline{\ensuremath{F_{\rm{tot}}}} = \overline{\ensuremath{F_{\rm{rad}}}} + \overline{\ensuremath{F_{\rm{conv}}}} = \int Q dz + \ensuremath{F_{\rm{bot}}},
\label{eqn:flux_definition}
\end{equation}
where $\ensuremath{F_{\rm{bot}}}$ is the flux carried at the bottom of the domain, and $\overline{\ensuremath{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 $\bar{\ensuremath{F_{\rm{rad}}}} = -k \vec{\nabla} \overline{T}$.
We note that $k$ and $-\partial_z \bar{T}$ fully specify $\bar{\ensuremath{F_{\rm{rad}}}}$ and in turn the convective flux, $\bar{\ensuremath{F_{\rm{conv}}}} = \bar{\ensuremath{F_{\rm{tot}}}} - \bar{\ensuremath{F_{\rm{rad}}}}$.
We define the temperature gradient and radiative temperature gradient 
\begin{equation}
\ensuremath{\nabla} \equiv -\partial_z \bar{T} \qquad
\ensuremath{\nabla_{\rm{rad}}} \equiv \frac{\bar{\ensuremath{F_{\rm{tot}}}}}{k}.
\label{eqn:gradrad_definition}
\end{equation}
We have defined the $\ensuremath{\nabla}$'s as positive quantities to align with stellar structure conventions and intuition.
Marginal stability is achieved when $\ensuremath{\nabla} = \ensuremath{\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 $\ensuremath{\nabla_{\rm{rad}}} = \ensuremath{\nabla_{\rm{ad}}}$ and $\bar{\ensuremath{F_{\rm{conv}}}} = 0$.

The addition of a nonzero $\ensuremath{\nabla_{\rm{ad}}}$ to Eqn.~\ref{eqn:temperature} was derived by \citet{spiegel_veronis_1960} and utilized by e.g., \citet{korre_etal_2019}.
In this work, we have decomposed the radiative diffusivity into a background portion ($\vec{\nabla}\dot \bar{\ensuremath{F_{\rm{rad}}}}$) and a fluctuating portion ($\chi \vec{\nabla}^2 T'$); by doing so, we have introduced a height-dependent $\ensuremath{\nabla_{\rm{rad}}}$ to the equation set while preserving the diffusive behavior on fluctuations felt by classical Rayleigh-B\'{e}nard convection.
Here, we will assume a model in which an unstable convection zone ($\ensuremath{\nabla_{\rm{rad}}} > \ensuremath{\nabla_{\rm{ad}}}$) sits below a stable radiative zone ($\ensuremath{\nabla_{\rm{rad}}} < \ensuremath{\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.

\subsection{Kinetic energy \& the dissipation-flux link}
\label{sec:theory_energy}

\begin{figure}[t!]
\centering
\includegraphics[width=\columnwidth]{theory_profiles.pdf}
\caption{
Temporally- and horizontally-averaged profiles from Eqn.~\ref{eqn:kinetic_energy_1D} in the simulation in Fig.~\ref{fig:vertical_dynamics_panels}.
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 $\bar{\mathcal{F}_{\editone{z}}}$, which go to zero at the bottom boundary and the top of the PZ.
    (bottom) Source terms from Eqn.~\ref{eqn:kinetic_energy_1D} normalized by the maximum of $\overline{\mathcal{B}}$ ($\bar{\mathcal{F}_{\editone{z}}}$ in the upper panel is similarly normalized).
The buoyancy source $\bar{\mathcal{B}}$ changes sign at the Schwarzschild boundary, and $\bar{\Phi}$ is positive-definite.
\label{fig:theory_profiles}
}
\end{figure}

Taking a dot product of the velocity and Eqn.~\ref{eqn:momentum} reveals the kinetic energy equation,
\begin{equation}
\frac{\partial \mathcal{K}}{\partial t}
    + \vec{\nabla}\dot\editone{\vec{\mathcal{F}}}
= \mathcal{B} - \Phi,
\label{eqn:kinetic_energy}
\end{equation}
where we define the kinetic energy $\mathcal{K} \equiv |\vec{u}|^2/2$, the fluxes of kinetic energy $\editone{\vec{\mathcal{F}}} \equiv \left[\vec{u}(\mathcal{K} + p/\rho_0) - \nu\vec{u}\vec{\times}\vec{\omega} \right]$, the buoyant energy generation rate $\mathcal{B} \equiv |\alpha| g w T'$, and the viscous dissipation rate $\Phi \equiv \nu |\vec{\omega}|^2$ where $\vec{\omega} = \vec{\nabla}\vec{\times}\vec{u}$ is the vorticity and $|\vec{u}|^2 = \vec{u}\dot\vec{u}$ \& $|\vec{\omega}|^2 = \vec{\omega}\dot\vec{\omega}$.
We next take a horizontal- and time-average of Eqn.~\ref{eqn:kinetic_energy} (we absorb the time-average into the horizontal-average $\bar{\,\cdot\,}$ notation for simplicity).
Assuming that $\bar{\mathcal{K}}$ reaches a statistically stationary state, convective motions satisfy
\begin{equation}
    \frac{d\bar{\mathcal{F}_{\editone{z}}}}{dz} = \bar{\mathcal{B}} - \bar{\Phi},
\label{eqn:kinetic_energy_1D}
\end{equation}
\editone{where $\mathcal{F}_z$ is the z-component of $\vec{\mathcal{F}}$.}
Each profile in Eqn.~\ref{eqn:kinetic_energy_1D} is shown in Fig.~\ref{fig:theory_profiles} for the simulation whose dynamics are displayed in Fig.~\ref{fig:vertical_dynamics_panels}.
As in Fig.~\ref{fig:grad_profiles}, 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 $\bar{\mathcal{F}_{\editone{z}}}$, neglecting the viscous flux term which is only nonzero in a small region above the bottom boundary.
We see that $\bar{\mathcal{F}_{\editone{z}}}$ is zero at the bottom boundary (left edge of plot) and at the top of the PZ.
In the bottom panel, we plot $\bar{\mathcal{B}}$ and $\bar{\Phi}$; we see that $\bar{\mathcal{B}}$ changes sign at the Schwarzschild CZ boundary, and that $\bar{\Phi}$ is positive-definite.

At the boundaries of the convecting region, $\bar{\mathcal{F}_{\editone{z}}}$ is zero (Fig.~\ref{fig:theory_profiles}, upper panel).
We integrate Eqn.~\ref{eqn:kinetic_energy_1D} vertically between these zeros to find
\begin{equation}
\int \bar{\mathcal{B}}\,dz = \int \bar{\Phi}\,dz.
\label{eqn:integral_constraint}
\end{equation}
Integral constraints of this form are the basis for a broad range of analyses in Boussinesq convection \citep[see e.g.,][]{ahlers_etal_2009, goluskin2016} and were considered in the context of penetrative stellar convection by \citet{roxburgh1989}.
Eqn.~\ref{eqn:integral_constraint} 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 $\ensuremath{\nabla} \rightarrow \ensuremath{\nabla_{\rm{ad}}}$ in both the CZ and PZ.
The buoyant energy generation is proportional to the convective flux, $\bar{\mathcal{B}} = |\alpha|g\bar{w T'} = |\alpha| g \bar{\ensuremath{F_{\rm{conv}}}}$, and is positive in the CZ and negative in the PZ (see Fig.~\ref{fig:theory_profiles}, bottom panel).
Breaking up Eqn.~\ref{eqn:integral_constraint}, we see that
\begin{equation}
\int_{\rm{CZ}} \bar{\mathcal{B}}\,dz \,\,=\,\,
\int_{\rm{CZ}} \bar{\Phi}\,dz + \int_{\rm{PZ}} \bar{\Phi}\,dz + \int_{\rm{PZ}}(-\bar{\mathcal{B}})\,dz.
\label{eqn:constraint_cz_pz_split}
\end{equation}
Eqn.~\ref{eqn:constraint_cz_pz_split} 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.~\ref{eqn:constraint_cz_pz_split} is satisfied.
We see that the viscous dissipation and \editone{buoyant deceleration} felt by flows in the PZ determine its size.
We now define
\begin{equation}
f \equiv \frac{\int_{\rm{CZ}} \bar{\Phi}\,dz}{\int_{\rm{CZ}}\bar{\mathcal{B}}\,dz},
\label{eqn:f_defn}
\end{equation}
the measurable fraction of the buoyant engine consumed by CZ dissipation.
Eqn.~\ref{eqn:constraint_cz_pz_split} can then be rewritten as
\begin{equation}
\frac{\int_{\rm{PZ}}(-\bar{\mathcal{B}})\,dz}{\int_{\rm{CZ}} \bar{\mathcal{B}}\,dz}
+ \frac{\int_{\rm{PZ}} \bar{\Phi}\,dz }{\int_{\rm{CZ}} \bar{\mathcal{B}}\,dz}
= (1 - f).
\label{eqn:first_pz_parameterization}
\end{equation}
We will measure and report the values of $f$ achieved in our simulations in this work.
Eqn.~\ref{eqn:first_pz_parameterization} provides two limits on a hypothetical PZ:
\begin{enumerate}
\item In the limit that $f \rightarrow 0$, viscous dissipation is inefficient.
Reasonably if we also assume that $\int_{\rm{PZ}}\bar{\Phi}\,dz \rightarrow 0$, Eqn.~\ref{eqn:first_pz_parameterization} 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 \citet{roxburgh1989} derived.
\item In the limit that $f \rightarrow 1$, viscous dissipation efficiently counteracts the buoyancy work in the CZ.
Per Eqn.~\ref{eqn:first_pz_parameterization}, 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.
\end{enumerate}
In general, we anticipate from the results of e.g., \citet{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.~\ref{fig:theory_profiles}, bottom panel\footnote{the bulk dynamics suggest by eye $f \sim 0.5$, but due to e.g., the height dependence of $\bar{\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 $\ensuremath{\delta_{\rm{p}}}$ develops above a CZ of depth $L_{\rm{CZ}}$, we model the PZ dissipation as
\begin{equation}
\int_{\rm{PZ}} \bar{\Phi}\,dz = \xi\frac{\ensuremath{\delta_{\rm{p}}}}{L_{\rm{CZ}}}\int_{\rm{CZ}}\bar{\Phi}\,dz = \xi \ensuremath{\delta_{\rm{p}}} \Phi_{\rm{CZ}}.
\label{eqn:xi_defn}
\end{equation}
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 $\bar{\Phi}(z = L_s) \approx \Phi_{\rm{CZ}}$ at the CZ-PZ boundary and that $\bar{\Phi}$ decreases with height in the PZ.
The shape of $\bar{\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 $\bar{\mathcal{B}} \propto \bar{\ensuremath{F_{\rm{conv}}}}$, we rewrite Eqn.~\ref{eqn:first_pz_parameterization},
\begin{equation}
-\frac{\int_{\rm{PZ}}\bar{\ensuremath{F_{\rm{conv}}}}\,dz}{\int_{\rm{CZ}}\bar{\ensuremath{F_{\rm{conv}}}}\,dz} + f\xi\frac{\ensuremath{\delta_{\rm{p}}}}{L_{\rm{CZ}}}
= (1 - f).
\label{eqn:theory_fraction}
\end{equation}
The fundamental result of this theory is Eqn.~\ref{eqn:theory_fraction}, which is a parameterized and generalized form of \citet{roxburgh1989}'s integral constraint.
This equation is also reminiscent of \citet{zahn1991}'s theory, and says that the size of a PZ is set by the profile of $\ensuremath{\nabla_{\rm{rad}}}$ near the convective boundary.
A parameterization like Eqn.~\ref{eqn:theory_fraction} can be implemented in stellar structure codes and used to find the extent of penetration zones under the specification of $f$ and $\xi$.
\editone{
    We note that an implementation of Eqn.~\ref{eqn:theory_fraction} likely requires an \emph{iterative} solve, as the penetration zone depth ($\ensuremath{\delta_{\rm{p}}}$) and thus the PZ integral of the flux, are not known \emph{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{\ensuremath{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 $\ensuremath{\nabla} = \ensuremath{\nabla_{\rm{ad}}}$ in the CZ and PZ),
\begin{equation}
\ensuremath{\mathcal{P}} \equiv -\frac{\bar{\ensuremath{F_{\rm{conv}}}}_{\rm{CZ}}}{\bar{\ensuremath{F_{\rm{conv}}}}_{\rm{PZ}}}.
\label{eqn:theory_P_defn}
\end{equation}
Since $\ensuremath{F_{\rm{conv}}} < 0$ in the PZ, the sign of $\ensuremath{\mathcal{P}}$ is positive.
Intuitively, $\ensuremath{\mathcal{P}}$ describes which terms are important in Eqn.~\ref{eqn:first_pz_parameterization}.
When $\ensuremath{\mathcal{P}} \ll 1$, the buoyancy term dominates in the PZ and dissipation can be neglected there.
When $\ensuremath{\mathcal{P}} \gg 1$, buoyancy is negligible and dissipation constrains the size of the PZ.
When $\ensuremath{\mathcal{P}} \sim 1$, both terms matter.
\editone{
    In this work, we have assumed that $\ensuremath{\mathcal{P}}$ and $\xi$ are fully independent parameters.
    We make this choice because $\ensuremath{\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 $\ensuremath{\mathcal{P}}$ increases, so too does the extent of the PZ, which likely in turn modifies the value of $\xi$).
}


\subsubsection{Case I: Discontinuous flux}
\label{sec:discontinuous_theory}
We first consider a model which satisfies
\begin{equation}
\overline{\ensuremath{F_{\rm{conv}}}}(z) = \ensuremath{F_{\rm{cz}}} \begin{cases}
1			&	z \leq L_s,\\
-\ensuremath{\mathcal{P}}_D^{-1}  & 	z > L_s 
\end{cases}.
\end{equation}
Here, $\ensuremath{F_{\rm{cz}}}$ is a constant value of flux carried in the convection zone and $\ensuremath{\mathcal{P}}_D$ is the penetration parameter (subscript D for discontinuous case).
Plugging this functional form of the flux into Eqn.~\ref{eqn:theory_fraction}, and integrating the CZ over a depth $L_{\rm{CZ}}$ below $L_s$ and the PZ over a height $\ensuremath{\delta_{\rm{p}}}$ above $L_s$, we predict
\begin{equation}
\frac{\ensuremath{\delta_{\rm{p}}}}{\ensuremath{L_{\rm{CZ}}}} = \ensuremath{\mathcal{P}}_D \frac{1 - f}{1 + \xi f \ensuremath{\mathcal{P}}_D}.
\label{eqn:discontinuous_prediction}
\end{equation}
Assuming that $f$ and $\xi$ are weak functions of $\ensuremath{\mathcal{P}}_D$, we see that, for small $\ensuremath{\mathcal{P}}_D$, the size of the penetration region is linearly proportional to $\ensuremath{\mathcal{P}}_D$, but saturates as $\ensuremath{\mathcal{P}}_D \rightarrow \infty$ due to dissipation.
Intuitively, this result makes sense: as $\ensuremath{\mathcal{P}}_D$ grows, the magnitude of $\overline{\ensuremath{F_{\rm{conv}}}}$ and the \editone{deceleration caused by} buoyancy in the PZ shrink, resulting in larger penetrative regions (but this growth cannot extend indefinitely).

\subsubsection{Case II: Piecewise linear flux}
\label{sec:linear_theory}
We next examine a model where the derivative of $\overline{\ensuremath{F_{\rm{conv}}}}(z)$ may be discontinuous at the CZ-PZ boundary,
\begin{equation}
\overline{\ensuremath{F_{\rm{conv}}}}(z) = 
\frac{\partial \ensuremath{F_{\rm{rad}}}}{\partial z}\bigg|_{\rm{CZ}}
\begin{cases}
(L_s - z) & z \leq L_s \\
-\ensuremath{\mathcal{P}}_L^{-1} (z - L_s) & z > L_s
\end{cases},
\label{eqn:linear_flux_theory}
\end{equation}
where $(\partial \ensuremath{F_{\rm{rad}}} / \partial z)|_{\rm{CZ}}$ is a constant and $\ensuremath{\mathcal{P}}_L$ is the penetration parameter (subscript L for linear case).
When $\ensuremath{\mathcal{P}}_L = 1$, $\bar{\ensuremath{F_{\rm{conv}}}}$ is a linear profile that crosses through zero at $z = L_s$.
Solving Eqn.~\ref{eqn:theory_fraction} with Eqn.~\ref{eqn:linear_flux_theory} and integrating over $L_{\rm{CZ}}$ in the CZ and $\ensuremath{\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 $\ensuremath{\delta_{\rm{p}}}$.
On that branch, we find
\begin{equation}
\frac{\ensuremath{\delta_{\rm{p}}}}{\ensuremath{L_{\rm{CZ}}}} = \sqrt{\ensuremath{\mathcal{P}}_L (1 - f)} \,\,(\sqrt{\zeta^2 + 1} - \zeta),
\label{eqn:linear_prediction}
\end{equation}
where $\zeta \equiv (\xi f/2)\sqrt{\ensuremath{\mathcal{P}}_L/(1-f)}$.
We expect the penetration height to be proportional to $\sqrt{\ensuremath{\mathcal{P}}_L}$ for small values of $\ensuremath{\mathcal{P}}_L$, and to again saturate at large values of $\ensuremath{\mathcal{P}}_L$ \editone{(as $\ensuremath{\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.~\ref{eqn:theory_fraction} through the predictions of Eqns.~\ref{eqn:discontinuous_prediction} and \ref{eqn:linear_prediction}.
Our goals are to see if the predicted scalings with the penetration parameter $\ensuremath{\mathcal{P}}$ are realized in simulations, and to measure the values of $f$ and $\xi$.

\section{Simulation Details}
We will now describe a set of simulations that test the predictions in Sec.~\ref{sec:theory}.
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 $\ensuremath{\mathcal{P}}$ and $\ensuremath{\mathcal{S}}$ can be varied separately \editone{by driving convection with internal heating, thus avoiding strongly superadiabatic boundary layers where $\ensuremath{\nabla} \rightarrow \ensuremath{\nabla_{\rm{rad}}}$}.
$\ensuremath{\mathcal{P}}$ is the ``Penetration Parameter,'' defined in Eqn.~\ref{eqn:theory_P_defn}, which compares the magnitude of the convective flux in the CZ and PZ; $\ensuremath{\mathcal{S}}$ is the ``stiffness,'' defined in Eqn.~\ref{eqn:stiffness}, and compares the buoyancy frequency in the stable radiative zone to the convective frequency. 
We suspect that some past experiments have implicitly set $\ensuremath{\mathcal{P}} \approx \ensuremath{\mathcal{S}}^{-1}$, which would result in negligible penetration for high stiffness \editone{(see discussion following Eqn.~\ref{eqn:stiffness})}.
Second, as we will show in Sec.~\ref{sec:results}, 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.

\editone{
    Appealing to the Buckingham $\pi$ theorem \citep{buckingham_1914}, we count nine fundamental input parameters in Eqns.~\ref{eqn:incompressible}-\ref{eqn:boussinesq}: $\rho_0$, $\alpha g$, $L_s$, $\nu$, $\chi$, $Q$, $\ensuremath{\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\'{e}nard convection.
    The remaining three parameters are $\ensuremath{\mathcal{S}}$, $\ensuremath{\mathcal{P}}$, and an additional parameter $\mu$, which we will hold constant and which sets the ratio between $\ensuremath{\nabla_{\rm{rad}}}$ and $\ensuremath{\nabla_{\rm{ad}}}$ in the convection zone.
}

We nondimensionalize Eqns.~\ref{eqn:incompressible}-\ref{eqn:boussinesq} on the length scale of the Schwarzschild-unstable convection zone $L_s$, the timescale of freefall across that convection zone 

\begin{equation}
    \editone{
    \tau_{\rm{ff}} = \left(\frac{L_s}{|\alpha| g Q_0}\right)^{1/3},
}
\end{equation}
and the temperature scale of the internal heating over that freefall time $\Delta T$; \editone{mass is nondimensionalized so that the freefall ram pressure $\rho_0(L_s/\tau_{\rm{ff}})^2 = 1$,}
\begin{equation}
\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\,\,\,
\vec{\nabla}^* = L_s^{-1} \vec{\nabla},\,\,\,\,\,
\\
    &\vec{u}^* = u_{\rm{ff}}\vec{u} = \frac{L_s}{\tau_{\rm{ff}}} \vec{u}, \qquad\qquad
p^* = \rho_0 u_{\rm{ff}}^2\varpi,
\\
&k^* = (L_s^2 \tau_{\rm{ff}}^{-1})k,\qquad
\ensuremath{\mathcal{R}} = \frac{u_{\rm{ff}} L_s}{\nu},\qquad
\Pran = \frac{\nu}{\chi}.
\end{split}
\end{equation}
For convenience, here we define quantities with $*$ (e.g., $T^*$) as being the ``dimensionful'' quantities of Eqns.~\ref{eqn:incompressible}-\ref{eqn:boussinesq}.
Henceforth, quantities without $*$ (e.g., $T$) are dimensionless.
The dimensionless equations of motion are
\label{sec:simulation_details}
\begin{align}
&\vec{\nabla}\dot\vec{u} = 0 
\label{eqn:nondim_incompressible} \\
&\partial_t \vec{u} + \vec{u}\dot\vec{\nabla}\vec{u} = -\vec{\nabla} \varpi + T \hat{z} + \ensuremath{\mathcal{R}}^{-1}\vec{\nabla}^2 \vec{u}
\label{eqn:nondim_momentum} \\
\begin{split}
\partial_t T + \vec{u}\dot\vec{\nabla} T + w \vec{\nabla}_{\rm{ad}}  + \vec{\nabla}\dot[-k \vec{\nabla} \overline{T}]\qquad\qquad 
\\
\qquad\qquad\qquad\qquad\qquad= (\Pran\ensuremath{\mathcal{R}})^{-1}\vec{\nabla}^2 T' + Q.
\label{eqn:nondim_temperature}
\end{split}
\end{align}
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 $\ensuremath{\nabla} = \ensuremath{\nabla_{\rm{ad}}}$ for $z \leq L_s$, and $\ensuremath{\nabla} = \ensuremath{\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 $\vec{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,
\begin{equation}
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}.
\label{eqn:sim_Q}
\end{equation}
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.
\editone{
We specify the flux using
\begin{equation}
    \mu \equiv \frac{\ensuremath{F_{\rm{bot}}}}{F_H}
    \label{eqn:mu_defn}
\end{equation}
}
and we choose $\mu = 10^{-3}$ so that most of the flux in the convection zone is carried by the convection.

\editone{
    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 $\ensuremath{\nabla_{\rm{rad}}}$ in the convection zone.
}
The average convective velocity depends on the magnitude of the convective flux, $\angles{\vec{|u|}} \approx F_H^{1/3} = (Q_{\rm{mag}}\Delta_H)^{1/3}$.
The characteristic convective frequency is ${f_{\rm{conv}} = \angles{\vec{|u|}} / L_s}$.
Empirically we find that for our choice of parameters, $\angles{\vec{|u|}} \approx 1$, so going forward we define $f_{\rm{conv}} = 1$.
The stiffness is defined,
\begin{equation}
\ensuremath{\mathcal{S}} \equiv \frac{N^2}{f_{\rm{conv}}^2} = N^2,
\label{eqn:stiffness}
\end{equation}
where $N^2$ is the Brunt-V\"{a}is\"{a}l\"{a}$\,$frequency in the radiative zone.
In our nondimensionalization, $N^2 = \ensuremath{\nabla_{\rm{ad}}} - \ensuremath{\nabla_{\rm{rad}}}$ in the radiative zone.
We use $\ensuremath{\mathcal{S}}$ as a control parameter.

\editone{
    In many prior studies, the stiffness has been set by the ratio of the subadiabaticity of $\ensuremath{\nabla_{\rm{rad}}}$ in the RZ to the superadiabaticity of $\ensuremath{\nabla_{\rm{rad}}}$ in the CZ,
    \begin{equation}
        \tilde{\ensuremath{\mathcal{S}}} = \frac{|\ensuremath{\nabla_{\rm{rad}}} - \ensuremath{\nabla_{\rm{ad}}}|_{\rm{RZ}}}{|\ensuremath{\nabla_{\rm{rad}}} - \ensuremath{\nabla_{\rm{ad}}}|_{\rm{CZ}}} = \frac{N^2}{|\ensuremath{\nabla_{\rm{rad}}} - \ensuremath{\nabla_{\rm{ad}}}|_{\rm{CZ}}}.
    \end{equation}
    In those studies, $\tilde{\ensuremath{\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.
}

\editone{
    Previous work has not defined $\ensuremath{\mathcal{P}}$, but its definition in our current study should apply to previous studies,
    \begin{equation}
        \ensuremath{\mathcal{P}} = -\frac{k_{\rm{CZ}}(\ensuremath{\nabla_{\rm{rad}}} - \ensuremath{\nabla_{\rm{ad}}})_{\rm{CZ}}}{k_{\rm{RZ}}(\ensuremath{\nabla_{\rm{rad}}} - \ensuremath{\nabla_{\rm{ad}}})_{\rm{RZ}}}.
    \end{equation}
    We note that $\ensuremath{\mathcal{P}}$ can be related to $\ensuremath{\mathcal{S}}$ and $\tilde{\ensuremath{\mathcal{S}}}$, $\ensuremath{\mathcal{P}} = (k_{\rm{CZ}}/k_{\rm{RZ}}) \tilde{\ensuremath{\mathcal{S}}}^{-1} = (k_{\rm{CZ}}/k_{\rm{RZ}}) (\ensuremath{\nabla_{\rm{rad}}} - \ensuremath{\nabla_{\rm{ad}}})_{\rm{CZ}} \ensuremath{\mathcal{S}}^{-1}$.
    Our use of internal heating to decouple convective perturbations from $\ensuremath{\nabla_{\rm{rad}}}$ in the CZ allows us to separately specify these nondimensional parameters.
    The distinction between $\ensuremath{\mathcal{S}}$ and $\ensuremath{\mathcal{P}}$ is perhaps clearer in the language of stellar evolution, where $\ensuremath{\mathcal{S}}$ is roughly the inverse square Mach number of the convection while $\ensuremath{\mathcal{P}}$ is set by the ratio of $\ensuremath{\nabla_{\rm{rad}}}$ and $\ensuremath{\nabla_{\rm{ad}}}$.
}

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

We now describe the two types of simulations conducted in this work (Case I and Case II).
We provide Fig.~\ref{fig:parameter_space} 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.


\subsection{Case I: Discontinuous flux}
\label{sec:numerics_case1}
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,
\begin{equation}
k = \begin{cases}
k_{\rm{CZ}}	&	z < 1 \\
k_{\rm{RZ}} &	z \geq 1
\end{cases},
\label{eqn:sim_discontinuous_k}
\end{equation}
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 $\ensuremath{\mathcal{S}}$ and $\ensuremath{\mathcal{P}}_D$ as inputs and specifying the radiative flux at the bottom boundary and in the RZ defines this system\footnote{
We solve the system of equations $\ensuremath{\mathcal{S}} = (\ensuremath{\nabla_{\rm{ad}}} - \ensuremath{\nabla_{\rm{rad}}})/f_{\rm{conv}}^2$, $\ensuremath{\mathcal{P}}_D = F_H/(k_{\rm{RZ}}[\ensuremath{\nabla_{\rm{ad}}} - \ensuremath{\nabla_{\rm{rad}}}])$, $F_{\rm{bot}} = k_{\rm{CZ}}\ensuremath{\nabla_{\rm{ad}}}$, and ${F_{\rm{bot}} + F_H = k_{\rm{RZ}}\ensuremath{\nabla_{\rm{rad}}}}$ to arrive at Eqns.~\ref{eqn:discontinuous_constants}.
},
\begin{equation}
\begin{split}
&k_{\rm{RZ}} = \frac{F_H}{f_{\rm{conv}}^2\ensuremath{\mathcal{S}}\ensuremath{\mathcal{P}}_D},\\
&k_{\rm{CZ}} = k_{\rm{RZ}}\frac{\mu}{1 + \mu + \ensuremath{\mathcal{P}}_D^{-1}},\\
&\ensuremath{\nabla_{\rm{ad}}} = f_{\rm{conv}}^{2}\ensuremath{\mathcal{S}}\ensuremath{\mathcal{P}}_D(1 + \mu + \ensuremath{\mathcal{P}}_D^{-1}),\\
&\ensuremath{\nabla_{\rm{rad}}} = \ensuremath{\nabla_{\rm{ad}}} - f_{\rm{conv}}^{2}\ensuremath{\mathcal{S}}.
\label{eqn:discontinuous_constants}
\end{split}
\end{equation}


\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{parameter_space.pdf}
\caption{
Each marker denotes a simulation conducted in this work in the $\ensuremath{\mathcal{R}}-\ensuremath{\mathcal{P}}$ parameter space at $\ensuremath{\mathcal{S}} = 10^3$. 
Purple circles represent Case I (Sec.~\ref{sec:numerics_case1}) simulations and orange diamonds represent Case II (Sec.~\ref{sec:numerics_case2}) simulations; empty circular markers have stress-free (SF) boundary conditions and all other simulations have no-slip boundaries.
The green ``x'' at $\ensuremath{\mathcal{P}} = 4$ and $\ensuremath{\mathcal{R}} = 400$ denotes the location in $\ensuremath{\mathcal{R}}-\ensuremath{\mathcal{P}}$ parameter space where we vary $\ensuremath{\mathcal{S}}$ in select Case I simulations.
Boxes denote the two ``landmark'' simulations.
The landmark Case I simulation has $\ensuremath{\mathcal{R}} = 400$ and $\ensuremath{\mathcal{P}}_D = 4$.
The landmark Case II simulation has $\ensuremath{\mathcal{R}} = 800$ and $\ensuremath{\mathcal{P}}_L = 4$.
Both landmark simulations have $\ensuremath{\mathcal{S}} = 10^3$ and no-slip boundary conditions.
\label{fig:parameter_space}
}
\end{figure}





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


\subsection{Case II: Piecewise linear flux}
\label{sec:numerics_case2}

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,
\begin{equation}
\partial_z k = \partial_z k_0
\begin{cases}
1	&	z < 1 \\
\ensuremath{\mathcal{P}}_L^{-1} &	z \geq 1
\end{cases}
\label{eqn:sim_linear_k}
\end{equation}
Since $k$ varies with height, formally the values of $\ensuremath{\mathcal{S}}$ and $\ensuremath{\mathcal{P}}$ also vary with height; we specify their values at $z = 2$.
By this choice\footnote{
We solve the system of equations where ${F_{\rm{bot}} = k_{\rm{bot}}\ensuremath{\nabla_{\rm{ad}}}}$, ${F_{\rm{bot}} + F_H = k_{\rm{ad}}\ensuremath{\nabla_{\rm{ad}}}}$, ${k_{\rm{ad}} = k_{\rm{bot}} + \partial_z k_0 L_s}$, ${\ensuremath{\mathcal{S}} = (\ensuremath{\nabla_{\rm{ad}}} - \ensuremath{\nabla}_{\mathrm{rad},z=2L_s})/f_{\rm{conv}}^2}$, and ${\ensuremath{\nabla_{\rm{rad}}} = F_{\rm{tot}}/k(z)}$.
}, we require
\begin{equation}
\partial_z k_0 = \frac{F_{H}}{f_{\rm{conv}}^2L_s \ensuremath{\mathcal{S}} \psi},\,\,
k_b = \frac{F_H \mu}{f_{\rm{conv}}^2\ensuremath{\mathcal{S}}\psi},\,\,
\ensuremath{\nabla_{\rm{ad}}} = f_{\rm{conv}}^2 \ensuremath{\mathcal{S}} \psi,
\end{equation}
where $\psi \equiv 1 + \ensuremath{\mathcal{P}}_L(1 + \mu)$.
We will study one sweep through $\ensuremath{\mathcal{P}}_L$ space at fixed $\ensuremath{\mathcal{R}}$ and $\ensuremath{\mathcal{S}}$ (see Fig.~\ref{fig:parameter_space}).
According to Eqn.~\ref{eqn:linear_prediction}, we expect $\delta_p \propto \ensuremath{\mathcal{P}}_L^{1/2}$.

\subsection{Numerics}
\label{sct:numerics}
We time-evolve equations \ref{eqn:nondim_incompressible}-\ref{eqn:nondim_temperature} using the Dedalus pseudospectral solver \citep{burns_etal_2020}\footnote{we use commit \href{https://github.com/DedalusProject/dedalus/commit/efb13bdaa09816dde3eee897bc2a15fc284ea2f1}{efb13bd}; the closest stable release to this commit is \href{https://github.com/DedalusProject/dedalus/releases/tag/v2.2006}{v2.2006}.} using timestepper SBDF2 \citep{wang&ruuth2008} 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 = 2 L_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 \ref{app:accelerated_evolution}.
In some simulations we start with $\bar{T} = T_0$, described above, and in others we impose an established penetrative zone in the initial state $\bar{T}$ according to Eqn.~\ref{eqn:initial_grad}.

Spectral methods with finite coefficient expansions cannot capture true discontinuities.
In order to approximate discontinuous functions such as Eqns.~\ref{eqn:sim_Q}, \ref{eqn:sim_discontinuous_k}, and \ref{eqn:sim_linear_k}, we must use smooth transitions.
We therefore define a smooth Heaviside step function,
\begin{equation}
H(z; z_0, d_w) = \frac{1}{2}\left(1 + \mathrm{erf}\left[\frac{z - z_0}{d_w}\right]\right).
\label{eqn:heaviside}
\end{equation}
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.~\ref{eqn:sim_Q} and Eqn.~\ref{eqn:sim_linear_k}, we use $d_w = 0.02$; while for Eqn.~\ref{eqn:sim_discontinuous_k} 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~\ref{app:simulation_table}.
We produce the figures in this paper using matplotlib \citep{hunter2007, mpl3.3.4}.
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 repository\footnote{\url{https://github.com/evanhanders/convective_penetration_paper}}, and in a Zenodo repository \citep{supp}.

\subsection{Penetration height measurements}
In our evolved simulations, the penetrative region has a nearly adiabatic stratification $\ensuremath{\nabla} \approx \ensuremath{\nabla_{\rm{ad}}}$.
To characterize the height of the penetrative region, we measure how drastically $\ensuremath{\nabla}$ has departed from $\ensuremath{\nabla_{\rm{ad}}}$.
We define the difference between the adiabatic and radiative gradient,
\begin{equation}
\Delta \equiv \ensuremath{\nabla_{\rm{ad}}} - \ensuremath{\nabla_{\rm{rad}}}(z).
\end{equation}
We measure penetration heights in terms of ``departure points,'' or heights at which the realized temperature gradient $\ensuremath{\nabla}$ has evolved away from the adiabatic $\ensuremath{\nabla_{\rm{ad}}}$ by some fraction $h < 1$ of $\Delta$.
Specifically,
\begin{equation}
L_s + \delta_{h} = \mathrm{max}(z) \,\,\mid\,\, \ensuremath{\nabla} > (\ensuremath{\nabla_{\rm{ad}}} - h\,\Delta).
\label{eqn:delta_p_measures}
\end{equation}
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 \citet{zahn1991}'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 \editone{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 \citep[for a discussion of alternate measurement choices, see][]{pratt_etal_2017}.


\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{lessturb_dynamics_panels.pdf}
\caption{
Temperature anomalies in vertical slices through the landmark simulations.
(top) Case I landmark ($\ensuremath{\mathcal{R}} = 400$, $\ensuremath{\mathcal{P}}_D = 4$, $\ensuremath{\mathcal{S}} = 10^3$) and (bottom) Case II landmark ($\ensuremath{\mathcal{R}} = 800$, $\ensuremath{\mathcal{P}}_L = 4$, $\ensuremath{\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.~\ref{fig:vertical_dynamics_panels}, temperature anomalies have different signs in the bulk CZ and PZ.
\label{fig:lessturb_dynamics_panels}
}
\end{figure}


\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{time_evolution.pdf}
\caption{
Time evolution of the landmark Case I simulation ($\ensuremath{\mathcal{R}} = 400$, $\ensuremath{\mathcal{P}}_D = 4$, $\ensuremath{\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.~\ref{eqn:f_defn}) and $\xi$ (bottom left panel, defined in Eqn.~\ref{eqn:xi_defn}).
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 $\ensuremath{\nabla}/\ensuremath{\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 $\ensuremath{\nabla_{\rm{ad}}}$.
The solid grey curve denotes the profile of $\ensuremath{\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.~\ref{eqn:delta_p_measures}) are plotted over the profile evolution ($\delta_{0.1}$ and $\delta_{0.9}$ as red lines, $\delta_{0.5}$ as a black line).
\label{fig:time_evolution}
}
\end{figure*}


\section{Results}
\label{sec:results}

We now describe the results of the 3D dynamical simulations described in the previous section.
Fig.~\ref{fig:vertical_dynamics_panels} 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.~\ref{sec:theory} using temporally averaged measures.

\subsection{Dynamics}
In Fig.~\ref{fig:lessturb_dynamics_panels} we display snapshots of the temperature anomalies in the two ``landmark'' simulations denoted by boxes in Fig.~\ref{fig:parameter_space}.
We display the temperature anomaly in the top panel of the Case I simulation with $\ensuremath{\mathcal{R}} = 400$, $\ensuremath{\mathcal{P}}_D = 4$, and $\ensuremath{\mathcal{S}} = 10^3$; this simulation is included in all three of our parameter space sweeps and represents the point where our ($\ensuremath{\mathcal{R}}, \ensuremath{\mathcal{P}}, \ensuremath{\mathcal{S}}$) cuts converge  in Fig.~\ref{fig:parameter_space}.
We display the temperature anomaly in the bottom panel of the Case II simulation with $\ensuremath{\mathcal{R}} = 800$, $\ensuremath{\mathcal{P}}_L = 4$, and $\ensuremath{\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.~\ref{fig:vertical_dynamics_panels} (bulk Re $\sim$ 5000).
Aside from the degree of turbulence, the dynamics are very similar in Figs.~\ref{fig:vertical_dynamics_panels} \& \ref{fig:lessturb_dynamics_panels}.
In particular, we observe that \editone{relatively} hot plumes in the CZ turn into \editone{relatively} cold plumes in the PZ (as they cross the dashed horizontal lines), and \editone{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 $\ensuremath{\nabla} \approx \ensuremath{\nabla_{\rm{ad}}} > \ensuremath{\nabla_{\rm{rad}}}$ at the base of the domain below the internal heating layer (lower dotted line in bottom panel of Fig.~\ref{fig:lessturb_dynamics_panels}).

While the landmark simulations in Fig.~\ref{fig:lessturb_dynamics_panels} are not as turbulent as the dynamics in Fig.~\ref{fig:vertical_dynamics_panels}, 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.

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{parameters_vs_p.pdf}
\caption{
Simulation results vs.~$\ensuremath{\mathcal{P}}$ for both Case I (left panels; solid purple circles in Fig.~\ref{fig:parameter_space}) and Case II (right panels; solid orange diamonds in Fig.~\ref{fig:parameter_space}).
Boxed data points denote landmark simulations from Fig.~\ref{fig:parameter_space}.
The top panels show the penetration height according to Eqn.~\ref{eqn:delta_p_measures}.
The Case I penetration heights (upper left) vary linearly with $\ensuremath{\mathcal{P}}$, in line with the prediction of Eqn.~\ref{eqn:discontinuous_prediction}.
The Case II penetration heights (upper right) vary like $\sqrt{\ensuremath{\mathcal{P}}}$, in line with the prediction of Eqn.~\ref{eqn:linear_prediction}.
In the middle panels, we measure $f$ according to Eqn.~\ref{eqn:f_defn}.
We find values of $f \in [0.6, 0.9]$, and changes in $f$ are secondary to changes in $\ensuremath{\mathcal{P}}$ for determining penetration heights.
In the bottom panels, we measure $\xi$ according to Eqn.~\ref{eqn:xi_defn}.
We find characteristic values of $\xi \in [0.5, 0.75]$, suggesting that the falloff of the $\bar{\Phi}$ in the PZ is well described by a linear function (at high $\ensuremath{\mathcal{P}}$ when $\xi \approx 1/2$), or by a cubic function (at low $\ensuremath{\mathcal{P}}$ when $\xi \approx 3/4$).
\label{fig:parameters_vs_p}
}
\end{figure*}


\subsection{Qualitative description of simulation evolution}

In Fig.~\ref{fig:time_evolution}, we show the time evolution of the landmark Case I simulation ($\ensuremath{\mathcal{R}} = 400$, $\ensuremath{\mathcal{S}} = 10^3$, and $\ensuremath{\mathcal{P}}_D = 4$) whose initial temperature profile sets $\ensuremath{\nabla} = \ensuremath{\nabla_{\rm{ad}}}$ in the convection zone ($z \lesssim 1$) and $\ensuremath{\nabla} = \ensuremath{\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.~\ref{fig:time_evolution}.
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.

\begin{figure*}[t]
\centering
\includegraphics{parameters_vs_s.pdf}
\caption{
Case I simulations results vs.~$\ensuremath{\mathcal{S}}$ at $\ensuremath{\mathcal{R}} = 400$, $\ensuremath{\mathcal{P}} = 4$.
Boxed data points denote the landmark simulation from Fig.~\ref{fig:parameter_space}.
(Left panel) Penetration heights vs.~$\ensuremath{\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.~$\ensuremath{\mathcal{S}}$.
We roughly observe a $\ensuremath{\mathcal{S}}^{-1/2}$ scaling.
\label{fig:parameters_vs_s}
}
\end{figure*}

In the right panel of Fig.~\ref{fig:time_evolution}, we plot the profile of $\ensuremath{\nabla}/\ensuremath{\nabla_{\rm{ad}}}$ in our simulation at regular time intervals, where the color of the profile corresponds to time, as in the left panels.
$\ensuremath{\nabla_{\rm{ad}}}$ is plotted as a dashed horizontal line while $\ensuremath{\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 $\ensuremath{\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 (256x64$^2$ 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 $\ensuremath{\delta_{\rm{p}}}$ vs.~time displayed in Fig.~\ref{fig:time_evolution}.
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 \editone{to evolve $\ensuremath{\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 $\ensuremath{\nabla} = \ensuremath{\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 \ref{app:accelerated_evolution}.


\subsection{Dependence on $\ensuremath{\mathcal{P}}$}
We find that the height of the penetration zone is strongly dependent on $\ensuremath{\mathcal{P}}$.
In the upper two panels of Fig.~\ref{fig:parameters_vs_p}, we plot the penetration height ($\delta_{0.1}$, $\delta_{0.5}$, $\delta_{0.9}$ from Eq.~\ref{eqn:delta_p_measures}) from Case I simulations (discontinuous $k$, upper left) and Case II simulations (discontinuous $\partial_z k$, upper right).
The fixed values of $\ensuremath{\mathcal{R}}$ and $\ensuremath{\mathcal{S}}$ are shown above these panels.
We find that the leading-order $\ensuremath{\mathcal{P}}$ scaling predictions of Eqns.~\ref{eqn:discontinuous_prediction} \& \ref{eqn:linear_prediction} describe the data \editone{well at intermediate values of $\ensuremath{\mathcal{P}}$} (orange lines).
At small values of $\ensuremath{\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 $\ensuremath{\nabla_{\rm{rad}}}$ profile in Figs.~\ref{fig:grad_profiles} \& \ref{fig:time_evolution} and Sec.~\ref{sct:numerics}).
At large values of $\ensuremath{\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.~\ref{fig:parameters_vs_p} demonstrate that that $f$ and $\xi$ are to leading order constant with $\ensuremath{\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 $\ensuremath{\mathcal{P}}$ increases.
As $\ensuremath{\mathcal{P}}$ increases, we find that CZ velocities decrease, leading to a decrease in the dissipation rate.
When $\ensuremath{\mathcal{P}}$ is small, the PZ-RZ boundary (which acts like a wall, left panel of Fig.~\ref{fig:vertical_dynamics_panels}) efficiently deflects convective velocities sideways resulting in increased bulk-CZ velocities.
As $\ensuremath{\mathcal{P}}$ grows, the velocities have access to an extended PZ in which to buoyantly \editone{decelerate} before deflection, resulting in slightly lower bulk velocities.
A similar trend of $\xi$ decreasing as $\ensuremath{\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.

\subsection{Dependence on $\ensuremath{\mathcal{S}}$}

We find that the height of the penetration zone is weakly dependent on $\ensuremath{\mathcal{S}}$.
In the left panel of Fig.~\ref{fig:parameters_vs_s}, we plot the penetration height of a few Case I simulations with $\ensuremath{\mathcal{P}}_D = 4$ and $\ensuremath{\mathcal{R}} = 400$ but with different values of $\ensuremath{\mathcal{S}}$.
The mean penetration height $\delta_{0.5}$ varies only weakly with changing $\ensuremath{\mathcal{S}}$, but that the values of $\delta_{0.1}$ and $\delta_{0.9}$ vary more strongly.
The \editone{PZ-RZ boundary layer} in which $\ensuremath{\nabla}$ changes from $\ensuremath{\nabla_{\rm{ad}}}$ to $\ensuremath{\nabla_{\rm{rad}}}$ becomes narrower as $\ensuremath{\mathcal{S}}$ increases.
To quantify this effect, we plot $\delta_{0.9} - \delta_{0.1}$ in the righthand panel of Fig.~\ref{fig:parameters_vs_s}.
We find that the width of this region varies roughly according to a $\ensuremath{\mathcal{S}}^{-1/2}$ scaling law, reminiscent of the pure-overshoot law described by \citet{korre_etal_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 large\footnote{These values are large for nonlinear simulations, but modest compared to astrophysical values. \editone{While there is observational uncertainty about the magnitude of deep convective velocities in the Sun, in the MESA model presented in Sct.~\ref{sec:solar_model}, $f_{\rm{conv}} \approx 10^{-6}$ s$^{-1}$ and $N \approx 10^{-3}$ s$^{-1}$, so $\ensuremath{\mathcal{S}} \approx 10^6$.}} values of $10^{2} \leq \ensuremath{\mathcal{S}} < 10^4$ in order to ensure $N^2 > \omega^2$ even in our highest enstrophy simulations.



\subsection{Dependence on $\ensuremath{\mathcal{R}}$}

We find that the height of the penetration zone is weakly dependent on $\ensuremath{\mathcal{R}}$.
In the upper left panel of Fig.~\ref{fig:parameters_vs_re}, we find a logarithmic decrease in the penetration height with the Reynolds number.
In order to understand how this could happen at fixed $\ensuremath{\mathcal{P}}$, we also plot the output values of $f$ (upper middle) and $\xi$ (upper right).
We find that $f$ increases with increasing $\ensuremath{\mathcal{R}}$, but is perhaps leveling off as $\ensuremath{\mathcal{R}}$ becomes large.
We find that $\xi$ does not increase strongly with $\ensuremath{\mathcal{R}}$ except for in the case of laminar simulations with $\ensuremath{\mathcal{R}} < 200$.
Eqn.~\ref{eqn:discontinuous_prediction} predicts that $\ensuremath{\delta_{\rm{p}}}$ should change at fixed $\ensuremath{\mathcal{P}}$ and $\xi$ if $f$ is changing.
In the bottom left panel, we show that the change in $\ensuremath{\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).

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{parameters_vs_re.pdf}
\caption{
(Upper left panel) Penetration heights vs.~$\ensuremath{\mathcal{R}}$ for Case I simulations (vertical cuts in Fig.~\ref{fig:parameter_space}).
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 $\ensuremath{\delta_{\rm{p}}}$ vs.~$\ensuremath{\mathcal{R}}$.
(Upper middle panel) $f$ increases with $\ensuremath{\mathcal{R}}$.
(Upper right panel) $\xi$ does not change appreciably with $\ensuremath{\mathcal{R}}$ for turbulent simulations with $\ensuremath{\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.~\ref{eqn:discontinuous_prediction}.
(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 $\ensuremath{\mathcal{R}} \rightarrow \infty$ and $\ell_{\nu} \rightarrow 0$.
Boxed data points denote the landmark simulation from Fig.~\ref{fig:parameter_space}.
\label{fig:parameters_vs_re}
}
\end{figure*}





We now examine why $f$ increases as $\ensuremath{\mathcal{R}}$ increases.
In the SF simulations, within the CZ, we can reasonably approximate $\bar{\Phi}$ as a constant $\Phi_{\rm{CZ}}$ in the bulk and zero within the viscous boundary layer,
\begin{equation}
\bar{\Phi}(z) = 
\begin{cases}
\Phi_{\rm{CZ}} 	& z > \ell_\nu \\
0				& z \leq \ell_\nu
\end{cases},
\label{eqn:cz_dissipation_model}
\end{equation}
where $\ell_\nu$ is the viscous boundary layer depth.
We have visualized a NS dissipation profile in the bottom panel of Fig.~\ref{fig:theory_profiles}; SF simulations look similar in the bulk, but drop towards zero at the bottom boundary rather than reaching a maximum.
Then, we have
\begin{equation}
\int_{\rm{CZ}} \bar{\Phi}\,dz \approx \Phi_{\rm{CZ}} \left( L_s - \ell_\nu \right),
\end{equation}
and so per Eqn.~\ref{eqn:f_defn},
\begin{equation}
f = f_\infty\left(1 - \frac{\ell_\nu}{L_s}\right),
\label{eqn:f_infty}
\end{equation}
where $f_\infty$ is the expected value of $f$ at $\ensuremath{\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.~\ref{fig:parameters_vs_re}, we find that Eqn.~\ref{eqn:f_infty} with $f_\infty = 0.755$ captures the high-$\ensuremath{\mathcal{R}}$ behavior.
To measure $\ell_\nu$, we first measure the height of the extremum of the viscous portion of the kinetic energy flux $\bar{\mathcal{F}}$ near the boundary, and take $\ell_\nu$ to be the twice that height.
We find that Eqn.~\ref{eqn:f_infty} 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.~\ref{eqn:cz_dissipation_model} is a poor model for $z \leq \ell_\nu$.
In the bottom right panel of Fig.~\ref{fig:parameters_vs_re}, we demonstrate that the depth of the viscous boundary layer follows classical scaling laws from Rayleigh-B\'{e}nard convection\footnote{
If you assume the Nusselt Number dependence on the Rayleigh number is throttled by the boundaries, Nu $\propto$ Ra$^{1/3}$ (as is frequently measured), and the Reynolds number is Re $\propto$ Ra$^{1/2}$, you retrieve Nu $\propto$ Re$^{2/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 \ensuremath{\mathcal{R}}^{-2/3}$.
} \citep{ahlers_etal_2009, goluskin2016}.
Combining these trends, we expect 
\begin{equation}
f = f_\infty(1 - C \ensuremath{\mathcal{R}}^{-2/3})
\end{equation}
for a constant $C$.
Thus as $\ensuremath{\mathcal{R}} \rightarrow \infty$, $f \rightarrow f_\infty$.

We use the fitted function of $f$ from the bottom middle panel, along with Eqn.~\ref{eqn:discontinuous_prediction}, 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.~\ref{eqn:discontinuous_prediction} to be smaller than 1.
Furthermore, the theory is derived in the limit of an instantaneous transition from $\ensuremath{\nabla_{\rm{ad}}}$ to $\ensuremath{\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 $\ensuremath{\mathcal{R}} \rightarrow \infty$ and plot this as a horizontal orange line on the upper left panel of Fig.~\ref{fig:parameters_vs_re}.
This value is coincidentally very near the value of $\delta_{0.5}$ achieved in our highest-$\ensuremath{\mathcal{R}}$ simulations.
Unfortunately, we cannot probe more turbulent simulations.
We can only run the $\ensuremath{\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.~\ref{fig:time_evolution} for similar evolution in a less turbulent, $\ensuremath{\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 $\ensuremath{\delta_{\rm{p}}}$ and $f$, $f$ and $\ell_\nu$, and $\ell_\nu$ and $\ensuremath{\mathcal{R}}$ --- all of which are effects we largely understand --- suggest that PZ heights should saturate at high $\ensuremath{\mathcal{R}}$.

In summary, we find that $\ensuremath{\delta_{\rm{p}}}$ decreases as $\ensuremath{\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.~\ref{eqn:f_infty}.
Stellar convection zones are not adjacent to hard walls\footnote{
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 $\ensuremath{\mathcal{P}} = 4$ here, we expect the simulation with $\ensuremath{\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.~\ref{eqn:f_infty}, we estimate that
\begin{equation}
f_\infty = 0.86\qquad\rm{and}\qquad
\xi = 0.6
\label{eqn:f_xi_estimates}
\end{equation}
are good first estimates for $f$ and $\xi$ when applying our theory of penetrative convection to stellar models.

\section{Testing our parameterization in a simple stellar model of the Sun}
\label{sec:solar_model}

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{mesa_profiles.pdf}
\caption{
(top) Profiles of $\ensuremath{\nabla}$ (green), $\ensuremath{\nabla_{\rm{ad}}}$ (purple), and $\ensuremath{\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.
\label{fig:mesa_profiles}
}
\end{figure}

Our simulation results present a strong case for a flux- and dissipation-based model of convective penetration, similar to those considered by \citet{zahn1991} and \citet{roxburgh1989}.
In this section, we discuss a simple stellar model of the Sun which we have created by implementing our parameterization into MESA (see Appendix~\ref{app:mesa}).
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.~\ref{eqn:theory_fraction} to spherical geometry.
To do so, we replace horizontal averages in Eqn.~\ref{eqn:integral_constraint} with integrals over latitude and longitude, and find that the relevant integral constraint contains the convective luminosity,
\begin{equation}
\int |\alpha| g L_{\rm{conv}}\,dr =   \int_V \rho_0 \Phi\,dV,
\label{eqn:mesa_integral_constraint}
\end{equation}
where $L_{\rm{conv}} = 4\pi\rho_0 r^2 \bar{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.~\ref{eqn:f_defn} and define $\xi$ similarly to Eqn.~\ref{eqn:xi_defn},
\begin{equation}
\int_{\rm{PZ}} \rho_0 \Phi \,dV = \xi \frac{V_{\rm{PZ}}}{V_{\rm{CZ}}}\int_{\rm{CZ}}\rho_0 \Phi \, dV,
\label{eqn:sphere_xi_defn}
\end{equation}
where $V_{\rm{PZ}}$ and $V_{\rm{CZ}}$ are the volumes of the PZ and CZ respectively.
Eqn.~\ref{eqn:sphere_xi_defn} generalizes Eqn.~\ref{eqn:xi_defn} outside of the assumption of a plane-parallel atmosphere.
Thus Eqn.~\ref{eqn:theory_fraction} in spherical geometry is
\begin{equation}
-\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),
\label{eqn:mesa_eqn}
\end{equation}

We implemented Eqn.~\ref{eqn:mesa_eqn} in MESA (see Appendix \ref{app:mesa} for details) and evolved a $1 M_\odot$ model to an age of 4.56 Gyr with $f = 0.86$ and $\xi = 0.6$ (Eqn.~\ref{eqn:f_xi_estimates}) to qualitatively understand how our penetration parameterization modifies a stellar model.
In the top panel of Fig.~\ref{fig:mesa_profiles} we display $\ensuremath{\nabla} \equiv d\ln T/d\ln P$ from the model which includes convective penetration.
Note that $\ensuremath{\nabla}$ (green) remains close to $\ensuremath{\nabla_{\rm{ad}}}$ (purple) below the Schwarzschild convective boundary ($\ensuremath{\nabla_{\rm{ad}}} = \ensuremath{\nabla_{\rm{rad}}}$) in a penetration zone.
After some depth $\ensuremath{\nabla} \rightarrow \ensuremath{\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.~\ref{fig:mesa_profiles}.
When a PZ is present beneath a CZ, $\ensuremath{\nabla}$ experiences a sharp jump from $\ensuremath{\nabla_{\rm{ad}}}$ to $\ensuremath{\nabla_{\rm{rad}}}$ (Fig.~\ref{fig:mesa_profiles}, top panel), resulting in an acoustic ``glitch'' in the sound speed profile.

In the model shown in Fig.~\ref{fig:mesa_profiles}, we find $H_p \approx 0.082R_\odot$ at the Schwarzschild CZ boundary, and the depth of the penetration zone in Fig.~\ref{fig:mesa_profiles} is $0.042R_\odot \sim 0.5 H_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 \citep[e.g.,][their Fig. 17]{christensen-dalsgaard_etal_2011}.
The difference $\Delta c = c_{\rm{PZ}} - c_{\rm{std}}$ that we see in this stellar model of the Sun (Fig.~\ref{fig:mesa_profiles}) 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\%)$ \citep{bergemann_serenelli_2014} and $\Delta c^2 / c^2 \approx \mathcal{O}(0.4\%)$ \citep{christensen-dalsgaard_etal_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  \citep[see Sct.~7.2.1 of][for a nice review]{basu2016}.
It is interesting, however, that the width of the PZ in Fig.~\ref{fig:mesa_profiles} is strikingly similar to the inferred width of the tachocline $(0.039 \pm 0.013)R_\odot$ that is reported by \citet{charbonneau_etal_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 \citep[see Fig.~1 in][where the Schwarzschild base of the CZ is $r/R_\odot \approx 0.712$, whereas the one in Fig.~\ref{fig:mesa_profiles} is at $r/R_{\odot} \approx 0.75$]{christensen-dalsgaard_etal_2011}.
Despite the limitations of this minimal proof of concept, Fig.~\ref{fig:mesa_profiles} 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.
\editone{
    We note briefly that the theory in e.g., Eqn.~\ref{eqn:mesa_eqn} 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.
}


\section{Discussion}
\label{sec:discussion}
In this work, we presented dynamical simulations of convective penetration, in which convection mixes ${\ensuremath{\nabla} \rightarrow \ensuremath{\nabla_{\rm{ad}}}}$ beyond the Schwarzschild boundary.
To understand these simulations, we used an integral constraint \citep[reminiscent of][]{roxburgh1989} and flux-based arguments \citep[similar to][]{zahn1991} 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.~\ref{eqn:theory_fraction} \& \ref{eqn:mesa_eqn}) 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 $\ensuremath{\delta_{\rm{p}}}$ should scale with these measurable parameters and a new flux-based ``penetration parameter'' $\ensuremath{\mathcal{P}}$.
We designed and analyzed two sets of simulations which showed good agreement with these theoretical predictions.
\editone{
    These simulations differ from past studies because we separately specify $\ensuremath{\mathcal{P}}$ and the stiffness $\ensuremath{\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 ($\ensuremath{\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.~\ref{fig:parameters_vs_p} and Eqn.~\ref{eqn:f_xi_estimates} here.
Future work should aim to understand how these parameters and the theory presented in e.g., Eqn.~\ref{eqn:mesa_eqn} 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 $\ensuremath{\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 \emph{both} occur in realistic settings, and their combined effects should be studied.
Future work should follow e.g., \citet{kapyla_etal_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 \citep{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$\,\ll1$ \citep{garaud2021}.
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, \citet{kapyla2021} 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 $\Pran$ 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 \citet{roxburgh1989}'s integral constraint with \citet{zahn1991}'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.



\begin{acknowledgments}
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.
\end{acknowledgments}


