# The shape of convection in 2D and 3D global simulations of stellar interiors

M.-G Dethero, J Pratt, D.G Vlaykov, I Baraffe, T Guillet, T Goffrey, A Le Saux, A Morison

► To cite this version:

M.-G Dethero, J Pratt, D.G Vlaykov, I Baraffe, T Guillet, et al.. The shape of convection in 2D and 3D global simulations of stellar interiors. Astronomy & Astrophysics - A&A, 2024, 692, pp.A46. 10.1051/0004-6361/202451814 . hal-04725562

HAL Id: hal-04725562

https://hal.science/hal-04725562v1

Submitted on 8 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The shape of convection in 2D and 3D global simulations of stellar interiors

M.-G. Dethero<fnref n="1" />, J. Pratt<fnref n="2" />,<fnref n="1" />,<fnref n="3" />,*, D. G. Vlaykov<fnref n="3" />, I. Baraffe<fnref n="3" />,<fnref n="4" />, T. Guillet<fnref n="3" />, T. Goffrey<fnref n="5" />, A. Le Saux<fnref n="6" />, and A. Morison<fnref n="3" />

Received 6 August 2024 / Accepted 14 September 2024

## ABSTRACT

*Context.* Theoretical descriptions of convective overshooting in stellar interiors often rely on a basic one-dimensional parameterization of the flow called the filling factor for convection. Several different definitions of the filling factor have been developed for this purpose, based on: (1) the percentage of the volume, (2) the mass flux, and (3) the convective flux that moves through the boundary.

*Aims.* We examine these definitions of the filling factor with the goal of establishing their ability to explain differences between 2D and 3D global simulations of stellar interiors that include fully compressible hydrodynamics and realistic microphysics for stars.

*Methods.* We study convection and overshooting in pairs of identical two-dimensional (2D) and three-dimensional (3D) global simulations of stars produced with MUSIC, a fully compressible, time-implicit hydrodynamics code. We examine pairs of simulations for (1) a 3 $M_\odot$ red giant star near the first dredge-up point, (2) a 1 $M_\odot$ pre-main-sequence star with a large convection zone, (3) the current sun, and (4) a 20 $M_\odot$ main-sequence star with a large convective core.

*Results.* Our calculations of the filling factor based on the volume percentage and the mass flux indicate asymmetrical convection near the surface for each star with an outer convection zone. However, near the convective boundary, convective flows achieve inward-outward symmetry for each star that we study; for 2D and 3D simulations, these filling factors are indistinguishable. A filling factor based on the convective flux is contaminated by boundary-layer-like flows, making a theoretical interpretation difficult. We present two possible new alternatives to these frequently used definitions of a filling factor, which instead compare flows at two different radial points. The first alternative is the penetration parameter of Anders et al. (2022, ApJ, 926, 169). The second alternative is a new statistic that we call the plume interaction parameter. We demonstrate that both of these parameters captures systematic differences between 2D and 3D simulations around the convective boundary.

**Key words.** convection – hydrodynamics – stars: evolution – stars: interiors

## 1. Introduction

Predicting the evolution of stars requires a prescription for the amount of mixing caused by convection at different values of the internal radius of the star. This mixing has been widely linked to the concept of a filling factor for the convective inflows (e.g., Schmitt et al. 1984; Stein & Nordlund 1989; Cattaneo et al. 1989, 1991; Toomre et al. 1990; Zahn 1991; Canuto & Dubovikov 1998; Brummell et al. 2002; Rempel 2004; Browning et al. 2004; Rogers et al. 2006; Käpylä et al. 2017; Cai 2020). Conceptually, a filling factor quantifies, based on fluid motion or heat, how much of the star is moving inward at a given radius. Several early works on stellar convection (e.g., Schmitt et al. 1984; Stein & Nordlund 1989; Hurlburt et al. 1984) have predicted that the large stratification present in stellar convection zones should produce convective motions that have a pronounced asymmetry between inflows and outflows, corresponding to a filling factor less than one-half. Indeed observations of the solar surface find

strong, coherent inflows surrounded by weaker, more diffuse outflows (e.g., as discussed in Nordlund et al. 2009). Where present, such an asymmetry would alter the shear interaction between opposing flow structures, and thus change the amount of fluid mixing. The filling factor has thus been viewed as a fundamental measure of how the stellar convection zone is structured.

Evaluating the filling factor at the boundary of convective instability measures how much fluid overshoots the convection zone and enters the radiative zone. This provides a link between the structure of convective flows in the convection zone and the depth of convective overshooting or penetration. Several works (Zahn 1991; Brummell et al. 2002) have suggested that a filling factor could play a more significant role in determining the overshooting depth than the radial velocity at the convective boundary. The convective velocities typically differ between 2D and 3D, and, as noted in these works, the convective flow structures are visually different; the filling factor has been discussed as a plausible source of difference between 2D and 3D stellar simulations (e.g., Schmitt et al. 1984; Cai 2020; Pratt et al. 2020).

<footnote n="1">¹ Department of Physics and Astronomy, Georgia State University, Atlanta, GA 30303, USA</footnote>

<footnote n="2">² Lawrence Livermore National Laboratory, 7000 East Ave, Livermore, CA 94550, USA</footnote>

<footnote n="3">³ Astrophysics, College of Engineering, Mathematics and Physical Sciences, University of Exeter, EX4 4QL Exeter, United Kingdom</footnote>

<footnote n="4">⁴ École Normale Supérieure de Lyon, CRAL (UMR CNRS 5574), Université de Lyon 1, 69007 Lyon, France</footnote>

<footnote n="5">⁵ Centre for Fusion, Space and Astrophysics, Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom</footnote>

<footnote n="6">⁶ Laboratoire de Météorologie Dynamique (IPSL), Sorbonne University, CNRS, Ecole Polytechnique, Ecole Normale Superieure, Paris, France</footnote>

<footnote>* Corresponding author; pratt34@lnl.gov</footnote>
The question of how to best define a filling factor remains open because it depends on how to accurately define the edges of structures in a convective flow. Without loss of generality, we call such structures “plumes” in this work. We also refer to convective flows moving inward toward the core as “inflows” and those moving outward toward the surface as “outflows”. Many works (e.g., Schmitt et al. 1984; Canuto & Dubovikov 1998; Brummell et al. 2002; Andrassy 2015) have defined the filling factor to be the fractional area occupied by inflows at a given radius interior to the star. These works have thus defined the edges of convective plumes using the radial velocity. We call this definition the volume-percentage filling factor. In contrast, Rempel (2004) has defined a filling factor as the fraction of the volume at a given radius with an inward mass flux; we call this the mass-flux filling factor. Zahn (1991) has defined a plume based on the convective flux, and has used the fraction of the convective flux carried by inflows to define a filling factor; we call this the convective-flux filling factor. In a similar vein, Anders et al. (2022) has proposed a “penetration parameter”, based on the change in the convective flux between the convection zone and the overshooting zone, as a predictive quantity.

These definitions of the filling factor have been evaluated in ideal box-in-a-star type simulations, which often use a moderate stratification and an ideal gas equation of state. However, they have not been examined in the kinds of global stellar simulations that we perform in this work. Our simulations are performed in a spherical shell that includes a large portion of the stellar radius and uses the stratification, temperature gradient, equation of state, and opacity tables that are extracted directly from stellar structure models accurate to the current state-of-the-art in stellar modeling. Our simulations also solve the equations for fully compressible convection; no additional assumptions are made that could impact the asymmetry of the convective flow<fnref n="1" />. We evaluate each definition of a filling factor for suitability to use in theoretical models of convection and overshooting based on two criteria, namely whether they: (1) are different for different stellar models and for 2D and 3D simulations, and (2) are correlated with the measured overshooting depth in the simulations.

This work is structured as follows. In Section 2, we discuss the global fluid simulations of stars that we use to study the filling factor. In Section 3 we describe observed differences between 2D and 3D simulations of stars. In Section 4 we describe how the volume-percentage, mass-flux, and convective-flux filling factors are calculated and present results from our simulations. We also evaluate calculations of the Anders penetration parameter. In Section 5 we present statistics based on the width and numbers of inflows and outflows. We use these diagnostics to build a new nondimensional parameter, which we call the plume interaction parameter, and we show that this parameter is systematically different in 2D and 3D simulations for the stars we examine. In Section 6, we discuss the implications of our results for future stellar evolution calculations.

## 2. Simulations

We produce simulations of stellar convection in four different stellar structures. One star that we examine is a 3 $M_\odot$ star that is

ascending the red giant branch, produced with the open-source Modules for Experiments in Stellar Astrophysics (MESA) code (Paxton et al. 2010); this star has a large outer convection zone because it is near the first dredge-up point. The second star that we examine is a one $M_\odot$ pre-main sequence star, called the young sun, produced with the Lyon stellar evolution code (Baraffe & El Eid 1991; Baraffe et al. 1997, 1998). This model has been described and examined in Pratt et al. (2016, 2017, 2020); we refer to that earlier work for additional details of the young sun, beyond the brief summary of relevant points given here. The third star is a one $M_\odot$ main-sequence star with a moderate outer convective layer, similar to the current sun, produced with the Lyon stellar evolution code. The fourth star is a 20 $M_\odot$ star at the zero-age main sequence (ZAMS) with a large convective core, produced with the Lyon stellar evolution code. These stars represent the different stages of evolution where convection occurs: pre-main sequence, main sequence, and evolved stars. They also represent different types of convection zones: outer convective envelopes that range between deep and shallow, as well as convective cores. Visualizations of the radial velocities in each of these stars are provided in Fig. 1.

In these simulations, each star is assumed to have a homogeneous chemical composition. We perform pairs of 2D and 3D Implicit Large Eddy Simulations (ILES) (Grinstein et al. 2007; Ritos et al. 2018; Margolin 2019) of these stars using the Multi-dimensional Stellar Implicit Code (MUSIC). In these pairs, the simulation volume and grid for the 3D simulation are identical to the 2D simulation except for the dimensionality, as in Pratt et al. (2020). Our simulations in this work only take convection into account; the possibility of studying additional physical effects such as rotation, a tachocline, chemical mixing, and magnetic fields is omitted from the current study.

We use the MUSIC code to solve the inviscid compressible hydrodynamic equations for density $\rho$, momentum $\rho u$, and internal energy $\rho e$:

$$ \frac{\partial}{\partial t} \rho = -\nabla \cdot (\rho u), \qquad (1) $$

$$ \frac{\partial}{\partial t}(\rho u) = -\nabla \cdot (\rho u \otimes u) - \nabla p + \rho g, \qquad (2) $$

$$ \frac{\partial}{\partial t}(\rho e) = -\nabla \cdot (\rho e u) - p \nabla \cdot u + \nabla \cdot (\chi \nabla T), \qquad (3) $$

using a second-order finite volume method, a MUSCL method (Van Leer 1997; Thornber et al. 2008) of interpolation on a staggered grid, and a van Leer flux limiter (as described in Van Leer 1974; Roe 1986; LeVeque et al. 2006). For 2D simulations, the finite volume method assumes azimuthal symmetry. Time integration in the MUSIC code is fully implicit and uses a Jacobian free Newton-Krylov (JFNK) solver (Knoll & Keyes 2004) with physics-based preconditioning (Viallet et al. 2016; Newman & Knoll 2013; Mousseau et al. 2000; Chen et al. 2014; Holod et al. 2021). The MUSIC code uses an adaptive time step, which is constrained identically for each pair of 2D and 3D simulations. MUSIC simulations use the same tabulated equation of state and opacity that are used by the 1D stellar evolution code that produced the stellar structure. In Eq. (2), $g$ is the gravitational acceleration, a spherically symmetric vector consistent with that used in the stellar evolution calculation. It is not evolved for any of the simulations with an outer convective envelope; the convective core simulation did allow the gravity to be recalculated, but changes in this gravity term were small (see Baraffe et al. 2023).

We study the pairs of MUSIC simulations in 2D and 3D described in Table 1. In all of our MUSIC simulations, the

<footnote n="1">¹ We generally refer to this combination of choices for modeling stellar interiors using a stratification and microphysics directly from stellar structure and evolution calculations, a spherical volume, and compressible hydrodynamic equations as realistic, as opposed to more idealized models.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Visualizations of radial velocity in 2D simulations (from left to right) bg2D, w2D, dcs2D, and cc2D. Outward flows are indicated in pink, while inward flows are in blue; the zero point in velocity is black. The maximum and minimum values of the color scale are defined by a radial velocity magnitude near the maximum for each simulation: bg2D (±6.5 km/s), w2D (±2.9 km/s), dcs2D (±0.63 km/s), cc2D (±1.5 km/s).</figcaption></figure>

<table><caption>**Table 1.** Parameters for compressible hydrodynamic simulations performed with MUSIC.</caption><thead><tr><th></th><th>Dim.</th><th>Evol.</th><th>M/M<sub>☉</sub></th><th>(R<sub>i</sub>, R<sub>CB</sub>, R<sub>o</sub>)/R</th><th>(θ, φ) (°)</th><th>(Δθ, Δφ) (°)</th><th>H<sub>p,CB</sub>/Δr</th><th>τ<sub>conv</sub> (10<sup>5</sup> s)</th><th>Time (τ<sub>conv</sub>)</th></tr></thead><tbody><tr><td>wide3D</td><td>3D</td><td>PMS</td><td>1</td><td>(0.21, 0.43, 1.00)</td><td>(140, 140)</td><td>(0.55, 0.55)</td><td>63</td><td>7.06 ± 0.3</td><td>3.8</td></tr><tr><td>wide2D</td><td>2D</td><td>PMS</td><td>1</td><td>(0.21, 0.43, 1.00)</td><td>(140, 0)</td><td>(0.55, -)</td><td>63</td><td>7.9 ± 0.7</td><td>104</td></tr><tr><td>wm3D</td><td>3D</td><td>PMS</td><td>1</td><td>(0.10, 0.43, 1.00)</td><td>(140, 9)</td><td>(0.14, 0.14)</td><td>258</td><td>3.91 ± 0.21</td><td>5.4</td></tr><tr><td>wm2D</td><td>2D</td><td>PMS</td><td>1</td><td>(0.10, 0.43, 1.00)</td><td>(140, 0)</td><td>(0.14, -)</td><td>258</td><td>3.76 ± 0.98</td><td>141</td></tr><tr><td>bg3D</td><td>3D</td><td>RGB</td><td>3</td><td>(0.02, 0.175, 0.90)</td><td>(140, 19)</td><td>(0.30, 0.30)</td><td>55</td><td>6.48 ± 0.07</td><td>22</td></tr><tr><td>bg2D</td><td>2D</td><td>RGB</td><td>3</td><td>(0.02, 0.175, 0.90)</td><td>(140, 0)</td><td>(0.30, -)</td><td>55</td><td>10.86 ± 2.16</td><td>163</td></tr><tr><td>cc3D</td><td>3D</td><td>MS</td><td>20</td><td>(0.194, 0.287, 0.38)</td><td>(90, 22)</td><td>(0.27, 0.52)</td><td>140</td><td>92.6 ± 7.7</td><td>3.77</td></tr><tr><td>cc2D</td><td>2D</td><td>MS</td><td>20</td><td>(0.194, 0.287, 0.38)</td><td>(90, 0)</td><td>(0.27, -)</td><td>140</td><td>36.5 ± 27.5</td><td>35.7</td></tr><tr><td>dcs3D</td><td>3D</td><td>MS</td><td>1</td><td>(0.40, 0.72, 0.97)</td><td>(120, 240)</td><td>(0.33, 0.64)</td><td>134</td><td>9.6 ± 1.0</td><td>3.9</td></tr><tr><td>dcs2D</td><td>2D</td><td>MS</td><td>1</td><td>(0.40, 0.72, 0.97)</td><td>(120, 0)</td><td>(0.33, -)</td><td>134</td><td>28.4 ± 5.5</td><td>84</td></tr></tbody></table>

**Notes.** Two different pairs of young sun simulations are included: (1) the wide2D and wide3D pair have lower radial resolution but the 3D simulation uses a large angle along the equator, while (2) w2D and w3D use a smaller equatorial angle but significantly higher radial resolution. The pair of red giant simulations is bg2D and bg3D; the pair of current sun simulations is dcs2D and dcs3D; the pair of convective core simulations is cc2D and cc3D. The simulation name, dimensionality, evolutionary state, and stellar mass *M* in units of the solar mass *M*<sub>☉</sub> are provided. The inner radius *R*<sub>i</sub> of the spherical shell, the radius of the convective boundary *R*<sub>CB</sub> determined by the Schwarzschild criterion, and the outer radius *R*<sub>o</sub> of the spherical shell for the simulation are given in units of the total radius of each star, *R*, as a triplet. The angular extent of the simulation in the polar and equatorial directions is given as (θ, φ), and the grid spacing in both angular directions is (Δθ, Δφ). The average global convective turnover time τ<sub>conv</sub> is provided as well as its standard deviation in time, and the total time span for each simulation is given in units of the convective turnover time.

compressible hydrodynamic Equations (1)–(3) are solved in a spherical shell using spherical coordinates: radius *r*, and angular variables θ and φ (in 3D). As we noted in Pratt et al. (2020), grid spacing is particularly important in determining the physics in ILES, because it is directly related to the effective numerical viscosity. In Table 1, the inner and outer radius of the spherical shell for each simulation is noted, and the radial and angular grid spacings are specified. The simulation volume and grid in the *r*–θ plane are identical for each pair of simulations. To obtain accurate statistics on overshooting, both radial resolution and a long-time series of data are critical. For that reason, it is ideal for us to examine convection in spherical shells. Herwig et al. (2023) have recently described a dipole mode that fills the entire convective core in their simulation of a 25 *M*<sub>☉</sub> main-sequence star; that mode is clearly not possible in either of the simulations

of the main-sequence star that we examine here. Nevertheless, the 2D and 3D simulations of this star both equally neglect the possibility of such a dipole mode dominating the flow. We leave further study of dipole modes in convective cores to future work.

To allow a clear comparison of the convective turnover time between 2D and 3D simulations, we define this fundamental parameter as in Pratt et al. (2016):

$$ \tau_{\text{conv}}(t) = \int_{\text{CZ}} dV \frac{H_p}{|v|} / \int_{\text{CZ}} dV. \quad (4) $$

In this expression, *H*<sub>p</sub> is the pressure scale height, and the magnitude of velocity in the denominator is calculated in either 2D or 3D, depending on the simulation. The integration covers the convection zone (CZ) and is volume-weighted using volume
element dV. The instantaneous value of $\tau_{\text{conv}}$ in Eq. (4) is averaged over time to produce the values in Table 1.

In Table 1, we introduce the quantity $H_{\text{P,GB}}/\Delta r$. This ratio shows how many grid spaces resolve each pressure scale height at the convective boundary. In this work, we often refer to the convective boundary; this is the boundary of the region of convective instability calculated by the Schwarzschild criterion, which does not evolve during the simulations in Table 1. Our highest resolution simulation pair is $\text{wm2D}/\text{wm3D}$. Simulation $\text{wm3D}$ has a grid of $r \times \theta \times \phi = 1312 \times 1024 \times 64$. Our simulations have sufficient radial resolution to produce a characteristic radial profile for velocity in 2D.

All simulations in this work are ILES, and convergence is not expected in the absolute sense that direct numerical simulations (DNS) converge. A universal shape of the velocity profiles can be observed with sufficiently small grid spacing, and the increase in the velocities becomes less as the grid spacing is progressively decreased (see also the discussion of ILES and convergence in Andrassy et al. 2024). A study of the effect of grid spacing was examined systematically for the young sun in Pratt et al. (2016), and for the current sun in Vlaykov et al. (2022). The main sequence core convection simulation was studied in Baraffe et al. (2023). Similar results are obtained for the red giant simulations. Because of the complications presented by ILES for convergence, we find the use of resolution criteria based on $H_{\text{P,GB}}/\Delta r$ to be more useful for convective overshooting than the traditional DNS-style convergence studies. Such resolution criteria allow for clear comparisons between simulations of different stars that use different grids.

All data studied here are produced during steady-state convection, a period where the time-averaged value of the total kinetic energy is well-defined and not evolving in time. Each 3D simulation includes at least 3 convective turnover times of data taken during steady-state convection; for each 2D simulation, the time span is more than $30 \tau_{\text{conv}}$. The uncertainties in the calculation of the average convective turnover time have a larger standard deviation for the 2D simulation than the 3D simulation; this is clearly impacted by the longer time series of data available for the 2D simulations (e.g., also observed in Pratt et al. 2020).

We examine simulations with two variations on energy boundary conditions that maintain the original radial profiles of density and temperature of the 1D stellar evolution model. For the wide2D/wide3D and $\text{wm2D}/\text{wm3D}$ simulations, which include the full stellar radius up to the photosphere, the surface radiates energy with the local surface temperature. In this case, the energy flux varies as $\sigma T_s^4$, where $\sigma$ is the Stefan-Boltzmann constant and $T_s(\theta, t)$ is the temperature along the surface. This boundary condition can only be effectively used when the near-surface layers are included in the simulation volume and the temperature gradient near the surface is sufficiently resolved; otherwise, it results in artificially high cooling rates. The young sun simulations have more than one grid space per pressure scale height in this region. For the bg2D/bg3D, cc2D/cc3D, and dcs2D/dcs3D simulations, which do not include the full stellar radius, we hold the energy flux and luminosity constant on the outer radial boundary, at values established by the stellar structure. For an examination of how these boundary conditions affect the dynamics, we refer to Pratt et al. (2016) and Vlaykov et al. (2022).

In all of the simulations studied in this work, we use the luminosity profiles accurate to the stellar structure models; we do not employ the tactic of luminosity boosting to shorten the convective turnover time and bring the thermal time-scale of the star closer to the convective turnover time. Luminosity boosting

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** RMS radial velocity $v_{r,RMS}$ vs. the interior radial coordinate of the star $r$ normalized by the star's radius $R$, for 2D and 3D simulations of the 3 $M_{\odot}$ red giant star. The solid and dashed lines indicate a time average, taken over the entire simulation time, of the horizontally averaged radial profile. The shaded regions represent one standard deviation above and below each time-averaged line. The radial position of the convective boundary, calculated by the Schwarzschild criterion, is indicated by a vertical black line.</figcaption></figure>

leads to a substantial reduction in computational costs and makes reaching thermal equilibrium feasible when a large enough boost factor is used. However, it can also distort the original background stratification of the star. Even if this is avoided, as discussed in Baraffe et al. (2021), luminosity boosting increases the overshooting depth, the local heating in the overshooting layer, and the shape of the spectrum of waves excited (Lecoanet et al. 2019; Le Saux et al. 2022). Our choice to use luminosities that have not been artificially boosted is motivated by the need to measure an overshooting depth for a specific star. The simulations that we study in this work are far from thermal equilibrium. However the temperature gradient is not evolving during the course of these simulations, and the statistics that we produce are meaningful.

For the velocity, we impose non-penetrative and stress-free boundary conditions in the radial directions for all simulations. The energy flux and luminosity are held constant at the inner radial boundary, at the value of the energy flux at that radius in the one-dimensional stellar evolution calculation. On the inner radial boundary of the spherical shell, we impose a constant radial derivative on the density for all simulations. At the outer radial boundary, we apply different boundary conditions that suit the local derivatives in density best. For simulations bg2D/bg3D, wide2D/wide3D, and $\text{wm2D}/\text{wm3D}$ this is a hydrostatic equilibrium boundary condition on the density that maintains hydrostatic equilibrium by assuming constant internal energy and constant radial acceleration due to gravity in the boundary cells (Grimm-Strele et al. 2015). For simulations dcs2D/dcs3D and cc2D/cc3D the outer radial boundary has a constant radial derivative imposed on the density. For simulations bg2D/bg3D, wide2D/wide3D, $\text{wm2D}/\text{wm3D}$, and dcs2D/dcs3D, we impose periodicity on all physical quantities at the boundaries in $\theta$ and $\phi$. For simulation cc2D the angular
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison of radial velocity in the 3 $M_\odot$ red giant star from (left) 2D simulation bg2D, and (right) 3D simulation bg3D. The visualization is zoomed in on a small region inside the convection zone. Outward flows are indicated in pink, while inward flows are in blue; the zero point in velocity is black. The maximum and minimum values of the color scale are defined by a radial velocity magnitude near the maximum for each simulation: bg2D (±6.5 km/s), bg3D (±2.9 km/s).</figcaption></figure>

direction is treated with reflective boundary conditions; for simulation cc3D both angular directions are treated with periodic boundary conditions.

**3. Differences between 2D and 3D stellar simulations**

For each of the pairs of simulations that we study, the root-mean-square (RMS) radial velocity profile in the 2D simulation is larger than in the 3D simulation (see for example Fig. 2). This is a common result for stellar convection also noted in Muthsam et al. (1995) and Pratt et al. (2020); the extent to which the 2D velocities are larger appears to be dependent on the stellar model examined. In addition, the radial velocities in 2D and 3D simulations look different; 3D simulations appear to be “rougher” and more small-scale structure is present. This is particularly visible at points where inflows and outflows interact (see Fig. 3).

The differences between 2D and 3D simulations also reach beyond velocity amplitudes into the structure of the flow. We examine the radial profile of the local entrophy in the $\phi$-direction, defined as the square of the $\phi$ component of the vorticity $\omega_\phi = \nabla \times \mathbf{v}|_\phi$. The local $\phi$ entrophy is larger in 2D than in 3D, with the largest differences occurring at, or near, the convective boundary (see Fig. 4). These plots reveal a different shape and structuring in the flow that occurs near convective boundaries based on the dimensionality.

We calculate the overshooting depth $\ell_{ov}$ by fitting the distribution of overshooting plumes calculated using the vertical kinetic energy flux with a generalized extreme value distribution, as described in Pratt et al. (2017). We adopt the location parameter from this fit as the overshooting depth $\ell_{ov}$. The values of $\ell_{ov}$ are provided in Table 2 for all of our simulations. For the simulation pairs bm2d/bm3D, bg2D/bg3D, and cc2D/cc3D

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Radial profiles of the local $\phi$ entrophy for 2D and 3D simulations of the 3 $M_\odot$ red giant star. The lines indicate a time average, taken over the entire simulation time, of the horizontally averaged radial profile. The shaded region represents one standard deviation above and below the time-averaged line. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

these numbers are extremely close. For the simulation pair dcs2D/dcs3D, the 2D simulation has somewhat deeper overshooting, while for the simulation pair wide2D/wide3D, the 3D simulation has somewhat deeper overshooting. Given the limited
amount of data for the 3D simulations, it is not clear whether these differences between the 2D and 3D simulations are statistically significant.

# 4. Filling factor calculations

## 4.1. Definition of the filling factor based on volume percentage

We define a volume-percentage filling factor to be the fractional volume occupied by either the inflows $\sigma_{vp,in}$ or the outflows $\sigma_{vp,out}$:

$$ \sigma_{vp,in} = \frac{V^{in}(r, t)}{V^{in}(r, t) + V^{out}(r, t)}, \qquad (5) $$

$$ \sigma_{vp,out} = \frac{V^{out}(r, t)}{V^{in}(r, t) + V^{out}(r, t)}. \qquad (6) $$

Here $V^{in}$ indicates the total volume of grid cells at a given radius that has an inward velocity, while $V^{out}$ indicates the total volume of grid cells that have an outward velocity. The natural consequence is that the sum of the filling factors of inflows and outflows must be one $\sigma_{vp,in}+\sigma_{vp,out} = 1$. Because of this relation, we generally use the notation $\sigma_{vp}$ to indicate the filling factor for the plumes moving toward the convective boundary, dropping the “in” and “out” labels. Conceptually, the volume-percentage filling factor is equating the situation where there are many small plumes with an equivalent single large inflow and single large outflow. Many works (Schmitt et al. 1984; Brummell et al. 2002; Andrássy 2015; Käpylä 2024) have used this kind of definition for a filling factor. Some have chosen to evaluate the area at the cell surfaces rather than the volume of the cell to calculate this filling factor. These two alternatives converge toward the same value as the cell size is decreased; using the volume to calculate this filling factor is convenient because MUSIC is a finite-volume code. In the work of Zahn (1991), the implication is that the filling factor is a single number, independent of radius. However other authors, including Schmitt et al. (1984), Cai (2020), Canuto & Dubovikov (1998), Käpylä (2024) have formulated a filling factor that is a function of the radial depth in the star; doing so allows us to examine ideas about non-local convection throughout the stellar radius. Based on ideas acquired from early simulations of solar convection, Canuto & Dubovikov (1998) predicted a highly asymmetric convection pattern with fast, concentrated inflows and slow, broad outflows. Based on observational data at the solar surface and the idea of a stratified star, they thus asserted that the filling factor for inflows is always less than half: $\sigma_{vp,in} < 1/2$.

A calculation of $\sigma_{vp}$ in our simulation pair bg2D/bg3D is shown in Fig. 5a. Both of these simulations evidence a $\sigma_{vp}$ close to a third near the surface of the star, in rough agreement with observations of solar surface convection (e.g., Nordlund et al. 2009). The similarity between our simulations and observations of the solar surface is striking, given that the resolution of near-surface dynamics is challenging for simulations of the stellar interior. As we examine $\sigma_{vp}$ deeper in the convective envelope, we find that it grows to approximately one-half, indicating that the convection becomes symmetric at the point that plumes are overshooting the bottom of the convection zone. This result is interesting when considered in conjunction with the kinetic energy flux (see Fig. 5b). The kinetic energy flux in all of our simulations is negative in the upper and middle parts of the convection zone. It becomes positive near the convective boundary, a result not seen in early works (e.g., Hurlburt et al.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Diagnostics for convective asymmetry: radial profiles of (a) the volume-percentage filling factor of the inward moving plumes $\sigma_{vp,in}$, and (b) the normalized kinetic energy flux, $F_{kin}$, for the red giant simulations bg2D/bg3D. The solid and dashed lines indicate time-averaged radial profiles. Shaded areas indicate one standard deviation above and below these averaged lines. A thin vertical line indicates the convective boundary determined by the Schwarzschild criterion.</figcaption></figure>

1984). These combined results indicate that flows in the convection zone are more complex than the simple picture of thin, fast-moving inward plumes. In the highly stratified deep interior, inward moving plumes can be both faster and wider than plumes in the near-surface layers.

We calculate the value of the volume-percentage filling factor for plumes moving toward the convective boundary in each simulation. We then extract the value at the convective boundary, $\sigma_{vp,CB}$, for each of our simulations; these can be found in Table 2. Across all simulations, the mean value of $\sigma_{vp,CB}$ is 0.499 and the
<figure><img src="image_7.png" /><figcaption>Fig. 6. Time-averaged volume-percentage filling factor evaluated at the convective boundary determined by the Schwarzschild criterion, $\sigma_{vp,CB}$, vs the overshooting depth in units of the pressure scale height at the convective boundary $\ell_{ov}/H_{p,CB}$ for all simulations described in Table 1. Error bars indicate one standard deviation around the time-averaged value of $\sigma_{vp,CB}$. Error bars obtained in the calculation of the overshooting depth $\ell_{ov}$ are smaller than the size of the data points.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Time-averaged mass-flux filling factor evaluated at the convective boundary determined by the Schwarzschild criterion $\sigma_{mf,CB}$ vs the overshooting depth in units of the pressure scale height at the convective boundary $\ell_{ov}/H_{p,CB}$ for all simulations described in Table 1. Error bars indicate one standard deviation around the time-averaged value of $\sigma_{mf,CB}$. Error bars obtained in the calculation of the overshooting depth $\ell_{ov}$ are smaller than the size of the data points.</figcaption></figure>

median value is 0.50. For each simulation, the value of 0.5 is within one standard deviation of the time-averaged value. We find no clear difference between the values of $\sigma_{vp,CB}$ calculated from 2D and 3D simulations. The distributions of these values, calculated at different points in time during steady state convection, strongly overlap; for some simulations, the time-averaged 2D value is slightly larger, and for others, the time-averaged 3D value is slightly larger. Given the wide range in overshooting depths that we calculate for the four different stars that we examine in this work, we find no clear correlation between the volume-percentage filling factor at the convective boundary and the overshooting length (see Fig. 6).

## 4.2. Definition of the filling factor based on mass flux

We define filling factors based on the vertical mass flux $F_{\text{mass}} = \rho v_r$ so that

$$ \sigma_{\text{mf,in}} = \frac{|F_{\text{mass}}^{\text{inflows}}|}{|F_{\text{mass}}^{\text{inflows}}| + |F_{\text{mass}}^{\text{outflows}}|}, \quad (7) $$

$$ \sigma_{\text{mf,out}} = \frac{|F_{\text{mass}}^{\text{outflows}}|}{|F_{\text{mass}}^{\text{inflows}}| + |F_{\text{mass}}^{\text{outflows}}|}. \quad (8) $$

The use of absolute values in the denominator of these ratios is necessary because, due to mass conservation, the direct sum of these fluxes is small at any point in time. We note that $\sigma_{mf,in} + \sigma_{mf,out} = 1$, and we generally use the notation $\sigma_{mf}$ to indicate the filling factor for the plumes moving toward the convective boundary. This definition of a mass-flux filling factor represents a mass-weighted rather than volume-weighted version of the filling factor.

Table 2 includes values of the mass-flux filling factor at the convective boundary, $\sigma_{mf,CB}$, for each of our simulations. Like

the volume-percentage filling factor, the mass-flux filling factor is close to one-half at this point. Fig. 7 demonstrates that, like the volume-percentage filling factor, there is no clear trend between 2D and 3D results, with strongly overlapping distributions from 2D and 3D simulations. There is also no clear correlation between the mass-flux flux filling factor and the overshooting length. For the selection of stars that we study in this work, the mass-flux filling factor provides similar information to the volume-percentage filling factor.

Canuto & Dubovikov (1998) defined a filling factor based purely on the velocity (see Eq. (49a) of that work). For the simulations examined in this work, this velocity-based filling factor produces results similar to the mass-flux filling factor.

## 4.3. Definition of the filling factor based on convective flux

The analytical work of Zahn (1991) is based on the radial profile of the convective flux

$$ \bar{F}_{\text{conv}}(r, t) = -c_p(r, t)\rho(r, t)u_r(r, \theta, t)T_1(r, \theta, t). \quad (9) $$

The bar on the right-hand side of this equation indicates an average over the horizontal directions, here indicated by $\theta$. The temperature fluctuation $T_1$ is the difference between the local temperature field $T$ and the average temperature at a given radius, angle, and time. We conveniently measure this temperature fluctuation as a deviation from an initial radial profile, because the radial profile of the temperature does not evolve during our simulations. This formulation of the convective flux neglects any changes in the density and specific heat that could be dependent on an angle.

Examining the convective flux in Eq. (9), we find that the difference between this flux in the inflows and outflows is large around the convective boundary, resulting in a characteristic
<table><caption>**Table 2.** Filling factors for convection and related diagnostics.</caption><thead><tr><th></th><th>l<sub>ov,k</sub>/H<sub>p,CB</sub></th><th>σ<sub>vp,CB</sub></th><th>σ<sub>mf,CB</sub></th><th>f<sub>z,CB</sub></th><th>P<sub>A</sub></th><th>N<sub>OL</sub></th><th>N<sub>CZ</sub></th><th>W<sub>OL</sub></th><th>W<sub>CZ</sub></th><th>σ<sub>int</sub></th></tr></thead><tbody><tr><td>wide3D</td><td>0.51</td><td>0.51</td><td>0.49</td><td>-0.11</td><td>0.25</td><td>16.0</td><td>7.6</td><td>2.8</td><td>5.3</td><td>0.53</td></tr><tr><td>wide2D</td><td>0.35</td><td>0.50</td><td>0.50</td><td>-0.10</td><td>0.10</td><td>14.3</td><td>6.4</td><td>3.4</td><td>7.5</td><td>0.45</td></tr><tr><td>wm3D</td><td>1.20</td><td>0.50</td><td>0.50</td><td>0.00</td><td>1.52</td><td>31.3</td><td>19.9</td><td>1.6</td><td>2.0</td><td>0.79</td></tr><tr><td>wm2D</td><td>1.11</td><td>0.49</td><td>0.51</td><td>-0.02</td><td>0.27</td><td>26.0</td><td>12.9</td><td>1.9</td><td>3.8</td><td>0.50</td></tr><tr><td>bg3D</td><td>0.60</td><td>0.49</td><td>0.51</td><td>0.00</td><td>0.91</td><td>14.9</td><td>10.5</td><td>3.4</td><td>4.7</td><td>0.72</td></tr><tr><td>bg2D</td><td>0.61</td><td>0.49</td><td>0.51</td><td>-0.04</td><td>0.06</td><td>12.8</td><td>3.3</td><td>3.8</td><td>16.7</td><td>0.23</td></tr><tr><td>cc3D</td><td>0.15</td><td>0.48</td><td>0.52</td><td>0.00</td><td>25.8</td><td>17.7</td><td>8.7</td><td>2.6</td><td>5.0</td><td>0.51</td></tr><tr><td>cc2D</td><td>0.13</td><td>0.46</td><td>0.54</td><td>-0.07</td><td>0.52</td><td>21.7</td><td>5.3</td><td>2.1</td><td>8.7</td><td>0.25</td></tr><tr><td>dcs3D</td><td>0.06</td><td>0.48</td><td>0.52</td><td>0.00</td><td>1.88</td><td>28.4</td><td>16.0</td><td>1.5</td><td>2.3</td><td>0.67</td></tr><tr><td>dcs2D</td><td>0.14</td><td>0.45</td><td>0.55</td><td>-0.02</td><td>0.44</td><td>23.2</td><td>8.1</td><td>2.1</td><td>6.0</td><td>0.36</td></tr></tbody></table>

**Notes.** The table includes the overshooting depth $l_{ov}$ in units of the pressure scale height at the convective boundary, the volume-percentage filling factor $\sigma_{vp,CB}$, the mass-flux filling factor $\sigma_{mf,CB}$, and the incompressible convective-flux filling factor $f_{z,CB}$. The subscript CB indicates that the quantity is evaluated at the convective boundary (CB), as defined by the Schwarzschild criterion. Several additional quantities that are evaluated above or below this convective boundary are also shown. This includes the Anders penetration parameter $P_A$ and the plume interaction parameter $\sigma_{int}$. It also includes for convective envelopes (convective cores) the average number of inflows (outflows) in the overshooting layer $N_{OL}$ and in the convection zone $N_{CZ}$, and the average width of inflows (outflows) in the overshooting layer $W_{OL}$ and in the convection zone $W_{CZ}$. The average widths are displayed in units of the percentage of the simulation volume at the given radius.

negative peak when the total convective flux is calculated (see Fig. 8). This is consistent with earlier works; both Brun et al. (2011) and Browning et al. (2004) defined the top of the overshooting layer as the point where the time average of the convective flux first changes sign. They then defined the bottom of the overshooting layer as the radius in the radiative zone where the convective flux becomes small. As is clear in the figure, using this definition of an overshooting layer along with our fully compressible simulations would mean defining a sizable portion of the convection zone as part of the overshooting layer. This shape of the convective flux is not the result of changes in the temperature gradient of the star; the average temperature profile has not evolved during the course of our simulations. The large extent of the shaded area in the figure also demonstrates that the time-variation in the convective flux is substantial. These characteristics are universally present for the different kinds of stars we have simulated.

In their derivation of a filling factor based on the convective flux, Zahn (1991) defined two functions that are “horizontal structures” of the fluid $h_1(\theta, t)$ and $h_2(\theta, t)$ such that

$$u_r(r, \theta, t) = u_{r,RMS}^{\text{in}}(r, t)h_1(\theta, t), \quad (10)$$

$$T_1(r, \theta, t) = T_{1,RMS}^{\text{in}}(r, t)h_2(\theta, t). \quad (11)$$

Here $T^{\text{in}}$ and $u^{\text{in}}$ are the temperature fluctuations and radial velocity only in the volume where the fluid is moving radially inward. We derive the convective-flux filling factor for the inflows; the corresponding equations for outflows are fully analogous. The RMS operation takes the average in the horizontal directions: over $\theta$ if the simulation is 2D, and over both $\theta$ and $\phi$ if the simulation is 3D. Zahn (1991) assumed that $h_1 = h_2 = h$. To satisfy this assumption requires that the temperature fluctuation and radial velocity are strongly correlated for inflows, or

$$\frac{u_r(r, \theta, t)}{u_{r,RMS}^{\text{in}}(r, t)} = \frac{T_1(r, \theta, t)}{T_{1,RMS}^{\text{in}}(r, t)}. \quad (12)$$

If an average over a sufficiently long time period is taken, then the numerators on both sides of this equation will be zero. However, it is unclear how well this assumption is satisfied at any point in time. As the width of a plume changes with the radius,

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Radial profile of time-averaged convective flux, normalized to its maximum magnitude value, for the 3 $M_{\odot}$ red giant simulation bg2D. The shaded region represents one standard deviation in time above and below the time-averaged line. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

the ratios in Eq. (12) could change as well; this could happen if the inflows are less coherent, or if a convective boundary layer alters the flow near the bottom of the convection zone. The degree to which Eq. (12) is satisfied could also be different for different stars.

Zahn (1991) then defined the convective-flux filling factor $f_z$ as a horizontally averaged function from these horizontal structures

$$f_z(r, t) \equiv \overline{h_1(r, \theta, t)h_2(r, \theta, t)} \equiv \overline{h^2(r, \theta, t)}. \quad (13)$$

Solving Eqs. (10) and (11) for $h_1(\theta, t)$ and $h_2(\theta, t)$, we express Eq. (13) as
$$f_z(r, t) = \frac{\overline{u_r(r, \theta, t)T_1(r, \theta, t)}}{u_{r,RMS}^{\text{in}}(r, t)T_{1,RMS}^{\text{in}}(r, t)}. \quad (14)$$

This expression clarifies that the filling factor defined by Zahn (1991) is essentially a horizontal average of the convective flux, normalized by a proxy for the convective flux in the inflows. The formula in Eq. (14) is in a convenient form for direct calculations. Combining this expression for $f_z$ with Zahn’s expression for incompressible convective flux in Eq. (9), we represent the convective flux in terms of the convective flux filling factor

$$\overline{F}_{\text{conv}}(r, t) = -c_p(r, t)\rho(r, t)u_{r,RMS}^{\text{in}}(r, t)T_{1,RMS}^{\text{in}}(r, t)f_z(t). \quad (15)$$

We find that $f_z(t)$ is highly variable in time, and that a long-time average is required to produce a smooth profile; this is unsurprising because this filling factor is related to $F_{\text{conv}}$, which also requires a long-time average to converge. Our formulation of the convective-flux filling factor is dependent on radius, a departure from the formula written by Zahn (1991), where the horizontal structure functions are not radially dependent.

A characteristic result for the time-averaged radial profile of the convective-flux filling factor is shown in Fig. 9. Far from the convective boundary, this filling factor is larger and positive. However, approaching the convective boundary, the convective-flux filling factor becomes small and/or negative. This appears to be related to the negative peak in the convective flux at the convective boundary. The implication is that the complex flow patterns, and thus complex fluxes that occur around this boundary contaminate the convective-flux filling factor. A full description of those flows is beyond the scope of the present work, but will be pursued in the future. For a filling factor, which is conceptualized as a percentage of the flow moving inward, a negative number makes little sense. Nevertheless, we have also documented the time-averaged value of $f_z$ at the Schwarzschild boundary for convective instability for all of the simulations examined in this work in Table 2.

### 4.4. Eliminating compressibility as a source of error for the convective flux filling factor

One possible source of error in the convective-flux filling factor formulated by Zahn (1991) is that the assumption of incompressibility could be unphysical around the convective boundary of a star; we therefore expand on Zahn’s incompressible definition of a filling factor. A definition of the horizontally-averaged convective flux that includes compressibility is

$$\overline{F}_{\text{conv}}(r, t) = -\overline{c_p(r, \theta, t)\rho(r, \theta, t)u_r(r, \theta, t)T_1(r, \theta, t)}. \quad (16)$$

We note that Käpylä et al. (2017) use a similar definition of this flux where the horizontal average does not include the specific heat. In this case, we define four horizontal structures such that

$$u_r(r, \theta, t) = u_{r,RMS}^{\text{in}}(r, t)h_1(\theta, t), \quad (17)$$

$$T_1(r, \theta, t) = T_{1,RMS}^{\text{in}}(r, t)h_2(\theta, t), \quad (18)$$

$$\rho(r, \theta, t) = \rho_{\text{RMS}}^{\text{in}}(r, t)h_3(\theta, t), \quad (19)$$

$$c_p(r, \theta, t) = c_{p,RMS}^{\text{in}}(r, t)h_4(\theta, t). \quad (20)$$

As in the incompressible case, the superscript indicates that the RMS only includes contributions where the fluid is moving radially inward, the RMS operation takes the average in the horizontal direction, and we derive the filling factor for inward plumes.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Radial profile of the time-averaged filling factor $f_z$ as defined in Eq. (14) for 2D and 3D simulations of the 3 M$_{\odot}$ red giant star. The shaded regions represent one standard deviation above and below each time-averaged line. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

We then define a compressible-convective-flux filling factor

$$f_{\text{comp}}(t) = \overline{h_1(\theta, t)h_2(\theta, t)h_3(\theta, t)h_4(\theta, t)}, \quad (21)$$

$$= \frac{\overline{c_p(r, \theta, t)\rho(r, \theta, t)u_r(r, \theta, t)T_1(r, \theta, t)}}{c_{p,RMS}^{\text{in}}(r, t)\rho_{\text{RMS}}^{\text{in}}(r, t)u_{r,RMS}^{\text{in}}(r, t)T_{1,RMS}^{\text{in}}(r, t)}. \quad (22)$$

We calculate $f_{\text{comp}}$ for all of the simulations studied in this work. The absolute difference between $f_z$ and $f_{\text{comp}}$ is small near the convective boundary in all of our simulations (see Fig. 10). The compressible filling factor $f_{\text{comp}}$ is slightly larger near the surface in both 2D and 3D, where Mach numbers are expected to be slightly larger; this is typical of all of our simulations. Including the effect of compressibility does not resolve the problematic aspects of the convective-flux filling factor near convective boundaries.

### 4.5. The Anders penetration parameter

Although it is not strictly a filling factor, we also examine a nondimensional number formulated by Anders et al. (2022), called the “penetration parameter”. Similar to Zahn (1991), their motivation is to predict convective overshooting and penetration using a ratio of convective fluxes. It is therefore fitting to examine it here as a comparison to the filling factor in our realistic global simulations of stars. The Anders penetration parameter $P_A$ is defined as the ratio of the time-averaged convective flux on either side of the convective boundary, namely

$$P_A = -\frac{\hat{F}_{\text{conv}|\text{CZ}}}{\hat{F}_{\text{conv}|\text{OL}}}. \quad (23)$$

In Anders et al. (2022), $\hat{F}_{\text{conv}|\text{CZ}}$ is described as the time-averaged convective flux evaluated “slightly” inside the convection zone where both outflows and inflows exist, and $\hat{F}_{\text{conv}|\text{OL}}$ is
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Radial profile of the absolute difference $|f_z - f_{\text{comp}}|$ between the incompressible convective-flux filling factor and the compressible convective-flux filling factor for 2D and 3D simulations of the 3 $M_\odot$ red giant star. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

the time-averaged convective flux located “slightly” inside the radiative zone, in a layer where only inflows exist. In global simulations of stars, this language is not specific enough to determine where the convective flux should be evaluated. In addition, there is the complication of the complex boundary-layer-like flows that we find inside the convective boundary. To calculate the Anders penetration parameter across global simulations of different stars, we define $\hat{F}_{\text{conv}|\text{CZ}}$ as the maximum time-averaged value of $F_{\text{conv}}(r)$ in the convective zone, a point that is unaffected by boundary layer flows. Similarly, we define $\hat{F}_{\text{conv}|\text{OL}}$ as the minimum time-averaged value of $F_{\text{conv}}(r)$, which can be found near the convective boundary, in the overshooting layer. Values of $P_A$ for each of our simulations, calculated in this way, are shown in Table 2. The Anders penetration parameter has the useful characteristic of producing larger average values in 2D than in 3D for all of our simulation pairs (see Fig. 11). We use the standard deviation in time of $\hat{F}_{\text{conv}|\text{CZ}}$ and $\hat{F}_{\text{conv}|\text{OL}}$ to calculate uncertainties for the penetration parameter; the error bars are large, in most cases overlapping data points from both 2D and 3D simulations. Whether the differences between 2D and 3D simulations that we observe in the Anders penetration parameter can explain differences in the overshooting depth in 2D and 3D, however, remains unclear (see Fig. 11). A longer time period of data may also be required to successfully evaluate $P_A$.

Another useful result is that the Anders penetration parameter clearly produces different values for different stellar models. For the simulation pair cc2D/cc3D, which has small overshooting lengths, the minimum of the convective flux in the overshooting layer is also a small value; consequently, the value of $P_A$ becomes large, exceeding 20 for the 3D simulation. This large value, in addition to the large error bars shown in Fig. 11, is a result of how the ratio is formulated. Using $1/P_A$ would tend to produce values less than one, as the formulations of filling factors did.

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Anders penetration parameter $P_A$ vs the overshooting depth in units of the pressure scale height at the convective boundary determined by the Schwarzschild criterion $\ell_{\text{ov}}/H_{\text{p,CB}}$ for all 2D and 3D simulations described in Table 1. Error bars consider one standard deviation of the convective flux at each of the two radial points that contribute to this number. A logarithmic scale is applied to the vertical axis because this parameter, as formulated, can be larger than one.</figcaption></figure>

# 5. Plume interactions and overshooting

## 5.1. The width of inflows

The essence of the filling factor is as a diagnostic for quantitatively measuring how symmetric or asymmetric inflows and outflows are. The formulation of a filling factor in this way is linked to the early measurement that convective flows observed on the solar surface are structured into thinner regions with intense inflows and broader regions with slower outflows. We therefore consider this early idea, and we examine directly how plume widths change with radius in our simulations. As a number that sums up the widths of inflows, the filling factor is related to a low-order statistic. By instead examining the widths directly, we can pursue higher-order statistics that may be different in 2D and 3D convective flows.

As with the volume-percentage filling factor, we define a single inflow as a continuous set of cells in the $\theta$ direction, at a given radius, that all have a negative radial velocity; similarly, we define an outflow based on positive radial velocity. For our 3D simulations, we perform the same calculation for each angle $\phi$ in our grid. This allows for the 2D and 3D calculations to be directly compared; otherwise to define a two-dimensional perimeter for a convective plume would require a more involved calculation (e.g., Haller 2015; Balasuriya et al. 2018; Rempel & Chian 2023). For each of our simulations, a characteristic profile is produced for the average widths of plumes as a function of radius (see Fig. 12). Two important features emerge in each profile: a maximum average width occurs approximately in the middle of the convection zone, which we call $W_{\text{CZ}}$, and a minimum average width occurs near the bottom of the overshooting layer, which we call $W_{\text{OL}}$. In our 2D red giant simulation, we find that the difference in the width of plumes between these two points is more than 10%. Fig. 12 also demonstrates that this characteristic profile is present for
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Widths of inflowing plumes for (a) 2D and 3D simulations of the young sun, and (b) three 2D simulations of the red giant that have different grid sizes as labeled, but are otherwise identical. The shaded region represents one standard deviation above and below the time-averaged line. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

simulations with a sufficient resolution; as the grid spacing is decreased, the average width is naturally smaller. We find that there is a wider distribution of average widths in the convection zones in each of the 2D simulations that we study than in their 3D counterparts, which may be due to the larger amount of data generated in 2D.

The distribution of plume widths in the convection zone is a fundamental diagnostic pointing toward the multi-scale nature of stellar convection. In comparison, the controlled situation of Rayleigh-Bénard convection produces convection rolls that have roughly the same size, dictated by the size of the experiment.

In some thin stellar convection zones, this can also be the case; however, in a large convection zone that is defined by a significant density stratification, inflows can have a range of sizes at any given radius. The stellar simulations we examine in this work have convection zones that represent a wide range of stratifications. Our 20 $M_\odot$ convective core simulation has a density ratio of ~2 from the top to the bottom of the convection zone, the current sun simulations have a density ratio of ~55, the red giant simulations include a density ratio of ~200, and the density ratio in the young sun simulations is greater than $10^5$. Particularly for the young sun, the convection becomes a truly multi-scale flow, with plumes and convection rolls of different sizes frequently interacting with each other. This is evidenced by the examination of higher-order statistics for plume widths. The time-averaged radial profile of the skewness of inflow width is generally positive throughout the convection zone in our simulations (see for example Figs. 13 and 14), only dipping briefly into the negative for the 3D simulations, where there is a smaller time-series of data available to average. A normal distribution is defined by a skewness of zero; a positive skewness indicates that the distribution of inflowing plume widths is not symmetrical about the mean value, such that inflowing plumes larger than the average are more prevalent than smaller ones.

These figures also include the excess kurtosis, which is zero for a normal distribution. In comparison to a Gaussian distribution, a distribution with positive excess kurtosis indicates a greater prevalence of events in the wings of the distribution. Therefore, the excess kurtosis could give us an indication of the importance of different scales in a convection zone that involves many length scales. In the 2D red giant simulation bg2D shown in Fig. 13a, the excess kurtosis is positive throughout. However, the 3D simulation bg3D shown in Fig. 13b has an excess kurtosis that is both positive and negative. We find similar results for the young sun pair wl2D/wl3D; in Fig. 14b, the 3D simulation has negative excess kurtosis. This reinforces the idea that examination of the plume width can expose differences between 2D and 3D simulations.

## 5.2. Plume numbers

The width of inflows can be directly related to their number. The number of inflows does not capture information on asymmetry; however, the number of plumes decreases when inflowing plumes are wider and increases when inflowing plumes are thinner. Some theoretical predictions of overshooting have used the number of plumes (Rieutord & Zahn 1995; Pinçon et al. 2016). In addition to this relationship to the plume width, the number of plumes can be related to two possible paradigms for convection in stellar interiors (Käpylä et al. 2017; Brandenburg 2016; Spruit 1996). The first has been described as a “tree-like” structure, where the number of inflows is dependent on depth. The second has been described as a “forest-like” structure where the number of inflows is depth-independent. To connect these ideas to overshooting, and the differences between 2D and 3D convection, we examine numbers of plumes in our simulations.

We calculate the number of inflows, $N_{in}$ by counting up the continuous regions of inflowing cells at each radius; we then average over time. The number of outflows $N_{out}$ is calculated analogously. In each of our simulations with convective envelopes, we find the largest number of inflowing plumes at the surface. As radius decreases, the number of plumes continually decreases until the convective boundary. The average number of inflows increases again just beyond this convective boundary.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Time-averaged profiles of skewness and excess kurtosis of the width of inflows for the red giant simulations (a) in bg2D, and (b) bg3D. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Time-averaged profiles of skewness and excess kurtosis of the width of inflows for the young sun simulations (a) in wm2D, and (b) wm3D. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

(see Fig. 15a), a signature that could indicate the break-up or dissipation of convective flow structures or the interaction between convection and the waves that populate the radiative region. In the figure, the number of inflowing plumes increases rapidly for radii above $r/R = 0.7$; this could be described as a tree-like structure, with large plumes dominating at the bottom of this region and much smaller convective flow structures dominating at the top. Below $r/R = 0.7$, the number of plumes changes only mildly until the bottom of the convection zone; this region could be described as having a “forest-like” structure, with convective flows of similar size dominating. Following this analogy, a mildly “tree-like” structure exists in the overshooting layer, as the number of inflowing plumes increases with the depth into the

overshooting layer; to our knowledge, this has not been observed before.

In the case of core convection (see Fig. 15b) the largest number of plumes occurs at the convective boundary. Deeper in the convection core, the number of plumes decreases and then increases toward the inner radial boundary of the simulation. Because the convective core is comparatively small, convection appears to be entirely “tree-like” in this case. This structure in the number of plumes looks highly similar to the structure that we observe in shallow outer convection zones, such as the current sun.

The picture of tree-like or forest-like convection evokes the question of whether plumes split or merge. That question can

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15. Radial profile of the average number of (a) inflows for simulations of the 3 $M_\odot$ red giant star, and (b) outflows for simulations of the 20 $M_\odot$ main-sequence star with a convective core. The radial position of the convective boundary determined by the Schwarzschild criterion is indicated by a vertical black line.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 16. Nondimensional ratios (a) the plume interaction parameter, and (b) the ratio of average numbers of plumes $N_{CZ}/N_{OL}$ vs the overshooting depth in units of the pressure scale height at the convective boundary determined by the Schwarzschild criterion, $\ell_{ov}/H_{p,CB}$, for all simulations studied in this work. A logarithmic scale is applied to the vertical axis.</figcaption></figure>

only clearly be addressed from a Lagrangian point of view, and so it is beyond the scope of this work.

### 5.3. Using plume widths and numbers in a non-dimensional parameter

Inspired by our relative success in using the Anders penetration parameter to differentiate between 2D and 3D simulations, we construct non-dimensional numbers from the number of plumes and their widths, using values at two radial points. The characteristic radial profile of plume widths (see Fig. 12) indicates that there is a funneling effect on the plumes as they reach and

then pass the convective boundary. We therefore construct a non-dimensional parameter to indicate the strength of this funneling effect. Because both the widths and numbers of plumes scale with the number of shear interactions between inflows and outflows, we call this the plume interaction parameter, which we define

$$ \sigma_{\text{int}} = W_{\text{OL}}/W_{\text{CZ}}. \qquad (24) $$

For each of our simulations, the plume interaction parameter is included in Table 2; it is always less than one for our simulations. Resolution of the complex small-scale flows in the overshooting layer can be a challenge for diagnostics. For simulation
wide2D, we find $\sigma_{\text{int}} = 0.45$, while for simulation `wm2D`, where the young sun is resolved about 4 times better, $\sigma_{\text{int}} = 0.50$. This gives a clear indication that the plume interaction parameter is sufficiently resolved in our simulations, so that it changes only a small amount with increased resolution. The value of $\sigma_{\text{int}}$ is also always larger for our 3D simulations than our 2D simulations. This suggests that the plume interaction parameter can encapsulate a general difference between 2D and 3D stellar convection. Using the standard deviation of the plume widths to calculate an uncertainty for the average plume interaction parameter, we find that these uncertainties are smaller than the differences between the 2D and 3D values of the plume interaction parameter in about half of our simulation pairs. Because the plume interaction parameter is based solely on the velocity rather than the convective flux, this measure is also physically distinct from the Anders penetration parameter.

In addition to the widths of plumes, the number of plumes can be combined in a ratio to produce a second nondimensional number, $N_{\text{CZ}}/N_{\text{OL}}$. As with the plume interaction parameter, this number is meaningfully different for our 2D and 3D simulation pairs, with the 3D simulation always having a larger value. Although a ratio based on the numbers of plumes does include different information from the widths of plumes, for our simulations they appear to be remarkably similar. Because the plume interaction parameter is more clearly linked to a filling factor, we focus on that diagnostic. Fig. 16 shows how both of our nondimensional ratios relate to the overshooting length in our simulations. No clear trend emerges for these pairs of simulations. However, the stellar models that we selected here are very different; investigation of a more similar set of stellar models could more easily exhibit a trend between the plume interaction parameter and the overshooting length; we are pursuing that study.

## 6. Summary and conclusions

Because 2D simulations have higher radial velocities than 3D simulations, it has generally been assumed that 2D simulations have a larger overshooting depth. Examining convection in four different models of stars, we show that the overshooting depth is often similar in 2D and 3D simulations. Because 2D and 3D convective flows are visibly different, one possible explanation has been that a filling factor could explain these results.

In this work, we have studied different definitions of the filling factor for convection in realistic global simulations of stars to understand differences in 2D and 3D convection as well as their link to an overshooting depth. Our calculations of a filling factor based on the volume percentage or based on mass flux result in characteristic profiles. These profiles reveal that, for stars with outer convective envelopes, the inward and outward flows are highly asymmetrical near the stellar surface with a value of about one-third (for the volume-percentage filling factor) or two-thirds (for the mass-flux filling factor). However, at the bottom of the convection zone, the filling factor is about one-half, and the convection is nearly perfectly symmetrical. For our convective core simulation, these profiles are similar to those of a small convective envelope; at the convective boundary, there is symmetry between inflows and outflows. This is a significant new result for understanding stellar convection because it indicates that a filling factor calculated either using a volume percentage or the mass flux does not distinguish between 2D and 3D simulations. Nor are these filling factors, or any arguments based on asymmetrical convection, able to predict the overshooting depth for different stars.

We study a filling factor based on the convective flux, as suggested by Zahn (1991). This calculation reveals that boundary-

layer-like dynamics, as discussed by Kupka & Muthsam (2017), are significant in realistic global simulations of stars. These flows contaminate the convective flux, distorting the signal near the convective boundary, where it might be used to predict an overshooting length. Connected to the convective flux, we also examine the Anders penetration parameter. This is not a filling factor because it involves two different radial points rather than a separation of inflows and outflows at a single radial point. However, we find that this parameter can distinguish between 2D and 3D stellar convection, so it may be useful in explaining the amount of overshooting in 2D and 3D global simulations of stars, as well as for the box-in-a-star simulations where it was developed. The convective flux in the overshooting layer is particularly intermittent, making this diagnostic converge only slowly.

We proceed to examine statistics of the widths of inward-flowing plumes (for convective envelopes) and outward-flowing plumes (for core convection). We find clear statistical differences between 2D and 3D simulations in the standard deviation and the kurtosis. We identify a universal shape to the radial profile of the average plume width. We also examine the radial profile of the average number of plumes, which we link to pictures of tree-like and forest-like convection. Based on these profiles, we observe tree-like convection in overshooting layers. We construct a nondimensional number from the ratio of the average radial profile of the inward plume widths at two different radial points; we call this the plume interaction parameter. We demonstrate that the plume interaction parameter captures differences between 2D and 3D simulations.

Although both the Anders penetration parameter and our plume interaction parameter successfully and reliably indicate differences between 2D and 3D stellar convection, for the set of stars that we examine in this work, they do not clearly correlate with the overshooting length. However, this set of stars was selected because they are very different from each other: they represent stars of different sizes, with different kinds of convection zones, at different evolutionary phases. A more systematic study of how these parameters change with overshooting depth may successfully demonstrate their use in predicting overshooting; that work is underway.

**Acknowledgements.** This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. The research leading to these results is partly supported by the ERC grants 787361-COBOM and by the STFC Consolidated Grant ST/V000721/1. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. This work also used the University of Exeter local supercomputer ISCA. Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-JRNL-857704.

## References

* Anders, E. H., Jermyn, A. S., Lecoanet, D., & Brown, B. P. 2022, ApJ, 926, 169

* Andrássy, R. 2015, Ph.D. Thesis, Universiteit van Amsterdam, The Netherlands

* Andrassy, R., Leidi, G., Higl, J., et al. 2024, A&A, 683, A97

* Balasuriya, S., Ouellette, N. T., & Rypina, I. I. 2018, Phys. D: Nonlinear Phenom., 372, 31

* Baraffe, I., & El Eid, M. F. 1991, A&A, 245, 548

* Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. 1997, A&A, 327, 1054

* Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. 1998, A&A, 337, 403

* Baraffe, I., Pratt, J., Vlaykov, D., et al. 2021, A&A, 654, A126
Baraffe, I., Clarke, J., Morison, A., et al. 2023, MNRAS, 519, 5333

Brandenburg, A. 2016, ApJ, 832, 6

Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512

Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825

Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79

Cai, T. 2020, ApJ, 891, 49

Canuto, V., & Dubovikov, M. 1998, ApJ, 493, 834

Cattaneo, F., Hurlburt, N. E., & Toomre, J. 1989, Solar and Stellar Granulation (Springer), 415

Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282

Chen, G., Chacon, L., Leibs, C. A., Knoll, D. A., & Taitano, W. 2014, J. Comput. Phys., 258, 555

Grimm-Strele, H., Kupka, F., Löw-Baselli, B., et al. 2015, New Astron., 34, 278

Grinstein, F. F., Margolin, L. G., & Rider, W. J. 2007, Implicit Large Eddy Simulation (Cambridge: Cambridge University Press)

Haller, G. 2015, Annu. Rev. Fluid Mech., 47, 137

Herwig, F., Woodward, P. R., Mao, H., et al. 2023, MNRAS, 525, 1601

Holod, I., Hoelzl, M., Verma, P., et al. 2021, Plasma Phys. Control. Fusion, 63, 114002

Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1984, ApJ, 282, 557

Käpylä, P. J. 2024, A&A, 683, A221

Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23

Knoll, D. A., & Keyes, D. E. 2004, J. Comput. Phys., 193, 357

Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comput. Astrophys., 3, 1

Lecoanet, D., Cantiello, M., Quataert, E., et al. 2019, ApJ, 886, L15

Le Saux, A., Guillet, T., Baraffe, I., et al. 2022, A&A, 660, A51

LeVeque, R. J., Mihalas, D., Dorfi, E., & Müller, E. 2006, Computational Methods for Astrophysical Fluid Flow (Springer Science & Business Media), 27

Margolin, L. G. 2019, Shock Waves, 29, 27

Mousseau, V., Knoll, D., & Rider, W. 2000, J. Comput. Phys., 160, 743

Muthsam, H., Göb, W., Kupka, F., Liebich, W., & Zöchling, J. 1995, A&A, 293, 127

Newman, C., & Knoll, D. A. 2013, SIAM J. Sci. Comput., 35, S445

Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 1

Paxton, B., Bildsten, L., Dotter, A., et al. 2010, ApJS, 192, 3

Pinçon, C., Belkacem, K., & Goupil, M. 2016, A&A, 588, A122

Pratt, J., Baraffe, I., Goffrey, T., et al. 2016, A&A, 593, A121

Pratt, J., Baraffe, I., Goffrey, T., et al. 2017, A&A, 604, A125

Pratt, J., Baraffe, I., Goffrey, T., et al. 2020, A&A, 638, A15

Rempel, M. 2004, ApJ, 607, 1046

Rempel, E. L., Chian, A. C.-L., de SA Silva, S., et al. 2023, Rev. Mod. Plasma Phys., 7, 32

Rieutord, M., & Zahn, J.-P. 1995, A&A, 296, 127

Ritos, K., Kokkinakis, I. W., & Drikakis, D. 2018, Comput. Fluids, 173, 307

Roe, P. 1986, Annu. Rev. Fluid Mech., 18, 337

Rogers, T. M., Glatzmaier, G. A., & Jones, C. 2006, ApJ, 653, 765

Schmitt, J., Rosner, R., & Bohn, H. 1984, ApJ, 282, 316

Spruit, H. C. 1996, ArXiv e-prints [arXiv:astro-ph/9605020]

Stein, R., & Nordlund, Å. 1989, ApJ, 342, L95

Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., & Williams, R. 2008, J. Comput. Phys., 227, 4873

Toomre, J., Brummell, N., Cattaneo, F., & Hurlburt, N. E. 1990, Comput. Phys. Commun., 59, 105

Van Leer, B. 1974, J. Comput. Phys., 14, 361

Van Leer, B. 1997, J. Comput. Phys., 135, 229

Viallet, M., Goffrey, T., Baraffe, I., et al. 2016, A&A, 586, A153

Vlaykov, D., Baraffe, I., Constantino, T., et al. 2022, MNRAS, 514, 715

Zahn, J.-P. 1991, A&A, 252, 179