The minimum neutron star mass in neutrino-driven supernova explosions

Abstract

Supernova theory has struggled to explain the lightest known neutron star candidate with an accurate mass determination, the $1.174\,\mathrm{M}_{\odot}$ companion in the eccentric compact binary system J0453+1559. To improve the theoretical lower limit for neutron star birth masses, we perform 3D supernova simulations for five stellar models close to the minimum mass for iron core collapse. We obtain a record-low neutron star mass of $1.192\,\mathrm{M}_{\odot}$ and a substantial kick of $\mathord{\sim}100\,\mathrm{km}\,\mathrm{s}^{-1}$. Given residual uncertainties in stellar evolution, a neutron star origin for the $1.174\,\mathrm{M}_{\odot}$ object remains plausible.

Introduction.—Compact object masses are among the most critical astrophysical observables because of their implications for high-density nuclear physics, stellar evolution, and supernova explosion physics. In recent years, new records for the maximum neutron star mass [1, 2] together with complementary constraints on neutron star radii and tidal deformability from X-ray [3] and gravitational wave [4, 5, 6] observations, have provided important insights on the nuclear equation of state.

The lowest neutron star mass also has significant implications. In contrast to the upper mass limit, the lower limit is due the astrophysical formation path, and is not a limit inherent to neutron stars and their equation of state. The best known way to make the lowest possible neutron star masses is the gravitational collapse of the iron core of massive stars [7] or of the O-Ne-Mg core of Super Asymptotic Giant Branch (SAGB) stars [8, 9, 10]. The masses of young neutron stars will reflect the core sizes and compositions in massive stars, although the exact location of the “mass cut” will depend on details of the explosion dynamics. Hence, the lowest neutron star mass is a key parameter for testing our understanding of stellar evolution and nuclear physics [11] through advanced burning stages to core collapse [12]. In the alternative scenario of accretion-induced collapse [13, 14, 15], a higher electron fraction and rotational stabilization likely result in a higher mass of collapsing white dwarf and the resulting neutron star.

Precise measurements of neutron star masses in binary systems have pushed minimum gravitational neutron star mass below $1.2\,\mathrm{M}_{\odot}$. The double neutron star system J0453+1559 was found to contain a companion with only $1.174\,\mathrm{M}_{\odot}$ [16]. There is also a candidate neutron star of $1.2\,\mathrm{M}_{\odot}$ in J1807-2500B [17, 18], though the alternative interpretation of this object being a white dwarf cannot be excluded in this case. A white dwarf origin for the $1.174\,\mathrm{M}_{\odot}$ object in J0453+1559 has also been proposed [19]. There are further claims of even lower neutron star masses [20], but these are beset with large error bars

Such low neutron star masses present a challenge for current stellar evolution and supernova explosion models, which, to date, have failed to obtain such low masses. The record has long been held by simulations of electron-capture supernovae with a minimum baryonic neutron star mass of $1.36\,\mathrm{M}_{\odot}$ [21], which translates into a gravitational mass of about $1.24\,\mathrm{M}_{\odot}$, with only small uncertainties from the nuclear equation of state. In this narrow stellar evolution channel for stars around $\mathord{\sim}9\,\mathrm{M}_{\odot}$ in zero-age main sequence (ZAMS) mass, dynamical collapse is believed to occur after the formation of an O-Ne-Mg core in SAGB stars due to rapid electron capture (EC) on ${}^{20}\mathrm{Ne}$ and ${}^{24}\mathrm{Mg}$ [22, 8, 9, 23, 24, 10, 25] (although collapse is not certain [26]), different from “normal” core-collapse supernovae that undergo further hydrostatic burning stages and collapse only after iron core formation. For electron-capture supernovae, a uniform core structure and therefore a unique neutron star mass is expected. Accretion-induced collapse in binary stars with an O-Ne-Mg white dwarf is an alternative avenue to EC supernovae. Similarly low neutron star masses were obtained in simulations of the collapse of low-mass iron core progenitors [27, 28, 29]. Recent simulations have in fact pushed the lower limit for predicted baryonic neutron star masses down slightly to $1.347\,\mathrm{M}_{\odot}$ [30]. Even lower values were reported in 2D simulations [31], but these were obtained with simplified neutrino transport for stellar evolution models based on artificial bare C/O cores with unrealistically low final C/O core masses $<1.45\,\mathrm{M}_{\odot}$ for such objects. This introduces an artificial bias in the remnant masses. Fixing the C/O core mass eliminates essential processes for core growth or core shrinkage by dredge-up, dredge-out [8, 32] and mass transfer that determine the final fate and the ultimate core sizes at the boundary between SAGB stars and massive stars. Even for ultra-stripped stars that lose much of their He shell by binary interaction, further C/O core growth after the ignition of C burning is important, and brings the final C/O core mass to $\gtrsim 1.45\,\mathrm{M}_{\odot}$ in stars that evolve to iron core collapse.

This tension between observations and theory leaves open several interpretations. Are stellar evolution models not capturing the core structure of the least massive supernova progenitors correctly? Do supernova explosions develop faster than predicted in current models to allow a smaller mass cut? Or are the observed low-mass objects indeed white dwarfs instead of neutron stars [19]?

One key priority for supernova simulations in answering these questions is to better scan the mass range of progenitors that produce the lightest neutron stars. Even though a number of supernova simulations of electron-capture supernovae and low-mass iron core collapse supernovae have been carried out [33, 34, 27, 35, 28, 36, 37, 38, 39, 40, 41], relatively few progenitor models are available to date due to technical difficulties that plague pre-collapse stellar evolution in this region of parameter space. The final stages of the progenitors in this regime are characterized by highly degenerate conditions in and around the core and complex burning behavior that can involve off-center ignition, convectively bounded flames, and sometimes powerful burning flashes with pre-supernova mass ejection [7]. Crucially, the core structure does not depend monotonically on initial progenitor mass [42]. This means that low-mass iron core collapse supernovae could form lighter neutron stars than electron capture supernovae. A fine grid of stellar evolution models is required to scan the variations in core size just above the iron-core formation threshold.

In this paper, we improve the theoretical lower limit for the neutron star mass compatible with current stellar evolution and supernova explosion theory. Starting from a finely spaced grid of low-mass iron-core progenitor models, we select five suitable candidates for particularly low neutron star masses. We then simulate the collapse and, if applicable, the explosion of these progenitor stars in three dimensions (3D) and predict the associated neutron star birth properties.

Figure 1: Entropy profiles of all 25 low-mass iron core progenitor models near the potential mass cut. Five models that have been identified as suitable candidates for three-dimensional supernova simulations are shown as thick colored lines. The remainder of the set are shown as thin gray lines, with lighter shades indicating higher ZAMS mass.

Progenitor Models.—We consider 25 single-star solar-metallicity progenitor models covering the ZAMS mass range from $9.45\,\mathrm{M}_{\odot}$ to $9.95\,\mathrm{M}_{\odot}$ and C/O core masses from $1.481\,\mathrm{M}_{\odot}$ to $1.585\,\mathrm{M}_{\odot}$. The models have been calculated with the stellar evolution code Kepler [43, 44], including fixes to the pair neutrino rates as applied previously in [42, 45]. The resolution in ZAMS mass is as small as $0.01\,\mathrm{M}_{\odot}$ where possible, but due to the complications of highly degenerate core evolution, in particular at the low-mass end, only 25 models were followed until collapse. In the early evolution stages, a dynamic adaptive nuclear reaction network is used that follows all necessary stable and unstable isotopes up to polonium (typically, some 1,800 species for the pre-supernova model), silicon burning and evolution of the iron core use a $\mathord{\sim}125$ isotope intermediate and full nuclear statistical equilibrium network for efficient advection of conserved quantities while reasonably capturing deleptonization and neutrino emission as well as energy in the nuclear excited states.

Due to the computer time requirements of 3D core-collapse supernova simulations with neutrino transport, a subset of five progenitor models was selected for follow-up. Suitable candidates for a low neutron star mass were identified based on the entropy profiles of the progenitors (Figure 1). The onset of the explosion is typically associated with entropy and density jumps at shell interface, in particular the pronounced entropy jump at the base of the oxygen-burning shell in many models [46, 47, 48, 49]. Figure 1 shows the biggest entropy jump for the $9.90\,\mathrm{M}_{\odot}$ model at just above $1.3\,\mathrm{M}_{\odot}$. Sizable entropy jumps at smaller mass coordinate $m$ occur, but are less prominent, and may not be sufficient to trigger an explosion. We select the $9.9\,\mathrm{M}_{\odot}$, $9.73\,\mathrm{M}_{\odot}$, $9.71\,\mathrm{M}_{\odot}$, $9.67\,\mathrm{M}_{\odot}$, and $9.59\,\mathrm{M}_{\odot}$ progenitors as five candidates with progressively smaller iron core masses of $1.31\,\mathrm{M}_{\odot}$, $1.30\,\mathrm{M}_{\odot}$, $1.29\,\mathrm{M}_{\odot}$, $1.27\,\mathrm{M}_{\odot}$, and $1.23\,\mathrm{M}_{\odot}$, but progressively lower likelihood of an explosion triggered by the entropy jump.

Of particular interest is the final evolution that leads to the small core sizes. Post helium burning, all models undergo 4 episodes of carbon core and shell burning before neon/oxygen burning is ignited off-center at $\mathord{\sim}0.19\,\mathrm{M}_{\odot}\ldots 0.32\,\mathrm{M}_{\odot}$ for the most to the least massive models and burns all the way downward to the center. The $9.9\,\mathrm{M}_{\odot}$ model first ignites core silicon burning slightly off-center, the other models ignite silicon burning directly in the center. All models encounter an oxygen shell burning phase before igniting silicon shell burning. In the $9.90\,\mathrm{M}_{\odot}$ and $9.73\,\mathrm{M}_{\odot}$ models, it is quenched well before silicon shell ignition, and sufficient oxygen remains to re-ignite and drive a powerful convective oxygen burning shell minutes ($9.90\,\mathrm{M}_{\odot}$) to seconds ($9.73\,\mathrm{M}_{\odot}$ and below) before core collapse just at the outer edge of the iron core, giving rise to the huge jumps in entropy seen in Figure 1. In the lower-mass models, oxygen is increasingly depleted in the pre-silicon burning phase of oxygen shell burning such that the that pre-collapse oxygen shell results in smaller entropy jumps.

Numerical Methods.—We perform three-dimensional supernova simulations with the relativistic neutrino hydrodynamics code CoCoNuT-FMT code [50]. CoCoNuT-FMT solves the relativistic equations of hydrodynamics in spherical polar coordinates using the piecewise parabolic method [51] for reconstruction and a hybrid HLLC/HLLE Riemann solver [52, 53]. A mesh coarsening scheme is used to avoid time step constraints near the axis of the spherical polar grid [29], and the innermost region of the grid is treated in spherical symmetry. Neutrinos are treated with the fast multi-group transport (FMT) scheme of [50], and the equations for the space-time metric are solved in the extended conformal flatness approximation [54]. We use a grid resolution of $550\times 128\times 256$ zones in radius, latitude, and longitude (corresponding to $1.4^{\circ}$ in angle) and $21$ exponentially spaced energy group from $4\,\mathrm{MeV}$ to $240\,\mathrm{MeV}$. We use the SFHo equation of state of [55] in the high-density regime. At low densities, the equation of state includes nucleons and nuclei (treated as a perfect gas), leptons and radiation, and nuclear burning is treated with a “flashing” scheme [56].

Figure 2: Average shock radius (top), criticality ratio $\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}$ between the advection and heating time scale, and baryonic neutron star mass (bottom) as a function of time for the five core-collapse supernova simulations.

Results.—Shock trajectories for the five models are shown in Figure 2. The $9.9\,\mathrm{M}_{\odot}$ and $9.71\,\mathrm{M}_{\odot}$ models are the earliest to explode. The $9.59\,\mathrm{M}_{\odot}$ model fails to explode early when the shell interface at $1.23\,\mathrm{M}_{\odot}$ reaches the shock, though it may still explode later, perhaps triggered by the more pronounced entropy jump at $1.4\,\mathrm{M}_{\odot}$. The critical threshold for shock revival by neutrinos is often expressed via the ratio between the advection time scale $\tau_{\mathrm{adv}}$ and the heating time scale $\tau_{\mathrm{heat}}$ [57, 58]. The critical ratio $\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}$ shows the same hierarchy, with the $9.71\,\mathrm{M}_{\odot}$ model reaching $\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}=1$ the earliest. Even though the $9.67\,\mathrm{M}_{\odot}$, $9.71\,\mathrm{M}_{\odot}$, and $9.73\,\mathrm{M}_{\odot}$ models initially have slightly higher values for a brief period, the $9.9\,\mathrm{M}_{\odot}$ simulation eventually shows the fastest shock propagation, followed very closely by the $9.71\,\mathrm{M}_{\odot}$ case. Evidently, the more rapid drop of the accretion rate associated with the bigger entropy jump in the $9.9\,\mathrm{M}_{\odot}$ model is more critical for the rapid development of an explosion than the smaller mass coordinate of the jump in the $9.67\,\mathrm{M}_{\odot}$ and $9.73\,\mathrm{M}_{\odot}$ models.

By the end of the simulation, the $9.9\,\mathrm{M}_{\odot}$ model has reached an explosion energy of $1.5\times 10^{50}\,\mathrm{erg}$, which is significantly higher than for the $9.67\,\mathrm{M}_{\odot}$ model with $4.7\times 10^{49}\,\mathrm{erg}$, and somewhat higher than for the $9.71\,\mathrm{M}_{\odot}$ model at $10^{50}\,\mathrm{erg}$. The $9.73\,\mathrm{M}_{\odot}$ model has the lowest energy by the end of the simulation, but has been terminated earlier than the other simulations and is likely to end up between the $9.67\,\mathrm{M}_{\odot}$ and the $9.71\,\mathrm{M}_{\odot}$ case.

The key outcome for the purpose of this paper is the neutron star mass. Figure 2 (bottom panel) shows the evolution of the baryonic neutron star mass, $M_{\mathrm{by}}$, which is obtained by integrating the entire mass on the grid at densities above $10^{11}\,\mathrm{g}\,\mathrm{cm}^{-3}$. Reflecting the more precipitous explosion in the $9.9\,\mathrm{M}_{\odot}$ model, the neutron star mass quickly saturates at $1.313\,\mathrm{M}_{\odot}$ in this case, whereas $M_{\mathrm{by}}$ still continues to increase in the other three exploding models. By the end of the simulation, the neutron star is already losing mass at a very small rate in the $9.9\,\mathrm{M}_{\odot}$ model, so the value of $M_{\mathrm{by}}=1.313\,\mathrm{M}_{\odot}$ is a rather firm upper limit for the final neutron star mass, barring the possibility of later fallback.

Figure 3: Total kick velocity (thick solid line) and the $x$- (thin black lines), $y$- (red), and $z$-components (blue) of the kick due to the gravitational tug-boat mechanism (hydrodynamic kick, thin solid lines) and due to anisotropic neutrino emission (dashed lines) for the $9.9\,\mathrm{M}_{\odot}$ model.

Using the fit formula for the cold neutron star binding energy from [59], $M_{\mathrm{by}}$ can be translated into the final gravitational mass $M_{\mathrm{grav}}$,

$$M_{\mathrm{by}}=M_{\mathrm{grav}}-0.084\,\mathrm{M}_{\odot}\left(M_{\mathrm{grav}}/\mathrm{M}_{\odot}\right)^{2},$$ (1)

which results in $M_{\mathrm{grav}}=1.192\,\mathrm{M}_{\odot}$. This is still in conflict with the mass measurement of $M_{\mathrm{grav}}=1.174\,\mathrm{M}_{\odot}$ for the lighter companion in J0453+1559. If we assume $2\,\sigma$ error bars instead of $1\,\sigma$ error bars, the results of [16] may be compatible with a mass as high as $1.182\,\mathrm{M}_{\odot}$, but even then a tension with observations remains. It is conceivable that stripped stars may have even slightly more favorable structure for low-mass neutron stars as their expected higher carbon mass fraction at the end of core helium burning may lead to more neutronisation in the supernova progenitor [60]. Something similar may occur for more metal-rich progenitors. Nonetheless, the discrepancy is very small in absolute terms and may be within the range of physical uncertainties in supernova progenitor and explosion modeling (see Conclusions).

Compatibility with the system parameters of J0453+1559, specifically the eccentricity and orbital period, also requires a suitable kick and amount of mass loss in the explosion. A small kick, as typically expected for electron-capture supernovae, is required [16]. For their alternative scenario involving the formation of a massive white dwarf and the ejection of several $0.1\,\mathrm{M}_{\odot}$ in a thermonuclear electron-capture supernova, Tauris & Janka [19] estimate a required kick of $\mathord{\sim}70\,\mathrm{km}\,\mathrm{s}^{-1}$. With a single-star progenitor we cannot estimate the realistic amount of mass ejection in a binary scenario where the progenitor would have undergone envelope stripping earlier in its evolution. It is plausible, however, that the resulting envelope mass in a binary scenario and the required kick will be of a similar scale in case of formation of the light companion by a low-mass iron core progenitor.

We therefore evaluate the kick for the $9.9\,\mathrm{M}_{\odot}$ model to demonstrate that the model is roughly consistent with these requirements with a kick of order $\mathord{\sim}100\,\mathrm{km}\,\mathrm{s}^{-1}$. Since the kick velocity is subject to stochastic model variations by at least a factor of two in both directions, only the rough scale of the kick is relevant. We compute the neutron star momentum from the vectorial momentum of the ejecta $\mathbf{p}_{\mathrm{ej}}$, which arises due to the gravitational tug-boat mechanism [61], and the emitted neutrinos by invoking momentum conservation [61, 62]. In terms of the relativistic momentum density $\mathbf{S}$ and the neutrino energy flux $\mathbf{F}$ (measured at a radius of $400\,\mathrm{km}$), we compute the kick velocity $\mathbf{v}_{\mathrm{kick}}$ at time $t$ as

$$M_{\mathrm{grav}}\mathbf{v}_{\mathrm{kick}}(t)=-\int\limits_{\mathrm{ejecta}}\mathbf{S}\,\mathrm{d}V-\int\limits_{0}^{t}\!\!\oint\frac{\mathbf{F}(r=400\,\mathrm{km})}{c}\,\mathrm{d}A,$$ (2)

where the relativistic volume and surface elements are to be used. By the end of the simulation, the kick velocity $|\mathbf{v}_{\mathrm{kick}}|$ has reached a value of $100\,\mathrm{km}\,\mathrm{s}^{-1}$ and is still growing, albeit at a decelerating pace. The kick from asymmetric neutrino emission is subdominant, but slightly reduces the total kick as it points in the direction opposite to the hydrodynamic kick. The $9.71\,\mathrm{M}_{\odot}$ model shows a similar growth of the kick velocity, which demonstrates that kicks of this order are not predicted to be unusual in this progenitor mass range.

Conclusions.— Our simulations set a new record for the lowest neutron star mass obtained in 3D supernova simulations with multi-group neutrino transport. The lowest gravitational neutron star mass compatible with current stellar evolution models is now at least as low as $1.192\,\mathrm{M}_{\odot}$. Despite the rapid explosion of the $9.9\,\mathrm{M}_{\odot}$ progenitor and the small amount of accretion after shock revival, we are able to obtain a kick of $100\,\mathrm{km}\,\mathrm{s}^{-1}$, which is mostly due to the gravitational tug-boat mechanism [61] with a small additional contribution from anisotropic neutrino emission. There is still a small tension with the $2\,\sigma$ error bars for the $1.174\,\mathrm{M}_{\odot}$ object in J0453+1559, whose neutron star nature has recently been questioned [19]. The tension is so small, however, so that minor variations in the progenitor evolution and supernova dynamics may well push the minimum neutron star mass down by a further $\mathord{\sim}0.01\,\mathrm{M}_{\odot}$ and resolve this discrepancy. For example, convective seed perturbations in the progenitor could result in a slightly earlier onset of the explosion [63, 64, 65] and reduce the neutron star mass further. The intricacies of degenerate nuclear burning may well introduce uncertainties in the pre-collapse core structure that also result in slightly smaller neutron star masses. Even minor details in the pre-collapse models, such as the usual non-rigorous treatment of the small effect of relativistic gravity in current stellar evolution code, could affect the pre-collapse evolution on a level of accuracy that matters for the comparison to J0453+1559 and other future neutron star candidates with precise mass determination in this mass range.

Regardless of whether the $1.174\,\mathrm{M}_{\odot}$ object in J0453+1559 can eventually be explained as a neutron star, our new simulations hold two important insights. Contrary to long-standing expectations, the lightest neutron star masses appear not to be made by electron-capture supernovae, but by iron-core collapse supernovae. Furthermore, models predict considerable non-monotonicity in the explosion and remnant properties for iron-core collapse supernovae near the minimum progenitor mass. The lightest neutron stars may originate from stars several $0.1\,\mathrm{M}_{\odot}$ above the progenitor mass threshold (in terms of ZAMS mass) for core-collapse supernovae. Furthermore, there is no strict correlation between neutron star masses and kicks. We predict that the lightest neutron stars can still have substantial (albeit below-average) kicks. Combining increasingly detailed and precise predictions from 3D explosion models and more precise neutron star mass and kick determinations at the low-mass end of the distribution holds considerable promise for testing and validating our understanding of supernova physics and stellar evolution in this challenging progenitor mass regime.

Data availability— The data from our simulations will be made available upon reasonable requests made to the authors.

Acknowledgments— The authors acknowledge support by the Australian Research Council through grants FT160100035 (BM) DP240101786 (BM, AH) DE210101050 (JP) and DP240103174 (AH). The authors acknowledge support by the Australian Research Council (ARC) Centre of Excellence (CoE) for Gravitational Wave Discovery (OzGrav) project numbers CE170100004 and CE230100016. This project was carried out using computer time allocations from Astronomy Australia Limited’s ASTAC scheme, the National Computational Merit Allocation Scheme (NCMAS), and from an Australasian Leadership Computing Grant. Some of this work was performed on the Gadi supercomputer with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government, and through support by an Australasian Leadership Computing Grant.