THE INFLUENCE OF THERMAL PRESSURE ON EQUILIBRIUM MODELS OF HYPERMASSIVE NEUTRON STAR MERGER REMNANTS

We use a set of eight EOS in this study. All EOS produce cold NSs in β-equilibrium that can have gravitational masses M_g above 2 M_☉. These include two EOS from Lattimer & Swesty (1991), the K₀ = 220 MeV and K₀ = 375 MeV variants (where K₀ is the nuclear compressibility modulus), denoted LS220 and LS375; the relativistic mean field (RMF) model EOS from Shen et al. (2011c), denoted HShen; two RMF models based on the NL3 and the FSUGold parameter set (Shen et al. 2011a, 2011b) denoted GShen-NL3 and GShen-FSU2.1; an unpublished⁷ RMF model based on the DD2 interaction denoted HSDD2; and two recent RMF model EOS fit to astrophysical measurements of NS masses and radii (Steiner et al. 2013), denoted SFHo and SFHx. All of these EOS are available in a common format for download from http://www.stellarcollapse.org.

The EOS of finite-temperature nuclear matter in nuclear statistical equilibrium (NSE) has contributions from a baryonic component (nucleons and nuclei), a relativistic electron/positron Fermi gas, a photon gas, and, if neutrinos are trapped, a neutrino gas. The Helmholtz free energies of these components add linearly, and the pressure is then the sum of the partial pressures and a function of baryon density ρ, temperature T, and electron fraction Y_e,

$$\begin{equation} P = P_\mathrm{baryon} + P_e + P_\gamma + P_\nu. \end{equation} \tag{ 1 }$$

While P_baryon varies between the employed EOS, we add P_e and P_γ using the Timmes EOS (Timmes & Arnett 1999) available from http://cococubed.asu.edu. In hot HMNSs, like in protoNSs, neutrinos are trapped and in equilibrium with matter. We include their pressure contribution to the EOS by treating them as a non-interacting relativistic Fermi gas with chemical potential $\mu _{\nu _i}$. For a single species of neutrinos and antineutrinos, the neutrino pressure in equilibrium is

$$\begin{equation} P_{\nu _i} = \frac{4\pi (k_B T)^4}{3(hc)^3} [F_3(\eta _{\nu _i}) + F_3(-\eta _{\nu _i})] \times \exp {\left(-\frac{\rho _\mathrm{trap}}{\rho }\right)}\,\,, \end{equation} \tag{ 2 }$$

where $\eta _{\nu _i} = \mu _{\nu _i} / (k_B T)$ is the neutrino degeneracy parameter. For HMNS conditions, all neutrino species are present, but ν_μ and ν_τ have $\mu _{\nu _i}=0$, since they appear only in particle–anti-particle pairs that have equal and opposite chemical potentials. For electron neutrinos we use $\mu _{\nu _e} = \mu _e + \mu _p - \mu _n$, for electron antineutrinos we use $\mu _{\bar{\nu }_e} = -\mu _{\nu _e}$. We include an attenuation factor exp (− ρ_trap/ρ) to account for the fact that neutrinos decouple from matter at low densities. We set ρ_trap = 10^12.5 g cm⁻³, which is a fiducial trapping density for protoNSs (e.g., Liebendörfer 2005). Taking the exact expression for the difference of the Fermi integrals from Bludman & van Riper (1978), we have the total neutrino pressure summed over all three species,

$$\begin{eqnarray} P_\nu &=& \frac{4\pi (k_B T)^4}{3(hc)^3} \left[\frac{21\pi ^4}{60} + \frac{1}{2}\eta _{\nu _e}^2\left(\pi ^2 + \frac{1}{2}\eta _{\nu _e}^2\right)\right] \nonumber\\ &&\times \exp {\left(-\frac{\rho _\mathrm{trap}}{\rho }\right)}\,\,. \end{eqnarray} \tag{ 3 }$$

We note that due to the neutrino statistical weight g = 1, for a single species of relativistic non-degenerate $\nu \hbox{{--}}\bar{\nu }$ pairs, the pressure is a factor of two lower than for e⁻–e⁺ pairs, since e⁻ and e⁺ have statistical weight (spin degeneracy) 2.

Figure 1 illustrates the contributions of the partial pressures to the total pressure as a function of baryon density ρ_b for neutron-rich HMNS matter at two temperatures, 0.5 MeV (a representative "cold" temperature) and 20 MeV (a representative "hot" temperature for HMNSs). For the 0.5 MeV EOS, we set the electron fraction Y_e by solving for ν-less β-equilibrium ($\mu _{\nu _e} = 0$). The resulting EOS describes ordinary cold NSs (at 0.5 MeV any thermal effects are negligible). For the 20 MeV case, we solve for Y_e by assuming ν-full β-equilibrium. We do so by making the assumption that any neutrinos produced during the merger are immediately trapped in the HMNS core, but stream away from regions below trapping density. The procedure is discussed in the next Section 2.2 and detailed in Appendix B.

Zoom In Zoom Out Reset image size

Near and above nuclear saturation density, ρ_nuc ≃ 2.6 × 10¹⁴ g cm⁻³ for the LS220 EOS, the baryon pressure is due to the repulsive core of the nuclear force and dominates in both cold and hot regimes. The thermal enhancement above ρ_nuc remains small even at 20 MeV. In the cold case, relativistically degenerate electrons (Γ = (dln P)(dln ρ)⁻¹ = 4/3) dominate below ρ_nuc. At 20 MeV, relativistic non-degenerate electron/positron pairs and photons (for both, P∝T⁴, independent of ρ_b; see, e.g., van Riper & Bludman 1977) are the primary contributors at low densities, while the baryon pressure is significantly thermally enhanced below nuclear saturation density and dominates above ∼10¹² g cm⁻³. The neutrino pressure is comparable to the degenerate electron pressure between ∼10^12.5 and 10¹⁴ g cm⁻³, but still subdominant to the nuclear component. The contribution of pairs and photons gradually becomes more important at all densities as the temperature increases. We note that for T = 0.5 MeV, the neutrino chemical potentials are all zero and the pressure of trapped neutrinos is 3 × (7/8) × P_γ, thermodynamically insignificant at T = 0.5 MeV.

The hydrostatic and rotational equilibrium equations that we solve in this study assume a barotropic EOS (P = P(ρ)) and do not provide constraints on thermal structure and composition (Y_e is the only relevant compositional variable in NSE). We must make some assumptions to be able to proceed and obtain P = P(ρ, T(ρ), Y_e(ρ, T(ρ))) for our general finite-temperature microphysical EOS. Old NSs in isolation are nearly isothermal (e.g., Prakash et al. 2001) and so are coalescing NSs until tidal heating becomes significant (e.g., Kochanek 1992; Lai 1994). During merger, the NS matter is shock-heated to tens of MeV and results of the few merger simulations that have been carried out with temperature-dependent EOS (e.g., Sekiguchi et al. 2011; Bauswein et al. 2010; Oechslin et al. 2007; Rosswog & Liebendörfer 2003; Ruffert & Janka 2001) indicate that the HMNS is far from being isothermal or isentropic. It has a very hot dense core with T ∼ 20–40 MeV surrounded by a lower-density cooler envelope/torus of 5–20 MeV, which may also be almost Keplerian and, hence, centrifugally supported. This result appears to be robust for equal-mass or near equal-mass NSNS systems (which may dominate the population; e.g., Lattimer 2012 and references therein). Mergers of non-equal mass systems in which the lower-mass NS is tidally wrapped around its more massive companion reach similar temperatures, but generally tend to have more mass at lower densities in the disk/torus (Oechslin et al. 2007).

There is no unique model/EOS independent mapping T = T(ρ), thus we must explore a variety of possibilities. In Figure 2, we contrast our set of temperature parameterizations with a T(ρ) profile obtained from a 1.35–1.35 M_☉ simulation using the HShen EOS by Sekiguchi et al. (2011) at ∼12 ms after merger. We consider very hot cores at 20, 30, and 40 MeV with cold envelopes (parameterizations c20p0, c30p0, and c40p0) and two parameterizations with very hot cores at 30 MeV and cool envelopes at 10 MeV and 5 MeV, c30p10 and c30p5, respectively. Since low-density regions have shorter neutrino cooling times, the c30p10 and c30p5 parameterization may represent early HMNSs, while the cold-envelope parameterizations c20p0, c30p0, and c40p0 may correspond to late-time HMNSs. Note that the c30p10 parameterization fits the temperature profile from the Sekiguchi et al. (2011) simulation quite well. Details on the functional forms of our parameterizations can be found in Appendix A. For the TOV case we also consider isothermal configurations as a limiting case.

Zoom In Zoom Out Reset image size

The choice of Y_e(ρ, T(ρ)) is equally difficult. Before merger, the NSs are in ν-less β-equilibrium (μ_ν = μ_e + μ_p − μ_n = 0). After merger, neutrinos are present. They are trapped in hot dense matter (μ_ν ≠ 0) and are streaming away from low-density regions. The equilibrium Y_e will shift and mixing due to non-linear hydrodynamics in the HMNS phase will distort any initial Y_e(ρ, T(ρ)) profile.

We deem the following prescription for Y_e to be the physically most sensible. We assume that the NSNS merger occurs so rapidly that the electron fraction Y_e of the ν-less β-equilibrium in the NSs becomes the trapped postmerger lepton fraction $Y_\mathrm{lep} = Y_e + Y_{\nu _e} - Y_{\bar{\nu }_e}$ above ρ_trap. Using the β-equilibrium condition with nonzero μ_ν, we solve for Y_e. At densities below ρ_trap we transition to Y_e given by ν-less β-equilibrium. Details of this procedure are given in Appendix B. We refer to this parameterization of Y_e as ν-full β-equilibrium. In addition and for comparison, we consider choices of constant Y_e = 0.1 and Y_e set according to ν-less β-equilibrium. We note that our parameterization of Y_e is ad hoc and cannot account for mixing and neutrino transport effects in the merger process. The right panel of Figure 2 depicts Y_e(ρ, T(ρ)) as obtained from the simulation of Sekiguchi et al. (2011) contrasted with Y_e profiles computed under the assumption of ν-less and ν-full β-equilibrium for various temperature parameterizations and for the T(ρ) as given by the simulation. None of the prescriptions fit the simulation-Y_e particularly well, which indicates that mixing and neutrino transport effects are important (but cannot be included here). The Y_e obtained using the temperature data from the simulation naturally fits best, in particular at low densities where neutrinos have decoupled from the matter and ν-less β-equilibrium holds.

In the top panels of Figure 3, we show the fractional pressure increase due to thermal effects as a function of baryon density for our set of temperature parameterizations for the LS220 EOS (left panel) and the HShen EOS (right panel) as two representative example EOS. We also distinguish between the choices of Y_e parameterization. For the parameterizations with cold "mantles" (cXp0), thermal effects are most important at densities near ∼ρ_nuc and quickly lose significance at lower and higher densities in both EOS. The thermal pressure enhancement is at most a factor of three (for the HShen) to five (for the LS220 EOS) for these parameterizations. The situation is different for the cases with hot plateaus, c30p10 and c30p5. For these, the thermal pressure is up to 20 times larger at low densities than predicted by the cold EOS. The Y_e parameterizations corresponding to ν-full and ν-less β-equilibrium yield qualitatively and quantitatively very similar results for both EOS.

Zoom In Zoom Out Reset image size

At low densities, the ν-full and ν-less β-equilibrium cases both lead to Y_e > 0.1 (see Figure 2). As a consequence, the pressure in the unrealistic Y_e = const. = 0.1, cXp0 parameterizations is lower than in the cold ν-less case at ρ_b ≲10^12.2 g cm⁻³. Due to the logarithmic scale of Figure 3, the graphs of cXp0 with Y_e = 0.1 start only there and the predicted pressure enhancement is higher than in the β-equilibrium cases, which lead to lower Y_e above ∼10^12.2 g cm⁻³ and below ∼ρ_nuc. In the cases with hot plateau (c30p10 and c30p5), thermal effects dominate over differences in Y_e at low densities. Finally, at ρ > ρ_nuc, where temperature effects are smaller, differences in Y_e become important. Since the nuclear component dominates there, lower Y_e corresponds to higher pressure (e.g., Lattimer & Prakash 2001) and both β-equilibrium cases yield Y_e > 0.1.

The lower panels of Figure 3 depict the relative contribution of the neutrinos to the total (hot) pressure in the HMNS temperature and Y_e parameterizations considered in this study. While there are clear temperature (see Equation (3)) and Y_e (through $\mu _{\nu _e}$) dependences, neutrino pressure plays only a minor role, making up at most ∼2% of the total pressure of the LS220 EOS. This is true also for the HShen EOS with the exception of the unrealistic Y_e = 0.1 case in which the neutrino pressure contribution grows to ≳10% of the total pressure at supranuclear densities.

Finally, we note that the temperature and Y_e prescriptions discussed here lead to regions that may be unstable to convection if not stabilized by a positive specific angular momentum gradient (e.g., Tassoul 1978). The spherically and axially symmetric equilibrium models that we construct in this study cannot account for convection and we leave an analysis of convective instability to future work.

We solve the TOV equation (e.g., Shapiro & Teukolsky 1983),

$$\begin{eqnarray} \frac{dP}{dr} &=& - \frac{G}{r^2}\! \left[\rho _\mathrm{b}\!\! \left(1 + \frac{\epsilon }{c^2} + \frac{P}{\rho _\mathrm{b} c^2}\right)\right]\!\!\! \left[ M_\mathrm{g}(r) + 4\pi r^3 \frac{P}{c^2} \right]\nonumber\\ &&\times \left[1 - \frac{2 G M_\mathrm{g}(r)}{r c^2}\right]^{-1}\!, \end{eqnarray} \tag{ 4 }$$

where r is the areal (circumferential) radius, ρ_b is the baryon density, is the specific internal energy, and M_g(r) is the gravitational mass enclosed by radius r, determined via

$$\begin{equation} \frac{dM_\mathrm{g}}{dr} = 4\pi r^2 \rho _\mathrm{b} \left[1 + \frac{\epsilon }{c^2} \right]\,\,. \end{equation} \tag{ 5 }$$

The baryonic mass is larger and given by

$$\begin{equation} \frac{dM_\mathrm{b}}{dr} = 4\pi r^2 \rho _\mathrm{b} \left(1-\frac{2GM_\mathrm{g}(r)}{rc^2}\right)^{-1/2}\,\,. \end{equation} \tag{ 6 }$$

We construct the TOV solutions using a standard fourth-order Runge–Kutta integrator on an equidistant grid with δR = 10² cm zones. After each integration sub-step, the EOS P = P(ρ_b) is inverted to obtain ρ_b. We use a variety of P(ρ_b) parameterizations: (1) T = const. (isothermal) with ν-full β-equilibrium above ρ_trap and ν-less β-equilibrium below, (2) T = const. with ν-less β-equilibrium, (3) T = const. with constant Y_e = 0.1, and (4) the phenomenological cXpX temperature parameterizations with ν-full β-equilibrium above ρ_trap and ν-less equilibrium below. We compute TOV solutions for all EOS and define the surface of the NS as the areal radius at which one of the following two conditions is true: (1) the pressure equals 10⁻¹⁰ of the central pressure; (2) the pressure predicted by the integration of Equation (4) drops below the lowest pressure value available in the equation of state table. The latter is not a limitation, because the high-density TOV configurations considered here have steep density and pressure profiles near their surfaces. The pressure dropping to very small values thus indicates that the surface has been reached.

Besides the EOS, temperature, and Y_e prescription, the central baryon density ρ_{b, c} is the only other free parameter. Since we are interested in the maximum mass that can be supported, we compute sequences with varying ρ_{b, c} for each EOS, but limit ourselves to ρ_nuc < ρ_{b, c} ⩽ ρ_{max, EOS}, where the latter is just the maximum density entry in the respective EOS table. HMNSs with central densities below ρ_nuc are not realistic (see Sekiguchi et al. 2011).

We make our TOV solver, all P = P(ρ_b) tables, and the Python scripts used to create the results in this paper available on http://www.stellarcollapse.org.

We generate axisymmetric equilibrium models using the code originally presented in (Cook et al. 1992, hereafter CST; see also Cook et al. 1994a, 1994b), which is based on the relativistic self-consistent field method of Komatsu et al. (1989a). The axisymmetric equilibrium equations are solved iteratively on a grid in (s, μ), where s is a compactified radial coordinate and μ = cos θ, where θ is the usual spherical polar angle. Additionally, metric functions are solved using Green's functions integrals expanded in terms of N_l Legendre polynomials. Consequently, the total numerical resolution is specified via a tuple of (N_s, N_μ, N_l), which we set to (500, 300, 16). The resolution is chosen so that the resulting integral quantities of the equilibrium solution (e.g., its gravitational mass) are precise to about one part in 10³. The surface of the star is defined by an enthalpy contour which is specified in the code by setting a surface energy density. This energy density has a default value of 7.9 g cm⁻³, and we have checked that increasing its value by a factor of 10⁶ leaves the physical quantities of the solution unchanged to our stated general error level of 10⁻³.

An axisymmetric HMNS equilibrium configuration is constructed by the CST code based on choices of (1) a barotropic EOS, (2) a rotation law, (3) the rotation rate, and, (4) the maximum mass-energy density E_max = [ρ_b(1 + /c²)]_max of the configuration.

In order to keep the size of the parameter space manageable, we restrict rotating configurations to the LS220 and HShen EOS and set up barotropic versions using the temperature and composition parameterizations described in Section 2.2. Since the EOS obtained with ν-full and ν-less β-equilibrium differ only very mildly (see Figure 3), we construct rotating configurations under the simple assumption of ν-less β-equilibrium.

We employ the "j − const." rotation law (see, e.g., CST), which is commonly used in the literature for HMNS models (e.g., Baumgarte et al. 2000). The degree of differential rotation is parameterized by $\tilde{A}$.⁸ In the Newtonian limit, this rotation law becomes $\Omega = \Omega _c / (1 + \tilde{A}^2 \varpi ^2/r_{\rm e}^2)$, where ϖ is the cylindrical radius, r_e is the radius of the star at its equator, and Ω_c is the central angular velocity. For $\tilde{A} = 0$, one recovers uniform rotation, while for large $\tilde{A}$, the specific angular momentum becomes constant (i.e., Ω∝ϖ⁻² in the Newtonian limit). We explore values of $\tilde{A}$ between 0 and 1. The latter value of $\tilde{A}$ corresponds to roughly a factor of two decrease of the angular velocity from the center to the HMNS surface, which is in the ball park of what is found in merger simulations (e.g., Shibata et al. 2005). Once the rotation law is fixed, the rotation rate is determined by specifying the axis ratior_p/e, defined as the ratio of the HMNS radius along the pole r_p divided by the radius at the equator r_e.

The final parameter to be chosen is the maximum energy density of the configuration. For simplicity and consistency with the choice of variables for the TOV solutions discussed in Section 2.3, we set E_max by choosing a maximum baryon density ρ_{b, max} and obtain E(ρ_{b, max}) from the EOS.

For each choice of EOS, $\rho _\mathrm{b,\mathrm{max}}$, and $\tilde{A}$, we compute a sequence of models with increasing rotation rate, stepping down from r_p/e = 1 (the nonrotating TOV case) until we reach mass shedding or until the code fails to converge to an equilibrium solution. In the case of uniform rotation ($\tilde{A} = 0$) the sequence always ends at mass shedding, the resulting rotating NS has spheroidal shape, and the maximum and central density coincide (ρ_{b, max} = ρ_c). Differentially rotating sequences, on the other hand, can bifurcate into two branches: one with ρ_{b, max} = ρ_c and spheroidal geometry and one with an off-center location of ρ_{b, max} and quasitoroidal shape. For differentially rotating models, the CST solver generally fails to converge to a solution at r_p/e before mass shedding and, therefore, possibly before the maximum mass for a given configuration is reached. This limitation means that the maximum masses we state for differentially rotating models are to be interpreted as lower bounds on the true maximum masses. The code developed by Ansorg et al. (2003) is far more robust than CST for such extreme configurations and these authors have argued that with increasing degree of differential rotation, arbitrarily large masses could be supported in extremely extended tori, but such configurations are unlikely to be astrophysically relevant.

It has been widely recognized that uniform rotation can support a supramassive NS against gravitational collapse (see, e.g., Friedman et al. 1986; Friedman & Ipser 1987). A supramassive NS is defined as a stable NS with a mass greater than the maximum mass of a TOV star with the same EOS (CST). At a given central density, the mass that may be supported rises with increasing angular velocity until the material on the NS's equator becomes unbound (the mass-shedding limit). This leads to the supramassive limit, a well defined maximum mass for uniformly rotating NSs with a specified EOS.

In Figure 5, we plot the baryonic mass M_b as a function of maximum baryon density for TOV and uniformly rotating mass-shedding sequences obtained with the LS220 EOS (left panel) and the HShen EOS (right panel). Focusing first on the TOV sequences, one notes that at low central densities (ρ_b ≲ few × ρ_nuc), M_b is significantly increased by thermal effects. This is because the mean density $\bar{\rho }_{{\rm _b}}$ of such configurations is in the regime in which thermal pressure is of greatest relevance (see Figure 3) and can alter the structure of the bulk of the NS. This carries over to the uniformly rotating case. The extended hot configurations reach mass shedding at lower angular velocities than their cold counterparts, but the extended, low $\bar{\rho }_{{\rm b}}$ cores of hot configurations receive sufficient rotational support to yield a higher M_b. This, however, is the case only for centrally hot cXp0 configurations. Models with hot envelopes (with parameterizations c30p5 and c30p10) benefit less from rotational support.

Zoom In Zoom Out Reset image size

With increasing maximum density, the baryonic masses of the TOV models for different temperature parameterizations converge for a given EOS. Near the density at which the maximum mass is reached, the increase in M_b in hot configurations has turned into a slight decrease for models computed with the LS220 EOS and has dropped to ≲5% for the HShen EOS (see also Figure 4). The mass-shedding sequences show a more complex behavior with increasing maximum density. As in the TOV case, the mean density $\bar{\rho }_b$ of the NSs increases and less material is experiencing enhanced pressure support due to high temperatures in the cXp0 models. Hence, these models move toward the $M_\mathrm{b}^\mathrm{max}$ of the cold supramassive limit (see the inset plots in Figure 5). For both EOS, the $M_\mathrm{b}^\mathrm{max}$ of hot configurations are all lower than the cold value. The cXp0 models reach supramassive limits that are within less than 2% of the cold supramassive limit for both EOS. The c30p10 and c30p5 models, on the other hand, have $M_\mathrm{b}^\mathrm{max}$ that are ∼5%–10% lower than the cold supramassive limit for both EOS. Table 2 summarizes key parameters of the hot and cold configurations at the supramassive limit.

Table 2. Uniformly Rotating Neutron Stars at the Supramassive Limit

Model	ρb, max	$M_\mathrm{b}^\mathrm{max}$	$M_\mathrm{g}^\mathrm{max}$	re	rp/e	Ω	T/|W|
	(1015 g cm−3)	(M☉)	(M☉)	(km)		(103 rad s−1)
LS220 cold	1.653	2.823	2.419	14.429	0.566	10.096	0.118
LS220 c20p0	1.652	2.760	2.384	14.788	0.574	9.647	0.106
LS220 c30p0	1.652	2.737	2.382	15.000	0.576	9.441	0.103
LS220 c30p5	1.710	2.671	2.322	15.300	0.587	9.031	0.088
LS220 c30p10	1.769	2.587	2.247	16.130	0.599	8.215	0.066
LS220 c40p0	1.625	2.717	2.383	15.201	0.577	9.262	0.101
HShen cold	1.220	3.046	2.649	17.101	0.564	8.233	0.117
HShen c20p0	1.196	3.006	2.629	17.760	0.573	7.745	0.105
HShen c30p0	1.171	3.009	2.648	18.173	0.574	7.511	0.103
HShen c30p5	1.228	2.916	2.564	18.665	0.588	7.086	0.084
HShen c30p10	1.261	2.808	2.467	20.070	0.604	6.238	0.060
HShen c40p0	1.139	3.012	2.664	18.474	0.574	7.355	0.101

Notes. Summary of mass-shedding uniformly rotating supramassive neutron star configurations at the maximum mass for each EOS and temperature prescription. These models are in ν-less β-equilibrium (see Section 2.2). ρ_{b, max} is the central density of the model with the maximum baryonic mass $M_\mathrm{b}^\mathrm{max}$. $M_\mathrm{g}^\mathrm{max}$ is the gravitational mass at the ρ_{b, max} at which $M_\mathrm{b}^\mathrm{max}$ occurs. r_e is the equatorial radius, r_p/e is the axis ratio, Ω is the angular velocity, and T/|W| is the ratio of rotating kinetic energy T to gravitational energy |W|.

Download table as:  ASCIITypeset image

The systematics of the supramassive limit with temperature prescription becomes clear when considering Figure 6. This figure shows the baryonic mass M_b and gravitational mass M_g for uniformly rotating NSs as a function of angular velocity Ω for the LS220 and HShen EOS at fixed densities near the maximum of M_b(ρ_{b, max}) (see Table 2). At fixed angular velocity below mass shedding, hotter configurations always yield higher M_g than their colder counterparts. For the LS220 EOS, as in the TOV case discussed in the previous section 3, hotter configurations have lower M_b. In the case of the HShen EOS, which generally yields less compact equilibrium models, the opposite is true, but the increase in M_b caused by thermal support is smaller than the increase in M_g.

Zoom In Zoom Out Reset image size

With increasing Ω, the mass-shedding limit is approached and hotter configurations systematically reach the mass shedding limit at lower angular velocities. The reason for this is best illustrated by comparing c30p0 models with c30p10 and c30p5 models, which have a high-temperature plateau at low densities of 10 MeV and 5 MeV, respectively. At low angular velocities, all c30pX models show the same thermal increase in M_g. However, the high pressure at low densities in the c30p10 and c30p5 models leads to significantly larger radii compared to the model without temperature plateau. Consequently, as Ω is increased, the configurations with plateau reach the mass-shedding limit at lower Ω. For the LS220 EOS, the c30p10 sequence terminates at ∼8200 rad s⁻¹, the c30p5 sequence terminates at ∼9200 rad s⁻¹, while the c30p0 sequence does not terminate before ∼9800 rad s⁻¹. The HShen model sequences show the same qualitative trends.

Differential rotation can provide centrifugal support at small radii while allowing a NS configuration to stay below the mass-shedding limit at its equatorial surface. Differentially rotating equilibrium configurations have been shown to support masses well in excess of the supramassive limit (e.g., Ostriker et al. 1966; Baumgarte et al. 2000; Morrison et al. 2004). Such configurations are referred to as "hypermassive." However, since there is (mathematically speaking) an infinite number of possible differential rotation laws, it is impossible to define a formal "hypermassive limit" for the maximum mass of HMNSs in the way it is possible for uniformly rotating supramassive NSs. Nevertheless, we can study the systematics of the supported baryonic (and gravitational) masses with variations in the HMNS temperature profile, maximum baryon density, and degree and rate of differential rotation for the rotation law considered in this study, which is not drastically different from what is found in merger simulations (e.g., Shibata et al. 2005).

In Figure 7, we show the supported baryonic mass M_b as a function of maximum baryon density ρ_{b, max} for cold, c20p0, and c40p0 temperature prescriptions, both EOS, and for different choices of $\tilde{A}$ (see Table 3 for a summary of quantitative results). The curves represent configurations with the minimum r_p/e at which an equilibrium solution is found by the CST solver (i.e., the most rapidly spinning setup). Note that the peaks of these curves represent only lower limits on the maximum HMNS mass. In addition, we plot only solutions with ratios T/|W| of rotational kinetic energy T to gravitational energy |W| below 25%, since more rapidly spinning models would be dynamically nonaxisymmetrically unstable (Chandrasekhar 1969; Baiotti et al. 2007). It is this limit which defines the rising branch of the M_b(ρ_{b, max}) curve at the lowest densities in Figure 7 for $\tilde{A} = 1.0$. Note that many of these configurations may still be unstable to secular rotational instabilities or rotational shear instabilities (e.g., Watts et al. 2005; Ott et al. 2007; Corvino et al. 2010).

Zoom In Zoom Out Reset image size

Table 3. Differentially Rotating Hypermassive Neutron Stars

Model	ρb, max	$M_\mathrm{b}^\mathrm{max}$	$M_\mathrm{g}^\mathrm{max}$	re	rp/e	$\tilde{A}$	Ωc	T/|W|
	(1015 g cm−3)	(M☉)	(M☉)	(km)			(103 rad s−1)
LS220 cold	0.993	3.648	3.140	17.258	0.376	0.5	15.476	0.244
LS220 c20p0	0.852	3.573	3.124	18.538	0.364	0.6	15.047	0.243
LS220 c30p0	0.706	3.568	3.167	19.611	0.344	0.7	14.888	0.249
LS220 c30p5	0.600	3.413	3.064	21.870	0.320	0.9	14.461	0.250
LS220 c30p10	0.990	3.090	2.723	19.208	0.421	0.9	16.330	0.187
LS220 c40p0	0.692	3.597	3.211	19.931	0.344	0.7	14.677	0.249
HShen cold	0.766	4.101	3.562	19.800	0.372	0.5	13.450	0.245
HShen c20p0	0.641	4.076	3.585	21.352	0.360	0.6	13.042	0.245
HShen c30p0	0.532	4.099	3.650	22.305	0.344	0.7	13.131	0.249
HShen c30p5	0.517	3.942	3.527	24.371	0.340	0.8	12.426	0.243
HShen c30p10	0.646	3.529	3.141	23.521	0.400	1.0	13.934	0.196
HShen c40p0	0.514	4.148	3.708	22.701	0.344	0.7	12.888	0.249

Notes. Summary of the differentially rotating HMNS configurations with the largest baryonic masses for each EOS and temperature prescription. These configurations are obtained in a sequence from $\tilde{A} = 0$ to $\tilde{A} = 1$ with spacing $\delta \tilde{A} = 0.1$ and are to be seen as lower bounds on the maximum achievable masses. The sequences considered here exclude dynamically nonaxisymmetrically unstable models with ratio of rotational kinetic energy to gravitational energy T/|W| > 0.25. The quantities listed in the table are the following: ρ_{b, max} is the baryon density at which the maximum baryonic mass $M_\mathrm{b}^\mathrm{max}$ occurs, $M_\mathrm{g}^\mathrm{max}$ is the gravitational mass at that density, r_e is the equatorial radius of the configuration, r_p/e is its axis ratio, $\tilde{A}$ is the differential rotation parameter at which $M_\mathrm{b}^\mathrm{max}$ obtains. Ω_c is the central angular velocity of the configuration and T/|W| is its ratio of rotational kinetic energy to gravitational energy. We note that the accuracy of the results listed in this table is set by the step size in r_p/e, which we set to δr_p/e = 0.004.

Download table as:  ASCIITypeset image

The overall shape of the M_b(ρ_{b, max}) curves in Figure 7 is qualitatively similar to what is shown in Figure 1 of Baumgarte et al. (2000) for Γ = 2 polytropes and Figure 2 of Morrison et al. (2004) for the cold Friedman–Pandharipande EOS (Friedman & Pandharipande 1981). The LS220 and HShen EOS yield qualitatively very similar results, but the supported HMNS masses found by the CST solver are, as expected, systematically higher for models with the HShen EOS than for those using the LS220 EOS. One notes, however, interesting variations with temperature prescription. At low ρ_{b, max}, thermal pressure support leads to increased M_b and more differentially rotating configurations have higher M_b. Sequences with $\tilde{A} \lesssim 0.5$ show similar systematics with density and temperature prescription as the uniformly spinning ones discussed in Section 4.1: As the density increases, hot configurations converge toward the cold sequence and reach their maximum M_b near and below the maximum of the cold sequence. Models with $\tilde{A} \gtrsim 0.5$, on the other hand, have more steeply rising curves with ρ_{b, max} and are discontinuous (i.e., exhibit a "kink") at their global maxima. At these points quasitoroidal solutions appear. Furthermore, the slope of the curve describing (as a function of ρ_{b, max}) the axis ratios r_p/e at which the solver stops converging discontinuously changes sign. We attribute this behavior, which was also observed by Morrison et al. 2004, to a bifurcation of the sequence between models, which continue shrinking in axis ratio until they become completely toroidal (r_p/e = 0), and less extreme models that stay quasitoroidal or spheroidal. Beyond the "kink" in $\tilde{A}\gtrsim 0.5$ sequences, thermal effects play little role.

The lower bounds of the range of ρ_{b, max} shown in the two panels of Figure 7 (and also Figure 8) are chosen for the following reason: fully dynamical merger simulations by, e.g., Sekiguchi et al. (2011), Baiotti et al. (2008), Shibata et al. (2005), Kiuchi et al. (2009), Bauswein et al. (2012), and Thierfelder et al. (2011), all suggest a rule of thumb that the postmerger maximum baryon density of the HMNS is typically not less than ∼80% of the central density of the progenitor NSs. We can derive a rather solid EOS-dependent lower limit on ρ_{b, max} for HMNS remnants from (equal mass) NSNS mergers in the following way. In order to form an HMNS, constituent equal-mass NSs must at the very least have a mass that is 50% of the maximum mass in the cold TOV limit. Hence, the premerger central density must at least be that of a TOV solution with $M_\mathrm{b} = 0.5 M_\mathrm{b}^\mathrm{max,TOV}$. Using the aforementioned empirical result from merger simulations, we arrive at

$$\begin{equation} \rho _\mathrm{b,min} = 0.8 \rho _\mathrm{b,TOV}(M_\mathrm{b} = M_\mathrm{b,max}/2)\,\,. \end{equation} \tag{ 7 }$$

For the LS220 EOS, ρ_{b, TOV}(M_b = M_{b, max}/2) ∼ 5.8 × 10¹⁴ g cm⁻³ and occurs at M_b (M_g) of 1.19 M_☉ (1.10 M_☉). For the HShen EOS, ρ_{b, TOV}(M_b = M_{b, max}/2) ∼ 4.4 × 10¹⁴ g cm⁻³ and occurs at M_b (M_g) of 1.28 M_☉ (1.20 M_☉). Applying the density cut given by Equation (7) excludes most dynamically nonaxisymmetrically unstable configurations.

Zoom In Zoom Out Reset image size

Figure 8, like Figure 7, shows baryonic mass as a function of maximum baryon density for both EOS and a variety of $\tilde{A}$, but contrasts models c30p5 and c30p10, which have hot plateaus at low densities, with cold models. The qualitative features discussed in the following are identical for both EOS. In the TOV case and at low densities, M_b is enhanced primarily by the hot core, since nonrotating solutions are compact and dominated by ρ_b ≳ 10¹⁴ g cm⁻³, where the high-temperature plateaus play no role. At higher densities, the M_b curves of hot models converge to near or below the cold TOV maximum M_b. The situation is different for uniformly and moderately differentially rotating models ($\tilde{A} \lesssim 0.5$). Rotation shifts these configurations to lower mean densities and the hot plateaus lead to equatorially bloated solutions. These reach their minimum r_p/e for which a solution can be found at lower angular velocities. Hence, centrifugal support is weaker and the configuration with the hottest plateau has the lowest M_{b, max}. The behavior is different at high degrees of differential rotation ($\tilde{A} = 1$). The cold and the c30p5 models are HMNSs and quasitoroidal already at the lowest densities shown in Figure 8. The c30p5 sequence has slightly larger M_b than the cold sequence. The c30p10 sequence, however, is spheroidal at low ρ_{b, max} and then discontinuously transitions to the quasitoroidal branch, which is marked by a large jump in M_b.

In order to illustrate this discontinuous behavior further, we plot in Figure 9 the equatorial radius of equilibrium solutions as a function of central angular velocity at $\tilde{A} = 1$ and for three different fixed ρ_{b, max}. We show curves obtained with the LS220 EOS for the cold, c30p5, and c30p10 temperature prescriptions. The curves are parameterized by decreasing r_p/e and terminate at the smallest value at which the solver converges. The three densities are chosen so that the first two are below and the third is above the jump of the c30p10 curve in Figure 8. At all ρ_{b, max}, the hot configurations have significantly larger radii than the cold models, but decreasing r_p/e leads to increasing Ω_c and only modest radius changes for cold and c30p5 models. This is very different for the c30p10 sequence. At ρ_{b, max} = 7.11 × 10¹⁴ g cm⁻³ these models do not become quasitoroidal and the r_e − Ω_c mapping becomes double-valued as the decrease in r_p/e turns from a decrease of r_p at nearly fixed r_e and increasing Ω_c into a steep increase of r_e and a decrease of Ω_c. As ρ_{b, max} increases, less material is at low densities where thermal pressure support is strong in the c30p10 models. Consequently, the solutions are more compact and stay so to smaller r_p/e. ρ_{b, max} = 8.16 × 10¹⁴ g cm⁻³ is the critical density at which the very last point in the sequence of decreasing r_p/e (the one shown in Figure 8) jumps discontinuously to large r_e. At ρ_{b, max} = 9.21 × 10¹⁴ g cm⁻³, which is above the critical density for c30p10 in Figure 8, the c30p10 models become quasitoroidal as r_p/e decreases and Ω_c increases. They exhibit the same systematics as the c30p5 and cold models. We note that what we have described for the c30p10 models also occurs for the c30p5 models, although at significantly lower densities ρ_{b, max} ≲ 5 × 10¹⁴ g cm⁻³ and even the cold models show similar trends at low densities.

Zoom In Zoom Out Reset image size

The sequences shown in Figures 7 and 8 are extreme configurations in the sense that models with smaller r_p/e cannot be found by the CST solver and may not exist for the rotation law that we consider here. Real HMNS may not by such critical rotators. In Figure 10, we plot M_b and M_g for the LS220 EOS as a function of central angular velocity Ω_c and temperature prescription. We fix the degree of differential rotation to $\tilde{A}=1$ and show sequences in Ω_c for a fixed maximum density ρ_{b, max} = 9.21 × 10¹⁴ g cm⁻³, which is the highest density shown in Figure 9. The transition to quasitoroidal shape is smooth and quasitoroidal configurations are marked with symbols. The end points of the M_b curves shown in Figure 10 and in Figure 9 correspond to the M_b values of the $\tilde{A}=1$ curves in Figures 7 and 8 at 9.21 × 10¹⁴ g cm⁻³.

Zoom In Zoom Out Reset image size

Figure 10 shows that, as in the case of uniform rotation (see Figure 6), hotter subcritically differentially spinning configurations have higher M_g. At the density chosen for this plot, they also have higher M_b, but at the higher densities at which the masses of uniformly spinning models peak, the M_b of hotter configurations are smaller than those of colder ones. It is particularly remarkable that the models with the hot plateau at low densities show the greatest thermal enhancement. They also transition to a quasitoroidal shape last but terminate the earliest in Ω_c. Nevertheless, for the ρ_{b, max} chosen here, they can support slightly more mass at critical rotation than their counterparts without low-density temperature plateau.

The existence of a maximum mass for equilibrium sequences of nonrotating (TOV) NSs is one of the most important astrophysical consequences of general relativity and, hence, is well known in the study of compact objects. The parameter space of hot differentially rotating HMNS models studied here is vast and complex. In the following, we briefly review the classical results on the stability of stationary NSs and formulate how one may reason regarding the stability of HMNS equilibrium models.

A particular useful approach to the stability problem is the turning-point method of Sorkin (1982). The turning-point method allows one to reason about the stability of sequences of equilibrium solutions solely by examining the parameter space of equilibrium models without dynamical simulations or linear perturbation analysis. The turning-point method has been used extensively in previous work on the stability of cold and uniformly rotating NSs (e.g., CST; Friedman et al. 1988; Stergioulas & Friedman 1995; Read et al. 2009).

An equilibrium sequence is a one-dimensional slice from the space of equilibrium models indexed by some parameter. Here we use ρ_{b, max} as our sequence parameter. A model in the space of equilibrium models may be defined by the following conserved quantities: the gravitational mass M_g, baryonic mass M_b, total angular momentum J, and total entropy S. Generally, as one changes the sequence parameter, ρ_{b, max}, the quantities (M_g, M_b, J, S) will vary. A turning point in the sequence occurs when 3 out of 4 of the derivatives d/dρ_{b, max} of (M_g, M_b, J, S) vanish. For this point in ρ_{b, max}, the turning point theorem shows (1) that the derivative of the fourth quantity in the tuple also vanishes, and (2) that the sequence must have transitioned from stable to unstable (Sorkin 1982; Kaplan 2014). This characterization of the space of equilibrium models relies on the assumption that the change in M_g depends to first order only on the total changes in baryonic mass M_b, angular momentum J, and entropy S, and not on changes to their higher moments. That is, changes in the distribution of entropy, baryonic mass and angular momentum. In nature, this will generally not be the case, since cooling and angular momentum redistribution will change the entropy and angular momentum distributions, respectively. However, these changes will be slow and not drastic so that changes to the total energy due to changes in these higher order moments will be small. We account for such changes approximately by considering different degrees of differential rotation and a range of temperature prescriptions in the following.

If we are considering the special case of zero-temperature configurations, then the entropy S is no longer relevant to the equilibrium's stability, since the change to the configuration's energy due to a change in entropy is also zero. In this case, a turning point may be identified when two out of three of the set d/dρ_{b, max}(M_g, M_b, J) are zero. Zero temperature is a very good approximation for our cold equilibrium models. In Figure 11, we plot M_g along constant M_b sequences with M_b = 2.9 M_☉ for the HShen EOS (M_b = 2.9 M_☉ corresponds to M_b of an HMNS formed from two NSs of M_g = 1.35 M_☉, assuming no mass loss). All of these curves have a minimum located at ρ_{b, max} ≳ 1 × 10¹⁵g cm⁻³.

Zoom In Zoom Out Reset image size

For the cold sequences, these minima are turning points because dM_g/dρ_{b, max} and dM_b/dρ_{b, max} are both zero. Any models along those curves at densities in excess of ρ_{b, max} at the minima are secularly unstable to collapse. For the hot temperature parameterizations,¹⁰ the minima are only approximations to the turning point (which we shall call approximate turning points) because only two out of four (dM_g/dρ_{b, max} and dM_b/dρ_{b, max}) of the derivatives of (M_g, M_b, J, S) are zero. We argue that these approximate turning points are good indicators of the onset of instability for the equilibrium sequences for several reasons. (1) We find that the approximate turning points for all considered temperature parameterizations and measures of differential rotation ($\tilde{A} = 0$ to $\tilde{A} = 1.1$ with spacing $\delta \tilde{A} = 0.1$) lie within the same ∼25% range in ρ_{b, max} indicated by the blue lines in Figure 11 (similarly within a ∼25% range in ρ_{b, max} for the LS220 EOS). (2) In cold uniformly rotating NS models, approximate turning points occur where one out of three of d/dρ_{b, max}(M_g, M_b, J) vanish. The study of such models shows that the actual turning point density is within only ∼1% of the approximate turning point density (where dM_g/dρ_{b, max} = 0 along the mass-shed sequence; see Figure 10 of Stergioulas & Friedman 1995). (3) The turning-point condition is a sufficient, but not necessary, criterion for secular instability. Thus instability must set in at ρ_{b, max} greater than the turning-point ρ_{b, max}, but may set in already at lower densities (see, e.g., Takami et al. 2011 for an example). It is thus conservative to use the approximate turning point located at the highest ρ_{b, max} over all sequences for a given EOS as an upper bound for the maximum stable ρ_{b, max} of HMNS models for that EOS.

Further to the above, we have verified (see Section 9.1 of Kaplan 2014) that the same density ranges contain approximate turning points when examining alternate pairs of conserved variables: both J and M_b, and J and M_g (in contrast to Figure 11, where we examine M_b and M_g). This gives us confidence that the method of approximate turning points is self consistent with respect to choice of the vanishing derivatives. Unfortunately, since the CST code employs only barytropic EOSs, we lack the infrastructure necessary to study the total entropy of the configurations, and note that an examination of total entropy of these models is an important goal for future work.

An HMNS remnant resulting from the merger of two NSs that does not promptly collapse into a black hole will settle into a quasiequilibrium state. More precisely, this is a state in which the HMNS is no longer in dynamical evolution, measured, for example, by oscillations in the HMNS maximum density. This should occur several dynamical times after merger. From this point on, the HMNS will evolve secularly along some sequence of equilibrium models. A secular evolution is, by definition, a dissipative process that may involve energy loss¹¹ from the system. Consequently, we may parameterize the secular evolution of the HMNS toward a turning point via the change in its total mass-energy, which, in our case, is the change in gravitational mass of the equilibrium model. This occurs in HMNSs via neutrino cooling and the emission of gravitational radiation. In addition, the rotational energy of the HMNS may be reduced by angular momentum redistribution via the MRI, provided this occurs sufficiently slowly to be characterized as a secular process. This can lead to a build up of magnetic field, or dissipation of the free energy of differential rotation as heat (see, e.g., Thompson et al. 2005 for a detailed discussion), which may lead to increased neutrino cooling. Furthermore, specific angular momentum transported to the HMNS surface may unbind surface material, leading to a decrease in J and M_b. These changes of M_b and J may be significant, but cannot be taken into account by the approximate description of the HMNS's evolution we are considering here. Our results should thus be interpreted with these limitations in mind.

A secularly evolving HMNS will, in general, evolve in the direction of decreasing gravitational mass M_g while (at least approximately) conserving its total baryonic mass M_b. This results in an increasing density and compactness of the star. Figure 11 shows, for a fixed temperature prescription and differential rotation parameter, that the gravitational mass M_g of a sequence with fixed baryonic mass M_b = 2.9 M_☉ (using the HShen EOS; we find qualitatively the same for the LS220) is decreasing with increasing density. This continues until, M_g reaches a minimum at an approximate turning point for ρ_{b, max} ≳ 1 × 10¹⁵ g cm⁻³. Here, δM_g = 0, and δM_b vanishes by our choice of a constant M_b sequence.

The curves in Figure 11 are shown for constant differential rotation parameter $\tilde{A}$. However, an HMNS of M_b = 2.9 M_☉ is not necessarily constrained to a specific curve. One would expect the HMNS to evolve to neighboring curves of less extreme differential rotation (decreasing $\tilde{A}$), in accordance with its loss of angular momentum due to gravitational waves and its redistribution of angular momentum due to other secular processes. Nevertheless, consider the limit in which the HMNS is constrained to a curve of constant $\tilde{A}$. Then it would evolve secularly until reaching the curve's minimum. At this point, any further energy loss implies that the HMNS must either (1) secularly evolve to a nearby equilibrium sequence with lower temperature or lower degree of differential rotation and higher density (another curve on the plot) or (2) undergo collapse to a black hole. Note that the densities at which the minimum occurs for different $\tilde{A}$ and temperatures are remarkably close to each other. For the sequences using the HShen EOS shown in Figure 11, the approximate turning points lie in the range 1.05 × 10¹⁵ g cm⁻³ < ρ_{b, max} < 1.30 × 10¹⁵ g cm⁻³ for all considered $\tilde{A}$ and both shown temperature prescriptions. The constant-M_b curves for other temperature parameterizations (c20p0, c30p0, c30p5, c30p10) are all located in-between the curves for the c40p0 and cold cases shown. Thus, we expect that the point of collapse for an HMNS will be marked by its evolution to this density regime regardless of the temperature distribution of the model.

From the above findings, we conclude that thermal effects have little influence on the stability of HMNSs in rotational equilibrium against gravitational collapse. However, our results do imply that thermal support will affect at what density the HMNS first settles to its quasiequilibrium state. The discussion in Section 4.2 and, in particular, Figure 10, illustrates that at subcritical rotation rates and densities significantly below those of the approximate turning points, models with hot temperature profiles have a larger M_b compared to models with cooler temperatures at the same ρ_{b, max}. Thus an HMNS with greater thermal support will reach a quasiequilibrium at a lower ρ_{b, max}, and thus have more energy to lose before it can evolve to the critical density regime for collapse.

While thermal effects may be important in setting the initial conditions for the secular evolution of an HMNS, they appear to be of little consequence to the stability of an HMNS in quasiequilibrium. Once in a quasiequilibrium state, the energy lost by an HMNS during its secular evolution is the most robust indicator for its progress toward instability and collapse. Figure 11 shows that this is true regardless of the degree of differential rotation of the HMNS. For a fixed temperature parameterization, the difference in M_g between different degrees of differential rotation is at most ∼0.005 M_☉, corresponding to ≲ 10% of the total energy lost during the HMNS's secular evolution.

Sekiguchi et al. (2011) conducted simulations of NSNS mergers using the HShen EOS and included neutrino cooling via an approximate leakage scheme. They considered three equal-mass binaries with component NS gravitational (baryonic) masses of 1.35 M_☉ (1.45 M_☉), 1.50 M_☉ (1.64 M_☉), 1.60 M_☉ (1.77 M_☉) denoted as L, M, and H, respectively. The HMNS formed from their high-mass binary collapses to a black hole within ≲ 9 ms of merger. The low-mass and the intermediate-mass binaries, however, form hot (T ∼ 5–30 MeV) spheroidal quasiequilibrium HMNSs that remain stable for at least 25 ms, the duration of their postmerger simulations.

Sekiguchi et al. (2011) argue that thermal pressure support could increase the maximum mass of HMNSs with T ≳ 20 MeV by 20%–30%. The results that we lay out in Sections 3 and 4 of our study suggest that it is not straightforward to disentangle centrifugal and thermal effects for differentially rotating HMNS. Our findings show that critically spinning configurations (i.e., configurations at which the maximum M_b is obtained for a given $\tilde{A}$) of hot models do not lead to an increase in the maximum supported baryonic mass by more than a few percent and in most cases predict a lower maximum mass than in the cold case. We find it more useful to consider the results of Sekiguchi et al. (2011) in the context of the evolutionary scenario outlined in Section 5.2.

In Figure 12, we plot M_b as a function of ρ_{b, max} for select sequences of uniformly and differentially rotating models obtained with the HShen EOS with the cold and c40p0 temperature prescriptions. We also mark the immediate postmerger densities of the L, M, and H models of Sekiguchi et al. (2011) and their evolutionary tracks (in ρ_{b, max}). The high-mass model H never settles into a quasiequilibrium and collapses to a black hole during the dynamical early postmerger phase. Its ρ_{b, max} evolves within ∼9 ms from 0.58 × 10¹⁵ g cm⁻³ to values beyond the range of the plot. Our secular-evolution approach cannot be applied to this model since it never reaches a quasiequilibrium state. The lower-mass M and L models enter Figure 12 at successively lower densities. Their "ring-down" oscillations are damped by ∼9 ms after which the HMNSs evolve secularly with ρ_{b, max} that increase roughly at the same rate in both models, suggesting that their rate of energy loss is comparable. At such early times, gravitational waves are most likely dominating energy loss (see the discussion of timescales in Paschalidis et al. 2012), and, indeed, model M and L exhibit similar gravitational wave amplitudes and frequencies (Sekiguchi et al. 2011, Figure 4). Focusing on model L, we now consider Figure 11, which shows sequences of constant M_b (for model L with M_b ∼ 2.9 M_☉). As the HMNS loses energy, M_g decreases and the HMNS evolves to the right (toward higher ρ_{b, max}). Model L enters its secular evolution at a central density of ∼0.56 × 10¹⁵ g cm⁻³ and evolves secularly to ∼0.68 × 10¹⁵ g cm⁻³ within ∼16 ms. Largely independent of its specific angular momentum distribution and thermal structure, Figure 11 suggests that this model will reach its global minimum M_g and, thus, instability in a small density range of ∼1.05–1.30 × 10¹⁵ g cm⁻³.

Zoom In Zoom Out Reset image size

Using our approximate secular evolution model for HMNSs discussed in Section 5.2, we linearly extrapolate the density evolution of model L in Sekiguchi et al. (2011). We expect a possible onset of collapse at t ≳ 58 ms after merger (and ≳49 ms after the start of the secular evolution). These numbers should be regarded as very rough estimates, given the limitations and rather qualitative nature of our model. Depending on its angular momentum when entering its secular evolution, its cooling rate, angular momentum redistribution and loss, model L may alternatively evolve into a long-term stable supramassive NS, since a baryonic mass of ∼2.9 M_☉ can in principle be supported by the HShen EOS at the supramassive limit (see Table 2). Furthermore, we have also checked that model L of Sekiguchi et al. (2011) contains sufficient angular momentum to be represented by the sequences identified in Figures 11 and 12. At a time of ∼10–15 ms after merger, model L has an angular momentum of 6 × 10⁴⁹ g cm² s⁻¹(6.8 in $c = G = M_{_\odot }$ units). Plots of similar sequences can be found in Kaplan (2014). They are consistent with this value.

The role of thermal pressure effects in all of the above is relatively minor (see the very similar ρ_{b, max} locations of the M_g minima in hot and cold configurations shown in Figure 11). However, when first entering the secular regime as a subcritical HMNS, a configuration with higher temperature and stronger thermal pressure support will be less compact and will have a lower ρ_{b, max} at a fixed M_b than a colder one. Hence, in the picture of secular HMNS evolution discussed in Section 5.2, such a configuration would have to evolve "farther" in ρ_{b, max} to reach criticality and, thus, can survive longer at fixed energy loss rates.

Paschalidis et al. (2012) performed NSNS merger simulations of Γ = 2 polytropes in which they approximated a thermal pressure component with a Γ = 2 Γ-law. Their postmerger HMNS enters its secular evolution in a quasitoroidal configuration with two high-density, low-entropy cores, a central, lower-density, hot region and a high-entropy low-density envelope. The total mass of their model can be arbitrarily rescaled, but in order to estimate temperatures and thermal pressure contributions, the authors scaled their HMNS remnant to a gravitational mass of 2.69 M_☉. With this, they estimated in their quasitoroidal HMNS peak and rms temperatures of ∼20 MeV and ∼5 MeV, respectively. Paschalidis et al. (2012) studied the effect of neutrino cooling on the HMNS evolution by introducing an ad-hoc cooling function that removes energy proportional to the thermal internal energy (neglecting the stiff temperature dependence of neutrino cooling). In order to capture effects of cooling during the limited simulated physical postmerger time, they drained energy from their HMNS at rates ∼100–200 times higher than realistic cooling by neutrinos.

The authors considered cases without cooling and with two different accelerated cooling timescales. As cooling is turned on in their simulations, the slope of the maximum baryon density ρ_{b, max}(t) of the HMNS increases discontinuously and the higher the cooling rate, the faster the evolution to higher ρ_{b, max}(t). The HMNSs in both cases with cooling become unstable at different times, but roughly at the same ρ_{b, max}. This is consistent with the secular HMNS evolution picture laid out in Section 5.2. Cooling reduces the total energy of the system (M_g) and drives the HMNS to higher ρ_{b, max} at fixed M_b until the (approximate) turning point is reached and collapse ensues. However, losses due to gravitational wave emission and angular momentum redistribution and shedding will have the same effect and may dominate in nature, since they are likely to operate more rapidly than neutrino cooling (see the discussion of timescales by Paschalidis et al. 2012).

Bauswein et al. (2010) carried out smoothed-particle hydrodynamics simulations of HMNSs in the conformal-flatness approximation to general relativity. They compared simulations using the full temperature dependence of the HShen and LS180 EOS¹² with an approximate treatment of thermal pressure via a Γ-law, P_th = (Γ_th − 1)_thρ_b. Although Bauswein et al. (2010) do not provide a figure showing the evolution of maximum baryon density, they show (in their Figure 5) graphs of cumulative mass as a function of distance from the center of the LS180-EOS HMNS at 8 ms after merger, roughly the time when the dynamical early postmerger phase is over and the secular HMNS evolution begins. From this, it may be observed that the HMNS with the lower thermal Γ (Γ_th = 1.5) is more compact than the model with Γ_th = 2. The HMNS evolved with the fully temperature-dependent LS180 EOS is in between the two, but closer to the Γ_th = 2 model. Bauswein et al. (2010) found that the more compact HMNS with Γ_th = 1.5 collapses after 10 ms, while the less compact Γ_th = 2.0 and full-LS180 cases collapse after ∼20 ms. This is consistent with the picture of secular HMNS evolution drawn in Section 5.2. Given a fixed number of baryons, a less compact configuration has a lower maximum baryon density after merger and, therefore, begins its secular evolution (in the sense of Figures 11 and 12) at a lower density than a more compact configuration. Consequently, it must lose more energy before reaching the critical density for collapse.

The above illustrates how thermal pressure effects may increase the lifetime of an HMNS by affecting the initial conditions for its secular evolution. From Section 4 one notes that hot configurations, at densities below ≲ 10¹⁵ g cm⁻³ (the exact value being EOS dependent), may support significantly larger masses than their cold counterparts at the same ρ_{b, max}. Thus, during the dynamical settle-down of two merging NSs to a secularly evolving HMNS remnant, a configuration with lower thermal pressure will need to evolve to higher ρ_{b, max} to reach an equilibrium configuration.