A Nuclear Equation of State Inferred from Stellar r-process Abundances

This section describes the DNS populations we use for this study, which will combine with the EOS input to determine specific NSM yields. First, we use a theoretical formation model (population synthesis) to estimate the past history of NSMs in our Galaxy. We will apply this DNS model to two cases in this work. Then, we consider the case for the current DNS population, with the underlying assumption that the r-process producing NSMs were of systems similar to present-day DNS systems in the Milky Way.

The continuously updated StarTrack code estimates the evolution of single and interacting binary stars, both individually and as populations (Belczynski et al. 2008). Frequently applied to interpret (and calibrated against) astronomical observations of several types—including binary pulsars and gravitational-wave observations—this code provides a useful benchmark to explore plausible self-consistent models for the Galactic binary population. After reviewing several recent studies (Dominik et al. 2012, 2013, 2015; Belczynski et al. 2016a, 2020, 2016b; Wysocki et al. 2018; Drozda et al. 2020), we choose the M15 model (submodel B) in Belczynski et al. (2016b) with strong pair-instability (pulsation) SNe and modest NS natal kicks (σ = 130 km s⁻¹; Wysocki et al. 2018). This model successfully reproduces compact object merger rates from the third observing run of the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo collaboration, remains qualitatively consistent with known constraints on the maximum NS mass (M_TOV), and is consistent with (most) observed Galactic pulsar masses.

We translate the M15 simulation (of a fiducial cosmological volume, as in Dominik et al. 2015; Belczynski et al. 2016a) to our application of the Milky Way as follows. Since we are interested only in the distribution and not the (highly uncertain) overall scale factor, we disregard overall normalization. In addition, because DNSs typically merge fairly quickly, we can ignore small differences between cosmological and galactic star formation history. We can therefore identify the relative likelihood of NSMs by the weighted data (M_1,i , M_2,i , s_i ), where s_i represents the star formation, metallicity, and population-weighted likelihood of the M_1,i −M_2,i mass pair, as described in previous work by Dominik et al. (2015). For computational efficiency, we keep only those mass pairs that are in the top 99% of the most frequent merging systems in the population synthesis model, i.e., those with ${s}_{i}\geqslant 0.01\times \max ({\vec{s}}_{i})$. Of several thousand unique mass combinations, these top 99% account for roughly 50% of all unique merging DNS systems in the population synthesis model output. Although we are potentially ignoring cases that might eject significant amounts of r-process material, our choice to work in abundance ratios rather than absolute abundances effectively removes the sensitivity to total ejecta mass. At most, we are removing cases that would contribute at the 1% level in Sr/Dy−Th/Dy abundance space, at the benefit of doubling our computational efficiency.

Figure 4 shows the masses of the merging DNSs in the StarTrack model in purple, weighted by their likelihood (s_i ). The most likely merging system has M₁ ≈ 1.35 M_⊙ and M₂ ≈ 1.1 M_⊙. Note that there are no systems in this model with M_1,2 > 2.0 M_⊙, and none within the first mass gap, as a consequence of the model choices in Belczynski et al. (2016b) and Wysocki et al. (2018). We will test two cases using this DNS population: one with a constraint on the maximum mass of a nonrotating NS ("with M_TOV") and one without such a constraint ("no M_TOV"). These two cases will be discussed in more detail in Section 3.2.

Zoom In Zoom Out Reset image size

Next, instead of using a synthetic DNS population, we can also use the observed mass distribution of NSs in the Milky Way. As a third, alternative, case we reconstruct the EOS informed by inferences about the mass distribution of individual Galactic NSs. Specifically, we assume the primary NS mass is drawn from a mixture of two normal distributions, with a primary peak at 1.34 M_⊙ and a secondary at 1.78 M_⊙ (from Alsing et al. 2018; see also Özel et al. 2012). We draw the secondary mass based on the mass-ratio distribution from Farrow et al. (2019): essentially a Gaussian peaked at q = 1, where q > 0.69 with 99% confidence. We choose 500 random selections from these distributions to represent the primary and secondary NS masses in the binary system. This DNS population is equivalent to assuming that the NSM progenitors of the r-process material in the Galaxy were of DNS systems similar to existing DNS binaries currently known in the Milky Way. The primary and secondary masses for this empirical distribution are also shown in Figure 4 in yellow. On average, this hypothesized Galactic DNS distribution predicts a much broader total mass distribution for DNSs than the observed distribution of known Galactic DNSs (see, e.g., Farrow et al. 2019), but extends to total masses more consistent with GW190425 (Abbott et al. 2020). Except for the input DNS population, this third case ("Alsing+") will use the same method and constraints as the "with M_TOV" case that uses the M15 model for the DNS masses.

Considering that the metal-poor stars likely originated in ultrafaint dwarf galaxies like Reticulum II and Tucana III, the correlation between [Fe/H] and time becomes unclear and prevents imposing a delay-time distribution on these merger pairs. Rather, we instead consider the cumulative effect that the entire population of NSMs would have had on the r-process abundance ratios.

Within the framework of this paper, the EOS bridges the DNS masses with theoretical descriptions of NSM yields and, therefore, with a total abundance distribution that can be compared to observations of metal-poor stars. Since the EOS is the unknown parameter, we will use an MCMC method to build an EOS from the abundance observations themselves. First, let us discuss a basis for building an EOS. To smoothly explore an EOS parameter space, we consider a four-parameter spectral-decomposition representation, described by Lindblom (2010) and implemented in LALSuite (LIGO Scientific Collaboration 2018). The energy density as a function of the pressure, (p), is described by the adiabatic index, Γ:

$$\begin{eqnarray}&&{\rm{\Gamma }}(p)=\displaystyle \frac{\epsilon +p}{p}\displaystyle \frac{{dp}}{d\epsilon }.\end{eqnarray} \tag{ 1 }$$

In the spectral-decomposition representation, the EOS can be expressed as an expansion of the adiabatic index:

$$\begin{eqnarray}&&{\rm{\Gamma }}(x)=\exp \left(\sum _{k}{\gamma }_{k}{x}^{k}\right),\end{eqnarray} \tag{ 2 }$$

where $x=\mathrm{log}(p/{p}_{0})$, and p₀ is the minimum pressure. In the LALSuite implementation, k goes from 0 to 3, and the EOS is therefore defined by four γ_k values. With some constraints on these parameters, arbitrary EOSs and NS mass–radius curves can be constructed. LALSuite provides a robust set of built-in functions to construct a neutron star EOS from $\vec{\gamma }$ and find properties such as the radius and tidal deformability as functions of the NS masses. Our method utilizes this capability of LALSuite to smoothly explore the $\vec{\gamma }$ parameter space and find the EOS-dependent properties needed by the ejecta yield descriptions.

There are several constraints from both dense-matter theory and NS observations available to guide the EOS parameter choices. In addition, not all allowed combinations of $\vec{\gamma }$ lead to physical EOSs. First, rather than consider the entire $\vec{\gamma }$ parameter range nominally allowed by LALSuite, we use the transformation and corresponding limits from Appendix B of Wysocki et al. (2020) to explore only those combinations that are physically allowed. These transformed parameters will be referred to as ${\vec{\gamma }}^{{\prime} }$.

Next, we fold in constraints from observation and theory. There are no direct constraints on the $\vec{\gamma }$ parameters since these are merely the representations to build the EOS (but see Jiang et al. 2020). However, solving the TOV equations with the EOS parameters does allow direct comparison with astronomical observables. Of primary importance is the maximum nonrotating NS mass, M_TOV. We use the updated mass of PSR J0740+6620 as a lower limit that M_TOV must satisfy for the EOS (Fonseca et al. 2021; Riley et al. 2021). This limit is implemented by adopting a cumulative distribution function (CDF) of a normal distribution (${ \mathcal N }$) with the centroid and 1σ spread from the measurements of PSR J0740+6620, namely, $\mathrm{CDF}(x,{ \mathcal N }(\mu ,\sigma ))$, where μ = 2.072 and σ = 0.07. This CDF enforces M_TOV of at least 2.072 M_⊙—with some room for uncertainty at the lower end—but imposes no upper limit. An upper limit on M_TOV has been suggested from the remnant mass of the GW170817 merger event based on inferences that the remnant probably collapsed somewhat quickly (≲0.1 s) into a black hole (e.g., Margalit & Metzger 2017; Shibata et al. 2019; Kawaguchi et al. 2021). One of our model cases enforces M_TOV at a somewhat higher 2.30 ± 0.15 M_⊙ limit ("with M_TOV") to include this observational constraint. As a second case ("no M_TOV"), we release the upper limit on M_TOV and allow the model to explore solutions with arbitrarily high NS masses, accommodating the possibility that the GW170817 remnant did survive longer.

On the microphysics side, both theory and experiment offer some constraints on the EOS. These constraints place limits on, e.g., the symmetry pressure, L, and the nuclear incompressibility, K_∞, at saturation densities (ρ₀ ≈ 2.8 × 10¹⁷ kg m⁻³, or n₀ ≈ 0.17 fm⁻³). Tight constraints on K_∞ exist from experiment; however, these values are defined for symmetric nuclear matter, whereas the EOSs defined with LALSuite assume neutron star material in β-equilibrium. Therefore, we implement limits on the somewhat more loosely constrained, but more analytically straightforward, parameter L. For an EOS, the pressure can be expressed as

$$\begin{eqnarray}&&p(n)={n}^{2}\displaystyle \frac{d\epsilon }{{dn}},\end{eqnarray} \tag{ 3 }$$

where n is the nucleon number density instead of the mass density, ρ. On a plot of energy density per nucleon versus number density, L is proportional to the slope at nuclear saturation:

$$\begin{eqnarray}&&L=3{n}_{0}{\left.\displaystyle \frac{d\epsilon }{{dn}}\right|}_{{n}_{0}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{3}{{n}_{0}}p({n}_{0}).\end{eqnarray} \tag{ 5 }$$

In this form, L can be straightforwardly evaluated with the functions available in LALSuite. Energy density versus density curves can be calculated directly by LALSuite; however, we caution that these values are the volumetric energy and mass densities, respectively. Therefore, we take the pressure-versus-mass densities from LALSuite and convert the mass density to a number density by dividing by the nucleon mass plus the typical binding energy per nucleon (≈−16 MeV). Note that this approximation is only valid around nuclear saturation. In addition, this transformation to obtain L may be different from the "true" value by a few MeV. In this way, we can obtain L for an EOS using the functions available in LALSuite. When evaluating the likelihood of a particular EOS (discussed in Section 4), we calculate L from Equation (5) and prefer EOSs within a liberal range of 30 MeV ≲ L ≲ 100 MeV, as determined by theory and experiment (e.g., Dutra et al. 2012; Lattimer 2012; Tews et al. 2019; Piekarewicz 2021). These limits are imposed as a likelihood in the model, which will be discussed in Section 4.

With a population of merging DNSs in hand (two cases to be tested with the StarTrack M15 model and one with the Alsing et al. 2018 distribution), we can compute the total, EOS-dependent r-process yields for each merger in our DNS list. We use Equation (6) from Krüger & Foucart (2020) to quantify the NSM dynamical ejecta and Equation (S4) from Dietrich et al. (2020) for the disk mass. The amount of mass that is lost through winds (neutrino-driven or viscous) from the accretion disk is assumed to be a constant fraction of 40% of the total disk mass as in, e.g., Radice et al. (2018). (We will discuss the implications of this assumption in Section 5.) In addition to depending on the NS masses, these fit choices also depend directly on physical properties of NSs such as the radius, compactness, and tidal deformability. Since these properties are EOS-dependent, the ejecta fits can be expressed as functions of the NS masses and the EOS, so the net ejecta can be expressed as ${m}_{\mathrm{ej}}({M}_{1},{M}_{2},\vec{\gamma })={m}_{\mathrm{dyn}}+0.4\,{m}_{\mathrm{disk}}$, where $\vec{\gamma }$ represents the EOS parameters, described in Section 3.2.

The available fits to NSM ejecta come with significant uncertainties, mostly as a result of fitting across multiple groups and simulations that may differ in their implementations. To account for these uncertainties, we allow the computed ejecta masses and lifetimes to vary from (truncated) normal distributions. Each distribution is centered at the expected value given by each equation and has a spread corresponding to the uncertainty of the fit (for details, see Radice et al. 2018; Krüger & Foucart 2020; Lucca & Sagunski 2020). For each DNS merger pair in the population synthesis output, we draw N = 200 samples from this truncated-normal distribution such that each NSM gives some statistical spread of masses rather than one discrete value; i.e., in the expressions above, each m_dyn and m_disk is a set of length N that accounts for uncertainty in each corresponding relationship. We then multiply each of these ejecta masses by the appropriate trajectory-dependent nucleosynthesis output to produce overall outflow yields. We also allow for a similar spread in the corresponding nucleosynthetic compositions produced by the dynamical and post-merger outflows to obtain a distribution of outflow yields. All nucleosynthesis output is based on previous calculations in Holmbeck et al. (2019) and Holmbeck et al. (2021a) that use the FRDM2012 nuclear mass model (Möller et al. 2012) for theoretically calculated nuclear reaction and decay rates.

For the dynamical ejecta composition, we use the thermodynamic evolution of a trajectory from the 1.4–1.4 M_⊙ NSM simulations of S. Rosswog (Piran et al. 2013; Rosswog et al. 2013) as in Korobkin et al. (2012). The initial composition may vary based on the EOS and masses of the merging NSs, so we choose final abundances from a nucleosynthesis network calculation that starts with a merger-dependent Y_e as in Holmbeck et al. (2019), where the initial Y_e is calculated using Equation (10) from Nedora et al. (2021). Since a spread in Y_e may be expected in the dynamical ejecta, we apply some variation on the initial Y_e . The Gaussian is centered at the calculated Y_e and given a 1σ spread of 0.01, motivated by the values in Table 1 of Nedora et al. (2021). Therefore, we calculate the final abundances for a range of initial Y_e values and combine those abundances with weights from a normal distribution. We draw the same number of samples (N = 200) as for the ejecta masses from this normal distribution to obtain a set of abundances, A _dyn, of dimension N × 3, where there are N yields for each of the three elements we study here: Sr, Dy, and Th.

Similarly, the composition of the post-merger outflows may be influenced by the lifetime (τ) of the remnant massive NS before it collapses into a BH (Lippuner et al. 2017). We draw random samples of the lifetime of the merger remnant based on the uncertainty in the lifetime fit from Lucca & Sagunski (2020). Then, lifetime-dependent r-process yields from Holmbeck et al. (2021a) are used to find the ultimate merger-dependent disk wind abundances, A _disk.

Finally, these random samples drawn from Y_e and τ for the dynamical and disk compositions are multiplied with their corresponding ejecta masses and added to obtain the final yields:

$$\begin{eqnarray*}&&{\boldsymbol{Y}}=\left({{\boldsymbol{A}}}_{\mathrm{dyn}}{m}_{\mathrm{dyn}}\right)+0.4\left({{\boldsymbol{A}}}_{\mathrm{disk}}{m}_{\mathrm{disk}}\right).\end{eqnarray*}$$

Note that Y also has dimensions of N × 3. For each merging DNS pair, the Sr/Dy and Th/Dy abundance ratios are computed by an element-wise division of the corresponding columns in Y . The final N-sized set of two abundance ratios is added to a two-dimensional distribution in abundance space as in Figure 2, after being weighted by the likelihood factor, s_i , for that merging pair. (For the Alsing+ case, s_i = 1.) Including these systematic uncertainties of the ejecta masses, the remnant lifetime, and the dynamical ejecta Y_e significantly increases computation time, but it allows flexibility in the model that would otherwise be too highly constrained by the literature fits to the simulations, especially as those fits update with new available data.

For the two cases that use the M15 model for the DNS population, Figure 5 shows some agreement with the EOS derived for GW170817. The two sets of posterior samples are tightly constrained at low densities (perhaps due to the constraint on L), but diverge at higher densities. This divergence is unsurprising since NSs are reasonable probes of the low-density EOS, and only the most massive NSs offer insight into densities greater than 6ρ₀ (Lattimer 2012).

The mass–radius curves derived for these EOSs in Figure 6 are slightly softer than the NICER/XMM-Newton median value, though there is still reasonable agreement within 2σ of their derived radii. Including agreement with the NICER/XMM-Newton in the likelihood function does not significantly change these results. Both M15 cases lead to EOS solutions narrowly grouped in mass–radius at M ∼ 1.35 M_⊙, most likely do the high frequency of DNSs with M₁ = 1.35 M_⊙ in the M15 model (see Figure 4). However, for the "no M_TOV" case, the posterior solutions sometimes violate causality, as indicated by the pink curves extending into the gray region in Figure 6. The two families of solutions are nearly identical in the ejecta parameters they predict, indicating some degeneracy when the spectral-decomposition parameters are propagated to ejecta and abundance observables through the NS mass–radius relationships. For the most common case in the StarTrack model of a 1.35–1.1 M_⊙ merger, both cases predict a wind mass of about 9 × 10⁻² M_⊙ and a much smaller dynamical mass of about 6.5 × 10⁻³ M_⊙. The two cases also give the same dynamical ejecta Y_e of about 0.17. Therefore, a 1.35–1.10 M_⊙ merger would likely produce the same kilonova signature with either of the two maximum-likelihood solutions in the with-/no-M_TOV cases. For a 1.40–1.35 M_⊙ merger, the predicted wind mass differs by about 20% between the two solutions, which could leave a distinguishable signature in the associated kilonova.

Figure 7 shows the Th/Dy and Sr/Dy abundances and the difference between the model and observations produced by the maximum-likelihood solution for the "with M_TOV" case. The model reasonably reproduces the most prominent feature in the abundance space of Figure 2 occurring at Sr/Dy ≈ 1.25 and Th/Dy ≈ −1.15, though the peak Th/Dy is somewhat lower. A couple of distinguishable artifacts also appear in the output abundances: first, a "spur" at high Sr/Dy. This feature commonly occurs for quite low-mass, low-asymmetry mergers. Because of their low total mass, the wind outflows from the long-lived remnant accretion disk drive the Sr abundances to high values. That said, the dynamical masses for mergers with M₁ ∼ 1.2 M_⊙ and q ≳ 0.9 will be small, but not zero. Although the dynamical mass may be very low, actinide production can still be sufficiently high to contribute a significant fraction of the total Th that is ejected. This spur therefore represents a limiting case in which low asymmetries and low masses drive not only high Sr production, but also a low-mass tidal tail that is rich in actinides.

Zoom In Zoom Out Reset image size

Second, with the ejecta equations employed here, there appears to be a diagonal floor in which the ejecta cannot simultaneously be rich in Sr/Dy and deficient in Th/Dy. Conditions producing this lower edge are ones with q ≳ 0.9 and total masses M₁ + M₂ ≈ 2.6–2.9. At the very lowest Th/Dy are mergers with M_1,2 ∼ 1.45 M_⊙. Conditions close to the solar abundances are typically produced by mergers with M_1,2 ∼ 1.3 M_⊙. These conditions produce long-lived remnants that have the same wind compositions and relatively small wind ejecta masses. Going from the upper right of the edge to the lower left (i.e., increasing the total binary mass), the contribution by the dynamical ejecta is increased. In nearly symmetric cases, low tidal deformabilities lead to dynamical 〈Y_e 〉 values that struggle to produce significant amounts of Th, but have a high yield of Dy. Their low-Th and high-Dy abundances on average bring the total Th/Dy and Sr/Dy yields down from pure-wind abundances as the NS masses increase. Therefore, this edge can be thought of as the limiting wind composition case with an increasing amount of moderately low-Y_e dynamical ejecta contributing to the total outflows. Both the high-Sr/Dy spur and the low-Th/Dy floor indicate limits in the computational method from both the compositions and the total ejecta masses. However, it is worth noting that the extension of the "observations" into the low-Th/Dy, high-Sr/Dy region—below the diagonal floor—are not populated by an actual observational measurement in a metal-poor star. Rather, recall the distribution was given a random spread to account for observations that statistically could have these particular abundance combinations in our attempt to build a complete sample by removing some observational bias. Therefore, it may be possible that no such combination exists in metal-poor stars that have primarily r-process origins of their heavy element material. More observations of metal-poor stars (currently underway, e.g., through RPA efforts) could shed light on the reality of this limit.

The abundances show some sensitivity to the percentage of the total disk mass that escapes from the system (assumed to be 40% here). Since most of the Sr that is produced in our model is from the wind component, altering this disk-to-wind fraction has the largest effect on the Sr/Dy abundances. Increasing (decreasing) this fraction effectively moves the output abundances to the right (left) on Figure 7, up to the diagonal abundance floor. In order to compensate for the surplus (deficiency) of Sr from increasing (decreasing) the ejecta fraction, the EOS must lead the merger to eject relatively less (more) Sr compared to Dy. Changing the Sr/Dy ratio can be achieved most notably by changing the NS compactness. Soft EOSs with relatively smaller NS radii will eject more dynamical ejecta for q ≳ 0.9, moderately low-mass systems, such as those in the DNS populations that use the M15 population synthesis model. With more Dy, the large Sr/Dy values that are achieved by higher disk-to-wind fractions are effectively reduced. Therefore, we expect that increasing (decreasing) the disk-to-wind fraction will lead to inferred EOSs that are more (less) compact. However, exploring the sensitivity of our model to this fraction would require a separate study, and for now we recognize it as an area for improvement.

There is one feature that is virtually absent in the model output: the lobe at moderately high Th and low Sr. To produce low Sr, the dynamical ejecta must be very high. However, NS−NS mergers with massive dynamical ejecta require high asymmetries and masses. Such massive, asymmetric cases will correspondingly produce a high Sr abundance from the disk outflows as well. In other words, in this model, a massive dynamical ejecta component will always come with a similarly massive (or even more massive) wind component. Therefore, this low-Sr/Dy region could possibly be populated if there were very little to no wind ejecta at all. An r-process source that may realize this case—and one that this model neglects—is an NSBH merger. Theoretically, disk wind outflows would not exist in the NSBH merger case, and the mass outflows would be entirely dynamical: perhaps rich in Th and poor in lighter r-process elements like Sr. Therefore, it is possible that the most Sr-poor abundances that could not be reproduced by the NS−NS mergers in this work could have NSBH origins instead. The StarTrack output also provides these NSBH binaries, and we reserve computing and including their nucleosynthetic outflows for future work.

The two cases that use the M15 model for the DNS population both produce reasonable agreement with all available data: the EOS from GW170817, the NS radius and mass measurements from NICER/XMM-Newton, and the r-process abundances of metal-poor stars. Here, we look at how the (Alsing et al. 2018) DNS distribution fares with our MCMC model to see if we can obtain additional agreement with the current distribution of DNS systems in the Galaxy.

While the two M15 cases show reasonable behavior for their posterior EOS samples, the Alsing+ case rather proceeds to very extreme values for the EOS parameters, which can be seen by the large fluctuations in the slope of the EOS in Figure 5. To reconcile the Alsing et al. (2018) DNS distribution with the observed abundances, a very extreme EOS is needed, perhaps even beyond the limits and constraints supplied to the model. For this case, there is no agreement with the GW170817 EOS at low densities. On the other hand, the posterior mass–radius curves for the Alsing+ case agree rather well with the NICER/XMM-Newton values, preferring a stiffer EOS than the cases that use the M15 model. Finally, Figure 8 shows the resulting abundances for the best-fit EOS case using the Alsing et al. (2018) NS mass distribution. There is nearly no agreement with the r-process abundance data, and the results group tightly around high Th/Dy and Sr/Dy values. This particular feature consists of the symmetric mergers with M_1,2 ≈ 1.4, the main peak in the Alsing et al. (2018) distribution in Figure 4. The "wisps" of abundances around the data result from the slightly asymmetric cases that appear in the Alsing et al. (2018) distribution.

Zoom In Zoom Out Reset image size

The clear failure of this model indicates that this estimate of the present-day NS mass distribution (and the associated assumption that DNSs can form with primaries drawn from this distribution) cannot be reconciled with the r-process abundances in metal-poor stars, under all of the other constraints and assumptions supplied to the model. This conclusion is in contrast to earlier work (Holmbeck et al. 2021a), which claims agreement between r-process stars and the Galactic DNS distribution (e.g., for the DD2 EOS). The largest source of difference between that work and the present model might possibly be the description of the dynamical ejecta Y_e . In this model, we employ a relationship between the Y_e of the dynamical ejecta and the NS properties, including the EOS (i.e., from Equation (10) of Nedora et al. 2021). Since the results are sensitive to this description, a robust and systematic study of dynamical ejecta Y_e as a function of NS properties would greatly benefit this type of work.