A publishing partnership

The following article is Open access

R-process Nucleosynthesis of Subminimal Neutron Star Explosions

, , , and

Published 2023 October 13 © 2023. The Author(s). Published by the American Astronomical Society.
, , Citation Chun-Ming Yip et al 2023 ApJ 956 115DOI 10.3847/1538-4357/acf570

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/956/2/115

Abstract

We show that a minimum-mass neutron star undergoes delayed explosion after mass removal from its surface. We couple the Newtonian hydrodynamics to a nuclear reaction network of ∼4500 isotopes to study the nucleosynthesis and neutrino emission during the explosion. An electron antineutrino burst with a peak luminosity of ∼3 × 1050 erg s−1 is emitted while the ejecta is heated to ∼109 K. A robust r-process nucleosynthesis is realized in the ejecta. Lanthanides and heavy elements near the second and third r-process peaks are synthesized as end products of nucleosynthesis, suggesting that subminimal neutron star explosions could be an important source of solar chemical elements.

Export citation and abstractBibTeXRIS

1. Introduction

Neutron star binaries have been intensively investigated for several decades, both observationally and theoretically. The discovery of pulsar PSR B1913+16 in a neutron star binary by Hulse & Taylor (1975), and its orbital decay due to gravitational wave radiation (Weisberg et al. 1981) revealed the evolution of close neutron star binaries. The first-ever multimessenger detection of a binary neutron star merger in 2017 (Abbott et al. 2017; Goldstein et al. 2017; Savchenko et al. 2017) confirmed the encounter of neutron stars in close binaries. Apart from direct merger, Clark & Eardley (1977) discussed the mass transfer in a close neutron star binary with different initial masses. It was found that stable mass transfer can be established under certain circumstances, which would further enhance the asymmetry of the systems. Blinnikov et al. (1984) subsequently derived that such mass transfer caused by Roche lobe overflow is slow and stable until the mass of the lighter and larger neutron star reaches ∼0.15 M. After that, tidal disruption of the minimum-mass neutron star occurs with a characteristic timescale longer than its hydrodynamic timescale, so it evolves through a series of quasi-equilibrium states until reaching the minimum stable mass ∼0.1 M allowed by the equation of state (EoS), and it explodes (Blinnikov et al. 1990, 2021). A recent postprocessing nucleosynthesis calculation by Panov & Yudin (2020) revealed the production of heavy elements through this event.

The nuclear reactions involved in the expansion of neutron star matter were examined by Lattimer et al. (1977). Simulations of subminimal neutron star explosions by Colpi et al. (1989, 1991, 1993) showed the essential roles of nuclear reactions in the development of hydrodynamic instabilities. The timescale of weak interactions is considerably longer than the hydrodynamic timescale, and thus β-equilibrium cannot be always established dynamically when a neutron star is perturbed. The pressure of neutron star matter is not uniquely determined at a given density if the β-equilibrium is not achieved. Hence, the minimum mass of neutron stars determined from the EoS, in which the β-equilibrium condition is often assumed, does not necessarily correspond to the onset of instability. Hydrodynamic simulations of such a configuration with consideration of departure from β-equilibrium were performed. It was demonstrated that a sufficiently large ratio of mass removal, ΔS, at the surface of the minimum-mass neutron star, can lead to an instability. Sumiyoshi et al. (1998) further presented the explosion of the subminimal neutron star in a simulation with an initial mass removal ratio ΔS = 0.22. An electron antineutrino burst with peak luminosity ∼1052 erg s−1 was predicted. Moreover, r-process was realized in the expanding neutron star matter when its density dropped rapidly, and rare-earth elements were produced.

Despite the inspiring predictions of neutrino burst and nucleosynthesis, several modifications of the model used by Sumiyoshi et al. (1998) are necessary for a more rigorous investigation. In this article, we revisit the reactive-hydrodynamic simulations of subminimal neutron star explosions with a more modern nuclear matter EoS. A nuclear reaction network, instead of the single nucleus approximation, is employed to compute the relevant nuclear reactions, including neutron capture, β-decay, and spontaneous fission, that are essential for nucleosynthesis. We predict the yield of stable elements produced through this event, and we show that lanthanides and heavy elements near the second and third r-process peaks are synthesized. Our results suggest that the explosions of subminimal neutron stars could be a potentially new class of astronomical r-process events in addition to other known heavy element production channels (e.g., Kajino et al. 2019; Kobayashi et al. 2020).

This paper is organized as follows: In Section 2, we present the computational setup of the hydrodynamic simulations coupled with the nuclear reaction network. In Section 3, we present the numerical results of the simulations and the nucleosynthesis calculations. We discuss and conclude the findings in Section 4.

2. Methodology

2.1. Initial Configuration

We adopt the EoS by Schneider, Roberts, and Ott (SRO; Schneider et al. 2017, 2019), computed using the NRAPR Skyrme parameterization (Steiner et al. 2005), an effective nonrelativistic Skyrme-type nuclear interaction model (Lattimer & Douglas Swesty 1991), and nuclear statistical equilibrium (NSE) approximation at low densities and temperatures, available in CompOSE, 3 in our study. We construct the initial configuration of the minimum-mass neutron star with a uniform temperature T = 108 K under hydrostatic and β-equilibrium. A pure Newtonian gravity and spherical symmetry are assumed. The minimum-mass neutron star has a central density of about 1.211 × 1014 g cm−3, a total mass M of about 0.0896 M, and a radius R of more than 250 km (see Figure 1). The general relativistic correction to the minimum-mass neutron star structure is negligible as its compactness C = GM/Rc2 ∼ 5 × 10−4 is very small.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Initial configuration of the minimum-mass neutron star at hydrostatic and β-equilibrium for the EoS by Schneider, Roberts, and Ott (SRO; Schneider et al. 2017, 2019) at uniform temperature T = 108 K. The density (solid line) and electron fraction (dashed line) are plotted vs. the radial coordinates in logarithm scale.

Standard image High-resolution image

2.2. Hydrodynamic Simulation

The set of Euler equations in the Lagrangian formalism, in spherical symmetry, are:

where G is the gravitational constant. The baryon mass m contained within the radius r from the center of the neutron star is chosen to be the Lagrangian mass coordinates, with the outward radial direction from the center defined as the positive direction. u, ρ, p, and ε are the radial velocity, baryon mass density, pressure, and specific internal energy density of the fluid, respectively, as functions of time t and m. A standard first-order Lagrangian finite-difference scheme (Richtmyer et al. 1967) with 100 uniform mass grid zones is employed to solve the Euler equations numerically with the introduction of an artificial viscosity term q whenever compression occurs, where l = 1.0 G M/c2 = 1.477 km is the artificial viscosity coefficient having the dimension of a length. Explicit first-order time discretization subjected to the Courant–Friedrichs–Lewy (CFL) condition is chosen to update these hydrodynamic variables. In Equation (4), , , , , and denote the extra heating/cooling rates of the mass elements associated with thermonuclear reactions, β-decays, β-decays neutrino emissions, spontaneous fission, and thermal neutrino emissions, respectively (see Section 2.3).

2.3. Nuclear Reaction Network

We implement and adopt the open-source code for nuclear reaction network computation developed by Timmes et al. (2000). A network consisting of ∼4500 isotopes (see Figure 2) is chosen to study nucleosynthesis in this work. The atomic number Z of the isotopes included in the network ranges from for neutron to for uranium. The stable elements synthesized through standard r-process channels and their short-lived neutron-rich isotopes temporarily produced during the r-process are included. The nuclear masses, partition functions, and thermonuclear reaction rates of pair reactions with proton, neutron, α particle, and photon, namely (n, γ), (n, p), (p, γ), (α, n), (α, p), (α, γ) processes, and their inverse processes, are used (Cyburt et al. 2010). Some special reactions, such as deuterium fusion and carbon burning, are implemented additionally. The thermal energy generation/absorption rate due to the thermonuclear reactions is computed by

where NA is the Avogadro constant, and Yi (Mi ) is the number fraction (nuclear mass) of isotope i. We sum over all isotopes i in the network to find the net change in Yi for all isotopes after the thermonuclear reactions.

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Isotopes included in the nuclear reaction network. The colored region covers the isotopes considered in the network, and the grids in red further illustrate the stable isotopes in the solar system (Lodders 2019). Extremely neutron-rich isotopes with proton number Z > 80 are not included because some of the nuclear reaction rates are not available to form linkages with other isotopes included.

Standard image High-resolution image

The β-decay half-lives and β-delayed neutron emission probabilities of nuclei provided by Möller et al. (2003) are adopted. The energy available Δj through the β-decay reaction j of an isotope with nuclear mass M(A, Z), where A denotes the mass number of the isotope, is given by

The chemical potential of electron μe enters the equation, since the Fermi level of degenerate electrons is high inside a neutron star. The endothermic β-decays are forbidden whenever μe is higher than the energy released. The average energy carried away by each electron antineutrino can be expressed as (Sumiyoshi et al. 1998)

We take the approximation that the neutrinos produced through β-decays escape freely from the subminimal neutron star without any dissipation of energy (Meyer 1989). Hence, the thermal energy generation rate and the neutrino cooling rate with respect to β-decays are computed by

where summation over all β-decay reactions j, with Δj > 0 under the given μe is performed. After the weak interactions, the heating effect owing to β-delayed neutron emissions is also calculated, which is found to be an insignificant contribution and is counted in for convention.

The spontaneous fission half-lives of heavy nuclei τf (in a unit of yr) with atomic number Z ≥ 90 in the network are evaluated using the semiempirical formula by Santhosh et al. (2010):

where N is the neutron number of the nuclei, and a = −43.25203, b = 0.49192, c = 3674.3927, d = −9360.6, and e = 580.75058 are fitting parameters. We approximate the fission fragments after spontaneous fission by the equation

when a nucleus undergoes spontaneous fission to release two daughter nuclei, and . The fission fragment asymmetry parameter γ controls the asymmetry of fission fragments. We set the values of γ to be 0.50, 0.45, and 0.40 in different models to investigate the influence of the fission fragment asymmetry on the subminimal neutron star explosions. The models with γ = 0.50 execute symmetric fission fragment approximation, while the other two classes of model roughly resemble the peaks of fission product yield distribution from more sophisticated fission fragment calculations (e.g., Hao et al. 2022). After calculating the change in Yi of all isotopes i in the network caused by spontaneous fission, the thermal energy generation rate through spontaneous fission is computed by

The thermal neutrino energy loss rate through pair-, photo-, plasma-, bremsstrahlung, and recombination neutrino processes are calculated using the open-source subroutine for thermal neutrino emission, 4 which adopts the analytical fitting formulae by Itoh et al. (1996). It enters Equation (4) under the assumption that the thermal neutrinos escape freely from the neutron star. It is found that the cooling effect caused by thermal neutrino emission is negligible as the temperature reached by the exploding subminimal neutron star is not adequately high.

2.4. Composition Evolution

The chemical composition evolution of neutron star matter is computed individually for each mass element. The Maxwell–Boltzmann equation is solved under the given temperature, density, and electron fraction to assign the mass fraction of each isotope in the nuclear reaction network as the initial composition at NSE, which is shown in Figure 3. A threshold density ρth ≡ 1013 g cm−3 is defined. A mass element is regarded as core-like (crust-like) matter if its density is above (below) ρth. Above ρth, the neutron star matter is composed of predominantly free neutrons, with a mass fraction ≳0.9, and nuclei forming exotic nuclear structures, such as pasta-like configurations, are immersed. The β-decay channels of nuclei are all blocked by the high electron chemical potential (μe ≳ 30 MeV, see Equation (7)). Therefore, we do not determine the composition of mass elements with a density exceeding ρth, and we neglect the nuclear reactions involved. The network is not activated for the mass elements composed of core-like matter in the hydrodynamic simulations.

Figure 3. Refer to the following caption and surrounding text.

Figure 3. Initial mass fractions of isotopes at NSE for the neutron star crust-like matter. Ni62 and Fe56 are the major isotopes near the surface of the neutron star. The electron fraction decreases gradually as density increases, and neutron-rich isotopes appear correspondingly. Free neutrons are dripped when the density is above ∼4 × 1011 g cm−3 and become the predominant composition near the threshold density ρth ≡ 1013 g cm−3.

Standard image High-resolution image

The maximum temperature reached during the explosion of a subminimal neutron star is not high enough for efficient neutrino capture (Meyer et al. 1998) or positron capture (Ruffini et al. 2010) to contribute significantly to leptonization of the neutron star matter. The modified URCA processes are the major weak interactions between the free neutrons and protons, which are characterized by a timescale much longer than the relevant hydrodynamic timescale (Colpi et al. 1989). Hence, the β-decays of nuclei discussed in Section 2.3 are the only weak interactions that alter the electron fractions of the mass elements. As a result, the electron fractions of mass elements regarded as core-like matter are fixed to the initial values obtained assuming β-equilibrium. When the density of a mass element drops below ρth, it transforms from core-like matter to crust-like matter. It is assumed that NSE is established after the transition from an exotic nuclear structure, such as pasta-like configurations, to a mixture of nuclei, free neutrons, and electrons as crust-like matter at that time step. The mass fractions of isotopes at NSE are assigned to be the composition of the mass element at the temperature, density, and electron fraction after the transition (see Figure 4). After that, the composition of the mass element is obtained by solving the network coupled to the Euler equations.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. Mass fractions of isotopes at NSE and threshold density ρth ≡ 1013 g cm−3. The solutions at different temperatures and electron fractions Ye are shown. The mass fraction of the free neutron is ∼0.9 (not shown in the figure) for the conditions concerned. For Ye between 0.016 and 0.030, the solutions at NSE are alike.

Standard image High-resolution image

The full network with all nuclear reactions considered is solved together with the Euler equations in order to investigate their effects consistently. It is believed that the leptonization owing to β-decays of nuclei and an overall net thermal energy generation by nuclear reactions enhance the expansion of the crust-like matter significantly. Hence, the leptonization and thermal energy generation are crucially responsible for the onset of hydrodynamical instability. The subminimal neutron star becomes unstable and explodes ultimately when it can no longer adjust itself against radial oscillations, and all mass elements are ejected. We perform the hydrodynamic simulation until the density of all the mass elements drops near the minimum value available by the EoS table (∼103 g cm−3). After that, we extrapolate 5 the density and temperature of the ejecta assuming homologous and adiabatic expansion described by ρ(t) ∝ t−3 and T(t) ∝ t−1 as functions of time t. We update the network for 1 Gyr after the hydrodynamic simulation to determine the stable elements produced by the explosion of the subminimal neutron star.

3. Results

We have checked that the initial model of the minimum-mass neutron star constructed is stable in the hydrodynamic simulation coupled to the nuclear reaction network if no perturbation is imposed. To initiate an instability, we remove several mass elements from the surface of the star. We vary ΔS, the initial mass removal ratio to the total mass of the star, from 0.40 to 0.30 in order to mimic the unstable mass transfer in the final approach of the binary system. Also, we set the fission fragment asymmetry parameter γ to be 0.50, 0.45, and 0.40 to study the effect of asymmetric fission fragment approximation (see Section 2.3). The models are named according to ΔS and γ as L-ΔS-γ.

3.1. Explosion of Subminimal Neutron Star

The hydrodynamic simulation results of the model L-40-50 are displayed in Figure 5. By imposing the initial mass removal, all the neutron star crust-like matter and part of the neutron star core-like matter are removed. The pressure gradient on the surface of the exposed core-like matter increases suddenly, and hence hydrostatic equilibrium is no longer maintained. The mass elements expand sequentially from the surface to the center. The perturbation made here, however, does not lead to a direct explosion of the whole configuration. Most of the mass elements oscillate and settle down at a larger radial position and a lower density temporarily after ∼0.01 s. The subminimal neutron star remains at quasi-equilibrium. If no nuclear reaction is included, the star oscillates around the new equilibrium position, with most of the mass elements remaining gravitationally bounded.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Radial positions (upper left panel), densities (lower left panel), temperatures (upper right panel), and velocities (lower right panel) of mass elements vs. time in logarithm scale in the model L-40-50. The temperatures are indicated by dashed–dotted lines prior to the activation of the nuclear reaction network to signal potentially large uncertainty due to the EoS construction ambiguity (see footnote 6). After the mass removal, the mass elements near the surface of the exposed neutron star core-like matter start expanding promptly, though most of them remain gravitationally bound. They fall back onto the neutron star at ∼0.01 s and undergo radial oscillation. A new equilibrium of the subminimal neutron star with lower central density is temporarily found between ∼0.01 and 0.05 s. Nonetheless, the nuclear reaction network is activated as the density of the mass elements is lower than the threshold density ρth ≡ 1013 g cm−3. The heating effect and leptonization caused by the nuclear reactions alter the structure of the neutron star in a quasi-equilibrium state. The star loses stability after ∼0.05 s, leading to a delayed explosion.

Standard image High-resolution image

Meanwhile, the nuclear reaction network is activated, as the densities of these mass elements are now lower than ρth. The temperatures and electron fractions of the mass elements during the transition from core-like matter to crust-like matter are plotted in Figure 6. As discussed in Section 2.4, the NSE isotope abundances right after the transition are assigned. Note that the temperatures of some mass elements are much lower than 2 × 109 K, the typical minimum temperature that the NSE approximation is regarded as valid during the transition. We assume that the NSE approximation still holds for these mass elements. Figure 7 illustrates the evolution of the isotope abundance of a mass element in the 30th percentile in mass coordinates after the transition, as an example. The differences among the assigned compositions of the mass elements at different temperatures are eliminated by a robust r-process very soon. Therefore, the simulation results are not noticeably affected, even if the NSE approximation may not be accurate at low temperatures.

Figure 6. Refer to the following caption and surrounding text.

Figure 6. Temperatures and electron fractions of the mass elements during the transition from core-like matter to crust-like matter in the model L-40-50. When the density of mass elements drops from ∼1014 g cm−3 to ρth ≡ 1013 g cm−3 due to expansion, a transition from core-like matter to crust-like matter occurs.

Standard image High-resolution image
Figure 7. Refer to the following caption and surrounding text.

Figure 7. Evolution of the isotope abundances of the mass element in the 30th percentile in mass coordinate vs. time in the model L-40-50. After the transition at ∼0.0006 s, As115 is the major nucleus present at NSE. The nuclear reactions, predominately neutron capture and β-decay, promote the subsequent formations of Se116, Br117, Kr118, and heavier isotopes within a short timescale. The mass fraction of free neutron drops only by ∼0.5% (not shown in the figure) within the duration shown. For other mass elements with lower temperatures during the transition, either Se116 or Kr118 is the major nucleus instead. These two isotopes serve as seed nuclei that undergo r-process with very similar evolutionary tracks as the one shown for As115. As a result, whether Kr118, Se116, or As115 is the major nucleus at NSE does not significantly affect the r-process later on.

Standard image High-resolution image

The r-process proceeds efficiently in all mass elements with evolutionary tracks analogous to the one displayed in Figure 7. The nuclear reactions increase the internal energy and electron fractions of these mass elements, which provide higher thermal pressure and electron degeneracy pressure that alter the hydrodynamic properties of the neutron star. The cumulative effects promote the expansion of mass elements, and the subminimal neutron star loses stability gradually. A delayed explosion of the whole configuration begins at ∼0.05 s, with all mass elements becoming unbounded as ejecta. The radial velocity of the ejecta increases monotonically after the explosion due to the thermal energy generation and leptonization associated with the ongoing nuclear reactions. The total kinetic energy and the magnitude of the gravitational potential of the ejecta are ∼1.4 × 1050 erg and ∼6 × 1047 erg, respectively, before the density of the outermost mass element drops below the minimum density allowed by the EoS table (at ∼0.5 s after the mass removal).

When the subminimal neutron star expands, a rapid rise in temperature up to ∼1010 K, 6 while the network is not yet activated, is observed at densities above ρth (see Figure 5). Correspondingly, the mass elements reach ρth at different temperatures, as illustrated in Figure 6. The thermodynamic properties of neutron star matter, however, are insensitive to temperatures below ∼1010 K in our simulations. 7 Since the temperature rise may be due to ambiguities in the EoS (see footnote 6) and has no effect on the hydrodynamics, we only focus on the nuclear reactions involved below ρth. After the network is activated, the mass elements are powered by the nuclear reactions. The high temperature of the crust-like matter, up to about 6 × 109 K, reached before the neutron star explodes is attributed to not only the heating effect by the nuclear reactions, but also the adiabatic compression during radial oscillations.

The hydrodynamic evolution of the subminimal neutron star is insensitive to the values of γ chosen in the models presented. The explosion timescale, however, varies significantly among the models with different values of ΔS. With a lower ΔS, the explosion is postponed, and the pulsating phase extends substantially. It takes a long time for the nuclear reactions to initiate an instability and lead to the explosion, which is indicated by the rapidly decreasing central density. After that, the star explodes in a similar way as illustrated in Figure 5. The hydrodynamic simulation results of the models with smaller values of ΔS are presented in Appendix A for reference.

3.2. Neutrino Emission and Electromagnetic Radiation

After the explosion of the neutron star, the nuclear reactions, mainly β-decays of nuclei, power the ejecta continuously. The temperature of the ejecta is sustained at ∼109 K and decreases slowly because of adiabatic cooling. The time evolution of electron antineutrino luminosity and the net cumulative thermal energy deposited on the subminimal neutron star by the network in the models L-40-50, L-40-45, and L-40-40 are shown in Figure 8. The electron antineutrino luminosity associated with β-decays of nuclei is predicted using Equation (8) during the explosion of a subminimal neutron star. A peak luminosity of ∼3 × 1050 erg s−1 is reached shortly after the mass removal in the three models. The little trough in the luminosity at ∼0.01 s is caused by the transient rise in electron chemical potential (see Equation (7)) associated with adiabatic compression during radial oscillations of the star. After the peak in neutrino emission, the average luminosity is maintained at ∼1050 erg s−1. Meanwhile, the net cumulative thermal energy deposited on the subminimal neutron star by the network is ∼1050 erg in the three models before the density of the outermost mass element drops below the minimum density provided by the EoS table at ∼0.5 s. The hydrodynamic evolution of models L-40-45 and L-40-40 is very similar to that of the model L-40-50. Also, the neutrino luminosities of the three models are very similar in general. These results are, therefore, insensitive to the variation in fission fragment asymmetry parameter γ. The electron antineutrino luminosity versus time and the net cumulative thermal energy deposited on the neutron star by the network in the models L-40-50, L-38-50, L-36-50, L-34-50, L-32-50, and L-30-50 are also shown in Appendix B as a comparison.

Figure 8. Refer to the following caption and surrounding text.

Figure 8. Upper panel: electron antineutrino luminosity emitted through β-decays of nuclei in hydrodynamic simulation vs. time in the models L-40-50, L-40-45, and L-40-40. A subplot zooms in on the initial peak and shows that the peak luminosity is ∼3 × 1050 erg s−1 in the three models. Lower panel: net cumulative thermal energy injected into the neutron star by the nuclear reaction network vs. time in the same set of models.

Standard image High-resolution image

In the model L-40-50, we evaluate the cumulative energy source terms contributed by β-decays, thermonuclear reactions, spontaneous fission, and neutrino emissions through β-decay of nuclei versus time (see Figure 9). It is found that β-decay is the dominating thermal energy source throughout the simulation, while the thermonuclear reactions and spontaneous fission of nuclei are rather insignificant. From the nucleosynthesis calculations, we notice that the net thermal energy generation rate and neutrino luminosity are sustained until ∼1.4 s and plunge when the mass fraction of free neutron vanishes (see Figure 10 for the evolution of relevant quantities of the mass element in the 30th percentile in mass coordinate as an example). The temperature of the ejecta remains high, and the expansion is continuously powered by the nuclear reactions until the free neutrons are exhausted, and the production of neutron-rich nuclei with short β-decay half-lives is terminated. An electromagnetic peak near soft gamma-ray lasting for a few seconds in total is expected accordingly, assuming blackbody emission. Likewise, the neutrino emission should last for a few seconds and diminish rapidly following the termination of r-process.

Figure 9. Refer to the following caption and surrounding text.

Figure 9. Cumulative energy source terms by the nuclear reaction network vs. time in the model L-40-50. The absolute value of the cooling term corresponding to β-decay neutrino emissions is plotted as a brown dashed–dotted line.

Standard image High-resolution image
Figure 10. Refer to the following caption and surrounding text.

Figure 10. Net thermal energy generation rate, local electron antineutrino luminosity emitted through β-decays of nuclei, and the free neutron mass fraction of the mass element in the 30th percentile in mass coordinate vs. time in the model L-40-50. The solid curves are obtained directly from hydrodynamic simulations, and the dashed curves are obtained by extrapolation of the density and temperature.

Standard image High-resolution image

3.3. Nucleosynthesis

During the explosion of a subminimal neutron star, the neutron star matter in the core, with initially very low electron fraction, is decompressed without experiencing very strong heating that can be achieved in other astronomical r-process sites, such as core-collapse supernovae and binary neutron star mergers, where the temperature can reach the order of 1011 K. The maximum temperature experienced by the mass elements in the simulation is ∼1010 K, which is much lower than the typical temperature scale in the former two scenarios. A robust r-process occurs because of the initially high neutron excess and low electron fraction of the ejecta. Beyond the simulation time, we extend the nucleosynthesis calculations by extrapolating the density and temperature of the ejecta, assuming free expansion. It is noticed that all mass elements experience similar evolutionary tracks of composition in the models L-40-50, L-40-45, and L-40-40. The temperatures and electron fractions of the mass elements in the three models when the network starts operating are listed in Table 1. The inner ejecta, originating from the neutron star core region, has a relatively high temperature and electron fraction compared to the other two categories after the decompression. The intermediate ejecta experiences moderate heating and has a low initial electron fraction. The outer ejecta, with almost no adiabatic compression during the radial oscillations, is initially cold and has moderate electron fraction.

Table 1. Nuclear Reaction Network Properties

Percentile (th) Ye,ini Tini (GK) (GK)Ejecta Category
1–100.02784.0516.310inner
11–200.02443.9665.886inner
21–300.02173.8345.339intermediate
31–400.01913.6514.903intermediate
41–500.02051.9465.6intermediate
51–600.02340.1003.079outer

Notes. Percentile in mass coordinate, averaged initial electron fraction Ye,ini, and averaged initial temperature Tini of the mass elements when the nuclear reaction network starts operation at the threshold density ρth ≡ 1013 g cm−3 in the models L-40-50, L-40-45, and L-40-40. The averaged maximum temperature reached by the mass elements after the activation of the network is listed as well. In most of the mass elements, the peak temperature is reached when the first adiabatic compression wave arrives at ∼0.01 s, while nuclear fission takes place only after ∼0.03 s. Therefore, the maximum temperature reached by the mass elements in the three models is nearly identical.

Download table as:  ASCIITypeset image

The mass fractions of isotopes produced in the models L-40-50, L-40-45, and L-40-40 after 1 Gyr are plotted versus mass number and atomic number in Figures 1113, respectively. The solar abundance distribution by Lodders (2019), scaled to match with the averaged production peak of Xe132 among the first to tenth mass elements, is shown by the black scatters in the figures for comparison. With the occurrence of a robust and long-lasting r-process in all mass elements, the composition evolution becomes insensitive to the tiny difference among the initial conditions of these mass elements. The end products of nucleosynthesis in all mass elements from the same model are very similar.

Figure 11. Refer to the following caption and surrounding text.

Figure 11. Averaged fractions of isotopes vs. mass number A and atomic number Z of every 10 mass elements at the end of nucleosynthesis calculation in the model L-40-50. The solar abundance distribution (Lodders 2019), scaled to match with the averaged mass fraction of Xe132 of the mass elements from the 1st to 10th percentile in mass coordinate, is shown by the black scatters.

Standard image High-resolution image
Figure 12. Refer to the following caption and surrounding text.

Figure 12. Same as Figure 11, but for the model L-40-45.

Standard image High-resolution image
Figure 13. Refer to the following caption and surrounding text.

Figure 13. Same as Figure 11, but for the model L-40-40.

Standard image High-resolution image

On the other hand, the production curves of the three models are distinct because of the different fission fragment asymmetries. In the model L-40-50, abundant production of elements near the third r-process peak (A ≈ 195), comparable to the solar abundance observation, is found. Lanthanides (Z = 57–71) are insufficiently synthesized, especially for those before gadolinium (Z = 64). Due to the symmetric fission fragment approximation applied, elements with Z < 45 are not produced after the spontaneous fission of heavy nuclei. The end products of nucleosynthesis are lacking in elements before the second r-process peak (A ≈ 130). In the model L-40-45, the solar abundance distribution is well recovered from the second r-process peak to the third r-process peak. Moreover, an excessive production of elements at the third r-process peak is obtained. As for the production curve of the model L-40-40, elements slightly less massive than the second r-process peak are present. However, elements on the two sides of the second r-process peak are significantly underproduced when compared to the solar abundance distribution. No elements around the first r-process peak (A ≈ 80) are formed in the three models.

Furthermore, the production curves of the model L-30-45 are shown in Figure 14 to illustrate the robustness of the results. The general features of the production curves in the models L-30-45 and L-40-45 are very similar. The nucleosynthesis is thus generic against the variation in ΔS among different models.

Figure 14. Refer to the following caption and surrounding text.

Figure 14. Same as Figure 12, but for the model L-30-45. The production curves are very similar to those shown in Figure 12 in general.

Standard image High-resolution image

4. Conclusions

We present a study of the explosion of subminimal neutron stars in close neutron star binaries. Based on the analysis of mass transfer stability and timescale in asymmetric neutron star binaries (Clark & Eardley 1977; Blinnikov et al. 1984), we construct the initial subminimal neutron star model with the minimum mass allowed by the EoS at equilibrium. The hydrodynamic instability is initiated by removing a portion of mass from its surface, following the methodology adopted by Colpi et al. (1993) and Sumiyoshi et al. (1998). A modern EoS and a detailed nuclear reaction network are coupled with the Euler equations to study the evolution of the neutron star. The major factor responsible for the instability development is leptonization by the β-decays of nuclei, which results in a higher electron degeneracy pressure of the neutron star matter. The heating effect following the β-decays serves as a supplementary factor that further enhances the expansion of neutron star crust-like matter. The overall thermal energy generated through thermonuclear reactions and spontaneous fission is rather minor in comparison. The net cumulative thermal energy injected into the neutron star matter from the network is in the order of 1050 erg, which is only ∼0.1%–1% of its rest mass. While nuclear energy is a negligible contribution of energy in other astronomical simulations of typical neutron stars, it plays a significant role in a subminimal neutron star, which is just marginally stable.

The peak luminosity of electron antineutrino emission due to β-decays of nuclei predicted by our models is 1–2 order lower than the values reported by Colpi et al. (1993) and Sumiyoshi et al. (1998). A maximized choice of the coefficient of β-decay rate estimation formula with uncertainty (Lattimer et al. 1977) used in their studies is responsible for the difference, and Sumiyoshi et al. (1998) have already pointed out that the β-decay rate is probably overestimated in their simulations. Adopting realistic β-decay half-lives of nuclei, our calculations predict a neutrino burst with a lower peak in luminosity but lasting for a longer duration. The ejecta is powered by the nuclear reactions, and thus the temperature stays at ∼109 K despite adiabatic cooling. At the end of the simulation, the radius of the outermost mass element grows to ∼104 km. Radiative cooling may become effective when the surface area of the expanding ejecta is large, and the temperature should decline more quickly. Further investigation is required to predict whether a soft gamma-ray burst, with high enough luminosity to be observed from the earth, is associated with the explosion of the subminimal neutron star. A bolometric light curve powered by the nuclear reactions may be obtained for further research.

The explosion of subminimal neutron star releases ∼0.05 M of ejecta in total. The r-process occurs in all mass elements, according to our calculations. The excessive production of elements near the third r-process peak is observed, especially in the models L-40-45 and L-40-40. We attribute the synthesis of heavy elements of this event to the initially low electron fraction and the relatively low temperature of the ejecta. The two conditions allow the nucleosynthesis to operate under a high neutron excess environment for a sufficiently long duration. A similar scenario is also realized in neutron star mergers with high neutron excess. The observation of lanthanides and elements at the third r-process peak associated with a binary neutron star merger by Tanvir et al. (2017) suggests the crucial contribution from neutron star ejecta to the heavy elements production in the Universe. The neutron star merger as a promising astronomical r-process site is under extensive studies (e.g., Freiburghaus et al. 1999; Bovard et al. 2017; Côté et al. 2018; Kullmann et al. 2023).

Meanwhile, several recent researches on binary neutron star mergers and black hole-neutron star mergers suggested that cold dynamical ejecta, originating from the outer regions of neutron stars, may abundantly produce heavy elements up to the third r-process peak (e.g., Cowan et al. 2021). Our simulation results basically agree with the prediction that cold ejecta from neutron stars is favorable for the formation of elements near the third r-process peak (Korobkin et al. 2012). Furthermore, the recent investigation on the nucleosynthesis of subminimal neutron star explosion by Panov & Yudin (2020) found that heavy elements up to the third r-process peak can be excessively formed in the ejecta from the inner crust initially. We notice that such a feature is also present in our results regarding ejecta originating from deeper regions of the subminimal neutron star. While the nucleosynthesis calculation in their work is performed without coupling to the hydrodynamic, we may realize that the robustness of the r-process during the subminimal neutron star explosion is insensitive to the hydrodynamic evolution of the ejecta with sufficiently low initial electron fraction. The total dynamical ejecta mass from binary neutron star mergers is typically in the order of 10−3–10−2 M (Bovard et al. 2017; Radice et al. 2018; Kasliwal et al. 2022), less than what we have obtained from subminimal neutron star explosions. Therefore, the latter should also contribute to the heavy element abundances if their event rate is comparable to that of binary neutron star mergers.

We demonstrate that the fission fragment calculation can significantly affect the nucleosynthesis results. The fragment asymmetry parameter γ is introduced to mimic the fission fragment distribution. The study by Hao et al. (2022) implied that the models with γ = 0.45 or 0.40 should give more reliable results than those with symmetric fission fragments, such as in the model L-40-50, as well as previous studies by Colpi et al. (1993) and Sumiyoshi et al. (1998). Further study may be conducted by determining the fission fragment distribution for each isotope individually with more sophisticated methods (Mumpower et al. 2020; Vassh et al. 2020; Hao et al. 2022), which should improve the accuracy of the predicted heavy elements production. Moreover, neutron-induced (Panov et al. 2010) and β-delayed fission (Panov et al. 2016) may be included. Nevertheless, the fission properties, as well as other nuclear reactions involving neutron-rich and heavy nuclei, are not well determined experimentally. Theoretical modeling must be used in r-process simulation studies, and thus systematic uncertainties always remain in these nucleosynthesis calculations, including the results presented in this paper.

A crucial assumption of our study is the existence of peculiarly light neutron stars with mass ∼0.1 M in close neutron star binaries. Even though such a configuration is allowed by modern EoS, a neutron star with such a small mass has never been observed (Suwa et al. 2018). Nevertheless, recent research on asymmetric neutron star binaries (e.g., Ferdman et al. 2020; Vincent et al. 2020) opens the door to the dynamical formation of low-mass neutron stars through mass transfer. In our study, a large fraction of mass (ΔS ≥ 0.30) is removed from the surface of the minimum-mass neutron star initially. Whether this treatment is realistic in describing the tidal stripping of the low-mass neutron star is unclear. We find that the delayed explosion of the subminimal neutron star is considerably postponed if the initial mass removal ratio ΔS is lower than 0.40. Although the hydrodynamic evolution of the star after losing stability is rather insensitive to the exact value of ΔS, it experiences more periods of radial oscillation during the quasi-equilibrium state if ΔS is smaller. In addition, the temperature of the neutron star crust-like matter temporarily exceeds 1010 K whenever the compression wave arrives. Neutrino capture and positron capture may additionally give rise to leptonization alongside β-decays of nuclei. Whether the electron fraction remains as low as the value found in our study has to be justified. On the other hand, a small value of ΔS may contradict the conditions described by Clark & Eardley (1977) and Blinnikov et al. (1984) that the tidal disruption of the whole subminimal neutron star happens at most within seconds. While the massive companion is assumed to be stable during the mass transfer in this paper, the influence of the massive companion must be considered if the instability of the subminimal neutron star takes a long duration to develop.

A robust r-process nucleosynthesis is realized in the ejecta from the explosion of a subminimal neutron star. Lanthanides and heavy elements near the second and third r-process peaks are synthesized as end products of the nucleosynthesis, with similar mass fractions as the solar abundances.

Acknowledgments

We acknowledge F. X. Timmes for making the nuclear reaction subroutine Torch open-source and the explanation about its usage. We thank A. S. Schneider for making the EoS publicly available. We also thank P. Möller and K. P. Santhosh for providing publicly available nuclear reaction rates. This work is partially supported by a grant from the Research Grant Council of the Hong Kong Special Administrative Region, China (Project Nos. 14300320 and 14304322). This material is based upon work supported by the National Science Foundation under Grant AST-2316807.

Appendix A: Models with Smaller ΔS

The hydrodynamic evolution of the model L-40-50 is discussed in detail in Section 3.1. As mentioned previously, the explosion is further postponed for models with lower ΔS. Consequently, the temperature reached by the mass elements can be even higher during radial oscillations prior to the explosion. The hydrodynamic simulation results of the models L-36-50 and L-32-50 are shown in Figures 15 and 16, respectively, to illustrate the influence of a smaller ΔS in the hydrodynamic simulations.

Figure 15. Refer to the following caption and surrounding text.

Figure 15. Same as Figure 5, but for the model L-36-50.

Standard image High-resolution image
Figure 16. Refer to the following caption and surrounding text.

Figure 16. Same as Figure 5, but for the model L-32-50.

Standard image High-resolution image

The evolution of the central density, electron antineutrino luminosity, and net cumulative thermal energy deposited on the neutron star by the network in the models L-40-50, L-38-50, L-36-50, L-34-50, L-32-50, and L-30-50 are shown in Figures 17 and 18. As a consequence of the difference in explosion timescale discussed in Section 3.1, the neutrino luminosity of models with lower ΔS fluctuates periodically following the radial oscillations of the star prior to the explosion. A peak luminosity of ∼3 × 1050 erg s−1 is always found in these models when the explosion occurs. The net cumulative thermal energy deposited on the star by the network is similar, starting from ∼0.3 s after the occurrence of explosion in all these models.

Figure 17. Refer to the following caption and surrounding text.

Figure 17. Evolution of central densities in the models L-40-50, L-38-50, L-36-50, L-34-50, L-32-50, and L-30-50. The nuclear reaction network is updated for every δ tnet = 1 × 10−4 s (solid lines), 2 × 10−4 s (dashed lines), and 5 × 10−4 s (dotted lines) before the explosions.

Standard image High-resolution image
Figure 18. Refer to the following caption and surrounding text.

Figure 18. Same as Figure 8, but for the models L-40-50, L-38-50, L-36-50, L-34-50, L-32-50, and L-30-50. The nuclear reaction network is updated for every δ tnet = 1 × 10−4 s (solid lines), 2 × 10−4 s (dashed lines), and 5 × 10−4 s (dotted lines) before the explosions.

Standard image High-resolution image

Appendix B: Numerical Test

Owing to limited computational resources, we do not solve the nuclear reaction network for every time step determined by the CFL condition during the hydrodynamic simulations. The extra energy source terms in Equation (4) are updated by solving the network with time steps δ tnet shorter than the β-decay half-lives ∼10−3 s of the neutron-rich nuclei presented. In Figures 17 and 18, we update the network for every δ tnet = 1 × 10−4 s (solid lines), 2 × 10−4 s (dashed lines), and 5 × 10−4 s (dotted lines) before the explosions. We notice that the explosions are advanced in all models whenever a smaller δ tnet is picked as we approach the physical limit δ tnet = 0. In particular, the explosion timescale as well as the peak in neutrino emission have not fully converged even for the finest δ tnet adopted.

Despite the difference in the hydrodynamics, the overall and local net energy generated by the network are insensitive to the choice of δ tnet. Moreover, the features of chemical element productions discussed in Section 3.3 are robust against the value of δ tnet chosen. The analysis performed in this paper is based on the simulation results using δ tnet = 1 × 10−4 s. Further investigation is required to study the exact explosion time of subminimal mass neutron star explosions.

Footnotes

  • 3  
  • 4  
  • 5  

    The extrapolation starts when the density of the mass elements reaches ∼104 g cm−3. To trace the evolution of all mass elements until the extrapolation can be started, we remove the outermost mass element whenever its density reaches near the minimum density available and continue the hydrodynamic simulation of the remaining mass elements up to about 1–1.5 s.

  • 6  

    For the EoS we adopt, there is a first-order phase transition from uniform to nonuniform nuclear matter configurations at density near 6–8 × 1013 g cm−3. During the EoS construction, the phase with lower free energy is chosen if there are solutions for both uniform and nonuniform phases. In some rare cases, the phase with slightly higher free energy is chosen to avoid an unphysical adiabatic index. Hence, the magnitude of temperature rise in our simulations is uncertain because the EoS near the phase transition point is rather ambiguous at low temperatures. For more details, please refer to the original literature on the EoS construction (Schneider et al. 2017).

  • 7  

    We have verified that the pressure and the specific internal energy density are insensitive to temperatures below ∼1010 K at densities above ρth, so the hydrodynamics is also unaffected.

10.3847/1538-4357/acf570
undefined