Quark Deconfinement in Rotating Neutron Stars

Abstract

In this paper, we use a three flavor non-local Nambu–Jona-Lasinio (NJL) model, an improved effective model of Quantum Chromodynamics (QCD) at low energies, to investigate the existence of deconfined quarks in the cores of neutron stars. Particular emphasis is put on the possible existence of quark matter in the cores of rotating neutron stars (pulsars). In contrast to non-rotating neutron stars, whose particle compositions do not change with time (are frozen in), the type and structure of the matter in the cores of rotating neutron stars depends on the spin frequencies of these stars, which opens up a possible new window on the nature of matter deep in the cores of neutron stars. Our study shows that, depending on mass and rotational frequency, up to around 8% of the mass of a massive neutron star may be in the mixed quark-hadron phase, if the phase transition is treated as a Gibbs transition. We also find that the gravitational mass at which quark deconfinement occurs in rotating neutron stars varies quadratically with spin frequency, which can be fitted by a simple formula.

1. Introduction

Exploring the properties of compressed baryonic matter, or, more generally, strongly interacting matter at high densities and/or temperatures has become a forefront area of modern physics [1,2,3]. Experimentally, the properties of such matter are being probed with the Relativistic Heavy Ion Collider RHIC at Brookhaven and the Large Hadron Collider (LHC) at Cern. Great advances in our understanding of such matter are expected from the next generation of heavy-ion collision experiments at FAIR (Facility for Antiproton and Ion Research, GSI), NICA (Nucloton-based Ion Collider fAcility, JINR), and, at lower energies, from radioactive beam facilities such as FRIB at MSU.

The universe was filled with hot and dense baryonic matter shortly after the Big Bang. Today, such matter is being created in the universe in the final stages of catastrophic stellar events (e.g., core-collapse supernovae, gamma-ray bursts) and exists permanently inside of neutron stars. Depending on mass and rotational frequency, gravity compresses the matter in the core regions of neutron stars to densities that are several times higher than the density of ordinary atomic nuclei [4,5,6,7,8,9,10,11,12]. At such huge densities atoms themselves collapse, and atomic nuclei are squeezed so tightly together that new particle states may appear and novel states of matter, foremost quark matter, may be created. This feature makes neutron stars superb astrophysical laboratories for a wide range of physical studies [4,5,6,7,9,11,13,14]. With observational data accumulating rapidly from both orbiting and ground based observatories spanning the spectrum from X-rays to radio wavelengths, there has never been a more exiting time than today to study neutron stars. The Hubble Space Telescope and X-ray satellites such as Chandra and XMM-Newton, for instance, have proven especially valuable. New astrophysical instruments such as the Five hundred meter Aperture Spherical Telescope (FAST), the square kilometer Array (skA), Fermi Gamma-ray Space Telescope (formerly GLAST), Astrosat, ATHENA (Advanced Telescope for High ENergy Astrophysics), and the Neutron Star Interior Composition Explorer (NICER) promise the discovery of tens of thousands of new neutron stars. In particular, the NICER mission is expected to open up a window on the inner workings of neutron stars, as this mission is dedicated to the study of extraordinary gravitational, electromagnetic, and nuclear physics environments embodied by neutron stars.

In this paper, we use a non-local extension of the SU(3) Nambu–Jona-Lasinio (NJL) model to investigate the possible existence of deconfined quarks in the cores of neutron stars. Investigations of dense neutron star matter based on the non-local SU(2) version of the NJL model have been carried out in Refs. [15,16].

The possible existence of such matter inside of neutron stars has already been suggested in the 1970s [17] but has remained an open issue ever since, as QCD cannot be solved for dense astrophysical matter with the tools usually employed by relativistic quantum field theories. We put particular emphasis on the possible existence of quark matter inside of rotating neutron stars, known as pulsars. In contrast to non-rotating neutron stars, whose core compositions do not change with time (are frozen in), the type and structure of the matter in the cores of rotating neutron stars depends on the spin frequencies of these stars, which could give rise to observable astrophysical signals of quark deconfinement [18,19,20,21,22].

A histogram of the rotational frequencies of 2510 observed pulsars is shown in Figure 1. This number is expected to increase dramatically once new instruments such as FAST and skA come into operation.

The paper is organized as follows. In Section 2 we introduce the nuclear equations of state (EoS) that were used to model the properties of (non) rotating neutron stars. The treatment of rotating neutron stars in the framework of general relativity theory is discussed in Section 3. The results of our investigations are presented in Section 4. A brief discussion of the results in provided in Section 5.

2. Models for the Nuclear Equations of State

2.1. Hadronic Matter

In the most primitive conception, the matter in the core of a neutron star is constituted from neutrons. At a slightly more accurate representation, the cores consist of neutrons and protons whose electric charge is balanced by leptons ($\lambda =\{e^{-},{\mu }^{-}\}$). Other particles, like hyperons ($B=\{n,p,\mathrm{\Lambda },\mathrm{\Sigma },\mathrm{\Xi }\}$) and the Δ resonance, may be present if the Fermi energies of these particles become large enough so that the existing baryon populations can be rearranged and a lower energy state be reached. To model this hadronic phase, we make use of the relativistic mean-field (RMF) approximation, in which the interactions between baryons are described by the exchange of scalar (σ), vector (ω), and isovector (ρ) mesons [26]. The standard mean-field Lagrangian is given by [4,5,27,28,29,30].

$$\begin{array}{ccc}& & \mathcal{L}=\sum _B{\overline{\psi }}_B\left[\left(i{\gamma }^{\mu }{\partial }_{\mu }-m_B\right)+{\mathrm{\Gamma }}_{\sigma B}-{\gamma }_{\mu }{\mathrm{\Gamma }}_{\omega B}^{\mu }-{\gamma }_{\mu }{\mathrm{\Gamma }}_{\rho B}^{\mu }\mathit{\tau }\right]{\psi }_B+\frac{1}{2}({\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -m_{\sigma }^2{\sigma }^2)-\frac{1}{3}{b}_{\sigma }{m}_n{\left(g_{\sigma B}\sigma \right)}^3\hfill \\ & & -\frac{1}{4}{c}_{\sigma }{\left(g_{\sigma B}\sigma \right)}^4-\frac{1}{4}{\omega }_{\mu \nu }{\omega }^{\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{\rho }}_{\mu \nu }·{\mathit{\rho }}^{\mu \nu }+\sum _{\lambda }{\overline{\psi }}_{\lambda }(i{\gamma }^{\mu }{\partial }_{\mu }-m_{\lambda }){\psi }_{\lambda },\hfill \end{array}$$

(1)

with the meson-baryon vortices given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\sigma B}=g_{\sigma B}\sigma ,{\mathrm{\Gamma }}_{\omega B}^{\mu }=g_{\omega B}{\omega }^{\mu },{\mathrm{\Gamma }}_{\rho B}^{\mu }=g_{\rho B}{\mathit{\rho }}^{\mu }.\end{array}$$

(2)

The σ and ω mesons in Equation (1) are responsible for nuclear binding while the ρ meson is required to obtain the correct value for the empirical symmetry energy. The cubic and quartic σ terms in Equation (1) are necessary, at the relativistic mean-field level, to obtain the empirical incompressibility of nuclear matter [28,29]. The field tensors ${\omega }_{\mu \nu }$ and ${\mathit{\rho }}_{\mu \nu }$ are defined as ${\omega }_{\mu \nu }={\partial }_{\mu }{\omega }_{\nu }-{\partial }_{\nu }{\omega }_{\mu }$ and ${\mathit{\rho }}_{\mu \nu }={\partial }_{\mu }{\mathit{\rho }}_{\nu }-{\partial }_{\nu }{\mathit{\rho }}_{\mu }$.

The parameters (i.e., coupling constants) of the theory must reproduce the bulk properties of infinite nuclear matter at saturation density, ${\rho }_0=0.16{\mathrm{fm}}^{-3}$. These are the binding energy $E/N$, effective nucleon mass $m_N^{*}/m_N$, nuclear incompressibility K, and the symmetry energy $a_{\mathrm{sy}}$ and its density derivative L. Of the six, the values of K, $a_{\mathrm{sy}}$, and L carry some uncertainty. The K value is believed to lie in the range between about 220 and 260 MeV [31], or between 250 and 315 MeV, as recently suggested in Reference [32]. The values for $a_{\mathrm{sy}}$ and L are in the ranges of 29 to 35 MeV and 43 to 70 MeV, respectively [33,34,35]. We have chosen two parameter sets, denoted GM1 and DD2, which approximately cover the uncertainties in the nuclear matter properties just above. They are listed in

Table 1. Properties of infinite nuclear matter at saturation density computed for parameter sets GM1 [36] and DD2 [37]. Shown are the saturation density ${\rho }_0$, energy per nucleon $E/N$, nuclear incompressibility K, effective nucleon mass $m_N^{*}/m_N$, asymmetry energy $a_{sy}$, and the density derivative of the symmetry energy, L.

Nuclear Matter Property	Units	GM1	DD2
${\rho }_0$	fm${}^{-3}$	0.153	0.149
E/N	MeV	$-16.3$	$-16.0$
K	MeV	300	243
$m_N^{*}/m_N$		0.70	0.56
$a_{sy}$	MeV	32.5	32.7
L	MeV	91.96	55.04

.

The coupling constants associated with GM1 are $g_{\sigma B}=9.572$, $g_{\omega B}=10.618$, $g_{\rho B}=8.198$, $b_{\sigma }=0.002936$, and $c_{\sigma }=-0.00107$ [36].

The field equations for the baryon fields ${\psi }_B$, which follow from Equation (1), are given by [4,5,8,27].

$$\begin{array}{c}\hfill \left(i{\gamma }^{\mu }{\partial }_{\mu }-m_B\right){\psi }_B=-g_{\sigma B}\sigma {\psi }_B+g_{\omega B}{\gamma }^{\mu }{\omega }_{\mu }{\psi }_B+g_{\rho B}{\gamma }^{\mu }\mathit{\tau }\mathbf{·}{\mathit{\rho }}_{\mu }{\psi }_B.\end{array}$$

(3)

The meson fields in Equation (3) are solutions of the following field equations [4,5,8,27].

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

(4)

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

(5)

$$\begin{array}{ccc}\hfill {\partial }^{\mu }{\mathit{\rho }}_{\mu \nu }+m_{\rho }^2{\mathit{\rho }}_{\nu }& =& \sum _B{g}_{\rho B}{\overline{\psi }}_B\mathit{\tau }{\gamma }_{\nu }{\psi }_B.\hfill \end{array}$$

(6)

In the standard RMF limit, the meson fields of Equations (4)–(6) simplify to [4,5,8,27].

$$\begin{array}{ccc}\hfill m_{\sigma }^2\overline{\sigma }& =& \sum _B{g}_{\sigma B}\frac{2J_B+1}{2{\pi }^2}{\int }_0^{k_B}\frac{m_B^{*}\left(\sigma \right)}{\sqrt{k^2+m_B^{*2}\left(\sigma \right)}}{k}^{2}dk-b{m}_N{g}_{\sigma }{\left(g_{\sigma }\overline{\sigma }\right)}^2-c{g}_{\sigma }{\left(g_{\sigma }\overline{\sigma }\right)}^3,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\overline{\omega }}_0& =& \sum _B\frac{g_{\omega B}}{m_{\omega }^2}{\rho }_B,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\overline{\rho }}_{03}& =& \sum _B\frac{g_{\rho B}}{m_{\rho }^2}{I}_{3B}{\rho }_B,\hfill \end{array}$$

(9)

where $\overline{\sigma }$, $\overline{\omega }$, and ${\overline{\rho }}_{03}$ denote the mean-field limits of σ, ω, and $\overline{\rho }$, respectively, and the effective baryon masses are given by $m_B^{*}\left(\overline{\sigma }\right)=m_B-g_{\sigma B}\overline{\sigma }$. $J_B$ and $I_{3B}$ denote the spin respectively isospin of a baryon of type B, and ${\rho }_B$ stands for the number density of baryon B.

In addition to the standard RMF theory discussed just above, we also consider a lagrangian where the meson-baryon vortices ${\mathrm{\Gamma }}_{MB}$ (where $M=\sigma ,\omega ,\rho$) of Equation (2) are no longer constant but rather depend on density [38,39,40]. In that case the values of the vortices are derived from relativistic Dirac-Brueckner calculations of nuclear matter, which use one-boson-exchange interactions as an input. A characteristic feature of the density-dependent theory is the occurrence of rearrangement terms in the expression for the baryon chemical potential, which leads to a more complex condition for chemical equilibrium compared to the standard RMF approximation [40]. The parameter set of the density-dependent (DD) treatment adopted in this paper is denoted DD2 ($G_V=0$), where a vanishing vector coupling constant $G_V=0$ among quarks has been chosen [37]. For DD2, the coupling constants at saturation density are $g_{\sigma N}=10.687$, $g_{\omega N}=13.342$, and $g_{\rho N}=3.627$ [37], which lead to the saturation properties of infinite nuclear matter shown in Table 1.

The equation of state of the standard mean-field treatment (for the DD formalism, see References [39,40,41,42]) is obtained by solving Equations (7)–(9) together with the charge conservation conditions (baryonic, electric) given by [4,5,27].

$$\begin{array}{c}\hfill {\rho }_b-\sum _B{\rho }_B=0,\sum _B{\rho }_B{q}_B+\sum _{\lambda }{\rho }_{\lambda }{q}_{\lambda }=0,\end{array}$$

(10)

where ${\rho }_b$ is the total baryonic density and $q_B$ and $q_{\lambda }$ are the electric charges of baryons and leptons, respectively. Particles in the hadronic phase are subject to the chemical equilibrium condition

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

(11)

where ${\mu }_B$ is the chemical potential and $b_B$ is the baryon number of baryon B. ${\mu }_n$ and ${\mu }_e$ denote the linearly independent chemical potentials of neutrons and electrons, respectively, which reflect baryon number and electric charge conservation on neutron star matter. New baryon or lepton states are populated when the right side of Equation (11) is greater than the particle’s chemical potential. The baryonic and leptonic number densities (${\rho }_B$ and ${\rho }_{\lambda }$) are both given by ${\rho }_i=(2J_i+1)k_i^3/6{\pi }^2$. The unknowns of the theory are the meson mean-fields (σ, ω, ρ), and the neutron and electron fermi momenta ($k_n$, $k_e$). Finally, the energy density and pressure of the hadronic phase are given by [4,5,27].

$$\begin{array}{ccc}\hfill {ϵ}_H& =& \frac{1}{3}b{m}_N{\left(g_{\sigma N}\overline{\sigma }\right)}^3+\frac{1}{4}c{\left(g_{\sigma N}\overline{\sigma }\right)}^4+\frac{1}{2}{m}_{\sigma }^2{\overline{\sigma }}^2+\frac{1}{2}{m}_{\omega }^2{\overline{\omega }}_0^2+\frac{1}{2}{m}_{\rho }^2{\overline{\rho }}_{03}^2+\sum _B\frac{2J_B+1}{2{\pi }^2}{\int }_0^{k_B}\sqrt{k^2+m_B^{*2}}{k}^{2}dk\hfill \\ & & +\frac{1}{{\pi }^2}\sum _{\lambda }{\int }_0^{k_{\lambda }}\sqrt{k^2+m_{\lambda }^2}{k}^{2}dk,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill p_H& =& -\frac{1}{3}b{m}_N{\left(g_{\sigma N}\overline{\sigma }\right)}^3-\frac{1}{4}c{\left(g_{\sigma N}\overline{\sigma }\right)}^4-\frac{1}{2}{m}_{\sigma }^2{\overline{\sigma }}^2+\frac{1}{2}{m}_{\omega }^2{\overline{\omega }}_0^2+\frac{1}{2}{m}_{\rho }^2{\overline{\rho }}_{03}^2+\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^{*2}}}\hfill \\ & & +\frac{1}{3}\sum _{\lambda }\frac{1}{{\pi }^2}{\int }_0^{k_{\lambda }}\frac{k^{4}dk}{\sqrt{k^2+m_{\lambda }^2}}.\hfill \end{array}$$

(13)

2.2. Deconfined Quark Phase

A popular model widely used to describe deconfined 3-flavor (up, down, strange) quark matter is the Nambu–Jona-Lasinio model [43,44,45,46,47]. Here we use a non-local extension of this model (n3NJL) [48,49], whose effective action is given by

$$\begin{array}{ccc}\hfill S_E& =& \int d^{4}x\left\{\overline{\psi }\left(x\right)\left[-i\partial /+\widehat{m}\right]\psi \left(x\right)-\frac{G_S}{2}\left[j_a^S\left(x\right)j_a^S\left(x\right)+j_a^P\left(x\right)j_a^P\left(x\right)\right]\right.\hfill \\ & & -\frac{H}{4}{T}_{abc}\left[j_a^S\left(x\right)j_b^S\left(x\right)j_c^S\left(x\right)-3j_a^S\left(x\right)j_b^P\left(x\right)j_c^P\left(x\right)\right]\left.-\frac{G_V}{2}\left[j_V^{\mu }\left(x\right)j_V^{\mu }\left(x\right)\right]\right\},\hfill \end{array}$$

(14)

where $\psi \equiv {(u,d,s)}^T$, $\widehat{m}=\mathrm{diag}(m_u,m_d,m_s)$ is the current quark mass matrix, ${\lambda }_a$ ($a=1,...,8$) denote the Gell-Mann matrices–generators of SU(3), and . The currents $j_a^{S,P}\left(x\right)$ and $j_V^{\mu }\left(x\right)$ are given by

$$\begin{array}{ccc}\hfill j_a^S\left(x\right)& =& \int d^{4}z\tilde{g}\left(z\right)\overline{\psi }\left(x+\frac{z}{2}\right){\lambda }_a\psi \left(x-\frac{z}{2}\right),\hfill \\ \hfill j_a^P\left(x\right)& =& \int d^{4}z\tilde{g}\left(z\right)\overline{\psi }\left(x+\frac{z}{2}\right)i{\gamma }_5{\lambda }_a\psi \left(x-\frac{z}{2}\right),\hfill \\ \hfill j_V^{\mu }\left(x\right)& =& \int d^{4}z\tilde{g}\left(z\right)\overline{\psi }\left(x+\frac{z}{2}\right){\gamma }^{\mu }{\lambda }_a\psi \left(x-\frac{z}{2}\right),\hfill \end{array}$$

(15)

where $\tilde{g}\left(z\right)$ is a form factor responsible for the non-local character of the interaction. Finally, the constants $T_{abc}$ in the t’Hooft term accounting for flavor-mixing are defined by

$$T_{abc}=\frac{1}{3!}{ϵ}_{ijk}{ϵ}_{mnl}{\left({\lambda }_a\right)}_{im}{\left({\lambda }_b\right)}_{jn}{\left({\lambda }_c\right)}_{kl}.$$

(16)

The current quark mass $\overline{m}$ of up and down quarks and the coupling constants $G_S$ and H in Equation (14), have been fitted to the pion decay constant, $f_{\pi }$, and meson masses $m_{\pi }$, $m_{\eta }$, and $m_{{\eta }^{\prime }}$, as described in [50,51]. The result of this fit is $\overline{m}=6.2$ MeV, $\mathsf{\Lambda }=706.0$ MeV, $G_S{\mathsf{\Lambda }}^2=15.04$, $H{\mathsf{\Lambda }}^5=-337.71$. The strange quark current mass is treated as a free parameter and was set to $m_s=140.7$ MeV. The strength of the vector interaction $G_V$ is expressed in terms of the strong coupling constant $G_S$. To account for the uncertainty in the theoretical predictions for the ratio $G_V/G_S$, we treat the vector coupling constant as a free parameter [52,53,54], which varies from 0 to $0.09G_S$.

For the mean-field approximation, the thermodynamic potential following from $S_E$ of Equation (14) is given by

$$\begin{array}{ccc}\mathsf{\Omega }& =& -\frac{3}{{\pi }^3}\sum _{f=u,d,s}{\int }_0^{\infty }d{p}_0{\int }_0^{\infty }dp\ln \left\{\left[{\widehat{\omega }}_f^2+M_f^2\left({\omega }_f^2\right)\right]\frac{1}{{\omega }_f^2+m_f^2}\right\}-\frac{3}{{\pi }^2}\sum _{f=u,d,s}{\int }_0^{\sqrt{{\mu }_f^2-m_f^2}}dp{p}^2\times \hfill \\ & & \left[\left({\mu }_f-E_f\right)\theta \left({\mu }_f-m_f\right)\right]-\frac{1}{2}\left[\sum _{f=u,d,s}\left({\overline{\sigma }}_f{\overline{S}}_f+\frac{G_S}{2}{\overline{S}}_f^2\right)+\frac{H}{2}{\overline{S}}_u{\overline{S}}_d{\overline{S}}_s\right]-\sum _{f=u,d,s}\frac{{\overline{\omega }}_f^2}{4G_V},\hfill \end{array}$$

(17)

where ${\overline{\sigma }}_f$, ${\overline{\omega }}_f$, and ${\overline{S}}_f$ are the quark scalar, vector, and auxiliary mean fields, respectively. Moreover, $E_f$ is given by $E_f=\sqrt{{\mathit{p}}^2+m_f^2}$ and we have defined ${\omega }_f^2={(p_0+i{\mu }_f)}^2+{\mathit{p}}^2$. The momentum dependent quark masses are given by $M_f\left({\omega }_f^2\right)=m_f+{\overline{\sigma }}_{f}g\left({\omega }_f^2\right)$. The quantities $g\left({\omega }_f^2\right)=\exp (-{\omega }_f^2/{\mathsf{\Lambda }}^2)$ are the Gaussian form factors, responsible for the nonlocal nature of the interaction among quarks [49]. The auxiliary mean fields are given by

$${\overline{S}}_f=-48{\int }_0^{\infty }d{p}_0{\int }_0^{\infty }\frac{dp}{8{\pi }^3}g\left({\omega }_f^2\right)\frac{M_f\left({\omega }_f^2\right)}{{\widehat{\omega }}^2+M_f^2\left({\omega }_f^2\right)}.$$

(18)

The vector interactions taken into account in the treatment shift the quark chemical potentials as ${\widehat{\mu }}_f={\mu }_f-g\left(w_f^2\right){\overline{\omega }}_f$ and ${\widehat{\omega }}_f^2={(p_0+i{\widehat{\mu }}_f)}^2+p^2$. The scalar and vector mean fields are obtained by minimizing the grand thermodynamic potential, $\partial \mathsf{\Omega }/\partial {\overline{\sigma }}_f=0$ and $\partial \mathsf{\Omega }/\partial {\overline{\omega }}_f=0$. Finally, the quark number densities are obtained from ${\rho }_f=\partial \mathsf{\Omega }/\partial {\mu }_f$.

To determine the equation of state one must solve a nonlinear system of equations for the fields ${\overline{\sigma }}_f$ and ${\overline{\omega }}_f$, and the neutron and electron chemical potentials ${\mu }_n$ and ${\mu }_e$. This system of equations consists of the mean field equations

$${\overline{\sigma }}_i+G_S{\overline{S}}_i+H{\overline{S}}_j{\overline{S}}_k=0,$$

(19)

with cyclic permutations over the quark flavors, ${\overline{\omega }}_f-2G_V\partial \mathsf{\Omega }/\partial {\overline{\omega }}_f=0$, and the baryonic and electric charge conservation equations ${\sum }_{f=u,d,s}{\rho }_f-3{\rho }_b=0$ and ${\sum }_{f=u,d,s}{\rho }_f{q}_f+{\sum }_{\lambda =e^{-},{\mu }^{-}}{\rho }_{\lambda }{q}_{\lambda }=0$, respectively. Finally, the pressure $p_Q$ and energy density ${ϵ}_Q$ are given by

$$p_Q={\mathsf{\Omega }}_0-\mathsf{\Omega },$$

(20)

and

$${ϵ}_Q=-p_Q+\sum _{f=u,d,s}{\rho }_f{\mu }_f+\sum _{\lambda =e^{-},{\mu }^{-}}{\rho }_{\lambda }{\mu }_{\lambda },$$

(21)

where ${\mathsf{\Omega }}_0$ was chosen by the condition that $P_Q$ vanishes at $T=\mu =0$ [48,49].

2.3. Quark-Hadron Mixed Phase

To model the mixed phase region of quarks and hadrons in neutron stars, we use the Gibbs condition for phase equilibrium between hadronic (H) and quark (Q) matter,

$$\begin{array}{c}\hfill P_H({\mu }_n,{\mu }_e,\left\{\varphi \right\})=P_Q({\mu }_n,{\mu }_e),\end{array}$$

(22)

where $P_H$ and $P_Q$ denote the pressures of hadronic matter and quark matter, respectively [55,56]. The quantity $\left\{\varphi \right\}$ in Equation (22) stands collectively for the field variables ($\overline{\sigma }$, $\overline{\omega }$, $\overline{\rho }$) and Fermi momenta ($k_B$, $k_{\lambda }$) that characterize a solution to the equations of confined hadronic matter (Section 2). We use the symbol $\chi \equiv V_Q/V$ to denote the volume proportion of quark matter, $V_Q$, in the unknown volume V. By definition, χ varies between 0 and 1, depending on how much confined hadronic matter has been converted to quark matter. Equation (22) is to be supplemented with the conditions of global baryon charge conservation and global electric charge conservation. The global conservation of baryon charge is expressed as [55,56].

$$\begin{array}{c}\hfill {\rho }_b=\chi {\rho }_Q({\mu }_n,{\mu }_e)+(1-\chi ){\rho }_H({\mu }_n,{\mu }_e,\left\{\varphi \right\}),\end{array}$$

(23)

where ${\rho }_Q$ and ${\rho }_H$ denote the baryon number densities of the quark phase and hadronic phase, respectively. The global neutrality of electric charge is given by [55,56].

$$\begin{array}{c}\hfill 0=\chi q_Q({\mu }_n,{\mu }_e)+(1-\chi )q_H({\mu }_n,{\mu }_e,\left\{\varphi \right\}),\end{array}$$

(24)

with $q_Q$ and $q_H$ denoting the electric charge densities of the quark phase and hadronic phase, respectively. We have chosen global rather than local electric charge neutrality. Local NJL studies carried out for local electric charge neutrality have been reported recently in References [57,58,59,60].

In Figure 2, we show the GM1 and DD2 equations of state (EoS) used in this work to study the quark-hadron composition of rotating neutron stars. The solid dots mark the beginning and the end of the quark-hadron mixed phases for these equations of state. Since the Gibbs condition is used to model the quark-hadron phase transition, pressure varies monotonically with the proportion of the phases in equilibrium. This would not be the case if the Maxwell construction were used to model the phase equilibrium between quarks and hadrons, in which case the pressure throughout the mixed phase is constant. For that reason, the mixed phase would strictly be excluded from neutron stars, but small cores made entirely of quark matter may still be possible in neutron stars close the maximum-mass model (see, for instance, Refrences [61,62], and references therein).

A Maxwell-like phase transition is generally supported by larger surface tensions, σ, of quark matter (see, however, Reference [63]). Values of $\sigma \sim 30{\mathrm{MeV}/\mathrm{fm}}^2$ have recently been suggested in the literature [64,65,66,67,68], but its actual value is an open issue.

The GM1 and DD2 equations of state are compared in Figure 2 with models for the equation of state that have recently been suggested in the literature. ’HLPS’ and ’Neutron matter’ show the constraints on the equation of state established by Hebeler, Lattimer, Pethick, and Schwenk [33,69]. The curves labeled ’EoS I’, ’EoS II’, and ’EoS III’ show the compact star matter equations of state determined by Kurkela et al. [70], which are based on an interpolation between the regimes of low-energy chiral effective theory and high-density perturbative QCD. One sees that the GM1 and DD2 models are well within these limits. The only difference concerns the behavior of the equation of state at sub-nuclear densities (labeled “Neutron matter” in Figure 2), where our models provide slightly more pressure. This, however, does in no way impact the results for the quark-hadron compositions shown in Section 4, because of the large masses of these stars.

Figure 3 shows the gravitational mass versus central neutron stars density (left panel) and gravitational mass versus radius of non-rotating neutron stars for the equations of state discussed in this section. The maximum masses of these neutron stars are between 1.9 and $2.1M_{\odot }$. Stellar rotation, which will be discussed in the next section, increases the masses of these stars to values that are between 2.2 and $2.4M_{\odot }$, while their central densities drops by around 20%.

3. Treatment of Rotating Neutron Stars in General Relativity Theory

The fact that rotation deforms neutron stars, stabilizes them against collapse, and drags along the local inertial frames inside and outside of them so that they co-rotate with the stars, renders the construction of models of rotating neutron stars very complicated. A line element which accounts for these three features has the form [5,71,72,73,74,75]

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

(25)

where ν, ψ, μ and λ denote metric functions and ω is the angular velocity of the local inertial frames. All these quantities depend on the radial coordinate r, the polar angle θ, and implicitly on the star’s angular velocity Ω. The metric functions and the frame dragging frequencies are to be computed from Einstein’s field equation,

$$R^{\kappa \sigma }-\frac{1}{2}R{g}^{\kappa \sigma }=8\pi T^{\kappa \sigma },$$

(26)

where $T^{\kappa \sigma }=T^{\kappa \sigma }(ϵ,P\left(ϵ\right))$ denotes the energy momentum tensor of the stellar matter, whose equation of state is given by $P\left(ϵ\right)$. The other quantities in Equation (26) are the Ricci tensor $R^{\kappa \sigma }$, the curvature scalar R, and the metric tensor, $g^{\kappa \sigma }$. No simple stability criteria are known for rapidly rotating stellar configurations in general relativity. However, an absolute limit on rapid rotation is set by the onset of mass shedding from the equator of a rotating star. The corresponding rotational frequency is known as the Kepler frequency, ${\mathsf{\Omega }}_{\mathrm{K}}$. In classical mechanics, the expression for the Kepler frequency, determined by the equality between the centrifugal force and gravity, is readily obtained as ${\mathsf{\Omega }}_{\mathrm{K}}=\sqrt{M/R^3}$. Its general relativistic counterpart is given by [5,71].

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{\mathrm{K}}=\omega +\frac{{\omega }_{,r}}{2{\psi }_{,r}}+e^{\nu -\psi }\sqrt{\frac{{\nu }_{,r}}{{\psi }_{,r}}+{\left(\frac{{\omega }_{,r}}{2{\psi }_{,r}}{e}^{\psi -\nu }\right)}^2},\end{array}$$

(27)

which is to be evaluated self-consistently at the equator of a rotating neutron star. The Kepler period follows from Equation (27) as $P_{\mathrm{K}}=2\pi /{\mathsf{\Omega }}_{\mathrm{K}}$. For typical neutron star matter equations of state, the Kepler period obtained for $1.4M_{\odot }$ neutron stars is typically around 1 ms. [5,8,71,72]. An exception to this are strange quark matter stars. Since they are self-bound, they tend to possess smaller radii than conventional neutron stars, which are bound by gravity only. Because of their smaller radii, strange stars can withstand mass shedding from the equator down to periods of around 0.5 ms [76,77].

A mass increase of up to $\sim 20$% is typical for rotation at ${\mathsf{\Omega }}_{\mathrm{K}}$. Because of rotation, the equatorial radii increase by several kilometers, while the polar radii become smaller by several kilometers. The ratio between both radii is around 2/3, except for rotation close to the Kepler frequency. The most rapidly rotating, currently known neutron star is pulsar PSR J1748-2446ad, which rotates at a period of 1.39 ms (719 Hz) [25], well below the Kepler frequency for most neutron star equations of state. Examples of other rapidly rotating neutron stars are PSRs B1937+21 [78] and B1957+20 [79], whose rotational periods are 1.58 ms (633 Hz) and 1.61 ms (621 Hz), respectively.

The density change in the core of a neutron star whose frequency varies from $0\le \mathsf{\Omega }\le {\mathsf{\Omega }}_{\mathrm{K}}$ can be as large as 60% [5,8]. This suggests that rotation may drive phase transitions or cause significant compositional changes of the matter in the cores of neutron stars [4,5,8].

4. Results

For each equation of state, a numerical code, based on Hartle’s rotation formalism [80,81,82], was used to predict the percentage of the mass made up by deconfined quark matter over a wide range of frequencies and gravitational masses. Figure 4 and Figure 5 show the composition of rotating neutron stars based on the lagrangian of Equation (1). As can be seen, the quark-hadron mixed phase as well as several different hyperon species are successively spun out of the neutron star if the rotation rate increases toward the Kepler frequency. Non-rotating neutron stars posses the most complex compositions, since they are the most dense members of the rotational sequence. The compositions shown in Figure 5 are snapshots taken from movies (see ”Supplementary Materials” at the end of the paper) showing the entire rotational evolution of this neutron star all the way from zero frequency to ${\mathsf{\Omega }}_{\mathrm{K}}$.

A heat map showing the quark-hadron content of rotating neutron stars computed for the GM1 ($G_V=0.05G_S$) equation of state is shown in Figure 6. As can be seen in this figure, up to 8% of the total gravitational mass of these neutron stars exists in the form a mixed quark-hadron phase. Lines of constant baryon number are also depicted on the figures as white lines, labeled with the logarithm of the star’s baryon number. These lines were included to give a sense of the path that a secluded neutron star would be expected to take through as it spins down.

By parsing out the maximum frequency where deconfined quark matter is expected to exist at the center of the neutron star for a given mass it is possible to get a curve through the gravitational mass–frequency diagram for the threshold above which one can expect to find deconfined quarks. These threshold frequencies for each mass were fit to determine a (quadratic) function for the curve. The fit equation, depicted as a dashed white line in Figure 6, was found to have the form

$$M\left(\mathsf{\Omega }\right)=a{\mathsf{\Omega }}^2+c,$$

(28)

where M is the neutron star’s gravitational mass in solar masses, Ω its rotational frequency, and a and c are parameters determined through fitting. The values for a and c can be found in

Table 2. Parameters for the empirical deconfinement threshold curve for each equation of state (EoS) with the form shown in Equation (28).

EoS	a (${\mathit{M}}_{\odot }{\mathit{s}}^2$)	c (${\mathit{M}}_{\odot }$)
GM1 ($G_V=0.05G_S$)	$2.48\times {10}^{-7}$	$1.91$
GM1 ($G_V=0$)	$2.75\times {10}^{-7}$	$1.71$
DD2 ($G_V=0$)	$2.56\times {10}^{-7}$	$1.89$

.

Figure 6 allows one to estimate the amount of quark-hadron matter that may exist in the cores of neutron stars that have both a measured frequency and mass, as illustrated for pulsars PSR J1614-2230 ($M=1.928\pm 0.017M_{\odot }$, rotational frequency $f=318$ Hz) [83,84] and PSR J0348+0432 ($M=2.01\pm 0.04M_{\odot }$, $f=26$ Hz) [85,86] in Figure 6. According to this calculation, up to around 7% of the mass of PSR J0348+0432 could be in the mixed quark-hadron phase, while the core of the more rapidly rotating pulsars PSR J1614-2230 may be hovering right at the quark deconfinement density.

5. Discussion and Summary

The type and structure of the matter in the cores of rotating neutron stars (pulsars) depends on the spin frequencies of these stars [5,8,87], which opens up a new window on the nature of matter deep in their cores. We find that, depending on mass and rotational frequency, up to around 8% of the mass of massive neutron stars may be in the mixed quark-hadron phase, if the quark-hadron phase transition is Gibbs-like. Examples of such stars are pulsars PSR J1614-2230 with a gravitational mass of $1.928\pm 0.017M_{\odot }$ [84] and PSR J0348+0432 with a mass of $2.01\pm 0.04M_{\odot }$ [85,86] (Figure 6). Pure quark matter in the centers of neutron stars is not obtained for any of the models for the nuclear equation of state studied in this work. We also find that the gravitational mass at which quark deconfinement occurs in rotating neutron stars varies quadratically with spin frequency, which can be fitted by a simple quadratic formula.

Our view of the interior composition of pulsars has changed dramatically since their first discovery almost 50 years ago. It has also become clear during that time period that all the ambient conditions that characterize pulsars tend to the extreme, making pulsars almost ideal astrophysical laboratories for a broad range of physical studies. Owing to the unprecedented wealth of high-quality data on pulsars provided by radio telescopes, X-ray satellites–and soon the latest generation of gravitational-wave detectors–it seems within reach to decipher the inner workings of pulsars, and to explore the phase diagram of cold and ultra-dense hadronic matter from astrophysics.