A New View on the Maximum Mass of Differentially Rotating Neutron Stars

D. Gondek-Rosińska; I. Kowalska; L. Villain; M. Ansorg; M. Kucaba

doi:10.3847/1538-4357/aa56c1

1. Introduction

Binary neutron star (BNS) mergers are thought to be promising sources of gravitational waves and of neutrinos (Rosswog 2015; Bernuzzi et al. 2016), as well as the progenitors of short gamma-ray bursts (Blinnikov et al. 1984; Eichler et al. 1989). After the detection of gravitational waves from binary black holes by the LIGO experiment (Abbott et al. 2016a, 2016b), they are even one of the most expected next targets. The outcome of a BNS merger is the formation of a stellar black hole, but the latter can either occur promptly or be slightly delayed and include a stage during which a massive and warm differentially rotating neutron star is produced (Shibata & Uryū 2002). Whatever is the actual scenario, highly sophisticated and realistic numerical simulations are needed to ascertain the signals to be expected and consequently to enable us to extract information on both the gravitational and the high-energy underlying physics. Although a lot of progress has been made, building numerical codes that take into account all the pertinent physical ingredients is such a difficult task that it is still far from being achieved, making it useful to resort to simplified studies concentrating on some specific aspects. One basic but key issue is the maximum life span of the potential short-lived material remnant, a question that can be approached by focusing first on its maximum mass. The analysis of the involved timescales and of the results of numerical simulations (Shibata et al. 2005) shows that a rough but reasonable approximation to address this problem consists in modeling the central body as a stationary axisymmetric rotating relativistic star in differential rotation, neglecting complicated inner motions of the matter, nuclear reactions, thermal effects, magnetic fields, etc. Doing so, it is, for instance, possible to study the influence of the degree of differential rotation or of the stiffness of the equation of state (EOS) on the maximum mass, which, for rotating stars, can be much higher than for static stars (see, e.g., Cook et al. 1992, 1994a, 1994b; Baumgarte et al. 2000; Lyford et al. 2003).

In a previous article (Ansorg et al. 2009, hereafter Paper I), a new investigation was presented of the structure of differentially rotating neutron stars, modelized as constant-density stars or relativistic $N=1({\rm{\Gamma }}=2)$ polytropes. This study, extended for other polytropic EOSs in Studzińska et al. (2016), relied on a multidomain spectral code (based on the so-called AKM-method; Ansorg et al. 2003a) that was formerly shown to enable the calculation of very extremal configurations of rigidly rotating relativistic stars (Ansorg et al. 2003b; Schöbel & Ansorg 2003) or rings (Ansorg 2005; Ansorg & Petroff 2005). In Paper I, only star-like configurations, i.e., with a spheroidal topology (without a hole), were considered, but allowing what are sometimes called "quasi-toroidal" configurations, in which the maximal density is not the central one. The focus was put on the solution space, and a noticeable result was the discovery of four "types" of configurations that coexist with each other even for reasonable profiles of angular momentum.

The purpose of the present article is to extend the study of Paper I by calculating, for ${\rm{\Gamma }}=2$ polytropes, astrophysically relevant quantities, such as the maximal mass, the angular momentum, the ratio between kinetic and potential energies, etc. Our highly accurate and stable spectral code enables us, for the first time, to study in detail those properties of differentially rotating neutron stars, taking into account the whole solution space identified in Paper I. Moreover, the understanding of the global structure of the parameter space makes it possible to explain the results of preceding studies, especially the works by Baumgarte et al. (2000) and Lyford et al. (2003), showing that strange features of some sequences they had obtained arise from the fact that their codes were jumping from one type of configuration to another, due to numerical limitations. Finally, the configurations we have calculated could be used as initial data to perform in a systematic way the stability analysis of differentially rotating neutron stars and to determine stability criteria for such objects.

The plan of the article is as follows. In Section 2, we start by recapitulating the issue of differential rotation in relativistic stationary rotating stars, with the primary goal of describing the hypothesis made in our work and defining variables and notations. Then, in Section 3, we briefly review current knowledge concerning the maximal mass of rotating relativistic stars and present our results in the context of the existence of the various types that were introduced in Paper I. Section 4 focuses on angular momentum and other quantities related to rotation and stability, and Section 5 summarizes the results achieved, in contrast with those of previous studies. Finally, the Appendix reviews the features of the numerical code, displaying tests of convergence and accuracy of our results, putting emphasis on the specificities of the version used in this article, as well as in Paper I.

2. Models of Differentially Rotating Relativistic Stars

Stationary and axisymmetric configurations of rotating relativistic stars have already been the subject of numerous works, be they analytical or numerical, making this topic quite classical (see Stergioulas 2003; Friedman & Stergioulas 2013, for reviews). Since in this article we use exactly the same equations and notations as in Paper I, we send the reader to that article for more details. Here, we shall only review the main assumptions concerning matter and its motion, in the framework of differentially rotating relativistic stars.

As in Paper I, we work with units such as $c=G=K=1$ , where c is the speed of light in vacuum, G is Newton's constant, and K is the polytropic constant. The latter is defined by the polytropic EOS, $p=K\,{\epsilon }_{B}^{{\rm{\Gamma }}}\,$ , with ${\rm{\Gamma }}=2$ in this article, which relates the pressure p to the rest-energy (or baryonic energy) density ${\epsilon }_{B}$ . For reasons that will be reviewed further, the main thermodynamical variable that shall be used is H, the dimensionless relativistic enthalpy defined from the pressure p and the total energy density as

$\begin{eqnarray}&&H(p)={\int }_{0}^{p}\,\displaystyle \frac{{dp}}{\epsilon (p)+p}.\end{eqnarray} \tag{ 1 }$

To smooth the way to comparison with other codes and systems of units, we remind the reader that for a ${\rm{\Gamma }}=2$ polytropic EOS with null temperature, the enthalpy (1) verifies, in units $G=c=1$ (but without any fixed values of K for the time being),

$\begin{eqnarray}&&H=\mathrm{log}(1\,+\,2\,K\,{\epsilon }_{B}),\end{eqnarray} \tag{ 2 }$

while the total energy density , which will be used in the following, satisfies

$\begin{eqnarray}&&\epsilon ({\epsilon }_{B})={\epsilon }_{B}\,+\,K\,{\epsilon }_{B}^{2},\end{eqnarray} \tag{ 3 }$

so that

$\begin{eqnarray}&&\epsilon (H)=\displaystyle \frac{1}{4K}({e}^{2H}\,-\,1).\end{eqnarray} \tag{ 4 }$

Describing the rotational properties is made much easier by introducing the fluid angular velocity with respect to infinity, ${\rm{\Omega }}={u}^{\phi }/{u}^{t}$ , where the u-variables are the nonzero components of the four-velocity and an auxiliary variable that we shall write as Y, defined as $Y={u}^{t}\,{u}_{\phi }$ , and which is some kind of specific angular momentum of the fluid. While the rotation profile of stationary and axisymmetric rotating Newtonian barotropes has to be such that the angular velocity only depends on the distance from the axis of rotation (Greenspan 1973), it is a well-known fact (Bardeen 1970; Butterworth & Ipser 1976) that for a relativistic one, conservation of the energy-momentum tensor implies the condition $Y=F({\rm{\Omega }}),$ where F is an arbitrary function. Such a weak constraint allows for numerous possibilities, some of which were, for instance, explored in Galeazzi et al. (2012) and Uryū et al. (2016). Here we should, however, restrict ourselves to the now classical law proposed by Komatsu et al. (1989), already used by many authors (e.g., Cook et al. 1992; Bonazzola et al. 1993; Goussard et al. 1998; Baumgarte et al. 2000; Lyford et al. 2003; Morrison et al. 2004; Villain et al. 2004) and in Paper I:

$\begin{eqnarray}&&F({\rm{\Omega }})={A}^{2}({{\rm{\Omega }}}_{c}\,-\,{\rm{\Omega }}),\end{eqnarray} \tag{ 5 }$

where A is a parameter with the dimension of a length, while it can be shown that ${{\rm{\Omega }}}_{c}$ is the limit of ${\rm{\Omega }}$ on the rotation axis. More precisely, we shall put Equation (5) in the form

$\begin{eqnarray}&&{\rm{\Omega }}={A}_{{\rm{\Omega }}}(Y)={{\rm{\Omega }}}_{c}\,-\,\displaystyle \frac{Y}{{(\hat{A}{r}_{e})}^{2}},\end{eqnarray} \tag{ 6 }$

where we introduced the equatorial radius of the star r_e, which is a typical length scale of the problem, and $\hat{A}=A/{r}_{e}$ . Notice that since the rotation profile tends to be rigid (or uniform) when $A\to \infty$ (r_e keeping a finished value), we follow Baumgarte et al. (2000) and parameterize the sequences and the degree of differential rotation with

$\begin{eqnarray}&&\tilde{A}={\hat{A}}^{-1}={r}_{e}/A.\end{eqnarray} \tag{ 7 }$

We refer the reader to Paper I and references therein for more details about the basic (astro)physical consequences of this law (Newtonian limit, etc.).

For a stationary and axisymmetric barotropic perfect fluid in rotation, Bianchi identities (or equivalently the relativistic Euler's equation) imply that there is a first integral of motion, which we choose to write as

$\begin{eqnarray}&&H={V}_{0}\,-\,V\,-\,{B}_{{\rm{\Omega }}}(Y),\end{eqnarray} \tag{ 8 }$

where

(i)
H is the relativistic enthalpy (1);
(ii)
V is defined by
$\begin{eqnarray}&&{u}^{t}={e}^{-V},\end{eqnarray} \tag{ 9 }$
V₀ being its value at the north pole; and
(iii)

$\begin{eqnarray}&&{B}_{{\rm{\Omega }}}(Y)={\int }_{0}^{Y}\,x\,\displaystyle \frac{{{dA}}_{{\rm{\Omega }}}(x)}{{dx}}\,{dx},\end{eqnarray} \tag{ 10 }$
with ${A}_{{\rm{\Omega }}}$ introduced in Equation (6).

Working with the law given by Equation (5), a configuration of a star corresponds to a solution of the first integral of motion together with the Einstein equations and a given EOS ( ${\rm{\Gamma }}=2$ polytrope here), which is obtained by fixing three parameters: the maximal enthalpy H_m, the central angular velocity ${{\rm{\Omega }}}_{c}$ , and the $\tilde{A}$ quantity. Due to the nonlinear nature of this system of equations, there is nevertheless not always uniqueness of the solution if those are indeed the fixed quantities, a feature that makes the solution space even richer than it is for rigidly rotating stars, as was shown and discussed in Paper I. Furthermore, we note that a sufficiently high degree of differential rotation can make the enthalpy (or equivalently some density) have its maximal value outside of the symmetry axis, so that its central value is not an optimum parameter to fix. In the code used for this work, resorting to a Newton–Raphson scheme enables us to easily change the fixed parameters (see the Appendix) and to explore the solution space by freely varying any parameters. Among the convenient quantities adapted to uniquely specify a fast-rotating star are the ratio between its polar and equatorial radii, ${r}_{\mathrm{ratio}}={r}_{p}/{r}_{e}$ , with $0\lt {r}_{\mathrm{ratio}}\lt 1$ , and the rescaled shedding parameter $\tilde{\beta }$ , defined in Paper I by

$\begin{eqnarray}&&\tilde{\beta }=\displaystyle \frac{\beta }{\beta +1},\end{eqnarray} \tag{ 11 }$

with (see Ansorg et al. 2003b)

$\begin{eqnarray}&&\beta =-{\displaystyle \frac{{r}_{e}^{2}}{{r}_{p}^{2}}\left.\frac{d({z}_{b}^{2})}{d({\rho }^{2})}\right|}_{\rho ={r}_{e}},\end{eqnarray} \tag{ 12 }$

where the equation $z={z}_{b}(\rho )$ describes the surface of the star. It can be verified that $0\lt \tilde{\beta }\lt 1$ , that $\tilde{\beta }\to 1$ when a star enters into the toroidal regime (i.e., ${r}_{\mathrm{ratio}}\to 0$ ), that $\tilde{\beta }\to 0$ in the mass-shedding limit, and that $\tilde{\beta }=1/2$ for a nonrotating spherical star.

As stated in the Introduction, we shall focus in this article on the issue of the structure of differentially rotating relativistic stars with maximal mass, which will be the topic of the next section. Having in mind astrophysical scenarios such as the collapses or the mergers described earlier, we shall in the following also briefly discuss quantities related to rotation and instabilities, such as the angular momentum, the ratio between kinetic and gravitational energies, or the Kerr parameter.

3. Maximum Mass of Neutron Stars

For a given EOS, the maximum mass of nonrotating neutron stars ( ${M}_{\max ,\mathrm{stat}}$ ) is obtained by solving the Tolman–Oppenheimer–Volkoff (TOV) equations, varying the value of the central energy density ${\epsilon }_{{\rm{c}}}$ . Calculations show that, for realistic EOSs, it falls in the range of 2–2.5 ${M}_{\odot }$ . In the case of rigid rotation, the centrifugal force enables its value to be higher by 12%–20% for neutron stars (e.g., Cook et al. 1994a, 1994b) and by 42%–44% for strange stars (Gondek-Rosińska et al. 2000; Gourgoulhon et al. 1999). Other studies established that differential rotation can even be much more efficient, making possible an increase of the maximum mass by more than 60%, especially for moderately stiff polytropic EOSs (Baumgarte et al. 2000; Lyford et al. 2003). As we will explain, our work demonstrates that, in the same conditions as considered in those previous studies, the actual effect of differential rotation can in fact be stronger.

3.1. Comparison with Previous Studies

Following Baumgarte et al. (2000), we study the rest mass M₀ of axisymmetric differentially rotating ${\rm{\Gamma }}=2$ polytropes by building sequences of configurations for fixed values of the degree of differential rotation $\tilde{A}$ and of the maximal energy density ${\epsilon }_{\max }$ in the range of ∼0–0.6. Each sequence is parameterized by the ratio of the polar to the equatorial radius ${r}_{\mathrm{ratio}}={r}_{p}/{r}_{e}$ . We start with a static, spherically symmetric configuration, ${r}_{\mathrm{ratio}}=1$ , and then construct the sequence by decreasing ${r}_{\mathrm{ratio}}$ or, equivalently, by increasing the central angular velocity ${{\rm{\Omega }}}_{c}$ or the angular momentum J, until we reach the mass-shedding limit or a star with ${r}_{\mathrm{ratio}}=0$ .

As was found in Paper I, for each fixed value of the maximum energy density ${\epsilon }_{\max }$ , there is a critical degree of differential rotation ${\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })$ , below which such a sequence terminates at the mass-shedding limit. Sequences of this kind were assigned the type A. On the other hand, when the degree of differential rotation $\tilde{A}$ is larger than ${\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })$ , sequences with fixed ${\epsilon }_{\max }$ starting from a static body do not end at the Keplerian limit, but enter into the regime of toroidal bodies, whose topology is not simply connected: a hole appears at their center when ${r}_{\mathrm{ratio}}=0$ . We call them type C sequences.

For each sequence, be it of type A or C, we looked for the maximum mass. Figure 1 shows the rest mass M₀ versus ${r}_{\mathrm{ratio}}$ for three examples of sequences with fixed maximal densities ${\epsilon }_{\max }$ and $\tilde{A}$ , and for configurations close to the maximum mass. The left panel corresponds to stars with a modest degree of differential rotation $\tilde{A}=0.5$ , which end at the mass-shedding limit and are consequently of type A. On the curves, the mass-shedding limit is reached for the smallest non-null value of ${r}_{\mathrm{ratio}}$ that is displayed and that our code always enabled us to calculate. As stated earlier, a configuration at the Keplerian limit is characterized by $\tilde{\beta }=0$ . The vanishing of this rescaled shedding parameter is evidence of the presence of a cusp at the equator for this kind of configuration. One consequence is that once we have obtained configurations in the neighborhood of the mass-shedding limit, it is easier, with the Newton–Raphson scheme, to reach the Keplerian limit by decreasing $\tilde{\beta }$ rather than ${r}_{\mathrm{ratio}}$ . As is illustrated by the examples depicted in Figure 1, the maximum mass is found close to but not at the mass-shedding limit (except in the case of uniformly rotating stars). More precisely, we observed that the smaller $\tilde{A}$ is, the closer the configuration with maximum mass is to the mass-shedding limit.

As for the right panel of Figure 1, it shows M₀ versus ${r}_{\mathrm{ratio}}$ for type C sequences with $\tilde{A}=1\gt {\tilde{A}}_{\mathrm{crit}}$ . For such sequences, the maximum mass is always obtained for ${r}_{\mathrm{ratio}}=0$ , independently of the value of ${\epsilon }_{\max }$ or $\tilde{A}$ . Notice, however, that, as in Paper I, we consider in this article only stars without a hole, which implies that we arbitrarily terminate this type of sequence at ${r}_{\mathrm{ratio}}=0$ , while the mass can probably become even larger for other nonsimply connected configurations with the same values of ${\epsilon }_{\max }$ and $\tilde{A}$ .

In Figure 2, we illustrate some of our results by showing the maximum rest mass ${M}_{0\max }$ as a function of the maximum energy density ${\epsilon }_{\max }$ for sequences starting from a nonrotating configuration and with fixed degrees of differential rotation, $\tilde{A}=0.5$ , 0.7, and 1. In other words, each curve represents the upper limit on the rest mass for a given $\tilde{A}$ . For those values, lines marked by $\tilde{A}=0.5$ and 0.7 contain maxima for sequences classified as type A, while for $\tilde{A}=1.0$ they are of type C. We also displayed the results for static stars (lowest line) and for rigidly rotating stars at the mass-shedding limit (corresponding to $\tilde{A}=0$ ), and, for comparison, we included calculations by Baumgarte et al. (2000) (dot-dashed lines). It can easily be observed that the agreement with their results is very good for small and moderate degrees of differential rotation $\tilde{A}$ and/or small values of the maximal energy density ${\epsilon }_{\max }$ . However, for larger degrees of differential rotation or energy densities, the discrepancy between their results and ours becomes more and more visible. This illustrates that with the high level of accuracy and of stability of our AKM-method-based code (see the Appendix), we were able to reach solutions with higher masses than could be considered in previous works based on other algorithms. For each fixed value of the degree of differential rotation $\tilde{A}$ , we also indicated in Figure 2 with a cross the maximum allowed mass (the maximum of maxima).

To illustrate further the differences between the configurations with the maximum allowed mass depending on the degree of differential rotation, we show in Figure 3 their shape for $\tilde{A}=0$ (rigid rotation), 0.5, and 0.7 (both belonging to type A sequences), but also $\tilde{A}=1.0$ (a type C solution). Since for rigid rotation the maximum mass is obtained at the Keplerian limit, the surface of the star is not smooth, but exhibits cusps along the equator. On the contrary, as soon as $\tilde{A}\ne 0$ , the maximum mass corresponds to a configuration whose outer surface is regular. As can be seen, the higher $\tilde{A}$ is, the farther from the shape with cusps the maximum mass configuration is. Notably, the configuration with $\tilde{A}=1.0$ has a toroidal shape, but belongs to spheroidal topology (being the last simply connected object in the type C sequence).

**Figure 3.** Isocontours of the relativistic enthalpy H in meridional cross sections of stars with the maximum allowed mass for rigidly rotating neutron stars at the mass-shedding limit ( $\tilde{A}=0$ , top left panel) and three fixed degrees of differential rotation $\tilde{A}=0.5$ (top right panel, type A), $\tilde{A}=0.7$ (bottom left panel, type A), and $\tilde{A}=1.0$ (bottom right panel, type C).
Download figure:
Standard image High-resolution image

**Figure 3.** Isocontours of the relativistic enthalpy H in meridional cross sections of stars with the maximum allowed mass for rigidly rotating neutron stars at the mass-shedding limit ( $\tilde{A}=0$ , top left panel) and three fixed degrees of differential rotation $\tilde{A}=0.5$ (top right panel, type A), $\tilde{A}=0.7$ (bottom left panel, type A), and $\tilde{A}=1.0$ (bottom right panel, type C).
Download figure:
Standard image High-resolution image

In Table 1, we summarize the properties of differentially rotating stars with maximum allowed mass, for $\tilde{A}$ , the degree of differential rotation, ranging from 0 to 1.5. For type A solutions, the higher $\tilde{A}$ is, the higher are the allowed mass ${M}_{0\max }$ , the compactness parameter $M/{R}_{\mathrm{circ}}$ (where M and ${R}_{\mathrm{circ}}$ are the gravitational mass and the circumferential stellar radius, respectively), the angular momentum J, the ratio between the kinetic and the gravitational potential energies T/W, and the Kerr parameter $J/{M}^{2}$ (which are indicators of the onset of instabilities for rotating stars; see Section 4). The higher $\tilde{A}$ is, the lower the maximum energy density is. Note that for configurations with the maximum mass belonging to the type A solution, the maximal density is always located in the stellar center. In contrast, for type C sequences, the maximum allowed mass ${M}_{0\max }$ is always obtained for ${r}_{\mathrm{ratio}}=0$ , and it is a decreasing function of $\tilde{A}$ . For such configurations, its highest value is obtained for the smallest possible value of $\tilde{A}$ compatible with type C ( $\tilde{A}={\tilde{A}}_{\mathrm{crit}}$ ), which is ∼0.7 for the ${\rm{\Gamma }}=2$ polytropic EOS. However, despite what could have been expected, ${M}_{0\max }$ is not a continuous function of $\tilde{A}$ , due to the existence of the different types, and more specifically to the ambiguity of the definition of the types for configurations with $\tilde{A}={\tilde{A}}_{\mathrm{crit}}$ . This ambiguity, as we shall see in the next subsection, is related to the fact that, as was shown in Paper I, there are in the solution space other types of sequences than those discussed up to now (types A and C). Furthermore, as can already be seen in Table 1, those types can be associated with even higher masses than the most common type A.

Table 1. Properties of Stars with Maximal Rest Mass ${M}_{0\max }$ for All Types of Sequences and for $\tilde{A}$ in $[0;1.5]$

Type	$\tilde{A}$	${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$	${M}_{0\max }$	$T/\| W\|$	r_p/r_e	${H}_{\max }$	${\epsilon }_{\max }({\epsilon }_{c}/{\epsilon }_{\max })$	J	$J/{M}^{2}$	$M/{R}_{\mathrm{circ}}$
A	0.0	1.000	0.206941	0.0832	0.58479	0.43773	0.350 (1)	0.02017	0.5690	0.17373
	0.1	1.027	0.207924	0.0856	0.57966	0.43702	0.349 (1)	0.02067	0.5773	0.17393
	0.2	1.108	0.210921	0.0926	0.56477	0.43524	0.347 (1)	0.02215	0.6014	0.17448
	0.3	1.240	0.216118	0.1042	0.54132	0.43151	0.343 (1)	0.02470	0.6389	0.17557
	0.4	1.422	0.223975	0.1204	0.51059	0.42494	0.335 (1)	0.02849	0.6871	0.17759
	0.5	1.657	0.235568	0.1419	0.47306	0.41433	0.323 (1)	0.03406	0.7440	0.18122
	0.6	1.959	0.253800	0.1708	0.42686	0.39772	0.304 (1)	0.04286	0.8094	0.18834
	0.7	2.507	0.295169	0.2222	0.35240	0.37329	0.306 (1)	0.06243	0.8921	0.21434

C	0.8	2.999	0.46319	0.2937	0.005	0.16416	0.097 (2.e–4)	0.1758	1.023	0.2504
	0.9	3.382	0.43357	0.2854	0.002	0.16750	0.100 (2.e–5)	0.1526	1.008	0.2450
	1.0	3.805	0.40851	0.2771	0.005	0.1720	0.103 (2.e–4)	0.1338	0.989	0.2415
	1.5	6.420	0.32590	0.2379	0.01	0.1968	0.121 (6.e–4)	0.0783	0.897	0.2275

B	0.4	1.785	0.721 s	0.336	0.035	0.152	0.089 (0.016)	0.422	1.078	0.270
	0.5	2.006	0.639 s	0.335	0.114	0.145	0.084 (0.26)	0.340	1.082	0.246
	0.6	2.223	0.571 s	0.331	0.144	0.140	0.081 (0.51)	0.277	1.088	0.222
	0.7	2.443	0.510 s	0.324	0.164	0.140	0.081 (0.75)	0.225	1.091	0.201

D	0.8	2.6279	0.4485	0.3116	0.1825	0.14242	0.08239 (0.944)	0.178	1.096	0.177

Note. In addition to the mass, for each configuration are displayed the ratio between the central and equatorial angular velocities ( ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ ), the ratio between the kinetic and gravitational binding energies ( $T/| W|$ ), the ratio between the polar and equatorial radii (r_p/r_e), the maximum enthalpy ${H}_{\max }$ , the maximum energy density ${\epsilon }_{\max }$ (with the central energy density ${\epsilon }_{c}$ in units of ${\epsilon }_{\max }$ ), the angular momentum J, the Kerr parameter $J/{M}^{2}$ (with M the gravitational mass), and the compactness parameter $M/{R}_{\mathrm{circ}}$ (with ${R}_{\mathrm{circ}}$ the circumferential radius). For more details on the accuracy, see the Appendix.

Download table as: ASCII Typeset image

3.2. Maximum Mass for All Types of Solutions

In Paper I, with a fixed ${\epsilon }_{\max }$ and a moderate degree of differential rotation, there are sequences of stars, without a static limit, coexisting with either type A or type C sequences. They belong to new types of sequences, called type B and type D, respectively, and exist only thanks to differential rotation. The four types of one-dimensional parameter sequences, denoted A, B, C, and D, are illustrated in the $({r}_{\mathrm{ratio}},\tilde{\beta })$ plane for fixed ${\epsilon }_{\max }=0.12$ in Figure 4, in which it can be seen⁴ that there are two threshold values of A, ${\tilde{A}}_{B}$ and ${\tilde{A}}_{D}$ , such that types A and B coexist for ${\tilde{A}}_{B}\lt \tilde{A}\lt {\tilde{A}}_{\mathrm{crit}}$ , while types C and D are found when ${\tilde{A}}_{D}\gt \tilde{A}\gt {\tilde{A}}_{\mathrm{crit}}$ . For a given ${\epsilon }_{\max }$ , the minimal value of the degree of differential rotation, ${\tilde{A}}_{B}$ , for which the type B exists is the minimum of the $\tilde{A}(\tilde{\beta })$ function for r_ratio close to zero, e.g. 0.001, and similarly ${\tilde{A}}_{D}$ is the maximum of the $\tilde{A}({r}_{\mathrm{ratio}})$ function for fixed $\tilde{\beta }=0$ . On the other hand, the curve with $\tilde{A}={\tilde{A}}_{\mathrm{crit}}$ is a separatrix that divides the plan into four domains. For a given ${\epsilon }_{\max }$ , the value of ${\tilde{A}}_{\mathrm{crit}}$ can be determined thanks to the fact that the $\tilde{A}({r}_{\mathrm{ratio}},\tilde{\beta })$ function possesses a saddle point that belongs to the separatrix and at which the four types coexist (see Paper I for details).

**Figure 4.** Typical structure of the solution space illustrating, in the $({r}_{\mathrm{ratio}},\tilde{\beta })$ plane, the various types of sequences for several values of the degree of differential rotation $\tilde{A}$ . The curves show the dependency between the shedding parameter $\tilde{\beta }$ and the ratio between polar and equatorial radii ${r}_{\mathrm{ratio}}$ for ${\rm{\Gamma }}=2$ polytropic stars with fixed maximal energy density ${\epsilon }_{\max }=0.12$ ( ${H}_{\max }=0.2$ ). The thick curve corresponds to the separatrix sequence with $\tilde{A}={\tilde{A}}_{\mathrm{crit}}=0.75904$ , which divides the diagram into four regions containing sequences of types: A (lower right corner), B (lower left corner), C (above separatrix), and D (between types A and B). Results are taken from Paper I.
Download figure:
Standard image High-resolution image

**Figure 4.** Typical structure of the solution space illustrating, in the $({r}_{\mathrm{ratio}},\tilde{\beta })$ plane, the various types of sequences for several values of the degree of differential rotation $\tilde{A}$ . The curves show the dependency between the shedding parameter $\tilde{\beta }$ and the ratio between polar and equatorial radii ${r}_{\mathrm{ratio}}$ for ${\rm{\Gamma }}=2$ polytropic stars with fixed maximal energy density ${\epsilon }_{\max }=0.12$ ( ${H}_{\max }=0.2$ ). The thick curve corresponds to the separatrix sequence with $\tilde{A}={\tilde{A}}_{\mathrm{crit}}=0.75904$ , which divides the diagram into four regions containing sequences of types: A (lower right corner), B (lower left corner), C (above separatrix), and D (between types A and B). Results are taken from Paper I.
Download figure:
Standard image High-resolution image

Having this in mind, the two types of sequences without a static limit, types B and D, are defined as follows:

1.
Type B indicates one-dimensional sequences that start at the mass-shedding limit ( $\tilde{\beta }=0$ ) but continuously enter into the toroidal regime ( ${r}_{\mathrm{ratio}}\to 0$ ), when fixing ${\epsilon }_{\max }$ and $\tilde{A}$ but varying another suited parameter. As a consequence, they are always characterized by small values of ${r}_{\mathrm{ratio}}$ , and, as in type C sequences, in our study they arbitrarily end at $({r}_{\mathrm{ratio}}=0,\tilde{\beta }=1)$ . They are found for ${\tilde{A}}_{B}\leqslant \tilde{A}\leqslant {\tilde{A}}_{\mathrm{crit}}$ .
2.
Type D also starts at the mass-shedding limit ( $\tilde{\beta }=0$ ), but, unexpectedly, they terminate there as well. As illustrated by Figure 4, configurations of this type fill a smaller part of the solution space and are less easily found than those of all other types. They appear for ${\tilde{A}}_{D}\geqslant \tilde{A}\geqslant {\tilde{A}}_{\mathrm{crit}}$ .

The three threshold values of $\tilde{A}$ are functions of ${\epsilon }_{\max }$ : ${\tilde{A}}_{B}({\epsilon }_{\max })\lt {\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })\lt {\tilde{A}}_{D}({\epsilon }_{\max })$ . Hence, another useful way to depict the various domains of (co)existence of the types of sequences is to indicate them in the ( ${\epsilon }_{\max },\tilde{A}$ ) plane. This is done for the ${\rm{\Gamma }}=2$ EOS in Figure 5, which shows that, for all reasonable values of ${\epsilon }_{\max },$ the solution space has more or less the same structure.

**Figure 5.** Regions of existence of A, B, C, and D types of differentially rotating neutron star sequences in the plane ( $\tilde{A},{\epsilon }_{\max }$ ). The central (red) curve indicates the critical value ${\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })$ , which corresponds to the separatrix, on which all types (A, B, C, and D) of solutions coexist. The green dashed line and the blue dot-dashed line correspond to the lower limit of existence of type B, ${\tilde{A}}_{B}$ , and to the upper limit, ${\tilde{A}}_{D}$ , of existence of type D, respectively. For $\tilde{A}\lt {\tilde{A}}_{B}({\epsilon }_{\max })$ or $\tilde{A}\gt {\tilde{A}}_{D}({\epsilon }_{\max })$ , only one type remains, respectively A or C.
Download figure:
Standard image High-resolution image

**Figure 5.** Regions of existence of A, B, C, and D types of differentially rotating neutron star sequences in the plane ( $\tilde{A},{\epsilon }_{\max }$ ). The central (red) curve indicates the critical value ${\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })$ , which corresponds to the separatrix, on which all types (A, B, C, and D) of solutions coexist. The green dashed line and the blue dot-dashed line correspond to the lower limit of existence of type B, ${\tilde{A}}_{B}$ , and to the upper limit, ${\tilde{A}}_{D}$ , of existence of type D, respectively. For $\tilde{A}\lt {\tilde{A}}_{B}({\epsilon }_{\max })$ or $\tilde{A}\gt {\tilde{A}}_{D}({\epsilon }_{\max })$ , only one type remains, respectively A or C.
Download figure:
Standard image High-resolution image

Having understood the global structure of the solution space, we were able to explore it in detail thanks to the high level of flexibility of our Newton–Raphson-based code (see the Appendix), which allows us to fix or vary any parameter. As can be guessed from Figure 4, the shedding parameter $\tilde{\beta }$ is well suited for finding type B sequences with fixed $\tilde{A}$ and ${\epsilon }_{\max }$ , at least when we are not too close to the entrance into the toroidal regime $({r}_{\mathrm{ratio}}=0,\tilde{\beta }=1)$ . Hence, we looked for the maximum mass of this type of sequences by following such curves in the solution space. In Figure 6, we show typical results obtained when studying the rest mass M₀ as a function of $\tilde{\beta }$ . As explained earlier, all lines start at the mass-shedding limit ( $\tilde{\beta }=0$ ) and end when entering into the toroidal regime ( $\tilde{\beta }=1$ ). It was found that for type B, the configuration with maximum mass was always at one or the other of these two positions in the solution space (which can easily be seen for the examples depicted in Figure 6). More precisely, when the maximum energy density ${\epsilon }_{\max }$ was sufficiently low, the configuration with maximum mass was at the Keplerian limit, while with increasing ${\epsilon }_{\max }$ it jumped to ( ${r}_{\mathrm{ratio}}=0,\tilde{\beta }=1$ ). One example of the shape of a star belonging to a type B sequence at the Keplerian limit and with ${r}_{\mathrm{ratio}}\sim 0$ was shown in Figure 1 of Paper I. It illustrates that the stability of our numerical code allows calculations of extreme configurations simultaneously strongly pinched and oblate.

**Figure 6.** Rest mass M₀ as a function of $\tilde{\beta }$ along sequences of type B with fixed ${\epsilon }_{\max }$ (and for all of them $\tilde{A}=0.8$ ).
Download figure:
Standard image High-resolution image

**Figure 6.** Rest mass M₀ as a function of $\tilde{\beta }$ along sequences of type B with fixed ${\epsilon }_{\max }$ (and for all of them $\tilde{A}=0.8$ ).
Download figure:
Standard image High-resolution image

Before commenting further on our results, and especially the values of the maximum mass obtained, we mention that a quite similar procedure was applied to look for the maximum mass of type D sequences. It led to typical results as shown in Figure 7, which illustrates that the maximum mass was always reached for one of the two mass-shedding configurations belonging to the sequence (more precisely, the mass was always for the one with the smallest ${r}_{\mathrm{ratio}}$ ).

**Figure 7.** Rest mass M₀ as a function of $\tilde{\beta }$ along sequences of type D with fixed ${\epsilon }_{\max }$ (and for all of them $\tilde{A}=0.8$ ).
Download figure:
Standard image High-resolution image

As was already explained, the structure of the solution space is such that, for fixed ${\epsilon }_{\max }$ and $\tilde{A}$ , there can be several types of configurations that exist simultaneously. This means that there is not always a unique solution to the system of equations for given values of the three parameters needed to calculate one, depending on what are the chosen parameters. More precisely, we have seen (Figure 5) that there is some overlap between the domains of existence of types A and B, but also between C and D, if one fixes ${\epsilon }_{\max }$ and varies $\tilde{A}$ . Consequently, when looking for the maximum mass for fixed ${\epsilon }_{\max }$ and $\tilde{A}$ , as illustrated by curves similar to those in Figure 2, all types of sequences have to be taken into account, and one should not only consider configurations reached from a static limit. As a matter of fact, when two neutron stars merge, which is one of the situations of interest for studies of maximum masses, a naive expectation could even be that the shape of the material remnant would, at first, be closer to that of a more toroidal type B configuration than to the more spherical shape of a type A star.

In Figure 8, we come back to some of the results already shown in Figure 2 and display the maximum rest mass ${M}_{0\max }$ as a function of the maximum energy density ${\epsilon }_{\max }$ for sequences with a fixed degree of differential rotation, $\tilde{A}=0.7$ (to make the comparison easier, we kept the curves for rigid rotation and for static stars, as well as the dot-dashed line associated with the calculations of Baumgarte et al. 2000). However, this time, in addition to the results obtained for type A sequences, we also included those for type B (whose existence for $\tilde{A}=0.7$ is proven by Figures 4 and 5).

From this figure, one concludes that type B stars can sustain a much higher mass than the more common type A configurations. If one compares with rigid rotation, the increase of the maximum mass can even be as large as around 150% (more than 0.5, compared to ∼0.2). The same kinds of conclusions arise from a careful study of the solution space for all types. For instance, Figure 9 shows the maximum rest mass ${M}_{0\max }$ versus the maximum energy density ${\epsilon }_{\max }$ for $\tilde{A}=0.8$ stars. As can be noticed from Figure 5, this value of the degree of differential rotation $\tilde{A}$ is one of the few for which all four types of solution can exist. More precisely, for low values of ${\epsilon }_{\max }$ , one has ${\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })\gt 0.8$ , so that the configurations are of types A and B, while for large values of ${\epsilon }_{\max }$ , ${\tilde{A}}_{\mathrm{crit}}({\epsilon }_{\max })\lt 0.8$ , and the configurations are of types C and D, the transition being for ${\epsilon }_{\max }\sim 0.08$ such that ${\tilde{A}}_{\mathrm{crit}}=0.8$ (see Figure 5).

One of the important conclusions that can be drawn from Figure 9 is that type C and D stars also give rise to masses much larger than those of rigidly rotating stars, even for a reasonable degree of differential rotation and maximum energy density. The precise values that we obtained for the highest increase of the maximum rest mass (for given $\tilde{A}$ ) within sequences of types C and D are presented in Table 1, together with other physical properties of those configurations that we describe in Section 4. Due to the quite small domain of existence of type D, be it in the ( ${\epsilon }_{\max },\tilde{A}$ ) plane (visible in Figure 5) or in an ( ${r}_{\mathrm{ratio}},\tilde{\beta }$ ) plane with fixed ${\epsilon }_{\max }$ and $\tilde{A}$ (visible in Figure 4), we included only one typical value for this type. Finally, we represented the highest increase of the maximum rest mass with respect to static configurations for all types of sequences in Figure 10, in which we also displayed the results of Lyford et al. (2003) for comparison.

**Figure 10.** Highest increase of the maximum rest mass with respect to static configurations as a function of $\tilde{A}$ for all types of sequences. As in Table 1, only one value is given for type D, due to the narrowness of its domain of existence (see text for more details). For comparison, we also display the results of Lyford et al. (2003).
Download figure:
Standard image High-resolution image

**Figure 10.** Highest increase of the maximum rest mass with respect to static configurations as a function of $\tilde{A}$ for all types of sequences. As in Table 1, only one value is given for type D, due to the narrowness of its domain of existence (see text for more details). For comparison, we also display the results of Lyford et al. (2003).
Download figure:
Standard image High-resolution image

To sum up our results, we found that the maximum mass of differentially rotating neutron stars depends on both the degree of differential rotation and the type of solution. From Figure 10, one deduces that the maximum mass is an increasing function of $\tilde{A}$ for a type A solution (associated with a low and modest degree of differential rotation) and a decreasing function for types B and C. Furthermore, configurations from the newly discovered types B and C can possess masses much larger than those of type A. More precisely, the highest increase of the maximal mass, around 4 times the maximal static mass, was obtained for a modest degree of differential rotation and for configurations belonging to type B sequences, which were not taken into account in other studies, mainly due to numerical limitations. Stars of type B sequences are indeed very oblate objects, with toroidal shapes such that ${r}_{\mathrm{ratio}}\lt 0.25$ . The agreement of our results with those of Lyford et al. (2003) is very good for type A solutions, and good for a modest degree of differential rotation, $\tilde{A}=0.8\mbox{--}1.0$ , for type C solutions. Our calculations of the maximum allowed mass are consequently the first that take into account all types of solutions. However, the obtained value of the maximum mass is much higher than the mass of the heaviest stars known to date. It is naturally an open question whether the considered configurations could be stabilized by differential rotation.

4. Rotational Properties

For astrophysical purposes, the question of the maximal mass of differentially rotating neutron stars is strongly related to that of their stability. Even without doing dynamical simulations or stability analysis, some basic conclusions can be drawn from the study of rotational properties of the stars, some of which are listed in Table 1. In this section, we shall focus on some of those properties, restricting the discussion to

1.
the comparison of the different types of configurations that all coexist if $\tilde{A}$ is fixed at $\tilde{A}\,\equiv \,0.8$ (as can be seen in Figure 5), a value of the degree of differential rotation that is an intermediate one; and
2.
their evolution for stars with maximal mass when $\tilde{A}$ changes (see Table 1).

4.1. Angular Momentum

The most obvious quantity to start with is angular momentum J, which is depicted in Figure 11 as a function of the maximal energy density ${\epsilon }_{\max }$ for stars with $\tilde{A}=0.8$ . A straightforward remark to make is that for types A and C the minimal value of J for fixed ${\epsilon }_{\max }$ is 0, while it is not for types B and D, which results from the fact that types A and C admit nonrotating limits while types B and D do not (see Paper I or Figure 4). Then, one can notice that for a given value of ${\epsilon }_{\max }$ , the angular momentum stored in a type B star is always larger than that in a type A star, as was already the case for the mass. Again, in agreement with what was the situation for the mass, the angular momentum of a type C star of fixed maximal energy density can be either higher or lower than for a type D star.

If one no longer fixes the maximal energy density, one notices that the values taken by the angular momentum can be much higher for type B, C, or D stars than for type A stars. Since such high values are reached for very small ratios of the radii and possibly for stars that do not have a nonrotating limit, this implies that calculations made only by accelerating a TOV star shall miss most of those configurations. Also, one can see that when the maximal energy density is close to the transition value that corresponds to the separatrix (see Paper I, Section 3, and Figure 5), there seems to be a huge gap in the maximal angular momentum for the type A star of highest maximal energy density and for the type C star of smallest maximal energy density, both being equal to ${\varepsilon }_{\max ,{\rm{S}}0.8}$ such that ${\tilde{A}}_{\mathrm{crit}}({\varepsilon }_{\max ,{\rm{S}}0.8})=0.8$ (see Figure 5). However, as we briefly discussed in Section 3, there is here an ambiguity in the definition of the types linked to the fact that for configurations exactly on the separatrix, types A and C share a branch: in Figure 4, it is the right one among the two that are going from the point of intersection of the separatrix lines (thick curve) to the mass-shedding limit (the horizontal axis defined by $\tilde{\beta }=0$ ).

Then, one shall notice that for stars of type A and D (but not B or C), the configurations with the highest value of J are also those with the maximal mass, so that they are (see Section 3)

1.
for type A: close to but not at the mass-shedding limit (with the exception of the case of rigid rotation);
2.
for type D: for the smallest value of ${r}_{\mathrm{rat}}$ among the two mass-shedding limit configurations (see Figure 7).

As far as types B and C are concerned, the situation is slightly different, maybe due to the fact that those types include stars with a vanishing polar radius (since we decided not to study stars with a hole), which have complicated internal distributions of physical quantities. More specifically, we observed that

1.
For type B: the maximal J is found at the mass-shedding limit, which is also the maximal mass for maximal energy densities not too close to the transition value (see Section 3). In Figure 11, this explained the continuity between the maximal J of types B and D when ${\epsilon }_{\max }$ increases. Notice that in this figure, the configurations of maximal mass for type B are indicated by a dotted line.
2.
For type C: configurations with maximal J are very close to those of maximal mass (at vanishing polar radius; see Section 3), but they do not exactly coincide with them.

Finally, if one comes back to Table 1 to have a glimpse at the influence of $\tilde{A}$ on J, one can see that increasing the degree of differential rotation (i.e., increasing $\tilde{A}$ ) allows higher maximal values of J only for stars of type A. For types B and C, it is the opposite, as was already the case with the maximal mass (see Section 3).

4.2. Ratio of Kinetic Energy to Gravitational Energy ( $T/| W|$ ) and Instabilities

For observational purposes, a quantity that is more directly interesting than the angular momentum is $T/| W|$ , the ratio between the kinetic (rotational) energy T and the gravitational binding energy W. This ratio indeed plays the role of an order parameter (Bertin & Radicati 1976) and is a good indicator of the possible onset of instabilities that can be the source of gravitational waves (Andersson 2003). It is displayed in Figure 12 for stars with $\tilde{A}\,\equiv \,0.8$ .

**Figure 12.** Ranges of the ratio between the kinetic and the potential energies as a function of the maximal energy density ${\epsilon }_{\max }$ for types A, B, C, and D with $\tilde{A}=0.8$ . Notice that for type B and the highest ${\epsilon }_{\max }$ , configurations of maximal $T/| W|$ and maximal mass do not coincide. The latter appear in this figure as a dotted line. Again, there is no gap between the lines of maxima for types B and D since on the separatrix they coincide at the mass-shedding limit of lowest ${r}_{\mathrm{rat}}$ (see text for more details).
Download figure:
Standard image High-resolution image

This figure shows many similarities with the figure for angular momentum J (Figure 11), for instance, in the facts that

1.
type B allows for higher values than type A;
2.
for types A and C the minimal value is 0, while it is not for types B and D;
3.
there seems to be a discontinuity of the maximal value when one goes from type A to type C, which is again due to the ambiguity in the definition of the types on the separatrix;
4.
for types A and D, the maximal value of $T/| W|$ is obtained for the configuration with maximal mass, whereas for type C they almost coincide, and for type B it happens only for the lowest ${\epsilon }_{\max }.$

However, one shall also notice some changes with respect to the situation for J:

1.
the relative difference between the possible values for type A and for other types is not as large for $T/| W|$ as it was for J;
2.
type C does not allow larger values than type D;
3.
for all types, the maximal value depends much less on ${\epsilon }_{\max }$ than it did for J.

As far as the influence of the degree of differential rotation is concerned, Table 1 tells us that in a way similar to what happens for J and M₀, the maximal value of $T/| W|$ is an increasing function of this degree (i.e., of $\tilde{A}$ ) only for stars of type A, and it is a decreasing one for types B and C. Of course, this table also confirms that much higher values of $T/| W|$ can be reached for types other than A.

Indeed, another result illustrated by Figure 12 is that for all configurations of types B and D, $T/| W|$ is at least equal to 0.2, which leads to the legitimate question of the (dynamical) stability of such stars (see Shibata et al. 2000; Andersson 2003). The study of this stability is beyond the scope of this work since it requires making dynamical simulations (Baumgarte et al. 2000; Shibata et al. 2000) or analysis of perturbations, but referring to previous studies, one can expect such dynamical instabilities as the so-called bar-mode instability. Furthermore, as we are dealing here with differentially rotating stars, another kind of instability could be triggered, the so-called low $T/| W|$ instability (see, e.g., Krüger et al. 2010 and references therein). As a relation between the latter and the appearance of a corotation point has been suggested (see Passamonti & Andersson 2015 and references therein), an easy way to get some more information on the possible stability of the configurations under study here is by depicting the ratio between central and equatorial angular velocities, ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ , as is done in Figure 13 for stars with $\tilde{A}=0.8$ .

**Figure 13.** Ranges of the ratio between the central and equatorial angular velocities ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ as a function of the maximal energy density ${\epsilon }_{\max }$ for types A, B, C, and D with $\tilde{A}=0.8$ .
Download figure:
Standard image High-resolution image

**Figure 13.** Ranges of the ratio between the central and equatorial angular velocities ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ as a function of the maximal energy density ${\epsilon }_{\max }$ for types A, B, C, and D with $\tilde{A}=0.8$ .
Download figure:
Standard image High-resolution image

Again, the picture is quite similar to the previous ones, with

1.
type B, which corresponds to higher values than type A;
2.
types A and C, which have 0 as a minimal value, whereas types B and D do not;
3.
an apparent discontinuity of the maximal value when one goes from type A to type C, which is still due to the ambiguity in the definition of the types on the separatrix.

However, the ratio ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ differs from the other quantities by the fact that, for types B and D, its maximal value is obtained for the configuration with maximal mass, while for types A and C they only almost coincide. Notice that the value of this ratio for the configuration of maximal mass is so close to its maximal value that we do not indicate it in Figure 13.

As far as the quantity ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ itself is concerned, this figure shows that types B, C, and D (and especially type C) allow large values (up to 5 for $\tilde{A}=0.8$ , a moderate degree of differential rotation), which implies a large window of possible corotation for instabilities such as those studied in Krüger et al. (2010) to be triggered. As we stated earlier, we shall not enter more into the details of this topic, and we just note, to conclude, that Table 1 displays, as one can expect, a strong and positive correlation between the value of ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ for the star with maximal mass and the value of $\tilde{A}$ .

4.3. Kerr Parameter $J/{M}^{2}$

To conclude our brief study of rotational quantities for stars with maximal mass or with a rotation profile characterized by $\tilde{A}=0.8$ , we shall look at the so-called Kerr parameter $J/{M}^{2}$ , whose value cannot be larger than 1 for a rotating black hole in general relativity. Indeed, it has been suggested (see, e.g., Giacomazzo et al. 2011) that supra-Kerrian stars, whose collapse would lead to a naked singularity, are somehow stabilized so that the cosmic censorship conjecture is respected. As previously done for other rotational quantities, we show in Figure 14 the Kerr parameter $J/{M}^{2}$ as a function of ${\epsilon }_{\max }$ for stars from all types and with $\tilde{A}=0.8$ , and we show this parameter for all stars with maximal mass but various values of $\tilde{A}$ and types in Table 1.

**Figure 14.** Ranges of the ratio between the angular momentum and the square of the gravitational mass, $J/{M}^{2}$ , the so-called Kerr parameter, which cannot be larger than 1 (limit indicated by the black horizontal line), for a rotating black hole in general relativity. It is displayed here as a function of the maximal energy density ${\epsilon }_{\max }$ for types A, B, C, and D with $\tilde{A}=0.8$ .
Download figure:
Standard image High-resolution image

**Figure 14.** Ranges of the ratio between the angular momentum and the square of the gravitational mass, $J/{M}^{2}$ , the so-called Kerr parameter, which cannot be larger than 1 (limit indicated by the black horizontal line), for a rotating black hole in general relativity. It is displayed here as a function of the maximal energy density ${\epsilon }_{\max }$ for types A, B, C, and D with $\tilde{A}=0.8$ .
Download figure:
Standard image High-resolution image

From the table or the figure, it can easily be checked that as soon as the density is large enough, no supra-Kerrian configuration exists, a conclusion quite similar to what was found in Giacomazzo et al. (2011). More precisely, for types A, B, and C, the Kerr parameter is a decreasing function of ${\epsilon }_{\max }$ . Additionally, we observe from Table 1 that for stars with maximal mass, $J/{M}^{2}$ is an increasing function of the degree of differential rotation ( $\tilde{A}$ ) for types A and B, whereas it is a decreasing one for type C. Furthermore, the Kerr parameter is always larger than 1 for type B stars, and it can be so for type D stars whose maximal density is not too large (such as the star with maximal mass displayed in Table 1). Although we previously saw that quantities such as $T/| W|$ or ${{\rm{\Omega }}}_{c}/{{\rm{\Omega }}}_{e}$ strongly suggest that they are not stable, the fact that the Kerr parameter is larger than 1 for such configurations could be a possible indication of their quasi-stability. Indeed, the dynamical simulations of Giacomazzo et al. (2011) showed that supra-Kerrian stars seem to be stabilized by differential rotation. If other conclusions from Giacomazzo et al. (2011) are correct, the final collapse of such stars would be associated with the excitation of various modes, making them very interesting sources of gravitational waves. Nevertheless, as was already stated, to be properly dealt with, this issue would require some dynamical study or some perturbative analysis, which are far beyond the scope of the current work.

5. Discussion of the Results and Conclusion

Using a highly accurate spectral code based on the Newton–Raphson scheme, we calculated configurations of relativistic differentially rotating neutron stars modeled as ${\rm{\Gamma }}=2$ polytropes for broad ranges of maximal densities and of the degree of differential rotation. We were able to fully explore the solution space for stars with a rotation profile described by the law proposed in Komatsu et al. (1989), although we considered only models with spheroidal topology (without a hole).

For the first time, the maximum mass and various other astrophysical quantities were calculated for all types of differentially rotating neutron stars, as defined in Paper I. The maximum mass of differentially rotating neutron stars was shown to depend not only on the degree of differential rotation but also on the type of the solution. Its value is an increasing function of the degree of differential rotation for type A solutions (associated with a low or modest degree of differential rotation) and a decreasing function for types B (with a modest degree of differential rotation) and C (with a modest or high degree of differential rotation). The highest increase of the maximal mass, 3–4 times the maximal nonrotating mass, is obtained for intermediate degrees of differential rotation, indicating that the corresponding configurations could be relevant in some astrophysical scenarios. Those configurations belong to sequences of type B that were not taken into account in previous studies, mainly due to numerical difficulties. In addition, the thorough investigation performed in the present article allowed us to understand the partial results obtained with other codes (Baumgarte et al. 2000; Lyford et al. 2003) and to show that the maximum possible mass of a differentially rotating neutron star could be much higher than previously thought, even for astrophysically reasonable configurations.

In order to try to go a step further in deciding whether configurations of the new types have pertinence in actual situations, we performed a rough analysis of their rotational parameters, such as their angular momentum and other quantities linked to the possible appearance of instabilities. We observed that the ratio between the kinetic and potential energies was indeed quite large for many newly discovered configurations, but we also noticed that so is their Kerr parameter (always higher than 1 for stars with the maximum mass belonging to types B and D and for some of type C), which could imply that they are somehow stabilized. However, the definitive answer to that question has to come from other analyses, be they perturbative or fully dynamical. A few hydrodynamical studies (Baumgarte et al. 2000; Shibata et al. 2000; Giacomazzo et al. 2011) have already shown that supra-Kerr stellar models seem dynamically stable but are subject to various secular instabilities, leading to the emission of gravitational waves. Another complementary approach would naturally be to use the configurations we have calculated as initial data to perform dynamical evolutions of differentially rotating neutron stars and to study the stability criteria for such objects.

To be of physical interest, the conclusions drawn in our study should also be supported by further investigations with more relativistic descriptions of the microphysics, such as the EOS. In Studzińska et al. (2016), we present results concerning the influence of the stiffness of a polytropic EOS on the various types of configurations and on their properties to examine how robust our results are. In other articles, we shall study the maximum mass of strange stars (M. Szkudlarek et al. 2016, in preparation) and analyze in detail the rotational properties of neutron stars described by realistic EOSs.

This work was partially supported by the Polish MNiSW grant no. N203 511238; by the FOCUS Programme of Foundation for Polish Science; by POMOST/2012-6/11 from the Program of Foundation for Polish Science co-financed by the European Union within the European Regional Development Fund; by the National Science Centre grants UMO-2014/14/M/ST9/00707, UMO-2013/01/ASPERA/ST/0001, and DEC-2013/08/M/ST9/00664; by the NewCompStar COST Action MP1304; and by the HECOLS International Associated Laboratory Astrophysics France–Poland program. Calculations were performed on the PIRXGW computer cluster funded by the Foundation for Polish Science within the FOCUS program.

Appendix: The Numerical Scheme

As in Paper I, the numerical calculations are done using a pseudo-spectral collocation point method that utilizes two domains: (i) a domain that covers the fluid's interior, and (ii) a spatially compactified domain describing the fluid's exterior. In order to avoid Gibbs phenomena, we choose the common boundary between the two domains to coincide with the surface shape of the fluid configuration. This shape is not known a priori but forms part of the elliptic "free" boundary value problem to be solved. Each of the two subdomains is characterized by a mapping

$\begin{eqnarray}&&{\varrho }^{2}={\varrho }_{k}^{2}(s,t),\qquad {z}^{2}={z}_{k}^{2}(s,t),\quad (s,t)\in {[0,1]}^{2},\end{eqnarray} \tag{ 13 }$

where ϱ and z are the coordinates used in the metric (see Paper I) and where k labels the subdomains ( $k\in \{0;1\}$ ). Here we have used the fact that the solutions are axially and equatorially symmetric, from which it follows that the metric coefficients are functions of the coordinate squares, ϱ² and z².

For the coordinate transformation

$\begin{eqnarray*}&&(s,t)\mapsto ({\varrho }^{2},{z}^{2})\end{eqnarray*}$

we take care of the fluid's unknown surface shape by means of a one-dimensional function G,

$\begin{eqnarray*}&&G:[0,1]\to {\mathbb{R}}.\end{eqnarray*}$

In particular, we write the following:

1.
Exterior subdomain, k = 0:
$\begin{eqnarray}&&{\varrho }_{0}^{2}(s,t)=t[{r}_{{\rm{e}}}^{2}-{r}_{{\rm{p}}}^{2}+{\xi }^{2}(s)],\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}&&{z}_{0}^{2}(s,t)=(1-t)[{\xi }^{2}(s)-{r}_{{\rm{p}}}^{2}]+G(t)-{r}_{{\rm{e}}}^{2}t,\end{eqnarray} \tag{ 15 }$
with
$\begin{eqnarray}&&\xi (s)={r}_{{\rm{p}}}+{r}_{{\rm{e}}}\displaystyle \frac{1-\sigma (s)}{\sigma (s)},\end{eqnarray} \tag{ 16 }$

$\begin{eqnarray}&&\sigma (s)=1-\displaystyle \frac{\sinh [(1-s)\mathrm{log}{\varepsilon }_{s}]}{\sinh (\mathrm{log}{\varepsilon }_{s})},\end{eqnarray} \tag{ 17 }$

$\begin{eqnarray}&&G(t=0)={r}_{{\rm{p}}}^{2},\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray}&&G(t=1)={r}_{{\rm{e}}}^{2}.\end{eqnarray} \tag{ 19 }$
Here r_p and r_e describe the polar and equatorial radii, respectively. The boundaries of the exterior domain are described by
$\begin{eqnarray}&&s=0:\,{\rm{Spatial}}\,{\rm{infinity}},\sqrt{{\varrho }^{2}+{z}^{2}}\to \infty ,\end{eqnarray} \tag{ 20 }$

$\begin{eqnarray}&&s=1:\,{\rm{Surface}}\,{\rm{of}}\,{\rm{the}}\,{\rm{fluid,}}\,{\rm{given}}\,{\rm{by}}{\rm{:}}\end{eqnarray} \tag{ 21 }$

$\begin{eqnarray}&&\qquad \{({\varrho }^{2},{z}^{2})=({r}_{{\rm{e}}}^{2}t,[G(t)-{r}_{{\rm{e}}}^{2}t]\,),t\in [0,1]\},\end{eqnarray} \tag{ 22 }$

$\begin{eqnarray}&&t=0:\,{\rm{Rotation}}\,{\rm{axis}},\varrho =0,\end{eqnarray} \tag{ 23 }$

$\begin{eqnarray}&&t=1:\,{\rm{Equatorial}}\,{\rm{plane}},{\text{}}z\,=\,0.\end{eqnarray} \tag{ 24 }$
2.
Interior subdomain, k = 1:
$\begin{eqnarray}&&{\varrho }_{1}^{2}(s,t)={r}_{{\rm{e}}}^{2}t,\end{eqnarray} \tag{ 25 }$

$\begin{eqnarray}&&{z}_{1}^{2}(s,t)=s[G(t)-{r}_{{\rm{e}}}^{2}t].\end{eqnarray} \tag{ 26 }$
The boundaries of the interior domain are described by
$\begin{eqnarray}&&s=0:\,{\rm{Equatorial}}\,{\rm{plane}},{\text{}}z\,=\,0,\end{eqnarray} \tag{ 27 }$

$\begin{eqnarray}&&s=1:\,{\rm{Surface}}\,{\rm{of}}\,{\rm{the}}\,{\rm{fluid}},{\rm{as}}\,{\rm{in}}\,{\rm{}}(21),\end{eqnarray} \tag{ 28 }$

$\begin{eqnarray}&&t=0:\,{\rm{Rotation}}\,{\rm{axis}},\varrho =0,\end{eqnarray} \tag{ 29 }$

$\begin{eqnarray}&&t=1:\,{\rm{Equator}},\varrho ={r}_{{\rm{e}}},z=0.\end{eqnarray} \tag{ 30 }$

The mapping of the exterior subdomain is chosen to resemble oblate spheroidal coordinates in which the entire class of Maclaurin spheroids exhibits a rapid spectral convergence rate. For general highly flattened relativistic stars we find, however, that the spectral convergence rate can be improved considerably by refining the exterior spectral mesh in the vicinity of the fluid's surface. We achieve this by rescaling the coordinate s and adjusting the free parameter ${\varepsilon }_{s}$ introduced in Equation (17). For an illustrative example see Figure 15.

**Figure 15.** Example for mappings of interior and exterior domains; see Equation (13). The coordinate transformations being chosen are specifically suited to extremely flattened configurations.
Download figure:
Standard image High-resolution image

In our pseudo-spectral collocation point scheme, all functions ${U}^{\kappa }$ ( $\kappa =0\ldots {n}_{\mathrm{eq}}-1$ ) to be determined by the free boundary value problem are considered at specific grid points $({s}_{k,i};{t}_{k,j})$ in the subdomains $k\in \{0;1\}$ . These grid points are given through

$\begin{eqnarray}&&{s}_{k,i}={\sin }^{2}\left(\displaystyle \frac{\pi i}{2({n}_{k}^{(s)}-1)}\right);\quad i=0\ldots {n}_{k}^{(s)}-1,\end{eqnarray} \tag{ 31 }$

$\begin{eqnarray}&&{t}_{k,j}={\sin }^{2}\left(\displaystyle \frac{\pi j}{2({n}_{k}^{(t)}-1)}\right);\quad j=0\ldots {n}_{k}^{(t)}-1.\end{eqnarray} \tag{ 32 }$

The integers ${n}_{k}^{(s)}$ and ${n}_{k}^{(t)}$ describe the number of grid points in the domain k with respect to the s- and t-directions (i.e., the spectral expansion orders). While ${n}_{0}^{(s)}$ may be different from ${n}_{1}^{(s)}$ , we assume the same numbers ${n}_{0}^{(t)}={n}_{1}^{(t)}$ of grid points at the common domain boundary.

We collect all function values

$\begin{eqnarray}&&{U}_{k,{ij}}^{\kappa }={U}^{\kappa }({s}_{k,i},{t}_{k,j}),\end{eqnarray} \tag{ 33 }$

as well as the values of the unknown surface function $G,{G}_{j}=G({t}_{1,j})$ , in order to build up a vector ${\boldsymbol{f}}$ . In addition, this vector is filled with two physical parameters that characterize, for a given EOS, the configuration. In particular, we choose them to be V_c (value of V at the center; see definition (9)) and ${{\rm{\Omega }}}_{c}$ (see definitions (5) and (6)). Note that it is sometimes possible to find more than one solution to a given pair of parameters.

The collection of elliptic equations valid in the subdomains, transition conditions at the common domain boundary, the vanishing pressure boundary condition at the fluid's surface, and certain parameter relations that one wishes to fulfill yield a discrete nonlinear system of the form

$\begin{eqnarray}&&{{\boldsymbol{F}}}^{(n)}({{\boldsymbol{f}}}^{(n)})=0,\end{eqnarray} \tag{ 34 }$

where n stands for the collection of all ${n}_{k}^{(s)}$ and ${n}_{k}^{(t)}$ ,

$\begin{eqnarray*}&&n=\{({n}_{k}^{(s)},{n}_{k}^{(t)});k=0;1\}.\end{eqnarray*}$

The dimension of this system is given by

$\begin{eqnarray}&&{n}_{\mathrm{total}}={n}_{\mathrm{eq}}\displaystyle \sum _{k=0}^{1}{n}_{k}^{(s)}{n}_{k}^{(t)}+{n}_{{\rm{G}}}+{n}_{\mathrm{par}},\end{eqnarray} \tag{ 35 }$

with ${n}_{{\rm{G}}}={n}_{0}^{(t)}$ and ${n}_{\mathrm{par}}=2$ . In particular, the transition conditions require the ${U}^{\kappa }$ to be continuous and to possess continuous normal derivatives. At domain boundaries that correspond to portions of the rotation axis or the equatorial plane, we require regularity conditions, which follow from the elliptic equations when specialized to this boundary. Via the integrated Euler Equation (8), the vanishing pressure boundary condition restricts the potentials at the fluid's surface. It adds n_G equations to the system. Finally, we may include specific parameter relations that we wish to be satisfied. For example, we could just prescribe certain values for the physical parameters contained in ${\boldsymbol{f}}$ . However, we also might wish to prescribe other parameters instead, say, rest mass M₀ and angular momentum J of the objects. For this reason we include the n_par physical parameters in the vector ${\boldsymbol{f}}$ and add the specific parameter relations to the system.

The solution ${{\boldsymbol{f}}}^{(n)}$ of the discrete algebraic system (34) describes the spectral approximation of the solution to the free boundary value problem. We find the vector ${{\boldsymbol{f}}}^{(n)}$ using a Newton–Raphson scheme,

$\begin{eqnarray}&&{{\boldsymbol{f}}}^{(n)}=\mathop{\mathrm{lim}}\limits_{m\to \infty }{{\boldsymbol{f}}}_{m}^{(n)},\end{eqnarray} \tag{ 36 }$

$\begin{eqnarray}&&{{\boldsymbol{f}}}_{m+1}^{(n)}={{\boldsymbol{f}}}_{m}^{(n)}-{[{{\boldsymbol{J}}}^{(n)}({{\boldsymbol{f}}}_{m}^{(n)})]}^{-1}{{\boldsymbol{F}}}^{(n)}({{\boldsymbol{f}}}_{m}^{(n)}),\end{eqnarray} \tag{ 37 }$

where the Jacobian matrix is given by

$\begin{eqnarray}&&{{\boldsymbol{J}}}^{(n)}=\displaystyle \frac{\partial {{\boldsymbol{F}}}^{(n)}}{\partial {{\boldsymbol{f}}}^{(n)}}.\end{eqnarray} \tag{ 38 }$

Note that for the convergence of the scheme a "good" initial guess ${{\boldsymbol{f}}}_{0}^{(n)}$ is necessary, which we provide through a known nearby function.

The linear step inside the Newton–Raphson solver, i.e., the solution of

$\begin{eqnarray*}&&{{\boldsymbol{J}}}^{(n)}\cdot \delta f=-{{\boldsymbol{F}}}^{(n)},\end{eqnarray*}$

is performed with the preconditioned "Biconjugate Gradient Stabilized (Bi-CGSTAB)" method (Barrett et al. 1993). A good convergence of this method requires a so-called preconditioning, which we construct in complete analogy to Ansorg et al. (2004) and Ansorg (2005) through a second- or fourth-order finite-difference representation of the Jacobian matrix of the nonlinear system.

For most configurations considered in this article, numerical solutions with extremely high accuracy were obtained with a moderate computational effort. This is illustrated by Figure 16, which displays the accuracy reached for some astrophysical parameters, e.g., the baryon mass, the circumferential radius, and the angular momentum, for a typical example of a not-too-oblate ( ${r}_{\mathrm{ratio}}=0.36$ ) type A star with a moderate degree of differential rotation ( $\tilde{A}=0.7$ ) and a modest maximum enthalpy ( ${H}_{\max }={H}_{{\rm{c}}}=0.38$ ). More precisely, in this figure are represented the relative differences between the nth spectral approximation S_n of all these quantities (denoted S) and their approximant of order 36, as functions of the number n of spectral points. In the case of more extreme configurations (e.g., close to the Keplerian limit, $\tilde{\beta }=0$ , or to the entrance in the toroidal regime, ${r}_{\mathrm{ratio}}=0$ ), the number of points needed to get a similar accuracy was larger, but the code also appeared able to perform the calculations. For subcritical configurations, n = 24 was sufficient, while for extreme ones, up to n = 34 was sometimes necessary. In Figure 16 it can be observed that the accuracy reached for such resolutions is of the order of ${10}^{-7}$ .

**Figure 16.** Illustration of the geometrical convergence rate for the rest mass M₀, the angular momentum J, and the circumferential radius ${R}_{\mathrm{circ}}$ of a differentially rotating neutron star described by a polytropic EOS with ${\rm{\Gamma }}=2$ . This configuration is of type A and characterized by a degree of differential rotation $\tilde{A}=0.7$ , a ratio of the coordinate radii ${r}_{\mathrm{ratio}}=0.36$ , and a maximal enthalpy ${H}_{\max }=0.38$ . For all three quantities, the plot displays the accuracy of the nth spectral approximation S_n, with $n={n}_{0}^{(s)}={n}_{1}^{(s)}={n}_{0}^{(t)}={n}_{1}^{(t)}$ .
Download figure:
Standard image High-resolution image

**Figure 16.** Illustration of the geometrical convergence rate for the rest mass M₀, the angular momentum J, and the circumferential radius ${R}_{\mathrm{circ}}$ of a differentially rotating neutron star described by a polytropic EOS with ${\rm{\Gamma }}=2$ . This configuration is of type A and characterized by a degree of differential rotation $\tilde{A}=0.7$ , a ratio of the coordinate radii ${r}_{\mathrm{ratio}}=0.36$ , and a maximal enthalpy ${H}_{\max }=0.38$ . For all three quantities, the plot displays the accuracy of the nth spectral approximation S_n, with $n={n}_{0}^{(s)}={n}_{1}^{(s)}={n}_{0}^{(t)}={n}_{1}^{(t)}$ .
Download figure:
Standard image High-resolution image

To conclude this appendix on the numerical scheme, we briefly describe the method used and the precision reached to identify differentially rotating neutron stars with maximum mass among all types of solutions, which was the main goal of this article. Figure 17 illustrates the precision and the method in the $({r}_{\mathrm{ratio}},{H}_{\max })$ plane for type A stars with $\tilde{A}=0.7$ . Once a rough estimation of the position in the $({r}_{\mathrm{ratio}},{H}_{\max })$ plane was obtained (after building sequences as described in the main text of the article), the code was used to find the value of M₀ for each configuration associated with $({r}_{\mathrm{ratio}}{}_{i},{H}_{\max }{}_{j})$ , where $i=0,1..,{i}_{\max }$ and $j=0,1..,{j}_{\max }$ , with typical values of ${i}_{\max }$ and ${j}_{\max }$ in the range of 5–15. Then the maximum mass was determined as the extremum of the ${M}_{0}({r}_{\mathrm{ratio}},{H}_{\max })$ function in that region of the plane, the corresponding configuration being also identified. Notice that when looking for the extremum of a smooth function such as ${M}_{0}({r}_{\mathrm{ratio}},{H}_{\max })$ , a spectral algorithm could also be used. The accuracy reached for the nth spectral approximation of ${M}_{0\max }$ is shown in the right panel of Figure 17, while the left panel presents the corresponding ${r}_{\mathrm{ratio}}$ and ${H}_{\max }$ . All values were obtained for ${i}_{\max }={j}_{\max }\,=\,10$ .