Thermal Effects in Binary Neutron Star Mergers

We study the impact of finite-temperature effects in numerical-relativity simulations of binary neutron star mergers with microphysical equations of state and neutrino transport in which we vary the effective nucleon masses in a controlled way. We find that, as the specific heat is increased, the merger remnants become colder and more compact due to the reduced thermal pressure support. Using a full Bayesian analysis, we demonstrate that this effect will be measurable in the postmerger gravitational wave signal with next-generation observatories at signal-to-noise ratios of 15.


INTRODUCTION.
The extreme conditions found in neutron stars make them an ideal means for probing the nuclear equation of state (EOS). Electromagnetic (EM) observations of pulsars have provided valuable information about the mass distribution of neutron stars (Lattimer & Prakash 2001;Alsing et al. 2018), and recent results from NICER offer constraints on their radii (Riley et al. 2021;Ludlam et al. 2022;Salmi et al. 2022). Binary neutron star (BNS) mergers give additional astronomical constraints; the gravitational waves (GW) and EM counterpart of GW170817 contained details about the EOS via tidal deformability measurements and ejecta characteristics (Abbott et al. 2017a;Margalit & Metzger 2017;Radice et al. 2018b;De et al. 2018;Abbott et al. 2018a;Dietrich et al. 2020;Kashyap et al. 2022).
Due to the high Fermi temperature of matter in a neutron star, constraints obtained from pulsars and BNS inspirals are informative of the zero-temperature equation of state (EOS). On the other hand, temperatures as high as 100 MeV might be reached in the post-merger phase (Perego et al. 2019), making it a possible probe of the jmf6719@psu.edu * Alfred P. Sloan fellow finite-temperature EOS. Current GW detectors have not yet observed a BNS post-merger (Abbott et al. 2017a(Abbott et al. ,b, 2020. Nevertheless, future detectors, like the proposed Einstein Telescope (ET) (Maggiore et al. 2020) and Cosmic Explorer (CE) (Abbott et al. 2017c) detectors, will feature improved sensitivity at the higher frequencies necessary to detect BNS post-merger signals. Inference on the EOS using post-merger data is possible at (postmerger) signal-to-noise ratios (SNR) as low as 8 (Breschi et al. 2022a). Additionally, sensitivity upgrades to current instruments promise higher BNS detection counts with better sky localization (Abbott et al. 2018b).
State-of-the-art BNS merger simulations typically incorporate thermal effects via full finite-temperature EOSs, often in the form of a table (Bauswein et al. 2010;Sekiguchi et al. 2011a), and realistic neutrino transport, such as via elaborate moment approximations (Foucart et al. 2016;Radice et al. 2022) or Monte Carlo methods (Foucart et al. 2022). Many studies perform simulations with multiple finite-temperature EOSs to demonstrate sensitivity (or lack thereof, as the case may be) of BNS merger observables under different scenarios, but the different cold-temperature behavior of each EOS makes it difficult to attribute specific outcomes to finitetemperature behavior (Neilsen et al. 2014;Radice et al. 2018b;Most et al. 2019;Perego et al. 2019;Hammond et al. 2021;Camilletti et al. 2022). Some studies have explored systematic changes in thermal effects through a so-called "hybrid EOS", which extends a cold nuclear EOS to finite temperatures using an ideal gas component with a fixed adiabatic constant Γ th (Bauswein et al. 2010;Figura et al. 2020). However, this is only a very rough approximation, as the effective Γ th of a full finitetemperature EOS varies considerably with density, temperature, and composition (Carbone & Schwenk 2019). Raithel et al. (2021) more recently considered finitetemperature effects through a more sophisticated hybrid EOS based on an approximation to the effective mass, but they do not fully explore how their model affects post-merger GW signals, nor is this approach an entirely self-consistent model. Furthermore, none of these studies incorporate all the relevant physics for modeling thermal effects, particularly consistent neutrino transport.
In this Letter, we present a first GR neutrino-radiation hydrodynamics study of finite-temperature effects of a realistic nuclear EOS on BNS mergers through modifications to the specific heat capacity. Our simulations show that an increased heat capacity results in denser, cooler remnants. This leaves clear imprints on the GW signal in the post-merger phase, which we show can be recovered in a parameter estimation pipeline tuned to a next-generation GW observatory.

METHODS.
We select three non-relativistic Skyrme-type nucleonic EOSs built with the framework of Schneider et al. (2019) and parameterized to produce the same cold nuclear matter bulk properties but different specific heat content. In Skyrme EOSs, the specific heat is controlled by the temperature-independent effective masses of neutrons and protons, m * n and m * p , respectively. These have a simple phenomenological description (Constantinou et al. 2014;Schneider et al. 2019;Margueron et al. 2018) that depends only on two parameters and on the nucleonic number densities, n n and n p (or, alternatively, the number density n = n n + n p and the proton fraction Y p = n p /n of matter), and converge toward the vacuum nucleon masses at zero density. The parameters were chosen to reproduce two nuclear matter observables at saturation density: the effective mass for symmetric nuclear matter, m * = m * n ≃ m * p , and the neutron-proton effective mass splitting for pure neutron matter, ∆m * . Guided by theoretical and experimental efforts (Margueron et al. 2018;Li et al. 2018;Huth et al. 2021;Zhang et al. 2021), the selected EOSs probe the average and extreme, but still plausible, expected values for m * , m * = {0.55, 0.75, 0.95} m n , while fixing the yet poorly constrained ∆m * to 0.10 m n . The same EOSs were used to study GW signals from core-collapse supernovae (Andersen et al. 2021).
To first order, the baryonic contribution to the specific heat of degenerate matter found in the core of a neutron star, which dominates over the lepton contribution, Constantinou et al. (2014). Thus, all else being equal, increasing m * leads to a larger specific heat capacity for matter in the merger remnant, to which we attribute the differences seen across our simulations.
Using the pseudospectral code LORENE (Gourgoulhon et al. 2016), we construct initial data for equal-mass binary neutron star systems in quasicircular orbit with a gravitational mass of M = 1.35 M ⊙ per star. We evolve each binary using THC M1 (Radice et al. 2022;Zappa et al. 2023), an extension of the THC numerical relativity code (Radice et al. 2014a,b) incorporating neutrino transport via a moment formalism (Thorne 1981;Shibata et al. 2011). The implementation in THC M1 makes use of the Minerbo closure for the radiation pressure tensor, which is exact in the optically thick limit. Thus while our neutrino treatment is approximate overall, we can capture effects such as neutrino trapping and dissipative effects from out-of-equilibrium weak reactions in the BNS remnant exactly. Our runs, while not modeling the magnetic field explicitly, also account for the effects of heating and angular momentum transport from magnetohydrodynamic (MHD) turbulence via a general relativistic large eddy simulation (GRLES) formalism calibrated with high-resolution GRMHD BNS simulations (Radice 2017;Kiuchi et al. 2018;Radice 2020).
We perform each simulation at two resolutions, designated as LR and SR, respectively corresponding to a grid spacing of ∆x ≈ 250 m and ∆x ≈ 180 m in the finest refinement level, which covers both stars during the inspiral and merger phases. We also run identical SR simulations with our older M0 solver (Radice et al. 2016) to validate our results; though the solver is less accurate and does not properly capture effects such as neutrino trapping, these runs support the major results of the M1 runs.  at most post-merger times. As suggested in this plot, all of the models, independent of resolution or neutrino treatment, produce long-lived remnants.
The GW strain in Figure 2 (bottom panel) demonstrates identical behavior in the inspiral for all three models, but the waveforms begin to deviate in the postmerger due to finite-temperature effects in the EOS. The differences in morphology include frequency evolution, amplitude, damping times, and modulations. This is more quantitatively seen in the GW spectrum (see Figure 3), where there is a clear rightward shift in the peak postmerger frequency, f 2 , as m * increases. Table 1 contains f 2 for both the LR and SR simulations. Although the precision of the NR waveform is limited by finite resolution and step size, the shift of ∆f ≳ 10 Hz is robust across resolutions, which suggests it is an effect of the EOS and not an artifact of the simulations. The M0 SR simulation demonstrates a similar trend, although the m * = 0.75 run is much closer to the m * = 0.95 model than either of the M1 cases. This provides an important sanity check of our results, but we stress that it is only qualitative; the M1 results should be used for quantitative analysis of post-merger effects, as a self-consistent approach which accounts for neutrino trapping noticeably changing the remnant's evolution (see Zappa et al. (2023); the appendix also contains a more detailed comparison). To detect these thermal effects in the postmerger via next-generation GW detectors and possibly constrain the m * nuclear parameter, we perform full Bayesian inference on the postmerger GW signals. To compute injections, we extract the postmerger waveforms from the SR simulations by applying a Tukey window to suppress the inspiral and spline interpolate the GW waveforms to a sampling rate of 16384 Hz. We further zero pad the signals to a segment of 1 s. For brevity, we consider here only noise-less injections. For parameter estimation, we employ the publicly available code bajes (Breschi et al. 2021) and use the UltraNest (Buchner 2021) sampler available as part of the bajes pipeline. We recover the injections by using the postmerger model NRPMw from (Breschi et al. 2022a) and compute posteriors on its parameters. We inject all signals at a luminosity distance corresponding to a post-merger SNR of 15 (which corresponds to an inspiral SNR of O(100), see Maggiore et al. (2020)), assuming the power spectral density of ET-D (Hild et al. 2011) to simulate the detector response. The priors are set in accordance to (Breschi et al. 2022b), section-IIB. In Figure 3, we show the reconstructed waveforms from NRPMw evaluated on the parameter space of the recovered posterior samples. We list the recovered SNRs and f 2 values in Table 1. At SNR = 15, the injected spectrum's f 2 frequencies sit well within the 90% credible intervals of the distribution of reconstructed waveform's f 2 frequencies. Furthermore, these intervals do not overlap, indicating that all three Table 1. Peak post-merger frequencies (f2) of the gravitational wave spectrum for the LR (M1), M0, and SR (M1) simulations and the NRPMw model. For reference, the recovered matched-filter SNR values and corresponding luminosity distance DL are also provided. For consistency, all peaks are measured after suppressing the inspiral. We also provide the mismatch, M (see Lindblom et al. (2008); Damour et al. (2011)), between the f SR 2 runs, with one row measured against m * = 0.75 and the other against m * = 0.95. We have shown that m * significantly influences the GW signals in BNS mergers, and we have analyzed this effect in detail for f 2 . The relationship between m * and f 2 is easily explained in terms of the specific heat. In-creasing the specific heat appears to soften the equation of state; because it requires more energy to increase the temperature, there is less thermal pressure available to support the star, thus producing a more rapidly rotating and compact remnant with lower temperatures, in agreement with core-collapse supernovae studies (Andersen et al. 2021;Schneider et al. 2019;Yasin et al. 2020). The EOSs we use have a simple relationship between m * and the specific heat capacity which may not be representative of the true nuclear EOS, which is expected to be both density and temperature dependent (Carbone & Schwenk 2019). Nevertheless, this study serves as a proof of concept demonstrating that future detectors like ET and CE can use f 2 to constrain the finitetemperature EOS. We also reiterate that this effect only affects the finite-temperature evolution, which will not be observable until next-generation detectors sensitive to the post-merger phase come online. Furthermore, we also note that this study does not provide a method for measuring m * from f 2 ; we only assert that m * leaves an imprint on f 2 .
We do acknowledge some limitations in our work; most notably, the length of the m * = 0.95 M1 simulation is much shorter than both the m * = 0.55 and m * = 0.75 runs due to a limitation in the M1 neutrino solver which introduced unphysical effects past 5 ms post-merger. This short length may explain why f m * =0.95 2 increases between the LR and SR runs while both f m * =0.55 2 and f m * =0.75 2 instead decrease, as a shorter signal will introduce a extra spread of frequencies to the power spectrum and possibly shift the peak. Another possible limitation concerns the omission of magnetic fields. Our simulations include a GRLES model which can account for some, but not all MHD effects. On the other hand, other simulations suggest that MHD effects on the post-merger gravitational wave frequency are negligible for realistic initial magnetic field strengths (Palenzuela et al. 2022).
We may consider several avenues for future work. Longer simulations would allow us to investigate the ejecta and consider possible effects on EM counterparts, as well as possible thermal effects on the lifetime of the remnants. Additionally, Zappa et al. (2023) indicate that resolution has a prominent influence on post-merger evolution, including collapse time and disk formation. Although we validate our results here with two resolutions, accurately determining the precise value of f 2 for each model would require higher resolution calculations. To investigate our hypothesis that this study's results are general, future simulations could also explore other EOS models with tunable finite-temperature behavior.
NOTE ADDED While this manuscript was under review, Raithel & Paschalidis (2023) announced new results on the impact of thermal effects on f 2 . Their work uses a semi-analytic prescription for the EOS and neglects neutrinos, but it is in good agreement with our findings.  The presence of a trapped neutrino gas in optically-thick regions can affect the temperature and composition of the remnant (Perego et al. 2019;Zappa et al. 2023). Many BNS simulations today employ leakage schemes, which do not explicitly model neutrino radiation, but instead estimate a cooling rate based on the local optical depth (Sekiguchi et al. 2011b;Neilsen et al. 2014;Palenzuela et al. 2015;Radice et al. 2016Radice et al. , 2018a. Because they do not perform consistent radiation transport, leakage schemes cannot capture the trapped neutrino gas in the remnant or potential out-of- The average temperature as a function of density at t ≈ 5 ms post-merger for the M0 (left) and M1 (right) SR simulations. Each average was calculated by constructing two-dimensional histograms in temperature and density, then averaging over the temperature with a weight corresponding to the mass of each density-temperature bin. To reduce noise in the data, five equally-spaced timesteps (with ∆t ≈ 0.05 ms) centered on t ≈ 5 ms were averaged together. equilibrium effects. The more accurate M1 scheme, however, directly models neutrino transport and can therefore capture both these effects (Radice et al. 2022;Zappa et al. 2023). Because our study represents the first investigation of thermal effects in BNS mergers with self-consistent neutrino transport, we will demonstrate here the influence of neutrinos on our results by comparing our M0 (a leakage-style scheme) simulations to our M1 runs.
In Fig. 4, we show the GW strain for both the M0 and M1 simulations. All three values of m * show identical behavior in the inspiral, which is to be expected; neutrinos emissions are highly sensitive to the temperature, which is quite low prior to merger. In the post-merger signal, the m * = 0.55 and m * = 0.75 differ only very slightly in amplitude and frequency. The m * = 0.95 runs show somewhat stronger deviations but are still quite similar in overall morphology. Nevertheless, the shifts in f 2 (see Table I in the main text) for changes in the neutrino solver are within a factor of two or three of the shifts due to changing m * ; were these results used to calibrate some model for determining m * from f 2 , it may result in a fairly significant error.
In Fig. 5, we show the average temperature as a function of density at 5 ms post-merger for both the M0 and M1 SR simulations. We calculate this temperature by constructing two-dimensional histograms in temperature and density, then performing a weighted average over the temperature for each density bin. This data is smoothed by averaging five time steps centered around 5 ms. The M1 data demonstrates a clear trend at higher densities (ρ ≳ ρ sat ) where higher values of m * lead to lower temperatures, and there is limited evidence suggesting an inverted trend in the outer layers of the star (ρ ≲ ρ sat ), a result which is consistent with core-collapse supernovae simulations (Schneider et al. 2019). One possible explanation is due to the increased compactness for higher m * ; material near the surface (where the density is low enough that m * has a much weaker effect on the EOS) falls deeper into the gravitational potential and heats up more as m * increases. However, the high-density trend is obscured in the M0 data, and there is no evidence of the low-density trend, as the m * = 0.55 temperature is often higher than the m * = 0.75 temperature at lower densities (where the M0 scheme suppresses neutrino absorption).
In weak equilibrium where µ i is the chemical potential for particle species i, with p, e, n, and ν e respectively representing the proton, electron, neutron, and electron neutrino. On the other hand, for a trapped neutrino gas in thermal equilibrium with the nucleonic matter of density ρ, temperature T , and electron fraction Y e , we have where m b is the average baryon mass and F i (x) is the i th Fermi function. This suggests that one can monitor the deviation from equilibrium by calculating the quantity where µ T νe is µ νe calculated under the assumption of thermal equilibrium (Eq. A2) using the evolved neutrino fractions from M1. For M0 this quantity should be set to zero, as the trapped neutrino gas is not explicitly modeled. We plot this deviation for ν e in Fig. 6. M0 shows noticeable deviations from weak equilibrium throughout the bulk of the remnant. This shows that the M0 simulations are not correctly capturing the thermal equilibrium of matter. On the other hand, the M1 simulations show that matter and neutrinos are in equilibrium in most regions in the remnant. We therefore conclude that a proper investigation of thermal effects must include full neutrino transport.