Simulating Lateral H/He Flame Propagation in Type I X-ray Bursts

X-ray bursts are the thermonuclear runaway of a mixed H/He layer on the surface of a neutron star. Observations suggest that the burning begins locally and spreads across the surface of the star as a flame. Recent multidimensional work has looked in detail at pure He flames spreading across a neutron star. Here we report on progress in multidimensional modeling of mixed H/He flames and discuss the challenges.


Introduction
Type I X-ray bursts (XRBs) occur in accreting binary systems between a neutron star and a low-mass companion.As a layer of fuel builds up on the surface of the neutron star, gravitational compression heats it up until the conditions at the base lead to a thermonuclear runaway.A challenging aspect of modeling these systems is capturing how the flame spreads from the initial ignition point across the surface of the star.This is important for understanding oscillations observed during the initial rise of some bursts, which are believed to arise from the bright hotspot coming in and out of view as the neutron star rotates [1].Unfortunately, one-dimensional simulations cannot capture these dynamics, as they are inherently multi-dimensional.
Most burst sources accrete a mix of 1 H and 4 He from a main sequence companion, while some ultra-compact systems accrete nearly pure 4 He [2,3].Depending on the accretion rate, the hydrogen can burn stably into helium which will then ignite a pure helium burst, or the helium could ignite before all the hydrogen is consumed, leading to a mixed H/He burst.Previous 2D and 3D simulations by our group and others have been limited to pure helium flames [4][5][6][7][8][9][10] as the mixed H/He reaction networks are usually larger and more difficult to integrate, and the timescales are slower due to the weak rate waiting points [11].We recently showed [10] that the nucleosynthesis and evolution of a spreading flame in 3D matches 2D simulations very closely, so we will focus on 2D simulations here.

Simulations
We use the Castro [12] compressible hydrodynamics code to perform our simulations, following the simulation methodology developed in [7].Castro uses the AMReX framework [13] for adaptive mesh refinement and parallelization across CPUs and GPUs.We use the stellar equation of state based on [14] and conductivities from [15] for the thermal diffusion operator.Castro can use an arbitrary nuclear reaction network, and we discuss some of the options below.All simulations use a two-dimensional axisymmetric simulation domain with a physical size of 2.4576 × 10 5 cm × 3.072 × 10 4 cm and a base grid of 1024 × 128 zones.The grid is subdivided using two levels of refinement, refining first by a factor of 4, then by a factor of 2. This gives a resolution of 30 cm at the finest level.We note that this is a little coarser than the He simulations we have done previously.However, since the addition of H to the atmosphere increases the scale height and H burning is less strongly peaked than He burning, we believe we can relax the resolution requirements.
Our initial runs used rprox, a 10 isotope nuclear reaction network introduced in [16] that models hot CNO burning, triple-α, and rp-process breakout.We vary the proportions of 1 H and 4 He in the fuel layer, and also include 1% of 14 O and 15 O to help start CNO burning.The oxygen abundance is split between the two isotopes in proportion to their respective β + -decay lifetimes.We assume that the composition is completely mixed by a period of convection preceding the ignition of the flame, as explored in [17].
The initial model is set up as described in section 4 of [7].In summary, it consists of a hot region behind the flame at T = 1.4 × 10 9 K and a cool region ahead of it at T = 3 × 10 8 K, with a density of 2 × 10 6 g cm −3 at the base of the atmosphere (z = 2000 cm).These two regions are constructed individually in 1D to be in hydrostatic equilibrium vertically, then smoothly blended at the interface, giving us an equilibrium starting condition that will lead to the ignition of a burning front.Figure 1a shows the temperature profiles for the two regions for one of the runs.The other profiles with different compositions are similar, but those with more helium are squashed due to the greater mass of the atmosphere above the surface.The hot region extends to r pert = 6.144 × 10 4 cm (25% of the domain), and the transition interface between the two regions has a width of δ blend = 2.048 × 10 3 cm.Paired with the axisymmetric geometry, this effectively models a hotspot at the pole of the star spreading outwards towards the equator.The rotation rate of the neutron star is set to an artificially high 1000 Hz, which helps to confine the flame laterally and allows us to use a smaller simulation domain.All of our simulations were run on the OLCF Summit machine, with the entire calculation offloaded to the NVIDIA V100 GPUs there.
Figure 2 shows snapshots of the mass fraction of ash (everything heavier than 16 O) after 80 ms for each of the initial compositions we used.We have found that this diagnostic is much easier to compare between the different compositions than the mean molecular weight, Ā, used in our previous work.The runs in the first two panels never produced well-defined flames.The  third panel has more burning than the previous two, but the flame front never moves beyond 10 5 cm.The most interesting is the fourth panel, with 9% H and 90% He, which shows a flame front with burning up to 56 Ni.The final panel models a pure He burst using the same network for comparison.Many of our runs exhibited an unexpected high-temperature region at the base of the atmosphere across the entire simulation domain, starting from the outer radial edge of the domain and spreading toward the pole.This often ignited and consumed all the fuel before the flame could set up and start propagating.The magenta layer on the right of the first panel in figure 2, for example, is a result of this excess burning.We attributed this to our selection of the base temperature in the cool region, and used a lower temperature in later runs to try to avoid this.This is one of the main differences between a pure He burst-a pure He network will not burn much at those temperatures.However, we found that the simulation behavior was very sensitive to any changes in the initial conditions, and it was difficult to reproduce a steady burning front.To explore this behavior in more detail, we experimented with 1D flame simulations using different reaction networks, and found that aprox19 produced a much more robust propagating flame front than rprox.Due to the different composition in that network, we replaced the oxygen isotopes with 12 C.
For our next set of runs with aprox19, we increased the resolution to 20 cm (using a base grid of 1536 × 192) and used a new lower boundary condition with a well-balanced pressure reconstruction, as done in [10].We reduced the cool region temperature to 2 × 10 8 K and the material below the atmosphere (i.e. the star) to 1 × 10 8 K, as mentioned above.The new initial conditions are shown in figure 1b.Interestingly, increasing the resolution this way actually scaled slightly better than linear in the number of compute nodes, as the new grid size distributed more evenly onto Summit's 6 GPUs per node.
Figure 3 shows the same diagnostics as figure 2 for the aprox19 runs.We can see that the run with 25% He again doesn't produce a flame, although it does have more ash present in    the atmosphere than the corresponding rprox run.This is due to vigorous carbon burning in the first few milliseconds of the simulation, which consumes the majority of the initial 12 C in the lower atmosphere.In other two runs, the right side of the atmosphere heats up and starts burning before the flame front.An example of how this progresses for the run with 10% H and 88% He is shown in figure 4. Finally, figure 5 shows the evolution of each species in the network over time, for the same run.We can see that most of the initial 12 C (solid red line) is consumed in under 1 millisecond.An interesting detail in this plot is the step-like features at early times, most visible in the profiles for 28 Si and 32 S.These are fairly evenly spaced in time with a separation of ∼ 5 × 10 −4 s, which matches the width of the domain divided by the average sound speed near the base of the

Mass fraction
Figure 5: Total mass fractions for each species in the aprox19 run with 10% H, excluding those which never increased above 10 −10 .atmosphere: L x /c s ≈ 2.4576 × 10 5 cm/5.2 × 10 8 cm s −1 ≈ 4.7 × 10 −4 s.This suggests they may be related to acoustic oscillations travelling across the simulation domain and reflecting off of the boundaries.

Summary
We ran several two-dimensional simulations of mixed H/He flames in XRBs, varying the proportions of H and He in the atmosphere.We found that the flames are very sensitive to the composition and initial temperatures, which makes it difficult to set up a flame in the first place.Using the rprox network, we were only able to produce a stable, propagating flame with an initial composition of 9% H and 90% He.The aprox19 network gave much better results in 1D simulations, but we have not yet been able to reproduce a clean H/He flame in 2D.
An unexpected hot region would often develop at the surface of the star, away from the flame, which would start burning and disrupt the flame.The source of this heating is unclear, and needs to be investigated further.We are currently testing out different reaction networks with pynucastro [18,19], which allows us to construct and inspect arbitrary reaction networks, and generate C++ code for use in our Castro simulations.It also allows us to disable specific rates at runtime, which can help us understand which rates are relevant for H/He burning under these conditions.Finally, we may need to consider different initial models, where the pre-runaway convective mixing is confined only to the He layer at high column depth, as explored in [20].

Figure 1 :
Figure 1: Initial vertical temperature profiles for the cool (solid) and hot (dashed) regions.

Figure 2 :
Figure 2: Total mass fraction of all species heavier than 16 O for rprox runs at t = 80 ms.

Figure 3 :
Figure 3: Total ash mass fraction at t = 80 ms for aprox19 runs.

Figure 4 :
Figure 4: Ash mass fraction evolution for the aprox19 run with 10% H and 88% He.