COMPARISONS OF TWO- AND THREE-DIMENSIONAL CONVECTION IN TYPE I X-RAY BURSTS

X-ray bursts (XRBs) are the thermonuclear runaway in a H/He layer on the surface of a neutron star. These transient events can be used to probe the structure of neutron stars and the equation of state (EOS) of dense material (Özel et al. 2010; Steiner et al. 2010). Furthermore, they are also the sites of rp-process nucleosynthesis (Schatz et al. 2001). For these reasons, understanding the dynamics of the explosion has seen substantial research interest in the past years.

One-dimensional studies (Taam 1980; Taam et al. 1996; Woosley et al. 2004) can reproduce the observed energies, durations, and recurrence timescales for XRBs, but use a parameterized model for convection, namely mixing length theory (which likely does not describe turbulent convection accurately, e.g., Arnett et al. 2015). An open question is whether a fully turbulent convective velocity field can modify the nucleosynthesis. Additionally, the convection may dredge up heavy element ash to the photosphere (Bhattacharyya et al. 2010; in't Zand & Weinberg 2010) thereby altering the opacity of the atmosphere, which affects the inference of neutron star mass and radius from photospheric radius expansion bursts. These are inherently three-dimensional problems.

Previously, we performed two-dimensional simulations, focusing first on pure He bursts (Malone et al. 2011), and then later on mixed H/He bursts (Malone et al. 2014). The latter study used an approximate network to capture the hot-CNO, triple-α, and initial rp-process breakout burning. There we found that we needed a spatial resolution of about 6 cm zone⁻¹ in order to accurately model the burning; for comparison, the extremely temperature-sensitive burning of the pure He models of Malone et al. (2011) required 0.5 cm zone⁻¹ resolution. In this paper, we extend our studies by performing the first three-dimensional model of convective burning in a H/He XRB, using the reaction network from Malone et al. (2014). This initial study compares to our two-dimensional results, and discusses the computational requirements for a more extensive study.

We use the publicly available⁵ MAESTRO code (Nonaka et al. 2010), which solves the equations of low Mach number hydrodynamics by reformulating the reactive Euler equations to filter soundwaves while retaining compressibility effects due to stratification and local heat release. By filtering dynamically unimportant soundwaves, MAESTRO enables efficient simulation of slow convective flows, such as those in XRBs (Malone et al. 2011, 2014), various progenitors of SNe Ia (Zingale et al. 2011, 2012; Nonaka et al. 2012), and in the cores of massive stars (Gilet et al. 2013). Also important for simulations like these is that the low Mach number formulation analytically enforces hydrostatic equilibrium of the base state, allowing us to maintain a hydrostatic atmosphere in the simulation code without the development of large spurious velocities (see, e.g., Zingale et al. 2002).

All of the MAESTRO options and microphysics used in our two-dimensional study of XRBs in Malone et al. (2014) are retained for this study. In particular, we use the new energy formulation variant of MAESTRO, based on the ideas in Klein & Pauluis (2012) and Vasil et al. (2013), which improves energy conservation and our treatment of gravity waves. We use the Helmholtz EOS from Timmes & Swesty (2000), which includes an ideal gas of nuclei, a photon gas, and an electron/positron gas with arbitrary degeneracy and relativistic parameters, and Coulomb corrections. MAESTRO is under continuous development, and we have improved the advection portion of the code since the construction of the interface states was last discussed in Almgren et al. (2008). We take the opportunity to document those changes in Appendix A.

We use the same parametrized initial model as in our two-dimensional study. Briefly, the model consists of a $M=1.4\;{M}_{\odot }$, R = 10 km neutron star, of which we model the outer $\sim 1.4\times {10}^{3}$ cm as an isothermal ($T=3\times {10}^{8}$ K), pure ${}^{56}\mathrm{Ni}$ gas. On top of the neutron star is a warm accreted layer of mainly H/He fuel that is slightly metal-rich compared to solar, with CNO metals tied up in ${}^{14}{\rm{O}}$ and ${}^{15}{\rm{O}}$ in a ratio comparable to their respective β-decay lifetimes. A smooth transition is applied between the density ($\rho =2\times {10}^{6}$ g cm⁻³) and temperature ($T=9.5\times {10}^{8}$ K) at the base of the accreted layer and the surface of the neutron star. The accreted layer is given an isentropic profile, making it convectively unstable, and the temperature decreases until a cut off temperature is reached. The original extent of the convective region is $\lesssim 2\times {10}^{3}$ cm. Figure 1 shows the density and temperature profile, along with the values of the cut-off densities that are part of the MAESTRO algorithm. For the three-dimensional simulations present here, the anelastic and base cut-off densities have been increased slightly to $2\times {10}^{3}\;{\rm{g}}\;{\mathrm{cm}}^{-3}$ to better quench the dynamics above the atmosphere. The two-dimensional simulations used the same parameters as in Malone et al. (2014). The reader is referred to the appendix of Malone et al. (2014) for more details of our model construction procedure. Finally we note that all of the problem setup files, inputs, and initial models for the runs presented here have been copied into the main MAESTRO code repository in Exec/SCIENCE/xrb_mixed/, allowing anyone to rerun these simulations.

Zoom In Zoom Out Reset image size

In this paper, we perform two three-dimensional simulations to assess the dynamics of the convective flow. We model the XRB using a plane-parallel geometry. Our wide simulation uses a uniform grid of 512 × 512 × 768 and our narrow simulation uses a grid of 256 × 256 × 768 zones, both with 6 cm zone⁻¹ spatial resolution—the same resolution used in our two-dimensional study. As the simulation evolves, the one-dimensional hydrostatic base state that MAESTRO carries is allowed to expand due to the heating, following the procedure described in Almgren et al. (2006).

2.1. Correction to the Network

In order to understand how dimensionality affects our results, we compare to updated two-dimensional calculations based on Malone et al. (2014). In particular, we use a 6 cm resolution 1024 × 768 zone calculation. Figure 2 shows the standard deviation of temperature (compared to other zones at the same height) as a function of height for the two- and three-dimensional runs, both at t = 0.02 s. The overall trend is the same for the two calculations, with the magnitude of the temperature fluctuations in the convective region (∼1400–3550 cm) $\delta T/\langle T\rangle \;\sim \;$ 10⁻³–10⁻⁴.

Zoom In Zoom Out Reset image size

Figure 3 shows the peak temperature and peak Mach number as a function of time for the runs. We see that they closely track one another, but that in the wide three-dimensional simulation there more sporatic spikes to moderate Mach number throughout the simulation. At the start of the calculation, there is always a period of transient behavior as the heating needs to set up a consistent convective velocity field, but the flow quickly settles down. For both simulations, the average Mach number after the transient is about 0.1; in the longer-duration two-dimensional case, the Mach number asymptotes to 0.1. The temperature plots all track one another well. We did not run the three-dimensional calculation as long as the two-dimensional calculation, to conserve computational resources.

Zoom In Zoom Out Reset image size

It is interesting to note how the peak Mach number translates into a timestep improvement compared to a fully compressible code. For this problem, the sound speed in the atmosphere is greater the deeper one goes into the atmosphere, but the Mach number is highest at the top of the atmosphere. As a result, the timestep increase will actually be better than the naive $1/M$ one would expect if the domain were uniform. A further complication is that, in the compressible code, when we reach the low density material at the top of the atmosphere that buffers us from the boundary, it is radiation pressure that dominates here, articially increasing the hydrodynamic soundspeed. This is a common limitation that arises from using an EOS that includes radiation instead of modeling the radiation field itself. For comparison, we started the same XRB simulation (in two dimensions) in the Castro code (Almgren et al. 2010), and the average timestep after $2.76\times {10}^{-4}$ s of evolution was ${\rm{\Delta }}{t}_{\mathrm{comp}}=2.79\times {10}^{-10}$ s. For the main three-dimensional MAESTRO calculation, the average timestep over the course of the entire simulation was ${\rm{\Delta }}{t}_{\mathrm{LM}}=1.93\times {10}^{-7}$ s—a $\sim 700\times$ improvement.

The convective velocity structure of the wide three-dimensional simulation is shown in Figure 4, highlighting the vertical velocity. These two images are representative of the flow throughout the simulation. We do not see the tight layering that was apparent in the older two-dimensional simulations (especially for narrower domains; see Figures 6 and 7, and the discussion in Section 4.2.1 of Malone et al. 2014) because of the fix to the reaction network discussed in Section 2.1. To better understand the difference in the nature of the convective flow, we need to examine the turbulent structure.

Turbulence is known to behave differently between two and three dimensions (see, e.g., Ouellette 2012). To get a feel for the turbulent nature of the convection in these simulations, we look at the kinetic energy power spectrum. Following the discussion regarding turbulence in stratified flows in Nonaka et al. (2012) and references therein, we calculate a generalized kinetic energy density spectrum as

$$\begin{eqnarray}&&{E}_{n}(k)=\displaystyle \frac{1}{{\rm{\Omega }}}{\displaystyle \int }_{{\boldsymbol{S}}({\boldsymbol{k}})}\displaystyle \frac{1}{2}{\hat{{\boldsymbol{V}}}}_{n}({\boldsymbol{k}})\cdot {\hat{{\boldsymbol{V}}}}_{n}^{\star }({\boldsymbol{k}})d{\boldsymbol{S}},\end{eqnarray} \tag{ 1 }$$

where ${\hat{{\boldsymbol{V}}}}_{n}({\boldsymbol{k}})$ is the Fourier transform of ${{\boldsymbol{V}}}_{n}({\boldsymbol{x}})={\rho }^{n}({\boldsymbol{x}})\tilde{{\boldsymbol{U}}}({\boldsymbol{x}})$ with n specifying the density weighting, ${\boldsymbol{S}}({\boldsymbol{k}})$ is the surface defined by $| {\boldsymbol{k}}| =k$, and the ⋆ denoting complex conjugation. We note that here we use $\tilde{{\boldsymbol{U}}}$, the local velocity on the grid, instead of explicitly calculating the turbulent velocity fluctuations from the full velocity field, including the base state expansion, ${\boldsymbol{U}}=\tilde{{\boldsymbol{U}}}+{w}_{0}{{\boldsymbol{e}}}_{r}$, because $\tilde{{\boldsymbol{U}}}$ is essentially the velocity perturbations on top of an otherwise hydrostatic background state (Nonaka et al. 2010). The volume, Ω, and surface element, $d{\boldsymbol{S}}$, are based on the dimensionality of the problem. The goal is to find the proper scaling of the energy density spectrum with wavenumber for both two- and three-dimensional flow.

The units of ${\hat{{\boldsymbol{V}}}}_{n}$ are [gⁿ ${\mathrm{cm}}^{1+D-3n}$ s⁻¹], where the extra power of D on the length scale comes from the integral over $d{\boldsymbol{x}}$ in the definition of the Fourier transform. In Equation (1), the integral is done in k-space, such that $d{\boldsymbol{S}}\sim {d}^{D-1}k$ with units [cm${}^{1-D}$], whereas the normalization is in real-space, so that Ω has units of [cm^D]. Upon integration of Equation (1), the dimensionality, D, drops out of the equation, and the units of the generalized kinetic energy density spectrum become [g²ⁿ ${\mathrm{cm}}^{3-6n}$ s⁻²] for both two- and three-dimensional configurations. For turbulent flows that have density variation (i.e., compressible or stratified flows), the typical Kolmogorov energy dissipation rate, $\epsilon (l)$, at a given length scale l should be weighted by the mass density (see Fleck 1983, 1996, for example): $\epsilon (l)=\rho {U}^{3}(l)/l$, which has units of [g cm⁻¹ s⁻¹]. The arguments of Nonaka et al. (2012) then apply to any dimension: the only combination of ${\epsilon }^{\alpha }{k}^{\beta }{E}_{n}(k)$ that yields a dimensionless quantity is when $\alpha =-2/3$, $n=1/3$, and $\beta =5/3$. If the physics of two-dimensional and three-dimensional turbulence were the same (this is likely not the case), then the spectrum defined in Equation (1) should scale as ${k}^{-5/3}$ for both two- and three-dimensional flows.

In evaluating Equation (1), we create equally spaced radial wavenumber bins, k_i, ranging from the smallest physical wavenumber, $1/L$, to the highest meaningful wavenumber, $1/(2{\rm{\Delta }}x)$, where L is the domain width. The Fourier transform of the kinetic energy density gives us

$$\begin{eqnarray}&&\hat{K}({k}_{x},{k}_{y},{k}_{z})=\displaystyle \frac{1}{2}({\hat{{\boldsymbol{V}}}}_{n}\cdot {\hat{{\boldsymbol{V}}}}_{n}^{\star }).\end{eqnarray} \tag{ 2 }$$

For each of the points in the three-dimensional $\hat{K}$ array, we define $| k| =\sqrt{{k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}}$ and determine which of the radial bins, k_i, this falls into and add the value of $\hat{K}$ to that bin's sum. Done this way, we are integrating up in spherical shells in k-space, using our discrete bins. The same procedure is done in two dimensions, but now we are working in the k_x–k_y plane, and are integrating up over annular regions in that plane, again defined by our discrete bins, k_i. We do not worry about the $1/{\rm{\Omega }}$ normalization, since we will normalize each spectrum such that its peak value is 1.

Figure 5 shows the power spectrum of the two-dimensional and wide three-dimensional XRB simulations at t = 0.02 s. For this analysis, we restrict the domain to just the vicinity of the convective region, including only the vertical range $1300\;\mathrm{cm}\lt z\lt 3550\;\mathrm{cm}$. For the three-dimensional case, we see that we have more than a decade in wavenumber where we achieve a ${k}^{-5/3}$ power-law scaling, indicative of Kolmogorov turbulence. We note that the region we are studying is not periodic in the vertical direction, but an FFT assumes periodicity, so the discontinuity through the vertical boundary may affect the behavior at high wavenumbers, perhaps accounting for the slow fall in the three-dimensional spectrum. For the two-dimensional case, the spectrum starts off with a ${k}^{-5/3}$ scaling, but then becomes steeper at moderate wavenumbers. Such a break in the power law scaling for two-dimensional turbulence is predicted for very idealized turbulence where the steeper part of the curve has a ${k}^{-3}$ scaling attributed to a cascade of enstrophy (e.g., Kraichnan 1967; Leith 1968; Batchelor 1969). Numerical simulation cannot achieve the idealized conditions (e.g., infinite domain and infinite Reynolds number) assumed in the ${k}^{-3}$ derivation, and sometimes achieve a steeper power law (e.g., the review by Gkioulekas & Tung 2006, and references therein). In our two-dimensional simulation, we see a moderate range in the spectrum after the break consistent with ${k}^{-3}$.

Zoom In Zoom Out Reset image size

We have also seen such a difference in scaling between two- and three-dimensional turbulence on smaller scales in reactive Rayleigh–Taylor simulations (Zingale et al. 2005), where we saw a spectrum that appeared to follow the ${k}^{-11/5}$ scaling predicted by Bolgiano–Obukhov statistics for a two-dimensional cascade (Niemeyer & Kerstein 1997). In that study, we found that a wide domain, giving more statistics, was essential to see this scaling. The difference in the scaling we observe in the present simulations suggests that there is a fundamental difference in how the cascade takes place between the two- and three-dimensional convection in XRBs.

Figure 5 also shows that there is relatively more power in small scale (higher wavenumber, k) features for the three-dimensional simulation compared to the two-dimensional calculation. This is made more explicit by looking at a colormap plot of the enstrophy density $\eta =| {\rm{\nabla }}\times \tilde{{\boldsymbol{U}}}{| }^{2}/2$, as is shown in Figure 6, where the left (right) panel shows the two-dimensional (three-dimensional) simulation at t = 0.02 s. The plot for the two-dimensional simulation is for the wide domain, but only half of the domain is shown to keep the same scale for both plots. For the three-dimensional simulation, the plot shows a slice through the center of the domain. The two-dimensional simulation plot appears to be dominated by moderate-sized vortices throughout the domain, while in three dimensions, we see structure on a much wider range of scales. This is similar to the results seen in comparisons of two- and three-dimensional simulations of novae (Kercek et al. 1998, 1999), although our two-dimensional results do not show as severe of a dominance of vortices as reported there.

Zoom In Zoom Out Reset image size

The panels of Figure 6 also show that in two-dimensional flow the convective motions penetrate deeper into the underlying neutron star than in the three-dimensional case. This can have implications for the amount of metal-rich material that can be dredged up into the atmosphere, potentially polluting the photosphere and adjusting the opacity. We leave these details for a future paper.

Accurate analysis of the convection during a thermonuclear runaway is challenging. Most convective analysis in the literature are focused on stellar convection, which reaches a steady-state. In that case, one can drop the time derivative in the energy equation and simply compare the balance of energy fluxes. A thermonuclear runaway is far from steady state. One can assume things are in a quasi-steady state over a somewhat short timescale and perform a RANS-like averaging of the energy balance, but it is not a priori clear the exact duration of this averaging timescale. Our data dumps are roughly once every eddy turnover time, which would likely not give good enough statistics for this approach. Instead, we have, as in Malone et al. (2011, 2014), focused on the adiabatic excess, Δ∇:

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\nabla }}={\rm{\nabla }}-{{\rm{\nabla }}}_{s};\quad {\rm{\nabla }}\equiv \displaystyle \frac{d\mathrm{ln}T/{dr}}{d\mathrm{ln}P/{dr}},\end{eqnarray} \tag{ 3 }$$

where the subscript s indicates the profile along an adiabat. Figure 7 shows the horizontal average of Δ∇ as a function of radius for both two- and three-dimensional simulations at t = 0.02 s. This view of the convective region confirms that the extent of the convective overshoot region is less in three dimensions than in two dimensions, as was seen in the comparison of Figure 6; in this snapshot, the average overshoot region in two dimensions is roughly 50% larger than that of the three-dimensional simulation. Furthermore, the upper boundary in the two simulations is a bit different. The two-dimensional simulation has a stronger degree of superadiabaticity, implying the thermal gradient is steeper than that of an adiabat. The three-dimensional simulation also appears to have, on average, a convectively stabilizing gradient around r = 3300 cm where the adiabatic excess becomes negative before a small overshoot region extends the convection to nearly the same distance as in the two-dimensional case.

Zoom In Zoom Out Reset image size

We described the first three-dimensional models of convective burning in an XRB. While the peak temperature and Mach number behave qualitatively the same as our two-dimensional calculations, the structure of the convective velocity field differs substantially, both in the global appearance and in the turbulent statistics. This is illustrated well by the difference in appearance of the enstrophy density, which in three dimensions shows the typical cascade to small scales. Since convective mixing is expected to distribute the synthesized nuclei throughout the atmosphere, potentially bringing some to the photosphere, modeling the convection accurately is important. Based on the differences seen between the two-dimensional and three-dimensional flows, this suggests that three-dimensional models should be the focus of our future simulation efforts.

The calculations presented here pave the way for a more detailed study of convective burning in XRBs. We plan to do a more thorough analysis of the convective flow patterns in both two and three dimensions in the future paper by including tracer particles in the flow. The tracers will help visualize the trajectory of the flow, and to build a statistical analysis of the transport during a thermonuclear runaway. The wide three-dimensional simulation used 2.8 M CPU-hours on the OLCF titan system (running with 768 MPI tasks and 16 threads per task). While MAESTRO can use AMR, in these simulations we would refine the entire convective region, so the cost savings would be small. Modifying the simulation to advect and store tracer particle information should increase the computational cost by only a few percent.

Our future calculations will push for increased realism of the reaction network. As detailed in Malone et al. (2014), the approximate network used here reasonably captures the overall energy release, but we plan to both improve the nuclear reaction network with a more clever selection of isotopes for an approximate network, and to investigate using larger networks whose integration can be accelerated using highly parallel hardware accelerators, such as GPUs or Intel Xeon Phi processors.

Thus far, we have only explored a single initial model, constructed with a simple parameterization. The real state of the accreted layer can vary, and there are two potential changes worth exploring. First, we extend the isentropic layer all the way to the surface of the model, but accretion would likely cause heating at the top of the atmosphere, which could truncate the convection region before the surface. This may prevent the convective plumes from reaching the steep density gradient at the very top of the atmosphere. Second, our base density is on the higher end of likely models. We should explore lower density models as well. The burning in that case would not be as vigorous, but our timestep should increase as the convective velocity decreases, making these simulations feasible. An initial study of the initial model variations can be done in two dimensions, and then selected models can be run in three dimensions as needed.

We also plan to push our calculations to larger scales. For the near future, however, these sort of calculations will be limited to convection-in-a-box studies. Capturing the range of length scales necessary to follow a laterally propagating burning front, while resolving the energy-generation region, is not currently possible. The complementary approach to ours are the calculations by Cavecchi et al. (2012), which used wide-aspect ratio zones and did not perform hydrodynamics in the vertical direction. Ultimately these two methods can inform one-another to build a picture of nucleosynthesis and dynamics of the burning front in XRBs.

Visualization was done with yt (Turk et al. 2011). The power spectrum calculation followed the "Making a Turbulent Kinetic Energy Power Spectrum" recipe in the yt Cookbook. The git-hashes of the codes used for the main three-dimensional simulation are MAESTRO: afb7a1479b2b... and BoxLib: 3fcc394f2774....

We thank Frank Timmes for making his equation of state publicly available. The work at Stony Brook was supported by DOE/Office of Nuclear Physics grants Nos. DE-FG02-06ER41448 and DE-FG02-87ER40317 to Stony Brook. Work at LANL was done under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. The work at LBNL was supported by the Applied Mathematics Program of the DOE Office of Advance Scientific Computing Research under U.S. Department of Energy under contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

Here we describe the details of the second-order Godunov integration schemes used to predict face and time-centered quantities in various steps of the algorithm. In our overall algorithm, there are three variations that share most of the same discretizations, with small differences that will be described below. In summary,

• 1.  
  Case 1: Construction of the advective velocities. In Steps 3 and 7 of the algorithm flowchart in Nonaka et al. (2010), we compute face- and time-centered normal velocities (i.e., we only compute u at x-faces, v at y-faces, etc.), ${{\boldsymbol{U}}}^{\mathrm{ADV},\star }$, given ${{\boldsymbol{U}}}^{n}$ and an associated source term, ${{\mathcal{S}}}_{{\boldsymbol{U}}}$.
• 2.  
  Case 2: Construction of the final velocity edge states. In Step 11 of the MAESTRO flowchart, we compute face- and time-centered velocities, ${{\boldsymbol{U}}}^{n+1/2}$, given ${{\boldsymbol{U}}}^{n}$ and an associated source term, ${{\mathcal{S}}}_{{\boldsymbol{U}}}$. This case is different from Case 1 in that we leverage the availability of the projected velocity field, ${\tilde{{\boldsymbol{U}}}}^{\mathrm{ADV}}$, during the characteristic tracing and upwinding steps. Also, at each face we need to compute all components of velocity, rather than just the normal components.
• 3.  
  Case 3: Construction of the scalar edge states. In Steps 4 and 8 of the MAESTRO flowchart, we compute face- and time-centered scalars, ${(\rho {X}_{k},\rho h)}^{n+1/2,(l)}$, given ${(\rho {X}_{k},\rho h)}^{n}$ and the associated source terms, ${{\mathcal{S}}}_{\rho {X}_{k}},{{\mathcal{S}}}_{\rho h}$. This is different from Case 2 in that the evolution equations for the scalars are in conservative form rather than advective form.

Our Godunov integration strategy is based on the piecewise parabolic method (PPM) of Colella & Woodward (1984). We modify this algorithm to account for the fact that (i) our underlying velocity field is spatially varying, (ii) we require a multidimensionally unsplit discretization, (iii) we have governing equations in both advective and conservative form,

$$\begin{eqnarray}&&\displaystyle \frac{\partial q}{\partial t}=-{\boldsymbol{U}}\cdot {\rm{\nabla }}q+{{\mathcal{S}}}_{q};\quad q=u,v,\;\mathrm{or}\;w,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial q}{\partial t}=-{\rm{\nabla }}\cdot ({\boldsymbol{U}}q)+{{\mathcal{S}}}_{q};\quad q=\rho {X}_{k}\;\mathrm{or}\;\rho h.\end{eqnarray} \tag{ 5 }$$

Which form is used for the scalars, $(\rho {X}_{k})$ or $(\rho h)$, is determined at runtime based on how we chose to bring these states to the interface. MAESTRO offers several possibilities, e.g., bringing $(\rho {X}_{k})$ to the interface as a single quantity, bringing $\rho \prime$ and X_k to the interface separately, or bringing ρ and X_k to the interface separately. Even more variation is allowed for $(\rho h)$, where we can use T instead of h for the interface prediction. We document both forms of interface reconstruction here. The full list of possible states is provided in the MAESTRO User's Guide. For the present simulations, we predict ρ and X_k separately for form $(\rho {X}_{k})$ on interfaces, and T is predicted and converted to h on the interface via the EOS.

For all cases, the idea is to use the original PPM algorithm is to predict preliminary face and time-centered states, ${q}^{1{\rm{D}}}$, using a one-dimensional advection equation for each direction d,

$$\begin{eqnarray}&&\displaystyle \frac{\partial q}{\partial t}=-{U}_{d}\displaystyle \frac{\partial q}{\partial {x}_{d}},\end{eqnarray} \tag{ 6 }$$

and then use these states in a multidimensionally unsplit discretization of the full equations of motion based on the ideas in Colella (1990) and Saltzman (1994) to compute updated face- and time-centered states. We now provide the details of our method.

A.1. Case 1

Here we compute face- and time-centered estimates of the normal velocity. We begin by computing preliminary face- and time-centered estimates of all velocity components at every face. Here, q will represent an arbitrary component of velocity ($u,v$, or w). The following developments are for the x-direction and we omit the $j\text{and}k$ subscripts for ease of exposition; the equations for the y- and z-directions are analogous. The first step is to construct a parabolic profile, q(x), in each zone following the methodology discussed in Colella & Woodward (1984)

$$\begin{eqnarray}&&q(x)={q}_{i,-}+\xi (x)({q}_{i,+}-{q}_{i,-}+{q}_{6,i}(1-\xi (x))),\end{eqnarray} \tag{ 7 }$$

where ${q}_{i,-}$ and ${q}_{i,+}$ are the limited edge values of the parabola and ${q}_{6,i}=6{q}_{i}-3({q}_{i,-}+{q}_{i,+})$ is related to the curvature. The quantity $\xi (x)$ converts x into the fraction of the zone from the left edge.

The parabolic profiles are then integrated along characteristics to get the average value swept over the high and low faces over the time step. By defining

$$\begin{eqnarray}&&\sigma =| {u}_{i}| \displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x},\end{eqnarray} \tag{ 8 }$$

we obtain:

$$\begin{eqnarray}&&{({q}_{L}^{1{\rm{D}}})}_{i+1/2}\\ =\;\left\{\begin{array}{ll}{q}_{i,+}-\displaystyle \frac{\sigma }{2}\left[({q}_{i,+}-{q}_{i,-})-\left(1-\displaystyle \frac{2}{3}\sigma \right){q}_{6,i}\right], & {u}_{i}\gt 0,\\ {q}_{i}, & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{({q}_{R}^{1{\rm{D}}})}_{i-1/2}\\ =\left\{\begin{array}{ll}{q}_{i,-}+\displaystyle \frac{\sigma }{2}\left[({q}_{i,+}-{q}_{i,-})+\left(1-\displaystyle \frac{2}{3}\sigma \right){q}_{6,i}\right], & {u}_{i}\lt 0,\\ {q}_{i}, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 10 }$$

Then, for the normal velocity components, we solve a Riemann problem to obtain the preliminary state at the face,

$$\begin{eqnarray}{u}_{i+1/2}^{1{\rm{D}}} & = & {\mathcal{R}}\left({({u}_{L}^{1{\rm{D}}})}_{i+1/2},{({u}_{R}^{1{\rm{D}}})}_{i+1/2}\right)\\ & = & \;\left\{\begin{array}{ll}0, & {({u}_{L}^{1{\rm{D}}})}_{i+1/2}\leqslant 0\;\mathrm{AND}\;{({u}_{R}^{1{\rm{D}}})}_{i+1/2}\geqslant 0,\\ {({u}_{L}^{1{\rm{D}}})}_{i+1/2}, & {({u}_{L}^{1{\rm{D}}})}_{i+1/2}+{({u}_{R}^{1{\rm{D}}})}_{i+1/2}\gt 0,\\ {({u}_{R}^{1{\rm{D}}})}_{i+1/2}, & {({u}_{L}^{1{\rm{D}}})}_{i+1/2}+{({u}_{R}^{1{\rm{D}}})}_{i+1/2}\lt 0.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

Next, we obtain the preliminary face- and time-centered transverse velocities using upwinding,

$$\begin{eqnarray}{v}_{i+1/2}^{1D} & = & {\mathcal{U}}\left({({v}_{L}^{1{\rm{D}}})}_{i+1/2},{({v}_{R}^{1{\rm{D}}})}_{i+1/2},{u}_{i+1/2}^{1{\rm{D}}}\right)\\ & = & \left\{\begin{array}{ll}{({v}_{L})}_{i+1/2}^{1{\rm{D}}}, & {u}_{i+1/2}^{1{\rm{D}}}\gt 0,\\ {({v}_{R})}_{i+1/2}^{1{\rm{D}}}, & {u}_{i+1/2}^{1{\rm{D}}}\lt 0.\end{array}\right.\end{eqnarray} \tag{ 12 }$$

We now have ${u}^{1{\rm{D}}}$ and ${v}^{1{\rm{D}}}$ at each face. The next step is compute updated face- and time-centered normal velocities by accounting for the transverse derivative and source terms we have ignored so far. In two dimensions, first compute ${({u}_{L})}_{i+1/2,j}$ and ${({u}_{R})}_{i+1/2,j}$ using, e.g.,

$$\begin{eqnarray}{({u}_{L})}_{i+1/2,j} & = & {({u}_{L}^{1{\rm{D}}})}_{i+1/2,j}-\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{v}_{i,j+1/2}^{1{\rm{D}}}+{v}_{i,j-1/2}^{1{\rm{D}}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{u}_{i,j+1/2}^{1{\rm{D}}}-{u}_{i,j-1/2}^{1{\rm{D}}}}{{\rm{\Delta }}y}\right)+\displaystyle \frac{{\rm{\Delta }}t}{2}{({{\mathcal{S}}}_{u})}_{{\boldsymbol{i}}},\end{eqnarray} \tag{ 13 }$$

followed by a Riemann solver to obtain the final face- and time-centered state,

$$\begin{eqnarray}&&{u}_{i+1/2,j}^{\mathrm{ADV},\star }={\mathcal{R}}\left({({u}_{L})}_{i+1/2,j},{({u}_{R})}_{i+1/2,j}\right).\end{eqnarray} \tag{ 14 }$$

In three dimensions, instead of Equation (13) we use

$$\begin{eqnarray}{({u}_{L})}_{i+1/2,j,k} & = & {({u}_{L}^{1{\rm{D}}})}_{i+1/2,j,k}-\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{v}_{i,j+1/2,k}^{1{\rm{D}}}+{v}_{i,j-1/2,k}^{1{\rm{D}}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{u}_{i,j+1/2,k}^{y| z}-{u}_{i,j-1/2,k}^{y| z}}{{\rm{\Delta }}y}\right)\\ & & -\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{w}_{i,j,k+1/2}^{1{\rm{D}}}+{v}_{i,j,k-1/2}^{1{\rm{D}}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{u}_{i,j,k+1/2}^{z| y}-{u}_{i,j,k-1/2}^{z| y}}{{\rm{\Delta }}z}\right)+\displaystyle \frac{{\rm{\Delta }}t}{2}{({{\mathcal{S}}}_{u})}_{{\boldsymbol{i}}},\end{eqnarray} \tag{ 15 }$$

where to account for transverse corner coupling, we compute the intermediate states as in Colella (1990) and Saltzman (1994). For example,

$$\begin{eqnarray}{({u}_{L}^{y| z})}_{i,j+1/2,k} & = & {({u}_{L})}_{i,j+1/2,k}-\displaystyle \frac{{\rm{\Delta }}t}{3}\left(\displaystyle \frac{{w}_{i,j,k+1/2}^{1{\rm{D}}}+{w}_{i,j,k-1/2}^{1{\rm{D}}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{u}_{i,j,k+1/2}^{1{\rm{D}}}-{u}_{i,j,k-1/2}^{1{\rm{D}}}}{{\rm{\Delta }}z}\right),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{u}_{i,j+1/2,k}^{y| z}={\mathcal{U}}\left({({u}_{L}^{y| z})}_{i,j+1/2,k},{({u}_{R}^{y| z})}_{i,j+1/2,k},{v}_{i,j+1/2,k}^{1{\rm{D}}}\right).\end{eqnarray} \tag{ 17 }$$

A.2. Case 2

Here we compute face and time-centered estimates of each component of velocity at every face. The details are the same as Case 1 up until Equation (8). Now we can use ${\tilde{{\boldsymbol{U}}}}^{\mathrm{ADV}}$ for characteristic tracing and upwinding whenever possible. Specifically, we define

$$\begin{eqnarray}&&{\sigma }_{\pm }=\left|{u}_{i\pm 1/2}^{\mathrm{ADV}}\right|\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x},\end{eqnarray} \tag{ 18 }$$

and compute

$$\begin{eqnarray}&&{({q}_{L}^{1{\rm{D}}})}_{i+1/2}\\ &&=\left\{\begin{array}{l}{q}_{i,+}-\displaystyle \frac{{\sigma }_{+}}{2}\left[({q}_{i,+}-{q}_{i,-})-\left(1-\displaystyle \frac{2}{3}{\sigma }_{+}\right){q}_{6,i}\right],\\ {u}_{i+1/2}^{\mathrm{ADV}}\gt 0,\\ {q}_{i},\\ \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{({q}_{R}^{1{\rm{D}}})}_{i-1/2}\\ &&=\left\{\begin{array}{l}{q}_{i,-}+\displaystyle \frac{{\sigma }_{-}}{2}\left[({q}_{i,+}-{q}_{i,-})+\left(1-\displaystyle \frac{2}{3}{\sigma }_{-}\right){q}_{6,i}\right],\\ {u}_{i-1/2}^{\mathrm{ADV}}\lt 0,\\ {q}_{i},\\ \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{q}_{i+1/2}^{1{\rm{D}}}={\mathcal{U}}\left({({q}_{L}^{1{\rm{D}}})}_{i+1/2},{({q}_{R}^{1{\rm{D}}})}_{i+1/2},{u}_{i+1/2}^{\mathrm{ADV}}\right).\end{eqnarray} \tag{ 21 }$$

We now have ${q}^{1{\rm{D}}}$ at each face. We now account for the transverse derivative terms to compute updated face- and time-centered states ${({q}_{L})}_{i+1/2,j}$ and ${({q}_{R})}_{i+1/2,j}$ using, in two dimensions,

$$\begin{eqnarray}{({q}_{L})}_{i+1/2,j} & = & {({q}_{L}^{1{\rm{D}}})}_{i+1/2,j}-\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{v}_{i,j+1/2}^{\mathrm{ADV}}+{v}_{i,j-1/2}^{\mathrm{ADV}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{q}_{i,j+1/2}^{1{\rm{D}}}-{q}_{i,j-1/2}^{1{\rm{D}}}}{{\rm{\Delta }}y}\right)+\displaystyle \frac{{\rm{\Delta }}t}{2}{({{\mathcal{S}}}_{q})}_{{\boldsymbol{i}}},\end{eqnarray} \tag{ 22 }$$

followed by an upwinding step to obtain the final face- and time-centered state,

$$\begin{eqnarray}&&{q}_{i+1/2,j}^{n+1/2}={\mathcal{U}}\left({({q}_{L})}_{i+1/2,j},{({q}_{R})}_{i+1/2,j},{u}_{i+1/2,j}^{\mathrm{ADV}}\right).\end{eqnarray} \tag{ 23 }$$

In three dimensions we use

$$\begin{eqnarray}{({q}_{L})}_{i+1/2,j,k} & = & {({q}_{L}^{1{\rm{D}}})}_{i+1/2,j,k}-\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{v}_{i,j+1/2,k}^{\mathrm{ADV}}+{v}_{i,j-1/2,k}^{\mathrm{ADV}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{q}_{i,j+1/2,k}^{y| z}-{q}_{i,j-1/2,k}^{y| z}}{{\rm{\Delta }}y}\right)\\ & & -\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{w}_{i,j,k+1/2}^{\mathrm{ADV}}+{v}_{i,j,k-1/2}^{\mathrm{ADV}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{q}_{i,j,k+1/2}^{z| y}-{q}_{i,j,k-1/2}^{z| y}}{{\rm{\Delta }}z}\right)+\displaystyle \frac{{\rm{\Delta }}t}{2}{({{\mathcal{S}}}_{q})}_{{\boldsymbol{i}}}\end{eqnarray} \tag{ 24 }$$

where to account for transverse corner coupling, we compute the intermediate states as in Colella (1990) and Saltzman (1994). For example,

$$\begin{eqnarray}{({q}_{L}^{y| z})}_{i,j+1/2,k} & = & {({q}_{L})}_{i,j+1/2,k}-\displaystyle \frac{{\rm{\Delta }}t}{3}\left(\displaystyle \frac{{w}_{i,j,k+1/2}^{\mathrm{ADV}}+{w}_{i,j,k-1/2}^{\mathrm{ADV}}}{2}\right)\\ & & \times \left(\displaystyle \frac{{q}_{i,j,k+1/2}^{1{\rm{D}}}-{q}_{i,j,k-1/2}^{1{\rm{D}}}}{{\rm{\Delta }}z}\right),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{q}_{i,j+1/2,k}^{y| z}={\mathcal{U}}\left({({q}_{L}^{y| z})}_{i,j+1/2,k},{({q}_{R}^{y| z})}_{i,j+1/2,k},{v}_{i,j+1/2,k}^{\mathrm{ADV}}\right).\end{eqnarray} \tag{ 26 }$$

A.3. Case 3

Here we compute face and time-centered estimates of a q that obeys a conservative equation. The details are the same as Case 2 up through Equation (22). We also note that in Step 2A, we use ${{\boldsymbol{U}}}^{\mathrm{ADV},\mathrm{pred}}$ instead of ${\tilde{{\boldsymbol{U}}}}^{\mathrm{ADV}}$. The difference between Case 2 is in the form of the corner coupling and transverse derivatives:

$$\begin{eqnarray}{({q}_{L})}_{i+1/2,j} & = & {({q}_{L})}_{i+1/2,j}^{1{\rm{D}}}-\displaystyle \frac{{\rm{\Delta }}t}{2}\\ & & \times \left(\displaystyle \frac{{({q}^{1{\rm{D}}}{v}^{\mathrm{ADV}})}_{i,j+1/2}-{({q}^{1{\rm{D}}}{v}^{\mathrm{ADV}})}_{i,j-1/2}}{{\rm{\Delta }}y}\right)\\ & & -\displaystyle \frac{{\rm{\Delta }}t}{2}{q}_{{\boldsymbol{i}}}\left(\displaystyle \frac{{u}_{i+1/2,j}^{\mathrm{ADV}}-{u}_{i-1/2,j}^{\mathrm{ADV}}}{{\rm{\Delta }}x}\right)+\displaystyle \frac{{\rm{\Delta }}t}{2}{({{\mathcal{S}}}_{q})}_{{\boldsymbol{i}}}.\end{eqnarray} \tag{ 27 }$$

We apply the same upwinding procedure in Equation (23) to obtain the final face- and time-centered state, ${q}_{i+1/2,j}^{n+1/2}$. In three dimensions we use

$$\begin{eqnarray}{({q}_{L})}_{i+1/2,j,k} & = & {({q}_{L}^{1{\rm{D}}})}_{i+1/2,j,k}-\displaystyle \frac{{\rm{\Delta }}t}{2}\\ & & \times \left(\displaystyle \frac{{({q}^{y| z}{v}^{\mathrm{ADV}})}_{i,j+1/2,k}-{({q}^{y| z}{v}^{\mathrm{ADV}})}_{i,j-1/2,k}}{{\rm{\Delta }}y}\right)\\ & & -\displaystyle \frac{{\rm{\Delta }}t}{2}\left(\displaystyle \frac{{({q}^{z| y}{w}^{\mathrm{ADV}})}_{i,j,k+1/2}-{({q}^{z| y}{w}^{\mathrm{ADV}})}_{i,j,k-1/2}}{{\rm{\Delta }}z}\right)\\ & & -\displaystyle \frac{{\rm{\Delta }}t}{2}{q}_{{\boldsymbol{i}}}\left(\displaystyle \frac{{u}_{i+1/2,j,k}^{\mathrm{ADV}}-{u}_{i-1/2,j,k}^{\mathrm{ADV}}}{{\rm{\Delta }}x}\right)+\displaystyle \frac{{\rm{\Delta }}t}{2}{({{\mathcal{S}}}_{q})}_{{\boldsymbol{i}}},\end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray}&&{({q}_{L}^{y| z})}_{i,j+1/2,k}={({q}_{L})}_{i,j+1/2,k}-\displaystyle \frac{{\rm{\Delta }}t}{3}\\ &&\times \left(\displaystyle \frac{{({q}^{1{\rm{D}}}{w}^{\mathrm{ADV}})}_{i,j,k+1/2}-{({q}^{1{\rm{D}}}{w}^{\mathrm{ADV}})}_{i,j,k-1/2}}{{\rm{\Delta }}z}\right),\end{eqnarray} \tag{ 29 }$$

along with the upwinding procedure in Equation (26).

Based on feedback from collaborators and the anonymous referee, we have performed several two-dimensional XRB tests to see how either the initial conditions or domain size can alter the effects of the convection.

The first test we performed was to alter the strength of the initial velocity perturbations that act as a seed to the convection. Even though we expect the convection to "forget" how it was initiated during the ∼200 convective turnover times we simulate, this is an important check. The default velocity perturbations in the simulations of the main paper, as well as those of Malone et al. (2014), were 1 × 10⁵ cm s⁻¹. Figure 8 shows both the peak temperature and peak Mach number as a function of time for five simulations, each with different initial velocity perturbations as shown in the labels; here we only plot every 500 steps to cut down on the image size. All simulations track each other very nicely, and thus the development of the convection does not depend strongly on the strength of the initial perturbations. Figure 9 shows colormaps of the magnitude of velocity for each of the five runs. They likewise compare well with one another.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second test was designed to see if the domain size choice we had made for our default runs, including those of the three-dimensional simulation, were affecting the evolution of the convective flow pattern. To this end, we performed two additional simulations where we extended the default two-dimensional domain (3072 cm × 4608 cm) in both the vertical direction (to 3072 cm × 5760) and again in the horizontal direction (to 6144 cm × 5760 cm). The results are shown in Figure 10 where we plot the magnitude of velocity for each of the three runs at t = 0.01 s; the default domain is on the left and has been shifted up by the 1152 cm difference, the middle frame shows the tall domain, and the right shows the tall and wide domain. These plots show that the shape and strength of the flow in the convective region are only weakly affected by the domain size. We note that had we chosen a domain width less than the scale-height of the convection, this would have constrained the eddy size quite strongly.

Zoom In Zoom Out Reset image size