Three-dimensional Instability of Flame Fronts in Type I X-Ray Bursts

We present the first realistic 3D simulations of flame front instabilities during type I X-ray bursts. The unperturbed front is characterized by the balance between the pressure gradient and the Coriolis force of a spinning neutron star (ν = 450 Hz in our case). This balance leads to a fast horizontal velocity field parallel to the flame front. This flow is strongly sheared in the vertical direction. When we perturb the front an instability quickly corrugates the front. We identify this instability as the baroclinic instability. Most importantly, the flame is not disrupted by the instability and there are two major consequences: the overall flame propagation speed is ∼10 times faster than in the unperturbed case and distinct flame vortices appear. The speedup is due to the corrugation of the front and the dynamics of the vortices. These vortices may also be linked to the oscillations observed in the light curves of the bursts.


Introduction
The type I X-ray bursts are thermonuclear explosions that burn the material that neutron stars accrete from their companion in low-mass X-ray binaries. The matter is stripped from the companion via Roche-lobe overflow and forms an accretion disk around the neutron star. It eventually reaches the star and spreads across its surface. With time, the fresh material is compressed by newly incoming material, and it sinks to deeper, higher-pressure layers. In this process the fuel, which is rich in hydrogen and helium, undergoes nuclear burning. Depending on the composition and the mass accretion rate, the heating from the nuclear burning may become uncompensated by the cooling processes, and the temperature rapidly increases leading to a thermonuclear runaway: the ignition of type I bursts (Fujimoto et al. 1981;Bildsten 1998;Strohmayer & Bildsten 2006;Galloway & Keek 2017).
One important aspect of the bursts is flame propagation: the fuel is assumed to cover the whole star, but synchronous ignition is unlikely (Shara 1982), and the flame needs to propagate from the first ignition site to engulf the star. The most likely ignition site is supposed to be the equator (Spitkovsky et al. 2002). The flame propagation needs to be fast enough to be compatible with the ∼1 s rise time of the burst light curves and yet cannot be a detonation, since the necessary conditions for detonations are not reachable (Fryxell & Woosley 1982;Malone et al. 2011). This poses stringent constraints on the flame-spreading mechanism. Furthermore, the flame propagation (and the subsequent cooling of the surface) should explain the temperature distribution across the star surface and thus the observed features of the observed light curves.
The details of burst flames have been explored in a series of papers: Malone et al. (2011) and Zingale et al. (2015) explored the onset of nuclear-burning ignition and the convection driven by nuclear heating. Spitkovsky et al. (2002) presented a shallow water calculation suggesting the Coriolis force should be important for ignition and front propagation. Cavecchi et al. (2013Cavecchi et al. ( , 2015Cavecchi et al. ( , 2016 explored the vertical structure of the front and its spreading across the surface. The picture emerging from these works is the following: the hot burning fluid expands vertically, creating a horizontal pressure gradient that pushes the hot fluid toward the cold one. The Coriolis force balances the pressure gradient (geostrophic balance), generating a horizontal flow parallel to the flame front. This flow follows the thermal wind profile (e.g., Pedlosky 1987) showing a strong vertical shear. The hot fluid is thus being confined within a few Rossby radii. The Rossby radius is the radius over which the Coriolis force is capable of deflecting the horizontal motion, where Ω is the angular velocity of the star, g the gravitational acceleration, and H the scale height of the fluid.
The flame propagation is driven mostly by conduction across the hot-cold fluid interface at the front. The latter is almost parallel to the horizontal, since the Rossby radius is of ordeŕ 5 10 4 cm, while the layer height is of order 10 2 cm, and the large conducting surface guarantees a high speed in good agreement with the observations. The tension provided by a sufficiently strong magnetic field, ∼10 8 G, will compete with the Coriolis force, reducing the thermal wind and increasing the confinement length, thus making the front even more horizontal, increasing the conducting surface and speeding up the flame (Cavecchi et al. 2016).
One important question left unanswered by these studies is the stability of the flame front. The thermal wind can be highly unstable to the baroclinic instability (see, e.g., Pedlosky 1987;Vallis 2006). In a nutshell, baroclinic instability is a kind of quasi-horizontal convection that taps the available gravitational potential energy when the surfaces of constant pressure and the surfaces of constant density are not coincident. Baroclinic instability has been previously discussed within the burst framework by, e.g., Fujimoto (1988) and Cumming & Bildsten (2000). These authors addressed angular momentum and mixing in differentially rotating atmospheres, in particular following the expansion of the burning fluid. Cavecchi et al. (2013) were the first to identify this instability along the fronts of burst flames. However, their simulations were 2.5D, assuming symmetry along the y direction (parallel to the flame front), and the instability could only develop on the vertical plane, bringing cold unburnt fluid from the top of the flame front down to the highest burning region, thus helping propagation.
In this paper we study the full development of front instabilities in 3D. The description of the numerical setup is given in section Section 2. In Section 3 we describe the development of the instability, and we discuss the results in Section 4.

Numerical Setup
In order to study the instability we first ignite and let the flame propagate in a 2.5D simulation as in Cavecchi et al. (2013Cavecchi et al. ( , 2015Cavecchi et al. ( , 2016, where the y direction has to be thought of as the equator to the pole direction and the z direction is the radial direction. We then extend the simulation in the third dimension (x, west to east), perturb the velocity field, and let the instability set in.
The code we used is described in Braithwaite & Cavecchi (2012) and Cavecchi et al. (2013Cavecchi et al. ( , 2015Cavecchi et al. ( , 2016, and specific details can be found there. For this paper we used only the hydrodynamical modules, switching the magnetic field off. The code is a finite-difference, explicit time-stepping code, using a three-step Runge-Kutta scheme designed to save memory impact (Williamson 1980). Derivatives, interpolations, and integrals use six-point stencils on a staggered grid that reach spectral-like accuracy (Lele 1992). Velocities in the three directions are evaluated on the cell faces; state variables like density, temperature, and mass fraction of the different species are evaluated at the cell centers.
The equations solved assume hydrostatic equilibrium in the vertical direction (parallel to the radial direction, z) and use pressure as the vertical coordinate. We simulate a box with reflecting boundary conditions in the y direction and periodic boundaries in the x direction. Vertically, we use symmetrical conditions for horizontal velocities and the state variables like temperature. We further impose cooling losses from the top, approximating radiative conditions at the top of the atmosphere (see Equation ( Different from our previous simulations we add a layer of hydrogen at the top and thicken the layer of carbon at the bottom. Our reaction network only burns helium into carbon through triple alpha reactions, so that hydrogen is not involved in the flame propagation. However, it provides a layer of light fluid on top of the burning fluid, emulating the role of newly accreted matter and providing space for the instability to develop into. A similar role is played by the carbon layer, which also emulates the ashes of previous bursts. The initial profiles of the mass fraction X of hydrogen and Y of helium are initialized as a function of pressure P (see also Figure 1, left): To ensure numerical stability we have added two regions of higher viscosity at the equatorial and polar walls and at the bottom of the carbon layer. The flame front perturbation takes place far enough in the simulation domain that these regions do not affect the dynamics of the instability. Finally, in order to be able to see the effects of cooling more quickly we also enhance the cooling from the top by a factor of 5. When the flame has reached steady propagation, we extend the domain in the x direction and perturb the U x velocity field in the following way. We identify a location y z , 0 0 ( ) near the highest burning and extract the value of U x 0 there. We then subtract the value U f x y z , , x 0 ( ) from the velocity field. f (x, y, z) is given by where j is a random phase between 0 and 2π. We use = y 7.2 0 km, d = y 187.5 m, z 0 =1.51 m, δ z=10 cm, and Δz=10 cm. This perturbation resemble a smooth bell-shape on the y-z plane, whose center varies sinusoidally with x. The choice of 64 for the last sin wavenumber is arbitrary, but the wavenumbers after ∼10 do not really contribute to the final result.

The Flame Front Instability
The right panel of Figure 1 shows the initial background conditions for the perturbed run. In particular, we plot the azimuthal velocity profiles with filled contours and the density stratification with dotted lines. One can identify the thermal wind profile and its vertical shear. Also note that where the density contours (isopycnals) slope upward is where the misalignment with isobars peaks (and so does baroclinicity ∇P×∇ρ): this is where we seed our perturbation (at the location of the dashed vertical line: = y 7.2 km). The angle the isopycnals form with the horizontal is called the wedge of instability (Pedlosky 1987;Vallis 2006), because exchanging fluid elements within this angle brings lighter elements among denser fluid and vice versa, and buoyancy then sustains the perturbation leading to the instability. The overall center of mass of the fluid is lowered and gravitational potential energy is converted into kinetic energy. This angle is very shallow: note that the vertical length scale of Figure 1 is just a few meters while horizontally it is kilometers, so that the motion following the instability is almost horizontal. Another aspect to consider is the effect of the pure horizontal velocity shear, which can also contribute to the growth of instabilities, even though usually on smaller scales (this is called the barotropic instability ;Pedlosky 1987).
This configuration and the development of the baroclinic instability is a very active field of research in geophysical sciences. A direct comparison with geophysical results is, however, difficult, because there are no general models for the development of the instability (especially for the nonlinear regime). Also, the systems studied in the literature are either simplified models with a few layers or have important differences from our simulation. For example, we have an important internal heat source (the burning) that reinforces the conditions for instability and is not always present in the geophysical models. When heat is included in those models, as in the case with heat released by condensation, it usually takes places in the top layers and not at the bottom like in our case.
One important point in our simulations regards the Coriolis force. We model an f plane (Pedlosky 1987), which is a system with constant Coriolis force, while in reality its intensity should change with latitude. This effect is captured for example by the so-called β-plane models, but in this exploratory work we neglect it. Perhaps the most easily comparable common geophysical model is the Eady model (Eady 1949), which consists of a uniformly stratified ideal gas, confined by two rigid lids on an f plane, with uniform vertical shear. In this model, the growth rate of the most unstable mode is given by ) . U is the velocity of the top layer and L the width of the channel where the instability takes place. While simplified, this model captures most of the features of the instability (Vallis 2006) and we will briefly compare our results to it. Finally, in the Appendix we show snapshots in 3D in order to provide a global view of the evolution of the perturbation that we describe below.

The Early Stages of the Instability
We have measured the growth rate of the instability during the linear phase. This is shown in Figure 2, where we overplot the prediction of (3) with parameters, taken form the unperturbed simulation, =Ŕ 5 10 Ro 4 cm, U=10 8 cm s −1 , L=3R Ro . This choice is approximated, and somewhat ad hoc, but the qualitative picture agrees. The fastest-growing mode has a wavelength of ∼2 km, which corresponds to approximately 4 Rossby radii, and the shortest wavelengths are stable. 3 This compares qualitatively well with what is predicted by the Eady model if one considers that the analogy is not complete since our conditions are not exactly the same (for instance, the top layer is not a rigid lid). Thus, we can consider this a robust result and expect that the fronts on the surface of neutron stars will corrugate when they exceed a few Rossby radii in length. Figures 3 and 4 show the instability at the end of the linear phase-beginning of the nonlinear one. In the left panel of Figure 3 we plot vertical, azimuthal slices, at the location of the dashed line of Figure 1, right, for the temperature and the meridional and vertical velocities. For each field we subtract the initial profile to highlight the effect of the instability. The right panel shows streamlines on horizontal slices at various heights. The red lines follow the unperturbed flow, while the black lines show the residuals. To guide the eye, the thickness of the lines is proportional to the local speed of the streamline, but for visibility the scale is different for each group of lines. The most apparent feature is the typical leaning of the velocity contours against the background flow. Also, note that the north-upward perturbation flows correspond to the higher temperatures and the south-downward ones to the colder fluid: there is a net flux of heat toward the pole. Analogously, the flow is carrying combustible fuel south and more ash-rich fuel north. This configuration is typical for the instability and is compatible with what we saw in the 2.5D simulations of Cavecchi et al. (2013Cavecchi et al. ( , 2015. Also in the 3D run we identify this trend as a contribution to the speedup of the flame.
It can be seen that at these stages the dominating mode has wavenumber 5, but asymmetries in the peaks of the contours show that the mode with wavenumber 1 is quickly catching up while the instability is saturating. This begins to be evident already in Figure 4, which is a snapshot taken only 0.01 s later, and becomes manifest in Figure 5, which shows the evolution in the fully nonlinear stage. The development is reminiscent of the mushroom shape of the Rayleigh-Taylor instability, but the motion is almost along a horizontal direction. The fluid slides along the isopycnal surfaces and breaks into cyclonic and anticyclonic systems: see the systems respectively to the right and to the left of the dotted line corresponding to an isobar in the top panel of Figure 5 and compare to the streamlines drawn in the bottom panel. 4

Later Stages and the Development of Giant Vortices
The net effect of the flows described above is to trigger circulation on top of the original azimuthal flow and to cause ignition at various locations northward. Eventually, various systems of locally higher-temperature, higher-pressure and lower-temperature, lower-pressure form at alternate y coordinates (corresponding to higher latitudes) creating a complex combination of anticyclones and cyclones (see Figures 6 and 7).
Bigger systems drift and interact with smaller vortices absorbing or disrupting them. One large anticyclone, approximately 4 km across (∼8R Ro ), is fully formed near the left atỹ 4 km, t=4.25 s (∼0.5 s after the perturbation was triggered). A second anticyclone, which originated from a small system detached roughly at t=4.07 s from the first one, is clearly recognizable at t=4.48 andỹ 11 km. This one in turn launched another small system at t=4.38 s that eventually becomes a third large eddy identifiable at t=4.69 s andỹ 16 km. The "launching" happens when the the major systems enter in contact with a smaller cyclonic system. Due to the mixing effect of this interaction, the flame reaches the north boundary very quickly. The fact that we observe three such systems is probably due to the fact that only three could fit within the domain. Different stars spinning at different frequencies will most probably show a different number of such anticyclones.
By late times (bottom rows of Figures 6 and 7) the three large anticyclonic systems are stable and most of the burning is localized at their core. A few cold cyclonic systems are still present between them and at the northmost latitudes, but the flame rapidly fills the gaps between the burning systems. Their centers slightly drift, "guided" by the background flow and the interaction with other subsystems, mostly toward the east. A more detailed discussion of the late evolution of the vortices, their drifts, and the effects on the observed light curves will be presented in future work.

Discussion
The first conclusion to be drawn is that, despite the quick and highly nonlinear evolution of the instability, the flame is able to survive until the end of the domain, which, as we said, is comparable to the equator-pole distance of a somewhat large star of 15 km of radius. We can then compare the perturbed flame results to the 2.5D simulation ones.
The first thing to compare is the global velocity of the flame. To measure the speed of the flame we have written a fronttracking tool that identifies the position of surfaces with a given value of the temperature and then returns the most advanced position. This is needed because the presence of the various anticyclones makes other methods very noisy. We set the tracking temperature to7 10 8 K, since this is coincident with  Figure 8. The comparison is striking: the unperturbed velocity is1.65 10 5 cm s −1 , while the perturbed one is1.72 10 6 cm s −1 . This means a speedup of ∼10 times. This is a very interesting result, because it can explain fast rise times for the bursts in stars with a weak magnetic field. The flame speed is comparable to the speed obtained by Cavecchi et al. (2016) in the case of a star with the same spin and a magnetic field of 10 8 G,1.61 10 6 cm s −1 , which was the highest achieved for a star with this spin. Since in the simulations of Cavecchi et al. (2016) the interaction between magnetic field and the flow of the fluid lead to significant changes, exploring the interaction of the instability and the vortices with the magnetic field seems a natural next step. Another important point we will explore is the dependence of the speedup on the spin of the star alone. We can expect the factor to depend on Ω since the mechanism leading to the speedup is the baroclinic instability and the presence of the vortices. The scale for both is set by the Rossby radius, which scales with 1/Ω, so that we can expect a similar scaling for the speedup. If we compare the total energy produced within the simulation domain (Figure 8, right) we can see that the increase of speed is reflected in the total energy released. This is to be expected and easily explained by the early ignition of the most distant fluid thanks to the almost horizontal convection and mixing triggered by the instability and the vortices. In Figure 8, left, the effects of the boundary conditions can be seen when the flame front reaches y∼23 km, and there is a clear change of the slope of the front position versus time. However, the fit for the speed is done in the interval before, where the boundary has no effect.
Second, we can extract azimuthal averages from the 3D run, so to have 2D profiles to compare to the 2.5D simulation. These averages are shown in Figure 9, where we plot the helium mass fraction Y, the temperature T, and reaction yield Q n . We compare the perturbed flame simulation (first column) to two situations in the unperturbed case: first to the 2.5D run at the same time (middle column) and then to the situation at a (much later) time when the flame has reached a position comparable to the one of the 3D simulation (last column). Note that we cannot calculate averages at every height, since we do not have data points above the top of the simulation and in 3D the top layer is highly corrugated by the vortices (as seen in Figure 7). For this reason, we calculate the averages up to 5.8 m, where data are available for all the x domain. However, we overplot the profile of the average height as a black solid line and the profiles of the minimum and maximum heights as dashed black lines to guide the eye.
First of all, some of the differences between the first and last columns are due to the fact that the time for the flame to reach ∼12 km is longer, so that near the origin the fluid had time to   cool, while the fluid ahead had been slightly burning. Then, looking at the helium mass fraction, we see that the vortices produced a much more mixed situation in the perturbed run, with more combustible fuel still present behind the flame. While the temperature is on average lower (a fact visible also in the lower average height of the fluid), a larger fraction of the fluid is burning and at deeper depth too. This is an important fact, especially comparing the first and last columns of Figure 9, since a larger part of the surface is hot. In terms of modeling observations, the 3D runs would predict shorter, but brighter bursts, as can be deduced also from the total energy release in the right panel of Figure 8. This can have deep implications when trying to extract from the observations the information about the nuclear reactions at play. For example, using the unperturbed run profile to fit short bursts would require much more efficient burning, changing the conclusions about the composition, for instance, or seemingly requiring higher yields than predicted by nuclear theory. This could have implications for understanding the observed burst rates, especially depending on the conclusions about the composition of the fluid, as suggested in Cavecchi et al. (2017). Looking at the curve of the maximal height in the left column of Figure 9, it is possible to identify two major "bumps," one approximately at = y 5 km and another one at approximately = y 12 km. It is tempting to consider whether the 3D average profile could be derived from the 2.5D run, since this would avoid expensive simulations if only global values were of interest. However, apart from the naive reduction of the timescales by a factor of ∼10, this seems a nontrivial task that is beyond the scope of this paper, and thus we will not treat it here.
Another important aspect of our results is related to the presence of the vortices. The instability that we observed radically changes the paradigm of flame propagation we have so far: our simulations show that a ring of fire initiating at the equator and propagating toward the poles will break down in vortices approximately the size of a few Rossby radii (we speculated on this in Spitkovsky et al. 2002). Likewise, if ignition takes place at midlatitudes, it is likely that the flame will not propagate as an expanding, perhaps drifting, smooth hurricane, but will instead spawn smaller vortices after its circumference exceeds a few Rossby radii and the flame velocities should be again ∼10 times faster. Furthermore, the local anisotropies (see in particular Figure 7, where the height profile is shown for the whole surface, and remember that it reflects the temperature as well) could lead to detectable burst oscillations.
Burst oscillations are detected when performing fast Fourier transforms of some burst light curves as a strong spike in the frequency domain (see Galloway et al. 2008;Watts 2012). The mechanism producing these oscillations is still debated (Heyl 2004;Cumming 2005;Piro & Bildsten 2005;Berkhout & Levin 2008;Heng & Spitkovsky 2009;Chambers et al. 2019); however, the general idea is that the surface emission pattern is not axially symmetric and the rotation of the star produces a modulation of the light curve that can be detected in the frequency domain. The frequency of the oscillations is often a few Hertz below the spin of the star and it is seen increasing in the cooling phase of the bursts. It is thought that the emission pattern must be moving counter-rotation-wise on the surface of the star to produce these effects. The most probable candidates for the mechanism of the oscillations are a hot spot (more likely in the rising phase of the bursts, while the flame is still propagating) or surface oscillation modes of the burning ocean. The modeling of the pulse profiles of the fluctuations attracts much attention, because they could be used to measure the effects of the strong star gravity on photon propagation, as done for the NICER detections (Gendreau et al. 2016;Bogdanov 2019), and thus put further constraints on the neutron star equation of state that determines the gravitational potential of the star.
In our simulation the vortices seem to drift along the residual flow of the fluid, which is toward the east, i.e., in the "wrong" direction, and, if detectable, would lead to an observed frequency above the spin frequency. A first glance at the drift in our simulation seems to indicate that there would be a positive difference of 1 Hz. However, we will calculate proper light curves in a follow-up paper. Furthermore, this simulation describes the rising phase of the bursts, where no substantial drift is observed. 5 When the fluid cools the flow should invert direction. If the vortices survive up until then, they would drift in the "correct" direction: when the equator is cooling, and the pressure profile along the north direction is inverted, vortices may drift west (see Spitkovsky et al. 2002). Also, even if the vortices themselves are not the direct origin of the burst oscillations in the tail of the bursts, 6 they could well be what excites the specific surface modes that produce the oscillations.
Finally, in the case of midlatitude ignition, the breaking up of the expanding hurricane into subsystems may weaken its signature as burst oscillations, and this would explain why they are so difficult to detect during the rise of the bursts. The expected launch of new telescopes like eXTP (Zhang et al. 2016) and STROBE-X (Wilson-Hodge et al. 2017) present an unprecedented opportunity. With their high collecting area, which also allows for fast time resolution during the rising part of the light curves, they may be able to detect the traces of these structures in future observations of the type I bursts and shed further light on the mechanism of the burst oscillations.