Comparing Early Evolution of Flames in X-Ray Bursts in Two and Three Dimensions

X-ray bursts (XRBs) result from thermonuclear burning of an accreted H/He or He layer on a neutron star (Galloway & Keek 2021). Observations of brightness oscillations in the rise of the light curve suggest that the burning begins localized and then spreads across the neutron star (Bhattacharyya & Strohmayer 2007). This spreading is inherently a multidimensional phenomenon, requiring hydrodynamic simulations that resolve the reactive zone, realistic reaction networks, and domain sizes that capture the scales over which rotation is important in order to understand the flame propagation and nucleosynthesis. Detailed nuclear physics is especially important for understanding the rp-process in mixed H/He bursts, as discussed in Schatz et al. (2001).

Here we show a first attempt at modeling the early evolution of a spreading hotspot in three dimensions. This builds on our earlier two-dimensional work (Eiden et al. 2020; Harpole et al. 2021) that developed our simulation framework and explored how the acceleration of the burning front depends on the initial model structure. As with those simulations, we use the freely available Castro simulation code (Almgren et al. 2010, 2020), with all of the code needed to run these simulations in our public GitHub repositories.

Our simulations complement the two- and three-dimensional models of flames of Cavecchi et al. (2013, 2015, 2016) and Cavecchi & Spitkovsky (2019) by focusing on resolving the reaction zone of the flame and modeling the vertical structure hydrodynamically instead of hydrostatically. We also focus on modeling flames starting from a small hotspot. This however means that we are confined to smaller domains, with similar computational resources, so one of the goals of this study is to understand how far we can push resolved XRB flame simulations. The work by Goodwin et al. (2021) also explores multidimensional evolution, looking at the thermal transport and how it affects the location where the hotspot ignites. Together, all three approaches help build a comprehensive multidimensional picture of X-ray bursts.

The goal of this paper is to understand the challenges of resolved three-dimensional hydrodynamical simulation and explore what is possible with today's resources. This will inform our follow-up simulations.

In Eiden et al. (2020), we developed the simulation methodology in Castro for modeling laterally spreading flames through the accreted layer on the surface of a neutron star. For these models, Castro evolves the compressible Euler equations with a Coriolis force arising from rotation, constant gravity (the plane-parallel approximation), reactive sources, and thermal diffusion (integrated explicitly):

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {X}_{k})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}}{X}_{k})=\rho {\dot{\omega }}_{k},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{U}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}}{\boldsymbol{U}})+{\rm{\nabla }}p=\rho {\boldsymbol{g}}-2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{U}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho E)}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}}E+p{\boldsymbol{U}})=\rho {\boldsymbol{U}}\cdot {\boldsymbol{g}}+{\rm{\nabla }}\cdot {k}_{\mathrm{th}}{\rm{\nabla }}T+\rho \dot{\epsilon }.\end{eqnarray} \tag{ 4 }$$

Here, ρ is the mass density, X_k are the mass fractions of the nuclear species with associated creation rates ${\dot{\omega }}_{k}$, U is the velocity, E is the specific total energy, p is the pressure, g is the gravitational acceleration, Ω is the rotation vector, k_th is the thermal conductivity, T is the temperature, and $\dot{\epsilon }$ is the energy release from the reaction network.

We use a general stellar equation of state with radiation, ions as a perfect gas, and degenerate/relativistic electrons described in Timmes & Swesty (2000), which closes the system:

$$\begin{eqnarray}&&E=e+\displaystyle \frac{1}{2}| {\boldsymbol{U}}{| }^{2},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&p=p(\rho ,e,{X}_{k}),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&T=T(\rho ,e,{X}_{k}).\end{eqnarray} \tag{ 7 }$$

Building on these first simulations, in Harpole et al. (2021), we explored several different rotation rates and initial model thermal structures to understand how the flame speed depends on the thermal structure of the envelope. Both of these papers used a two-dimensional axisymmetric geometry. The three-dimensional simulations we present here build on these works. We use a neutron star crust temperature of 2 × 10⁸ K and rotation rate of 1000 Hz, as followed in Harpole et al. (2021). This model was chosen here since it gives a clean well-defined flame. The model is constructed following the same procedure described in Eiden et al. (2020) and Harpole et al. (2021). We generate a hot and a cool model, and blend them laterally to setup an initial perturbation that is in hydrostatic equilibrium. The structure of these models is shown in Figure 1.

Zoom In Zoom Out Reset image size

Moving to 3D is very expensive, so to further reduce the computational expense, we use an anisotropic resolution, with 32 cm resolution laterally and 16 cm resolution vertically. This is finer vertically than done in Harpole et al. (2021). Because the hotspot is placed in the center of the domain, the size of the domain also becomes a constraint, and as a result, the evolution time we can reach will become limited by the time it takes the burning to approach the edge of the domain. Consequently, we focus here on the early evolution.

Castro uses the the AMReX adaptive mesh refinement library (Zhang et al. 2019) to manage the discretization and domain decomposition. Our 3D simulation used a base grid of 768²×192 zones with two levels of refinement, the first jumping by a factor of 4 and the next by a factor of 2. We subcycle in time, so the coarser grids take a larger time step than the finer grids. We used static mesh refinement, to achieve better load balancing, fully refining the atmosphere below a height of 3600 cm. The domain has a size ${\left(1.96608\times {10}^{5}\,\mathrm{cm}\right)}^{2}$ × 2.4576 × 10⁴ cm. The simulation methodology and parameters are mostly identical to those used in Harpole et al. (2021). We use an unsplit piecewise parabolic method (Colella 1990; Miller & Colella 2002) for the advection and operator (Strang) splitting for the reactions, evolving an internal energy evolution equation during the burn, as described in Zingale et al. (2021). These simulations used the seven-isotope He-burning network described in Timmes et al. (2000). This was one of the networks used in Eiden et al. (2020). Thermal diffusion is treated explicitly, using the conductivities from Timmes (2000). The overall integration strategy is second-order accurate in space and time. The main difference with the previous simulations is that we switched from a hydrostatic to reflecting boundary at the bottom of the domain and used a simple well-balanced scheme (similar to Käppeli & Mishra (2016) but adapted to deal with the characteristic tracing in the piecewise parabolic method as described in Zingale et al. 2002). Tests demonstrated that this boundary does a better job than our previous approach at supporting the atmosphere as the flame propagates for long times.

The simulations were run on the Oak Ridge Leadership Computing Facility Summit supercomputer, on 342–1366 nodes, with six NVIDIA V100 GPUs per node. The entire computation was offloaded to GPUs using the AMReX C++ abstraction layer for parallel for loops, as described in Katz et al. (2020). GPU offloading was critical to being able to perform these simulations, as we run more than an order of magnitude faster using all the GPUs on a node as compared to all of the CPU cores on the node. Overall, about 250,000 node hours were used for the calculation. The full state of the calculation was output only every 0.005 s, as each output file is 2.4 TB in size (for single-precision data). We output a few fields more frequently for visualization, and global diagnostics were output every time step. We ran a 2D simulation with the same resolution and refinement strategy to compare with here.

Figure 2 shows the three-dimensional simulation via top-down volume rendering of the mean molecular weight,

$$\begin{eqnarray}&&\displaystyle \frac{1}{\bar{A}}=\displaystyle \sum _{k}\displaystyle \frac{{X}_{k}}{{A}_{k}},\end{eqnarray} \tag{ 8 }$$

where X_k is the mass fraction of species k and A_k is the atomic weight (in mass units). The structure is shown at 10, 20, and 40 ms of simulation time. While initially quite symmetric, the axial symmetry is broken, likely due to roundoff error, and the flame structure takes on a more tenuous form. The panels show the flame spreading considerably through the domain as time evolves, and the simulation is halted once the flame nears the boundaries.

Zoom In Zoom Out Reset image size

Figure 3 shows the corresponding two-dimensional axisymmetric simulation. In the axisymmetric geometry, each zone is essentially a torus, rotated about the vertical axis. The initial model and simulation parameters are identical to the three-dimensional case. We see that the overall structure and evolution looks identical to that reported in Harpole et al. 2021. We also see that there is a well-defined flame in the simulation plane that spreads outward with time. A vertical slice through the three-dimensional simulation (in the x–z plane) is shown in Figure 4 at the final time, 40 ms. Aside from the lack of reflection symmetry from the absence of axisymmetric geometry, it looks very similar to the two-dimensional case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because the 3D flame does not stay perfectly circular, it is difficult to define the flame speed using the same method employed in Eiden et al. (2020). Instead, we will look at the mass of the ash material to assess how quickly the burning is taking place. Figure 5 shows the mass of the different species as a function of time, for the 2D and 3D calculation. The axisymmetric 2D calculation has slightly less volume than the 3D domain, as rotating about the symmetry axis produces a cylindrical domain inscribed in the 3D domain. To compensate for this, we scale the masses by the total mass on the grid. This plot shows that the masses of the heavy species, especially ¹²C, grow quickly with time. The differences between the 2D and 3D simulations appear small—the burning is slightly faster in 2D, but the trends track very well. This slight difference is likely because the assumption of axisymmetry there does not allow for the complex structure we see in the evolution of composition in the 3D flame. This strong agreement suggests that using 2D axisymmetric calculations is a good model for the early flame evolution. As they are so much less computationally expensive, this will allow us to explore much more complex reaction networks and understand the nucleosynthesis in more detail.

Zoom In Zoom Out Reset image size

We have extended our simulation framework to three dimensions and explored the differences between a two-dimensional axisymmetric simulation and a fully three-dimensional hydrodynamical simulation. Overall, the early evolution of the flame spreading from a hotspot behaves similarly in two and three dimensions. As this three-dimensional simulation pushed the limits of available computational resources, this result suggests that we can continue to use two-dimensional axisymmetric simulations to explore the initial flame structure and spend our computational resources on larger networks and bigger domains.

These results also suggest that we can do initial explorations of full star bursts in two dimensions, perhaps using axisymmetric spherical coordinates (modeling an annular region in neutron star radius and colatitude). In this geometry, point ignition requires starting at the poles, but this would allow us to see the effect of varying gravity with latitude (due to centrifugal forces), as well as begin to consider light curve modeling.

There is still a role for three-dimensional simulations—the convective burning in the accreted layer should create turbulence that the flame will encounter as it propagates. Although we have shown that for the flame in a Type Ia supernova, the turbulence from the era of convection is not strong enough to disrupt the flame (Nonaka et al. 2011), the flame in an XRB is considerably slower and thicker, so it remains to be seen what effect it might have. This can be assessed in a different geometry, like a long channel, which would allow us to trade domain size for resolution to increase the numerical Reynolds number. We also will consider higher-order simulation methodologies, which could allow us to capture these effects at lower resolution. We have already demonstrated a fully fourth-order accurate (in space and time) algorithm for coupling hydro, diffusion, and reactions in Castro (Zingale et al. 2019) that could be used here. The main outstanding issue with applying that work to the present problem is the multilevel time integration.

We are currently exploring mixed H/He bursts as well as the sensitivity of the flame properties to the size of the reaction network used. Both of these initial studies use our two-dimensional axisymmetric geometry.

Castro is open-source and freely available at http://github.com/AMReX-Astro/Castro. The problem setup used here is available in the git repo as flame_wave. The metadata describing the build environment and the global diagnostic output files are available on Zenodo at Zingale (2023). We thank Alice Harpole for contributions to the AMReX Astrophysics suite.

The work at Stony Brook was supported by DOE/Office of Nuclear Physics grant DE-FG02-87ER40317. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U. S. Department of Energy under Contract No. DE-AC02-05CH11231. This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. 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, awarded through the DOE INCITE program. We thank NVIDIA Corporation for the donation of a Titan X and Titan V GPU through their academic grant program. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

Facilities: NERSC - , OLCF - .

Software: AMReX (Zhang et al. 2019), Castro (Almgren et al. 2010, 2020), GCC (https://gcc.gnu.org/), helmeos (Timmes & Swesty 2000), linux (https://www.kernel.org/), matplotlib (Hunter 2007, http://matplotlib.org/), NumPy (Oliphant 2007; van der Walt et al. 2011), python (https://www.python.org/), valgrind (Nethercote & Seward 2007), VODE (Brown et al. 1989), yt (Turk et al. 2010).