Models of binary neutron star remnants with tabulated equations of state

ABSTRACT

The emergence of novel differential rotation laws that can reproduce the rotational profile of binary neutron star merger remnants has opened the way for the construction of equilibrium models with properties that resemble those of remnants in numerical simulations. We construct models of merger remnants, using a recently introduced 4-parameter differential rotation law and three tabulated, zero-temperature equations of state. The models have angular momenta that are determined by empirical relations, constructed through numerical simulations. After a systematic exploration of the parameter space of merger remnant equilibrium sequences, which includes the determination of turning points along constant angular momentum sequences, we find that a particular rotation law can reproduce the threshold mass to prompt collapse to a black hole with a relative difference of only |$\sim 1{{\ \rm per\ cent}}$| with respect to numerical simulations, in all cases considered. Furthermore, our results indicate a possible correlation between the compactness of equilibrium models of remnants at the threshold mass and the compactness of maximum-mass non-rotating models. Another key prediction of binary neutron star merger simulations is a relatively slowly rotating inner region, where the angular velocity Ω (as measured by an observer at infinity) is mostly due to the frame dragging angular velocity ω. In our investigation of the parameter space of the adopted differential rotation law, we naturally find quasi-spherical (Type A) remnant models with this property. Our investigation clarifies the impact of the differential rotation law and of the equation of state on key properties of binary neutron star remnants and lays the groundwork for including thermal effects in future studies.

1 INTRODUCTION

The detection of gravitational waves (GW) from the inspiral phase of the GW170817 binary neutron star (BNS) merger event (Abbott et al. 2017a, b) combined with complementary information from its electromagnetic counterpart (Abbott et al. 2017b, c; Goldstein et al. 2017) have produced new constraints on the equation of state (EOS), see Bauswein et al. (2017), Abbott et al. (2018), De et al. (2018), Fattoyev, Piekarewicz & Horowitz (2018), Most et al. (2018), Abbott et al. (2019), Montaña et al. (2019), Capano et al. (2020), Landry, Essick & Chatziioannou (2020), Dietrich et al. (2020), Breschi et al. (2021), and references therein, as well as Chatziioannou (2020), Dietrich, Hinderer & Samajdar (2021) for recent reviews. A second likely BNS merger event, GW190425, was reported in Abbott et al. (2020b) and more are expected in the next years (Abbott et al. 2020a). Although the sensitivity of the LIGO and Virgo GW detectors was not sufficient to detect the post-merger phase in GW170817 (Abbott et al. 2017a, d), such a detection is likely to be achieved in the future, either with upgraded or next-generation detectors, see e.g. Clark et al. (2014, 2016), Chatziioannou et al. (2017), Bose et al. (2018), Yang et al. (2018), Torres-Rivas et al. (2019), Martynov et al. (2019), Oliver et al. (2019), Easter et al. (2019), Tsang, Dietrich & Van Den Broeck (2019), Breschi et al. (2019), Hall & Evans (2019), Easter et al. (2020), Ackley et al. (2020), Haster et al. (2020), Aggarwal et al. (2021), Ganapathy et al. (2021), Page et al. (2021).

The outcome of a BNS merger is closely tied to the EOS and the total mass M = m₁+m₂ of the system, where m₁ and m₂ are the binary’s components masses, see Shibata & Hotokezaka (2019), Bernuzzi (2020), Radice, Bernuzzi & Perego (2020), Friedman & Stergioulas (2020) for recent reviews. If M < M_thres (the threshold mass for prompt black hole formation), the merger results in a hot, massive, and differentially rotating, compact object with a substantial material disc around it. If, at the same time, M > M_{max, rot} (the maximum mass of a cold, uniformly rotating neutron star), the remnant will initially survive several tens of milliseconds (ms) due to the support of differential rotation and thermal pressure. However, the loss of angular momentum, due to GW emission, as well as dissipative effects (e.g. shear viscosity, magnetic breaking, and effective viscosity due to the development of the magnetorotational instability, see Shibata & Hotokezaka (2019), Ciolfi (2020), Sarin & Lasky (2021), Ruiz, Shapiro & Tsokaros (2021) for recent reviews and also Radice (2020)) will ultimately lead to a delayed collapse to a black hole. A remnant with mass M_{max, rot} > M > M_max, where M_max is the maximum mass of a non-rotating star, will be long-lived, spinning down on the time-scale of electromagnetic emission, before reaching the axisymmetric instability limit. Only if M < M_max, can a stable remnant form.

During the first few milliseconds after its formation, the remnant is still highly non-axisymmetric, featuring also strong non-linear oscillations and deformations away from equilibrium. Characteristic non-linear features are combination tones and spiral deformations (Stergioulas et al. 2011; Bauswein & Stergioulas 2015; Bauswein, Stergioulas & Janka 2016; Bauswein & Stergioulas 2019). On a somewhat longer time-scale, one can regard the remnant as a quasi-stationary, slowly drifting equilibrium state with the addition of linear oscillations. If one neglects some aspects of the state of the remnant (non-axisymmetric deformations, oscillations, time-dependence, and thermal structure), one can construct simplified, stationary axisymmetric models of its structure.

Merger remnants that survive for more than a few milliseconds before collapsing to black holes have been studied through numerical simulations (Hotokezaka et al. 2011; Sekiguchi et al. 2011; Bauswein et al. 2012; Bauswein & Janka 2012; Hotokezaka et al. 2013; Bernuzzi et al. 2014; Dietrich et al. 2015; De Pietri et al. 2016; Radice et al. 2018). The remnant’s structure, including its rotation profile was studied extensively in Kastaun & Galeazzi (2015), Paschalidis et al. (2015), Bauswein & Stergioulas (2015), Kastaun, Ciolfi & Giacomazzo (2016), East et al. (2016), Endrizzi et al. (2016), Kastaun et al. (2017), Ciolfi et al. (2017), Hanauske et al. (2017), Endrizzi et al. (2018), Kiuchi et al. (2018), Ciolfi et al. (2019), East et al. (2019), De Pietri et al. (2020), Kastaun & Ohme (2021). A common finding is that the remnant’s rotation profile exhibits a maximum away from the centre, which is in sharp disagreement with the differential rotation law by Komatsu, Eriguchi & Hachisu (1989), hereafter KEH, which was widely used in the context of differentially rotating neutron stars (see Ansorg, Gondek-Rosińska & Villain 2009; Espino & Paschalidis 2019; Espino et al. 2019 for different types of equilibrium models that can be constructed with the KEH rotation law).

A 3-parameter piecewise extension of the KEH rotation law was used in Bauswein & Stergioulas (2017), Bozzola, Stergioulas & Bauswein (2018), in order to allow the outer regions to rotate more slowly than the core, reaching high masses (typical of remnants) without encountering mass-shedding (see also Galeazzi, Yoshida & Eriguchi (2012), Uryū et al. (2016) for other rotation laws). Two different 3-parameter and 4-parameter rotation laws were proposed by Uryū et al. (2017), who presented selected example equilibrium models.

Zhou et al. (2019) constructed differentially rotating strange star models, using the 4-parameter rotation law of Uryū et al., whereas Passamonti & Andersson (2020) and Xie et al. (2020) studied the onset of the low T/|W| instability (Watts, Andersson & Jones 2005) in models constructed with the 3-parameter rotation law of Uryū et al. In Camelio et al. (2021) models of stationary remnants of a BNS merger at ∼10−50 ms after merger were presented, which were differentially rotating, hot, and baroclinic, using their own, 5-parameter rotation law. The models were constructed with the assumption of spatial conformal flatness (IWM-CFC approximation).

An important aspect of modelling post-merger remnants is to separate the effects of (i) the differential rotation law, (ii) the cold part of the EOS, and (iii) thermal effects on the structure of the remnant and on its dynamical properties (stability to axisymmetric collapse and oscillations). To do so, we embarked on a systematic study of each of these three effects in separation from the other two. Our first step was to present equilibrium sequences of rotating relativistic stars, constructed with the 4-parameter rotation law of Uryū et al., adopting a cold, relativistic N = 1 polytropic EOS and choosing rotational parameters motivated by simulations of binary neutron star merger remnants (Iosif & Stergioulas 2021). A distinctive feature of the Uryū et al. law is that it allows for the angular velocity to attain a maximum value Ω_max away from centre (as seen in simulations), which was not possible with the KEH law. We compared the sequences of equilibrium models to published sequences that used the KEH rotation law, revealing only a small influence of the choice of rotation law on the mass of the equilibrium models and a somewhat larger influence on their radius. Both Type A and Type C solutions (in correspondence to the classification of KEH-type models by Ansorg et al. 2009) were found. While our models were highly accurate solutions of the fully general relativistic structure equations, we also demonstrated that for models relevant to merger remnants the IWM-CFC approximation still maintains an acceptable accuracy.

Here, we take a second step in this program and construct sequences of models of post-merger remnants, using the 4-parameter rotation law of Uryū et al. and different tabulated EOS. In Bauswein & Stergioulas (2017), the threshold mass to black hole collapse, as determined in simulations, was reproduced in a semi-analytic way, using equilibrium models obeying a piecewise extension of the KEH rotation law. Following the analysis detailed in that work, our sequences are constructed using an empirical relation between angular momentum at merger and the radius and compactness of the progenitor stars (assuming equal masses). Again, we find both Type A and Type C sequences. For a particular combination of rotation-law parameters, we find that the sequence of merger remnants terminates very near the threshold mass to collapse (as obtained by numerical simulations) for all three representative EOS that we used in this study. For somewhat different combinations of rotation-law parameters, we find sequences of merger remnants with realistic rotation profiles, for which the angular velocity in the core is close to the angular velocity of frame dragging, reproducing a characteristic feature seen in binary neutron star merger simulations. The next step in this program will be the inclusion of thermal effects, which we are planning to present in the future.

The structure of the paper is as follows: in Section 2 we discuss the theoretical framework and numerical methods. In Section 3, we present the main results. In Section 4, we discuss our findings.

Throughout the text we set c = G = 1 in equations (except for equations where units are explicitly included) and choose appropriate physical units to report numerical results. We also denote with R_X the radius of non-rotating neutron stars with gravitational mass |$X \, \mathrm{M}_\odot$|⁠. E.g. R_1.4 stands for the radius of a |$1.4\, \mathrm{M}_\odot$| star.

2 THEORETICAL FRAMEWORK AND METHODS

2.1 Space-time metric and matter description

2.2 Rotation law

2.3 Numerical scheme

In order to build our equilibrium configurations, we use an extended version of the rns code (Stergioulas & Friedman 1995), which implements the iterative Komatsu et al. (1989) scheme with improvements by Cook, Shapiro & Teukolsky (1992). The initial rns code was updated to tackle differential rotation in Stergioulas, Apostolatos & Font (2004) and further extended for the 3-parameter, piecewise KEH law in Bauswein & Stergioulas (2017). In Iosif & Stergioulas (2021), we extended rns with the implementation of the 4-parameter Uryu + rotation law of equation (4). This allowed for the construction of models with realistic rotation profiles for BNS merger remnants that have off-centre maxima in the angular velocity profile. The solutions were shown to be highly accurate and converging at second order with an increasing number of grid points. A standard resolution was chosen that yields solutions with 3D virial theorem index (GRV3) of the order of 10⁻⁵. In this study, we employ a grid size of SDIV × MDIV = 401 × 201 (compactified radial times angular) for all models. We refer to Iosif & Stergioulas (2021) for further details on the numerical scheme.

2.4 Equations of state

Considerable uncertainty still exists in the determination of the EOS of dense nuclear matter. Fig. 1 shows the gravitational mass M versus the circumferential radius R for non-rotating models constructed with several different hadronic EOS that cover the large uncertainty range that existed before the historic detection of gravitational waves from the source GW170817. The initial analysis of GW170817 resulted in a constraint on neutron star radii |$R=11.9_{-1.4}^{+1.4}$| km (Abbott et al. 2018) for both stars involved in the merger, at the 90 per cent credible level. In the meantime, a large number of studies presented multimessenger constraints on the neutron star radius, taking into account observations in the electromagnetic spectrum as well as nuclear-theory computations using chiral effective field theory. Recent studies provide predictions for the radius of a |$1.4 \, \mathrm{M}_{\odot }$| neutron star with an uncertainty range of |$R_{1.4}=12.32_{-1.47}^{+1.09}$| km (90 per cent credible level) (Landry et al. 2020), |$R_{1.4}=11.0_{-0.6}^{+0.9}$| km (90 per cent credible level) (Capano et al. 2020), |$R_{1.4}=11.75_{-0.81}^{+0.86}$| km (90 per cent credible level) (Dietrich et al. 2020), |$R_{1.4}=12.2_{-0.5}^{+0.5}$| km (1σ level) (Breschi et al. 2021), and |$R_{1.4}=11.94_{-0.87}^{+0.76}$| km (90 per cent credible level) (Pang et al. 2021). Furthermore, the Neutron Star Interior Composition Explorer (NICER; Gendreau et al. 2016) measurements of PSR J0740+6620 have yielded radius estimates of |$R_{1.4}=12.45_{-0.65}^{+0.65}$| km (1σ level) (Miller et al. 2021), km (95 per cent credible level) and |$R_{1.4}=12.18_{-0.79}^{+0.56}$| km (95 per cent credible level) for different high-density EOS parametrizations (Raaijmakers et al. 2021; Riley et al. 2021).

Taking into account the above range of radii uncertainties, we selected three tabulated, zero-temperature, hadronic EOS that correspond to typical neutron star radii between 11 and 13 km. These are APR (Baym, Pethick & Sutherland 1971; Akmal, Pandharipande & Ravenhall 1998; Douchin & Haensel 2001), DD2 (Möller, Nix & Kratz 1997; Hempel & Schaffner-Bielich 2010; Typel et al. 2010) and MPA1 (Müther, Prakash & Ainsworth 1987), shown with darker colours in Fig. 1. All three EOS satisfy the current constraints for the maximum neutron star mass (Demorest et al. 2010; Antoniadis et al. 2013; Cromartie et al. 2020) as well as the minimum radius constraint, when combining causality and GW170817, of R_1.6 ≥ 10.68 km (Bauswein et al. 2017). EOS with strong phase transitions are not included in this study.

2.5 Construction of merger remnant sequences

Our aim is to construct sequences of equilibrium models that mimic characteristic properties of post-merger remnants and reach the threshold mass to prompt collapse.

We construct sequences of models of merger remnants with different combinations of {λ₁, λ₂} and with remnant masses of |$M_\text{tot} = \lbrace 2.2, 2.3, 2.4, 2.5, \dots \rbrace \, \mathrm{M}_\odot$|⁠. We continue to larger values with a step of |$0.1 \, \mathrm{M}_\odot$| up to the maximum possible M_tot for which we can construct an equilibrium sequence for the particular rotation law and EOS. For each value of M_tot, we compute the corresponding R_NS and M_NS of a non-rotating star. From (8) we compute the corresponding angular momentum of the remnant. The detailed properties of the equilibrium models of merger remnant sequences are reported and discussed in detail in Section 3.

2.6 Constant angular momentum sequences and turning points

The dynamical instability to prompt collapse to a black hole is detected through numerical simulations or by finding the models for which the frequency of the fundamental quasi-radial mode vanishes. For uniform rotation, the dynamical instability limit for prompt collapse is very close to the secular instability limit and it occurs slightly earlier (see Friedman & Stergioulas 2013) for a detailed discussion). In the case of differential rotation, Weih, Most & Rezzolla (2018) demonstrated (through numerical simulations) that for particular choices of the KEH law the dynamical instability also sets in very close to the secular instability limit (the central density of dynamically unstable models was at most several per cent smaller than the central density at the turning points).

Given the above findings for uniformly rotating as well as differentially rotating models with the KEH law and since we do not yet have dynamical or perturbative calculations for models constructed with the Uryu + law, we adopt the line connecting the turning points of constant angular momentum sequences as a reasonably approximate indicator of dynamical instability.

3 MAIN RESULTS

3.1 Sequences of type C models and threshold mass to prompt collapse

Initially, we focus on two combinations of {λ₁, λ₂} that were shown (Iosif & Stergioulas 2021) to yield Type C solutions according to the classification of Ansorg et al. (2009). These are sequences along which the models transition continuously from quasi-spherical to quasi-toroidal configurations, as the polar to equatorial axial ratio r_p/r_e is decreased. We find that setting parameters {λ₁, λ₂} equal to the pairs of values {2.0, 0.5} and {1.5, 0.5}, continues to result in Type C solutions for the tabulated EOS we examined, as was the case for the N = 1 polytropes in Iosif & Stergioulas (2021).

As can be seen in Fig. 2 for the APR EOS, the J-constant curves for both of these rotation laws are overlapping and the turning points we locate are quite close as well. This is to be expected, as these two particular rotation laws correspond to similar Ω(r) rotational profiles (fig. 8 Iosif & Stergioulas 2021). Asterisks in black and red denote the remnant models found for the respective J_merger values predicted by the empirical relation (8). The gravitational masses of these configurations for the APR EOS start at |$2.2 \, \mathrm{M}_\odot$| for the least massive model and end at |$2.84 \, \mathrm{M}_\odot$| for the most massive model. The picture is similar for the other two EOS we consider, DD2 (see Fig. 3) and MPA1 (see Fig. 4). The only difference is that higher masses (as well as higher angular momentum values) are reached for the most massive remnant models at |$3.28 \, \mathrm{M}_\odot$| for DD2 and at |$3.2 \, \mathrm{M}_\odot$| for MPA1.

Qualitatively, the remnant sequence rises to larger masses at an almost constant (steep) slope that after a point abruptly diminishes, allowing the remnant sequence to intersect with the turning point line. As in Bauswein & Stergioulas (2017), we find that this intersection is related to the threshold mass for prompt collapse, |$M_\text{thres}^\text{sim}$| (as determined by numerical simulations in Bauswein et al. 2020). This mass is indicated by a blue horizontal dash-dotted line in Figs 2, 3, and 4. The intersection can also be determined in Figs 5, 6, and 7, which show the angular momentum J as a function of the gravitational mass M of the remnant sequence (blue line) and the line connecting the turning points of J-constant sequences for {λ₁, λ₂} = {2.0, 0.5} (black line) and {λ₁, λ₂} = {1.5, 0.5} (red line), for the three EOS. Note that the data for the remnant sequences in Figs 5, 6, and 7, as obtained from the empirical relation (8), form a straight line, in agreement with the form of the original empirical relation (7).

Also, the agreement at the level of 1 per cent between our current models and numerical simulations is a significant improvement with respect to the agreement at the level of 3 per cent – 7 per cent in Bauswein & Stergioulas (2017), where a particular 3-parameter piecewise KEH-type rotation law was used. This is achieved without a direct reconstruction of particular merger remnants (i.e. without trying to match all properties of a remnant, as extracted from a numerical simulation), but using equilibrium models for particular {λ₁, λ₂} values, in conjunction with the empirical relation for the angular momentum at merger.

The above findings indicate that the Uryu + law is not simply qualitatively appropriate for merger remnants, in the sense that it allows for the maximum angular velocity to appear off-centre, but it can also yield precise numerical results, at least for certain properties of the remnants.

3.2 Domain of type A solutions

Merger remnants that do not collapse promptly, can evolve towards nearly axisymmetric, quasi-stationary configurations (at least before a possible delayed collapse sets in) that can be approximated with suitable equilibrium models. This involves Type A solutions,² i.e. sequences of models that remain quasi-spherical (the maximum density is always at the centre) as the axial ratio r_p/r_e is decreased (i.e. the rotation rate increases) until the mass-shedding limit is reached.

In our recent work Iosif & Stergioulas (2021), we found that the Uryu + rotation law with {λ₁, λ₂} = {2.0, 1.0} and with {λ₁, λ₂} = {1.5, 1.0} yields Type A solutions for the N = 1 polytropic EOS. In addition, we highlighted the fact that according to recent numerical simulations (Hanauske et al. 2017; De Pietri et al. 2020), a value of λ₂ = 1 seems to be favoured over λ₂ = 0.5, for the case of compact remnants from BNS mergers, while λ₁ ranges between 1.7 and 1.9 for realistic EOS. We therefore probe in more detail models with this range of parameters. To that end, we set the value of λ₂ to 1.0 and explore the range λ₁ ∈ [1.5, 2.0] with a step of 0.1.

A vertical ‘scan’ of the mass versus ϵ_max parameter space for specific Type A {λ₁, λ₂} pairs (fixing the value of the maximum energy density and gradually decreasing the axial ratio r_p/r_e), revealed that these sequences reached a point, after which it was not possible to further construct equilibrium solutions as the maximum density was increased. This behaviour is much more stark for the case of {λ₁, λ₂} = {2.0, 1.0} than for {λ₁, λ₂} = {1.5, 1.0}. We note that the terminal models encountered for each choice of {λ₁, λ₂} pairs are not close to mass-shedding. For the highest {λ₁, λ₂} = {2.0, 1.0}, the Ω_e/Ω_K ratio of the terminal models ranges between 0.6 and 0.7 as the angular momentum increases, whereas for {λ₁, λ₂} = {1.5, 1.0} it ranges between 0.6 and 0.8. However, the maximum density remained at the centre of the configuration even at the highest rotation rates that were achieved and we classify these models as Type A.

Fig. 8 displays the six Type A remnant sequences corresponding to the six {λ₁, λ₂} pairs that we investigated³ for the APR EOS. Remnant models are shown as asterisks (different colours correspond to different {λ₁, λ₂} values). The terminal models along each J-constant sequence for each {λ₁, λ₂} pair are shown as dots of matching colour. We find that the choice of λ₁ = 1.6 allows for the most massive Type A remnant model for this EOS, with a gravitational mass of |$M_\text{tot}=2.7 \, \mathrm{M}_\odot$|⁠.

For the {λ₁, λ₂} = {1.5, 1.0} and {λ₁, λ₂} = {1.6, 1.0} cases, we explicitly show the constant angular momentum sequences (as light blue and dark blue dashed lines, respectively). We note that the {λ₁, λ₂} = {1.6, 1.0} J-constant lines in Fig. 8 nearly merge into those of {λ₁, λ₂} = {1.5, 1.0} as the maximum density increases. J-constant lines constructed with other pairs of {λ₁, λ₂} also tend to merge with those of {λ₁, λ₂} = {1.5, 1.0}, but for clarity, we omit the J-constant lines for λ₁ = {1.7, 1.8, 1.9, 2.0} in Fig. 8. Moreover, as a rule of thumb the J-constant lines for {λ₁, λ₂} = {1.5, 1.0} reach higher maximum densities than the J-constant lines for {λ₁, λ₂} = {1.6, 1.0}. Note though, that for the highest J-constant sequence for the APR and DD2 EOS (Figs 8 and 9 respectively), this trend is reversed, i.e. it is the highest J-constant line for {λ₁, λ₂} = {1.6, 1.0} that actually reaches higher maximum densities. In our investigation of this behaviour, we found that an increasingly smaller step size in axial ratio was required (reaching as low as 10⁻⁴) to locate equilibrium solutions for {λ₁, λ₂} = {1.5, 1.0} (which yields the weakest differential rotation we consider) in the parameter space of increasing energy density values and high angular momentum values. This effect is indicating numerical stiffness. However, one would have to perform similar trials with other equilibrium codes to clarify whether this issue: (i) has some physical origin (e.g. it could be possible that weaker differential rotation cannot produce high enough angular velocities to accommodate arbitrarily high angular momentum values), (ii) or if it is due to numerical stiffness of the problem, (iii) or if it is due to the existence of different types of differentially rotating solutions that have the same central density and axial ratio.⁴

Gathering all the evidence, some interesting observations can be made in connection to earlier works in the literature, where the KEH rotation law was used. First of all, we note that Type C remnant models are able to reach higher masses than Type A models, in agreement with findings in Gondek-Rosińska et al. (2017), Studzińska et al. (2016), Espino & Paschalidis (2019) for the KEH rotation law. Concerning the terminal models encountered for the Type A J-constant sequences, they can be interpreted as a confirmation that the domain of Type A solutions shrinks for higher densities and stronger differential rotation. Specifically, with stronger differential rotation we do not find Type A solutions above a certain maximum energy density, whereas we can still find Type C solutions⁵ (or Type A solutions with a weak differential rotation).

The above property of differentially rotating models was highlighted for the KEH rotation law in Studzińska et al. (2016), Gondek-Rosińska et al. (2017) for polytropes, Espino & Paschalidis (2019) for realistic EOS and Szkudlarek et al. (2019) for strange quark stars. From our findings for the Uryu + rotation law, it seems that the different types of solutions are not tied to the particular KEH rotation law (for which they were originally discovered), but appear also for other, more general rotation laws, such as the one by Uryū et al. considered here.

Figs 9 and 10 show the same investigation of Type A models as in Fig. 8, but for the EOS DD2 and MPA1. We note that also for these EOS, the choice of λ₁ = 1.6 allows for the construction of the most massive Type A remnant model and also for a full remnant sequence (i.e. starting from the lowest remnant mass of |$M_\text{tot}=2.2 \, \mathrm{M}_\odot$| we consider). Within the parameter range that we investigated, the highest mass Type A remnant model reached with EOS DD2 was |$3.1 \, \mathrm{M}_\odot$|⁠, while for EOS MPA1 it was |$3 \, \mathrm{M}_\odot$|⁠.

In Iosif & Stergioulas (2021, fig. 8) we showed that for {λ₁, λ₂} = {1.5, 1.0} the differential rotation is quite weak, compared to other values of λ₁ in the range we consider here. For the stiffest EOS we consider (DD2) this leads to Type A models for which we can find turning points along the J-constant sequences (light blue dotted curve in Fig. 9). Moreover, these turning points are close to the axisymmetric instability limit for uniform rotation.

3.3 Frame dragging contribution to rotation

Figs 11, 12, and 13 show the angular velocity Ω(r_c) rotational profiles versus the circumferential radial coordinate r_c at the equatorial plane, for the three EOS under study, APR, DD2, and MPA1. For each EOS, a triplet of panels is shown, with each left-hand panel corresponding to rotation laws {λ₁, λ₂} = {2.0, 0.5}, each middle panel to {λ₁, λ₂} = {1.5, 0.5}, and each right-hand panel to {λ₁, λ₂} = {1.6, 1.0}. Every individual panel shows Ω(r_c), as well as the frame dragging angular velocity ω(r_c) for the most massive, the least massive and for an intermediate mass remnant model constructed with the particular choice of parameters {λ₁, λ₂}. A common finding shared among the three EOS explored is that models with {λ₁, λ₂} = {2.0, 0.5} reach the highest angular velocity compared to the other two options and models with {λ₁, λ₂} = {1.6, 1.0} reach the smallest angular velocity peaks in their profile. This is in agreement with corresponding rotational profiles for these parameter values for polytropic equilibrium models with N = 1 (fig. 8 Iosif & Stergioulas 2021).

Examining Type A remnants, we observe that for the most massive models the central part of the configuration (i.e. approximately up to |$r_c \sim 5 \, \text{km}$|⁠) rotates slowly compared to the rest of the configuration. However, this rotation rate is mostly due to the contribution of the frame dragging ω, which means that with respect to a ZAMO this part of the model is almost non-rotating. To our knowledge, this is the first time that a differential rotation law has been shown to reproduce this feature, already known from numerical simulations: similar rotation profiles have been presented and analysed in Kastaun & Galeazzi (2015), Endrizzi et al. (2016), Kastaun et al. (2016, 2017), Ciolfi et al. (2017).

Note that we did not observe the Ω ∼ ω behaviour near the centre of the star only for the most massive Type A models (i.e. the rightmost dark blue asterisks of the remnant sequences for {λ₁, λ₂} = {1.6, 1.0} in Figs 8, 9, and 10 corresponding to the light blue curves of the right-hand panels in Figs 11, 12, and 13 for each EOS considered). This feature also appears in less massive models that have higher central densities (i.e. models neighbouring the terminal models shown as dark blue dots for {λ₁, λ₂} = {1.6, 1.0} in Figs 8, 9, and 10). Specifically, for all the rotation law parameters yielding Type A configurations that we considered (i.e. the six {λ₁, λ₂} pairs with λ₂ = 1.0 in Figs 8, 9, and 10), we found that ω approaches Ω as the maximum density increases along a J-constant sequence, until the Ω ∼ ω feature appears as we reach the terminal model. To clarify even further, we found that

• along a J-constant sequence, the Ω ∼ ω behaviour appears at lower central densities as the strength of differential rotation is increased⁶ (see left-hand panel of Fig. 14).

• for a specified strength of differential rotation (i.e. fixing the differential rotation law’s parameters at certain values), the Ω ∼ ω behaviour appears at lower central densities as the angular momentum J is increased (see middle panel of Fig. 14).

• as the maximum energy density is increased along a J-constant sequence, the frame dragging at the centre ω_c increases faster than the central angular velocity Ω_c, thus giving rise to the Ω ∼ ω feature (see right-hand panel of Fig. 14).

Note that while Fig. 14 demonstrates the above for the MPA1 EOS, similar behaviour is observed for the other two EOS that we studied.

3.4 Ω(F) profiles

Figs 15, 16, and 17 follow the same organization as corresponding figures for the rotation profiles Ω(r_c) and present the angular velocity Ω(F) profiles that define each rotation law. We find that for each choice of rotation law parameters (e.g. if one observes the left-hand panels of Figs 15, 16, and 17 ‘vertically’, etc.) the resulting Ω(F) profiles are qualitatively similar between the three EOS that we consider. This is in agreement with analogous behaviour of the corresponding Ω(r_c) profiles in Figs 11, 12, and 13.

Note that in all cases, the inverse profile, F(Ω), would not be a one-to-one function. In the case of the KEH rotation law, one can simply integrate F(Ω) in the equation of hydrostationary equilibrium. However, for the Uryu + rotation law, one needs to express the equation of stationary equilibrium in terms of an integral Ω(F) (see Iosif & Stergioulas (2021) for details).

3.5 Structure of the remnants: surface and density distribution

Stellar surfaces in the x − z plane (with x and z normalized with the equatorial coordinate radius r_e) for the most massive, least massive, and an intermediate mass remnant model are shown in Figs 18, 19, and 20, for the three EOS employed and for three different choices of parameters {λ₁, λ₂} (similar to the corresponding figures of the rotation profiles). We note that for tabulated EOS we always have a finite minimum pressure value P_min at the minimum density ρ_min. For the three EOS employed, the range of ρ_min is sub-grid varying from |$\sim 8 \, \mathrm{g \, cm^{-3}}$| to |$\sim 10^3 \, \mathrm{g \, cm^{-3}}$|⁠, i.e. so small that not a single grid point is removed from the model. Furthermore, since the rns code constructs numerical tables in log enthalpy, a very small enthalpy value H_min (of the order of 10⁻¹⁶) is used in the tabulated EOS files. Then, the stellar surface is defined as the location where the pressure P = P_min and the enthalpy H = H_min. Interpolating the model’s enthalpy by demanding that its value matches H_min allows us to pinpoint the surface’s location. In addition, Figs 21, 22, and 23 display the 2D rest-mass density distributions for the most massive models in the meridional plane.

Note that the surfaces and meridional density profiles are qualitatively similar for the Type C solutions obtained with {λ₁, λ₂} = {2.0, 0.5} and {λ₁, λ₂} = {1.5, 0.5}. Both choices lead to a quasi-toroidal surface shape, typical for Type C solutions. For {λ₁, λ₂} = {2.0, 0.5}, we observe a stronger deformation close to the rotation axis than for {λ₁, λ₂} = {1.5, 0.5}, which is explained by the stronger differential rotation.

For the Type A solutions obtained with {λ₁, λ₂} = {1.6, 1.0}, the surfaces of all models (most massive to least massive) are quite similar to each other (when coordinates are scaled by r_e). Even though the remnant models have significant oblateness, they still retain their quasi-spherical shape (the central density is also the maximum density). It is interesting to note that for these selected Type A remnants, the axial ratio r_p/r_e ranges between ∼0.47 and 0.54, whereas for the Type C remnants the corresponding range is considerably lower at ∼0.3−0.42.

Recently, Kastaun et al. (2016) introduced a new measure that replaces density profiles, mass, and compactness in a way that can be used unambiguously for rapidly and differentially rotating merger remnants, without a clearly defined surface. Therefore, we stress that the profiles presented here serve simply as an indication about the different configurations possible for the different cold EOS employed and for the {λ₁, λ₂} options we considered. In a realistic remnant with a hot envelope and mass-shedding, the density distribution would not terminate at the same radius as in our models and one would need to define an approximate surface shape, based e.g. on the location where the density drops to a certain fraction of the maximum density.

3.6 Mass versus equatorial radius

Having presented the main characteristics of our models in the previous subsections, we further analyse our findings by constructing M(R_e) plots (i.e. gravitational mass versus equatorial circumferential radius) for the remnant models. We focus on the choice {λ₁, λ₂} = {2.0, 0.5} that represents the strongest differential rotation we consider, but note that a similar picture holds for the case {λ₁, λ₂} = {1.5, 0.5}.

Figs 24, 25, and 26 show M(R_e) for the EOS APR, DD2, and MPA1 respectively. In each of these figures, remnant sequences are shown as dots (denoting the equilibrium models), connected by solid interpolated lines. The dotted lines represent J-constant sequences, turning points are marked by crosses and a linear fit (dashed line) approximates the turning point sequence. The annotation in these figures means that while the empirical relation (8) by Lucca et al. (2021) provides a predicted J_merger value for a desired target value M_tot, the intersection of the remnant sequence and the turning point sequence has already taken place. Therefore, the target M_tot model for the specific J_merger value predicted by (8), does not exist, since it would exceed the value of M_thres.

Collecting data from all EOS in a single M(R_e) plot, Fig. 27, one could extrapolate the remnant sequences so that they intersect with the turning point sequences and obtain the intersection (M_thres, R_{e − thres}). We note that the M_thres values determined in this way are in good agreement with those reported in Table 2. Still, we regard the values in Table 2 as better estimates for |$M_\text{thres}^\text{eq}$|⁠, since they involve bulk quantities of the star, such as the angular momentum J and the mass M, with direct input from numerical simulations, via (8), to determine J. Moreover, the precise determination of R_{e − thres}, in actual remnants is affected by the thermal properties of low-density material and one needs to resort to a particular definition, based on an isodensity surface, where the density has become a certain small fraction of the maximum density. The notion of a ‘bulk’ region and the relative measure introduced in Kastaun et al. (2016) is relevant in this respect.

3.7 Mass versus equatorial compactness and a criterion for prompt collapse

The almost identical slope of the turning point lines in the M(R_e) plot of Fig. 27, implies that the threshold mass to collapse will be attained at about the same value of the ratio M/R_e, for all three EOS considered. In correspondence to the usual definition of the compactness of a non-rotating star C = M/R, we define the ratio C_e = M/R_e as the equatorial compactness.

It is interesting that in simulations with the SFHo EOS by Kastaun & Ohme (2021) a ‘bulk’ compactness of 0.31 for the remnant 1ms before collapse is reported, with |$C_\text{max}^{\rm TOV}=0.295$| for the maximum mass TOV model for this specific EOS. In Kastaun et al. (2016) it was found that the core of the remnant has a mass profile that resembles that of a TOV solution, with similar findings echoed by Ciolfi et al. (2017). The latter study went on to suggest the conjecture that merger remnants that do not admit a TOV core equivalent⁸ promptly collapse to a black hole. It will be interesting to explore these prospects with equilibrium modelling in future work, now that we have provided proof of concept that the Uryu + law can capture the Ω ∼ ω feature in the core of remnants, as it has been observed in simulations. We note that the physical quantities of our remnant models are collectively reported in Tables 5, 6, and 7 for EOS APR, DD2, and MPA1.

4 DISCUSSION AND OUTLOOK

In this study, we followed up on our recent work (Iosif & Stergioulas 2021) and expanded our investigation of relativistic equilibrium models, using the 4-parameter differential rotation law of Uryū et al. (2017) and cold, tabulated EOS of high-density matter. We constructed sequences of merger remnants, taking advantage of empirical relations for the angular momentum at merger, derived through numerical simulations. Comparing the rotational profiles Ω(r_c), the Ω(F) profiles, the stellar surfaces and the meridional rest mass density profiles of our remnant models between the three EOS we employed, we found that all respective profiles studied are qualitatively similar when the rotation law parameters are held fixed to different values. Furthermore, choices for the rotation law parameters that yield the same Type of solutions (i.e. either Type A or C) also show qualitative similarities in profiles even when varying the parameter λ₁.

In addition, we constructed constant angular momentum sequences. From the intersection of the line connecting the turning points of J-constant sequences and the merger remnant sequence, we were able to reproduce the threshold mass to prompt collapse, M_thres with a relative difference of only |$\sim 1{{\ \rm per\ cent}}$|⁠, compared to accurate binary neutron star merger simulations. We stress that this was achieved using the same rotation law, i.e. the same values for the rotation law parameters {λ₁, λ₂} for all three EOS used in this study.

Our investigation of the threshold mass to collapse, using equilibrium models, also points towards a possible connection between the equatorial compactness C_e = M/R_e at the threshold and the maximum compactness of non-rotating models |$C_\text{max}^\text{TOV}$|⁠. This lends support to the conjecture by Ciolfi et al. (2017), that merger remnants collapse if their relatively cold and slowly rotating inner region does not admit a stable TOV equivalent.

Another key prediction of binary neutron star merger simulations is a relatively slowly rotating inner region, where the angular velocity Ω (as measured by an observer at infinity) is mostly due to the frame dragging angular velocity ω. In our investigation of the parameter space of the Uryu + rotation law, we naturally find quasi-spherical (Type A) remnant models with this property. Both the density distribution and the angular velocity profile of these models have striking similarities with merger remnants produced in simulations (except for the low-density regions, since we neglect thermal effects).

In a forthcoming study, we plan to take the next step in this program and construct equilibrium models with finite-temperature EOS, using temperature and electron fraction profiles extracted from simulations of BNS mergers. This should allow us to isolate the effect of the thermal state of the remnant on the key properties that we discuss above. Based on Kaplan et al. (2014), Camelio et al. (2019), Camelio et al. (2021), Chakravarti & Andersson (2020) we do not expect dramatic deviations for bulk properties of the remnants (such as their mass and angular momentum) when thermal effects are included, whereas the corresponding angular velocity profiles could exhibit a larger effect depending on the different thermal treatments (Raithel, Paschalidis & Özel 2021). One can still expect an increased radius, due to the added thermal support (see e.g. Koliogiannis & Moustakidis 2021), and consequences on longer time-scales, such as convective instabilities.

In our study, we construct models of merger remnants that have angular momentum equal to the angular momentum at the time of merger, as extracted from simulations. Due to the excitation of non-axisymmetric oscillations in the remnant (as well as non-axisymmetric features, such as spiral arms), the GW emission in the post-merger phase will result in a reduction of the angular momentum with time. The remnant’s evolutionary path due to various dissipative processes taking place, could be tracked by constructing evolutionary sequences of equilibrium models. To a first approximation, one could neglect mass-losses and construct evolutionary sequences keeping the rest mass of the remnant, M₀, fixed in time.

With this viewpoint and following the related discussion in Kaplan et al. (2014), we show in Fig. 29 the rest mass M₀ versus the maximum energy density ϵ_max for the MPA1 EOS and the rotation law with {λ₁, λ₂} = {1.6, 1.0}. Analogously to Fig. 10, the blue dashed lines represent constant angular momentum sequences, with sequences for J = {3.0, 3.5, 4.0} as well as those reported in Table 7.

We note that no remnant models are located for these lower J-constant sequences, since their lower masses (implying a binary total mass |$M_\text{tot} \lt 2.2 \, \mathrm{M}_\odot$|⁠) would make them irrelevant for this study. The arrows in Fig. 29 approximate the evolutionary path of remnants as they lose angular momentum (assuming constant rest mass). A remnant will evolve towards larger densities and higher compactness, as its angular momentum decreases. In addition, the rotation law may also evolve in time. Three representative evolutionary paths are shown, resulting in a delayed collapse to a black hole on a short time-scale (when the remnant is hypermassive); a delayed collapse on a longer time-scale (when the remnant is supramassive and can still exist as a uniformly rotating model); and a stable remnant (when its mass is smaller than the maximum mass for non-rotating models).

The description of post-merger remnants as quasi-equilibrium models has applications in the interpretation of the post-merger GW spectrum; in the study of the threshold mass to prompt collapse; in the construction of universal or empirical relations between properties of the remnants and properties of non-rotating models; and in modelling of processes taking place on longer time-scales that are relevant for multimessenger follow-up studies of GW detections (Ciolfi et al. 2021; Rosati et al. 2021). We plan to elaborate further on these aspects in future work.

ACKNOWLEDGEMENTS

We thank Andreas Bauswein, John Friedman, Dimitrios Papadopoulos, and Christos Tsagas for helpful discussions and comments. We also thank the anonymous referee for providing useful comments that improved the final manuscript. PI gratefully acknowledges support by a Virgo-EGO Scientific Forum (VESF) PhD fellowship. NS acknowledges support by the COST actions CA16214 ‘PHAROS’, CA16104 ‘GWVerse’, and CA18108 ‘QG-MM’. The authors gratefully acknowledge the Italian Instituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS), and the Netherlands Organization for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium.

DATA AVAILABILITY

The data underlying this article are available in the article.

Footnotes

REFERENCES

APPENDIX: COMPARISON OF ANGULAR MOMENTUM EMPIRICAL RELATIONS

In this appendix, we present a more detailed comparison for the empirical relations (7) and (8) discussed in Section 2.5 that were used to estimate the angular momentum at the time of merger.

In Fig. A1 the two angular momentum empirical relations (7) and (A1) are shown for the DD2 EOS. An M_tot ranging from |$2.2 \, \mathrm{M}_\odot$| to |$3.3 \, \mathrm{M}_\odot$| is assumed. Note that the threshold mass for prompt collapse to a black hole for this particular EOS is |$M_\text{thres} \simeq 3.325 \, \mathrm{M}_\odot$| as reported from numerical simulations by Bauswein et al. (2020, Supplementary Material). In equation (A1) we have set R_NS = 13.1 km, which is a typical circumferential radius value for progenitor neutron stars constructed with the DD2 EOS (see Fig. 1). The two curves show very good agreement, deviating for the most part at the 2.5 per cent level and at the 5 per cent level only for very low mass remnants.

To investigate further the two angular momentum empirical formulas, we consider six zero-temperature, tabulated EOS, namely APR, DD2, and MPA1 already mentioned in Section 2.4 with the addition of SFHx (Steiner, Hempel & Fischer 2013), SLy4 (Douchin & Haensel 2001), and WFF1 EOS (Wiringa, Fiks & Fabrocini 1988). For each EOS, we consider a remnant mass range starting from |$2.2 \, \mathrm{M}_\odot$| and increasing with a step of |$0.1 \, \mathrm{M}_\odot$| up to approximately the threshold mass value M_thres, as the latter is reported in Bauswein et al. (2020). Constructing the corresponding non-rotating progenitor models with M_NS = M_tot/2, we then use the empirical relations (7) and (8) to get predictions for the value of J_merger from each formula.

Tables A1, A2, A3, A4, A5, and A6 report the J_merger estimates from the Bauswein & Stergioulas (2017) and Lucca et al. (2021) expressions. For the EOS DD2, MPA1, SFHx, and SLy4 excellent agreement is noticed, as the relative difference ranges from |$\sim 0.5{{\ \rm per\ cent}}$| up to |$\sim 5{{\ \rm per\ cent}}$|⁠. For the APR EOS, somewhat larger deviations (⁠|$\gt 6{{\ \rm per\ cent}}$|⁠) are noticed for remnants with masses over |$2.6 \, \mathrm{M}_\odot$|⁠, with the Bauswein & Stergioulas (2017) approximate formula predicting larger angular momentum values. The picture is similar for the WFF1 EOS, with the deviations augmented to the 10 per cent level.

Overall, the agreement between the empirical formulas from Bauswein & Stergioulas (2017) and Lucca et al. (2021) at a level of a few per cent is notable, given that (i) they are independent analyses of different numerical simulations data sets and (ii) the formula appearing in Bauswein & Stergioulas was a benchmark for the DD2 EOS with an added approximate treatment to describe other EOS. In fact, since one cannot be certain that the specifics of the angular momentum extraction were performed in an identical way between the two works, the actual agreement could prove to be even better. In any case, an accuracy of a few per cent is adequate to get useful insights into merger remnants properties with equilibrium models.