Nextgeneration gravitationalwave instruments are expected to constrain the equation of state of dense nuclear matter by observing binaries involving neutron stars. We highlight a problematic systematic error in finitetemperature merger simulations where shock heating associated with the neutronstar surface gives rise to elevated temperatures. We demonstrate the severe implications of this artificial heating by computing static and dynamical tidal parameters for neutron stars immersed in simulation temperature profiles. The induced systematic errors must be addressed if we want to build robust gravitationalwave signal models for neutronstar or indeed neutron starblack hole binaries.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena']
Problematic systematics in neutronstar merger simulations
{'Problematic systematics in neutron-star merger simulations': 'F. Gittins, R. Matur, N. Andersson, and I. Hawke Mathematical Sciences and STAG Research Centre, University of Southampton, Southampton SO17 1BJ, United Kingdom (Dated: September 23, 2024)  \n∗  \nNext-generation gravitational-wave instruments are expected to constrain the equation of state of dense nuclear matter by observing binaries involving neutron stars. We highlight a problematic systematic error in finite-temperature merger simulations, where shock heating associated with the neutron-star surface gives rise to elevated temperatures. We demonstrate the severe implications of this artificial heating by computing static and dynamical tidal parameters for neutron stars immersed in simulation temperature profiles. The induced systematic errors must be addressed if we want to build robust gravitational-wave signal models for neutron-star, or indeed neutron star-black hole, binaries.', 'I. INTRODUCTION': "Gravitational-wave astronomy has the potential to provide precise constraints on the extreme physics represented by neutron stars. The observation of the first binary neutron-star merger event, GW170817, in LIGOVirgo highlighted this promise through the absence of a distinguishable imprint of tidal interaction in the gravitational-wave signal [1-3]. This, in turn, placed constraints on the neutron-star tidal deformability and the permitted range of neutron-star radii. To date, GW170817 is the only observed event for which such constraints have been obtained. This is expected to change with the development of more sensitive gravitationalwave interferometers, like the Einstein Telescope [4] and Cosmic Explorer [5]. Neutron-star physics provides a key science driver for these advanced instruments.  \nDesigned to be about a factor of 50 more sensitive than the interferometers that were in operation in 2017, the next-generation detectors should allow us to put much tighter bounds on the equation of state of matter at supranuclear densities. This will involve precision measurements of the response of a neutron star to the tidal interaction with a companion during the late stages of compact binary inspiral. The tide induced on a neutron star has two aspects: At lower frequencies, the dominant feature is the static tide , represented by the tidal deformability [6-8]. As the binary orbit shrinks and the system approaches merger, the dynamical tide starts to play a role [9-12]. The dynamical tide has various features, but it is dominated by the contribution from the fundamental mode of oscillation of the star (the f -mode) [13]. Moreover, the tidal deformability and the f -mode frequency are known to be linked by a quasi-universal relation [14], which means that the two features may not be truly 'independent' from the parameter-extraction point of view. There are convincing arguments that we need to include both static and dynamical aspects in our tidal models if we want to maximise the information extracted from observations [15-17].  \nAdditional contributions to the dynamical tidee.g. , from low-frequency gravity g -modes arising from composition and entropy gradients [10, 18] or interface i -modes associated with sharp phase transitions [19, 20]-are expected to leave a weak imprint on the gravitational-wave signal. Detections of such fine-print contributions would shed light on the composition and state of the highdensity matter. This is an important issue, but not the focus of the discussion here.  \nGravitational-wave astronomy relies, in general, on a robust link between observations and the underlying theory. In particular, for the problem of neutron-star tides, waveform models need to be calibrated against two distinct bodies of work [21-27]. At low frequencies (early stages of inspiral), the problem is well described by perturbation theory coupled to a post-Newtonian model for orbital evolution. At higher frequencies-as the system approaches merger-the dynamics become nonlinear and require numerical simulations. Each approach to the problem comes with its own set of computational issues. The specific feature we highlight here relates to the late stages of evolution where numerical simulations are necessary.  \nA variety of approaches are currently used to build waveform models for binaries involving neutron stars. For example, the time-domain effective-one-body framework [21-23] calibrates the model using numericalrelativity simulations of merging black-hole binaries, while the frequency-domain closed-form tidal approximant method [24-27] normalises with respect to neutronstar simulations. There are also hybrid models that combine the two techniques, see for instance Ref. [28]. Arguably, the current state of the art relies on data from neutron-star merger simulations which implement fits of one-parameter nuclear-matter models [29], implicitly assuming the stellar material to be cold and the nuclear reactions to maintain equilibrium. However, realistic neutron-star mergers are hot, out-of-equilibrium events, so we ultimately need to work towards calibrations based on realistic, finite-temperature simulations [30-32] to capture the anticipated physics. Such simulations are essential if we want to connect the pre-merger behaviour with the violent dynamics, including electro-  \nmagnetic counterpart emission, of the post-merger remnant.  \nThe issue we bring to the fore relates to a (known, but somewhat ignored) systematic error present in all current (grid-based) finite-temperature neutron-star merger simulations involving shock-capturing schemes: The temperature in the simulated stars is much too high before merger. The problem is due to the simple fact that the employed numerical schemes treat the sharp drop in density at the neutron-star surface as a 'shock' leading to generation of (numerical) entropy and artificial heating. The level of heating varies between simulation codes, but it is generally the case that the temperature is ramped up to 5-10 MeV throughout the star very soon after the beginning of the simulation (see Fig. 1 and Appendix A). As a result, the simulated neutron star is several orders of magnitude too hot. In fact, at this level the thermal pressure is significant, effectively making the star much less compressed than it ought to be. Based on a comparison to proto-neutron stars, which typically reach similar temperatures, one would estimate that the neutron-star radius can change by as much as a factor of 2. This obviously impacts on both the tidal deformability and the f -mode frequency and hence may render debatable any waveform calibration involving such numerical simulation data. This is clearly problematic.  \nFIG. 1. The peak temperature reached in the neutron-star matter during a numerical simulation. Results correspond to the DD2 equation of state and are provided for both single isolated neutron stars (at three different resolutions; 369 m in black, 185 m in green and 92 m in purple) and one of the companions in a neutron-star binary simulation (with grid resolution 221 . 5 m in red). In the binary simulation, the initial separation is set to 45 km. The isolated neutron-star simulations suggest that the issue with the spurious heating remains for affordable levels of resolution. The comparison to the binary neutron-star case shows that the artificial heating is enhanced when the star is moving across the numerical grid. The temperatures indicated here are representative of those seen in all current finite-temperature (grid-based) neutron-star simulations.  \n<!-- image -->  \nAs a first step towards a better quantitative under-  \nstanding of the issue, we calculate both static and dynamical tidal parameters for two specific realistic finitetemperature equations of state. Our results provide a more precise illustration of the impact of the artificial heating and a clear indication of the systematics associated with current finite-temperature neutron-star simulations.", 'II. SYSTEMATICS IN FINITE-TEMPERATURE SIMULATIONS': 'The fluid equations that determine the neutron-star structure are closed by a thermodynamical equation of state . For our calculations, we assume a three-parameter model for the nuclear matter-for simplicity, only comprising neutrons, protons and electrons-with the temperature T , baryon-number density n b and electron fraction Y e as the natural thermodynamical variables. Such models represent the current state of the art in nuclear astrophysics [33] and are regularly implemented in numerical simulations [30-32].  \nFor the assumed three-parameter model, the first law of thermodynamics may be expressed as  \ndf = -n b s dT + f + p n b dn b + n b µ ∆ dY e , (1)  \nwhere s is the entropy per baryon, p is the isotropic pressure, µ ∆ = µ p + µ e -µ n and { µ n , µ p , µ e } are the neutron, proton and electron chemical potentials, respectively. The various quantities are as measured in the local inertial reference frame of the fluid. The free-energy density f = f ( T, n b , Y e ) represents the fundamental equation of state from which all the thermodynamical information about the stellar material can be derived.  \nIn order to illustrate the artificial surface heating in binary simulations, we use temperature results from two separate inspiral-merger simulations with different prescriptions for the finite-temperature matter: APR [34] and DD2 [35] (as implemented in the CompOSE library [33, 36, 37]). The first simulation is taken from Refs. [31, 38], which uses the Einstein Toolkit codebase to collide neutron stars described by the APR nuclear-matter equation of state. The simulation clearly exhibits the characteristic artificial surface heating, reaching temperatures of order 10 MeV, see Fig. 1 of Ref. [31]. In addition, mainly in order to demonstrate the behaviour with a different numerical set-up, we have performed simulations with the WhiskyTHC framework adopting the DD2 model for the microphysics. The details of the simulations are provided in Appendix A. We consider two physical situations: an isolated neutron star and a binary merger for which the two stars (obviously) move across the numerical grid.  \nThe results show that, at typical simulation resolutions, the bulk of an isolated neutron star heats up to temperatures of ∼ 1 -4 MeV (with the peak temperature slightly higher, see Figure 1) due to fluid shocks on  \nthe grid. The evolution of the peak temperature reached in each simulation is displayed in Fig. 1. For the DD2 merger simulation, the results show that the artificial heating is further enhanced when the stars move across the numerical grid.', 'III. IMPACT ON TIDES': "The main question we want to answer is: What is the impact of the observed artificial heating on the main tidal parameters that one would aim to extract from gravitational-wave observations? In order to explore this issue, we solve for the linear perturbations of a spherically symmetric, perfect-fluid relativistic star. Our perturbation calculation closely follows Ref. [39] and we provide further details in Appendix B. From the outset, we assume that the star is immersed in a temperature profile T = T ( n b ). (This will be either uniform or lifted from a numerical simulation.) To determine the corresponding matter composition, we assume that the background fluid is in chemical equilibrium such that  \nµ ∆ ( T, n b , Y e ) = 0 . (2)  \nThis fixes Y e = Y e ( n b ) and reduces the thermodynamical state to depend solely on n b . The structure of the equilibrium neutron star is then straightforwardly obtained by solving the standard equations supplemented with the functions ε = ε ( n b ) and p = p ( n b ) obtained from the equation of state. This provides the background on which the linear perturbations are computed. Although our setup is similar to Ref. [39], here we are concerned with different temperature profiles, in particular those extracted from merger simulations.  \nAs already mentioned, there are two regimes in a binary inspiral: the early regime, which is characterised by the static tide, and the late regime, where the dynamics become important (and non-linear aspects come into play). The static tide is commonly represented by the neutron star's tidal deformability Λ. This quantity enters the inspiral waveform through its mass-weighted average for the binary at fifth order in the post-Newtonian approximation [1, 40]. The tidal deformability provides a dimensionless measure of the star's susceptibility to a companion's gravitational field. If the nuclear matter is particularly stiff, the neutron star can support large deformations (large values of Λ). Meanwhile, if the stellar material is comparatively soft, the neutron star will be more compact and have a small Λ. These qualities make Λ an attractive observable with which to constrain the properties of dense nuclear matter [40], as demonstrated in the case of GW170817 [1-3].  \nWe determine the static, quadrupolar deformations using the formalism detailed in Ref. [6] to obtain Λ. The required static perturbations depend on an additional aspect of the nuclear matter, the equilibrium adiabatic in-  \nex  \nΓ = ε + p p dp dε , (3)  \nwhich is determined from partial derivatives of the thermodynamical functions (see Ref. [39]).  \nDuring the later, dynamical regime, the tidal driving frequency increases as the binary inspirals, radiating gravitational waves. As the orbital frequency increases, it will eventually become resonant with the low-frequency g -modes of the neutron star, as well as implicate the other higher frequency oscillations, like the f -and p -modes [9, 10]. Indeed, the f -mode is expected to provide the dominant contribution to the tide [13]. At this point, the assumption of a static tide breaks down and the problem becomes dynamical. In the absence of dissipation, the Newtonian oscillation problem is self-adjoint, implying that the modes form a complete basis in terms of which the tide can be decomposed. This mode-sum representation elegantly addresses the mathematical challenges of solving the full time-dependent, tidal response problem. Although motion in general relativity is inherently dissipative due to gravitational radiation, the natural vibrational modes of neutron stars are still expected to dominate the dynamical tide.  \nTo represent the dynamical tide, we calculate the quadrupolar, quasi-normal oscillation modes of a neutron star using the equations provided in Refs. [41, 42]. The equation of state now enters the fluid perturbations through the adiabatic index  \nΓ 1 = ε + p p ( ∂p ∂ε ) s,Y e , (4)  \nassuming frozen composition during an oscillation ( i.e. , that the mode dynamics are fast enough that we may ignore nuclear reactions). Similarly to Γ, the determination of Γ 1 involves thermodynamical derivatives (for more detail, see Ref. [39]).  \nThe two adiabatic indices Γ and Γ 1 characterise stratification in the fluid. In general, stellar material can support both entropy and composition gradients. When Γ 1 ≥ Γ, the star is convectively stable and supports lowfrequency g -mode oscillations. Meanwhile, when Γ 1 = Γ, the g -modes vanish. For realistic nuclear-matter models, the neutron star is expected to be stably stratified throughout, except possibly at low densities, and will therefore support g -mode oscillations. Indeed, a recent study has shown that these low-frequency fluid oscillations, which likely become resonant with the tide as the binary inspirals, may be within reach of next-generation observatories [43].  \nWe present a sample of numerical results in Table I. For each of the two nuclear-matter models we consider, we show the tidal parameters (static and dynamical) when the star has cold, uniform temperature and when the star is placed in a temperature profile obtained from a numerical-relativity merger simulation (modelled on the  \nresults in Fig. 4 in Appendix A for DD2), assuming the same matter model (and stars with the same baryon mass). This provides an immediate quantitative measure of the impact of the artificial simulation temperatures on the tidal parameters.  \nTABLE I. The effects on the tidal parameters due to artificial temperatures in numerical-relativity merger simulations. The table lists the stellar radius R , tidal deformability Λ, f -mode frequency ω f and first three g -mode frequencies ω g 1 , ω g 2 and ω g 3 . For each nuclear-matter model, APR and DD2, the parameters are shown for two neutron stars with different temperatures; one neutron star has uniform temperature, whereas the other is immersed in a merger simulation profile. The stars described by APR have baryon mass M b = 1 . 40 M ⊙ and those for DD2 have M b = 1 . 55 M ⊙ . It is evident that the neutron stars with simulation-level temperatures have significantly altered tidal parameters compared to the cold stars.  \nFor the APR matter model, the numerical results show that Λ is 16% larger due to the artificially high temperature in the numerical simulation than for the corresponding cold neutron star. The difference is even starker with DD2, for which the deformability increases by more than 25% with the simulation temperature profile. It should be noted that the two equation-of-state cases are not intended to be directly compared to each other. Indeed, the involved neutron stars have different masses, as well as different numerical treatments for the hydrodynamics and the spacetime evolution. Rather, the point of the analysis is to illustrate that the artificial temperatures in two separate merger simulations lead to the same qualitative behaviour: The tidal parameters are significantly distorted.  \nThe results for the tidal deformability are quite intuitive. The artificially high simulation temperatures contribute to a strong thermal pressure, which leads to larger stellar radii. In fact, from Table I, we see that the neutron-star radii approximately double. This leads to an effective stiffening of the nuclear matter and means that Λ-which scales with the areal radius R as ∝ R 5 -increases. It is not particularly surprising that higher temperatures lead to larger tidal deformabilities (this feature is generic). What is notable is the extent to which the values change. If we were to calibrate gravitationalwaveform models with these kinds of simulations, these aspects would inevitably manifest as systematic errors that bias parameter inference.  \nWe now turn our attention to the effect on the dynamics. Presented in Table I are the (real) oscillation  \nfrequencies of the f -mode and first three g -modes. The results show that these are also drastically altered by the artificial temperature. We see that, when the neutron star is placed in the temperature profile of a simulation, the spectrum shifts considerably. The f -mode is the most weakly affected of the oscillations, decreasing by 16% for the APR model, but only 5% for DD2. In contrast, the low-frequency g -modes are substantially altered, increasing in frequency due to enhanced entropy gradients caused by the high temperatures.  \nFinally, we calculate the tidal parameters' evolution with respect to uniform temperature. The results are presented in Fig. 2 for APR nuclear matter and are consistent with the behaviour we have just discussed: The tidal deformability and g -mode frequencies increase with temperature, while the f -mode gradually oscillates slower. We see that the f -mode is the least affected by temperature. In each panel of Fig. 2, we show the tidal values from Table I using the simulation profile for comparison. The tidal deformability and g -mode frequencies in the simulation are very roughly similar to that of a star with uniform T ∼ 10 MeV. However, even at T = 15MeV, the frequency of the f -mode remains higher than that of the simulated neutron star.", 'IV. WORDS OF CAUTION': "We have demonstrated that the key tidal parameters of neutron stars in finite-temperature, merger simulations are severely distorted due to artificial heating arising from numerical shocks close to the stellar surface. Using data from two numerical-relativity simulations, with different numerical methodologies and matter prescriptions, we showed how quantities associated with both the static and dynamical tide are substantially shifted by the artificially high temperatures. Specifically, this includes the neutron-star tidal deformability, already constrained by the GW170817 observation, as well as the oscillation modes that dominate the dynamical tide, which we hope to probe with upcoming, next-generation gravitationalwave interferometers, the Einstein Telescope and Cosmic Explorer. Although we considered only two simulations, we expect these features to be generic for all grid-based numerical-relativity codes, as they all reach unphysically high temperatures.  \nThe main take-away message from this work is cautionary. If we were to use results from finite-temperature non-linear simulations to calibrate gravitational waveforms, we may incur considerable systematic error in the parameter inference. Current instruments are not sensitive enough to the tide for these issues to manifest, but the problem will need to be addressed if we want to make precision observations with the Einstein Telescope and Cosmic Explorer.  \nFuture work will need to be dedicated to either reducing the temperature systematics or correcting for them in the gravitational-wave analyses. As a first step, we  \n<!-- image -->  \nFIG. 2. The tidal deformability (left panel) and mode eigenfrequencies (right panel) against uniform temperature of an M b = 1 . 40 M ⊙ neutron star described by the APR equation of state. The frequencies correspond to the f - and first three g -modes, as indicated by the legend. The horizontal, dashed lines correspond to the parameters listed in Table I calculated using the simulation temperature profile. As the neutron star is heated up, the f -mode frequency decreases with larger neutron-star radius. The other quantities, including the tidal deformability and g -mode frequencies, increase.  \n<!-- image -->  \nneed to better understand the origin of the problem. This includes reconciling the results we have presented here with previous work using simpler equation-of-state models (like polytropes) and a phenomenological thermal (Γ-law) representation. For example, the results in Ref. [44] demonstrate much better agreement between simulations and perturbation theory. To some extent, the difference in the results likely stems from the fact that the phenomenological equation-of-state models assume chemical equilibrium. To illustrate this, consider Eq. (1) as a relation for the numerical errors. By enforcing chemical equilibrium we have µ ∆ = 0 and the final term does not contribute. Away from chemical equilibrium we would need the error dY e to 'match' the other errors to give a sensible temperature. However, in current simulations the evolved quantity is (roughly) n b Y e , and hence the accuracy will be poor at low densities. This is likely related to the observation-see for example the results presented in Fig. 3 of [31]-that finite-temperature equation-of-state simulations are prone to large errors at low temperatures. Of course, these arguments only hint at the origin of the problem we are discussing. We do not yet have a solution.  \nIn order to make progress, it would be helpful to understand to what extent the issue we have raised is relevant for particle-based simulations, like those in Refs. [45, 46]. In principle, such simulations provide a better representation of the neutron-star surface, but it remains to be established to what extent this alleviates the artificial heating issue.  \nLet us make two final remarks. First, and this is important in order to keep the discussion in the proper context, the problem we have discussed does not influence current gravitational-waveform models for neutron-star  \ntides. These models, like the work in Ref. [27], do not (yet) involve thermodynamically consistent matter models and hence do not suffer the artifical heating problem (at least not at the level we have indicated here). Having said that, if we want to do better in the future (and we do!), then we need to get a handle on the temperature problem. Second, the artificial heating of low-density matter may have considerable impact on the related problem of matter ejecta, for which finite-temperature simulations are being used [47]. In this problem, the additional thermal pressure may tend to unbind matter and an elevated temperature will affect the composition (the electron fraction) of the outflows, which may in turn alter the nuclear reaction rates and the associated kilonova signature. Similarly, an artifically high temperature will impact on neutrino opacities, and hence need to be carefully considered in efforts to implement realistic neutrino transport [48]. If this concern is justified remains to be explored.", 'ACKNOWLEDGMENTS': 'NA and IH acknowledge support from STFC via grant number ST/R00045X/1. The authors thank Peter Hammond for sharing temperature data from a numericalrelativity merger simulation and Tim Dietrich and Nick Stergioulas for useful discussions. We acknowledge the use of the IRIDIS High Performance Computing Facility and associated support services at the University of Southampton. The software developed to support this article is available in a GitHub repository [49] and is written in the Julia programming language [50-54]. The figures were generated using matplotlib [55, 56].', 'Appendix A: Numerical simulations': 'For the finite-temperature inspiral-merger simulations we considered two distinct models. The first simulation used the APR matter model [31, 38] and was performed using the Einstein Toolkit [57]. The initial data was created using Lorene [58], while the hydrodynamical and spacetime evolution was performed using GRHydro [59-62] and McLachlan [63, 64], respectively. McLachlan uses the BSSN formulation [65-67] of the Einstein equations. The simulation was performed by setting the atmosphere temperature and rest-mass density to 0 . 02 MeV and 6 . 2 × 10 6 g cm -3 , respectively. A fourth-order Runge-Kutta method was employed, with the Courant-Friedrichs-Lewy condition set to 0 . 25. Neutron stars were tracked using NSTracker .  \nThe second set of simulations, for the DD2 matter model, involved initial data generated using the Fuka solver [68] and the evolution of the hydrodynamics and spacetime are performed with WhiskyTHC [69-72] and CTGamma [73], respectively. We used the z4c formulation [74] of the Einstein equations, the Carpet adaptive mesh refinement driver (as in the APR model) [75] of Cactus [76], and track the neutron stars with the BNSTrackerGen component of WhiskyTHC . The atmosphere temperature and rest-mass density are 0 . 02 MeV and 6 . 2 × 10 3 g cm -3 , respectively. The time evolution was carried out using a fourth-order Runge-Kutta method and the Courant-Friedrichs-Lewy condition is set to 0 . 15. We used the Sophie Kowalevski release of the Einstein Toolkit [57, 77].  \nFirst of all, the artificial heating for simulations of isolated (static) neutron stars and three different numerical resolutions, corresponding to scales of ∼ 369 m, ∼ 185 m and ∼ 92 m, respectively, is provided in Fig. 3. The  \n- [1] B. P. Abbott et al. , GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119 , 161101 (2017).\n- [2] B. P. Abbott et al. , GW170817: Measurements of Neutron Star Radii and Equation of State, Phys. Rev. Lett. 121 , 161101 (2018).\n- [3] B. P. Abbott et al. , Properties of the Binary Neutron Star Merger GW170817, Phys. Rev. X 9 , 011001 (2019).\n- [4] M. Punturo et al. , The Einstein Telescope: a thirdgeneration gravitational wave observatory, Classical Quant. Grav. 27 , 194002 (2010).\n- [5] D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Barsotti, G. Billingsley, D. A. Brown, Y. Chen, D. Coyne, R. Eisenstein, M. Evans, P. Fritschel, E. D. Hall, A. Lazzarini, G. Lovelace, J. Read, B. S. Sathyaprakash, D. Shoemaker, J. Smith, C. Torrie, S. Vitale, R. Weiss, C. Wipf, and M. Zucker, Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, in Bull. Am. Astron. Soc. , Vol. 51 (2019) p. 35,  \nresults show that the star heats up to temperatures of ∼ 1 -4 MeV due to fluid shocks on the grid. It is notable that, while the shock originate at the surface, the heat rapidly propagates to the stellar interior as the simulation progresses.  \nThe inspiral-merger simulation shows that the artificial heating is enhanced when the stars move across the numerical grid. The typical temperature profile extracted from our DD2 simulation, demonstrating that the peak temperatures reaches well above 10 MeV, is shown in Fig. 4.', 'Appendix B: Perturbative mode calculations': "The mode calculations closely follow the analysis in Ref. [39]. Specifically, to address noise with lowfrequency oscillations, we use an augment developed in Ref. [78] and perturbations in the exterior are calculated using the method from Ref. [79]. A mode with discrete (complex) frequency ω corresponds to a solution such that the ingoing gravitational waves in the exterior, with amplitude ˜ A in , vanish.  \nThe mode-calculations show that, at low temperatures in the APR neutron star, the first g -mode-the g 1 -mode-is an interfacial i -mode that arises due to the core-crust phase transition [39]. Such oscillations appear when there is a discontinuity in the energy density of the equation of state. This identification is made clear by examining the mode's eigenfunctions, which possess a sharp cusp at the point in the star where the discontinuity lies. It is notable that this mode oscillates at a relatively high frequency of Re( ω/ 2 π ) = 508 . 5 Hz. To illustrate the changes, we show in Fig. 5 the neutron-star oscillation spectra for the finite-temperature DD2 equation of state. The qualitative results for APR, summarised in Table I, are the same.  \narXiv:1907.04833 [astro-ph.IM].  \n- [6] T. Hinderer, Tidal Love Numbers of Neutron Stars, Astrophys. J. 677 , 1216 (2008).\n- [7] T. Binnington and E. Poisson, Relativistic theory of tidal Love numbers, Phys. Rev. D 80 , 084018 (2009).\n- [8] T. Damour and A. Nagar, Relativistic tidal properties of neutron stars, Phys. Rev. D 80 , 084035 (2009).\n- [9] A. Reisenegger and P. Goldreich, Excitation of Neutron Star Normal Modes during Binary Inspiral, Astrophys. J. 426 , 688 (1994).\n- [10] D. Lai, Resonant Oscillations and Tidal Heating in Coalescing Binary Neutron Stars, Mon. Not. R. Astron. Soc. 270 , 611 (1994).\n- [11] T. Pitre and E. Poisson, General relativistic dynamical tides in binary inspirals without modes, Phys. Rev. D 109 , 064004 (2024).\n- [12] A. Hegade K. R., J. L. Ripley, and N. Yunes, Dynamical tidal response of nonrotating relativistic stars, Phys. Rev. D 109 , 104064 (2024).  \nFIG. 3. An illustration of the artificial heating associated with shocks at the neutron-star surface. The results are obtained from the evolution of a single isolated neutron star and the three panels represent different simulation resolutions. The left, middle and right panels show the low-resolution (LR), medium-resolution (MR) and high-resolution (HR) simulation results, with a grid spacing of ∼ 369 m, ∼ 185 m and ∼ 92 m in the finest grid, respectively. The neutron stars are modelled using the finite-temperature, composition-dependent DD2 equation of state. While the temperature in the core of the NS reaches ∼ 4 MeV in the low-resolution simulation, it reduces to ∼ 1 -2 MeV for the high-resolution simulation within nearly ∼ 1 . 6 ms.  \n<!-- image -->  \nFIG. 4. The typical temperature profile (in the equatorial plane) extracted from the DD2 binary neutron-star simulation. The histogram illustrates the distribution of baryon mass in the density-temperature plane. The data is extracted after 3 . 9 ms of evolution.  \n<!-- image -->  \n- [13] N. Andersson and P. Pnigouras, Exploring the effective tidal deformability of neutron stars, Phys. Rev. D 101 , 083001 (2020).\n- [14] T. K. Chan, Y. H. Sham, P. T. Leung, and L. M. Lin, Multipolar universal relations between f -mode frequency and tidal deformability of compact stars, Phys. Rev. D", '90 , 124023 (2014).': '- [15] G. Pratten, P. Schmidt, and T. Hinderer, Gravitationalwave asteroseismology with fundamental modes from compact binary inspirals, Nat. Commun. 11 , 2553 (2020).\n- [16] N. Williams, G. Pratten, and P. Schmidt, Prospects for distinguishing dynamical tides in inspiralling binary neutron stars with third generation gravitational-wave detectors, Phys. Rev. D 105 , 123032 (2022).\n- [17] G. Pratten, P. Schmidt, and N. Williams, Impact of Dynamical Tides on the Reconstruction of the Neutron Star Equation of State, Phys. Rev. Lett. 129 , 081102 (2022).\n- [18] N. Andersson and P. Pnigouras, The g-mode spectrum of reactive neutron star cores, Mon. Not. R. Astron. Soc. 489 , 4043 (2019).\n- [19] D. Tsang, J. S. Read, T. Hinderer, A. L. Piro, and R. Bondarescu, Resonant Shattering of Neutron Star Crusts, Phys. Rev. Lett. 108 , 011102 (2012).\n- [20] D. Neill, W. G. Newton, and D. Tsang, Resonant shattering flares as multimessenger probes of the nuclear symmetry energy, Mon. Not. R. Astron. Soc. 504 , 1129 (2021).\n- [21] S. Bernuzzi, A. Nagar, T. Dietrich, and T. Damour, Modeling the Dynamics of Tidally Interacting Binary Neutron Stars up to the Merger, Phys. Rev. Lett. 114 , 161103 (2015).\n- [22] T. Hinderer, A. Taracchini, F. Foucart, A. Buonanno, J. Steinhoff, M. Duez, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, B. Szilagyi, K. Hotokezaka, K. Kyutoku, M. Shibata, and C. W. Carpenter, Effects of NeutronStar Dynamic Tides on Gravitational Waveforms within the Effective-One-Body Approach, Phys. Rev. Lett. 116 ,  \n<!-- image -->  \nFIG. 5. The oscillation spectra of an M b = 1 . 55 M ⊙ neutron star described by the finite-temperature DD2 equation of state for two different temperature profiles; uniform temperature T = 0 . 2 MeV (left panel) and the simulation profile (right panel). Here, ˜ A in represents the amplitude of ingoing gravitational radiation and ω is the (real) frequency of the perturbations. A mode of the neutron star corresponds to when ˜ A in = 0. The weakly damped oscillations can be identified by the singularities in the spectrum. (Rapidly damped, spacetime w -modes are not visible.) The vertical, dotted lines indicate the mode eigenfrequencies listed in Table I. The high temperatures in the merger simulation lead to a substantial distortion of the mode spectrum. The f -and p -mode frequencies decrease. We see in the right panel how the first p -mode becomes visible. The impact of the simulation temperatures on the g -modes is to increase their frequencies.  \n<!-- image -->', '181101 (2016).': '- [23] J. Steinhoff, T. Hinderer, A. Buonanno, and A. Taracchini, Dynamical tides in general relativity: Effective action and effective-one-body Hamiltonian, Phys. Rev. D 94 , 104028 (2016).\n- [24] T. Dietrich, S. Bernuzzi, and W. Tichy, Closed-form tidal approximants for binary neutron star gravitational waveforms constructed from high-resolution numerical relativity simulations, Phys. Rev. D 96 , 121501 (2017).\n- [25] T. Dietrich, S. Khan, R. Dudi, S. J. Kapadia, P. Kumar, A. Nagar, F. Ohme, F. Pannarale, A. Samajdar, S. Bernuzzi, G. Carullo, W. Del Pozzo, M. Haney, C. Markakis, M. Purrer, G. Riemenschneider, Y. E. Setyawati, K. W. Tsang, and C. Van Den Broeck, Matter imprints in waveform models for neutron star binaries: Tidal and self-spin effects, Phys. Rev. D 99 , 024029 (2019).\n- [26] T. Dietrich, A. Samajdar, S. Khan, N. K. JohnsonMcDaniel, R. Dudi, and W. Tichy, Improving the NRTidal model for binary neutron star systems, Phys. Rev. D 100 , 044003 (2019).\n- [27] A. Abac, T. Dietrich, A. Buonanno, J. Steinhoff, and M. Ujevic, New and robust gravitational-waveform model for high-mass-ratio binary neutron star systems with dynamical tidal effects, Phys. Rev. D 109 , 024062 (2024).\n- [28] K. Kawaguchi, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Frequency-domain gravitational waveform models for inspiraling binary neutron stars, Phys. Rev. D 97 , 044044 (2018).\n- [29] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Constraints on a phenomenologically parametrized neutron-star equation of state, Phys. Rev. D 79 , 124032 (2009).\n- [30] A. Perego, S. Bernuzzi, and D. Radice, Thermodynamics conditions of matter in neutron star mergers, Eur. Phys.', 'J. A 55 , 124 (2019).': "- [31] P. Hammond, I. Hawke, and N. Andersson, Thermal aspects of neutron star mergers, Phys. Rev. D 104 , 103006 (2021).\n- [32] J. Fields, A. Prakash, M. Breschi, D. Radice, S. Bernuzzi, and A. d. S. Schneider, Thermal Effects in Binary Neutron Star Mergers, Astrophys. J. 952 , L36 (2023).\n- [33] M. Oertel, M. Hempel, T. Klahn, and S. Typel, Equations of state for supernovae and compact stars, Rev. Mod. Phys. 89 , 015007 (2017).\n- [34] A. S. Schneider, C. Constantinou, B. Muccioli, and M. Prakash, Akmal-Pandharipande-Ravenhall equation of state for simulations of supernovae, neutron stars, and binary mergers, Phys. Rev. C 100 , 025803 (2019).\n- [35] M. Hempel and J. Schaffner-Bielich, A statistical model for a complete supernova equation of state, Nucl. Phys. A837 , 210 (2010).\n- [36] S. Typel, M. Oertel, and T. Klahn, CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose.obspm.fr, Physics of Particles and Nuclei 46 , 633 (2015).\n- [37] S. Typel, M. Oertel, T. Klahn, D. Chatterjee, V. Dexheimer, C. Ishizuka, M. Mancini, J. Novak, H. Pais, C. Providencia, A. Raduta, M. Servillat, and L. Tolos, CompOSE Reference Manual, arXiv:2203.03209 [astroph.HE] (2022).\n- [38] P. Hammond, Source, parfile, and plotting scripts relating to 'Thermal aspects of neutron star mergers' (arXiv 2108.08649) (2021).\n- [39] F. Gittins and N. Andersson, Neutron-star seismology with realistic, finite-temperature nuclear matter, arXiv:2406.05177 [gr-qc] (2024).\n- [40] ' E. Flanagan and T. Hinderer, Constraining neutron-star  \n- tidal Love numbers with gravitational-wave detectors, Phys. Rev. D 77 , 021502 (2008).\n- [41] L. Lindblom and S. L. Detweiler, The quadrupole oscillations of neutron stars., Astrophys. J. Suppl. Ser. 53 , 73 (1983).\n- [42] S. Detweiler and L. Lindblom, On the nonradial pulsations of general relativistic stellar models, Astrophys. J. 292 , 12 (1985).\n- [43] W. C. G. Ho and N. Andersson, New dynamical tide constraints from current and future gravitational wave detections of inspiralling neutron stars, Phys. Rev. D 108 , 043003 (2023).\n- [44] J. A. Font, N. Stergioulas, and K. D. Kokkotas, Nonlinear hydrodynamical evolution of rotating relativistic stars: numerical methods and code tests, Mon. Not. R. Astron. Soc. 313 , 678 (2000).\n- [45] A. Bauswein, N. Stergioulas, and H. T. Janka, Revealing the high-density equation of state through binary neutron star mergers, Phys. Rev. D 90 , 023002 (2014).\n- [46] S. Rosswog and P. Diener, SPHINCS BSSN: Numerical Relativity with Particles (2024), arXiv:2404.15952 [grqc].\n- [47] V. Nedora, F. Schianchi, S. Bernuzzi, D. Radice, B. Daszuta, A. Endrizzi, A. Perego, A. Prakash, and F. Zappa, Mapping dynamical ejecta and disk masses from numerical relativity simulations of neutron star mergers, Classical Quant. Grav. 39 , 015008 (2022).\n- [48] P. L. Espino, P. Hammond, D. Radice, S. Bernuzzi, R. Gamba, F. Zappa, L. F. L. Micchi, and A. Perego, Neutrino Trapping and Out-of-Equilibrium Effects in Binary Neutron-Star Merger Remnants, Phys. Rev. Lett. 132 , 211001 (2024).\n- [49] F. Gittins, RealisticSeismology (2024).\n- [50] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Rev. 59 , 65 (2017).\n- [51] C. Rackauckas and Q. Nie, DifferentialEquations.jl-A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia, J. Open Res. Software 5 , 15 (2017).\n- [52] M. Kittisopikul, T. E. Holy, and T. Aschan, Interpolations.jl (2022).\n- [53] A. Pal, F. Holtorf, A. Larsson, T. Loman, Utkarsh, F. Schaefer, Q. Qu, A. Edelman, and C. Rackauckas, NonlinearSolve.jl: High-Performance and Robust Solvers for Systems of Nonlinear Equations in Julia, arXiv:2403.16341 [math.NA] (2024).\n- [54] P. K. Mogensen and A. N. Riseth, Optim: A mathematical optimization package for Julia, J. Open Source Software 3 , 615 (2018).\n- [55] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng. 9 , 90 (2007).\n- [56] C. Rowley, PythonCall.jl (2022).\n- [57] F. Loffler et al. , The Einstein Toolkit: a community computational infrastructure for relativistic astrophysics, Classical Quant. Grav. 29 , 115001 (2012).\n- [58] E. Gourgoulhon, P. Grandcl'ement, and J. Novak, LORENE.\n- [59] I. Hawke, F. Loffler, and A. Nerozzi, Excision methods for high resolution shock capturing schemes applied to general relativistic hydrodynamics, Phys. Rev. D 71 , 104006 (2005).\n- [60] L. Baiotti, I. Hawke, P. J. Montero, F. Loffler, L. Rezzolla, N. Stergioulas, J. A. Font, and E. Seidel, Three-d  \n- mensional relativistic simulations of rotating neutronstar collapse to a Kerr black hole, Phys. Rev. D 71 , 024035 (2005).\n- [61] B. Giacomazzo and L. Rezzolla, WhiskyMHD: a new numerical code for general relativistic magnetohydrodynamics, Classical Quant. Grav. 24 , S235 (2007).\n- [62] P. Mosta, B. C. Mundim, J. A. Faber, R. Haas, S. C. Noble, T. Bode, F. Loffler, C. D. Ott, C. Reisswig, and E. Schnetter, GRHydro: a new open-source generalrelativistic magnetohydrodynamics code for the Einstein toolkit, Classical Quant. Grav. 31 , 015005 (2014).\n- [63] P. Diener, E. Schnetter, and J. Tao, McLachlan, a Public BSSN Code.\n- [64] D. Brown, P. Diener, O. Sarbach, E. Schnetter, and M. Tiglio, Turduckening black holes: An analytical and computational study, Phys. Rev. D 79 , 044023 (2009).\n- [65] T. Nakamura, K. Oohara, and Y. Kojima, General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes, Prog. Theor. Phys. Suppl. 90 , 1 (1987).\n- [66] M. Shibata and T. Nakamura, Evolution of threedimensional gravitational waves: Harmonic slicing case, Phys. Rev. D 52 , 5428 (1995).\n- [67] T. W. Baumgarte and S. L. Shapiro, Numerical integration of Einstein's field equations, Phys. Rev. D 59 , 024007 (1998).\n- [68] L. J. Papenfort, S. D. Tootle, P. Grandcl'ement, R. Most, E, and L. Rezzolla, New public code for initial data of unequal-mass, spinning compact-object binaries, Phys. Rev. D 104 , 024057 (2021).\n- [69] D. Radice and L. Rezzolla, THC: a new high-order finite-difference high-resolution shock-capturing code for special-relativistic hydrodynamics, Astron. Astrophys. 547 , A26 (2012).\n- [70] D. Radice, L. Rezzolla, and F. Galeazzi, Beyond secondorder convergence in simulations of binary neutron stars in full general relativity, Mon. Not. R. Astron. Soc. 437 , L46 (2014).\n- [71] D. Radice, L. Rezzolla, and F. Galeazzi, High-order fully general-relativistic hydrodynamics: new approaches and tests, Classical Quant. Grav. 31 , 075012 (2014).\n- [72] D. Radice, L. Rezzolla, and F. Galeazzi, High-Order Numerical-Relativity Simulations of Binary Neutron Stars, in Numerical Modeling of Space Plasma Flows ASTRONUM-2014 , Astronomical Society of the Pacific Conference Series, Vol. 498, edited by N. V. Pogorelov, E. Audit, and G. P. Zank (2015) p. 121, arXiv:1502.00551 [gr-qc].\n- [73] D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, and P. Diener, High accuracy binary black hole simulations with an extended wave zone, Phys. Rev. D 83 , 044045 (2011).\n- [74] S. Bernuzzi and D. Hilditch, Constraint violation in free evolution schemes: Comparing the BSSNOK formulation with a conformal decomposition of the Z4 formulation, Phys. Rev. D 81 , 084003 (2010).\n- [75] E. Schnetter, S. H. Hawley, and I. Hawke, Evolutions in 3D numerical relativity using fixed mesh refinement, Classical Quant. Grav. 21 , 1465 (2004).\n- [76] G. Allen, T. Goodale, J. Masso, and E. Seidel, The cactus computational toolkit and using distributed computing to collide neutron stars, in Proceedings. The Eighth International Symposium on High Performance Distributed Computing (Cat. No.99TH8469) (1999) p. 57.\n- [77] R. Haas et al. , The Einstein Toolkit (2022).\n- [78] C. J. Kruger, W. C. G. Ho, and N. Andersson, Seismology of adolescent neutron stars: Accounting for thermal effects and crust elasticity, Phys. Rev. D 92 , 063009 (2015).\n- [79] N. Andersson, K. D. Kokkotas, and B. F. Schutz, A new  \nnumerical approach to the oscillation modes of relativistic stars, Mon. Not. R. Astron. Soc. 274 , 1039 (1995)."}