Axisymmetric models for neutron star merger remnants with realistic thermal and rotational profiles

Merging neutron stars are expected to produce hot, metastable remnants in rapid differential rotation, which subsequently cool and evolve into rigidly rotating neutron stars or collapse to black holes. Studying this metastable phase and its further evolution is essential for the prediction and interpretation of the electromagnetic, neutrino, and gravitational signals from such a merger. In this work, we model binary neutron star merger remnants and propose new rotation and thermal laws that describe post-merger remnants. Our framework is capable to reproduce quasi-equilibrium configurations for generic equations of state, rotation and temperature profiles, including nonbarotropic ones. We demonstrate that our results are in agreement with numerical relativity simulations concerning bulk remnant properties like the mass, angular momentum, and the formation of a massive accretion disk. Because of the low computational cost for our axisymmetric code compared to full 3+1-dimensional simulations, we can perform an extensive exploration of the binary neutron star remnant parameter space studying several hundred thousand configurations for different equation of states.


I. INTRODUCTION
With densities substantially exceeding those in atomic nuclei, neutron stars (NSs) provide an interesting 'astrophysical laboratory' to probe matter under the most extreme conditions and they can deliver physical information that complements other ongoing efforts to understand nuclear matter [1,2]. NSs originate in supernova explosions or binary neutron star (BNS) mergers [3]. In either case, they are hot and differentially rotating in the first minute of their lives [4,5]. Because of the growing possibilities of detecting them via gravitational wave interferometers and in the whole electromagnetic spectrum (from radio to gamma rays [6]), and because they involve nuclear matter at densities and temperatures that cannot be probed in terrestrial experiments, BNS remnants have been carefully investigated in a number of recent studies, e.g., [7][8][9][10][11][12].
The physical realism of 3+1 numerical relativity simulations has enormously increased over recent years, but realistic simulations come at the price of several hundred thousand core hours on supercomputers per physical millisecond, which makes an efficient exploration of the remnant parameters impossible. Moreover, studies that focus on the exploration of the microphysics, such as the effects of neutrino oscillations [13][14][15], need physically motivated background models but usually cannot afford at the same time a 3+1 numerical relativity approach. For these reasons very fast, yet still physically reliable, axisymmetric models of newly formed merger remnants are needed.
The vast majority of NS studies neglect differential rotation and assume rigid rotation. The first model of a NS in differential rotation made use of the so-called j-constant rotation law 1 [16,17], which is a good qualitative description of the proto-neutron star formed in a core-collapse supernova, where the core rotates faster than the envelope. In order to improve on these approximations, Uryū et al. [18] proposed a new model for the rotation profile of a BNS merger remnant that mimics the output of dynamical simulations [5,[19][20][21], where the angular velocity reaches a maximum in the envelope and approaches Keplerian rotation at large radii. Since then, other authors used Uryū and collaborators' model [22][23][24]. However, a proper inclusion of the thermal profile of the BNS merger remnant, that can reach temperatures up to a hundred MeV, has not been done yet. Moreover, until recently, hot NS models had been obtained through the so-called effectively barotropic approximation, where all thermodynamical quantities were put in a one-to-one relation [25]. This is a strong assumption for a remnant, that is expected to be baroclinic, i.e., not effectively barotropic [4]. Recently, Camelio et al. [26] developed a technique to obtain a stationary, hot, differentially rotating, baroclinic NS model, opening the way to a larger class of thermal and rotational profiles.
Modeling BNS merger remnants with stationary codes is an important complementary approach to full hydrodynamical simulations, since it allows for a much faster and wider exploration of the possible parameter space. In addition, stationary configurations can be used as initial profiles for dynamical simulations. Last but not least, the study of stationary configurations provides important indications on stellar stability [25,[27][28][29][30][31][32]. This is impor-tant because unstable stars are more likely to be observed through gravitational, neutrino, and electrodynamic radiation [e.g., 3, 6,33,34], allowing for an in-depth study of the involved physics.
In this work, we first develop a model for the stationary remnant of a BNS system at ∼ 10-50 ms after merger, which is differentially rotating, hot, and baroclinic (Sec. II). In particular, we propose new rotation and thermal laws for the remnant and apply the baroclinic formalism developed in Camelio et al. [26]. We then explore the model parameter space and discuss the remnant stability with simple heuristics (Sec. III). We conclude in Sec. IV. In Appendix A, we provide details of our numerical implementation. The parameter space exploration results and the profiles of the most realistic stellar models found are available to the community on Zenodo [35].
Unless otherwise specified, we set c = G = M = k B = 1. Our code unit for lengths approximately corresponds to 1.477 km, that for angular velocity to 32.31 kHz 2π rad, that for energy to 1.115 × 10 60 MeV, and that for time to 4.925 µs. Moreover, the saturation density is ρ n = 4.339×10 −4 and the neutron mass is m n = 8.423×10 −58 .

A. Equation of state
The equations of state (EOSs) adopted in this work are piecewise polytropes with a crust [36] and a thermal component [26]: where , ρ, s are respectively the total energy density, the rest mass density, and the entropy per baryon, a i , k i , Γ i are cold piecewise polytropic parameters valid in a given density range ρ i−1 < ρ < ρ i and are obtained by fits [36], Γ th = 1.75 is the thermal exponent and we set its value so that it is in the range expected for the high-density part of the EOS [37,38], and k th is the thermal constant and its value is determined for each EOS so that the thermal pressure at ρ = 2ρ n and s = 2 k B is 30% of the cold pressure. This value has been chosen after inspecting tabulated EOSs and could be easily adjusted for further studies, if needed. We consider a subset of EOSs from Read et al. [36] that fulfill the most recent radius and maximum mass constraints obtained from nuclear physics and astrophysical observations [39]: ALF2 [40], SLy [41], APR4 [42], and ENG [43], see Appendix A 1.

B. Euler equation
We determine the NS configuration with our version [26,32] of the XNS code [44,45]. The code assumes stationarity (and hence axisymmetry), circularity (and hence the absence of meridional currents), and conformal flatness [46,47]. The conformal flatness assumption does not change the theory of the modeling of the neutron star and its stability described in this section; however, the exact values of the total stellar quantities like mass and angular momentum may vary at most up to a few percent with respect to the values obtained in full General Relativity [24,48]. This level of precision is acceptable for this initial study.
It is possible [26] to cast the Euler equation in a form that is reminiscent of thermodynamical equations dQ(p, F ) = dp by defining the potential where p is the pressure, h = + p the total enthalpy density, r, θ are respectively the quasi-isotropic radius and polar angle coordinates, Ω = u φ /u t is the fluid angular velocity seen from infinity (u is the fluid four-velocity), F = u t u φ is the redshifted angular momentum per unit enthalpy and unit rest mass [24], α is the lapse function, and Γ = αu t is the Lorentz factor with respect to the zero angular momentum observer. From Eqs.
(2)-(3) it follows that the angular velocity Ω and the enthalpy density h can be obtained by differentiation: The advantage of using the potential Q to define the stellar model is that in this way we can obtain "baroclinic" configurations [26], that allow for a more realistic representation of merger remnants [4] than the commonly used "effectively barotropic" approximation. In an effectively barotropic model, one thermodynamical variable fixes all the other ones, while this is not true in a baroclinic model. Note that we choose a version of the potential Q that depends on F instead of Ω since in a BNS merger remnant the profile of the angular velocity is not monotonic [18,21,26]. Our model for a BNS merger remnant is defined by the following potential: where b is the "baroclinic" parameter, H(p) is the "heat function", G(F ) is the "rotation law" 2 , ρ 0 is a parameter equivalent to the central density 3 ,s(ρ) is (one version of) the "thermal law", namely a one-to-one relationship between the thermodynamical quantities, and p and ρ(p) the total derivative ofp(ρ) = p ρ,s(ρ) and its inverse, respectively.  [49].
To complete the definition of our model, we must choose the rotation and thermal laws. For the rotation law, we propose where G 0 , Ω 0 , Ω M , F 0 , and σ are free parameters. This rotation law is smooth (its second derivative is continuous), it has an easy analytical form, a minimum (resp. maximum) at the center (resp. at F 0 ), and it is Keplerian at large radii 4 . When the star is effectively barotropic (b = 0), the derivativeΩ = −G (F ) is equal to the angular velocity profile, see Fig. 1a and cf. Eqs. (4) and (6). In this case, Ω 0 and Ω M are the axial and maximum angular velocities, the latter reached off-axis exactly at F 0 . F 0 (resp. σ −1 ) is the scale of the variation for the low (resp. high) angular momentum part of the rotation law. To reduce the number of free parameters, we assume that G 0 = 0, which implies that ρ 0 and Ω 0 are the central density and axial angular velocity 5 also in the baroclinic case [26], and that G is continuous in F 0 , which leads to During the first tens of milliseconds after the merger, the remnant is not isentropic [4]: temperatures and entropy increase for decreasing density up to a critical value that in the nonbarotropic case they consider (b = 0) is equivalent to the angular velocity Ω.
where the temperature peaks. At lower densities the temperature decreases adiabatically, while the entropy per baryon keeps increasing, but with a lower rate. This behavior can be reproduced with our EOS assuming the following thermal law (see Fig. 1b-c): with Γ th − 1 > Γ T > 0, which implies where ρ M is approximately the peak density for the temperature and ρ L is a density scale, k s is a multiplicative constant that sets the scale of the entropy, and Γ T is the temperature polytropic index at lower density. Following the description of Perego et al. [4] of the BNS merger remnant at ∼ 10 ms after the merger, we set ρ m = ρ L = ρ n and Γ T = 2/3 (i.e., adiabatic expansion in the envelope). At later times (20-30 ms) and in the low-density region ρ < 10 −4 ρ n , this value is expected to decrease to Γ T = 1/3 [4].

A. Search
For each EOS, we run about 100,000 simulations in order to explore the parameter space, varying the following six parameters with a uniform distribution: • central density ρ 0 = [2; 10] ρ n .
• maximal-to-axial angular velocity ratio Ω M /Ω 0 = [1; 2]. • rotation law scale σ −1 = [1.5, 10] km (we report it in km because it can be interpreted, approximately, as the radial scale of the rotation distribution at large distance from the rotation axis).
As already discussed, see Sec. II B, the other parameters are set as follows: G 0 = 0, F 0 from Eq. (9), ρ M = ρ L = ρ n , and Γ T = 2/3. The values and ranges of the parameters are chosen to approximately reproduce the models evolved by Hanauske et al. [5] and Perego et al. [4]. In particular, we set b < 0 so that there is a hot ring in the equatorial plane instead of two hot caps in the polar regions and its range is set to resemble the models of Perego et al. [4] and Kastaun et al. [21], and to include the effectively barotropic model as special case (b → 0). The numerical details of how we find our solutions with a modified version of the XNS code are reported in Appendix A 2.
We remark that time evolution of the BNS remnant can be mimicked by varying the free input parameters of our model. In fact, shortly after merger, the central density ρ 0 increases on a t 10 ms timescale [4,5], b → 0 on a t 50 ms timescale [4], and Ω M /Ω 0 increases and Ω 0 decreases on a t 10 ms timescale due to angular momentum transfer from the central region outward [5]. On a longer timescale, if the remnant does not collapse to a black hole, we expect that k s → 0 on a t 10 s timescale due to the loss of entropy caused by neutrino emission, so that the central density ρ 0 keeps increasing due to cooling, and that Ω M /Ω 0 → 1 and σ → 0 as the star approaches rigid rotation due to neutrino diffusion and magnetic viscosity (both of which we do not include). Whether the axial angular velocity Ω 0 increases or decreases on a longer timescale depends on the total angular momentum loss by neutrino emission and magnetic braking [33,34], on that gained by accretion, and by the evolution of the stellar moment of inertia. However, it is plausible to assume that the BNS merger remnant is spun up as it happens for a proto-neutron star [50].
Before discussing the results, we remark that only ∼ 7-8% of the parameter combinations in the searches gives a valid solution of the Einstein and Euler equations. The failure of a particular parameter combination may be due to the physics (e.g., the mass shedding limit has been exceeded) or to the numerics (i.e., the code is not stable enough; a "false negative"). We increased the stability of the code by choosing physically motivated parameters and by slowly increasing the rotational and thermal content of the star at the beginning of the iterative process (see Appendix A 2). However, it is unavoidable that a fraction of the unsuccessful runs might consist of false negatives. On the other hand, the successful configuration are physical in the sense that they are solutions of the Einstein and Euler equations, but despite our efforts of realistic modeling we cannot be sure that they all approximate the result of dynamical evolution of mergers. For example, a small number of successful parameter combinations (of the order 10) results in stellar models with gravitational mass M > 4. Considering BNS population scenarios, such high masses are astrophysically unlikely (but not impossible for rapidly rotating models with extremely stiff EOSs [51]) for BNS merger remnants and we will exclude these configurations from the following analysis.
Unless otherwise stated, we will consider the ALF2 EOS in this section. The reason is that we can reliably invert the EOS and obtain the rest mass density ρ and entropy per baryon s from the pressure p and the enthalpy density h only for this EOS (see Appendix A 1 for details). We checked that the other quantities follow the same qualitative trends of the ALF2 EOS, see for example Fig. 3.
The parameters and stellar quantities of the successful configurations found in the search can be downloaded from Zenodo [35].

B. Stellar properties
By exploring several thousand configurations, we find some physically motivated correlations in the parameters of the successful configurations, see Fig. 2: The maximum of the axial angular velocity Ω 0 increases with density ρ 0 (Fig. 2a), since gravity is stronger and it is possible to reach faster rotation without mass shedding, cf. Fig. 4b. Similarly, the maximum of the entropy scale k s increases with increasing density ρ 0 (Fig. 2b; the other EOSs reach k s = 9 with a similar trend of ALF2). The maximum of the axial angular velocity Ω 0 and of the rotation scale σ −1 are greater for smaller rotation ratio Ω M /Ω 0 (Figs. 2c-d) due to the necessity for the equatorial angular velocity to be lower than the Keplerian frequency in order to avoid mass shedding. Remarkably, the position of the maximum of the angular velocity F 0 (which also serves as a scale for the inner part of the rotation law) is correlated with the baroclinic constant b, while it is not correlated with the central density ρ 0 (Figs. 2e-f) and with k s .
For a given central density, the gravitational mass M of our BNS merger remnant model is larger than the nonrotating mass and can even be larger than the cold rigidly rotating Keplerian one, see Fig. 3. In Fig. 4 we show some trends of other stellar quantities. Some trends are obvious and expected: the stellar angular momentum J grows with the axial angular velocity Ω 0 (Fig. 4a), the maximum of the Keplerian angular velocity Ω kep grows with the central density ρ 0 , and the average entropy per baryon S/M 0 grows with the entropy scale k s . The circumferential radius R cir grows with entropy scale k s (due to an increasing thermal pressure, Fig. 4e) and when the equatorial angular velocity Ω eq approaches the Keplerian one Ω kep (since the configuration approach mass shedding, Fig. 4f). Assuming that the configuration collapses to a black hole, one can estimate the mass of the accretion disk that remains outside of the innermost stable circular orbit (see discussion in Sec. III D). The disk mass M disk also grows when the equatorial angular velocity approaches the Keplerian one, as expected. It is worth pointing out that, when present, the estimated mass of the accretion disk is of the order of that expected from simulations, e.g., Refs. [5,39,52,53], and much larger than the one that is expected from a single rotating NS [31,32].

C. Stability
The solutions that are found for a given EOS are not necessarily dynamically stable. There are many types of instabilities that may be present [for a review see e.g. 54], and whether or not a particular one is relevant depends on how its associated timescale compares with the timescale of the viscous processes at work. In this paper, we will consider a non-comprehensive set of possible instabilities that may be present in the remnant.
Low T /|W |-instability: The dynamical study of dif-ferentially rotating configurations allowed the discovery of the so-called "low-T /|W | instability" [22,55,56]. The low-T /|W | instability sets in when an oscillation mode co-rotates with the matter in a point of the star. Since in a BNS merger remnant the angular velocity is not monotonic with the radius, it is possible for an oscillation mode to co-rotate with the matter in two points [22,23]. Performing the numerical evolution in General Relativity of an initially cold remnant with a rotation law from Uryū et al. [18], Xie et al. [23] found that this instability is present for the relatively low value of T /|W | = 0.16, where T is the kinetic energy (not to be confused with the temperature) and W = M p + T − M is the gravitational binding energy (M p is the proper mass). Similarly, making use of Newtonian gravity, assuming the Cowling approximation, and exploring a larger number of remnant configurations, Passamonti and Andersson [22] found that this instability may set in for T /|W | 0.02 and as T /|W | grows it initially becomes more relevant until the mode stabilizes to a specific value of T /|W |. We find that T /|W | grows for increasing axial angular velocity Ω 0 (Fig. 5a), and that the maximum of T /|W | decreases with increasing entropy scale k s (Fig. 5b) and increasing rotational scale σ −1 (Fig. 5c). The anti-correlation between T /|W | and k s is not mediated through the central density ρ 0 since both Ω 0 and k s increase with increasing ρ 0 , cf. Figs. 2a-b. We interpret the anti-correlation of k s and T /|W | with the fact that, increasing the thermal pressure, the star is less strongly bound. We conclude that a larger entropy content contributes in stabilizing the star against the low-T /|W | instability.
Convective instability: The convective instability has a local character and sets in when a displaced fluid element is accelerated away from its equilibrium position. In a hot, rotating star, the forces that are applied on a fluid element are gravity, buoyancy (due to the pressure gradient), and the centrifugal force [see e.g. 57].
In a non-rotating and hence spherical star, convective instability is driven by buoyancy. In this case, necessary conditions for convective instability are a non-barotropic EOS and entropy (or composition) gradients. For nonrotating NSs, the onset of convective instability is controlled by the Schwarzchild criterion [58], that is, a star is convectively unstable when the Schwarzschild discrim-inantS(r) is negative, wherer is the Schwarzschild radius and c 2 s the speed of sound. As pointed out in Camelio et al. [26], for our EOS this is identical to or, equivalently (since in our case Γ th < Γ i ), a star is unstable against convection if sgn(ρ c − ρ) · sgn ds dr < 0, where the critical density ρ c happens to be the same critical density for the EOS inversion (see Appendix A 1; the value of ρ c for each EOS is reported in Table II and marked with a vertical line in Fig. 3). In our case, since the thermal law for the effective barotropic case is such that ds/dρ > 0 [we set Γ T = 2/3 in Eq. (10)], Eq. (14) tells us that if the density decreases monotonically from the center outward, then the star is stable in the region with ρ < ρ c and unstable for ρ > ρ c . On the other hand, in isentropic stars the driver for convective instability is the centrifugal force. When the isentropic star is differentially rotating, a necessary criterion for convective stability is [54,59] dj 2 (r, π/2) dr > 0, (15) where the square of the specific (per unit mass) angular momentum j = h u φ /ρ is differentiated along the equatorial plane, θ = π/2. As shown in the top panel of Fig. 7, this criterion is generally respected with our differential rotation law.
Having a configuration that is differentially rotating, non-isentropic, and baroclinic (namely non effectively barotropic) at the same time means that the simple criteria (12) and (15) are no more valid. This is due to the fact that not only the gravitational force is no more balanced by the buoyant force alone, but also to the fact that the three forces are not necessary parallel. However, due to the qualitative nature of our discussion, we will still make use of criteria (14) and (15) to allow for this simple remark on the remnant stability: that an increase of ρ 0 favors a buoyancy-driven convective instability, because configurations with ρ 0 > ρ c (right of the dotted line in Fig. 3) are convectively unstable (at least in a part of the star), while configurations ρ 0 < ρ c (left of the dotted line in Fig. 3) are convectively stable if the density decreases monotonically with the radius everywhere (and the entropy increases with the radius). Note that convective instability has been found 30-50 ms after merge in numerical simulations [60,61], and some of our models do have negative entropy gradients at intermediate radii. We remark moreover that this already approximate consideration is valid only for the simplified EOSs we are considering. In a more realistic EOS, the value of Γ th may be density-dependent and the simple non-rotating convective instability criterion we derived, see Eqs. (13)- (14), should be revised.
Axisymmetric secular stability: Another type of instability is the "secular instability". It sets in when a configuration evolves to a similar one with lower energy (i.e. a stabler one) on a timescale longer than the hydrodynamical one. Here we are concerned with axisymmetric instabilities, which can be determined simply by studying the stationary (axisymmetric) configurations. In practice, a configuration defined by a set of parameters is stable if all configurations close in the parameter space with the same baryon mass M 0 , angular momentum J, and total entropy S, have a greater grav-itational mass M . For rigidly rotating, isentropic NSs, secular stability can be checked with the turning point criterion [25,[27][28][29], that is, a star becomes secularly unstable when where ρ 0 is the central density. In the case of a cold and non-rotating NS, secular instability implies and is implied by instability against dynamical perturbations [62], while in general, for a rigidly rotating and isentropic star, secular instability is a sufficient but not necessary condition for instability against dynamical perturbations [e.g., 30].
In the general case we are interested in, namely differentially rotating, non-isentropic, and nonbatropic NSs, the turning point criterion [Eq. (16)] cannot be applied because the number of free parameters is greater than the number of conserved quantities (i.e., M 0 , S, J [28]). Moreover, two configurations close in the parameter space may be not connected by any dynamical evolution, and therefore we cannot apply the definition of secular instability. In the case of the ALF2 EOS, for which we can compute M 0 and S, we tried anyway to look at the variation of M with respect to the parameters, keeping M 0 , J, S constant. Unfortunately, we were not able to draw clear conclusions from this analysis and further work is required 6 .

D. Selected models
In this section we show the stellar profiles for three selected models, chosen according to the following criteria: A: A model that is a plausible outcome of a BNS merger. B: The model with M < 4 that is closest to an effectively barotropic configuration, namely that with the baroclinic parameter b closest to zero, to be compared with model A. C: The model with the greatest disk mass and M < 4. The specific model parameters and properties are summarized in Table I and their rest mass density, temperature, and angular velocity distribution are shown in Figs. 6-7. Model A is a realistic BNS merger remnant. Its mass M , angular momentum J, and central density ρ 0 are in the expected range, e.g., [63,64]. Its shape is qualitatively similar to that obtained in simulations, (cf. Fig. 13 6 It would be interesting, in the future, to find a physically motivated parametrization of the star with a number of free parameters less or equal to the number of global constraints on the stellar evolution (i.e., 3), so that the the turning point criterion [27,28] would be applicable, and to compare this secular stability with dynamical simulations as done in Takami et al. [30] and Camelio et al. [32]. in Perego et al. [4]). The temperature forms a hot ring in the equatorial plane (cf. Fig. 7 of Perego et al. [4]); the temperature profile is continuous but not smooth due to the fact that the EOS is piecewise defined. The angular velocity curve peaks 3-4 km from the rotational axis (cf. Fig. 5 of Hanauske et al. [5]).
Baroclinicity is fundamental in order to obtain the right thermal distribution. This can be realized comparing the profiles of models A and B (the latter being almost effectively barotropic): in model A the density and temperature isocontours are not parallel and this permits the existence of the hot equatorial ring, while in model B they are parallel and as a consequence the temperature profile has an onion-like shape. This is a consequence of baroclinicity [26].
Between the three chosen models, model C has the biggest circumferential radius and its equatorial angular velocity is the closest to the Keplerian one. Unsurprisingly, a significant amount of matter with large specific angular momentum is present. In case of black hole formation, this matter could form a disk. We followed the approach of Margalit et al. [31] (see also Shapiro [65] and Camelio et al. [32]), namely we computed the baryon mass of the matter whose angular momentum per unit mass j is larger than that of the innermost stable circular orbit of a black hole with the same mass and angular momentum of the original system 7 . This is equivalent to assume that there is no angular momentum transfer or loss during the collapse and that dynamical effects like shocks play no role, which is clearly not true [e.g., 33,34], but at the same time it is a first order approximation that allows us to make a semi-quantitative estimations of the expected disk mass M disk . We found it to be M disk ≈ 0.4 M , which is substantially larger than that expected from the collapse of a marginally stable, Keplerian, rigidly rotating, cold NS [31,32]. The disk mass is in the range of what is expected from dynamical simulations (actually at the upper end of the expectations) [5,12,39,52,53]. With a large potential energy reservoir of 3.6 × 10 52 erg ( /0.05)(M disk /0.4M ), this configuration is a good candidate for launching a powerful short GRB, provided that the energy can be deposited in a low enough density environment. For the latter, one usually assumes that the central object needs to collapse to a black hole, but see [66] for the possibility to launch relativistic outflows in the presence of a central neutron star.
We provide the community with 10 realistic BNS remnant models (including model A) with the ALF2 EOS I: Parameters and stellar quantities of the models shown in Sec. III D. The EOS is ALF2. The quantities are: central density ρ0, axial angular velocity Ω0, entropy scale ks, angular velocity ratio ΩM /Ω0 and rotation law scale σ −1 , baroclinic parameter b, gravitational and baryon mass M and M0, stellar angular momentum J, average entropy per baryon mass S/M0, kinetic-to-gravitational energy ratio T /|W |, circumferential radius Rcir, Keplerian angular velocity at the equator Ω kep , angular velocity at the equator Ωeq, expected disk mass M disk , maximum temperature Tmax. and stellar properties in the ranges 2.4 < M < 3.5, 5.5 < J < 7.5, 1.5 < ρ 0 /ρ n < 4, Ω M /Ω 0 > 1.05, compatible to numerical relativity simulations, e.g., [67,68], plus models B and C. The stellar profiles and a python script to read them can be downloaded from Zenodo [35]. This dataset can be used as background models for microphysical studies or as initial conditions for dynamical evolution.

IV. CONCLUSION
In this paper we studied realistic stationary models for post-merger configurations after a BNS merger. We modeled the EOS with cold polytropic pieces [36] plus a thermal component as described in more detail in [26]. Our remnant model is controlled by the central density and other parameters that fix the rotational and thermal distributions. We explored a broad range of postmerger configurations and discussed their stability based on qualitative criteria.
In particular we have • introduced new rotation and thermal laws, Eqs. (8) and (10), that are motivated by numerical simulations [4,5]. • applied the technique recently developed by Camelio et al. [26] to BNS merger remnants. We obtained baroclinic (i.e., not effectively barotropic) configurations, which are more suitable to model merger remnants than the effectively barotropic ones [4].
• performed an extensive parameter space study in which we included the effects of differential rotation, temperature, and baroclinicity. Our main results are: • the central density ρ 0 , the axial angular velocity Ω 0 , and the thermal scale k s are the parameters that have the largest impact on the global remnant properties, see Fig. 4. • baroclinicity (implemented with the parameter b) is necessary to reproduce the thermodynamical profile of BNS merger remnants, in particular the existence of a hot ring in the equatorial plane [4,21], compare models A and B in Sec. III D. • the collapse of a BNS merger remnant to a black hole may generate a massive disk which could provide the central engine to launch a short gamma ray burst [31,32], see Figs. 4d and 7. • the increase of the central density ρ 0 may cause convective instabilities and the increase of the axial angular velocity Ω 0 may cause low-T /|W | instability. If no convection is present, an increased thermal content (k s ) seems to increase stability by reducing the maximal T /|W | that can be reached by the model. • we make available to the community the results of our parameter search and a set of realistic models [35]. The approach described here can be extended further: • an even more realistic description of the remnant physics, namely (i) the inclusion of composition in the model, (ii) the adoption of more realistic EOSs (for example the new piecewise parametrization of O'Boyle et al. [69] or a tabulated EOS), (iii) the addition of the magnetic field [see Ref. 70, for an example of proto neutron star studied in Newtonian gravity], and (iv) the use of a rotation curve that is truly Keplerian by construction, not only because it approaches the Keplerian trend at large radii (like in this work and in Uryū et al. [18]), but also the Keplerian frequency. • addition of physically motivated restrictions on the free parameters of the stationary remnant model to simplify the study of its stability and the test of these predictions with dynamical simulations [5,32] and/or a perturbative study [71,72]. der grant number Dnr. 107/16, the research environment grant "Gravitational Radiation and Electromagnetic Astrophysical Transients (GREAT)" funded by the Swedish Research council (VR) under Dnr 2016-06012 and by the Knut and Alice Wallenberg Foundation. We acknowledge stimulating interactions with the COST Actions CA16104 "Gravitational waves, black holes and fundamental physics" (GWverse) and CA16214 "The multimessenger physics and astrophysics of neutron stars" (PHAROS). The authors gratefully acknowledge the Italian Istituto 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.  [36] to implement them, including the SLy crust (we use only one index i = 1, . . . , 7 running from the crust to high density). We summarize the EOS properties in Table II. In order to recover ρ, s from h , p, we note that Eq. (B14) of [26] can be generalized to This means that one can find a valid stellar configuration without knowing the rest mass and entropy distributions inside the star.
In the case of a piecewise polytropic EOS such as Eq. (1), the degeneracy can in principle be more problematic than for a 1-piece polytropic one because there can be more than 2 valid ρ, s couples for given h , p. An analysis of the piecewise EOSs considered in this paper shows that this is not the case: the degeneracy is the same of the one-piece polytropic EOS of Camelio et al. [26], i.e., there is only one critical density ρ c in the range of the last high-density piece 8 , ρ c > ρ 6 . In Table II 8 In practice we studied Eq. (A1) in the range of each polytropic piece and checked whether the critical density ρ c,i , given by [cf.
Eq. (B15) of 26]: (i) exists (i.e., the RHS is positive) and (ii) lays in the correct range (i.e., ρ i−1 ≤ ρ c,i ≤ ρ i ). We found that, for the EOSs considered, the only piece with a critical density (respecting conditions i-ii) is the last one, and we therefore write ρc ≡ ρ c,7 . Note that if the RHS of Eq. (A1) is negative and ρc exists the only valid solution is the high-density one with ρ > Γ 7 −1 √ Γ 7 ρc. On the other hand, it can be that ρc < ρ < Γ 7 −1 √ Γ 7 ρc (i.e., the RHS of Eq. (A1) is positive), but the recovered s 2 is negative, that is the entropy is not physical and the only valid solution is the low-density one. Unfortunately, in general one cannot exclude one specific branch.   [26]. ρc is the critical density for the inversion.
we report the critical density for each EOS considered; apart for the ALF2 EOS, for which ρ c 45ρ n , the other critical densities lie in the range 4ρ n ρ c 6ρ n and are even lower than the central density of the maximal mass configuration of the spherical (nonrotating) model. Since the ALF2 has such a high value of ρ c , we can safely choose the low density branch of the solution like we did in Camelio et al. [26], while we cannot do the same for the other EOSs. For this reason, we are able to compute the total stellar rest mass M 0 , entropy S, and disk mass M disk only for ALF2. All other quantities, such the stellar gravitational mass M and angular momentum J, can be computed also for the other EOSs.

Neutron star
We used a modified version [26,32] of the XNSv2 code [44,45]; we refer the reader to the original papers for details on the implementation. The only difference with respect to our previous work [26] is that, at the beginning of the iterative procedure to determine the stellar configuration, we have slowly increased the thermal and rotational content of the star by varying k s and Ω 0 , in order to increase the stability of the numerical scheme. We set the following parameters in our code: • inner radial grid: boundary at r = 30 M , 3000 evenly spaced points, • outer radial grid: boundary at 1000 M , 3000 increasingly spaced points, • absolute tolerance of max(h ) for convergence: 10 −11 , • planar symmetry, • 50 points in the angular grid (in one of the hemispheres), • 500 relaxing iterations (see discussion above), • 30 Legendre polynomials. In order to implement the rotation law (8), we define two functions G 1 and G 2 : where Θ is the step function and We start by solving the system (3)-(4) with G ≡ G 1 . If F < F 0 , then we solve again the system with G ≡ G 2 . We found that in this way we increase the precision of the solution close to the maximum, F F 0 . Moreover, instead of solving the Newton-Raphson with F as independent variable, we found that it is numerically more stable to solve the equations using Ω as independent variable, even if the rotation law is defined in F . The heat function [Eq. (7)] is determined by numerical integration.