A"Hyperburst"in the MAXI J0556-332 Neutron Star: Evidence for a New Type of Thermonuclear Explosion

The study of transiently accreting neutron stars provides a powerful means to elucidate the properties of neutron star crusts. We present extensive numerical simulations of the evolution of the neutron star in the transient low-mass X-ray binary MAXI J0556--332. We model nearly twenty observations obtained during the quiescence phases after four different outbursts of the source in the past decade, considering the heating of the star during accretion by the deep crustal heating mechanism complemented by some shallow heating source. We show that cooling data are consistent with a single source of shallow heating acting during the last three outbursts, while a very different and powerful energy source is required to explain the extremely high effective temperature of the neutron star, ~350 eV, when it exited the first observed outburst. We propose that a gigantic thermonuclear explosion, a"hyperburst"from unstable burning of neutron rich isotopes of oxygen or neon, occurred a few weeks before the end of the first outburst, releasing 10^44 ergs at densities of the order of 10^11 g/cm^3. This would be the first observation of a hyperburst and these would be extremely rare events as the build up of the exploding layer requires about a millennium of accretion history. Despite its large energy output, the hyperburst did not produce, due to its depth, any noticeable increase in luminosity during the accretion phase and is only identifiable by its imprint on the later cooling of the neutron star.

as "deep crustal heating" (Brown et al. 1998) which, over a long time, will lead to an equilibrium between this heating and photon and neutrino cooling mechanisms (Miralda-Escude et al. 1990;Brown et al. 1998;Colpi et al. 2001). Depending on the star's core temperature, the timescale to establish this equilibrium can range from a few decades for an initially very cold core, up to millions of years for the hottest stars (Wijnands et al. 2013). Regarding evolution on short timescales, theoretical modeling found that in the case of a long and strong enough accretion outburst the crust can be driven out of thermal equilibrium with the core (Rutledge et al. 2002), which led to the prediction that subsequent cooling of the crust, once accretion has stopped and the surface temperature is not anymore controlled by the mass accretion but rather by the internal evolution of the star, should be observable on a timescale of a few years. This prediction has been amply confirmed and to date crust cooling after an accretion outburst has been observed in almost a dozen cases, which we reproduce in Figure 1.
After observations of the cooling of the first quasipersistent source that went into quiescence (KS 1732-260: Wijnands et al. 2001;Cackett et al. 2006), theoretical modeling found it was not possible to reproduce the high temperature of the first data point, obtained two months after the end of the outburst, within the deep crustal heating scenario (Shternin et al. 2007) and the introduction of another energy source, located at low densities and dubbed "shallow heating", was found to be necessary in the models (Brown & Cumming 2009). The deep crustal heating generates an energy Q dc of about 1.5 to 2 MeV per accreted nucleon (Gupta et al. 2008;Haensel & Zdunik 2008;Fantina et al. 2018;Shchechilin et al. 2021), most of it through pycnonuclear fusions in the inner crust at densities above 10 12 g cm −3 . In contrast, shallow heating has been found to deposit energy at densities well below 10 11 g cm −3 but its strength Q sh , as well as the depth at which it is deposited, appear to vary significantly from one source to another. Typical values needed for Q sh are of the order of 1-3 MeV, but in some cases (XTE J1701-462, Page &Reddy 2013, andSwift J174805.3-244637, Degenaar et al. 2015) it was not found to be required, while in the case of MAXI J0556-332 (Deibel et al. 2015) it could be well above 10 MeV.

MAXI J0556-332
In the present work we focus on the thermal evolution of the neutron star in MAXI J0556-322, an X-ray transient that was discovered in early 2011 (Matsumura et al. 2011) when it had started a major outburst. It was in outburst for 16 months before it returned to quiescence in May of 2012. Although no thermonuclear bursts or pulsations were detected, the behavior of the source in an X-ray hardness-intensity diagram strongly suggested that the accreting compact object in MAXI J0556-332 is a neutron star: the tracks traced out at the highest fluxes showed a strong resemblance to those of the most luminous neutron star low-mass X-ray binaries, the Z sources. Based on a flux comparison with other Z sources, Homan et al. (2011) suggested that the source is a very distant halo source. MAXI J0556-332 has shown three smaller outbursts since its discovery outburst, in late 2012, 2016, and 2020. Homan et al. (2014) studied the cooling of the neutron star after the first two outbursts. They found that during the first few years in quiescence after the 2011/2012 outburst, the neutron star in MAXI J0556-332 was exceptionally hot compared to the other cooling neutron stars that have been studied, as can be seen in Figure 1. The smaller second outburst that started in late 2012 and lasted ∼115 days did not produce detectable deviations from the cooling trend seen after the first outburst. Deibel et al. (2015) showed that to produce the extremely high temperatures observed after the first outburst, a high amount of shallow heating was required (∼6-16 MeV per accreted nucleon). They fur-ther concluded that the shallow heating mechanism did not operate during the second outburst. Parikh et al. (2017a) analyzed the neutron star cooling observed after the first three outbursts. Reheating of the crust was observed after the third outburst. It was concluded that the strength of the shallow heating in MAXI J0556-332 varied from outburst to outburst. For two other sources the cooling has been studied after multiple outbursts as well. In MXB 1659-29 the strength of the shallow heating was found to be constant in the two outbursts after which cooling was studied. For Aql X-1, with numerous but short outbursts, the shallow heating was found to differ in both strength and depth between outbursts (Degenaar et al. 2019). A definitive explanation for a varying shallow heating in MAXI J0556-332 and Aql X-1 has not been provided yet.
Here we present a re-modelling of the crustal cooling data of MAXI J0556-332, including data taken after the end of the most recent outburst from 2020. We argue that the high crustal temperatures caused by the 2011/2012 outbursts were not the result of anomalously strong shallow heating, but were instead caused by a gigantic thermonuclear explosion that occurred at some time during the last three weeks of that outburst. As this explosion was much more energetic than any other thermonuclear X-ray burst observed to date, even about a hundred times more powerful than a superburst, we find it appropriate to call it a "hyperburst". We infer it must have been produced by unstable thermonuclear burning of neutron rich isotopes of oxygen or neon.
The paper is organized as follows. In Section 2 we present the analysis of the data taken during and soon after the 2020 outburst of MAXI J0556-332. In Section 3 we briefly describe how we model the neutron star temperature evolution and in Section 4 we compare two different scenarios of shallow heating. In section Section 5 we propose the new scenario of a hyperburst as the cause of the hot crust of MAXI J0556-332 when it exited the first 2011-2012 observed outburst and in Section 6 we try to identify the source of this explosion. We discuss our results in Section 7 and conclude in Section 8.

NICER observations
The 2020 outburst of MAXI J0556-332 was monitored extensively with the X-ray Timing Instrument (XTI) on board the Neutron Star Interior Composition Explorer (NICER; Gendreau et al. 2016). The XTI provides coverage in the 0.2-12 keV band and consists of 56 nonimaging X-ray concentrators that are coupled to Focal Plane Modules (FPMs), each containing a silicon drift Each data point represents a single GTI. Background was not subtracted. The start of the outburst (as observed with MAXI) and the occurrence of a small X-ray burst are marked with dashed orange and green lines, respectively. (b) NICER (black) and Swift (red) coverage of the decay into quiescence. Each data point represents the background subtracted 0.3-10 keV count rate in Crab units. The blue line shows the estimate of the end of the outburst (MJD 59165.6), obtained by finding the intersection between an exponential fit to the last four NICER and first two Swift data points, and a constant fit to the last five Swift data points.
detector. At the time of the observations 52 of the 56 FPMs were functional. A total of 107 ObsIDs are available for the 2020 outburst, each with one or more goodtime intervals (GTI). Data were reprocessed with the nicerl2 tool that is part of HEASOFT v6.29. Default filtering criteria 1 were used. We additionally used the count rate in the 13-15 keV band, where no source contribution is expected, to filter out episodes of increased background (13-15 keV count rates >1.0 counts s −1 per 52 FPMs). After filtering, a total exposure of ∼332 ks remained.
A full 0.5-10 keV outburst light curve was made with one data point per GTI. The GTIs varied in length from 16 s to ∼2.1 ks. The outburst light curve is shown in Figure 2a. NICER observations started six days after the first MAXI/GSC detection of source activity (Negoro et al. 2020, see orange dashed line in Figure 2) and caught the source during a fast rise towards the peak of the outburst. After the peak of the outburst the NICER count rate dropped rapidly for ∼20 days, after which the count rate decreased more slowly until the end of the outburst, about half a year later.
We also produced 0.5-10 keV light curves for each ObsID with a time resolution of 1 s to search for possible X-ray bursts. One candidate X-ray burst was found in ObsID 3201400198 (MJD 59147, 2020 October 25). A segment of the light curve of the observation is shown in Figure 3. The burst is very short (<10 s from start to end) and peaks at a count rate a factor ∼3 higher than the persistent emission. Fits to the spectrum of the persistent emission slightly favor a soft spectral state at the time of the burst. Fitting a "double-thermal" model (Lin et al. 2007) with n H fixed to 3.2×10 20 cm −2 (Parikh et al. 2017b) gives black-body and disk-black-body temperatures of 1.45(1) keV and 0.42(1), respectively, with a reduced χ 2 of 1.00 (a typical hard state model, blackbody with power-law, resulted in a slightly worse reduced χ 2 of 1.07). The unabsorbed 0.5-10 keV flux was 5.8×10 −11 erg cm −2 s −1 , corresponding to a luminosity of 1.3×10 37 erg s −1 for a distance of 43.6 kpc (Parikh et al. 2017b). Careful inspection of the data suggests that the burst is not the result of a short, sudden increase in the background. The burst is also present in various subsets of detectors, indicating that the burst was not instrumental. The count rates during the burst are too low to perform a study of the spectral evolution, but the burst was more pronounced at higher energies (>2 keV), which is consistent with a thermonuclear (type I) X-ray burst. We used WebPIMMS 2 to convert the peak count rate (75 counts s −1 ) into a 0.01-100 keV luminosity (assuming a black-body temperature of 2.5 keV, a distance of 43.6 kpc, and an n H of 3.2×10 20 cm −2 ): ∼ 3.8 × 10 38 erg s −1 , which is very close to the empirical Eddington limit of neutron stars (Kuulkers et al. 2003). This would constitute the first detection of a type I X-ray burst from MAXI J0556-332 and would add further proof that the compact object is a neutron star.
Background subtracted spectra were produced for each ObsID using the nicerbackgen3C50 tool (Remillard et al. 2021), which uses several empirical parameters to construct background spectra from a library of observations of NICER background fields. FPMs 14 and 34 were excluded, since they tend to be affected by noise more often than the other detectors. We note that back- ground spectra could only be extracted for 76 of the 107 ObsIDs; for the remaining ones the empirical parameters fell outside the range covered by the library background observations. From the resulting spectra a backgroundsubtracted 0.3-10 keV light curve was produced. In §2.2 this light curve will be used in conjunction with Swift data to estimate the end of the outburst and the start of quiescence.

Swift observations
The Neil Gehrels Swift Observatory X-Ray Telescope (Swift-XRT; Burrows et al. 2005) observed MAXI J0556-332 eight times near the end of the 2020 outburst. A 0.3-10 keV light curve with one data point per ObsID was produced by the online Swift-XRT products generator 3 (Evans et al. 2007(Evans et al. , 2009). Background was not subtracted, since the corrections were expected to be only minor (∼5%). The Swift light curve (red data points) is shown together with the NICER backgroundsubtracted data (black) in Figure 2b. Both the Swift and NICER data were normalized to Crab units.
The NICER data in Figure 2b cover a switch from a slow to a rapid decay, around MJD 59153. The Swift data covered the end of the decay and the start of quiescence, as indicated by the levelling off of the count rates. To estimate the start of quiescence (t 0 ) we determined the intersection of the decay (modeled with an exponential) and the quiescent level (modeled with a constant). A similar method was previously employed for other sources (see Fridriksson et al. 2010;Waterhouse et al. 2016;Parikh et al. 2017bParikh et al. , 2019. The exponential fit was made to the last four NICER and the first two Swift data points. The constant fit was made to the last five Swift data points. The third Swift data point was excluded since the higher count rate indicated a possible short episode of enhanced accretion, which was also ob-served several times soon after the end of the 2011/2012 outburst (Homan et al. 2014). From the fits we find t 0 =MJD 59165.6. This time is marked with a blue line in Figure 2b, along with the exponential (dashed line) and constant fits (dotted line).
Again using the online Swift-XRT products generator we extracted a single background-subtracted quiescent spectrum for the last five Swift observations. The spectrum had an exposure time of ∼12.6 ks, contained ∼30 counts, and was fit with XSPEC v12.12 (Arnaud 1996) using C-statistics. We used a similar absorbed neutron star atmosphere model as in Homan et al. (2014) and Parikh et al. (2017a): constant × tbabs × nsa. Note that these authors employed nsa instead of the frequently used neutron star atmosphere model nsatmos (Heinke et al. 2006), since the latter model does not allow temperatures above log 10 (kT ) = 6.5. constant is set to a value of 0.92 to account for cross-calibration with other X-ray detectors. Although only a Swift spectrum is analyzed in our work, the n H and distance used in our fits are based on the Chandra and XMM-Newton spectra fitted by Parikh et al. (2017a); not using the cross-calibration constant would result in a temperature that is too low. The absorption component tbabs is an updated version of the tbnew feo component used by Homan et al. (2014), but it yields the same column densities. We set the abundances to WILM and cross sections to VERN, and fixed the n H to the value obtained by Parikh et al. (2017a): 3.2×10 20 cm −2 . For the neutron star atmosphere component nsa we fixed the neutron star mass (M ) and radius 4 (R) to 1.4 M and 10 km, respectively, and for the normalization we used a distance of 43.6 kpc. The only free parameter in the fit was the temperature of the nsa component. When converted to the effective temperature measured at infinity 5 we find a value of 166±8 eV. This data point (for which we use the midpoint of the five observations as date: MJD 59176.6) was added to the cooling data obtained by Parikh et al. (2017a) from prior Swift, Chandra, and XMM-Newton observations.

MODELING THE EVOLUTION OF THE MAXI J0556-332 NEUTRON STAR
We build on our previous experience in modeling accretion heated neutron stars using an updated version of the code NSCool (Page 2016), which solves the gen-4 In this paper by radius we always mean the coordinate radius in Schwarzschild coordinates, i.e., not the "radius at infinity". 5 T ∞ eff = T eff /(1 + z), where (1 + z) ≡ (1 − R S /R) −1/2 is the gravitational redshift factor, with R S = 2GM/c 2 being the Schwarzschild radius. eral relativistic equations of stellar structure and thermal evolution. The physics we employ is standard and briefly described in Appendix A.
To model the mass accretion rate during the four observed outbursts we apply a similar method as previously described in Ootes et al. (2016Ootes et al. ( , 2018 and Parikh et al. (2017a). However, instead of using a mix of Swift, NICER and MAXI data to estimate the evolution of the bolometric flux, we opt to only use the MAXI data, since they provide a single consistent set that covers all four outbursts and requires just one count-rate-to-flux conversion factor. For the 2-20 keV MAXI count rate to bolometric flux conversion factor we used the value from Parikh et al. (2017a): 2.353×10 −8 erg cm −2 count −1 .
From the bolometric flux F bol we calculate the daily average mass accretion rate usinġ where D is the source distance, fixed here at D = 43.6 kpc, c the speed of light and η ≈ 0.2 the accretion efficiency factor (Shapiro & Teukolsky 1983). Heating was assumed to only take place during the outburst intervals, which were fixed by hand and are shown as gray shaded areas in Figure 6. The temporal evolution of the heating luminosities, H dc (t) for the deep crustal heating and H sh (t) for the shallow heating, follow the instantaneous mass accretion rate as with their respective strengths Q dc and Q sh . We follow Haensel & Zdunik (2008) and use Q dc = 1.93 MeV that is distributed in density at the locations of each reaction (see figure 3 in Haensel & Zdunik 2008) for the model starting with 56 Fe at low densities. For the shallow heating we assume Q sh is distributed uniformly, per unit volume, within a shell covering a density range from ρ sh up to 5ρ sh . Both Q sh and ρ sh are considered, in a purely phenomenological manner, as free parameters that are adjusted to fit the data. In our model of a thermonuclear hyperburst we assume an energy X hb × 10 18 erg g −1 is deposited almost instantaneously at time t hb in a shell going from our outer boundary layer (at ρ = ρ b = 10 8 g cm −3 ) up to a density ρ hb . We envision the layer at density ρ hb as a "critical layer" and the latter as a "critical density", being the point where the hyperburst is triggered. The temperature of this critical layer at the moment when the hyperburst is triggered, T hb , is not fixed in any way but will be an output of our simulations and will be referred to as a "critical temperature". X hb , t hb , and ρ hb , are also treated phenomenologically as free parameters to fit the data. An energy of 10 18 erg g −1 corresponds to about 1 MeV per nucleon and is a typical energy released by fusion of C, O or Ne into iron-peak nuclei (see Table 1) so that X hb roughly represents the mass fraction of the exploding nuclear species.
In the previous study of the evolution of MAXI J0556-332 in Parikh et al. (2017a) a χ 2 minimisation technique was used to obtain the best fit to the data. Here, however, we apply the more robust method, for models with many parameters, of Markov Chain Monte Carlo (MCMC) using our recently developed MCMC driver MXMX (see Lin et al. 2018;Ootes et al. 2019;Degenaar et al. 2021). The parameters we vary are the mass M and radius R of the star, which determine the crust thickness that is a dominant factor in controlling, on one hand, the time scale for heat transport and, on the other hand, the amount of matter present and thus the amount of energy needed to heat it 6 , the initial redshifted core temperature T 0 , the depth of the light element layer in the stellar envelope, y L , which affects the observed effective temperature T eff , the impurity parameter Q imp that can strongly reduce the thermal conductivity in the solid crust, the strength and depth of the shallow heating, Q sh and ρ sh , and, when considered, the properties of the hyperburst, X hb , t hb , and ρ hb . In Appendix B we describe the various settings we apply in the various scenarios presented in the following sections.

SHALLOW HEATING IN THE FOUR OBSERVED ACCRETION OUTBURSTS
To extend the work of Parikh et al. (2017a) and incorporate the new information from the 2020 accretion outburst we set up a first scenario, "A", with 20 parameters: M , R, T 0 , 5 zones with different Q sh can take any value between 8.2 and 10.2 (in g cm −3 ). Under these premises, we ran our MCMC (see Appendix B for details) and generated about 2 millions samples. We present in Figure 4, continuous lines, the posterior distribution of the eight parameters Q sh . There is a clear dichotomy between outburst 1 on one hand and the next three ones 6 Notice that our spectral fits to deduce T eff were performed assuming M = 1.4M and R = 10 km: for full consistency with our MCMC runs these should be performed for a range of M and R. These, we expect, would result in shifts of T eff of the order of 10% and should not have major impact on our conclusions.
sh show no upper limit because, as found by Deibel et al. (2015), when it is much above 10 MeV the outer crust becomes so hot that most of the shallow heating energy is lost to neutrinos and increasing Q sh has no further effect on the temperature profile so that an arbitrary cut-off has to be introduced (we choose 20 MeV). However, the distributions of both Q (j) sh and log 10 ρ (j) sh in outbursts 2, 3, and 4 point to the possibility that shallow heating had the same properties in these three outbursts. To evaluate this possibility we consider a second scenario, "B", in which Q for j = 2, 3, 4. The resulting posterior distributions, from a second MCMC run (see Appendix B for details) generating also about 2 millios samples, are plotted in Figure 4 as dotted lines. To compare these two scenarios we plot in Figure 5 the resulting χ 2 distributions: it is clear that there is no really significant difference between them and considering that scenario "B" has 16 parameters versus 20 in scenario "A", it seems there is no gain in considering that the shallow heating had different properties in the last three outbursts.
(Scenario "C" is described in the next section.) It is important to notice that in scenario "B" the distribution of Q (234) sh is essentially determined by the third outburst in which both scenarios result in the same posterior distribution. Hopefully further observations of the relaxation of MAXI J0556-332 in the near future will constrain more strongly the properties of shallow heating in the fourth outburst and strengthen our present conclusion that data from outbursts 2, 3, and 4, are compatible with having identical Q sh and ρ sh .

A GIGANTIC THERMONUCLEAR EXPLOSION IN THE FIRST ACCRETION OUTBURST
Having shown that the properties of shallow heating were likely identical in outburst 2, 3, and 4 we now postulate that this was also the case during the first 2011-2012 outburst and that the extremely high temperature of the neutron star at the end of that outburst is due to another phenomenon. We explore here the possibility that this event may have been a gigantic thermonuclear explosion, triggered deep in the outer crust. We, thus, present a third scenario, "C", with a single parametrization of the shallow heating, Q sh and ρ sh , operating in all four accretion outbursts and add a sudden heat injection with parameters X hb , t hb , and ρ hb , as presented in Section 3, with t hb adjustable to any time during the first outburst, whose duration was 1.31 year, and ρ hb allowed to take values between 10 8 up to 10 12 g cm −3 and X hb in the range 0 − 0.05.
To explore scenario "C" we ran again our MCMC (see Appendix B for details) and in Figure 5 the resulting distribution of the χ 2 of this scenario is also exhibited. There is no significant difference between the χ 2 of the three scenarios, "A" and "C" have slightly larger distri-bution below 20 but this could simply reflect the fact they have more parameters than "B" (20 and 17 versus 16, respectively). This simple comparison indicates that scenario C appears statistically at least as good as the first two.
In Figure 6 we illustrate in detail one of the good models found in this manner. The two upper panels display the overall evolution of both the daily mass accretion rateṀ and the star's red-shifted effective temperature T ∞ eff while the four central ones exhibit the excellent fit to the 19 data points that we obtain for this more than a decade long evolution. In the lower left panel of this figure one sees the constant rise in the temperature at the bottom of the exploding layer, at density ρ hb , which could be consistent with this temperature reaching a critical value at which some nuclear burning becomes unstable resulting in a thermonuclear runaway. We attempt to identify the fuel of the explosion in the next section. In the lower right panel we present the temperature profile in the whole star at the time just before the explosion, compared with a model where shallow heating was not implemented, showing that the occurrence of shallow heating was crucial in setting this profile in the outer crust. In Figure 7 we detail the evolution of the internal temperature, before (left panel) and after (right panel) the hyperburst starting from the beginning of the first accretion outburst, at time t start . The profile at 48 days corresponds to the highest temperature in the region where the shallow heating is operating and is occurring when the mass accretion rate is at its maximum. In the subsequent profiles one notices the decreasing temperature of this low density region due to both the decrease of the mass accretion rate and the flow of heat inward, the latter resulting in a rising temperature in the density range 10 10 to 10 12 g cm −3 : this explains the constant temperature rise at the point where the explosion will be triggered, at ρ hb = 10 11 g cm −3 , as seen in the lower left panel of Figure 6. Notice that in the low density layers, due to their low heat capacity and resulting short thermal timescales, the temperature follows the variations of the mass accretion and heating rates. At higher densities the heating has a more cumulative effect and the temperature keeps rising while the mass accretion rate decreases (a similar evolution can be clearly seen in the "Supporting Information" movie of Ootes et al. 2016). The right panel of Figure 7 illustrates the relaxation of the temperature after the heat injection from the hyperburst: due to a fixed energy injection per gram the initial profile peaks at the highest density ρ hb . The initial rapid drop in temperature is due to neutrino losses (Deibel et al. 2015) and during the relaxation under accretion, i.e., at t before t 0 , the Temperature profile at explosion Figure 6. Evolution of a 1.6 Modot model, with a radius of 11.2 km and an initial red-shifted core temperature T0 = 1.2 × 10 7 K. Its shallow heating parameters are Q sh = 0.390 MeV distributed in the density range [1.15 − 5.7] × 10 9 g cm −3 . The two upper panels show the whole temporal evolution, spanning more than ten years and four accretion outbursts, of the daily mass accretion rate,Ṁ , and red-shifted effective temperature, T ∞ eff , while the four central panels show details of the relaxation after each accretion outburst. (In parentheses are the initial and final dates of the outburst.) With 19 measurements of T ∞ eff this model has a χ 2 of 24.6. In all these panels the grey shaded areas delineate the periods of accretion. In the "Post-outburst 4" panel the two dotted lines show the 3σ range of predictions from the whole set of our scenario "C" models. The lower left panel displays the evolution of the local temperature in the critical layer at density ρ hb = 10 11 g cm −3 during the first accretion outburst where the hyperburst occurred at time of 56030 MJD., i.e., three weeks before the end of the outburst. The lower right panel displays the local temperature profile in the whole star, continuous line, just before the time when the hyperburst will be triggered, where the red dot signals the position of the critical layer at ρ hb , while the dashed curve shows the temperature profile at the same time in the absence of shallow heating. The blue shaded density range shows the region where the shallow heating is assumed to take place. The dotted curves show the melting temperature of the main matter (upper dotted curve) and of 28 Ne bubbles in pressure equilibrium with the main matter (lower dotted curve).  Figure 7. Evolution of the neutron star internal temperature profile for the same model as in Figure 6. Left panel: profiles starting at the beginning of the first accretion outburst, tstart, and five posterior times till just before the hyperburst explosion, t − hb . Right panel: profiles starting just before, t − hb , and after, t + hb , the hyperburst explosion, with three posterior times till the end of the accretion outburst, t0, and two later times during the quiescent phase. In both panels the blue shaded region marks the density range where the shallow heating is operating (during accretion). (Numerically, the hyperburst energy is smoothly injected during 200 seconds, the actual time interval from t − hb to t + hb . Moreover, no energy is injected close to the outer boundary at 10 8 g cm −3 but this anomaly is corrected by diffusion in less than a day as seen from the t hb + 1 day profile.) 35.5 36.0 36.5 log 10 L γ max (erg s −1 ) (f) Figure 8. Histograms of the distributions of some results from our scenario C MCMC: (a) total hyperburst energy, E hb , (b) total heating (shallow + deep crust), E h , during the four outbursts, total energy lost to (c) neutrinos, Eν , and (d) photons, Eγ, from the beginning of the first outburst till the end of the fourth one. Panel (e) shows the ratio of energy lost through photon emission to the total injected energy E hb + E h . The last panel, (f), shows the histogram of the maximum luminosity coming out of the stellar interior during the hyperburst. peak moves back toward the shallow heating region. After the end of the accretion the absence of heating and the heat flow toward the surface result in the maximum temperature peak moving toward higher densities. As a result of the inward heat flow, a large part of the energy released by the explosion is being stored at higher densities and will result in a very long cooling time: one sees in the "Cooling curve" panel of Figure 6 that even after the perturbation from the third outburst the cooling trajectory is a continuation of the trajectory from the first outburst (Parikh et al. 2017a).
The global energetics of scenario "C" resulting from our MCMC are presented in Figure 8. We find that a typical hyperburst energy E hb is of the order of 10 44 ergs, which is about 2 orders of magnitude larger than the one of a superburst (Cumming et al. 2006) and is comparable to the energy output of a magnetar Giant Flare (Kaspi & Beloborodov 2017). However, the total energy injected into the star from both the shallow and the deep crustal heating, E h , is larger than E hb by a factor of a few and the energy lost through photon emission, E γ , is of comparable magnitude: in most cases, however, E γ is somewhat smaller than the total heating energy E hb + E h . The energy not lost to photons is stored into the core and neutrino losses will contribute to the global energy balance in the long term (Brown et al. 1998;Colpi et al. 2001), but a detailed study of this issue is beyond the scope of the present work. Finally, the last panel of Figure 8 shows the distribution of the peak luminosity L γ max coming from the stellar interior during the hyperburst that is radiated as photons from the surface: unfortunately this luminosity L γ max is always significantly smaller than the X-ray luminosity inferred from the observed flux, ∼ 10 38 erg s −1 for a distance of 43.6 kpcs as we assume here, and thus unlikely to be noticeable in the data. However, since we do not model the explosion at low densities, our outer boundary being at density ρ b = 10 8 g cm −3 , we cannot exclude that a hyperburst may trigger a lower density X-ray burst, in a manner similar to the superburst precursors (see, e.g., Galloway & Keek 2021), which could be detectable.
We also show in Figure 9 the posterior distribution of the strength, Q sh , and lower density, ρ sh , of the shallow heating in our scenario "C". The Q sh distribution peaks at 0.6 MeV, which is slightly higher that in scenario "B" where the peak was at 0.4 MeV. As in scenario "B" the distribution of ρ sh is leaning on its lower permitted value: exploring this peak at lower densities would, however, imply lowering the outer boundary density of our models in NSCool and require a significant extension of 0.0 0.5 1.0 Q sh (MeV) 8 9 10 log 10 ρ sh (g cm −3 ) Figure 9. Histograms of the distribution of the shallow heating strength, Q sh , and lower density, ρ sh in scenario C.
the code's numerics that will be implemented in a future paper 7 .
We emphasize that the MCMC process was driven to fit the data using as parameters for the hyperburst only its occurrence time, t hb , depth, ρ hb , and total energy through X hb . We made no further assumption about the nature of the explosion, its fuel and the temperature at which it was triggered. To characterize the physical conditions of the hyperburst trigger, we display in Figure 10 a summary of the most important results of this scenario "C". Beside the posterior distribution of 5 of our most relevant MCMC parameters we also report the posterior distribution of two quantities that are output of our calculation: 1) the temperature of the critical layer just before ignition, T hb , and 2) the cooling sensitivity (Eqs. 4 and 3) that will be of interest in the next section. The first noticeable result is a slightly bimodal distribution in terms of mass and radius: the dominant peak is located at ∼ (2M , 10.5 km), and the second smaller one around ∼ (1.5M , 12.5 km). These peaks also exhibit clear differences in critical density and temperature: for the first quantity, there is a threshold at 10 11 g cm −3 above (below) which the dominant (lower) peak can be found, while for the latter quantity their respective temperatures are ∼ 10 8.6 ≈ 4 × 10 8 K for the dominant and ∼ 10 8.4 ≈ 2.5 × 10 8 K for the lower one. Of relevance also is that the dominant peak favors a trigger time t hb very close to the end of the accretion outburst while the second peak favors t hb 1.25 year, i.e., about 25 days before the end of the outburst. Notice, finally, that the column density y hb corresponding to ρ hb has a narrow, symmetric distribution centered at 10 15 g cm −2 . log 10 T hb (K) log 10 dε cool /dT 2 0 3 0 4 0 5 0 6 0 χ 2 Figure 10. Corner plot of 1D and 2D histograms of the posterior distributions of five MCMC parameters, mass and radius of the model and the three hyperburst control parameters, X hb , t hb , and ρ hb , as well as its equivalent column density y hb , and two deduced quantities, the critical temperature, T hb , and the cooling sensitivity, log 10 dε cool /dT (see Eqs. 4 and 3 with ε cool measured in erg g −1 s −1 ), calculated at the critical density at time t hb , as well as the χ 2 of the model fit to the 19 data points. Neutron drip thresholds Figure 11. Evolution of the charge number Z, under the effect of double electron captures, of the six lightest α-nuclei, 12 C, 16 O, 20 Ne, 24 Mg, 28 Si, and 32 S, from their production site, at densities below 10 9 g cm −3 , when pushed to increasing densities. Continuous lines allow forbidden electron captures while dashed ones only consider allowed transitions. Vertical dotted lines show the neutron drip density thresholds for each isobaric sequence. The grey shaded region corresponds to the predicted location of the origin of the hyperburst and in this density range excluding or including forbidden transitions makes no difference. Details of the calculation are provided in Appendix C.

THE FUEL OF THE HYPERBURST
In the classical short H/He X-ray bursts, as well as the long C superbursts, matter in the stellar envelope is pushed by accretion toward higher densities and temperatures till it reaches a critical point where the nuclear burning becomes unstable and a thermonuclear explosion ensues. The time scale for this compression is hours to days in the case of the short bursts and months to years for the superbursts. In our scenario, an explosion is triggered at densities ∼ 10 11 g cm −3 , corresponding to column densities of the order of 10 15 g cm −2 : under an Eddington accretion rate (ṁ Edd ≈ 10 5 g cm −2 s −1 ) at such densities matter is barely progressing by one millimeter per day and the time scale to reach this point is centuries to millennia. Moreover, matter is being pushed toward lower temperatures because of the inverted temperature gradient (right lower panel of Figure 6). However, we are not in a steady state situation but rather in a year and a half long strong outburst during which the critical layer is still warming up (left lower panel of Figure 6 and left panel of Figure 7): the physical origin of the explosion appears thus to be the result of the temperature at the critical layer to be rising with time rather than this matter being pushed toward higher densities.
A standard criterion in X-ray burst modeling to identify the location of the ignition layer (Fujimoto et al. 1981;Bildsten 1998) is to compare the temperature dependence of the nuclear burning rate, ε nucl , with the cooling rate where F is the heat flux. As long as the inequality is satisfied the burning is stable. Histograms of the values of log 10 dε cool /dT , calculated with the temperature profile of each model just before the time t hb at ρ hb , were presented in Figure 10: we found that triggering of the hyperburst occurred at values log 10 dε cool /dT ∼ 2 (with ε cool in erg g −1 s −1 ). The candidate nucleus to trigger the hyperburst must have been produced by nuclear burning, either stable or explosive, at lower densities and then pushed to high densities by accretion. On its journey this nucleus will undergo electron captures as described, e.g., by Haensel & Zdunik (1990). The most abundant light nuclei produced by the nuclear burning in the envelope are αnuclei and in Figure 11 we show how the six lightest ones (excluding 8 Be which is unstable) evolve under double electron captures when pushed to increasing densities: these are our pre-candidates for the triggering of the hyperburst.
One interesting feature of the physical conditions where the hyperburst is triggered, which may affect the conclusions about the triggering species, is that the crust main matter, in the temperature and density ranges where the explosion starts, is in a solid (likely crystalized) state: the lower right panel of Figure 6 displays the melting temperature curve of the main matter which shows that on this temperature profile solidification occurs close to 10 10 g cm −3 . However, it has been predicted (Horowitz et al. 2007) that when the high Z nuclei solidify there is a phase separation with the low Z ones that remain in the liquid phase: the emerging state could thus be one of bubbles of liquid O and/or Ne embedded in a solid medium. We also display in the lower right panel of Figure 6 the melting temperature for 28 Ne bubbles in these conditions which shows that they remain in a liquid state, given this temperature profile, up to densities well above 10 11 g cm −3 . Oxygen bubbles  Figure 12. Upper panels: contour lines where log 10 dεnuc/dT = 2 (ε nucl being in erg g −1 s −1 and T in K) for the fusion of 16 C+ 16 C (labelled "Z = 6"), 20 O+ 20 O (labelled "Z = 8"), and 28 Ne+ 28 Ne (labelled "Z = 10"), and for each case with seven mass fractions X = 10 −n with n = 0, 1, 2, . . . , 6, as dash-n-dot lines. Lower panels: contour lines of life-time τ = 10 3 yrs of the 16 C, 20 O, and 28 Ne nuclei under the same fusion reactions as in the upper panels. For the treatment of screening, these nuclei are assumed immersed in a background of 56 Ca, with mass fraction 1 − X, which is the nuclear species predicted to be present at ρ ∼ 10 11 g cm −3 in the model of Haensel & Zdunik (2008). In the left panels we apply the minimum reaction rates while in the right panel the maximum ones are used. Details of the calculations are presented in Appendix D. In all panels the background contours reproduce the 1σ, 2σ, and 3σ, confidence ranges of the 2D distribution of log 10 ρ hb − log 10 T hb already displayed in Figure 10.
have even lower melting temperatures than neon ones because of their lower Z.
In Figure 12 we present the critical conditions for the burning of α-nuclei with Z = 6, 8, and 10 (details are presented in Appendix D) in the ρ − T plane where the 1σ, 2σ, and 3σ, critical confidence regions for the values of (ρ hb , T hb ) from Figure 10 are shown in the background. In the upper panels we plot contour lines of log 10 dε nucl /dT = 2, since log 10 dε cool /dT = 2 is a typical value found, see Figure 10, in our scenario "C" MCMC. In the lower panels we plot contour lines of the lifetimes τ = 1, 000 yrs of the same nuclei, since 1,000 yrs is a typical time for matter to be pushed to these densities by accretion. Even if the medium average mass fraction X of the triggering nucleus is small, its mass fraction in the bubbles where it is concentrated during medium crystalization, if concentration did occur, could be close to unity. To encompass this range of possibilities all critical contour lines in Figure 12 consider a range of X from 10 −6 , no bubble formation, up to unity, perfect bubble formation. Depending on the actual mass fraction of the possible exploding nuclei we have a variation of a factor of a few in the predicted density where the explosion started. At these densities and relatively low temperatures, uncertainties on the burning rates are largely dominated by uncertainties in the screening and we apply either the estimated minimum rates, in the left panels, or the maximum ones, in the right panels. Notice that the burning rates depend strongly on the charge number, Z, of the nuclei and have much less sensitivity to their mass number A so that considering different isotopes has little impact on this part of our inquiry. From this Figure 12 we see that, clearly, carbon is not the trigger nucleus of our hyperburst (it is most likely the trigger nucleus of superbursts acting at lower densities ∼ 10 9 g cm −3 , Bildsten 2001 andCumming et al. 2006): it is only at very low concentrations, X 10 −4 , that carbon burning can be unstable within the critical region of the ρ−T plane but, unless something beyond our simple calculations happened, it would have been consumed by stable burning before reaching such high densities as seen from the lower panels of Figure 12. Oxygen appears to be the best candidate but neon cannot be excluded particularly if actual reaction rates are close to the maximum value we apply and the hyperburst critical density was on the high values found by the MCMC. Oxygen could be 20 O coming from 20 Ne after a pair of electron captures or 24 O coming from 24 Mg after two pairs of electron captures, as pictured in Figure 11. Similarly, neon could be 28 Ne coming from 28 Si or 32 Ne coming from 32 S. The lower panels of Figure 12 show that both oxygen and neon nuclei need to be pushed to densities higher than their explosion densities before they are exhausted and, thus, exhaustion by stable burning is not an issue.

DISCUSSION
Having now data for a fourth outburst in MAXI J0556-332, our first finding is that shallow heating during outbursts 2, 3, and 4 can be described by identical parameters, Q sh and ρ sh (see Eq. 2 and the description following it). The thermal evolution of this neutron star when exiting its first outburst implies, however, a very different physical condition, as is easily intuited by a look at Figure 1. Such condition had previously been described as an extremely strong continuous shallow heating, with Q sh ∼ 10-20 MeV compared to ∼ 1 MeV in outbursts 2 and 3 (Deibel et al. 2015;Parikh et al. 2017a).
Our second finding is that the high temperature of the MAXI J0556-332 neutron star when exiting outburst 1, as well as its subsequent years long thermal evolution, can be very well modeled assuming a very strong and sudden energy release, a "hyperburst", instead of a continuous energy release as in the shallow heating scenario. From our extensive MCMC simulation, scenario "C", we find that this event would have occurred during the last 3 weeks of the accretion outburst and released an energy of the order of 10 44 ergs. In this scenario the energy was deposited in a layer extending from the lowest density up to ∼ 10 11 g cm −3 : global energetics as well as other characteristics are presented in Figs. 8 and 10. Shallow heating was also assumed to occur during the first outburst, and described by the same parameters as in the next three outbursts (see Figure 9 for the parameter posterior distributions).
To identify candidates for the triggering of the hyperburst we simply applied the standard one-zone model criterion of Equation (4). Values of dε cool /dT were obtained directly from perturbation of the T -profile in each one the models of our MCMC and we found that log 10 dε cool /dT ∼ 2, see Figure 10, with ε cool measured in erg g −1 s −1 . Values of log 10 dε nucl /dT = 2 were calculated and displayed in the two upper panels of Figure 12: considering the enormous uncertainties on the fusion rate due to screening in, and near to, the pycnonuclear regime, the left upper panel employed minimum rates and the right upper panel maximum ones. Lower panels of the same figure show "exhaustion lines", i.e., contour lines where the burning time scale is 1,000 yrs, implying that the described nucleus would be depleted by stable burning when reaching the corresponding density since 1,000 yrs is roughly the time it takes, under MAXI J0556-332's average mass accretion rate, for mat-ter to reach 10 11 g cm −3 . The conclusion of this is that C could be the triggering nucleus, but it is unlikely to be able to reach the required density before being depleted, leaving O and Ne as the best candidates. Isotopes of O and Ne have slightly different burning rates, but differences are small enough that the results of Figure 12 are practically only dependent on the charge Z of the nucleus and not on their mass number A. Once produced at low densities, nuclei are compressed by accretion and undergo double electron capture reactions that gradually reduce their charge Z, see Figure 11, but without changing the mass A: the triggering oxygen is likely 20 O or even 24 O, descendants of 20 Ne and 24 Mg, while the triggering neon would be 28 Ne or 32 Ne, descendants of 28 Si and 32 S.
Fusion of low Z nuclei liberates about 1 MeV per nucleon (see Table 1), i.e. about 10 18 erg g −1 , and we parametrized the energy injected during the explosion as X hb ×10 18 erg g −1 so that X hb approximately reflects the average mass fraction, X , of the exploding nuclei. Our MCMC found a peak at X hb ∼ 0.02, exhibited in Figure 10, while the minimum value found was 0.011. Thus, a small mass fraction of a few percents of low Z nuclei is needed to provide the required energy. While accretion is ongoing, in the upper layers of the accreted material H burns into He through the hot CNO cycle and He into C through the triple-α reaction. At high accretion rates, as is the case in MAXI J0556-332, breakout reactions from the CNO cycle add the rp-process (Wallace & Woosley 1981) to H burning and lead to the production of high Z nuclei, possibly up to the SnSbTe cycle at Z = 50-52 and A = 103-107 (Schatz et al. 2001). In the absence of X-ray bursts, the models of Schatz et al. (1999) with continuous burning in the ocean at accretion rates ∼Ṁ Edd produced A = 12 nuclei at a mass fraction ∼ 4% and A = 24 at ∼ 2−3% with other low A nuclei in much smaller abundances and heavy nuclei up to A ∼ 80 through the rp-process. This is likely representative of the burning in most of the observed accretion phases of MAXI J0556-332 that were close toṀ Edd most of the time (see, e.g., upper panel of Figure 6) as explosive burning is quenched at high mass accretion rates (e.g. Bildsten 1998;Galloway & Keek 2021). However, during low mass accretion moments, as, e.g., toward the end of an accretion outburst, bursting behavior may appear as was observed in the case of the end of the outburst of XTE J1701-462 whenṀ was down to about 10%Ṁ Edd (Lin et al. 2009). We notice that MAXI J0556-332's first outburst, in 2011-2012, was almost a carbon copy of XTE J1701-462's 2006-2007 outburst but, nevertheless, no X-ray bursts were detected while we may have found a small X-ray burst toward the end of the 2020 outburst as presented in § 2.1. Within this context of bursting behaviour, we can consider the ashes of the three cases used by Lau et al. (2018): in the extreme rpprocess X-ray burst case (based on Schatz et al. 2001) one finds negligible amounts of C and O left but about 1% of A = 20 nuclei, while in their KEPLER X-ray burst case (based on Cyburt et al. 2016) one has to go up to A = 28 and 32 to find nuclei with a significant mass fraction, about 5% and 10%, respectively: in both cases we have enough seed nuclei which, after multiple double electron captures when pushed to high densities, could be the trigger nuclei for a hyperburst. The third case we contemplate is the occurrence of a superburst (made possible, e.g., by C ashes from long term continuous burning as seen above): the resulting ashes (based on Keek & Heger 2011) contain negligible amount of C, O, and Ne, but almost 1% of A = 28 nuclei which could lead to a seed nucleus for hyperburst triggering after electron captures. Since the hyperburst is triggered at densities ∼ 10 11 g cm −3 , a region whose chemical composition is the result of many centuries of accretion, of which we only have a decade long snapshot, all of the above burning possibilities may be partial realities and show that having ∼ 1 − 3% of low Z nuclei present at such densities is a realistic possibility.
We see from Figure 10 that values of the column density at the explosion point, y hb , range from 5 × 10 14 up to 2 × 10 15 g cm −2 . Under an Eddington rate it takes then approximately from 150 to 600 years, respectively, of continuous accretion for matter to be pushed to the explosion point. From MAXI J0556-332's inferred bolometric flux and Equation (1) we deduce a total accreted mass, in the four observed outbursts, of almost 10 26 g resulting in a 10 years average Ṁ 2.8 × 10 17 g s −1 4.4 × 10 −9 M yr −1 , i.e., about 30% of the Eddington rate. Since MAXI J0556-332 was never seen previously to its 2011-2012 outburst the long term actual fraction of time it spends in quiescence is possibly much larger than indicated by the last decade of observations and its long term Ṁ smaller that estimated above. Hyperbursts in MAXI J0556-332 are thus likely a once in a few millennia events and the probability to have witnessed one in 10 years since its discovery, or in 25 years of having continuous all sky monitoring in X-rays, is of the order of 1%: a small but not too small value. The seven known persistent Z sources (Hasinger & van der Klis 1989;Fridriksson et al. 2015), also estimated to be accreting close to, or above, the Eddington limit most of the time, are likely to experience hyperbursts more frequently, once every few centuries, but unfortunately these are not detectable under continuous accretion due to the small increase in surface luminosity during the ex-plosion, see panel (f) in Figure 8. From what we found in the present work, a hyperburst can only be recognized if accretion stops shortly afterwards and we detect an anomalously hot neutron star. It could be considered too much of a coincidence that the hyperburst need to have occurred only a few days/weeks before the end of the outburst. Notice, however, that even if it had happened several months earlier, the star would still have exited the outburst with a very high effective temperature and allowed the identification of the occurrence of the hyperburst: Figure 1 shows that it took more than 100 days after the end of the outburst, i.e., after the occurrence of the hyperburst, for the crust to cool down and T ∞ eff to drop below 250 eVs. Had the hyperburst occurred 3-4 months before the end of the outburst, the star would have entered quiescence with a temperature in excess of 250 eVs which would already have singled out MAXI J0556-332 as an anomalously hot star. At least in the case of MAXI J0556-332 Nature appears to have been cooperating! We now have a handful of binary systems under transient accretion in which post outburst cooling data are available. Most of them have much lower mass accretion rates than MAXI J0556-332 which implies that their crust experience lower temperatures and that their matter would need to be compressed to still higher densities for a hyperburst to be triggered. Hyperbursts in these system are, thus, even rarer events than in high accretors and occurrence, and identification, of one in any of these systems appears to be, unfortunately, highly unlikely, unless many more such systems are discovered. XTE J1701-462 is the exception within this sample as it exhibited an outburst very similar to MAXI J0556-332 first observed outburst, while its neutron star exited it with a much lower surface temperature of the order of 150 eV: it may need many more such outbursts till a hyperburst is triggered.

CONCLUSIONS
We have shown that the highly anomalous thermal evolution of the MAXI J0556-332 neutron star when exiting its 2011-2012 accretion outburst and during the following years can be very well modeled as due to the occurrence of a hyperburst, i.e., explosive nuclear burning of some neutron rich isotope of O or Ne in a region of density ∼ 10 11 g cm −3 . Such a deep explosion, depositing ∼ 10 44 ergs, results in a hot outer crust whose cooling takes several years and provides an excellent fit to the data. Published models of nuclear, either stable or explosive, burning during accretion show that the subsistence of a small amount, a few percents in mass fraction, of the needed low Z nuclei to these densities is likely. However, it takes many centuries of accretion to accumulate this material in transient systems and hyperbursts as postulated here are, thus, very rare events and we are unlikely, unfortunately, to witness another one.
Our complementary results that modeling of the cooling of the neutron star after the subsequent three accretion outbursts, in 2012, 2015-2016, and 2020, can be realized employing a single parametrization of the still enigmatic shallow heating, lead to a prediction about the future cooling of the neutron star which will be monitored in the coming years. The example presented in Figure 6 shows a future smooth evolution, see the "Postoutburst 4" panel, for at least 1,000 days, and the 3σ range deduced from our scenario "C" MCMC run (containing more than 4 millions cooling curves) shows that this simple behavior is a strong prediction of the scenario.
In the global modeling of 12 years of evolution of MAXI J0556-332 in our scenario "C" with the hyperburst occurring at the end of the first outburst, parameters describing the shallow heating during this outburst were also assumed to have identical values as in the next three outbursts. We can, hence, describe the whole evolution of the MAXI J0556-332 neutron star during a period of 12 years through 4 very different accretion outbursts with a single consistent parametrization of the shallow heating. This result is in agreement with the study of the MXB 1659-29 neutron star that could describe its evolution through two accretion outbursts, spanning a period of almost two decades (1999 till 2018), also with a single shallow heating parametrization (Parikh et al. 2019). These results are, however, in contrast with the modeling of the evolution of the Aquila X-1 neutron star, which exhibits frequent but short accretion outbursts, that required a variation of the shallow heating parameters between different outbursts (Degenaar et al. 2019).
We would like to thank the NICER and Swift teams for approving and executing the observations made during and after the 2020 outburst. DP and MB acknowledge financial support by the Mexican Consejo Nacional de Ciencia y Tecnología with a CB-2014-1 grant #240512 and the Universidad Nacional Autónoma de México through an UNAM-PAPIIT grant #109520. MB also acknowledges support from a postdoctoral fellowship from UNAM-DGAPA and support from a grant of the Ministry of Research, Innovation and Digitization, CNCS/CCCDI -UEFISCDI, Project No. PN-III-P4-ID-PCE-2020-0293, within PNCDI III. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We take the equation of state and chemical composition of the crust from Haensel & Zdunik (2008) for the model assuming 56 Fe at low densities, and fix the crust-core transition density at ρ cc ≈ 1.5 × 10 14 g cm −3 . Reasonable variation of this transition density is known to have very little effect on the crust relaxation modeling (Lalit et al. 2019).
The specific heat is obtained by adding the contribution of the degenerate gas of electrons, the degenerate gas of dripped neutrons in the inner crust and the nuclei, including the Coulomb interaction contribution in the liquid phase from Slattery et al. (1982) and in the solid phase from Baiko et al. (2001) with a liquid-solid phase transition occurring at a Coulomb coupling parameter Γ = 180. The only strong interaction modifications to the neutron specific heat we include is from the effect of pairing for which we follow Levenfish & Yakovlev (1994) while the pairing phase transition critical temperature is taken from Schwenk et al. (2003).
For the thermal conductivity, dominated by electrons, we follow Potekhin et al. (1999) and Gnedin et al. (2001) for electron-phonon scattering in the solid phase to which we add electron-impurity scattering following the simple treatment of Yakovlev & Urpin (1980). When ions are in the liquid phase we apply the results of Yakovlev & Urpin (1980). We neglect the very small contribution of electron-electron scattering (Shternin & Yakovlev 2006). We also neglect heat transport from neutrons in the inner crust that is a very minor contribution (Schmitt & Shternin 2018).
In a very hot crust as in MAXI J0556-332, neutrino emission from plasmon decay can become a significant sink of energy and we employ the result of Itoh et al. (1996). We also include neutrino emission from pair annihilation from the same authors, a process that may contribute in cases of extremely hot outer crusts. For the small contribution of electron-ion bremsstrahlung we follow Kaminker et al. (1999) and we also include neutrino emission for the Cooper breaking and formation process in 1 S 0 neutron superfluid phase transition in the inner crust following the treatment described in Page et al. (2009). We do not include neutrino emission from the Urca cycles (Schatz et al. 2014) which may have a significant effect (Deibel et al. 2015) but would force us to introduce more parameters in our MCMC runs to explore the variations in chemical composition that control it.
Our handling of the heating, both from deep crustal heating and shallow heating, and the thermonuclear explosion was described in § 3. We neglect the effect of nuclear burning in the envelope after the accretion phase, in the form of diffusive nuclear burning, as it has only a small effect for a few days (Wijngaarden et al. 2020).
We do not focus on the physics of the stellar core in the present work and we simply follow the minimal scenario as described in Page et al. (2004) with no pairing taken into account.
B. SOME DETAILS OF OUR MONTE CARLO RUNS Our priors for the distributions of the MCMC parameters are as follows. For the mass M and radius R of the neutron stars we limit the range of the former from 1.2 to 2.4 M and of the latter from 8 to 16 km with a joint probability distribution in these ranges that is displayed in Figure 13. This distribution encompasses the deduced posteriors from the two classes of models, "PP" and "CS", of Raaijmakers et al. (2021). The initial, red-shifted, core temperature T 0 is taken as uniformly distributed in the range (0.1 − 1.5) × 10 8 K. About the column density of light elements in the envelope, y L , measured in g cm −2 , we assume a uniform distribution of log 10 y L in the range 6 to 10, and choose it independently at the beginning of each accretion outburst, thus having four values y L (k), k = 1, . . . , 4. The resulting T 0 − T ∞ eff relationship is illustrated in the right panel of Figure 3 in Degenaar et al. (2017). Our next set of parameters are the impurity parameters, Q imp , and we divide the crust in five density ranges, in g cm −3 : [10 8 , 10 11 ], [10 11 , 10 12 ], [10 12 , 10 13 ], [10 13 , 10 14 ], [10 14 , ρ cc ], with a value Q (i) imp , i = 1, · · · , 5, in each range that is uniformly distributed between 0 and 100. The parametrization of the heating, either shallow or from the hyperburst has been described in the main text in Section 3.
When varying M and R we recalculate the structure of the crust by integrating the Tolman-Oppenheimer-Volkoff equation of hydrostatic equilibrium from the outer boundary at ρ = ρ b = 10 8 g cm −3 and radius r = R inward till we reach the crust core density ρ cc , thus giving us the core's mass and radius, m c and r c . This procedure gives us a self-consistent structure of the crust. In this work, we are not interested in the response of the core but still need to define its density and chemical composition profile to employ NSCool. For this purpose we start with a core structure calculated using the APR EOS (Akmal et al. 1998) and a stellar mass of 1.4M (that has a core mass m  (Christen 2007;Goodman & Weare 2010). For the present work we found the optimal configuration was having 5 tempered chains, with "temperatures" T = 1, 2, 5, 10, and 100, each having 100 walkers. The basic chain applied stretches and the other chains simple walks. After initial burn-in, the chains of scenarios "A" and "B" had more than 2 millions points while the one of scenario "C" had more than four millions points. We calculated the integrated autocorrelation lengths (Sokal 1997) τ of each parameter of each walker: typical values are 50 − 100 and the longest ones do not exceed 200 in all three scenarios. We thus have, in each scenario, more than a hundred effectively independent samples from each walker.
For completeness we display in Figure 14 the posterior distributions of the remaining 10 parameters of our scenario "C" MCMC that were not displayed in Figure 9 and 10 for not being crucial to the purpose of this work. Interesting to notice are the preference for cold cores, T 0 , and the contrast in the distribution of Q imp in region 1 and 2, i.e., essentially below neutron drip, versus regions 3, 4, 5, i.e., above neutron drip. One expects a hyperburst to significantly reduce the impurity content as all low Z nuclei are burned into iron peak nuclei while the MCMC prefers high impurity in the region (1), i.e., at densities below 10 11 g cm −3 , precisely where the hyperbust happened (and in region (2) whose matter should have also been processed by previous hyperbursts). However, low Q imp are not excluded in regions (1) and (2), only disfavored, and low Q imp values actually just favor a thick crust, i.e., not too high neutron star masses and not too small radii. There is, thus, nothing contradictory, nor conclusive, in these impurity content posteriors. They would favor the secondary peak, low M and large R, in the M and R distributions discussed in Section 5, if one added the new prior that the hyperburst would result in low Q imp , at least in region 1. (For this reason we choose our illustrative model for Figure 6 and 7 as a 1.6 M star with an 11.2 km radius.)  Figure 14. Posterior distribution of the remaining parameters in our scenario "C" which are not displayed in Figure 9 and 10.

C. ELECTRON CAPTURES ON α-NUCLEI
For studying electron captures, we have followed the scheme of Sato (1979) and Haensel & Zdunik (1990). We choose an initial nucleus (Z, N ), N = A − Z being the neutron number of the nucleus, and follow its evolution as pressure is increased. We only consider α nuclei, i.e., Z is even and Z = N . When reaching a critical point where the electron chemical potential µ e is large enough a first electron capture occurs resulting in an odd-odd nucleus (Z − 1, N + 1) at which point, due to pairing, a second electron capture immediately happens leading to an even-even nucleus (Z − 2, N + 2). The criterion for the first electron capture is that W N (Z, N ) + W C (Z, n e ) + µ e (n e ) = W N (Z − 1, N + 1) + W C (Z − 1, n e ) where W N (Z, N ) = M N (Z, N )c 2 is the energy, M N being its mass, of the nucleus (Z, N ), W C (Z, n e ) the Coulomb energy of the nucleus' Wigner-Seitz cell immersed in a medium with electron density n e , and µ e is the electron chemical potential. For the Coulomb energy we apply the ion sphere value W C (Z, n e ) = −0.9Z 5/3 e 2 /a e where e is the elementary charge unit and a e ≡ (3/4πn e ) 1/3 (Shapiro & Teukolsky 1983). We deduce M N (Z, N ) from the atomic mass M atom (Z, N ) of the AME2020 table  as M N (Z, N ) = M atom (Z, N ) − Zm e + B el where m e is the electron mass and the electrons binding energy is approximated as B el = 14.4381 Z 2.39 + 1.55468 × 10 −6 Z 5.35 eV (Lunney et al. 2003). We then keep increasing the pressure looking for further electron captures leading to (Z−4, N +4), a.s.o., until the neutron drip point is reached. The resulting evolutions are presented in Fig. 11.

D. NUCLEAR FUSION PROCESSES AND SCREENING
We consider nuclei of charge and mass numbers Z i and A i and mass M i undergoing a fusion reaction (A 1 , Z 1 ) + (A 2 , Z 2 ) → (A c , Z c ). For this we need to know the fusion cross section σ and calculate the fusion rate including a proper treatment of plasma screening. The cross section is commonly expressed as where S(E) is the "astrophysical S-factor", E being the center of mass energy, and η = Z 1 Z 2 e 2 / v is the Gamow parameter, with v = 2E/µ the relative velocity, µ = M 1 M 2 /(M 1 + M 2 ) being the reduced mass. For our purpose we need an extended set of fusion reactions involving neutron rich nuclei and we employ the results of Afanasjev et al. (2012) that provide a consistent scheme covering thousand of possible reactions. These authors calculated S(E) using the São Paulo potential with the barrier penetration formalism (SP-BP: Chamon et al. 2002;Gasques et al. 2005) and provide simple analytical fits to their results.
In the purely thermonuclear regime, i.e., at high temperatures where screening can be neglected, the fusion rate is given by the standard result (e.g., Kippenhahn et al. 2012) where n i is the number density of nucleus i, r B = 2 /(2µZ 1 Z 2 e 2 ), and S(E) is evaluated at the Gamow peak energy E pk = (πk B T /2) 2/3 E 1/3 a with E a = Z 1 Z 2 e 2 /r B . F = exp[−τ ] describes the Coulomb barrier penetration with τ = 3E pk /k B T and the prefactor P = 16τ 2 /(3 5/2 π). (The 1 + δ 12 term avoids double counting in the case nuclei 1 and 2 are identical.) Based on their previous work (Gasques et al. 2007), Afanasjev et al. (2012) state that typical uncertainties on their S(E) calculations are of the order of 2 to 4 for stable nuclei but can be about a factor 10 in the cases of neutron rich nuclei and even up to 100 at low energies. For example, in the case of the three typical fusion reactions of 12 C+ 12 C, 12 C+ 16 O, and 16 O+ 16 O we find that the thermonuclear rates from Eq. D3 using S(E) from Afanasjev et al. (2012) can be up to 4-5 times smaller that the classical values from Caughlan & Fowler (1988). For another comparison, including important cases of reactions with no experimental measurements, Umar et al. (2012) performed dynamical density-constrained time dependent Hartree-Fock (DC-TDHF) calculations of S(E) for several C and O reactions and compared them to the static SP-BP results: in the cases of 12 C+ 16 O, and 16 O+ 16 O they obtain excellent agreements with experimental values for the DC-TDHF results (while the SP-BP results are about 4 time smaller as seen above), in the cases of 12 C+ 24 O and 16 O+ 24 O they find that the SP-BP results are also smaller that the DC-TDHF ones, by about one order of magnitud, while in the case of the 24 O+ 24 O the situation reversed with SP-BP values being larger that DC-TDHF ones. (No experimental values are available for theses last three reactions.) As seen below, variations of even a factor 10 in S(E) have very little impact on our results. Notice that resonances are not included in any of these models that only provide an average S(E) but their effect is small as long as their width is much smaller than the width of the Gamow peak  .
The second ingredient needed to calculate realistic fusion rates in the dense medium of the neutron star crust is proper inclusion of the screening of the Coulomb repulsion. In our cases this effect increases the rates by many tens of orders of magnitude. We employ the treatment proposed by Gasques et al. (2005) and Yakovlev et al. (2006) which allows to cover the whole range from the T = 0 pycnonuclear regime up to the thermonuclear one by adjusting E pk and appropriately modifying the two functions P and F in Equation (D3). In the thermonuclear regime screening can be accurately incorporated while in the opposite pycnonuclear regime uncertainties are enormous, possibly up to 10 orders of magnitude. We will apply the three representative case of Yakovlev et al. (2006) that give minimal and maximal screening effects as well at their optimal model. From the rate per unit volume R one deduces the energy Table 1. Energy released per reaction, Q, and per nucleon, Q/A, both in MeV, for fusions between nuclei with same charge Z and mass numbers A1 and A2. Nuclear binding energies are taken from the 2020 Atomic Mass Evaluation Wang et al. 2021 generation rate, per unit mass, ε nucl = QR/ρ where Q is the energy released per reaction. Representative Q-values are listed in Table 1. Relevant screening factors, energy generation rates ε nucl and their temperature sensitivity dε nucl /dT are plotted in Figure 15. The flat portion of the ε nucl curve exhibit the pycnonuclear regime and comparison with Figure 12 show that the triggering of the hyperburst occurred in the regime of transition from purely pycnonuclear to thermonuclear.  Figure 15. Screening factors, energy generation rates, ε nucl , and their temperature sensitivity, dε nucl /dT , for fusions of several isotopes of C, O, and Ne, at chosen densities. We used the "optimal" reaction rates and in all cases assumed a mass fraction X = 10 −3 within a background of 56 Ca.