Effects of nuclear matter properties in neutron star mergers

The dynamics in mergers of binary neutron star (BNS) systems depend sensitively on the equation of state (EOS) of dense matter. This has profound implications on the emission of gravitational waves (GWs) and the ejection of matter in the merger and post-merger phases and is thus of high interest for multi-messenger astronomy. Today, a variety of nuclear EOSs are available with various underlying microphysical models. This calls for a study to focus on EOS effects from different physical nuclear matter properties and their influence on BNS mergers. We perform simulations of equal-mass BNS mergers with a set of 9 different EOSs based on Skyrme density functionals. In the models, we systematically vary the effective nucleon mass, incompressibility, and symmetry energy at saturation density. This allows us to investigate the influence of specific nuclear matter properties on the dynamics of BNS mergers. We analyze the impact of these properties on the merger dynamics, the fate of the remnant, disk formation, ejection of matter, and gravitational wave emission. Our results indicate that some aspects of the merger are sensitive to the EOS around saturation density while others are sensitive to the behavior towards higher densities, e.g., characterized by the slope of the pressure as a function of density. The detailed density dependence of the EOS thus needs to be taken into account to describe its influence on BNS mergers.

The nuclear EOS is a key ingredient in BNS merger simulations.Experiments (see, e.g., Danielewicz et al. 2002;Tsang et al. 2012;Lattimer & Lim 2013;Le Fèvre et al. 2016;Roca-Maza et al. 2015;Russotto et al. 2016;Birkhan et al. 2017;Adhikari et al. 2021Adhikari et al. , 2022) ) and theoretical calculations based on chiral effective field theory (see, e.g., Hebeler et al. 2013;Tews et al. 2013;Lynn et al. 2016;Drischler et al. 2019Drischler et al. , 2020;;Keller et al. 2023) constrain the EOS up to around nuclear saturation density  0 = 0.16 fm −3 .However, for densities inside NSs (two times  0 or higher) the EOS represents one of the largest sources of uncertainties in numerical models of BNS mergers.Observations of NSs can provide excellent constraints of cold matter at supra-nuclear densities.The discovery of two-solar-mass neutron stars (Demorest et al. 2010;Antoniadis et al. 2013;Fonseca et al. 2021) has significantly reduced the uncertainties in the neutron star mass-radius relation.These measurements suggest that in order to support neutron stars with such high mass, the EOS cannot be too soft.Furthermore, observations from NASA's NICER mission can constrain the mass and radius of pulsars simultaneously (Riley et al. 2019;Miller et al. 2019;Riley et al. 2021;Miller et al. 2021) leading to novel constraints on the EOS (see, e.g., Raaĳmakers et al. 2019;Dietrich et al. 2020;Legred et al. 2021;Raaĳmakers et al. 2021;Huth et al. 2022).Moreover, multi-messenger observations of binary compact object mergers like GW170817 / AT2017gfo can be used to probe the EOS at finite temperatures and even higher densities.However, to extract constraints from the observations of the merger of compact objects, detailed numerical relativity hydrodynamics simulations exploring the influence of many different EOSs are necessary.
The EOS directly influences the dynamics of the merger and determines the fate and structure of the remnant and the GW emission in the post-merger phase.Shock heating and tidal torques can eject up to 10 −2  ⊙ of material on dynamical timescales (1 − 10 ms) (Rosswog et al. 1999;Korobkin et al. 2012;Hotokezaka et al. 2013a;Rosswog 2013;Sekiguchi et al. 2015Sekiguchi et al. , 2016)).The masses of the very neutron-rich tidal ejecta and moderately neutron-rich shock-heated ejecta depend sensitively on the EOS.After the coalescence, a large fraction of the material ejected from the NSs stays gravitationally bound and forms an accretion disk around the central object.The density in the disk is typically between 10 10 − 10 13 g cm −3 so in a regime where the EOS is better constrained.However, the disk mass, structure, and composition are determined by the merger dynamics.Furthermore, as long as the remnant does not collapse to a black hole (BH) it interacts with the disk.Therefore, the dynamics and mass ejection in the disk phase are indirectly influenced by the nuclear EOS.In the early disk phase (∼ 100 ms after the merger), mass is ejected by oscillations in the deformed remnant (Bauswein et al. 2013;Radice et al. 2018;Nedora et al. 2019Nedora et al. , 2021c)).Furthermore, neutrinos emitted from the hot disk remnant interface can be reabsorbed in the outer layers of the disk supporting mass ejection (Perego et al. 2014).After ∼ 1 s the disk cools down to a point where further cooling by neutrino emission becomes inefficient.As a result, 10% -50% of the disk is ejected due to heating by an effective viscosity originating from the magnetorotational instability (MRI) (see, e.g., Fernández & Metzger 2013;Just et al. 2015;Fujibayashi et al. 2018Fujibayashi et al. , 2020;;Fahlman & Fernández 2022;Hayashi et al. 2022).The ratio of dynamical ejecta to disk outflow mass depends mainly on the fate of the central object.If the remnant NS does not collapse to a BH the mass of the disk ejecta is one to two orders of magnitude larger than that of the dynamical ejecta (Fujibayashi et al. 2020).However, if a BH is formed within 10 − 20 ms after the merger the ratio is closer to one (Fujibayashi et al. 2023).
Many aspects of the nuclear EOS influence the merger dynamics such as the detailed density dependence of the pressure, thermal effects, and the dependence on the composition.Many studies aim to quantify nuclear physics uncertainties in numerical simulations by employing a small sample of different EOS models Bovard et al. (2017); Sekiguchi et al. (2016); Radice et al. (2018); Nedora et al. (2021b), thereby varying many of the above-mentioned aspects at the same time.While this approach is useful to reveal correlations with the general stiffness of the EOS, the impact of its specific features cannot be studied this way.Only a few works have focused on the dependence of BNS mergers on individual aspects of the EOS.For example, Most & Raithel (2021) investigated the impact of the slope of the symmetry energy on the post-merger dynamics and the temperature dependence has been investigated by changing the thermal pressure independently from the cold EOS in several works (see, e.g., Bauswein et al. 2010;Hotokezaka et al. 2013a;Raithel et al. 2021;Fields et al. 2023).
In this work, we perform 3D general-relativistic simulations of merging NSs based on various EOSs with different nuclear matter properties.The employed EOS models are based on the work of Yasin et al. (2020) using Skyrme energy-density functionals combined with the liquid-drop model and are created with the SROEOS code (Schneider et al. 2017(Schneider et al. , 2018)).By adjusting the parameters of the Skyrme functional, each nuclear matter property can be varied independently.A similar approach was used in numerical simulations of supernovae (Schneider et al. 2019;Yasin et al. 2020;Andersen et al. 2021) and very recently also in BNS-merger simulations (Fields et al. 2023).This enables us to study the impact of the nucleon effective mass, the incompressibility, the symmetry energy, and the saturation density and energy separately.We specifically focus on the effective mass and incompressibility in this work.Furthermore, we systematically vary all the nuclear matter properties mentioned from the values of the LS220 (Lattimer & Swesty 1991) to the ones of the Shen EOS (Shen et al. 1998) as in the work of Yasin et al. (2020).
The paper is organized as follows.Section 2 describes our setup for performing BNS merger simulations and explains the postprocessing techniques used to analyze them.In Section 3, we introduce the method to create the EOS models as well as their key properties.The results of the simulations are reported in Sections 4 to 6, which describe the dynamics of the merger and post-merger phases, the post-merger gravitational wave emission, and the properties of the ejected material, respectively.Finally, we summarize our conclusions in Section 7.

SIMULATION SETUP
The simulations we employ to model BNS systems have been carried out employing the framework of the EinsteinToolkit suite (Haas et al. 2020;Löffler et al. 2012), itself based on the Cactus computational toolkit (Goodale et al. 2003).
To handle general relativistic hydrodynamics we employ the WhiskyTHC code (Radice & Rezzolla 2012;Radice et al. 2014a,b).Modeling the NS matter as a perfect fluid, it solves Euler's equations for the balance of energy and momentum, coupled to conservation laws for the neutron and proton densities: In Eqs. ( 1) and (2),  p,n are respectively the proton and neutron number densities,   is the fluid four-velocity and   is the fluid stress-energy tensor.In this formulation, due to charge neutrality, the electron fraction  e is directly related to the proton number density, i.e.,  e ≡  e /(  +   ) =  p /(  +   ).
Equations ( 1) and ( 2) include source terms due to neutrino interactions:  p,n are the net lepton number exchange rates, due to the absorption and emission of electron neutrinos and antineutrinos, while  is the net energy deposition rate due to the absorption and emission of (anti-) neutrinos of all flavors.Neutrino emission and absorption are modeled with a leakage (Galeazzi et al. 2013;Neilsen et al. 2014) scheme and the so-called "M0" scheme of Radice et al. (2018), respectively.These are both "grey" (i.e., energy integrated) schemes evolving three neutrino species: electron neutrinos  e , electron antineutrinos νe and heavy neutrinos  x , which accounts for all others (anti-) neutrino flavors.Furthermore, the M0 scheme evolves the distribution function of neutrinos on a ray-by-ray grid, which we set up consisting of 2048 rays uniformly spaced in latitude and longitude with a radial resolution Δ ≈ 244 m.
In our setup, WhiskyTHC employs a finite-volume scheme for the discretization of the hydrodynamic quantities.The scheme reconstructs primitive variables with the fifth-order MP5 method (Suresh & Huynh 1997), from which numerical fluxes are computed with the HLLE flux formula (Harten et al. 1983), augmented with a positivitypreserving flux limiter (Hu et al. 2013;Radice et al. 2014a) in order to handle the transition to vacuum regions (which we fill with an atmosphere of density  atmo ≈ 1.24 × 10 3 g cm −3 ).The hydrodynamics is coupled to a dynamically evolved spacetime.Einstein's equations are written in the BSSNOK formulation (Shibata & Nakamura 1995;Baumgarte & Shapiro 1998;Brown 2009), and discretized with fourth-order finite-differences stencils by the McLachlan code (Brown et al. 2009).We furthermore choose the "1+log" and "Gamma-driver" gauge conditions (see, e.g., Baumgarte & Shapiro 2021).
The time evolution is performed by the strong-stability-preserving RK3 integrator (Gottlieb & Shu 1998) using a method-of-lines scheme.The time step is determined by the Courant-Friedrich-Lewy (CFL) criterion taking the speed of light as the maximum propagation speed.The CFL factor is chosen to be 0.15.
The mesh for our simulations is handled by the Carpet code (Schnetter et al. 2004), which employs a "moving boxes" Berger-Oliger adaptive mesh refinement scheme (Berger & Oliger 1984;Berger & Colella 1989).We employ a Cartesian-coordinates computational domain made of 7 refinement levels.The resolution on the finest level is 0.128 M ⊙ ≈ 189 m and the full grid spans 1024 M ⊙ ≈ 1512 km in each direction.To reduce the computational cost, we only evolve the  ≥ 0 part of the domain, and impose reflecting boundary conditions at  = 0.

Postprocessing methods
We define the NS as the region where the rest mass density is larger than 10 13 g cm −3 .Since the gauge is evolved dynamically during the simulation we need to compare simulation data in a gaugeindependent way.For this, we divide all grid cells in a snapshot of a simulation into 100 uniformly spaced density bins up to 10 15 g cm −3 .Any quantity of interest  is then averaged in each bin, weighted by the conserved rest mass density √ : where  bin is the rest mass density of the bin and ∫  bin d 3  indicates the volume integral over all cells in the density bin.
Furthermore, we make use of the complex azimuthal mode decomposition given by (Paschalidis et al. 2015;East et al. 2016a,b;Radice et al. 2016;Nedora et al. 2021c) to study the deformation of the remnant.We estimate the relevance of neutrino heating based on the net timescale of the neutrino heating   .It is given by the ratio of the conserved internal energy density  =  −  and the local net neutrino heating rate   =  M0  −  Leak  : where  is the lapse function.
The properties of the ejecta are extracted on a detection sphere with a radius of 300 km.The Bernoulli criterion is used to determine whether a fluid element is unbound.It is defined by where ℎ ∞ = lim ,→0 ℎ is the asymptotic specific enthalpy Foucart et al. (2021).Accordingly, the asymptotic velocity of an ejected fluid element is defined as . This criterion is less restrictive than the geodesic criterion −  > 1, which does not take the ejecta's thermal and binding energy into account.Typically, the ejection rate of matter meeting the geodesic criterion stops after ∼ 10 ms post-merger and roughly corresponds to the fluid elements ejected by dynamical processes (see Section 5 for further discussion).Therefore, we define fluid elements fulfilling the geodesic criterion as dynamical ejecta (the same convention is used by Nedora et al. 2021c;Combi & Siegel 2023).Ejecta that only satisfy the Bernoulli, but not the geodesic criterion are defined as disk ejecta.The extraction radius of 300 km ensures that the remnant accretion-disk system is fully contained (the disk typically extends to ∼ 100 km) but is small enough for different bursts of ejecta to be resolved.If instead a much larger extraction radius is chosen, the time at which ejected material is registered depends more on the ejecta velocity than on the time of ejection, which hinders the identification of different ejecta components.For an in-depth discussion of different ejection criteria and their impact see, e.g., Foucart et al. (2021).
To extract GW waveforms and spectra we employ the Newman-Penrose scalar Ψ 4 (Newman & Penrose 1962), following Sect.5 of Hinder et al. (2013) for its practical implementation.The second time derivative of the GW strain polarization components ℎ + and ℎ × can be related to Ψ 4 by: where Ψ 4 is expanded in weighted spherical harmonics of spin weight -2, −2   (, ).The expansion coefficients   are extracted at multiple finite coordinate radii inside the simulation domain and extrapolated to null infinity along outgoing radial null geodesics.From them, the strain components ℎ + and ℎ × are obtained by performing the time integration with the "fixed frequency integration" (FFI) method (Reisswig & Pollney 2011).Finally, the power spectral density (PSD) of the signal is given by with the frequency-domain strain components

Overview of the BNS models
We perform simulations for 9 EOSs in Table 2 (see Section 3 for more details).Each model initially consists of two irrotational identical  = 1.365 ⊙ NSs on quasi-circular orbits with an initial separation of 45 km.This combination corresponds to a chirp mass of 1.188 ⊙ and is thus compatible with the GW source of GW170817.This orbital separation corresponds to an inspiral phase of 2 − 3 orbits before the merger.The initial data for all the selected simulations are constructed using the spectral elliptic solver LORENE (Gourgoulhon et al. 2001).In the construction of the initial data, the minimum temperature slice of the EOS table is used and the composition is determined by neutrinoless beta-equilibrium.An overview of the simulations is compiled in Table 1.Each simulation is run for at least 40 ms post-merger.The only exceptions are the models LS175 † and LS220 † .In LS175 † , a BH is formed Table 1.Overview of all EOS models and their key results, including the simulated time  end , time until collapse to a BH  BH , disk mass at the end of the simulation  disk , masses of different ejecta components  ej , frequencies of the two most prominent peaks in the post-merger GW spectrum  1 and  2 , radius  NS and reduced tidal deformability Λ of the NSs in the initial data of our simulations.The EOS models and labels are based on Yasin et al. (2020) almost immediately after the merger so barely any mass is ejected.In LS220 † , the remnant collapses after ∼ 8 ms which is enough time to eject matter and form an accretion disk.However, the mass ejection and GW emission stop completely after 30 ms.In all other models, the remnant is a massive NS which stays stable for the duration of the simulation and is surrounded by an accretion disk.

EOS MODELS
In BNS merger simulations, the influence of the EOS on the dynamics of the merger is often discussed only in terms of the softness or stiffness of the EOS, as encapsulated, e.g., in the radius of a cold non-rotating NS or the reduced tidal deformability of the initial binary system (see, e.g., Bauswein et al. 2016;Takami et al. 2015;Rezzolla & Takami 2016;Nedora et al. 2021a).This approach has been successful in describing the peak frequencies of the post-merger GW emission, typically within a 10% error, as well as the threshold mass to prompt collapse (see, e.g., Bauswein et al. 2021).The same approach has been extended to features of the remnant-disk system and the dynamical ejecta (Dietrich & Ujevic 2017;Krüger & Foucart 2020;Nedora et al. 2021c;Henkel et al. 2023).In this work, in order to enable a more detailed description of the EOS' features, we employ the expansion parameters of the energy density of nuclear matter at nuclear saturation density.
The energy per particle of nuclear matter   at zero temperature is a function of number density  =   +   and isospin asymmetry Here,   and   are the neutron and proton number densities, respectively, and  =   / is the proton fraction, which is equal to the electron fraction  e due to charge neutrality.The energy per particle can be expanded around saturation density  =  0 and symmetric matter  = 0: where  is the baryonic energy density and  = is the relative density difference to saturation density.The expansion coefficients (called nuclear matter properties) for symmetric nuclear matter with  = 0 are the binding energy  and the incompressibility .The third term in Eq. ( 9) depends on the isospin asymmetry and is called symmetry energy.It is defined via the second derivative of the energy per particle with respect to  and can be expanded as with the symmetry energy at saturation density  sym and the slope parameter , which define two nuclear matter properties for neutron matter ( = 1).In this quadratic approximation, the symmetry energy is the difference between the energy per particle of neutron matter and symmetric nuclear matter: It has been shown in microscopic calculations that the quadratic expansion of the energy per particle provides a good description of the EOS around saturation density (see, e.g., Drischler et al. 2016;Wellenhofer et al. 2016;Somasundaram et al. 2021;Keller et al. 2023), with non-quadratic contributions arising mainly from the kinetic energy density that is fully included in EOS functionals.
As in Yasin et al. (2020), we use the open-source SROEOS code (Schneider et al. 2017(Schneider et al. , 2018)), which is based on a Skyrme density functional plus liquid drop model to describe the low-density EOS.The nuclear matter parameters , ,  sym , and  0 are systematically varied by adjusting the parameters of the Skyrme density functional while the effective mass  * at saturation density is varied directly.To explore the impact of these parameters we use nine EOS tables.Eight of the tables were created with the SROEOS code, six of which were developed and used in Yasin et al. (2020) while two are used for the first time in this work.Finally, the ninth EOS is the Shen EOS (Shen et al. 1998).The nuclear matter properties of these EOS models are summarized in Table 2.
With each of the EOS models, we systematically vary one nuclear matter property at a time.Following Yasin et al. (2020), we use the model LS220 † as the fiducial model which is based on the same Skyrme parametrization as the LS220 EOS.Its incompressibility parameter is  = 220 MeV and its effective mass is not density-dependent and simply given by the neutron mass  * =   = 939.5654MeV.First, we vary the  parameter to 175 and 255 MeV, resulting in the EOS LS175 † and LS255 † , respectively.These values represent the upper and lower bounds predicted by chiral effective field theory calculations (Hebeler et al. 2011;Drischler et al. 2016Drischler et al. , 2019)).Second, we change the nucleon effective mass at saturation density to  * = 0.8  ( * 0.8 ) and  * = 0.634  ( * S ).The value of  * = 0.634  is chosen because it is used in the Shen EOS.Finally, starting from  * S , we also change the remaining nuclear matter properties to match the values of the Shen EOS, starting with  = 281 MeV, followed by the symmetry energy  sym = 36.9MeV, and the value of the saturation density  0 = 0.145 fm −3 and the binding energy  = 16.3MeV.This results in the EOSs ( * ) S , ( *   sym ) S , and SkShen, respectively.We include the original Shen EOS for comparison.
Figure 1 shows the mass-radius relation of cold NSs for all of Table 2. Nuclear matter properties for all employed EOSs models, following the notation in Yasin et al. (2020) Note that  * is the effective nuclear mass at saturation density.Also given are the pressure and the dimensionless slope of the pressure at saturation density for  = 0 and   = 0.1.the EOSs.Changing the incompressibility has a large influence on the radius of NSs with high masses (compare LS175 † , LS220 † , and LS255 † or  * S and ( * ) S in Fig. 1).The effective mass impacts intermediate to high mass NSs but has less of an effect at high masses compared to the incompressibility.Varying the symmetry energy mostly has a large effect on the radii of low mass NSs.The horizontal dashed line marks the NSs with a mass of 1.365 ⊙ .These are the radii of the NSs in the initial data of the simulations (see Section 2.2).
Figure 2 shows the pressure versus density at  = 0 and  e = 0.1 for all the employed EOSs.The pressure is related to the derivative of the energy density with respect to the particle density.Since the pressure of symmetric matter vanishes at saturation density, the slope of the pressure at  = 0 is governed by the incompressibility.Increasing  leads to larger pressure at densities above  0 , but also smaller pressures below  0 (compare LS175 † , LS220 † , and LS255 † in Fig. 2).Changing the symmetry energy also changes the slope parameter , which characterizes the pressure of neutron matter.Thus, it has a relatively large effect close to saturation density but smaller at higher densities (compare ( * ) S and ( *   sym ) S in Fig. 2).Finally, the effective mass influences the degeneracy pressure from the kinetic energy density, which for our Skyrme density functional is given by (Lattimer & Swesty 1991) Since the effective mass enters the degeneracy pressure inversely, decreasing the effective mass increases the cold pressure.Furthermore, for the EOSs used here reducing the effective mass implicitly increases the  parameter due to the way that the Skyrme parameters are determined.Therefore, the pressure is increased for the full range of densities in the NS ( > 10 14 g cm −3 ).While increasing  and reducing  * both generally increase the pressure (i.e., the stiffness of the EOS), the two parameters have different effects around saturation density.This can be visualized by considering the pressure and the dimensionless slope of the pressure  above,  only affects the slope of the pressure at saturation density but not the value of the pressure itself.Reducing the effective mass (LS220 † ,  * 0.8 ,  * S ), increases the pressure at saturation density and its slope but to a lesser degree than increasing .This is due to an increase in the next-order term in the symmetry energy  sym .Finally, ( *   sym ) S , SkShen, and Shen exhibit a higher pressure at saturation density, but their dimensionless pressure slope is comparable to that of  * S .The radii of the initial NSs are shown in the lower panel of Fig. 3.They are mostly correlated with the pressure at saturation density.
The total energy density and pressure are often divided into a cold and thermal part (indicated by the subscripts cold and th, respectively).

𝑃(𝑛
where  th (,  = 0) = 0.The temperature dependence is often described by the thermal index Γ th given by (see, e.g., Bauswein et al. 2010) For a non-interacting gas of non-relativistic fermions with densitydependent effective mass, the thermal index is determined by the density dependence of the effective mass (see, e.g., Constantinou et al. 2015;Lattimer & Prakash 2016;Keller et al. 2021): Chiral effective field theory calculations have found that Eq. ( 15) is a good approximation for the full thermal index via Eq.( 14) (Carbone & Schwenk 2019;Keller et al. 2021Keller et al. , 2023)).In Fig. 4, the thermal index inside the NSs is shown versus density for the three values of the effective mass used in our simulations, and for the Shen EOS.For the non-relativistic Skyrme density functional based EOSs, the thermal index is the same using Eq. ( 14) or Eq. ( 15).Therefore, lowering the effective mass at saturation density increases the thermal index.We also note that this behavior is due to the parametrization of the density-dependence of the effective mass, which decreases with increasing density, leading to an increasing thermal index.However, the results from Carbone & Schwenk (2019); Keller et al. (2021Keller et al. ( , 2023) ) show that Γ th () has a maximum and starts to decrease at higher density, which originates from a minimum of the effective mass around saturation density (see also Huth et al. (2021)).Note that the thermal index of the Shen EOS does not follow the expression from Eq. ( 15) due to relativistic mean field functional used.
Several previous works investigated thermal effects in binary compact object mergers in isolation by changing the thermal pressure independently from the cold EOS (see, e.g., Bauswein et al. 2010;Hotokezaka et al. 2013a;Raithel et al. 2021;Fields et al. 2023).The setup of Fields et al. (2023) is also based on variations of the effective mass with the SROEOS code and is very similar to ours.However, by employing a functional with more free parameters than ours, they also constrain the  parameter and the pressure at four times  0 .Therefore, the variations of the effective mass in their work only change the thermal pressure and do not affect the EOS at zero temperature, while in our EOS function both the cold and the thermal pressure are changed by varying the effective mass.

MERGER AND POST-MERGER DYNAMICS
Many of the observable features of a BNS merger depend sensitively on the dynamics during the first 10−20 ms after the merger.Therefore, we investigate this early post-merger phase in this section.We specifically focus on the effect of the EOS on the collapsing behavior, the initial contraction and deformation of the massive NSs, the role of thermal effects, and the formation of the disc.

Evolution of the central object
Figure 5 shows the time evolution of the maximum density inside the remnant NS.As the two cores of the original NSs combine after the merger, the maximum density rises.Consequently, the pressure in the NS core increases as well.If the pressure given by the EOS is relatively low (i.e., for a soft EOS) a prompt collapse occurs, as is the case of model LS175 † .The other models initially form a massive NS remnant.Note that the pressure at high densities (several times the saturation density) is relevant for the collapse.In the context of the quantities introduced in Section 3 and shown in Fig. 3, the pressure at high density is determined by a combination of the pressure and its slope at saturation density.
All models except LS175 † do not promptly collapse.In these cases, the cores of the initial NSs bounce off each other initially.Thus, in the center of the newly formed massive NS, two density maximums are present and form an oscillating bar shape.The oscillations can be seen in the maximum density in Fig. 5.It typically takes 10−20 ms for the NS cores to merge.During this process, the central density and consequently the pressure in the NS core keeps increasing.The contraction ends either due to the delayed formation of a BH or when the pressure becomes large enough to balance the central object's self-gravity.In this case, a meta-stable NS remains in the center.For LS220 † , several oscillations with increasing amplitude occur before the central object collapses to a BH, at roughly 8 ms after the merger.Because these models exhibit a larger slope of the pressure as a function of density, the pressure rises faster as the central density increases.In Fig. 5, one can see that the central density of models with larger incompressibilities (LS255 † , ( * ) S , ( *   sym ) S ) stops rising already after 10−15 ms while it keeps increasing for more than 20 ms in  * 0.8 ,  * S .However, this trend does not hold for the models SkShen and Shen in which the contraction of the remnant continues for longer times as well.The central density at the end of the non-collapsing simulations seems to be determined by the pressure at 6 − 7 × 10 14 g cm −3 (compare Fig. 2 with Fig. 5).
During this initial phase of contraction, the remnant is highly deformed and oscillates, where the dominant deformation mode is a  = 2 bar-shaped deformation.However, after some time, the  = 1 deformation becomes dominant.This change influences the gravita- tional wave emission, disk formation, and matter ejection (Bernuzzi et al. 2014;Kastaun & Galeazzi 2015;Paschalidis et al. 2015;East et al. 2016a,b;Lehner et al. 2016;Radice et al. 2016;Nedora et al. 2019Nedora et al. , 2021c)).In Fig. 6, the  = 1 and 2 deformation calculated with Eq. ( 4) is shown.The models undergo the change from a dominating  = 2 to a dominant  = 1 mode at different times.In LS255 † , ( * ) S , and ( *   sym ) S , the transition occurs already after roughly 9, 15, and 17 ms, respectively, which is significantly earlier compared to the other models.This roughly matches the time at which the contraction of the massive NSs in these models stops, indicating that the stability of the initial bar-shaped mode might be linked to the contraction of the remnant and therefore to the slope of the pressure.

Thermal effects
During the merger, the initially cold NS matter is heated to several tens of MeV or more by shocks originating from the interface of the collision.The high-density NS cores remain comparatively cool because the shocks do not penetrate them.As the colder cores merge, the hot matter is redistributed into a ring shape at densities of approximately 1−5×10 14 g cm −3 .This evolution can be seen in Fig. 7, which shows temperature profiles inside the remnant for simulation  * 0.8 at approximately 0.5, 3, and 10 ms post-merger.Because the center of the NS is cold and very dense, the cold pressure dominates and thermal contributions to the pressure are not relevant in this region.However, the thermal pressure plays a significant role in the hot ring close to saturation density and even becomes larger than the cold pressure in the outer layers of the NS as can be seen in Fig. 8. Thermal effects are therefore especially important for disk formation and matter ejection but their impact on the post-merger gravitational wave emission is small (see, e.g., Bauswein et al. 2010;Hotokezaka et al. 2013b).
Generally, a stiffer EOS leads to less shock-heating due to two effects.Larger initial NS radii lead to an earlier merger because the NSs come into contact at larger separations.This leads to smaller orbital and radial velocities (see, e.g., Radice et al. 2020) and thus a less intense initial bounce and less shock-heating.Furthermore, with a softer EOS, the remnant will contract more before the first bounce occurs, which leads to stronger oscillations in the early post-merger phase.
Apart from the temperature, the thermal pressure depends on the thermal index which is determined by the density dependence of the effective mass in our EOS models.Reducing the effective mass at saturation density increases the thermal index and therefore the effectiveness of shock heating (Hotokezaka et al. 2013a).However, reducing the effective mass in our EOS setup also increases the  and  sym parameters and thus the cold pressure, which in turn leads to a less violent merger.Therefore, the effect of changing the thermal index cannot be observed independently in our simulations.
The lower panel of Fig. 5 shows the evolution of the average temperature inside the remnant as a function of time.In LS220 † the remnant is heated up much more than in the non-collapsing cases due to the increasingly violent oscillations.Because stiffer EOSs result in a less violent merger and reduced shock heating, non-collapsing remnants with larger maximum densities generally also have larger temperatures.However, there are some deviations from this general trend.
First, the average remnant temperature depends mostly on the outer layers of the NS (and thus on the EOS close to saturation density) since this is the region that is heated by shocks.Furthermore, the initial NS radii (and therefore the radial infall velocities of the NSs at merger) depend on the EOS at lower densities as well.The maximum density, however, is more sensitive to the pressure at high densities.This can be seen by comparing  * S and ( * ) S , which have the same pressure at saturation but differ at high density.Consequently, the maximum density is lower in ( * ) S but the average temperatures are very similar.Similarly, ( * ) S and ( *   sym ) S show a very similar evolution of the central density but the average temperature of ( *   sym ) S is lower because the two EOSs differ only at low densities.
Second, the thermal index influences the effectiveness of shock heating which can be seen by comparing LS255 † and  * S .The LS255 † EOS exhibits lower pressures for all densities compared to  * S (see Fig. 2).Therefore, the central density in LS255 † is higher than in  * S .Nonetheless, its average remnant temperature is lower because  * S has a higher thermal index which increases the efficiency of the shock heating.

Comparison of Shen and SkShen
The SkShen EOS matches the effective mass, incompressibility, symmetry energy, binding energy, and saturation density of the Shen EOS.The two EOS models are thus similar close to saturation density but differ more and more at higher densities.Furthermore, the thermal index of the Shen EOS is lower compared to SkShen.This is expected because matching the effective mass at saturation density does not reproduce the thermal index if the density dependence of the effective mass is different.
Due to the increased pressure at high densities, the central density in the model Shen is larger than in SkShen.The average remnant temperature, however, is larger in SkShen because it has a larger thermal index, resulting in more effective shock heating.Nonetheless, considering the large deviation of the two EOSs at high densities, the evolution of the remnants in the models SkShen and Shen are remarkably similar.This indicates, that the post-merger evolution is largely determined by the EOS around saturation density, even though the maximum density in SkShen reaches up to 3 times the saturation density.

Accretion disk evolution
We define the disk mass as the total baryon rest mass outside of the NS, i.e., the region with  <  disk = 10 13 g cm −3 .However, this definition for the NS surface is relatively arbitrary.Specifically, it includes the transition region between the disk and NS.Thus, we also calculate all disk-related quantities for  disk = 10 12 g cm −3 which excludes this transition region.Therefore, the second definition offers a more conservative estimate of which matter contributes to the disk.Figure 9 shows the evolution of the disk mass for all simulations except LS175 † .The left and right panels correspond to  disk = 10 12 and  disk = 10 13 g cm −3 , respectively.The disk masses at the end of the simulations are listed in Table 1.
The formation of the accretion disk depends on the EOS properties due to multiple effects.The most important factor for the disk mass is the fate of the merger remnant.In simulation LS175 † , the central object collapses to BH almost immediately after the merger, so no disk is formed, while the BH formation in LS220 † is delayed long enough for a disk to form.However, roughly half of the disk mass is swallowed upon collapse.For non-collapsing mergers, the relation of the disk mass to the EOS is more complicated.On one hand, the disk mass originating from the tidal disruption of the NSs is larger for stiffer EOSs.For a softer EOS, on the other hand, matter ejection due to shock heating is enhanced.Furthermore, the remnant is more compact and rotates faster (see, e.g., Bernuzzi 2020;Nedora et al. 2021c), which also increases the disk mass.In addition to the EOS at zero temperature, thermal effects have to be considered.A larger thermal index enhances the amount of matter ejected to the disk because of shock heating (Hotokezaka et al. 2013a).It is thus hard to find a correlation between the disk mass and specific nuclear matter properties.
At the end of the simulations that do not form a BH the disk mass is close to 0.25 ⊙ .The only exception is  * 0.8 with a disk mass of 0.14 ⊙ .The mass of the disk in LS255 † rises very fast in the first 10 ms after the merger.Furthermore, the disk is particularly massive in comparison to the other models, considering that LS255 † is comparable to  * 0.8 in terms of stiffness.This changes if the inner boundary of the disk is defined by  < 10 12 g cm −3 .
Even though the evolution of the disk mass and angular momentum varies for the different EOSs, the composition and structure of the disk are relatively similar in all non-collapsing simulations.The mass-weighted histogram of the electron fraction, entropy per baryon, and temperature in the disk at  = 30 ms is shown in Fig. 10.All histograms show a triple peak structure, typical for equal mass BNS mergers (Nedora et al. 2021c).The peak at the lowest entropy and electron fraction corresponds to the interface between the remnant and the disk.The dotted lines in Fig. 9 show the histograms for the matter at densities below 10 12 g cm −3 which lack this high temperature, low entropy, and low electron fraction tail.
The second peak at / ≈ 4 − 8  B and  e ≈ 0.1 − 0.2 represents the bulk of the disk.The disk bulk has densities between 10 10 − 10 12 g cm −3 and extends to roughly 100 km.For most simulations, the main peak in the electron fraction and entropy are located at  e = 0.16 and / = 7.5  B .The largest deviations from this trend are LS220 † and LS255 † .Since the central object in LS220 † is a BH, the hot and dense interface to the NS is missing, so there are no spiral arms in the disk.Furthermore, the strong neutrino irradiation from the central NS is also missing.Therefore, the disk has lower temperatures, entropies, and electron fractions.Due to the higher incompressibility of LS255 † , a larger fraction of the disk is located closer to the massive NS.This effectively shifts the average bulk entropy and energy toward lower values.The same is true for ( * ) S and ( *   sym ) S , however, to a lesser extent.
The third peak with high electron fractions and entropies is due to the matter outside of the main bulk of the disk.During the disk evolution, this matter becomes partially unbound.Thus, its composition matches that of the disk wind ejecta in Fig. 12 (see Section 5 for an extended discussion).

EJECTA PROPERTIES
The mass ejection in BNS mergers occurs through multiple channels.These can be separated into two categories: dynamic and disk (or secular) ejecta.The former are expelled within a few milliseconds after the merger, while the latter are ejected over timescales of 10 ms up to 10 s (see, e.g., Fujibayashi et al. 2023).
As described in Section 2.1, we define the dynamical ejecta as matter that fulfills the geodesic criterion, while matter that fulfills the Bernoulli criterion but not the geodesic criterion is associated with the disk ejecta.The same division is used in Nedora et al. (2021c); Combi & Siegel (2023).Note, however, that this division scheme is only approximate as the transition from dynamical ejecta to disk ejecta happens continuously between 5-10 ms after the merger.Furthermore, we separate the dynamical ejecta into the tidal and shock-heated components.For this, we simply count dynamically ejected fluid elements toward the tidal component if their electron fraction is smaller than 0.1 and toward the shock-heated component otherwise.The classification into dynamical and disk ejecta and tidal and shock-heated ejecta is made on the detection sphere at 300 km radius.Note that this distinction based on the electron fraction only works, if the ejecta extraction radius is small enough.At later times, neutrino captures increase the electron fraction of the tidal ejecta making it more difficult to distinguish them from the shock-heated ejecta.
In    fluid elements are clearly visible in the lower left corner of both panels (i.e., at low  e , early times, and in the equatorial direction).However, the distinction between the shock-heated and disk ejecta is less clear.We observe roughly separate bursts of mass ejection at  −  merg ≈ 8, 14, 20, . . .ms in the left panel.The first burst is the shock-heated ejecta, while the later bursts are caused by the oscillations following the initial bounce of the newly formed massive NS.Especially in the case of  * S , the division based on the geodesic criterion does not properly separate the first burst from the subsequent ones and slightly underestimates the shock-heated ejecta mass.
Figure 12 shows the accumulated mass histogram of the ejecta electron fraction, entropy per baryon, and asymptotic velocity.Solid lines show the dynamical ejecta and dotted lines the disk ejecta.If the tidal ejecta do not experience significant shock heating, they consist of almost pure cold NS matter ( e < 0.1) and have low entropies (/ < 10 B ).This part forms a peak at low electron fractions and entropies in Fig. 12.However, a fraction of the tidally removed matter can be reprocessed by shocks and thus counts towards the shock-heated ejecta component in our classification scheme.In ( *   sym ) S , SkShen, and Shen, this is the main part of the tidal ejecta.
The shock-heated ejecta reach much higher temperatures and entropies compared to the tidal ejecta.Therefore, positron and electronneutrino captures increase the electron fraction of the matter to  e ≈ 0.1 − 0.4.A small part of the shock-heated ejecta is responsible for the high-velocity and high-entropy tail in the right panels of Fig. 12.These ejecta might lead to additional observable features in the kilonova and is larger for softer EOSs (see, e.g., Dean et al. 2021, and references therein).
The evolution of the mass ejection rate and the integrated total ejected mass for the disk component is shown in Fig. 13.On the timescales of our simulations, the mass ejection in the early post-merger phase is primarily driven by the oscillating double-core structure in the massive NS.With each bounce, matter becomes unbound as the central density reaches a minimum (Nedora et al. 2019(Nedora et al. , 2021c;;Combi & Siegel 2023).The left panel of Fig. 11 shows several such bursts of matter ejection after  ≈ 10 ms.Furthermore, energy deposition by neutrino absorption and scattering (Dessart et al. 2009;Perego et al. 2014;Just et al. 2015;Radice et al. 2018) enhances the ejection matter from the early post-merger remnant.In LS220 † , the hypermassive NS collapses and thus no longer provides an efficient mechanism for the ejection of matter through neutrino irradiation and the oscillation of the NS.Thus the ejection of matter stops after  > 20 ms.For all non-collapsing models, however, the ejection of matter is still ongoing at the end of the simulation and viscous effects will eject a 10%−50% of the disk mass on timescales of seconds (see, e.g., Fernández & Metzger 2013;Just et al. 2015;Fujibayashi et al. 2018Fujibayashi et al. , 2020;;Fahlman & Fernández 2022;Hayashi et al. 2022).Thus, the total disk ejecta mass is not converged in our simulations.
We find that the mass ejection in the disk phase is more effective while the  = 2 bar-shaped deformation of the remnant drives two spiral arms into the disk (see Fig. 6).The model LS255 † exhibits a period of low mass ejection after ∼ 15−20 ms that lasts approximately 20 ms.A similar feature is visible for ( *   sym ) S after ∼ 25 ms, lasting approximately 10 ms.Both correlate well with periods of low  = 2 deformations of the remnant with a delay of ∼ 10 ms which is roughly the time it takes ejecta to travel to the detection radius at 300 km.Later, a one-armed spiral wave, driven by the  = 1 deformation, appears, resulting in an increased mass ejection rate.Similar results were found by Nedora et al. (2019Nedora et al. ( , 2021c)).
Directly above the massive NS the neutrino heating is very strong.However, the baryonic matter density in the polar direction is low, so the neutrino heating unbinds only a small amount of matter in this region.This matter is ejected toward the polar direction with relatively high velocities and electron fractions.
The neutrino heating in the diagonal direction enhances the mass outflow significantly and increases its electron fraction.This is visible in Fig. 11: in the right panel, the angles  ≈ 45 − 75 • are marked by dotted lines and encompass the largest ejecta component (the roughly horizontal bar-shaped region).In the left panel, this component can be seen after 15 − 20 ms at  e ≈ 0.3.This part of the ejecta forms a peak at  e ≈ 0.3 in Fig. 12 for all models except LS220 † , where it is missing due to the lack of neutrino emission from the collapsed massive NS. Figure 11 also shows how this component evolves.With time the disk becomes thinner and more transparent to neutrinos.As a result, the neutrino-irradiation intensifies and the outflow becomes more proton-rich and moves closer to the equatorial plane.
Initially, the oscillations of the massive NS ejects matter also in the equatorial direction.The dense bulk of the disk shields the matter in the equatorial plane outside the disk from neutrino irradiation.Therefore, the equatorial ejecta have lower electron fractions ( e ≈ 0.2 − 0.25).Furthermore, the energy deposition by neutrinos is reduced, so the amount of mass ejected in this direction decreases continuously as the oscillations of the massive NS die down.

Correlation of ejecta masses to EOS properties
The masses of the ejecta components at  −  merg = 40 ms are given in Table 1 for all simulations.In the following analysis, we will focus on the correlation of the dynamical ejecta masses with the properties of the employed EOSs.Since the mechanism of ejection is entirely different for the tidal and shock-heated dynamical ejecta components, it is not surprising that they depend differently on the properties of the EOS.This is visible in Fig. 14, which shows the masses of the tidal and shock-heated ejecta versus the dimensionless slope of the pressure and the absolute pressure at saturation density of the corresponding EOS (as defined in Section 3), respectively.It is well-known in the literature that the shock-heated ejecta mass is correlated with the EOS stiffness (see, e.g., Dietrich et al. 2017; Masses of the tidal and shock-heated ejecta components as a function of the dimensionless slope of the pressure and the pressure at saturation density, respectively.Radice et al. 2018;Nedora et al. 2021a;Henkel et al. 2023).This is often quantified by properties of cold NSs (like the reduced tidal deformability), which depend on the full range of densities present in NSs and thus also on the slope of the pressure at saturation density.We find that the shock-heated ejecta mass is correlated with pressure at saturation density and entirely unaffected by the slope of the pressure at saturation density.Specifically, we find that pairs of models that only differ in the incompressibility (LS220 † and LS255 † or  * S and ( * ) S ) exhibit almost the same amount of shock-heated ejecta.This points towards the shock-heated ejection mechanism being most sensitive to the EOS close to the saturation density and entirely independent of the high-density part of the EOS.
The mass of the tidal ejecta varies between 10 −6 − 4.3 × 10 −3  ⊙ .Even though it is relatively small in comparison to the other ejecta components, it plays a significant role because it can produce actinides and fissioning isotopes due to its very low electron fraction.This is important for galactic chemical evolution and the kilonova light curve, which might show signatures of fission reactions if very heavy elements are produced.Furthermore, the tidal ejecta component is expected to be larger in the merger of asymmetric binaries.In the models ( *   sym ) S , SkShen, and Shen, the tidal ejecta are caught by shocks early on, so almost no material with  e < 0.1 remains.While this implies that almost no tidal ejecta component exists in these models by our definition, it does not mean that no matter is ejected by tidal torques.For the other models, the tidal ejecta mass shows a correlation with the dimensionless slope of the pressure and is significantly larger for LS255 † and ( * ) S compared to the other models.Figure 15 shows the distribution of the electron fraction in the equatorial plane for LS220 † ,  * S , LS255 † , and ( * ) S .The tidal arms are visible as the low electron-fraction (i.e., red) regions, which are larger for models with larger slopes of the pressure.
The fact that a larger slope of the pressure influences the tidal ejecta mass might be related to the fact that it affects the distribution of the mass inside the NS.For a large slope of the pressure with respect to the density, the central density in the NS is reduced but at the same time, the density in the outer layers is increased.This effect can be seen by comparing the mass distribution in the NSs in the initial data of our simulations.Figure 16 shows the spatial mass distribution d d = 4 2 versus the radius for the models LS175 † , LS220 † , and LS255 † .For the sake of comparison, the radius is normalized to the outermost radius of the NS .Since there is more mass located in the outer layers of the NSs, the ejection of matter by tidal forces is more efficient.

GRAVITATIONAL WAVE EMISSION
We extract the gravitational waveform as outlined in Section 2.1.Figure 17 shows the + mode of the GW strain for all simulations except LS175 † , since no significant post-merger signal is produced after the prompt collapse.We only show the  = 2,  = 2 mode as it is by far the most dominant mode.It shows the low-frequency pre-merger phase ( < 0), as well as the post-merger phase ( > 0).The latter consists of two periods.A transition period ( ≲ 3 ms), during which the system readjusts from the inspiral of two NSs to one rotating and oscillating massive NS, and the longer ring-down phase that follows, during which the amplitude gradually decreases (see, e.g., Baiotti 2019).In LS220 † , the GW signal stops abruptly after the collapse of the hypermassive NS.In LS255 † , ( *   sym ) S , and ( * ) S , the GW amplitude drops within the first 10-15 ms while it decreases more slowly within 20-30 ms in  * 0.8 ,  * S , SkShen, and Shen.This decrease in the amplitude is directly linked to the  = 2 deformation of the massive NS (see Fig. 6).As described in Section 4.4, the  = 2 deformation decreases much faster in the simulations with higher pressure slopes.Future detections of post-merger GWs could therefore be used to constrain the slope of the pressure at high densities.
The post-merger GW time-frequency spectrograms (upper panels) and the corresponding Fourier spectra (lower panels) are shown in Fig. 18.The thin lines in the Fourier spectra represent the full spectrum, while the thick lines are produced by excluding the inspiral GW signal.Initially, several frequencies are present but they decay within approximately 5 milliseconds (Takami et al. 2015;Rezzolla & Takami 2016).Afterward, the spectrum is dominated by a single frequency, often called  2 .
This frequency has been identified and studied in many previous works (Stergioulas et al. 2011;Bauswein & Janka 2012;Bauswein et al. 2012;Hotokezaka et al. 2013b;Takami et al. 2014Takami et al. , 2015;;Rezzolla & Takami 2016;Fields et al. 2023) and is attributed to the quadrupole mode of the NS (Bauswein & Janka 2012;Bauswein et al. 2012).The solid vertical lines show the maximum of the Fourier spectra, while the dotted lines follow the time-dependent maximum frequency.The following effects can be seen in the time-frequency spectrograms: • In LS220 † , the post-merger signal shows a "chirp-like" behavior (the peak frequency rises quickly) shortly before the remnant collapses to a BH.This is because the rotational velocity of the massive NS increases significantly as it collapses.
• In the first 10 ms of the post-merger GW emission, the  2 frequency can vary slightly.This is especially visible for the stiffest EOSs: ( * ) S , ( *   sym ) S , SkShen, and Shen.For those, the  2 frequency increases until  ≈ 5 ms and subsequently decreases again until  ≈ 10 ms.A similar effect is described by Rezzolla & Takami (2016).They determine the frequency during the transient phase separately from the ring-down phase and label it  2, .
• After the initial transient phase, a continuous shift of the  2 emission towards higher frequencies is visible for  * 0.8 and  * S .A similar but weaker increase is visible for  * S , SkShen, and Shen.In these models, the gravitational wave amplitude stays high for an extended period.The GWs carry away angular momentum, the NS contracts leading to a larger  2 (Maione et al. 2017).
The second most prominent peak is the so-called  1 peak (sometimes also called  − ), which always lies at lower frequencies than the  2 peak and disappears after ∼ 5 ms (Stergioulas et al. 2011;Takami et al. 2014Takami et al. , 2015;;Rezzolla & Takami 2016).Its origin has been attributed to the interaction of the  2 and the quasi-radial  0 mode (named  2−0 ) (Stergioulas et al. 2011) as well as the orbital motion of antipodal bulges rotating around the central remnant with a slower frequency (thus called  spiral ) (Bauswein & Stergioulas 2015).Depending on the remnant compactness, either one or both frequencies might be present (Bauswein & Stergioulas 2015;Bauswein et al. 2016;Rezzolla & Takami 2016;Maione et al. 2017;Kiuchi et al. 2020).We define  1 independently of its origin as the second highest peak with a frequency at least 400 Hz below the  2 peak.The extracted  1 peaks are marked by a dashed vertical line in Fig. 18.
The  1 and  2 frequencies for all simulations are listed in Table 1.They span from  2 = 2.24 kHz for Shen to  2 = 3.10 kHz for LS220 † and  1 = 1.66 kHz for SkShen to  1 = 2.14 kHz for LS220 † where softer EOSs generally produce larger frequencies.Decreasing the effective mass and increasing the incompressibility and symmetry energy lowers the pressure in the center of the NS.This decreases the density, which in turn reduces  2 (Bauswein & Janka 2012;Bauswein et al. 2012).By changing the nuclear matter properties to the values of the Shen EOS (i.e., the progression  * S , ( * ) S , ( *   sym ) S , SkShen), both peak frequencies approach those of the Shen simulation.
The Fourier spectra and the time-frequency spectrograms of models Shen and SkShen are very similar.Especially the position and width of the  1 and  2 peaks as well as the amplitude and time dependence of the spectrogram match almost perfectly.This implies, that the EOS-impact on the post-merger GW emission is well described by the nuclear matter properties of the EOS, while the details of the microphysics and models (i.e., relativistic mean field and Skyrme density functionals) only play a minor role.It might thus be possible to directly constrain the properties of nuclear matter with future detections of the post-merger GW spectra.Furthermore, this implies that the post-merger GW emission is mostly sensitive to the EOS around saturation density, since this is the region where SkShen and Shen match.
Many works have found that the  2 frequency is correlated with the properties of the EOS.Takami et al. (2015); Rezzolla & Takami (2016); Bauswein et al. (2016); Kiuchi et al. (2020) have provided fit formulae relating  1 and  2 with various properties of cold isolated NSs.We compare their fit formulae with our results.Figure 19 shows the residuals of the fits to our extracted frequencies.
All fit formulae predict the correct frequencies within an uncertainty of 5−10%.However, some trends in the difference can be identified.Fit formulae for  2 (solid symbols connected with solid lines) underestimate the value for  2 in LS220 † and overestimate it in LS255 † ,  * S , ( * ) S , ( *   sym ) S , SkShen, and Shen.The fit from Kiuchi et al. (2020) performs better, but also uses the most fit parameters.The fits for  1 do relatively well for the models with  * / < 1 but overestimate the frequency for LS255 † and LS220 † .Our set of simulations is too small to make quantitative predictions for the peak frequencies.However, we plan to increase the number of models in future works.Therefore, we will potentially be able to enhance the accuracy of universal relations in the future by incorporating additional parameters describing the EOS into the fits.

CONCLUSIONS
In this work, we systematically investigated the influence of the nuclear EOS on BNS mergers.To this end, we performed simulations of merging equal mass NS binaries while individually varying the nuclear matter properties in the employed EOSs from the fiducial model LS220 † to those of the Shen EOS.This study focused on the influence of the nucleon effective mass and the incompressibility parameter.We provide a comprehensive analysis of the remnant and disc dynamics, the post-merger GW spectrum, as well as the proper-ties of the dynamical and early disk ejecta.The key findings can be summarized as follows: • The incompressibility governs the slope of the cold pressure at saturation density.Therefore, it has a large impact when the density is much larger than saturation density, i.e., in the core of the merger remnant, which typically reaches three to five times saturation density.This has implications for the emission of gravitational waves, as well as for the fate of the remnant.Increasing the incompressibility decreases the compactness of the remnant core and therefore lowers the dominant post-merger GW frequency.For our choice of initial conditions, varying the incompressibility  within theoretical uncertainties results in a prompt collapse for low  and a NS remnant that is stable for the duration of the simulation for high .In the latter case, the steep rise of the pressure halts the contraction of the  20) and ( 23) (labeled RT), Kiuchi et al. (2020), Eq. (5.4) (labeled K), Bauswein et al. (2016), Eq. ( 2) (labeled B), and Takami et al. (2015), Eq. ( 25) (labeled T).Shown are the relative deviations of the fits to the measured frequencies Δ  =  NR −  fit .Solid symbols connected with solid lines represent fits for  2 , while crosses connected with dashed lines represent fits for  1 .
newly formed remnant shortly after collapse and therefore leads to less oscillations in the remnant compared to other non-collapsing simulations.
• The effective nucleon mass is important for both the pressure at zero temperature as well as thermal contributions in our EOS models.Lowering the effective mass increases the cold pressure at all densities inside the NS.Similar to the model with increased incompressibility, this results in a less compact remnant and a reduction of the peak frequency in the post-merger GW spectrum.However, a significant difference is that the pressure density dependence is less steep leading to a longer period of contraction and oscillations as imprinted in the post-merger GW amplitude.Another important aspect of the effective mass is its influence on the thermal pressure.Decreasing the effective mass results in a larger thermal index.At the same time, it also reduces shock heating during and after the plunge because the remnant is less compact.Thus, the total thermal pressure is reduced for lower effective masses, even though the thermal index is higher.
• The dynamical ejecta mass is correlated with the effective mass and the incompressibility.On the one hand, larger incompressibilities and lower effective masses reduce the compactness of the initial NS and therefore increase the amount of tidal ejecta.This effect is especially strong for higher incompressibilities.The mass of the shock-heated ejecta, on the other hand, are mostly correlated by the effective mass.Mass ejection in the early disk phase is more complex because multiple effects, such as remnant deformation, oscillations, and neutrino emission play a role.
• We matched all nuclear matter properties to those of Shen EOS.The resulting SkShen EOS is mostly similar to the Shen EOS except for larger pressures at high densities and a slightly lower thermal index.In the corresponding simulations, the evolution of the remnant and the accretion disk, as well as the amounts of ejected mass are relatively similar, and the post-merger GW spectrum of the two models is remarkably similar.Therefore, we conclude that nuclear matter properties are a useful measure to quantify EOS effects in BNS mergers.
Nuclear matter properties are a promising way to quantify EOS effects in BNS mergers that is independent of the microphysical frame-work used.Furthermore, they can be constrained by nuclear physics experiments and theoretical calculations, enabling the combination of results from nuclear physics, astrophysics, and multi-messenger astronomy.However, more studies investigating the role of EOS properties in compact object mergers are needed in the future.

Figure 1 .
Figure 1.The mass-radius relation of cold non-rotating NSs for the EOSs used in this paper.The horizontal dashed line marks 1.365 ⊙ which is the mass of the NSs explored in the simulations in this work.

Figure 2 .
Figure 2. Pressure versus density at  = 0 and  e = 0.1 for the EOSs used in this work.The lower panel shows the curves normalized to the pressure in the fiducial model LS220 † for comparison.The vertical dashed lines mark the saturation density 0.155 fm −3 .
= 0 ,=0,  =0.1 at saturation density for  = 0 and   = 0.1.Both quantities are shown in the top panel of Fig. 3 and given in Table 2 for all EOS models considered in this work.As discussed P| n = n0, T = 0, Ye = 0.1 [MeV fm 3 ]

Figure 3 .
Figure 3. Top panel: the pressure and the dimensionless slope of the pressure at saturation density for  = 0 and   = 0.1.Arrows show the relevant nuclear matter property changes.Lower panel: the pressure at saturation density for  = 0 and   = 0.1 and the radius of the initial NSs.

Figure 4 .
Figure 4. Thermal index Γ th for the LS220 † ,  * 0.8 ,  * S , SkShen, and Shen EOSs.SkShen exhibits a slightly different Γ th compared to  * S due to the different saturation density.Although the saturation density and the effective mass at saturation density are the same in SkShen and Shen, their thermal index is different due to the underlying relativistic mean field functional in Shen.

Figure 5 .
Figure 5. Maximum density (upper panel) and average remnant temperature (lower panel) versus time.

Figure 6 .
Figure 6.Complex azimuthal  = 1 and 2 decomposition of the density in the  -plane calculated with Eq. (4).Solid and dashed lines correspond to the  = 2 and  = 1 modes, respectively.The colors denote the different models and have the same meaning as in Fig.5.The transition to a dominant  = 1 mode happens first for LS175 † , followed by ( *   sym ) S , and ( *  ) S .

Figure 9 .
Figure 9. Evolution of the disk mass for all simulations except LS175 † .The top and bottom panels correspond to different definitions of the disk boundary.
Fig. 11, we examine the detailed ejecta composition for the simulation  * S .We use  * S as reference here but all other simulations show similar trends to the ones discussed below.The left panel shows the 2D histogram of the electron fraction of the ejected fluid elements versus the time at which they reach 300 km.The right panel shows the histogram of their electron fraction versus the polar angle of ejection .The division into dynamical and disk ejecta components is approximately shown by the dashed contour in the left panel.The dashed horizontal line shows the cut at  e = 0.1.Tidally ejected

Figure 10 .
Figure 10.Mass-weighted histogram of the electron fraction, entropy per baryon, and temperature in the disk at  = 30 ms.Dotted lines show the histograms when matter with  > 10 12 g cm −3 is excluded.The EOS models are split into two panels for better readability.

Figure 11 .
Figure 11.Mass-weighted histograms of ejecta properties in model  * S .Top panel: Two-dimensional histogram of the electron fraction versus the time after the merger at which the ejecta cross the detection sphere at 300 km.The dashed contour marks the transition point from dynamical ejecta to shockheated ejecta and the horizontal dashed line shows the cut at  e = 0.1 used to distinguish the tidal ejecta.Bottom panel: Two-dimensional histogram of the electron fraction versus the polar angle.

Figure 12 .
Figure12.Mass-weighted histogram of the electron fraction, entropy per baryon, and asymptotic velocity of the ejected matter.Solid and dotted lines represent dynamical ejecta and disk ejecta, respectively.The EOS models are split into two panels for better readability.

Figure 13 .
Figure 13.Mass ejection rate and total ejected mass of the disk ejecta component as a function of time.
Figure14.Masses of the tidal and shock-heated ejecta components as a function of the dimensionless slope of the pressure and the pressure at saturation density, respectively.

Figure 15 .Figure 16 .
Figure 15.Distribution of  e in the   plane for the models LS220 † ,  * S , LS175 † , and ( *  ) S 15 ms after the merger.The tidal ejecta are visible as large, smooth areas with low  e .The dimensionless slope of the pressure and accordingly the size of the tidal arms increases from the left to the right panels.