Gravitational-Wave Instabilities in Rotating Compact Stars

Abstract

It is generally accepted that the limit on the stable rotation of neutron stars is set by gravitational-radiation reaction (GRR) driven instabilities, which cause the stars to emit gravitational waves that carry angular momentum away from them. The instability modes are moderated by the shear viscosity and the bulk viscosity of neutron star matter. Among the GRR instabilities, the f-mode instability plays a historically predominant role. In this work, we determine the instability periods of this mode for three different relativistic models for the nuclear equation of state (EoS) named DD2, ACB4, and GM1L. The ACB4 model for the EoS accounts for a strong first-order phase transition that predicts a new branch of compact objects known as mass-twin stars. DD2 and GM1L are relativistic mean field (RMF) models that describe the meson-baryon coupling constants to be dependent on the local baryon number density. Our results show that the f-mode instability associated with $m=2$ sets the limit of stable rotation for cold neutron stars ($T\lesssim {10}^{10}$ K) with masses between $1M_{\odot }$ and $2M_{\odot }$. This mode is excited at rotation periods between 1 and 1.4 ms (∼20% to ∼40% higher than the Kepler periods of these stars). For cold hypothetical mass-twin compact stars with masses between $1.96M_{\odot }$ and $2.10M_{\odot }$, the $m=2$ instability sets in at rotational stellar periods between 0.8 and 1 millisecond (i.e., ∼25% to ∼30% above the Kepler period).

1. Introduction

Neutron stars are compact stellar objects formed when massive stars explode in so-called supernova explosions [1]. They possess very dense cores with densities up to $\sim {10}^{15}{\mathrm{g}/\mathrm{cm}}^3$ and canonical masses of around $1.4M_{\odot }$ [2,3]. However, these objects can reach masses larger than $2M_{\odot }$ and approach what is known as the Tolman-Oppenheimer-Volkoff (TOV) mass limit, an upper bound mass limit for cold, non-rotating neutron stars, analogous to the Chandrasekhar limit for white dwarfs. The typical radii of neutron stars are between 10 to 15 km. Neutron stars are frequently observed in the radio wavelength as pulsars. Pulsars are rapidly rotating neutron stars that emit electromagnetic radiation from their poles. The rotation frequencies of so-called fast pulsars are several hundred revolutions per second. Pulsars can emit radiation from large wavelengths like radio waves to small wavelengths like high-energy gamma rays, depending on their rate of rotation. Rotating neutron stars are referred to as pulsars. They can spin very rapidly, often making many hundred rotations per second. Many neutron stars form radio pulsars, emitting radio waves that appear from the Earth to pulse on and off like a lighthouse beacon as the star rotates at very high speeds. Neutron stars in X-ray binaries accrete material from a companion star and flare to life with tremendous bursts of X-rays [2]. Neutron stars can occur isolated or be part of binary stellar systems. One such configuration of stars is the Hulse-Taylor binary stellar system that consists of a neutron star and a pulsar [4,5]. Due to the emission of gravitational waves, the system loses energy and the orbits change exactly as predicted by Einstein’s general theory of relativity [5]. Einstein’s field equation is given by $G^{\mu }{}_{\nu }=8\pi G{T}^{\mu }{}_{\nu }$, where $G^{\mu }{}_{\nu }=R^{\mu }{}_{\nu }-1/2{\delta }^{\mu }{}_{\nu }R$ denotes the Einstein tensor, $T^{\mu }{}_{\nu }$ is the energy-momentum tensor of the stellar matter, and G is the gravitational constant. The quantity $R^{\mu }{}_{\nu }$ is the Ricci curvature tensor and R is the scalar curvature.

The thermodynamic variables that describe neutron star matter are energy density $\epsilon$, baryonic density n and pressure P at a given temperature T or entropy. To achieve hydrostatic stellar stability, the matter of neutron stars must have a neutral electric charge. Another important condition on the matter is $\beta$-equilibrium, which guarantees that the matter is in the lowest possible energy state where there are no net particle reactions. The equation of state (EoS) of neutron star matter can be expressed as $P(\epsilon ,T)$, where $\epsilon =\epsilon (n,T)$ and $P=P(n,T)$ acquire parametric forms. The EoS enters Einstein’s field equation through the energy-momentum tensor $T^{\mu }{}_{\nu }=(\epsilon +P)u^{\mu }{u}_{\nu }+{\delta }^{\mu }{}_{\nu }P$, where $u^{\mu }$ and $u_{\nu }$ denote four velocities. Despite many years of research, the composition of matter in the core of a neutron star is known only very vaguely. In the simplest conception, the matter consists of neutrons and a smaller number of protons whose charge is balanced by electrons and muons. However, because of the enormously high densities in the core, the neutron chemical potential may well exceed the in-medium masses of hyperons and even of the $\Delta \left(1232\right)$ baryon, in which case neutron stars would contain significant populations of hyperons and even of $\Delta \left(1232\right)$ particles [6]. Theoretical studies have also shown that quark matter, made of up, down, and strange quarks could be present in the cores of neutron stars (see [7] and references therein). Interestingly, successive first-order phase transitions to new states of matter can give rise to new branches of compact stars which are separated from the branch of conventional neutron stars by regions of instability [8]. One such branch comprises the mass twin stars studied in Ref. [9]. In this paper, we shall study rotational instabilities of mass twin stars and compare them with those of conventional neutron stars, which may help to distinguish both types of stars from each other observationally.

An absolute limit on rapid stellar rotation is set by the Kepler Frequency, at which mass shedding from the equator sets in, which destabilized the star. More stringent limits on rapid rotation are believed to be set by gravitational-radiation reaction (GRR) driven instabilities (in particular by f-modes or by r-modes), which lead to loss of angular momentum via the emission of gravitational radiation from a rotating star, if the star spins at the critical frequency at which these modes are excited [10,11,12,13]. This situation is shown schematically in Figure 1. We investigate the GRR driven f-modes [11,14,15]. The f-mode instability is expected to compete with the r-mode instability [12,13], which could set an even tighter limit on rapid stable rotation than the f-mode. Here we focus on the f-modes since they are currently better understood theoretically than the r-modes.

A comprehensive overview of rotational instabilities in rotating neutron stars can be found in the papers by [13,17]. It is noteworthy to mention that as yet no fully general relativistic calculations of the pulsation modes of rapidly rotating compact stars exist. In the absence of such studies, we look to follow the standard strategy [13] adopted to estimate the instability of timescales [11,18], which assumes that the true mode solution is well represented by the solution to the non-dissipative perturbation equations. These solutions are then used to evaluate the effects of different dissipation mechanisms (bulk and viscosity) and add their contributions to the rate of change of the mode energy according to

$$\frac{dE}{dt}=-2\frac{E}{\tau },$$

(1)

where $\tau$ is the timescale of the dissipative process and $dE/dt$ is determined by the equations for the evolution of a viscous fluid coupled to gravitational radiation [11,14,19]. The instability modes are taken to have the time dependence ($m=2,3,4,\dots$)

$$\propto \exp [i{\omega }_m\left(\mathrm{\Omega }\right)t+im\varphi -t/{\tau }_m\left(\mathrm{\Omega }\right)],$$

(2)

where ${\omega }_m$ is the frequency of the mode that depends on the rotation frequency $\mathrm{\Omega }$ of the star, $\varphi$ denotes the azimuthal angle, and ${\tau }_m$ is the time scale of the mode that determines its growth or damping. The problem of determining the maximum stable frequency of a compact star is then reduced to finding the root of [11,16]

$$0=\frac{1}{\tau \left(\mathrm{\Omega }\right)}\equiv \frac{1}{{\tau }_{\mathrm{GR}}\left(\mathrm{\Omega }\right)}+\frac{1}{{\tau }_{\eta }\left(\mathrm{\Omega }\right)}+\frac{1}{{\tau }_{\zeta }\left(\mathrm{\Omega }\right)},$$

(3)

where ${\tau }_{\mathrm{GR}}\left(\mathrm{\Omega }\right)$, ${\tau }_{\eta }\left(\mathrm{\Omega }\right)$, and ${\tau }_{\zeta }\left(\mathrm{\Omega }\right)$ are the dissipation time scales related to the gravitational-radiation reaction instability that drives the instability modes, and are related to the shear ($\eta$) and bulk ($\zeta$) viscosities that damp the modes. The critical frequencies obtained from Equation (3) depend on m and constitute the maximum rotational frequency (the smallest rotational period) of the stable compact star.

Our paper is organized as follows. In Section 2 we introduce the models for the EoS used in this work. The models can be divided into two different categories, namely relativistic field-theoretical models and a polytropic model which leads to the existence of mass-twin stars. The parameters of the models are provided in Section 3. Section 4 deals with the general relativistic stellar structure equations that follow Einstein’s field equation. The equations that are to be solved to determine the GRR-driven instabilities in rotating compact stars are reviewed in Section 5. The role of the shear and bulk viscosity for these instabilities are discussed in Section 6. A discussion along with a brief summary of the results is provided in Section 7.

2. Equation of State of Dense Neutron Star Matter

2.1. Choice of Lagrangian

Hadronic neutron star matter is modeled using the relativistic nuclear field theory that describes the baryon-baryon interaction through the exchange of mesons. The equations of the theory are solved in the framework of the standard relativistic mean-field (RMF) model. The standard RMF model describes the matter in the core of a neutron star as a relativistic quantum gas consisting of leptons ($\lambda \in \{e^{-},{\mu }^{-}\}$) and different types of baryons like nucleons, hyperons and the spin-$\frac{3}{2}$ delta isobar quartet ($B\in \{n,p,\mathrm{\Lambda },\mathrm{\Xi },\mathrm{\Sigma },\mathrm{\Delta }\}$) that interact via the exchange of scalar and vector mesons ($M\in \{\sigma ,\omega ,\rho \}$). The Lagrangian for interacting baryons is given by [15,20,21,22]

$${\mathcal{L}}_{\mathrm{B}}=\sum _B{\overline{\psi }}_B[{\gamma }_{\mu }(i{\partial }^{\mu }-g_{\omega B}\left(n\right){\omega }^{\mu }-\frac{1}{2}{g}_{\rho B}\left(n\right)\tau ·{\rho }^{\mu })-(m_B-g_{\sigma B}\left(n\right)\sigma )]{\psi }_B$$

(4)

where the scalar meson ($\sigma$) describes the attraction between baryons, the vector meson ($\omega$) describes the repulsion between baryons and the isovector meson ($\rho$) describes the baryon-baryon interactions in isospin asymmetric systems. Also, $g_{\omega B}\left(n\right)$, $g_{\rho B}\left(n\right)$ and $g_{\sigma B}\left(n\right)$ are the meson-baryon coupling constants, ${\psi }_B$ is the baryon field operator, n is the baryon number density, and $\mathit{\tau }$ = (${\tau }_1,{\tau }_2,{\tau }_3$) are the Pauli isospin matrices. The meson Lagrangian is given by [15,20,21,23,24]

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{M}}=\hfill & \hfill \frac{1}{2}({\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -m_{\sigma }^2{\sigma }^2)-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }+\frac{1}{2}{m}_{\omega }^2{\omega }_{\mu }{\omega }^{\mu }+\frac{1}{2}{m}_{\rho }^2{\mathit{\rho }}_{\mu }·{\mathit{\rho }}^{\mu }-\frac{1}{4}{\mathit{G}}_{\mu \nu }·{\mathit{G}}^{\mu \nu },\hfill \end{array}$$

(5)

where $F^{\mu \nu }={\partial }^{\mu }{\omega }^{\nu }-{\partial }^{\nu }{\omega }^{\mu }$ and ${\mathit{G}}^{\mu \nu }={\partial }^{\mu }{\mathit{\rho }}^{\nu }-{\partial }^{\nu }{\mathit{\rho }}^{\mu }$ denote the field tensors of $\omega$ and $\rho$ mesons, respectively. Additional nonlinear scalar self-interactions that contribute to a Lagrangian, allowing the standard RMF model to reproduce empirical values for the nuclear incompressibility ($K_0$) and effective nucleon mass ($m^{\ast }/m_N$) at $n_0$ have been considered in the literature [25,26],

$${\mathcal{L}}_{\mathrm{NL}\sigma }=-\frac{1}{3}{b}_{\sigma }{m}_N{\left[g_{\sigma N}\left(n\right)\sigma \right]}^3-\frac{1}{4}{c}_{\sigma }{\left[g_{\sigma N}\left(n\right)\sigma \right]}^4,$$

(6)

where $b_{\sigma }$ and $c_{\sigma }$ are constants fixed by the symmetric nuclear matter (SNM) parametrization. The parametrization determines the meson-nucleon coupling constants in the RMF model that reproduce properties of isospin SNM at the nuclear saturation density ($n_0$), which is given by

$$g_{iN}\left(n\right)=g_{iN}\left(n_0\right),i\in \{\sigma ,\omega ,\rho \}.$$

(7)

In order to soften the high-density EoS to satisfy the nuclear matter constraints and reproduce empirical values for the slope of the asymmetry energy ($L_0$) at $n_0$, the following nonlinear vector self-interactions and cross-interaction have also been suggested in the literature [27,28,29,30]

$$\begin{array}{ccc}{\mathcal{L}}_{\mathrm{NL}\omega }& =& \frac{1}{4}{g}_{{\omega }^4}\left(g_{\omega N}^2{\omega }_{\mu }{\omega }^{\mu }\right),\end{array}$$

(8)

$$\begin{array}{ccc}{\mathcal{L}}_{\sigma \omega \rho }& =& \frac{1}{2}{g}_{\rho N}^2{\mathit{\rho }}_{\nu }·{\mathit{\rho }}^{\nu }(g_{\sigma \rho }{g}_{\sigma N}^2{\sigma }^2+g_{\omega \rho }{g}_{\omega N}^2{\omega }_{\mu }{\omega }^{\mu }),\end{array}$$

(9)

where $g_{{\omega }^4}$, $g_{\sigma \rho }$ and $g_{\omega \rho }$ are coupling constants fixed by the SNM parametrization. These terms are not considered in this paper. The leptons present in the core of the neutron star along with the hadronic matter contribute due to the condition of the chemical equilibrium of neutron star matter, so the Lagrangian for non-interacting leptons is given by [15]

$${\mathcal{L}}_{\mathrm{L}}=\sum _{\lambda }{\overline{\psi }}_{\lambda }(i{\gamma }_{\mu }{\partial }^{\mu }-m_{\lambda }){\psi }_{\lambda }.$$

(10)

The total standard RMF Lagrangian considered in this paper is thus given by [15,20,21,22,31]

$${\mathcal{L}}_{\mathrm{RMF}}={\mathcal{L}}_{\mathrm{B}}+{\mathcal{L}}_{\mathrm{M}}+{\mathcal{L}}_{\mathrm{NL}\sigma }+{\mathcal{L}}_{\mathrm{L}}.$$

(11)

2.2. Dyson Equation and Baryon Self-Energies

The baryon and meson field equations are derived by evaluating the Euler-Lagrange equation for the particle fields ${\psi }_B,{\psi }_{\lambda },\sigma ,\omega$, and $\rho$ in ${\mathcal{L}}_{\mathrm{RMF}}$ [15,20,21]. For ${\psi }_B$ the result is the inhomogeneous Dirac equation given by

$$(i{\gamma }^{\mu }{\partial }_{\mu }-m_B){\psi }_B\left(x\right)=\left\{-g_{\sigma B}\sigma \left(x\right)+g_{\omega B}{\gamma }^{\mu }{\omega }_{\mu }\left(x\right)+g_{\rho B}{\gamma }^{\mu }\tau ·{\rho }_{\mu }\left(x\right)\right\}{\psi }_B\left(x\right).$$

(12)

The equation of motion for the $\sigma$-meson field follows as

$$({\partial }_{\mu }{\partial }^{\mu }+m_{\sigma }^2)\sigma \left(x\right)=\sum _B{g}_{\sigma B}{\overline{\psi }}_B\left(x\right){\psi }_B\left(x\right)-m_N{b}_N{g}_{\sigma N}{\left(g_{\sigma N}\sigma \left(x\right)\right)}^2-c_N{g}_{\sigma N}{\left(g_{\sigma N}\sigma \left(x\right)\right)}^3,$$

(13)

which constitutes an inhomogeneous Klein-Gordon equation. The equations of motion for the $\omega$ and $\rho$ fields read

$$\begin{array}{ccc}\hfill {\partial }^{\mu }{F}_{\mu \nu }\left(x\right)+m_{\omega }^2{\omega }_{\nu }\left(x\right)\hfill & =& \hfill \sum _B{g}_{\omega B}{\overline{\psi }}_B\left(x\right){\gamma }_{\nu }{\psi }_B\left(x\right),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill {\partial }^{\mu }{\mathit{G}}_{\mu \nu }\left(x\right)+m_{\rho }^2{\mathit{\rho }}_{\nu }\left(x\right)\hfill & =& \hfill \sum _B{g}_{\rho B}{\overline{\psi }}_B\left(x\right)\mathit{\tau }{\gamma }_{\nu }{\psi }_B\left(x\right),\hfill \end{array}$$

(15)

which both constitute homogeneous Proca equations. Next, we introduce the equations obeyed by the propagators ${\Delta }^{0\sigma }$, ${\Delta }^{0\omega }$, and ${\Delta }^{0\rho }$ of free $\sigma$, $\omega$, and $\rho$ meson, respectively. They are given by [15]

$$\begin{array}{ccc}\hfill ({\partial }^{\mu }{\partial }_{\mu }+m_{\sigma }^2)D^{0\sigma }(x,x^{\prime })\hfill & =& \hfill {\delta }^4(x-x^{\prime }),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill ({\partial }^{\lambda }{\partial }_{\lambda }{\delta }^{\mu }{}_{\nu }-{\partial }^{\mu }{\partial }_{\nu }+m_{\omega }^2{\delta }^{\mu }{}_{\nu })\left(g_{\mu \kappa }+\frac{{\partial }_{\mu }{\partial }_{\kappa }}{m_{\omega }^2}\right)D^{0\omega }(x,x^{\prime })\hfill & =& \hfill g_{\nu \kappa }{\delta }^4(x-x^{\prime }),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill ({\partial }^{\lambda }{\partial }_{\lambda }{\delta }^{\mu }{}_{\nu }-{\partial }^{\mu }{\partial }_{\nu }+m_{\rho }^2{\delta }^{\mu }{}_{\nu })\left(g_{\mu \kappa }+\frac{{\partial }_{\mu }{\partial }_{\kappa }}{m_{\rho ,r}^2}\right)D^{0\rho }(x,x^{\prime };r,r^{\prime })\hfill & =& \hfill g_{\nu \kappa }{\delta }^4(x-x^{\prime }){\delta }_{r{r}^{\prime }},\hfill \end{array}$$

(18)

where r$(=1,2,3)$ accounts for the electric charge carried by a $\rho$ meson. With the help of Equations (16) through (18), the meson field Equations (13) through (15) can be written as [15]

$$\begin{array}{ccc}\hfill \sigma \left(x\right)\hfill & =& \hfill \sum _B{g}_{\sigma B}\int d^4{x}^{\prime }{D}^{0\sigma }(x,x^{\prime }){\overline{\psi }}_B\left(x^{\prime }\right){\psi }_B\left(x^{\prime }\right),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\omega }_{\mu }\left(x\right)\hfill & =& \hfill \sum _B{g}_{\omega B}\int d^4{x}^{\prime }\left(g_{\mu \nu }+\frac{{\partial }_{\mu }{\partial }_{\nu }}{m_{\omega }^2}\right)D^{0\omega }(x,x^{\prime }){\overline{\psi }}_B\left(x^{\prime }\right){\gamma }^{\nu }{\psi }_B\left(x^{\prime }\right),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\rho }_{\mu }^r\left(x\right)\hfill & =& \hfill \sum _{r^{\prime }}\sum _B{g}_{\rho B}\int d^4{x}^{\prime }\left(g_{\mu \nu }+\frac{{\partial }_{\mu }{\partial }_{\nu }}{m_{\rho ,r}^2}\right)D^{0\rho }(x,x^{\prime };r,r^{\prime }){\overline{\psi }}_B\left(x^{\prime }\right){\tau }^{r^{\prime }}{\gamma }^{\nu }{\psi }_B\left(x^{\prime }\right).\hfill \end{array}$$

(21)

The Dyson equation of the baryon two-point Green function, $g_1\left(p\right)$, can then be written as

$$\begin{array}{c}\hfill \left({\gamma }^{\mu }{p}_{\mu }-m_B\right)g_1^B\left(p\right)=-1+{\mathrm{\Sigma }}^B\left(p\right)g_1^B\left(p\right),\hfill \end{array}$$

(22)

where the baryon self-energy ${\mathrm{\Sigma }}^B\left(p\right)$ is given by [15]

$$\begin{array}{ccc}\hfill {\mathrm{\Sigma }}^B\left(p\right)\hfill & =& \hfill i\sum _{B^{\prime }}\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{\mathrm{e}}^{i\eta q^0}\{g_{\sigma B}{D}^{0\sigma }\left(0\right)g_{\sigma B^{\prime }}+\left({\gamma }^{\mu }{g}_{\omega B}\right)D_{\mu \nu }^{0\omega }\left(0\right)\left({\gamma }^{\nu }{g}_{\omega B^{\prime }}\right)\hfill \\ & & \hfill +\left({\gamma }^{\mu }{\tau }^r{g}_{\rho B}\right)D_{\mu \nu }^{0\rho }(0;r,r^{\prime })\left({\gamma }^{\nu }{\tau }^{r^{\prime }}{g}_{\rho B^{\prime }}\right)\}{g}_1^{B^{\prime }}\left(q\right).\hfill \end{array}$$

(23)

Using the spectral representation of the 2-point Green function [15], one obtains for the self-energy of $\sigma$-mesons the expression

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\sigma }^B=\frac{g_{\sigma B}}{m_{\sigma }^2}\sum _{B^{\prime }}{g}_{\sigma B^{\prime }}{n}_{B^{\prime }}^S,\hfill \end{array}$$

(24)

where $n_{B^{\prime }}^S\equiv \langle {\overline{\psi }}_{B^{\prime }}{\psi }_{B^{\prime }}\rangle$ denotes the baryon scalar density which is given by

$$\begin{array}{ccc}\hfill n_B^S\hfill & =& \hfill \frac{2J_B+1}{2{\pi }^2}{\int }_0^{k_B}dk{k}^2\frac{m_B^{\ast }}{\sqrt{m_B^{\ast 2}+k^2}}\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& =& \hfill \sum _B\frac{2J_B+1}{2{\pi }^2}\left(\frac{m_B^{\ast }}{2}\right)\left[k_B{E}_B^{\ast }-m_B^{\ast 2}ln\left(\frac{E_B^{\ast }+k_B}{m_B^{\ast }}\right)\right].\hfill \end{array}$$

(26)

The quantity $m_B^{\ast }$ in Equations (25) and (26) denotes the effective (in-medium) mass of a baryon given by

$$m_B^{\ast }=m_B-{\mathrm{\Sigma }}_{\sigma }^B,$$

(27)

and $E_B^{\ast }$ in those equations stands for the effective single-particle energy of a baryon

$$\begin{array}{c}\hfill E_B^{\ast }\left(k\right)=\sqrt{m_B^{\ast 2}+k^2},\hfill \end{array}$$

(28)

evaluated at the Fermi surface, i.e., $E_B^{\ast }=\sqrt{m_B^{\ast 2}+k_B^2}$, where $k_B$ is the baryon Fermi momentum of a baryon. The quantity $J_B$ is the spin degeneracy factor of a baryon.

Next, we turn to the self energies of $\omega$ and $\rho$ mesons in Equation (23). They are given by

$$\begin{array}{ccc}\hfill {\mathrm{\Sigma }}_{\omega }^B\hfill & =& \hfill \frac{g_{\omega B}}{m_{\omega }^2}\sum _{B^{\prime }}{g}_{\omega B^{\prime }}{n}_{B^{\prime }},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\mathrm{\Sigma }}_{\rho }^B\hfill & =& \hfill \frac{g_{\rho B}}{m_{\rho }^2}\sum _{B^{\prime }}{g}_{\rho B^{\prime }}{I}_{3B^{\prime }}{n}_{B^{\prime }},\hfill \end{array}$$

(30)

where $I_{3B^{\prime }}$ is the 3-component of isospin. The quantity $n_{B^{\prime }}\equiv \langle {\psi }_{B^{\prime }}^{†}{\psi }_{B^{\prime }}\rangle$ stands for the baryon number density given by

$$\begin{array}{c}\hfill n_B=(2J_B+1)\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\mathrm{\Theta }(k_B-k),\hfill \end{array}$$

(31)

which leads to

$$n_B=\frac{2J_B+1}{6{\pi }^2}{k}_B^3,$$

(32)

where $k_B$ is the Fermi momentum of a baryon of type B.

The self energies shown in Equations (24), (29) and (30) are to be computed under the condition that neutron star matter is electrically charge neutral, and in chemical equilibrium. The first condition leads to

$$\sum _B{n}_B{q}_B+\sum _{\lambda }{n}_{\lambda }{q}_{\lambda }=0,$$

(33)

where $q_B$ and $q_{\lambda }$ are the baryon number and electric charge in units of elementary charge, respectively. The condition of chemical equilibrium reads

$${\mu }_B={\mu }_n-q_B{\mu }_e,$$

(34)

which expresses the chemical potential ${\mu }_B$ of a baryon in terms of the independent chemical potentials of neutrons, ${\mu }_n$, and electrons ${\mu }_e$. The latter is given by ${\mu }_e=E_e$, where $E_e=\sqrt{m_e^2+k_e^2}$ is the energy-momentum relation of free, relativistic electrons.

The densities shown in Equations (25) and (32) are for compact star matter at zero temperature. In the case of the matter at finite temperature, these densities are given by [32]

$$n_B^s=(2J_B+1)\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\left[f_{B-}\left(k\right)-f_{B+}\left(k\right)\right]\frac{m_B^{\ast }}{E_B^{\ast }\left(k\right)}$$

(35)

and

$$n_B=(2J_B+1)\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\left[f_{B-}\left(k\right)-f_{B+}\left(k\right)\right],$$

(36)

where the quantity $f_B^{\pm }$ denotes the Fermi-Dirac distribution function,

$$f_B^{\pm }\left(k\right)=\frac{1}{{\mathrm{e}}^{\mp (E_B^{\ast }\left(k\right)-{\mu }_B^{\ast })/T}+1}.$$

(37)

The effective baryon chemical potential, ${\mu }_B^{\ast }$, in Equation (37) is defines as

$${\mu }_B^{\ast }={\mu }_B-{\mathrm{\Sigma }}_{\omega }^B-{\mathrm{\Sigma }}_{\rho }^B{I}_{03}=\sqrt{m_B^{\ast 2}+k_B^2}$$

(38)

and ${\mu }_B$ is given by

$${\mu }_B={\mathrm{\Sigma }}_{\omega }^B+{\mathrm{\Sigma }}_{\rho }^B{I}_{03}+\sqrt{m_B^{\ast 2}+k_B^2}.$$

(39)

The total baryon number, n, of neutron star matter is given by

$$n=\sum _B{n}_B.$$

(40)

2.3. Equation of State in Standard RMF Theory

The RMF energy density ($\epsilon$) and pressure (P) are determined from the energy-momentum tensor which is given in terms of the RMF Lagrangian,

$$T^{\mu \nu }=g^{\mu \nu }\mathcal{L}+\sum _B\frac{\partial \mathcal{L}}{\partial \left({\partial }_{\mu }{\psi }_B\right)}{\partial }^{\nu }{\psi }_B,$$

(41)

where ${\psi }_B$ represents all baryons present in the neutron star matter. Spacetime is flat on the length scale of particle interactions inside a neutron star (∼1 fm), so the flat space Minkowski metric can be used. This leads to the following expressions for the energy density $\epsilon =\langle T^{00}\rangle$ to [15,22,32,33]

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{RMF}}\hfill & \hfill =\sum _B\frac{2J_B+1}{2{\pi }^2}{\int }_0^{k_B}dk{k}^2\sqrt{k^2+m_B^{\ast 2}}+\sum _{\lambda }\frac{1}{{\pi }^2}{\int }_0^{k_{\lambda }}dk{k}^2\sqrt{k^2+m_{\lambda }^2}\hfill \\ & \hfill +\frac{1}{2}\sum _B\left({\mathrm{\Sigma }}_{\sigma }^B{n}_B^S\right)+\frac{1}{2}\sum _B\left({\mathrm{\Sigma }}_{\omega }^B{n}_B\right)+\frac{1}{2}\sum _B\left({\mathrm{\Sigma }}_{\rho }^B{I}_{3B}{n}_B\right)-{\mathcal{L}}_{\mathrm{NL}\sigma }\left({\mathrm{\Sigma }}_{\sigma }^N\right),\hfill \end{array}$$

(42)

where ${\mathrm{\Sigma }}_{\sigma }^B$, ${\mathrm{\Sigma }}_{\omega }^B$, and ${\mathrm{\Sigma }}_{\rho }^B$ are the baryon self-energies given in Equations (24), (29), and (30), respectively. The expression for the total pressure $P_{\mathrm{RMF}}=\frac{1}{3}{\sum }_{i=1}^3\langle T^{ii}\rangle$ has the form [15,22,32,33]

$$\begin{array}{cc}\hfill P_{\mathrm{RMF}}\hfill & \hfill =\frac{1}{3}\sum _B\frac{2J_B+1}{2{\pi }^2}{\int }_0^{k_B}\frac{k^{4}dk}{\sqrt{k^2+m_B^{\ast 2}}}+\frac{1}{3}\sum _{\lambda }\frac{1}{{\pi }^2}{\int }_0^{k_{\lambda }}\frac{k^{4}dk}{\sqrt{k^2+m_{\lambda }^2}}\hfill \\ \hfill \hfill & \hfill -\frac{1}{2}\sum _B\left({\mathrm{\Sigma }}_{\sigma }^B{n}_B^S\right)+\frac{1}{2}\sum _B\left({\mathrm{\Sigma }}_{\omega }^B{n}_B\right)+\frac{1}{2}\sum _B\left({\mathrm{\Sigma }}_{\rho }^B{I}_{3B}{n}_B\right)+{\mathcal{L}}_{\mathrm{NL}\sigma }\left({\mathrm{\Sigma }}_{\sigma }^N\right).\hfill \end{array}$$

(43)

2.4. Density-Dependent RMF Models

Density-dependent RMF (DDRMF) models are an extension of the standard RMF model approach that accounts for meson-baryon coupling constants which are dependent on the local baryon number density n [23,24,34],

$$g_{iB}\left(n\right)=g_{iB}\left(n_0\right)f_i\left(x\right),$$

(44)

where i∈$\{\sigma ,\omega ,\rho \},$$x=n/n_0$, and $f_i\left(x\right)$ provides the functional form for the density dependence. The quantity $n_0=0.149$ fm${}^{-3}$ denotes the baryon number density of nuclear matter at saturation density. The most commonly utilized ansatz for $f_i\left(x\right)$ are given by

$$f_i\left(x\right)=a_i\frac{1+b_i{(x+d_i)}^2}{1+c_i{(x+d_i)}^2}$$

(45)

for i∈$\{\sigma ,\omega \}$, and to

$$f_{\rho }=exp[-a_{\rho }(x-1)]$$

(46)

for the $\rho$ meson. The density dependence of the meson-baryon couplings eliminate the need for the additional scalar nonlinear self-interactions in Equation (6) [22,23,24,33,35,36,37]. The Lagrangian of the model can thus formally written as

$${\mathcal{L}}_{\mathrm{DDRMF}}={\mathcal{L}}_{\mathrm{RMF}}-{\mathcal{L}}_{\mathrm{NL}\sigma }.$$

(47)

The novel term entering the DDRMF model is the so-called rearrangement contribution $\tilde{R}$ [35],

$$\tilde{R}=\sum _B\left(\frac{1}{g_{\omega B}}\frac{\partial g_{\omega B}}{\partial n}{n}_B{\mathrm{\Sigma }}_{\omega }^B+\frac{1}{g_{\rho B}}\frac{\partial g_{\rho B}}{\partial n}{I}_{3B}{n}_B{\mathrm{\Sigma }}_{\omega }^B-\frac{1}{g_{\sigma B}}\frac{\partial g_{\sigma B}}{\partial n}{n}_B^S{\mathrm{\Sigma }}_{\sigma }^B\right).$$

(48)

The rearrangement term modifies the expressions of the chemical potential as shown below [22,23,33],

$${\mu }_B={\mathrm{\Sigma }}_{\omega }^B+{\mathrm{\Sigma }}_{\rho }^B{I}_{03}+\sqrt{m_B^{\ast 2}+k_B^2}+\tilde{R}.$$

(49)

The expressions for the energy density and pressure are derived from the standard RMF by the removal of the nonlinear self-interactions and the contribution of the rearrangement energy. They are given by [22,33],

$$\begin{array}{ccc}\hfill {\epsilon }_{\mathrm{DDRMF}}\hfill & =& \hfill {\epsilon }_{\mathrm{RMF}}+{\mathcal{L}}_{\mathrm{NL}\sigma },\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill P_{\mathrm{DDRMF}}\hfill & =& \hfill P_{\mathrm{RMF}}-{\mathcal{L}}_{\mathrm{NL}\sigma }+n\tilde{R},\hfill \end{array}$$

(51)

where ${\epsilon }_{\mathrm{RMF}}$ and $P_{\mathrm{RMF}}$ are given by Equations (42) and (43), respectively.

In this paper, we use two parametrizations named GM1L and DD2 [22,36] for the RMF models. The equations of state of these parametrizations are shown in Figure 2, together with the polytropic model ACB4 that will be discussed in Section 2.5. The ACB4 EoS differs from DD2 and GM1L, in that, it contains a phase transition that leads to an additional stable configuration of massive objects called mass-twin stars. GM1L is an RMF model where only the $\rho$ meson coupling depends on density. In contrast to this, the couplings of all meson ($\sigma$, $\omega$, $\rho$) are treated as density-dependent. GM1L is also different from a standard RMF model because it is parametrized to fix the slope of the asymmetry energy, $L_0$, at $n_0$ [33]. The value of this slope has become constrained in recent years and may have a significant effect on neutron star composition and properties [38,39,40,41,42,43,44]. The isovector meson-baryon coupling constant of the GM1L model becomes density-dependent in the same way as the DDRMF approach given by Equation (25), which is convenient because $L_0$ can then be tailored by adjusting the coefficient $a_p$ without affecting the other saturation properties and does not impair the preexisting core parametrization [45].

2.5. A Polytropic Model for the EoS of Ultradense Matter

The ACB4 model exhibits a fundamentally different approach to a neutron star’s EoS than that of DD2 and GM1L, which are RMF-derived models. The ACB4 model consists of a piece-wise ($i=1,2,3,4$) polytropic representation of a neutron star’s EoS at super-nuclear densities [9,46]

$$P\left(n\right)=k_i{(n/n_0)}^{{\mathrm{\Gamma }}_i},n_i<n<n_{i+1},i=1\dots 4,$$

(52)

where ${\mathrm{\Gamma }}_i$ is the polytropic index in each of the density regions labeled by $i=1\dots 4$. The quantity ${\mathrm{\Gamma }}_i$ represents a stiff nucleonic EoS [47], the second polytrope corresponds to a first-order phase transition where pressure is constant in this region, $P_c={\kappa }_2$ (${\mathrm{\Gamma }}_2=0$), and polytropes 3 and 4 represent regions that lie above the phase transition with that correspond to high-density matter, like stiff quark matter [9]. The pressure, energy density, and number density of the model are given by [48],

$$P\left(n\right)=n^2\frac{d\left(\epsilon \right(n)/n)}{dn},$$

(53)

$$\epsilon \left(n\right)/n=\int dn\frac{P\left(n\right)}{n^2}=\frac{1}{n_0^{{\mathrm{\Gamma }}_i}}\int dn\kappa n^{{\mathrm{\Gamma }}_i-2}=\frac{1}{n_0^{{\mathrm{\Gamma }}_i}}\frac{\kappa n^{{\mathrm{\Gamma }}_i-1}}{{\mathrm{\Gamma }}_i-1}+C,$$

(54)

$$\mu \left(n\right)=\frac{P\left(n\right)+\epsilon \left(n\right)}{n}\frac{1}{n_0^{{\mathrm{\Gamma }}_i}}\frac{\kappa {\mathrm{\Gamma }}_i}{{\mathrm{\Gamma }}_i-1}{n}^{{\mathrm{\Gamma }}_i-1}+m_0,$$

(55)

where integration constant C is fixed by the condition that $\epsilon (n\to 0)=m_{0}n$. The above expressions are inverted to obtain the equations for number density and pressure as a function of chemical potential, $\mu$,

$$n\left(\mu \right)={\left(n_0^{{\mathrm{\Gamma }}_i}(\mu -m_0)\frac{{\mathrm{\Gamma }}_i-1}{\kappa {\mathrm{\Gamma }}_i}\right)}^{1/({\mathrm{\Gamma }}_i-1)},$$

(56)

$$P\left(\mu \right)=\kappa {\left(n_0^{{\mathrm{\Gamma }}_i}(\mu -m_0)\frac{{\mathrm{\Gamma }}_i-1}{\kappa {\mathrm{\Gamma }}_i}\right)}^{{\mathrm{\Gamma }}_i/({\mathrm{\Gamma }}_i-1)}.$$

(57)

This ACB4 EoS predicts a new family (branch) of stable compact stars, known as mass-twin stars, at densities that are higher than those of neutron stars, which will be discussed in more detail Section 4.1. The new branch has its origin in the phase transition highlighted in Figure 3. This type of phase transition is not present in the relativistic mean-field models of Section 2.3 which, therefore, do not lead to mass-twin stars. The parameters of the ACB4 model are shown in

Table 1. The parameters of the ACB4 model [9,46].

i	${\mathbf{\Gamma }}_{\mathit{i}}$	${\mathit{\kappa }}_{\mathit{i}}$ (MeV/fm${}^3$)	n${}_{\mathit{i}}$ (fm${}^{-3}$)	m${}_{0,\mathit{i}}$ (MeV)
1	4.921	2.1680	0.1650	939.56
2	0.0	63.178	0.3174	939.56
3	4.000	0.5075	0.5344	1031.2
4	2.800	3.2401	0.7500	958.55

[46].

3. The Model Parameters of the RMF Models

The values of the density-dependent parameters ($a_M$, $b_M$, $c_M$, $d_M$), meson masses ($m_{\sigma }$, $m_{\omega }$, $m_{\rho }$) and the meson-nucleon coupling constants ($g_{\sigma \mathrm{N}}\left(n_0\right)$, $g_{\omega \mathrm{N}}\left(n_0\right)$, $g_{\rho \mathrm{N}}\left(n_0\right)$) are adjusted to the charge and diffraction radii, spin-orbit splittings, the neutron skin thickness of finite nuclei and the bulk properties of symmetric nuclear matter (SNM) at $n_0$ [23,49,50]. The bulk properties of SNM at $n_0$ include the energy per nucleon, $E_0$, effective nucleon mass $m_N^{\ast }/m_N$, nuclear incompressibility $K_0$, nucleon potential $U_N$, symmetry energy $J_0$ and its derivative $L_0$. The values of these quantities for the nuclear parametrizations DD2 and GM1L are compiled in

Table 2. Properties of nuclear matter at nuclear saturation density for DD2 and GM1L [36].

Saturation Properties	Units	DD2 Model	GM1L Model
$E_0$	MeV	$-16.02$	$-16.30$
$K_0$	MeV	242.7	300.0
$m^{\ast }/m_N$	−	0.56	0.70
J	MeV	32.8	32.5
$L_0$	MeV	55.3	55.0
$n_0$	${\mathrm{fm}}^{-3}$	0.149	0.153
$-U_N$	MeV	75.2	65.6

.

The numerical values shown in Table 2 are obtained with the model parameters listed in

Table 3. Parameters of the DD2 and GM1L models for the nuclear EOS used in our work [24,36].

Parameter	Units	DD2 Model	GM1L Model
$m_{\sigma }$	GeV	0.5462	0.550
$m_{\omega }$	GeV	0.783	0.783
$m_{\rho }$	GeV	0.763	0.770
$g_{\sigma N}$	−	10.687	9.5722
$g_{\omega N}$	−	13.342	10.618
$g_{\rho N}$	−	3.6269	8.983
${\tilde{b}}_{\sigma }$	−	−	0.0029
${\tilde{c}}_{\sigma }$	−	−	$-0.0011$
$a_{\sigma }$	−	1.3576	−
$b_{\sigma }$	−	0.6344	−
$c_{\sigma }$	−	1.0054	−
$d_{\sigma }$	−	0.5758	−
$a_{\omega }$	−	1.3697	−
$b_{\omega }$	−	0.4964	−
$c_{\omega }$	−	0.8177	−
$d_{\omega }$	−	0.6384	−
$a_{\rho }$	−	0.5189	0.3898

[22,36].

In

Table 4. Particle compositions and coupling constants investigated with the GM1L model.

Model	p	n	$\mathbf{\Sigma }$	$\mathbf{\Lambda }$	$\mathbf{\Xi }$	$\mathbf{\Delta }$	e${}^{-}$	${\mathit{\mu }}^{-}$	${\mathit{g}}_{\mathit{\sigma }\mathit{H}}$	${\mathit{g}}_{\mathit{\omega }\mathit{H}}$	$\frac{{\mathit{g}}_{\mathit{\rho }\mathit{H}}}{{\mathit{g}}_{\mathit{\rho }\mathit{N}}}$	$\frac{{\mathit{g}}_{\mathit{\sigma }\mathbf{\Delta }}}{{\mathit{g}}_{\mathit{\sigma }\mathit{N}}}$	$\frac{{\mathit{g}}_{\mathit{\omega }\mathbf{\Delta }}}{{\mathit{g}}_{\mathit{\omega }\mathit{N}}}$	$\frac{{\mathit{g}}_{\mathit{\rho }\mathbf{\Delta }}}{{\mathit{g}}_{\mathit{\rho }\mathit{N}}}$
GM1L_b	•	•	•	•	•	•	•	•	HPf	SU(3)	1	1.2	1.2	1.2
GM1L_c	•	•	−	−	−	−	•	•	−	−	−	−	−	−
GM1L_f	•	•	•	•	•	•	•	•	HPf	SU(3)	1	0.8	0.8	0.8

we show the particle compositions and the coupling schemes that were chosen to fix the coupling constants of the GM1L model. Models GM1L_b and GM1L_f include all particles of the baryon octet plus all four electrically charged states of the $\Delta \left(1232\right)$ isobar, while model GM1L_c account only for protons and neutrons in chemical equilibrium with each other. Hypernuclear potential fits (HPf) were used to determine the $\sigma$-hyperon ($g_{\sigma H}$) coupling constants. The $\omega$-hyperon ($g_{\omega H}$) coupling values were chosen on the basis of a broken SU(3) flavor symmetry [51]. For the $\rho$ meson, we have chosen the same coupling value as for the $\rho$-nucleon ($g_{\rho N}$) coupling, and the relative coupling strength of the $\Delta \left(1232\right)$ isobar varies between 0.8 and 1.2. Values larger than 1.2 are excluded from our coupling scheme, since the maximum mass of a neutron star calculated for this case falls below the required $2M_{\odot }$ constraint set by the masses of PSR J1614–2230 and PSR J0343+0432 which are $M=1.928\pm 0.017$ M${}_{\odot }$ [52,53] and $M=2.01\pm 0.04$ M${}_{\odot }$ [54], respectively. A detailed discussion of these coupling scenarios can be found in Refs. [22,33].

4. Stellar Structure Equations

4.1. Non-Rotating Compact Stellar Objects

The treatment of compact stellar objects is relatively straightforward as long as these objects are spherically symmetric [55,56], which is assumed here. The components of the element

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }$$

(58)

are then given by

$$d{s}^2=-{\mathrm{e}}^{2\nu \left(r\right)}d{t}^2+{\mathrm{e}}^{2\lambda \left(r\right)}d{r}^2+r^{2}d{\theta }^2+r^2{sin}^2\theta d{\varphi }^2,$$

(59)

with the metric functions $\nu \left(r\right)$ and $\lambda \left(r\right)$ to be determined from Einstein’s field equation (we use the geometric unit system where $G=c=1)$,

$$R^{\kappa \sigma }-\frac{1}{2}{g}^{\kappa \sigma }R=8\pi T^{\kappa \sigma }(\epsilon ,P\left(\epsilon \right)),$$

(60)

and the conservation of energy-momentum

$${\nabla }_{\kappa }{T}^{\kappa \sigma }=0.$$

(61)

The quantity ${\nabla }_{\kappa }$ in Equation (61) denotes the covariant derivative of the energy-momentum tensor given by ${\nabla }_{\kappa }{T}^{\kappa \sigma }={\partial }_{\kappa }{T}^{\kappa \sigma }+{\mathrm{\Gamma }}_{\kappa \mu }^{\kappa }{T}^{\mu \sigma }+{\mathrm{\Gamma }}_{\kappa \mu }^{\sigma }{T}^{\kappa \mu }$. Treating the stellar matter as a perfect fluid leads for $T^{\mu \nu }$ to

$$T^{\mu \nu }(\epsilon ,P\left(\epsilon \right))=u^{\mu }{u}^{\nu }(\epsilon +P\left(\epsilon \right))+g^{\mu \nu }P\left(\epsilon \right).$$

(62)

The four velocities $u^{\mu }$ and $u^{\nu }$ in (62) are give by

$$u^{\mu }\equiv \frac{d{x}^{\mu }}{d\tau },u^{\nu }\equiv \frac{d{x}^{\nu }}{d\tau }.$$

(63)

The covariant components of the metric tensor in Equation (59) read

$$g_{tt}=-{\mathrm{e}}^{2\nu },g_{rr}={\mathrm{e}}^{2\lambda },g_{\theta \theta }=r^2,g_{\varphi \varphi }=r^2{sin}^2\theta .$$

(64)

Due to the underlying symmetries, the metric functions are only dependent on the radial distance r measured from the star’s origin.

The Einstein Equation (60) and the conservation of energy-momentum (61) are used in deriving the Tolman-Oppenheimer-Volkoff (TOV) equation for a spherically symmetric, non-rotating mass distribution (i.e., a compact star). This equation is given by [15,57,58]

$$\frac{dP}{dr}=-\frac{(\epsilon \left(r\right)+P\left(r\right))(m\left(r\right)+4\pi r^{3}P\left(r\right))}{r^2\left(1-\frac{2m\left(r\right)}{r}\right)},$$

(65)

where $P\left(r\right)$ and $ϵ\left(r\right)$ are the pressure and energy density, respectively, of matter at a radial distance r from the star’s center. The TOV equation is solved by specifying the energy density at the center of the star, which, via the EoS, fixes the pressure at the center of the object. The gravitational mass contained in a sphere of radius r is obtained from $dm\left(r\right)=4\pi r^2\epsilon \left(r\right)$, from which the gravitational mass of a star follows as

$$M=4\pi {\int }_0^{R}dr{r}^2\epsilon \left(r\right).$$

(66)

The TOV equation is solved for a range of central energy densities that results in a mass-radius curve that describes a given family of compact stars. As a by-product, the TOV equation also provides information about the baryon number and energy density profiles of compact stars, which are needed to analyze GRR driven instabilities, as will be discussed in Section 5.

In Figure 4, we show the gravitational masses of compact stars as a function of central energy density computed for the EoS models described in Section 2 and Section 3. By construction, the ACB4 model for the EoS which comprises a piece wise polytropic representation, leads to two stable (a neutron star as well as a mass-twin star) branches of compact stars. In contrast, the relativistic mean-field models for the EoS, GM1L$\_b$, GM1L$\_c$, GM1L$\_f$, and DD2, produce only one stable branch of compact (neutron) stars. The neutron stars of the ACB4 model are up to around four times denser than SNM. The mass-twin stars have central densities that are between around five to nine times higher than the density of SNM. The radii of all these stars are shown in Figure 5. As can be seen, the maximum-mass neutron stars obtained with the mean-field models are in the narrow range between 11.5 km and 11.8 km. In comparison, the radius of the most massive twin star is around 11 km, while the radius of the most mass neutron stars of the stellar ACB4 sequence is 13.9 km.

4.2. Rotating Compact Stellar Objects

It is considerably more complicated to determine the properties of rotating compact objects compared to their non-rotating counterparts. The rotating object experiences rotational deformation where there is flattening at the pole and a radial blowup in the equatorial direction. These complications lead to a dependence on the polar coordinate, $\theta$, and the radial coordinate, r, in the star’s metric. The star’s rotational stability against gravity allows it to have more mass than its non-rotating counterpart, affecting the geometry of spacetime and causes the star’s metric to become dependent on the star’s rotational frequency. The general relativistic effect of the dragging of local inertial frames adds a non-diagonal term, $g^{t\varphi }$, to the metric tensor. These additional aspects lead to a line element that has the form [15,59]

$$d{s}^2=-{\mathrm{e}}^{2\nu }d{t}^2+{\mathrm{e}}^{2\psi }{(d\varphi -\omega dt)}^2+{\mathrm{e}}^{2\mu }d{\theta }^2+{\mathrm{e}}^{2\lambda }d{r}^2,$$

(67)

where the quantities $\nu$, $\psi$, $\mu$, and $\lambda$ denote the metric functions which are dependent on the radial coordinate r, the polar angle $\theta$, the star’s angular velocity $\mathrm{\Omega }$, and implicitly on the frame-dragging frequency $\omega$. The quantity $\omega$ describes the angular velocity of the local inertial frames which are dragged along by the rotational motion of a compact object in general relativity theory, known as Lense-Thirring frame dragging [60,61,62,63].

The numerical solution of Einstein’s equations for compact rotating objects is known to be a cumbersome and complicated task [59,64,65]. An alternative treatment, which is easier to implement and has proven to be a practical tool for the construction of models of general relativistic, rapidly rotating neutron stars is Hartle’s perturbative method [66,67]. Within the latter, a perturbative solution of the rotating stellar structure equations, based on the Schwarzschild metric, is developed. The method leads to results that are in good agreement with those obtained by a numerically exact treatment of Einstein’s field equation [68]. This is especially true for the mass increase due to the fast rotation. The values of the Kepler frequencies resulting from the Hartle method differ between 10 and 15% from the numerically exact values, depending on which equation of state has been used for the calculations [68]. We perform our investigations in the framework of this method.

Specifying the conditions for stable rapid rotation is a non-trivial issue in the framework of Einstein’s general relativity [11,15,16]. Gravitational-radiation reaction driven instabilities within a rapidly rotating star may act as a criterion for an upper bound on the rotation frequency due to the onset of emission of gravitational waves carrying away angular momentum from the star. This instability will be discussed in Section 5. A second criterion is given by the onset of mass shedding from the equator of the rapidly spinning compact object, which occurs at the Kepler frequency, $f_{\mathrm{K}}$. If a star spins as fast as $f_{\mathrm{K}}$, mass loss at the equator drives the star out of hydrostatic equilibrium, which leads to a loss of stability. Since no star can spin more rapidly than at the mass-shedding frequency, $f_{\mathrm{K}}$ sets an absolute bound on rapid rotation [15,59]. The general relativistic expression for the Kepler frequency is derived by applying the extremal principle to the circular orbit of a point mass rotating at the star’s equator [15]. From this, the orbital velocity of a comoving observer rotating at the star’s equator (relative to a locally non-rotating observer with zero angular momentum in the $\varphi$-direction) is found as

$$V_{\pm }=\frac{\partial \omega /\partial r}{2\partial \psi /\partial r}{\mathrm{e}}^{\psi -\nu }\pm \sqrt{{\left(\frac{\partial \omega /\partial r}{2\partial \psi /\partial r}\right)}^2{\mathrm{e}}^{2(\psi -\nu )}+\frac{\partial \nu /\partial r}{\partial \psi /\partial r}}.$$

(68)

The general relativistic Kepler frequency is then given by

$$f_{\mathrm{K}}=\frac{1}{2\pi }\left({\mathrm{e}}^{\nu -\psi }{V}_{+}+\omega \right),$$

(69)

where $\nu$, $\psi$, and $\omega$ denote the metric functions and the frame-dragging frequency shown in the line element (67). The quantity $V_{+}$ is the equatorial velocity in the rotational “forward” direction. Equations (68) and (69) are to be computed self-consistently at the star’s equator via a recursive numerical method that solves Einstein’s field equation for a given EoS [15,59]. A guess value for $f_{\mathrm{K}}$ is given for the first iteration of Einstein’s field equations for the metric functions and the frame dragging frequency. Sufficiently accurate self-consistent solutions are typically obtained after 10 to 15 steps of iteration of Einstein’s equations.

In Figure 6 we show the Kepler periods $P_{\mathrm{K}}(=1/f_{\mathrm{K}})$ obtained from Equation (69) for compact stars computed with the models of the EoS of this work. The mass-twin stars of the ACB4 EoS have the smallest rotational periods since these stars have the smallest radii of all the stellar models. The gray rectangle in the upper part of the figure marks the location of pulsars assumed to have masses between 1.3 and $1.6M_{\odot }$ and rotation periods greater than 1.4 ms (720 Hz). The latter corresponds to pulsar PSR J1748-2446ad, which is the most rapidly spinning neutron star observed to date [69]. As can be seen from this figure, all pulsars observed to date rotate well above the Kepler periods that follow Einstein’s field equations for the EoS of this work.

Figure 7 shows the gravitational masses of non-rotating as well as rotating compact stars. The latter rotate at their respective Kepler frequencies $f_{\mathrm{K}}$.

The lines marked 1′–1, 2′–2, 3′–3, and 4′–4 show the evolutionary paths of isolated pulsars (whose baryon numbers are constant) spinning down from high to low frequencies. This leads to rotation-driven re-population effects and could even trigger phase transitions to new types of matter in the stellar cores. For more details, we refer to Ref. [31].

5. Gravitational-Radiation Reaction Driven Instabilities

Neutron stars can emit gravitational waves due to the instabilities driven by their rapid rotation. This gravitational-radiation reaction offers insight into the effects of different dissipation mechanisms within the star, like viscosity, and insight into the possible speed limit on this compact object’s rate of rotation dependent on its EoS. Most importantly, the GRR drives the f-mode instability within these rapidly rotating objects. These modes propagate in the opposite direction to the star’s rotation and have no radial nodes in the non-rotating limit of uniformly rotating polytropes [70]. The root of Equation (3) determines the maximum stable rotational frequency of a compact star. The dissipation timescales in the equation all depend on m, so they can be written as ${\tau }_{\mathrm{GR},m}$, ${\tau }_{\eta ,m}$ and ${\tau }_{\zeta ,m}$.

The dependence of the eigenfrequency $\omega$, which is not to be confused with the frame dragging frequency, on the star’s angular velocity $\mathrm{\Omega }$ is expressed in terms of the function [11]

$${\alpha }_m\left(\mathrm{\Omega }\right)=\frac{\omega \left(\mathrm{\Omega }\right)-m\mathrm{\Omega }}{\omega \left(0\right)}.$$

(70)

The ${\alpha }_m\left(\mathrm{\Omega }\right)$ are very slowly varying functions of $\mathrm{\Omega }$ with $a_m\approx 1$ over the entire range of velocities [11]. The angular velocity dependence of the damping times are expressed as dimensionless functions,

$$\begin{array}{ccc}& \gamma \left(\mathrm{\Omega }\right)& \hfill =\frac{{\omega }_m\left(\mathrm{\Omega }\right)}{{\omega }_m\left(0\right)}{\left(\frac{{\tau }_{\eta }\left(0\right)}{{\tau }_{\mathrm{GR}}\left(0\right)}\frac{{\tau }_{\mathrm{GR}}\left(\mathrm{\Omega }\right)}{{\tau }_{\eta }\left(\mathrm{\Omega }\right)}\right)}^{1/(2m+1)},\hfill \end{array}$$

(71)

$$\begin{array}{ccc}& {ϵ}_m\left(\mathrm{\Omega }\right)& \hfill =\frac{{\tau }_{\zeta }\left(0\right)}{{\tau }_{\eta }\left(0\right)}\frac{{\tau }_{\eta }\left(0\right)}{{\tau }_{\zeta }\left(0\right)}.\hfill \end{array}$$

(72)

These functions are also very slowly varying except for the highest angular velocities $\mathrm{\Omega }$. The functions are difficult to determine, so they are substituted with their corresponding Maclaurin spheroid functions [11,71], where ${\alpha }_m\left(\mathrm{\Omega }\right)$ and ${\gamma }_m\left(\mathrm{\Omega }\right)$ are calculated by courtesy of [72].

The critical rotational stellar frequency at which a given mode becomes unstable is given by the transformation of Equation (70) using Equations (71) and (72) [11],

$$\begin{array}{cc}& {\mathrm{\Omega }}_{\mathrm{GRR}}\equiv {\mathrm{\Omega }}_m=\frac{{\omega }_m\left(0\right)}{m}\\ & \times \left\{{\alpha }_m\left({\mathrm{\Omega }}_m\right)+{\gamma }_m\left({\mathrm{\Omega }}_m\right){\left(\frac{{\tau }_{\mathrm{GR},m}\left(0\right)}{{\tau }_{\eta ,m}\left(0\right)}\frac{{\tau }_{\mathrm{GR},m}\left(0\right)}{{\tau }_{\zeta ,m}\left(0\right)}{ϵ}_m\left({\mathrm{\Omega }}_m\right)\right)}^{1/(2m+1)}\right\},\end{array}$$

(73)

where

$${\omega }_m\left(0\right)=\sqrt{\frac{2m(m-1)}{2m+1}\frac{GM}{R^3}}$$

(74)

is the frequency of the vibrational mode in a non-rotating star. The individual damping times are given by [11]

$$\begin{array}{ccc}\hfill \frac{1}{{\tau }_{\zeta }}\hfill & =& \hfill \frac{1}{2E}\int d^3\mathit{x}\zeta \delta \sigma \delta {\sigma }^{\ast },\hfill \end{array}$$

(75)

$$\begin{array}{ccc}\hfill \frac{1}{{\tau }_{\eta }}\hfill & =& \hfill \frac{1}{E}\int d^3\mathit{x}\eta \delta {\sigma }^{ij}\delta {\sigma }_{ij}^{\ast },\hfill \end{array}$$

(76)

$$\begin{array}{ccc}\hfill \frac{1}{{\tau }_{\mathrm{GR}}}\hfill & =& \hfill \frac{\widehat{\omega }}{2E}\sum _{l=l_{\min }}^{\infty }{N}_l{\omega }^{2l+l}\delta D_l^m\delta D_l^{\ast m}.\hfill \end{array}$$

(77)

The details for the derivations of the above equations can be found in Chapter 16 in [15]. The timescale for gravitational-radiation reaction ${\tau }_{\mathrm{GR},m}$ [73] and the dissipation timescale ${\tau }_{\eta ,m}$ related to the shear viscosity, denoted by $\eta$, due to particle collisions in the dense neutron star core is given by [14]

$$\begin{array}{ccc}\hfill \frac{1}{{\tau }_{\mathrm{GR}}}\hfill & =& \hfill \frac{-32\pi G{\mathrm{\Omega }}^{2m+2}}{c^{2m+3}}\frac{{(m-1)}^{2m}}{{[(2m+1)!!]}^2}{\left(\frac{m+2}{m+1}\right)}^{2m+2}{\int }_0^R\rho r^{2m+2}dr,\hfill \end{array}$$

(78)

$$\begin{array}{ccc}\hfill \frac{1}{{\tau }_{\eta }}\hfill & =& \hfill (m-1)(2m+1){\int }_0^R\eta \left(r\right)r^{2m}dr{\left[{\int }_0^R\rho r^{2m+2}dr\right]}^{-1},\hfill \end{array}$$

(79)

where the expression for the shear viscosity $\eta$ is given by [74,75]

$$\eta \left(\rho \right)=347{\rho }^{9/4}{T}^{-2}.$$

(80)

The dissipation timescale ${\tau }_{\zeta ,m}$ related to the bulk viscosity is given by [14]

$$\frac{1}{{\tau }_{\zeta ,m}}=\frac{1}{2E}\int \zeta \frac{\delta \rho \delta {\rho }^{\ast }}{{\rho }^2}{\mathrm{d}}^{3}x,$$

(81)

where the Eulerian density perturbations $\delta \rho$ are given by [14]

$$\frac{\delta \rho }{\rho }=\alpha R^2{\mathrm{\Omega }}^2\frac{d\rho }{dp}\left[\frac{2l}{2l+1}\sqrt{\frac{l}{l+1}}{\left(\frac{r}{R}\right)}^{l+1}+\delta \mathrm{\Psi }\left(r\right)\right]Y_{l+1l}{\mathrm{e}}^{i\omega t}.$$

(82)

The quantity $\delta \mathrm{\Psi }\left(r\right)$ in the above equation is proportional to the gravitational potential and is given by [14]

$$\begin{array}{c}\hfill \frac{d^2\delta \mathrm{\Psi }}{d{r}^2}+\frac{2}{r}\frac{d\delta \mathrm{\Psi }}{dr}+\left[4\pi G\rho \frac{d\rho }{dp}-\frac{(l+1)(l+2)}{r^2}\right]\delta \mathrm{\Psi }=-\frac{8\pi Gl}{2l+1}\sqrt{\frac{l}{l+1}}\rho \frac{d\rho }{dp}{\left(\frac{r}{R}\right)}^{l+1}.\hfill \end{array}$$

(83)

The lowest order expression for the energy of a given l-mode reads

$$\frac{1}{E}=2{\alpha }^{-2}{\mathrm{\Omega }}^{-2}{R}^{2l-2}{\left[{\int }_0^R\rho r^{2l+2}dr\right]}^{-1},$$

(84)

and the bulk viscosity $\zeta$ is given by [76]

$$\zeta \left(\rho \right)=6.0\times {10}^{-59}{\rho }^2{T}^6.$$

(85)

Finally, the bulk viscosity timescale for a given stellar mass M, radius R, density profile $\rho$ and pressure P are obtained from [14]

$$\begin{array}{c}\hfill \frac{1}{{\tau }_{\zeta ,m}}={\mathrm{\Omega }}^2{R}^{2l+2}{\left[{\int }_0^R\rho r^{2l+2}dr\right]}^{-1}\times \int \zeta {\left[\frac{d\rho }{dp}\left(\frac{2l}{2l+1}\sqrt{\frac{l}{l+1}}{\left(\frac{r}{R}\right)}^{l+1}+\delta \mathrm{\Psi }\left(r\right)\right)Y_{l+1l}{e}^{i\omega t}\right]}^2{d}^{3}x.\hfill \end{array}$$

(86)

The GRR frequencies ${\mathrm{\Omega }}_m$ from Equation (73) are then converted into the GRR-driven instability periods for cold stars ($T\sim {10}^{10}$ K) with different masses and radii along with the harmonic indices $l=m=2,3,4,5$. Figure 8 illustrates the GRR periods for each model and harmonic index compared to each model’s Kepler period.

The limit on stable rotation for all models is set by the $m=2$ mode because it is excited firstly and occurs at a slower frequency, and thus occurs much earlier, than the Kepler frequency where the instability due to mass shedding begins. For example, Figure 8 illustrates that the set of $m=2$$P_{\mathrm{GRR}}$ ($=2\pi /{\mathrm{\Omega }}_{\mathrm{GRR}}$) periods with masses between $M=1.0M_{\odot }$ and $2.0M_{\odot }$ for DD2 have values starting from $1.285$ ms to $0.930$ ms. The set of calculated Kepler periods with masses between $M=1.0M_{\odot }$ and $2.0M_{\odot }$ for DD2 have values from 1.096 ms to 0.671 ms; the $m=2$ f-mode periods for DD2 are ∼20% to ∼40% higher than the measured Kepler periods.

6. Shear and Bulk Viscosity

The bulk viscosity in neutron star matter is due to pressure and density perturbations that become slightly out of phase since the weak interaction needs a long timescale to reestablish local thermodynamic equilibrium [76]. The shear viscosity arises in neutron star matter due to neutron-neutron scattering, under the condition that the temperature exceeds the superfluid transition temperature [74,75]. The bulk viscosity has a greater effect on the GRR-driven instability periods of hot stars ($T\gtrsim {10}^9$ K) than the shear viscosity due to the bulk viscosity’s dependence on the temperature being $T^6$ compared to the shear viscosity’s dependence, $T^{-2}$, as respectively seen in Equations (85) and (80). The bulk viscosity can regulate the GRR driven instability in rapidly rotating neutron stars, pushing their critical rotational periods toward smaller values, possibly even as small as the Kepler period as it dominates over the shear viscosity [15].

There is a notable difference in value among the set of $P_{\mathrm{GRR}}$ illustrated in Figure 9 compared to Figure 8. Without the bulk viscosity function contributing to the calculation of $P_{\mathrm{GRR}}$, the values become higher for the GRR instability periods. The highest $P_{\mathrm{GRR}}$ value for DD2 with bulk viscosity neglected is $P_{\mathrm{GRR}}=$$1.521$ ms, which is $18\%$ higher compared to the highest $P_{\mathrm{GRR}}$ value for DD2 with bulk viscosity contributing, $P_{\mathrm{GRR}}=$$1.285$ ms. The value $P_{\mathrm{GRR}}=$$1.521$ ms is the GRR instability period that corresponds to the the $m=4$ mode for a $1.0$$M_{\odot }$ star. There is $39\%$ difference between the $m=4$ instability period value and the measured Kepler period for DD2 that corresponds to a $1.0$$M_{\odot }$ star, $P_{\mathrm{GRR}}=$$1.096$ ms, compared to the $17\%$ difference between the Kepler period and the $m=2$, $1.0$ $M_{\odot }$ star instability period ($P_{\mathrm{GRR}}=$$1.285$ ms).

7. Discussion

Neutron stars are compact objects that can rapidly rotate and emit radiation due to such rotation. Not only can these objects emit radio waves as a pulsar, but they can also emit gravitational waves due to gravitational-radiation reaction driven f-mode instabilities once the star spins fast enough. As seen by the results, this phenomenon occurs before the star approaches an unstable rotation rate known as the Kepler frequency, the frequency at which mass shedding at the star’s equator occurs. The $m=2$ mode is excited before the Kepler frequency of the star is reached and the system emits such gravitational waves that carry away angular momentum, effectively slowing down the star’s rotation speed. The f-mode instability limit on rotation is consistent for both relativistic mean field models, DD2 and GM1L, as well as the piece-wise polytropic ACB4 model.

There are ongoing disputes concerning the constraints on the EoS of these compact objects. The data from observed gravitational waves and the data from NICER have discrepancies among them regarding the possible masses and radii a star can have. It is expected that as more research has progressed on gravitational waves and neutron star mergers, and when NICER collects more data from its observations, the constraints between the two camps will cover a tighter range and ultimately overlap.

Further study on the temperature dependence of the bulk viscosity and shear viscosity by calculating the EoS of each model with varying temperatures (from 1 MeV to around 50 MeV) to retrieve new GRR driven f-mode instabilities would be beneficial as it would be a practical replication of real compact stars, considering, in reality, most of them vary in temperature. It will also give insight on whether these GRR frequencies are a viable limit for these rapidly rotating objects, or if the contribution of the bulk viscosity on the neutron star matter will require the Kepler frequency to be the set limit on rotation due to the star having to spin faster in order to excite one of the modes—the value of $P_{\mathrm{GRR}}$ decreases below the value of the Kepler period.

We note that compact stars with sharp phase transitions in their interiors can have quasi-normal modes (see, for example, Refs. [78,79,80,81,82,83] and references therein). These modes could be clear indicators of the presence of quark matter in the inner cores of compact objects. Generally, g-modes associated with sharp hadron-quark phase transitions in the inner cores of cold hybrid stars have frequencies ranging from ∼500 Hz to ∼1.5 kHz, and have always lower frequencies than the fundamental mode.

The secular instability of a wide variety of low-frequency g-modes, but not including those due to a hadron-quark phase transition, has been estimated in Ref. [84]. The main conclusion of this paper is that the secular instability in rotating hybrid stars sets in when the compact object’s rotation frequency is ∼30% less than the g-mode frequency of the non-rotating configuration. If this were the case, the g-mode instability might set in at rotation rates that are lower than those of the f-mode. Despite this, in Ref. [84] it was also shown that dissipation due to viscosity might completely suppress the instability outside of an extremely narrow temperature window of $\sim {10}^9$ K. In Ref. [85] the secular instability of g-modes produced by the presence of a hadron-quark mixed phase was studied. The main result is that the appearance of a mixed phase produces a sharp increase in the frequency of the g-modes that might prevent the secular instability from setting in. In this sense, more research related to this subject is needed to better understand this matter.

We note that according to the turning point method of Friedman, Ipser, and Paker [86] axisymmetric stellar configurations can exist which are secularly unstable. It would be very interesting to investigate this instability for the stellar models of our work and to compare the results with the f-mode instability of our study. Such a comparison is planned for future work.

Finally, though the r-mode instability is not explored in this study, it is expected to compete with the f-mode instability [12,13]. Determining which of these two modes places a more stringent limit on the stable rotation periods of neutron stars require substantial additional study of the very nonlinear equations describing the r-modes.