The imprint of convection on Type I X-ray bursts: Pauses in photospheric radius expansion lightcurves

Motivated by the recent observation by NICER of a type I X-ray burst from SAX J1808.4-3658 with a distinct"pause"feature during its rise (Bult et al. 2019), we show that bursts which ignite in a helium layer underneath a hydrogen-rich shell naturally give rise to such pauses, as long as enough energy is produced to eject the outer layers of the envelope by super-Eddington winds. The length of the pause is determined by the extent of the convection generated after ignition, while the rate of change of luminosity following the pause is set by the hydrogen gradient left behind by convection. Using the MESA stellar evolution code, we simulate the accumulation, nuclear burning and convective mixing prior to and throughout the ignition of the burst, followed by the hydrodynamic wind. We show that the results are sensitive to the treatment of convection adopted within the code. In particular, the efficiency of mixing at the H/He interface plays a key role in determining the shape of the lightcurve. The data from SAX J1808.4-3658 favors strong mixing scenarios. Multidimensional simulations will be needed to properly model the interaction between convection and nuclear burning during these bursts, which will then enable a new way to use X-ray burst lightcurves to study neutron star surfaces.


INTRODUCTION
As per the recent MINBAR catalogue (Galloway et al. 2020), about one fifth of type I X-ray bursts from accreting neutron stars (Lewin et al. 1993;Galloway & Keek 2021) reach high enough luminosities to provoke a radiatively-driven expansion of the neutron star envelope.In these "photospheric radius expansion" (PRE) bursts, the star's photosphere moves outward and appears 10-100 times larger, for a few to tens of seconds.PRE bursts offer a unique opportunity to study not only the surface but also the interior of neutron stars.Indeed, they have been used to place joint constraints on both neutron star mass and radius and the dense matter equation of state ( Özel et al. 2016, and references therein).Another way to constrain the mass is to measure the gravitational redshift of spectral features from heavy elements being ejected during the burst (Li et al. 2018;Strohmayer et al. 2019).These techniques rely on theoretical models which describe the expansion of the star's envelope and the winds that drive it.
The recent deployment of the Neutron Star Interior Composition Explorer (NICER) telescope has drastically improved observations of PRE bursts, since the instrument's soft X-ray response allows spectral evolution of the PRE to be followed as the blackbody temperature drops to ≲1 keV (Keek et al. 2018).A most interesting PRE burst was recently observed by NICER from the millisecond pulsar SAX J1808.4-3658(Bult et al. 2019).During the burst rise, the luminosity briefly "paused" for ≈0.7 s before reaching its peak.The ratio between the bolometric luminosity at the peak and pause was ≈1.68, which is very similar to the ratio between the pure helium and solar composition (X ≈ 0.7) Eddington luminosities, given by where M is the neutron star mass, σ T is the Thomson scattering cross-section, m p is the proton mass, and X is the mass fraction of ionized hydrogen (free protons).Bult et al. (2019) interpreted this as an observation of the rapid ejection of a solar or hydrogen-rich layer, followed by the usual helium PRE phase.This is consistent with the observed burst recurrence times and energetics which indicate that SAX J1808.4-3658 is in the burst regime where hydrogen is depleted by hot CNO burning well before unstable ignition of helium (Galloway & Cumming 2006;Goodwin et al. 2019).
A similar idea for a two-staged mixed H/He PRE burst was put forward by Sugimoto et al. (1984) to explain the bursting behaviour of 4U/MXB 1636-53, which showed a bimodal distribution of peak luminosity (see also Galloway et al. 2006).Following this suggestion, Kato (1986) computed steady-state solutions of outflows from a small H layer on top of a He-layer, finding the timescale for the ejection of the H layer (and jump in luminosity) to be on the order of 0.1 to 1 s, inversely proportional to the luminosity of the model.However, this strongly depends on not only the mass of the H layer, assumed to be 10 −16 M ⊙ , but also on the steady-state assumption, which cannot reproduce the actual ejection of the H layer.It is clear that in order to understand this type of burst, we need time-dependent hydrodynamic simulations combined with realistic neutron star envelopes.
X-ray bursts are challenging to model because of the many different types of physics involved.Previous studies have followed the time-dependent nuclear burning and convection during the thermonuclear runaway, using stellar evolution codes such as KEPLER (Woosley et al. 2004;Cyburt et al. 2010), or Modules for Experiments in Stellar Astrophysics, MESA, (Paxton et al. 2011;Meisel 2018), but did not resolve the formation of the wind in PRE bursts.Yu & Weinberg (2018) first demonstrated the ability of MESA to resolve both the nuclear burning and convective mixing at the onset of bursts from pure He accretion, followed by the hydrodynamic ejection of a super-Eddington wind and the PRE phase.This opened up the possibility of simulating time-dependent PRE bursts with an emphasis on the role of composition in the resulting wind.
In this paper, we use MESA (version 15140; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019) ) to simulate the accumulation phase, ignition, super-Eddington wind, and decay of a mixed H/He PRE burst.This represents the first full simulation of PRE bursts resulting from accretion of H and He.Our main result is that the rise in luminosity at the start of the burst pauses temporarily once the luminosity reaches the Eddington luminosity.As a wind develops and mass is ejected, deeper layers are eventually exposed that have been depleted in hydrogen by a combination of convective mixing and nuclear burning.This ends the pause and the luminosity begins to rise again as the outflowing material becomes less hydrogen rich and therefore has a larger Eddington luminosity.The resulting lightcurve has a distinct shape in contrast to pure helium bursts, and it depends on the gradient of hydrogen left behind by convection.However, as we will show, these results are very sensitive to the treatment of convection within the code.Indeed, the entrainment of protons by convection from the H layer into the He burning zone below leads to short timescale nuclear burning, which feeds back into the convection.This is a regime that cannot be adequately simulated in one dimension, making the detailed shape of light curve uncertain.
We begin in Section 2 with a simple model of mixed H/He PRE bursts, and explain the connection between the hydrogen profile in the envelope post-convection and the shape of the lightcurve.The main stages and important parameters of this model are also summarized in Figure 1.In Section 3, we describe our MESA simulations and show detailed results using the simplest prescription for convection.In Section 4, we vary the prescription for convection and assess how the results are affected.In Section 5, we summarize our findings, elaborate on issues related to the treatment of convection in one-dimensional simulations, and give an interpretation for the SAX J1808.4-3658data.

EVOLUTION OF THE COMPOSITION PROFILE AND LIGHTCURVE SHAPE
The structure of the neutron star envelope at ignition is determined mainly by the accretion rate onto the neutron star and composition of the infalling gas.For Ṁacc ≳ 2 × 10 −10 M ⊙ yr −1 , hydrogen burns via the hot CNO cycle at a constant rate (Bildsten 1998).Then, it can be shown that the column depth, y(r) ≡ r ∞ ρ(r ′ )dr ′ , at which hydrogen is depleted is (2) where X 0 and Z CNO are the initial hydrogen and CNO nuclei mass fractions (Cumming & Bildsten 2000) 1 .We scale the accretion rate to the Eddington accretion rate corresponding to a neutron star with R = 12 km and accreted hydrogen mass fraction X 0 = 0.7, giving ṀEdd = 8πRm p c/(1 Therefore at ignition, the envelope consists of two layers: an outer H-rich layer of depth y d in which the hydrogen abundance drops from the accreted value to zero, and an inner layer of pure helium where ignition of the burst  2).We are interested in bursts that ignite at y > y d (red color indicates nuclear burning).B) Heat from nuclear burning creates a growing convection zone (blue semicircles) which penetrates into the H-rich layer, resulting in additional nuclear burning.C) As the convection zone retreats, it leaves behind a layer of constant hydrogen fraction at a column yc,min < y d , a mixed H-He layer and nuclear ashes at depth.Winds progressively eject the layers, up to a column yw > yc,min, during which the observed luminosity at infinity is the Eddington luminosity, which depends on the hydrogen fraction X of the material.Since X is initially constant, we first observe a pause after the initial burst rise.After the H layer is ejected, the luminosity once again rises in a manner that depends on the hydrogen gradient dX/dy.See text for further details.
occurs.This initial state is illustrated in column A of Figure 1.
From an energetics standpoint, we know that only a column y w ∼ 10 6 -10 7 g cm −2 can be ejected by winds (Weinberg et al. 2006), which is smaller than y d .However, Weinberg et al. (2006) also showed that prior to the wind, a convection zone will grow and extend to a column smaller than y w (Figure 1 column B), which will mix the H and He and result in a change in the composition of the ejecta as a function of time.We define the minimum column depth reached by convection as y c,min , below which hydrogen is not mixed and X = X 0 is roughly constant 2 .
Wind models for PRE bursts (Ebisuzaki et al. 1983;Paczynski & Proszynski 1986;Joss & Melia 1987;Guichandut et al. 2021) show that the luminosity at infinity is always very close to the Eddington luminosity 2 X in fact decreases linearly with y, but since y c,min ends up being ∼1% of y d or less, the variation in X over this column is negligible.
L Edd (Equation 1), as any "extra" energy in the form of a super-Eddington flux gets used up to drive mass-loss.Therefore, during the initial ejection of y c,min , L ≈ L Edd will be constant.This is the observed pause, and its duration is where Ṁ18 = Ṁ /(10 18 g s −1 ) is the mass-loss rate, which we assume to be constant during the pause.After the pause, the ejection of the mixed layers will begin, and the luminosity will rise as the hydrogen fraction X in the ejecta decreases.The rate dL/dt at which the luminosity will increase, assuming it stays near Eddington, will be proportional to dX/dt ∝ Ṁ (dX/dy), thus linking the shape of the lightcurve to the hydrogen gradient in the envelope.Column C of Figure 1 illustrates the compositional nature of these two stages, the pause and the rise.Note that if X = 0 at columns y < y w , a third stage will appear where the luminosity peaks at the helium L Edd and remains there for the rest of the PRE (until all of y w has been ejected).
The burst lightcurve is therefore determined by 1) the nuclear burning and convection that occur during the rising phase, setting the hydrogen profile, and 2) the mass-loss rate during the wind phase.In the next section, we describe simulations with MESA to investigate both of these factors.

MESA SIMULATIONS
We model a single burst with several separate MESA runs 3 , in an approach similar to Yu & Weinberg (2018).First, we follow the ignition and convective rise of the burst under the assumption of hydrostatic equilibrium, then use MESA's hydrodynamic solver to follow the super-Eddington wind phase.However, unlike Yu & Weinberg (2018), we leave nuclear burning on during the wind, as much of the energy in bursts with hydrogen comes from slower reactions that continue into the wind phase.

Accumulation and formation of the layers
The basic physical setup is the same for all simulations: we assume a non-rotating neutron star of mass M = 1.4M ⊙ and radius R = 12 km, and ignore general relativistic corrections.The envelope initially consists of an 56 Fe substrate 4 with a column depth y = 10 11 g cm −2 .The outer boundary is initially set at an optical depth τ = 100 to avoid numerical issues caused by radiation-dominated regions becoming convective (Paxton et al. 2013) 5 .We begin accreting a solar-like composition ( 1 H, 4 He and 12 C with mass fractions X = 0.7, Y = 0.28 and Z = Z CNO = 0.02 respectively) at a constant rate of Ṁacc = 3 × 10 −10 M ⊙ yr −1 = 0.014 ṀEdd .We assume carbon to be the only metal being accreted for simplicity.What matters for the ignition of this type of burst is to achieve hydrogen depletion, and any isotope part of the CNO cycle would work, because the CNO abundances quickly adjust to the equilibrium ratio 3 Our MESA inlists, models and simulation setup are available at 10.5281/zenodo.8048553. 4To build the starting model, we took the the ns env model file from the ns he problem provided as part of MESA's test suite, then relaxed the neutron star radius from 10 to 12 km, and finally accreted additional 56 Fe to the target column depth.The point of increasing the mass of the iron substrate is to build a large enough buffer between the flashing zone and the inner boundary, allowing heat to diffuse inward without reflecting back. 5In this first part of the simulation, the outer layers do not expand significantly and we mainly focus on the effects of convection at large depths.We later relax this outer boundary to a more appropriate value to model the wind and lightcurve (Section 3.3).
of 14 O and 15 O in the hot CNO cycle (Bildsten 1998) 6 .Throughout accretion, the luminosity at the base of the substrate is fixed to 1.8×10 34 erg s −1 , equivalent to Ṁacc times 1 MeV per nucleon which is roughly the expected crust heating at low accretion rates (Brown 2000).As discussed in Yu & Weinberg (2018), Ṁacc and L base determine the ignition depth, which we have chosen to be at y ≈ 3 × 10 8 g cm −2 , greater than y d so that the burst ignites in a pure He layer.
To reduce the complexity of the computations, especially during the wind phase, we use MESA's cno extras.netnuclear network to model the nuclear burning, which includes 17 isotopes up to 24 Mg (we also add 56 Fe as an additional inert element).This is fewer and lighter isotopes than the approx21.netnetwork used by Yu & Weinberg (2018) for He bursts (which has α-capture reactions up to iron-group elements), but contains a similar number of reactions due to the addition of hot CNO.It is also a much smaller network than that used by Woosley et al. (2004), who studied the energy generation carefully but did not model the hydrodynamic wind.Since we are focused primarily on the wind, this limited network is adequate because most of the energy generated during the burst comes from hydrogen and helium burning in CNO and triple-α reactions.However, this assumption means that our calculations do not accurately predict the nuclear burning ashes, both ejected in the wind and leftover afterwards.
The left panel of Figure 2 shows the composition profile of the envelope after 5.9 days of accretion.The infalling CNO species convert to equilbrium oxygen ratios in under an hour.The hydrogen depletion column is y d = 3.7 × 10 7 g cm −2 , consistent with Equation 2 for the chosen Ṁacc .At larger column depths, the CNO cycle is starved of protons, and the abundances are determined by β-decay rates only.Prior to ignition, some helium has already started stably burning and converting to 12 C.

Ignition and convective rise
After having built up the He layer, unstable triple-α burning triggers the thermonuclear runaway.We define "ignition" as the moment t = t 0 when the He layer is convective and has a larger maximum nuclear energy generation rate than the H layer.The convection zone begins growing in the He layer as it would in a pure He burst, then hits the H layer about 1.1s later.The mixing of H into the convection zone leads to a sudden increase in the nuclear energy generation rate at the top of the zone.The outcome of this mixing event depends strongly on the prescription used in the code for convection.Here, we use the prescription of Henyey et al. (1965) for mixing length theory, with the dimensionless parameter α mlt = 1.5 dictating the ratio between the mixing length and local pressure scale height.In this section, we present results using the Schwarzschild criterion to determine convective boundaries (as assumed also by Yu & Weinberg 2018).This ignores the effects of composition gradients on the convective stability, but simplifies the interpretation of our results.We explore the effect of changing the prescription for convection in Section 4.
In Figure 3, we show Kippenhahn diagrams for the history of convection and nuclear burning as a function of depth, throughout the flash phase and beginning of the wind.The bottom panel is zoomed in to show short timescales following the collision between the He convection zone and the H layer.In this collision, fresh protons are brought in below the depletion depth y d where they can capture onto seed nuclei, causing a rapid injection of energy and a local steepening of the temperature gradient, which in turns drives further expansion of the convection zone.The first proton captures to happen are 18 O(p, α) 15 N and 15 N(p, α) 12 C.The remaining protons (and fresh ones coming from the top) then quickly capture onto the carbon and build up 13 N.The nuclear reactions in this first stage do not produce enough heat to generate large scale convection.Instead, the convection splits into many zones, as radiative gaps as small as 0.1% of the scale height appear.These zones and gaps are clearly seen in the bottom panel of Figure 3.The maximum number of individual convective zones is 66, and it occurs 0.1 ms after the collision.The maximal extent of the convection during this initial stage is to a column ≈10 6 g cm −2 .We later refer to this stage as the precursor.
About 0.4 ms after the initial collision, enough nitrogen has built up to trigger a "second ignition" via the 13 N(p, γ) 14 O reaction.This time, so much heat is released that the convection zone grows massively in less in 0.1 ms, its column depth extent decreasing by a factor of ≈30.The convection is still split but the radiative gaps now maintain over time, resulting in a period of layered convection, which extends down to a minimum column y c,min = 1.1 × 10 4 g cm −2 .The remnants of this period of rapid burning appear in the final composition profiles (see e.g. the magenta and orange lines for 14 O and 13 N respectively in the right panel of Figure 2), and will later be partly ejected by the wind.Further oxygen, neon and finally sodium-burning reaches the end-point of our network at 24 Mg.The final composition profiles that go into the hydrodynamic calculation are shown in the right panel of Figure 2. If a more complete network were to be used, we would expect the ashes to proton and α capture to heavier elements, namely Ca and Si, over the following ∼tens of seconds (Woosley et al. 2004).
At the onset of both of these stages, the precursor and 13 N ignition, convective velocities on the order of Figure 3. Kippenhahn diagrams for the Schwarzschild run, centered on the moment of collision between the convection zone and the H layer.The bottom panel is zoomed into a 1.5 ms window following the collision.The color scale traces the energy generation or loss (nuclear burning minus all neutrino losses), while the green hatches mark the convection zones.The solid black lines show the location of the depletion depth y d .In the top panel, the dashed line shows the luminosity coming out of the atmosphere, normalized by the Eddington luminosity of the accreted gas (scale on the right-hand side).In the bottom panel, it shows the integrated nuclear power ( dm ϵnuc).Moments where certain reactions dominate are labeled (see Section 3.2). 10 7 cm s −1 are generated, which allow for the overall convective envelope (disregarding gaps) to expand by meters in tens of microseconds.In most of the envelope, the large temperatures are such that these velocities remain well below the sound speed.However, in some convective zones, the local Mach number does reach up to ∼15%.This implies that the rapid burning is approaching a dynamical regime, but not so much as to make the hydrostatic assumption invalid.We return to the question of convective velocities in Section 4.
As we have explained in Section 2, the main predictor for the shape of the lightcurve is the hydrogen gradient left over by convection.We now investigate what creates this gradient.Figure 4 is another Kippenhahn diagram of the same simulation that shows the evolution of the hydrogen abundance after the collision.A few ms after the collision, once the convection retreats, the hydrogen gradient is already set.The dashed line traces how much hydrogen has burned away since ignition.At the collision, a significant amount of protons capture onto seed nuclei, and this continues throughout the burst even as the convection zone retreats.At ignition, the total mass of the envelope (minus the iron substrate) is ∼10 22 g, and ∼2% of it is H.By the time the gradient is set, about 25% of the hydrogen has burned away (another 10% burns in the rest of the flash, mostly at depths y > y w , before ejection by winds).What is important in determining the gradient is not the details of the nuclear burning; once convection subsides, hydrogen burning at shallow depths (near y c,min ) is slow and does not substantially affect the gradient and therefore the lightcurve.Instead, it is the efficiency of the convective mixing which determines how fast and how far protons can be brought downwards to regions of high temperatures where they can quickly burn away.We further investigate the point about the efficiency of mixing in Section 4.

Wind and collapse
As the outgoing luminosity rises and approaches the Eddington limit a short time after ignition, the outer layers become radiation pressure dominated and the envelope begins to expand.In Newtonian gravity, we know from previous work that appreciable expansion of ∼100 m above the stellar surface occurs at ≈90% of L Edd (see Figure 12 of Guichandut et al. 2021).At this point, we turn off convection7 and accretion, turn on MESA's hydrodynamics calculation, and relax the outer boundary to an optical depth τ = 2/3 in order to resolve the photosphere, as in Yu & Weinberg (2018).Mass-loss is then done by removing any grid points with a density ρ < ρ thresh = 10 −7 g cm −3 .To avoid issues caused by going off the opacity tables at low density, we switch to an interpolation formula for electron scattering opacity as a function of temperature from Paczynski (1983).This is a good approximation at the high temperatures of the wind where electron scattering dominates.
The presence of jumps in the 1 H mass fractions as a result of convection (see right panel of Figure 2) added some numerical difficulties in the simulation of the winds.Since the acceleration of a fluid parcel due to radiation is proportional to its opacity, a jump in the hydrogen fraction X will result in a density inversion, as the uppermost fluid element is ejected faster than the bottom one can follow.We found that these density inversions tended to expand as they moved outwards, and caused the MESA integration to diverge as they approached ρ thresh .A solution that worked in all cases was to manually soften those composition jumps using a smoothing spline on the 1 H mass fractions prior to running the wind simulation, as shown in the top panel of  .In order to preserve the overall gradient, this smoothing was done using monotonic functions.The 4 He profiles were then adjusted such that the sum of mass fractions of all species remained 1 everywhere.
During the wind, the rate of change of composition in the atmosphere is determined by the mass-loss rate where ρ and v are the gas density and velocity at a radial distance r from the center of the star.We evaluated Ṁ at three different locations: a) the "sonic point", i.e.where the velocity v = kT /µm p where µ is the mean molecular weight of the gas, b) the photosphere, i.e. the location r = r ph where the luminosity L = 4πr 2 ph σT 4 , and c) the surface of the model or outer boundary of the grid.As shown in Figure 5, despite some small variations, the mass-loss rate is overall constant across all locations.This is no surprise, as we expect these winds to reach a steady-state, characterized by a constant Ṁ (r) at fixed t, in a time much shorter than the evolution timescale of the burst (Joss & Melia 1987;Guichandut et al. 2021).We can then write the total column ejected as a function of time, independently of location.The total ejected column y w is the final value of y ej , equal to 5.83 × 10 5 g cm −2 in this simulation (see dashed line in Figure 5).This is roughly consistent with results obtained by Yu & Weinberg (2018) for the total ejected mass of pure He bursts igniting at a similar column depth as we have here.As L Edd (X = 0) L Edd (X = 0.7) Figure 6.Results of the Schwarzschild run.Top: Hydrogen profile after convective mixing.To avoid numerical difficulties during the hydrodynamic wind (see Section 3.3), the profile was smoothed using a monotonic cubic spline.The dashed lines show the minimum column reached by convection yc,min, and the total column ejected by winds yw.Bottom: Lightcurve of the burst.The pause occurs after surpassing L Edd of the accreted material (X = 0.7, bottom dotted line).The following rise takes place over a much longer timescale since yw≫ yc,min.In general, the outgoing luminosity follows the Eddington luminosity of the material which is currently being ejected (blue dotted line), which we can track using yej(t) (see Section 3.4).
shown by these authors, the total ejected column, and therefore the total duration of the PRE, would increase (decrease) for bursts which ignite at larger (smaller) column depths.
Approaching the end of the super-Eddington phase of the burst, the nuclear luminosity which is driving the wind begins to die down.The outflow then separates into two regions, an outer unbound wind of high velocity which is being ejected, and an inner atmosphere which is collapsing back into, eventually, a hydrostatic configuration.This can be seen from Figure 5 at t − t 0 > 9 s, where the mass-loss rate first drops sequentially from the inside (lower radii first).The timescale for the collapse to reach the surface from the sonic point is ≲1 s, which is roughly the sound crossing time between those two lo-cations (Guichandut et al. 2021).In some other simulations (Section 4), the collapse results in numerical issues as the infalling gas becomes supersonic, which causes the time-step to drop and the evolution of the model to come to a stop.The moment where these numerical issues arise coincides with the sonic point crossing the photosphere and the wind effectively becoming optically thin.Then, the implicit optically thick assumption made by MESA to treat the radiative transfer becomes incorrect.Future work is needed to investigate the infall phase in more detail.

Lightcurve
In order to plot the observed lightcurve, we evaluate the radiative luminosity of our models as a function of time at the photosphere.We show the lightcurve for the main Schwarzschild run in the bottom panel of Figure 6.Its shape is consistent with the hydrogen profile in the envelope post-convection, shown in the top panel.On the luminosity axis, the ratio between the peak and pause luminosities is ≈1.6, consistent with the ratio of Eddington luminosities (1 + X 0 )/(1 + X(y w )) with X(y w ) = 0.07 (Figure 6 top panel).On the time axis, we also have to take into account how the mass-loss rate varies throughout the burst.The pause duration is 0.6 s which, for y c,min = 1.1 × 10 4 g cm −2 , corresponds to a mass-loss rate Ṁ18 = 0.33 according to Equation 3.This is in good agreement with the time-averaged value of Ṁ during the pause, 3.5 × 10 17 g s −1 (note however that Ṁ changes significantly during the pause, from ≈6.5 × 10 16 g s −1 to ≈7.6 × 10 17 g s −1 , see Figure 5).If the mass-loss rate remained at this value throughout the remainder of the wind, the total duration of the PRE phase would be (y w /y c,min )∆t pause ≈ 32 s.But since Ṁ increases by a factor of ∼3 after the pause (Figure 5), the ejection is much faster and the PRE only lasts ≈9 s.
As expected, the lightcurve can be obtained by tracking L Edd as the hydrogen mass fraction in the ejecta X(y ej (t)) evolves in time.This is shown by the blue dotted line in Figure 6 which closely follows the luminosity from the simulation (black line).Comparing the two, we see that two features of the observed lightcurve are unexplained by composition changes.First, the luminosity during the pause is not exactly flat, but instead slowly decreases throughout its duration.This effect can be understood by considering the energetics of the expansion.At the beginning of the pause, the wind is not yet established.To do so, it needs to both lift material out of the gravitational potential and expand it (effectively raising its enthalpy).These two contributions account for the observed decrease in luminosity8 .Second, near the end of the super-Eddington phase and before the decay, a bump in luminosity appears.This is related to the wind dying down, which "returns" the gravitational energy and enthalpy that was required to sustain it back to the radiation.However, as discussed in the previous section, this part of the lightcurve is uncertain because of the wind becoming optically thin.
From the observational perspective, it could in principle be possible to infer the shape of the hydrogen profile from the lightcurve only.Assuming that the energy used to eject mass (GM M w /R where M w = 4πR 2 y w ) is equal to a fraction η of the observed burst energy (integrated luminosity which can be determined from the fluence if the distance to the source is known), one could infer y w from the lightcurve only.Then, given the total duration of the PRE, one could find the average Ṁ , and finally use Equation 3 to obtain y c,min .In our simulations, we find η ≈ 0.31 − 0.37, but this is likely to change for different ignition depths.We plan to study the energy budget of PRE bursts in more detail in future work.

IMPACT OF CHANGING THE TREATMENT OF CONVECTION
It is unlikely that the implementation of convection in Section 3 with the Schwarzschild criterion is an accurate representation of the true hydrodynamic phenomena.For one, the rapid nuclear burning induces local changes in composition towards heavier species, and locally increases the mean molecular weight µ.The creation of µ-gradients can either have a stabilizing or de-stabilizing effect on the thermal profile.This is especially relevant starting at the collision, where the growth of the He-rich convection zone will be inhibited by the H-rich material on top, as pointed out by Weinberg et al. (2006).They showed that a jump in temperature between the convection zone and the overlying radiative layer would develop in order to overcome the stabilization of the boundary due to composition.However, the effectiveness of the composition jump may be decreased by entrainment of fluid at the convective-radiative boundary which would erode the stabilizing composition gradient (e.g.Anders et al. 2022).Second, we may expect that the tiny radiative gaps that appear in between convection zones (see Section 3.2) will be destroyed by some form of overshoot mixing.This is because, at the Schwarzschild convective .Results for all convective prescriptions and spatial resolutions tested in this work.Left: Mass fraction of 1 H in the atmosphere after the convective and before the ejection by winds.The circles label the locations of yw resulting from each simulations.The black line shows the hydrogen gradient at ignition, from which all runs start.Right: Lightcurve of the bursts.The dotted lines mark the locations of L Edd for X = 0.7 and X = 0, as in Figure 6.Some simulations could not integrate through the decay phase and were stopped short at the "x" symbols.The inset zooms in on the pauses.boundary, the fluid parcel has zero acceleration by definition, but it may have enough inertia to continue rising and mix the fluid all the way to the next convective zone.
To better understand the importance of these effects, we ran additional simulations, starting from the same model at ignition (Figure 2, left panel), but changing the prescription for convection during the thermonuclear flash.We first ran a simulation using the Ledoux criterion instead of Schwarzschild for locating convective boundaries, which takes into account compositional gradients.When this is used, semiconvective and thermohaline mixing become available.For these, we set the dimensionless parameters α sc = 0.1, and α th = 2 (MESA uses these to determine diffusion coefficients, see Paxton et al. 2013), following the ns he test suite problem in MESA.Then, we ran simulations using both the Schwarzschild and Ledoux criteria where we also forced all radiative gaps of radial extent less than 10% of the scale height to close and become convective instead9 .This is meant to, in a simplified way, mimic the overshoot at the top of each convective zone.Finally, we tested the convergence of our simulations by running a high resolution version of every prescription.For these, the number of grid points during the flash phase was increased from ∼5000 to ∼15000 (the exact number varies in time and is adjusted by MESA using the mesh delta coeff control).
We show in Figure 7 the hydrogen composition profiles and lightcurves for all the simulations mentioned above.In each case, the main two features of our simplified model hold: 1) a larger y c,min results in a longer pause, and 2) a steeper hydrogen gradient leads to a faster rise in luminosity.The peak luminosity of the burst is then inversely proportional to X(y w ), where y w ≲ 10 6 g cm −2 is roughly constant across different bursts, and reaches the helium L Edd in half of the simulations.Moreover, lightcurves with higher peak luminosities have shorter PRE times -this is because the fluence is conserved for a given ignition depth, no matter the exact the shape of the lightcurve.The slow decrease of the luminosity during the pause is also observed in every case.For runs which managed to integrate through the decay phase (ones which do not end in an "x" in Figure 7), the bump at the end of the super-Eddington phase is also present.
Once y w has been completely ejected, all models join on the same exponential cooling track.
While Figure 7 demonstrates agreement with the basic model described in Section 2, it also clearly shows that the results, and in particular the observable lightcurve, depend on the prescription for convection.We cannot confidently claim that one prescription is more realistic than another, given the complex interactions between nuclear burning and mixing in this inherently multidimensional process, and so cannot predict lightcurves using these simulations.Nevertheless, we can assess the impact that different convective prescriptions have on the overall simulation.First, using the Ledoux criterion instead of the Schwarzschild criterion means that formerly convective regions become semi-convective instead, as the composition gradients stabilize the thermal profile.This effect inhibits the mixing, and, therefore (Section 3.2), reduces the amount of hydrogen burned away.This can be seen in the left panel of Figure 7 by comparing the orange and blue lines.Closing the radiative gaps between convective zones ("CG" in the figure legend) naturally has the opposite effect, with the mixing becoming stronger.The impact of spatial resolution is also interesting (compare solid and dashed lines in Figure 7).In the Schwarzschild and Ledoux runs, we obtain an increase in the total number of convective zones when increasing the resolution.Indeed, the addition of grid points in the convection zone allows it to split even more.Therefore, these runs are clearly not converged.However, in the CG runs, this splitting is effectively cancelled, as tiny zones are merged together at the end of every step, and we find good agreement between the models with different resolutions.
Another aspect to consider is the criterion for locating convective boundaries.In the presence of composition discontinuities, as is the case in our simulations, simplified applications of the Schwarzschild or Ledoux criteria can lead to unphysical scenarios and ultimately impede the growth of convective zones (see Gabriel et al. 2014 for an extensive discussion).In Paxton et al. (2018), an optional "predictive mixing" scheme was implemented in MESA to address this issue.This algorithm iteratively tests if the convective boundary is still correctly located once the material becomes fully mixed on the convective side, applying corrections otherwise.This may not be accurate in our situation where mixing also leads to rapid burning10 .Nonetheless, we ran some exploratory simulations of the flash including predictive mixing.The results were very similar in all cases, except when using the Ledoux criterion without closing gaps.Then, the hydrogen abundances were reduced (mixing was enhanced) compared to the original model without predictive mixing, but not all the way to the models with the Schwarzschild criterion.
Finally, as mentioned in section 3.2, we found some large convective velocities resulting from the rapid nuclear burning, implying a departure from a hydrostatic envelope.We also find similar velocities in the other simulations presented in this section.However, the usual implementation of mixing-length theory does not factor in the time required to accelerate to these velocities.To explore this, we evaluated the timescale for convective acceleration at any given step n, given by Wood (1974) as ), where ℓ is the mixing length and v c is the convective velocity.If the actual timestep taken by the code ∆t = t n+1 − t n is shorter than τ a , then velocities have increased too quickly.In our simulations, we found a few instances where ∆t/τ a ∼ 10 −3 , but only during the precursor phase.This casts doubt on the existence and behaviour of this precursor.Alternatively, MESA has the option to enable acceleration-limited convection (Paxton et al. 2015).The most recent version also includes a more complete theory of time-dependent convection (Jermyn et al. 2023) which will also limit accelerations.However, our attempts at running simulations with either option were unfruitful, with the evolution eventually being driven to prohibitively short timesteps in trying to converge.Getting these models to converge and reach Eddington is a potential avenue for future work, as they are likely the most accurate way to model convection with rapid burning while remaining in one dimension.

SUMMARY AND DISCUSSION
We have shown that variations in chemical composition in the envelope of neutron stars accreting mixed H/He fuel are reflected in the lightcurves of their PRE bursts.After the ignition of the thermonuclear runaway in the He-rich layer, a convective zone expands outwards and mixes the fuel with the overlying H-rich shell (Figure 2).The resulting H abundance profile determines the shape of the lightcurve, namely the duration of initial pause and the subsequent slope in luminosity (Figure 6).Due to convection, the mass of the layer with solar composition, y c,min ∼ 10 4 − 10 5 g cm −2 , is much reduced compared to the initial y d ∼ 10 7 g cm −2 set by stable hydrogen burning during accretion (Equation 2).This results in a rapid ejection of a hydrogen-rich shell and a short observed pause on the order of 1 s or less (Equation 3).Subsequently, the luminosity rises toward the helium Eddington luminosity as hydrogen depleted layers are exposed by the wind.
We find that the hydrogen profile in the envelope is sensitive to the details of convection and mixing following the collision with the H layer (Figures 3 and 4).As a result, the exact shape of the lightcurve of a given event is uncertain as it depends on the choice of convective prescription and spatial resolution (Figure 7).The critical factor in setting the hydrogen profile is the efficiency of the mixing within the convective regions.However, this mixing is inhibited when the convection splits into many zones interspersed with radiative gaps.This splitting occurs even when ignoring compositional gradients (Schwarzschild criterion), suggesting that the culprit is the local energy deposition from rapid nuclear burning.Moreover, we found that increasing the spatial resolution of the simulations led to an increase in the number of zones and gaps, significantly reducing the efficiency of mixing, such that our simulations are not converged.This non-convergence however is mitigated by overshoot mixing at the top of convective zones, which we modeled in a simplified way by closing radiative gaps less than 10% of the scale height.
Even disregarding problems related to splitting of the convection, a more fundamental issue stems from the approximate treatment of convection with mixing-length theory.In the collision event, nuclear burning releases energy on tens of microsecond timescales, which is close to or even shorter than local convective turnover times.This is in violation of the standard assumptions of mixing length theory.It also amplifies the differences between the Schwarzschild and Ledoux criterion, in contrast to situations with long dynamical timescales such as the main-sequence, where both criteria should lead to similar outcomes (Anders et al. 2022), although we did find that an improved implementation of convective boundaries using predictive-mixing brought Ledoux closer to Schwarzschild.This timescale problem has been noted before in the context of late-stage evolution of massive population III stars.There, a helium burning convective region encroaches upon a hydrogen shell, mixing in protons which burn on timescales of hours to days, which is short compared to the monthlong convective turnover times (Marigo et al. 2001).The proper modeling of these situations, which are also known as level-3 mixing or convective-reactive phases (Herwig et al. 2011), continues to be an active area of research (e.g.Davis et al. 2019;Clarkson & Herwig 2021), with a particular focus on multidimensional hydrodynamics simulations (Woodward et al. 2014;Stephens et al. 2021).
Although the general shape of the lightcurves in our simulations agrees with the burst from SAX J1808.4-3658 reported by Bult et al. (2019), there are some differences.In the observed burst, the pause is ∼0.7 s long, similar to our principal Schwarzschild run, suggesting a similar extent of the convection.However, the subsequent rise is rapid, reaching the helium Eddington luminosity in just ∼1.3 s.This would imply a mixing event which is strong enough to produce a very steep hydrogen gradient.Since the total pause plus rise duration is ∼3 times that of the pause, the hydrogen profile would have to go from 11 X ≈ 0.7 to X = 0 in the span of y c,min ≈ 10 4 g cm −2 to ≈3 y c,min .Or, if for example the mass-loss rate increases by a factor of 3 from the pause to the rise, as it does in our simulations (Section 3.3), then in the span of y c,min to ≈9 y c,min ; in any case, the whole hydrogen gradient spans a decade in column depth at most.None of our simulations achieve this -our fastest rise time is ≳6 s for the Schwarzschild+CG model (in fact, factoring redshift, these times should be ∼20−30% longer).Moreover, our Figure 7 shows a general trend that steep hydrogen gradients are also associated with smaller y c,min values; to reproduce the rise seen in SAX J1808.4-3658,we may need such strong mixing that it would push y c,min to very small columns and dissolve the pause entirely.
One way to match the rapid rise but keep the pause duration the same as observed in SAX J1808.4-3658 could be to burn hydrogen more effectively with the same mixing efficiency and convective extent.In fact, our simulations do not model hydrogen burning completely, because we were limited to a small nuclear network which reached its end at 24 Mg prior to the wind launch (see right panel of Figure 2).We investigated the effect of a larger network by running a simulation equivalent to our Schwarzschild run, but using MESA's rp 153.net nuclear network, which includes isotopes up to 56 Ni, until the wind launch.We found that the outer hydrogen profile was unchanged, with y c,min and initial hydrogen gradient staying the same.The effects of the additional hydrogen burning were limited to large columns ≳ 0.3y w .At these depths, hydrogen completely burned away, whereas with the smaller network a small amount (≲ 0.1 for the well-mixed models) remains.The lightcurve for such a burst would initially look the same as in our original Schwarzschild run (Figure 6 bottom 11 Goodwin et al. (2019) inferred a hydrogen mass-fraction X 0 ≈ 0.57 +0.13 −0.14 for the companion, based on an analysis of Type I Xray burst recurrence times and energetics.The observed ratio of peak to pause luminosities in Bult et al. (2019) favors the upper end of this range.panel), but would continue rising all the way to the helium Eddington luminosity, instead of levelling of to ∼90% of it.This suggests that additional burning is not the explanation for the rapid rise after the pause.
More observations of PRE bursts in the pure helium ignition regime would help to further understand and constrain the hydrogen ejection model.Note that previous PRE bursts from SAX J1808.4-3658,observed with the Rossi X-ray Timing Explorer, have shown an increase in luminosity during the Eddington phase, but only by ∼20% and on ∼7 s timescales (see Figure 3 in Galloway et al. 2017).In these bursts, pauses are not clearly seen, which would indicate small values of y c,min , although this could also be due to the choice of time bins used for the analysis.Such variations in the shape of the lightcurve (slow or fast luminosity increases) across different bursts from a single source may also imply that the dynamics of convection are very sensitive to initial conditions at ignition.Furthermore, a puzzling aspect of the burst reported in Bult et al. (2019) is the secondary peak following the PRE phase.This is unexplained by our hydrogen ejection model, and could instead require multidimensional effects.
On the theoretical side, the obvious next step in order to refine predictions for these bursts will be to improve the treatment of convection during the thermonuclear flash, in particular for the collision between the He and H layers. Due to the timescales involved and the limitations of mixing length theory, we know that only multidimensional hydrodynamical simulations can yield accurate results.This may pose a significant numerical challenge, although recent works by Malone et al. (2014) and Zingale et al. (2015) have shown promising results in this direction, demonstrating the use of low-mach number hydrodynamics to model two and three-dimensional convection in thermonuclear explosions.
Finally, other improvements need also be made in the hydrodynamical part of the simulation in order to correctly model the super-Eddington wind.First, we faced numerical problems with "staircases" in mass fractions leading to density inversions in the wind, which we simply smoothed out in this work.It would be interesting to investigate such density inversions as they propagate outward in future work.We also had issues at the end of the super-Eddington phase and collapse of the atmosphere.To properly model this part of the PRE, we will likely need hydrodynamical simulations which can handle optically thin radiative transfer as well as shocks (if our findings that infall velocities can be supersonic are correct).Hydrodynamical simulations would also be useful to model the super-Eddington winds in multiple dimensions, where the effects of rotation and magnetic fields could be taken into account.Lastly, for accurate observational predictions, it would be pertinent to include general relativistic corrections to the hydrodynamic equations, as they are known to result in larger photospheric radii (Paczynski & Proszynski 1986;Guichandut et al. 2021).

Figure 1 .
Figure 1.Diagram illustrating the three main stages of the mixed H/He burst.A) As hydrogen burns stably throughout accretion, distinct hydrogen and helium-rich layers build up in the envelope, their boundary being at a known column depth of y d (Equation2).We are interested in bursts that ignite at y > y d (red color indicates nuclear burning).B) Heat from nuclear burning creates a growing convection zone (blue semicircles) which penetrates into the H-rich layer, resulting in additional nuclear burning.C) As the convection zone retreats, it leaves behind a layer of constant hydrogen fraction at a column yc,min < y d , a mixed H-He layer and nuclear ashes at depth.Winds progressively eject the layers, up to a column yw > yc,min, during which the observed luminosity at infinity is the Eddington luminosity, which depends on the hydrogen fraction X of the material.Since X is initially constant, we first observe a pause after the initial burst rise.After the H layer is ejected, the luminosity once again rises in a manner that depends on the hydrogen gradient dX/dy.See text for further details.

Figure 2 .
Figure 2. Composition profiles before (left) and after (right) the thermonuclear flash, which ignites at t = t0.The lines for different elements have the same color in left and right panels.In the ∼1.4 s to reach the Eddington luminosity, convection has significantly mixed the He and H layers. Dotted lines show the depletion column y d (left) and the minimum extent of convection yc,min (right).All runs shown in this paper begin at the ignited model t = t0, which the left panel leads up to.The panel on the right shows the result of mixing using the Schwarzschild criterion prescription for convective boundaries (see Section 3.2).

Figure 4 .
Figure 4. Kippenhahn diagram for the Schwarzschild run as in Figure 3, now with the color scale representing the hydrogen abundance.The solid black and gray lines show y d and yc,min (same values as in Figure 2).The dashed line shows the relative change in the total mass of hydrogen at a given time compared to the initial amount at ignition (scale on the right-hand side).

Figure 5 .
Figure5.The mass-loss rate of the wind as a function of time for the Schwarzschild run, in solid lines.The different colors show Ṁ evaluated with Eq.(4) at different locations (see text).The dashed lines show the ejected column yej (scale on the right-hand side).The red shaded region marks the pause in the lightcurve (bottom panel of Figure6).

Figure 6
Figure6.In order to preserve the overall gradient, this smoothing was done using monotonic functions.The 4 He profiles were then adjusted such that the sum of mass fractions of all species remained 1 everywhere.During the wind, the rate of change of composition in the atmosphere is determined by the mass-loss rate

L
Edd (X(y ej (t))) Figure7.Results for all convective prescriptions and spatial resolutions tested in this work.Left: Mass fraction of 1 H in the atmosphere after the convective and before the ejection by winds.The circles label the locations of yw resulting from each simulations.The black line shows the hydrogen gradient at ignition, from which all runs start.Right: Lightcurve of the bursts.The dotted lines mark the locations of L Edd for X = 0.7 and X = 0, as in Figure6.Some simulations could not integrate through the decay phase and were stopped short at the "x" symbols.The inset zooms in on the pauses.