Massive Neutron Stars with a Color Superconducting Quark Matter Core

TL;DR

This paper develops an equation of state for neutron star matter using the NJL model, demonstrating that neutron stars with color superconducting quark cores can exceed 2 solar masses, aligning with recent observations.

Contribution

It shows that the same NJL model parameters can describe both nuclear and quark matter phases, predicting massive neutron stars with quark cores.

Findings

Maximum neutron star mass exceeds 2.01 solar masses.

Equation of state aligns with gravitational wave data.

Color superconducting quark matter can be stable at high densities.

Abstract

We construct the equation of state for high density neutron star matter at zero temperature using the two-flavor Nambu--Jona-Lasinio (NJL) model as an effective theory of QCD. We build nuclear matter, quark matter, and the mixed phases from the same NJL Lagrangian, which has been used to model free and in-medium hadrons as well as nuclear systems. A focus here is to determine if the same coupling constants in the scalar diquark and vector meson channels, which give a good description of nucleon structure and nuclear matter, can also be used for the color superconducting high density quark matter phase. We find that this is possible for the scalar diquark (pairing) interaction, but the vector meson interaction has to be reduced so that superconducting quark matter becomes the stable phase at high densities. We compare our equation of state with recent phenomenological parametrizations based on generic stability conditions for neutron stars. We find that the maximum mass of a neutron star, with a color superconducting quark matter core, exceeds $2.01 \pm 0.04\,M_\odot$ which is the value of the recently observed massive neutron star PSR J0348+0432. The mass-radius relation is also consistent with gravitational wave observations (GW170817).

Figures

c in Ref. Alford *et al.* ([2013](#bib.bib56)), i.e., the pure neutron star and hybrid star configurations form a connected sequence. (As already mentioned, for case A this is not immediately evident on the scale of Fig. [IV](#S4)a.) A connected sequence of this type has also been obtained in the calculations of Ref. Abgaryan *et al.* ([2018](#bib.bib42)), see their Fig. 7 for the case of $\Delta_{P}>5\%$.

c in Ref. Alford *et al.* ([2013](#bib.bib56)), i.e., the pure neutron star and hybrid star configurations form a connected sequence. (As already mentioned, for case A this is not immediately evident on the scale of Fig. [IV](#S4)a.) A connected sequence of this type has also been obtained in the calculations of Ref. Abgaryan *et al.* ([2018](#bib.bib42)), see their Fig. 7 for the case of $\Delta_{P}>5\%$.

Tables

Table 1. Table 1: Values for the model parameters which are determined in the vacuum, single hadron, and nuclear matter sectors. The regularization parameters and the current quark mass are given in units of GeV, and the coupling constants in units of GeV -2 .

${\Lambda }_{\mathrm{IR}}$	${\Lambda }_{\mathrm{UV}}$	$G_{\pi }$	$G_S$	$G_A$	$G_{\omega }$	$G_{\rho }$	$m$
0.240	0.645	19.04	11.24	4.20	6.03	14.17	0.016

Table 2. Table 2: Transition densities ρ tr subscript 𝜌 tr \rho_{\mathrm{tr}} and maximum star masses M star max superscript subscript 𝑀 star max M_{\rm star}^{\rm max} for the cases A, B, and C, in comparison to the case R of Fig.  1 .

Case	$c_s$	$c_v$	${\rho }_{\mathrm{tr}}\left[{\mathrm{fm}}^{-3}\right]$	$M_{\mathrm{star}}^{\max }\left[M_{\odot }\right]$
$A$	0.80	0.50	0.643	2.078
$B$	0.90	0.60	0.584	2.055
$C$	0.98	0.68	0.496	2.071
$R$	0.50	0.90	0.379	1.918

Equations

$$\mathcal{L}$$

$$\displaystyle-G_{\omega}\big{(}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\gamma^{\mu}\psi\big{)}^{2}-G_{\rho}\big{(}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\gamma^{\mu}\vec{\tau}\psi\big{)}^{2}$$

$${\mu }_p$$

$${\mu }_u$$

$$V^{(\mathrm{NM})}(M,{\omega }_0,{\rho }_0)$$

$$M$$

$${\omega }_0$$

$$V_{\mathrm{vac}}=12i\int \frac{{\mathrm{d}}^{4}k}{(2\pi {)}^4}\ln \frac{k^2-M^2}{k^2-M_0^2}+\frac{(M-m{)}^2}{4G_{\pi }}-\frac{(M_0-m{)}^2}{4G_{\pi }}.$$

$$V_N=-2\sum _{\alpha =p,n}\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}\left({\mu }_{\alpha }^{\ast }-E_N(k)\right)\Theta \left({\mu }_{\alpha }^{\ast }-E_N(k)\right),$$

$${\mu }_{\alpha }^{\ast }={\mu }_{\alpha }-3{\omega }_0\mp {\rho }_0(\alpha =p,n),$$

$$V^{(\mathrm{QM})}(M,\Delta ,{\omega }_0,{\rho }_0)=V_{\mathrm{vac}}+V_Q+V_{\Delta }-\frac{{\omega }_0^2}{4G_{\omega }}-\frac{{\rho }_0^2}{4G_{\rho }}-\frac{{\mu }_e^4}{12{\pi }^2},$$

$$\Delta =-2G_S⟨{\psi }^{T}C{\gamma }_5{\tau }_2{\lambda }_2\psi ⟩.$$

$$V_{\mathrm{Q}}=-6\sum _{\alpha =u,d}\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}\left({\mu }_{\alpha }^{\ast }-E(k)\right)\Theta \left({\mu }_{\alpha }^{\ast }-E(k)\right).$$

$${\mu }_{\alpha }^{\ast }={\mu }_{\alpha }-{\omega }_0\mp {\rho }_0(\alpha =u,d).$$

$$V_{\Delta }$$

$$+\ln \frac{k_0^2-{\left({ϵ}_{\alpha }-{\mu }_I^{\ast }\right)}^2}{k_0^2-{\left(E_{\alpha }-{\mu }_I^{\ast }\right)}^2}]+\frac{{\Delta }^2}{4G_S},$$

$${\mu }_B^{\ast }$$

$$P$$

$$\mathcal{E}$$

$$x_1({\mu }_B,{\mu }_I)$$

Full text

I INTRODUCTION

The study of strongly interacting cold matter at high baryon densities — neutron star matter — has been a very important and active area of research for several decades Itoh (1970); Baym et al. (1971); Glendenning (1997); Heiselberg and Hjorth-Jensen (2000); Dean and Hjorth-Jensen (2003); Baym et al. (2018). This field has attracted increased attention recently however, because of the observation of massive neutron stars exceeding two solar masses Demorest et al. (2010); Antoniadis et al. (2013) and gravitational wave measurements of a binary neutron star merger event Abbott et al. (2018, 2017); Evans et al. (2017). Among the many theoretical tools used to study dense hadronic matter are nonrelativistic potential models Gandolfi et al. (2012, 2014), effective field theories Hebeler et al. (2011); Kaiser (2012), and relativistic theories based on a mean field description of nucleons interacting via meson exchange Muller and Serot (1995); Chin and Walecka (1974). The role of hyperons in dense hadronic matter has also been of interest Kaplan and Nelson (1986); Dapo et al. (2008). On the other hand, studies based on effective quark theories of QCD — often using the framework provided by the Nambu–Jona-Lasinio (NJL) model Nambu and Jona-Lasinio (1961a, b); Vogl and Weise (1991); Hatsuda and Kunihiro (1994) — have focused on the existence of a possible quark matter phase at high densities Buballa (2005); Shovkovy (2005); Alford et al. (2008); Fukushima and Hatsuda (2011) and the role of color superconductivity Bailin and Love (1984); Alford et al. (1998); Rapp et al. (1998). In these approaches the transition between the hadronic and quark matter phases has been described by using either the Maxwell or Gibbs constructions Glendenning (1992); Wu and Shen (2017), or by interpolation methods based on hadron-quark continuity Schäfer and Wilczek (1999); Macher and Schaffner-Bielich (2005); Masuda et al. (2013). Finite size effects in the mixed phase caused by the surface tension of various geometrical shapes, Coulomb interactions, and charge screening have also been studied Heiselberg et al. (1993); Endo et al. (2006); Lugones et al. (2013).

The aim of this paper is to present results for the equation of state of high density matter and the properties of neutron stars, using effective quark degrees of freedom for the description of the hadronic as well as the quark phase. Most of the calculations done so far on hybrid star matter — hadronic matter which converts into quark matter at high baryon densities — combined a relativistic hadronic mean field theory for nuclear matter with some version of the NJL model for quark matter, see for example Refs. Baym et al. (2018); Contrera et al. (2014); Benic et al. (2015); Li et al. (2018). More phenomenological approaches used a parametrization of the nuclear and quark matter phases in terms of polytropes Abgaryan et al. (2018) or relativistic density functionals Kaltenborn et al. (2017). In Ref. Hempel et al. (2013) an interesting attempt was made to describe a mixture of elementary hadrons and quarks by using the same Lagrangian, the mixture being determined by a scalar field which increases the hadron masses and decreases the quark masses as the baryon density increases. An application of this approach to neutron stars can be found in Ref. Dexheimer et al. (2018).

The purpose of our present work differs from the above mentioned approaches in the following respect: We wish to investigate whether a single effective quark theory of QCD, which can describe the quark structure of free hadrons Cloët et al. (2005a, 2008, 2014); Hutauruk et al. (2016); Ninomiya et al. (2017) as well the role of quarks in bound nucleons and nuclear systems Cloët et al. (2005b, 2006, 2016), can also produce reasonable scenarios for the hadron-quark phase transition and the properties of neutron stars. Our work is based on the NJL model with the proper-time regularization scheme Hellstern et al. (1997), which incorporates important aspects of confinement in hadronic systems Bentz and Thomas (2001). The crucial point which leads to saturation of the nuclear matter binding energy in the mean field approximation is the scalar polarizability of the in-medium nucleons Bentz and Thomas (2001), i.e., the non-linear dependence of the nucleon mass on the constituent quark mass, which arises naturally in the Faddeev approach based on the quark-quark correlations in the scalar and axial-vector diquark channels Cloët et al. (2014). A key question that we will address in this work is whether the same strength of the scalar diquark interaction — which is required to reproduce the nucleon mass and other quark-quark correlation effects in baryons — can also be used as the pairing interaction in color superconducting quark matter.

In the course of our investigation, we will formulate a few conditions for a reasonable scenario of the hadron-quark phase transition and the resulting properties of neutron stars, and we will show how these conditions can be satisfied in our model. An important question which we wish to address is whether the resulting picture is consistent with the generic stability conditions for a hybrid star, formulated in Ref. Alford et al. (2013) and generalized in Ref. Ranea-Sandoval et al. (2016). Closely related to this is the role of the repulsive effects in quark matter which arise from the interaction in the vector meson channels. The importance of the isoscalar-vector repulsion in producing a sufficiently stiff quark matter equation of state has been pointed out in several recent papers Baym et al. (2018); Whittenbury et al. (2016); Hell et al. (2013), and here we will also include the isovector-vector repulsion which is very important in nuclei Cloët et al. (2009, 2012). We will compare the required strengths to those adjusted to the saturation density and symmetry energy of nuclear matter.

In order to gain insight into these and other related questions, we will confine ourselves to two light quark flavors in both the nuclear matter and the quark matter phases, and neglect finite size effects in the hadron-quark mixed phase. It is believed on rather general grounds that at very high baryon densities a three-flavor color superconducting phase of quark matter is realized Alford et al. (2008), but whether nuclear matter directly jumps to this phase or to an intermediate two-flavor color superconducting quark phase depends on model details like the in-medium strange quark mass and the strength of the interaction in the relevant diquark channels Buballa (2005); Ruester et al. (2005); Abuki and Kunihiro (2006); Steiner et al. (2002). In any case, the role of strangeness should be investigated consistently both in the hadronic phase, based for example on the quark-diquark description of the baryon octet as given in Ref. Carrillo-Serrano et al. (2016), as well as in the quark phase, and we wish to leave such extensions to a future work. Concerning finite size effects, it has been argued in several papers Endo (2011); Neumann et al. (2003); Voskresensky et al. (2003); Pinto et al. (2012); Fraga et al. (2019) that surface tension and charge screening tend to work against spatially extended mixed phases, regaining qualitatively the simple picture of a Maxwellian first-order phase transition. We will come back to this point in later sections.

The outline of the paper is as follows: Sec. II discusses the NJL model and the parameters that enter the calculations; Sec. III introduces the equations of state for nuclear matter and quark matter; Sec. IV presents a discussion of our results and a comparison with data and other related work; and a summary is given in Sec. V.

II LAGRANGIAN AND MODEL PARAMETERS

The two-flavor NJL Lagrangian relevant for this study reads Vogl and Weise (1991); Hatsuda and Kunihiro (1994):

[TABLE]

where $\psi$ is the quark field, $m$ is the current quark mass, $C$ is the charge conjugation matrix, $\lambda^{A}\,\,\,(A=2,5,7)$ are the antisymmetric color Gell-Mann matrices, and $\bm{\tau}$ are the Pauli isospin matrices. The 4-fermion coupling constants in the scalar $\bar{q}q$ channel, the isoscalar and isovector vector $\bar{q}q$ channels, and the scalar $qq$ channel are denoted by $G_{\pi}$, $G_{\omega}$, $G_{\rho}$ and $G_{S}$, respectively.111The coupling constant $G_{S}$ of Eq. (1) is the same as $G_{D}$ used in most works on the NJL model for high density quark matter. It is related to $G_{s}$ used in our previous works on the NJL model (see for example Eq. (7) of Ref. Cloët et al. (2009)) by $G_{S}=\frac{3}{2}G_{s}$. The other model parameters are the 4-fermion coupling constant in the axial-vector $qq$ channel Cloët et al. (2014), and the infrared (IR) and ultraviolet (UV) regularization parameters, which are used with the proper-time regularization scheme Hellstern et al. (1997); Bentz and Thomas (2001).222The interaction Lagrangian in the axial-vector diquark channel has the form $G_{A}\left(\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\gamma_{\mu}C\tau_{i}\tau_{2}\lambda^{A}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}^{T}\right)\left(\psi^{T}C^{-1}\gamma^{\mu}\tau_{2}\tau_{i}\lambda^{A}\psi\right)$, where $G_{A}$ is related to $G_{a}$ of our previous works (see for example Eq. (7) of Ref. Cloët et al. (2009)) by $G_{A}=\frac{3}{2}G_{a}$. Because we do not consider spin triplet pairing here, this interaction is used only in the quark-diquark bound state equation. We also mention that chiral symmetry of the interaction Lagrangian requires additional terms in the vector $qq$ and the axial-vector $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{q}q$ channels, which however are not directly related to our calculation.

We stress that in this work we use the same model parameters as in several previous calculations which focused on the structure of nuclear matter systems Cloët et al. (2006, 2009, 2012). That is, the parameters of the model are determined in the vacuum, the single hadron sector, and nuclear matter sector as follows: We fix the IR cut-off $\Lambda_{\mathrm{IR}}=0.24\,$GeV, and choose the UV cut-off $\Lambda_{\mathrm{UV}}$, $m$, and $G_{\pi}$ so as to give a constituent quark mass in vacuum of $M_{0}=0.4\,$GeV, the pion decay constant $f_{\pi}=0.93\,$GeV, and the pion mass $m_{\pi}=0.14\,$GeV using the standard methods based on the gap and Bethe-Salpeter equations Cloët et al. (2014). The scalar diquark coupling $G_{S}$ and its axial-vector counterpart $G_{A}$ are then determined in the Faddeev equation approach to reproduce the vacuum values of the nucleon mass ($M_{N0}=0.94\,$GeV) and the nucleon axial coupling ($g_{A0}=1.267$) Cloët et al. (2005b). Finally, by using the model description for nuclear matter explained in Sec. III, the vector couplings $G_{\omega}$ and $G_{\rho}$ are determined from the binding energy per-nucleon in symmetric nuclear matter ($E_{B}=-15.7\,$MeV) and the symmetry energy ($a_{4}=32.0\,$MeV) at the saturation density Cloët et al. (2009). We note that in this framework the saturation density of symmetric nuclear matter becomes $\rho_{0}=0.16\,$fm*-3*. The resulting values for the model parameters are given in Tab. 1.333Note that $G_{\pi}$ obtained in the present proper-time regularization scheme is almost three times larger than the value obtained in the 3-momentum cut-off scheme (see for example Ref. Buballa (2005)). As a consequence, for the same vacuum value $M_{0}=0.4$ GeV, our vacuum chiral condensate is roughly three times smaller in magnitude than in the 3-momentum cut-off scheme. Nevertheless, our value $\langle\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\psi\rangle_{0}^{1/3}=-0.216$ GeV is still close to the upper limit of the range $[-0.33,-0.24]$ derived from QCD sum rules at a renormalization scale of 1 GeV  Dosch and Narison (1998). By choosing smaller values for the input $M_{0}$ the magnitude of the chiral condensate increases Ninomiya et al. (2015). Using the nuclear matter equation of state presented in Sec. III, and the parameters of Tab. 1, gives an effective nucleon mass of $M_{N}=0.744\,$GeV at nuclear matter saturation density, and an incompressibility of $K=0.370\,$GeV.

In the quark matter phase there should be no effects from color confinement so we set $\Lambda_{\rm IR}=0$ in this phase. For the parameters $\Lambda_{\rm UV}$, $m$ and $G_{\pi}$, which follow from the properties of the vacuum and the pion, we use the same values as in nuclear matter (see Tab. 1). As discussed in Sec. I, we wish to address the question whether the same coupling constants $G_{S}$, $G_{\omega}$ and $G_{\rho}$ as determined by the free nucleon mass and the properties of nuclear matter can also be used to describe the phase transition to color superconducting quark matter and neutron stars. We therefore introduce two scaling parameters $c_{s}$ and $c_{v}$ such that in the quark matter phase $G_{S}\to c_{s}\,G_{S}$, $G_{\omega}\rightarrow c_{v}\,G_{\omega}$, and $G_{\rho}\rightarrow c_{v}\,G_{\rho}$. (We use a common scaling factor for the vector-isoscalar and vector-isovector interaction in order to reduce the number of parameters.) The dependence of our results on the scaling parameters $(c_{s},c_{v})$ will be investigated in detail in Sec. IV.

III NUCLEAR MATTER AND QUARK MATTER

In this section we present expressions for the effective potential ($V$) of nuclear matter (NM) and quark matter (QM) in the mean field approximation of the two-flavor NJL model for given values of the two independent chemical potentials $\mu_{B}$ and $\mu_{I}$ for baryon number and isospin. The corresponding chemical potentials for nucleons and quarks are444In principle a further chemical potential ($\mu_{8}$) is needed in QM to guarantee that the mean value of the color is zero, but it turns out to be very small for the two-flavor case Buballa (2005); Aguilera et al. (2005) and we neglect it here for simplicity. We also note that in many works the chemical potentials of $u,d$ quarks are expressed as $\mu_{q}=\mu+Q_{q}\,\mu_{Q}$ ($q=u,d$), where $Q_{q}$ is the electric charge. Physical quantities as functions of baryon and charge density of course do not change with different definitions of chemical potentials for baryon number and charge.

[TABLE]

The electron Fermi gas terms are also included, with the chemical potential $\mu_{e}=-2\mu_{I}$ determined by $\beta$ equilibrium. Muon contributions will not be included for simplicity. Below we summarize the unregularized expressions, and refer to Ref. Bentz and Thomas (2001) for a detailed discussion on the proper-time regularization scheme.

The effective potential of NM for fixed chemical potentials $\mu_{B},\mu_{I}$ is given by Bentz et al. (2003); Lawley et al. (2006)

[TABLE]

where $M$ is the (in-medium) constituent quark mass and $\omega_{0}$ and $\rho_{0}$ are the isoscalar- and isovector-vector mean fields, which must be determined by minimization of the effective potential. Those minimization conditions are equivalent to the following definitions in terms of the in-medium quark condensates:

[TABLE]

The vacuum (Mexican hat shaped) contribution is

[TABLE]

The Fermi motion of nucleons in the scalar and vector mean fields is described by the term

[TABLE]

where the effective chemical potential for protons and neutrons is given by

[TABLE]

and $E_{N}(k)=\sqrt{\bm{k}^{2}+M_{N}^{2}}$. The nucleon mass $M_{N}(M)$ is determined as a function of the constituent quark mass $M$ from the relativistic Faddeev equation for the nucleon, which is approximated as a quark-diquark bound state Cloët et al. (2014). As mentioned in Sec. I, the function $M_{N}(M)$ develops a (positive) curvature with decreasing $M$, which is crucial for the saturation of the NM binding energy Bentz and Thomas (2001).

The effective potential for QM for fixed chemical potentials $\mu_{B}$ and $\mu_{I}$ is given by Bentz et al. (2003); Lawley et al. (2006)

[TABLE]

where the constituent quark mass $M$ and the mean vector fields $\omega_{0}$ and $\rho_{0}$ in QM are given in terms of quark condensates as in Eqs. (5) and (6), while the energy gap $\Delta$ arising from the pairing in the scalar diquark channel is related to the BCS-type quark condensate (order parameter of the broken color symmetry) by

[TABLE]

It is well known Shovkovy (2005); Alford et al. (2008) that this type of 2-flavor pairing leaves chiral symmetry intact, and in the limit $m=0$ color superconducting quark matter may therefore exist in the chiral symmetric phase ($M=0$). All four quantities $M,\Delta,\omega_{0},\rho_{0}$ are determined by minimization of the effective potential for fixed chemical potentials.

The vacuum part $V_{\mathrm{vac}}$ in Eq. (12) is the same as given in Eq. (7). The difference between its value at $M$ calculated in QM and its value at $M$ calculated in NM corresponds to the bag constant in the NJL model. The term $V_{Q}$ describes the Fermi motion of quarks moving in the scalar and vector mean fields:

[TABLE]

Here $E(k)=\sqrt{\bm{k}^{2}+M^{2}}$ and the effective up and down quark chemical potentials are defined by

[TABLE]

The term $V_{\Delta}$ in Eq. (12) arises from the pairing in the scalar diquark channel, and is given by (see the papers on color superconductivity cited in Sec. I and also Refs. Kiriyama et al. (2001); Huang and Shovkovy (2003))

[TABLE]

where $\epsilon_{\pm}=\sqrt{\left[E(k)\pm\mu_{B}^{*}/3\right]^{2}+\Delta^{2}}$ is the quark dispersion relation555As can be seen from our results in Sec. IV, we always have $\Delta>|\mu_{I}^{*}|$, i.e., the “gapless modes” Shovkovy (2005) (poles at zero frequency for finite $\Delta$) are not realized in our calculation, which is consistent with the findings of other works Ruester et al. (2005) for intermediate or strong pairing strength. and $E_{\pm}=\left|E(k)\pm\mu_{B}^{*}/3\right|$. The effective chemical potentials for baryon number and isospin are

[TABLE]

Using the above forms for the effective potentials, the pressure ($P$), baryon and isospin densities ($\rho_{B}$, $\rho_{I}$), and the energy density (${\cal E}$) are obtained by

[TABLE]

The charge neutrality condition $\rho_{C}=\left(\rho_{B}+\rho_{I}\right)/2=0$ then implies a relation between the two chemical potentials $\mu_{B}$ and $\mu_{I}$. The charge neutral equation of state is then a function of only one variable, which we take to be the baryon density $\rho_{B}$.

The equation of state for the globally charge neutral mixed phase is then calculated by using the Gibbs construction Glendenning (1992) as follows: If there is a line in the $(\mu_{B},\mu_{I})$ plane between two points (called X and Y) along which the NM and QM phases have equal effective potentials (denoted as $V^{\rm(mixed)}$) and opposite charges, then along this line we have trivially $V^{\rm(mixed)}=x_{1}\,V^{\rm(NM)}+(1-x_{1})\,V^{\rm(QM)}$ for any number $x_{1}$. We can calculate the densities and the energy density along this line by differentiation of $V^{\rm(mixed)}$ according to Eq. (18), noting that $\partial V^{\rm(mixed)}/\partial x_{1}=0$. The requirement that the charge density along this line vanishes, i.e., $\rho_{C}=x_{1}\rho_{C}^{\rm(NM)}+(1-x_{1})\rho_{C}^{\rm(QM)}=0$, determines $x_{1}(\mu_{B},\mu_{I})$ as the volume fraction of NM in the mixed phase as

[TABLE]

If we approach the point X along a charge neutral line of the NM phase, then $x_{1}=1$ at point X, and if we leave point Y towards a charge neutral line of the QM phase, then $x_{1}=0$ at point Y, i.e., $x_{1}$ decreases from $1$ to [math] along the line X $\rightarrow$ Y. The baryon density in the mixed phase is then given by $\rho_{B}^{\rm(mixed)}=x_{1}\rho_{B}^{\rm(NM)}+(1-x_{1})\rho_{B}^{\rm(QM)}$, and a similar expression holds also for the energy density.

IV RESULTS AND DISCUSSION

Before presenting our results, we would like to explain how they qualitatively depend on the choice of parameters. The equation of state in the QM phase is largely controlled by three parameters in our model: the strength of the attractive isoscalar-scalar pairing interaction ($G_{S}$) and the strengths of the repulsive isoscalar- and isovector-vector interactions $G_{\omega}$, $G_{\rho}$. If we increase the pairing attraction, the baryon density $\rho_{\rm tr}$ where the transition NM $\rightarrow$ QM occurs decreases, resulting in an overall softening of the equation of state and lower masses of neutron stars. If the pairing interaction is increased beyond a certain limit, QM becomes the stable phase of the system also at low densities, which we regard as unphysical. On the other hand, if we increase the vector repulsion the transition density $\rho_{\rm tr}$ increases, then the equation of state becomes stiffer and the masses of the neutron stars increase. If the vector repulsion is increased beyond a certain limit, NM becomes the stable phase also for high densities, which we again regard as unphysical.

As explained at the end of Sec. II, in order to show these dependencies quantitatively we scale the value of $G_{S}$ given in Tab. 1 by a factor $c_{s}$, and the values of $G_{\omega}$ and $G_{\rho}$ by a common factor $c_{v}$, in the calculation of the QM sector. To characterize the results obtained for different values of $c_{s}$ and $c_{v}$ we define a “physically reasonable” scenario by the following three conditions: (1) The phase transition NM $\rightarrow$ QM occurs in the range $2\,\rho_{0}\leq\rho_{\rm tr}\leq 4\,\rho_{0}$;666There is no fundamental reason for this choice of limits and our qualitative results do not change much if this condition is relaxed. (2) the maximum mass of the star satisfies $M_{\rm star}^{\rm(max)}\geq 1.97\,M_{\odot}$ to be compatible with recent observations Demorest et al. (2010); Antoniadis et al. (2013); and (3) the hybrid star with a superconducting QM core is stable against density fluctuations, i.e., ${\rm d}M_{\rm star}/{\rm d}\rho_{B}(r=0)>0$ where $\rho_{B}(r=0)$ is the central baryon density, in a region of densities above the transition density.

We find that all three conditions can be satisfied in the rather narrow parameter region marked by yellow in Fig. 1. In the blue region of Fig. 1, which continues towards smaller values of $(c_{s},c_{v})$ not included in the figure, at least one of the conditions (2) or (3) is not satisfied. Because the yellow region extends up to $c_{s}=0.98$, we find that in practice it is possible to use the same value of $G_{S}$ for the scalar pairing strength in QM as for the scalar diquark interaction in the single nucleon sector. On the other hand, if we would use the same value of the vector couplings $G_{\omega}$, $G_{\rho}$ as determined in the nuclear matter sector (see Tab. 1), there would be no phase transition to QM. In order to satisfy the three conditions explained above, the vector coupling in QM must be smaller than in NM by a factor of $0.45\leq c_{v}\leq 0.68$.777Interestingly, using the same NJL framework to explore possible explanations of the EMC effect Aubert et al. (1983); Geesaman et al. (1995); Malace et al. (2014); Cloët et al. (2019) in nuclear structure functions, we found that the coupling of the vector mean field to the struck quark must also be substantially reduced Cloët et al. (2006). In this case the reduction was associated with asymptotic freedom in QCD Detmold et al. (2006). In the white region marked as “NM only”, either the transition density is too high or only the NM phase is realized, while in the white region marked as “QM only”, either the transition density is too low or only the QM phase is realized.

Figure 2: (Color online) Phase diagrams for the point B of Fig. 1 ($c_{s}=0.9,c_{v}=0.6$) in: (a) the plane of chemical potentials, and (b) in the plane of densities. In the top panel (a) the “$+$” and “$-$” signs refer to the electric charge state of nuclear matter (NM) or quark matter (QM). The black dashed line markes electrically neutral matter, with the arrows indicating the increase of baryon density. The segment X $\rightarrow$ Y markes the mixed phase.

In order to discuss our results for the phase structure, we select point B of Fig. 1 as a representative example where all three conditions are satisfied. Figs. IV shows the phase diagrams in: (a) the plane of chemical potentials, and (b) the plane of densities, focusing on the region of the phase transition for electrically neutral matter. In Fig. IVa the stable phase (phase with the larger pressure) is marked as NM or QM for each point in the $(\mu_{B},\mu_{I})$ plane, together with the sign of the electric charge density. The dashed line marks electrically neutral matter. The baryon density increases as we follow this line from the left end in the NM phase to the right end in the QM phase. We see that in the section of the NM-QM mixed phase — the line X $\rightarrow$ Y in Fig. IVa — $\mu_{B}$ stays almost constant (within our numerical accuracy of $\pm 1$ MeV), while $|\mu_{I}\!|$ increases by about 10 MeV.888Note, however, that the constancy (almost zero width) of $\mu_{B}$ during the phase transition depends somewhat on the definitions of the chemical potentials, see Footnote 4. With our present definitions, we find very small widths ($<3$ MeV) for most of the parameters $(c_{s},c_{v})$ in the colored regions of Fig. 1. Large widths occur only for weak pairing strength and weak vector interaction ($c_{s}\lesssim 0.2$, $c_{v}\simeq 0$), which is consistent with the results of Refs. Bentz et al. (2003); Lawley et al. (2006), but in this region the conditions (2) and (3) cannot be satisfied at the same time.

Fig. IVb shows the phase diagram in the plane of the densities $(\rho_{B},\rho_{C})$ in the vicinity of $\rho_{C}=0$. The dashed line, which marks electrically neutral matter, passes through the region of the mixed phase (shown in white). Figs. IV illustrates that along the phase transition line X $\rightarrow$ Y there is only small change of the chemical potentials, while the change in baryon density (about $0.1$ fm*-3*) is appreciable but not too large, so as to keep our resulting hybrid star gravitationally stable as will be shown below. We found qualitatively very similar results for the phase diagrams for all parameter sets $(c_{s},c_{v})$ in the colored regions of Fig. 1.

Figure 3: (Color online) The dependence of physical quantities on the baryon chemical potential in electrically neutral matter for the case of point B in Fig. 1 ($c_{s}=0.9,c_{v}=0.6$). Figure (a) shows the baryon density and (b) shows the constituent quark mass $M$ and the energy gap $\Delta$.

In Fig. IVa we show the baryon density of electrically neutral matter as function of the baryon chemical potential for the parameter set B. As can be anticipated from Figs. IV, the change of the baryon density during the phase transition appears as a jump, i.e., within our numerical accuracy (of $\pm 1$ MeV) it occurs at constant $\mu_{B}$. In Fig. IVb we show the constituent quark mass $M$ and the energy gap $\Delta$ as functions of $\mu_{B}$ for the same parameter set B. The values of $M$ in the QM phase are small compared to those in the vacuum or the NM phase, which indicates that QM in the range of $\mu_{B}$ above the phase transition is already reasonably close to a phase where chiral symmetry is restored. Our values of $\Delta$ are larger than in most of the previous works done in the 3-momentum cut-off scheme (see for example Refs. Buballa (2005); Kiriyama et al. (2001)), which is mainly because of the larger coupling constant $G_{S}$ in Eq. (13) for the proper-time regularization scheme, see Tab. 1 for the reference value. (Our values for the BCS-type condensate in Eq. (13) are very similar to the values obtained with the 3-momentum cut-off scheme.) Nevertheless, in the region of $\mu_{B}$ just above the phase transition, we observe a qualitative agreement of our results shown in Fig. IVb with the results for the two-flavor quark phase shown in Fig. 33 of Ref. Baym et al. (2018), or Fig. 6.8 of Ref. Buballa (2005). We found that all parameter sets $(c_{s},c_{v})$ shown by the yellow region in Fig. 1 give results which are qualitatively very similar to those shown in Figs. IV.

Figure 4: (Color online) Pressure as function of the baryon density for electrically neutral matter. Figure (a) shows the results obtained for the points A and B in comparison to point R of Fig. 1, and (b) shows the result obtained for the point C. The black solid line is the result of the pure NM case.

In order to separately show the counteracting effects of pairing and vector repulsion in the QM phase, we show in Fig. IVa the pressure of electrically neutral matter as a function of the baryon density for the points A, B and R in Fig. 1, as well as the result for the pure NM case. Here the point R, which satisfies the above conditions (1) and (3) but not (2), is used as a reference point. Starting with case A, we see that by increasing the pairing strength (${\rm A}\rightarrow{\rm R}$) the transition density decreases substantially, resulting in an overall softer equation of state. If the vector coupling is then increased (${\rm R}\rightarrow{\rm B}$) the transition density increases again, but stays below case A. As a result, there is a net decrease of the transition density and a resulting softening of the equation of state as we go along the solid line of Fig. 1 in the direction of increasing coupling constants (A $\rightarrow$ B). This trend continues as we extend the line $c_{s}=c_{v}+0.3$ up to the point C, for which the results are shown in Fig. IVb. In this case $c_{s}=0.98$, i.e., the coupling constant $G_{S}$ is practically the same as determined from the free nucleon mass and the vector couplings $G_{\omega}$, $G_{\rho}$ are smaller than the nuclear matter values by a factor of $c_{v}=0.68$. The transition density is about $3\,\rho_{0}$ for this case.

The results shown in Figs. IV indicate that the phase transition, which we described here by the Gibbs construction, is in fact very similar to the usual Maxwellian first order phase transition. This result, which for the case B can be anticipated from Figs. IV and IV, holds for the whole yellow region of parameters shown in Fig. 1. As a consequence, the spatial extension of the mixed phase in a hybrid star will be small compared to the overall size, i.e., if the central density is sufficiently high the QM phase will begin to form at the center of the star with an almost sharp boundary to the surrounding NM phase, which is the scenario depicted schematically in Fig. 25 of Ref. Heiselberg and Hjorth-Jensen (2000). Our results for the phase transition also suggest that the inclusion of finite size effects (surface tension of various geometrical shapes, Coulomb interactions and charge screening) in the mixed phase will not lead to qualitative changes of the overall physical picture: As mentioned in Sec. I, these effects work against spatially extended mixed phases, regaining qualitatively the simple picture of a Maxwellian first-order phase transition  Endo (2011); Neumann et al. (2003); Voskresensky et al. (2003).

In order to see whether our results are consistent with phenomenological parametrizations based on generic stability criteria for hybrid stars Alford et al. (2013); Ranea-Sandoval et al. (2016), we show our equation of state (energy density vs. pressure) for case B in Fig. 5. We immediately see that the physical picture is consistent with those parametrizations: The increase of the energy density during the phase transition ($\Delta{\cal E}=136.7$ MeV/fm3) is smaller than the critical value determined from Eq. (2) of Ref. Alford et al. (2013) ($\Delta{\cal E}_{\rm crit}=545.9$ MeV/fm3 in our case), which means that the pressure of QM can counteract the additional downward gravitational pull — exerted by the additional energy in the core — on the rest of the star. We can see also consistency with the more general parametrizations based on a constant speed of sound in QM Ranea-Sandoval et al. (2016), and our value $c_{\rm QM}^{2}=0.56$ is sufficiently high to get a heavy hybrid star Alford et al. (2013).

Figure 6: (Color online) Neutron star masses as functions of the central baryon density $\rho_{B}(r=0)$. Figure (a) shows the results obtained for the points A and B in comparison to point R of Fig. 1, and the (b) shows the result obtained for the point C. The black solid line shows the pure NM case.

By using our equations of state as an input to solve the Tolman-Oppenheimer-Volkoff (TOV) equations Tolman (1939); Oppenheimer and Volkoff (1939), we can calculate the properties of neutron stars. The results for the neutron star mass as a function of the central baryon density are shown in Fig. IVa for the cases of points A, B, and R of Fig. 1, and in Fig. IVb for case C. For case A (high transition density) one can obtain the largest star masses, but the range of stability against density fluctuations (${\rm d}M_{\rm star}/{\rm d}\rho_{B}(r=0)>0$) is narrow, i.e., the hybrid star tends to be unstable. If we increase the pairing strength (A $\rightarrow$ R), the overall equation of state becomes softer and the neutron star masses become smaller, however the stability of the hybrid star against density fluctuations is substantially improved. Increasing the vector coupling (R $\rightarrow$ B), we can obtain stars which satisfy $M_{\rm star}^{\rm(max)}\geq 1.97\,M_{\odot}$ and show a reasonable range of stability. Fig. IVb shows the results obtained for the point C in Fig. 1. Here the star mass near the onset of the phase transition is smaller than for the cases A and B, which expresses the softening of the equation of state near the phase transition density along the line A $\rightarrow$ B $\rightarrow$ C (see Fig. IV), but for very high central densities the QM equation of state is stiff enough to support a heavy hybrid star within a reasonably broad range of stability.

In Figs. IV we also indicate the recently observed values of neutron star masses. GW170817 denotes the neutron star coalescence event observed by the gravitational wave measurements Abbott et al. (2018, 2017); Evans et al. (2017), where neutron stars with masses in the range $M_{\rm star}=1.17\sim 1.60\,M_{\odot}$ were observed and PSR J0348+0432 refers to the observation of massive neutron stars (pulsars) of mass $M_{\rm star}=2.01\pm 0.04\,M_{\odot}$ Antoniadis et al. (2013). The result of PSR can be considered as a lower limit for calculations of the maximum mass of neutron stars. We see that our parameter sets A, B, and C satisfy this constraint, but the maximum star mass for case R is too small. In Tab. 2 we list the transition densities and maximum neutron star masses for the cases A, B, and C, in comparison to the reference case R.

Figure 7: (Color online) Mass–radius relation for neutron stars. Figure (a) shows the results obtained for the points A and B in comparison to point R of Fig. 1, and (b) shows the result obtained for the point C. The black solid line shows the pure NM case.

In Fig. IVa we show the relation between the neutron star masses and radii for the cases A, B, and R of Fig. 1, and Fig. IVb shows the results for the case C. The low density part of the NM curve (lower part of the solid line in Figs. IV) shows the characteristics of a case where the pressure drops to zero (or nearly to zero) at a finite value of the baryon density, as the density decreases,999If the pressure drops to zero at $\rho_{B}=0$, the lower part of the NM curves in Figs. IV would turn to the right side, instead of turning to the left. which is indicative of a bound state in the absence of gravity. Fig. IVa clearly shows that QM in case A (highest value of $\rho_{\rm tr}$) tends to be unstable. (The narrow region of stability shown in Fig. IVa for case A becomes invisible in Fig. IVa.) Similar to the previous figures, the effects of the softening of the QM equation of state caused by increasing the pairing strength (A $\rightarrow$ R) and the improvement of the stability caused by increasing the vector couplings (R $\rightarrow$ B) are clearly seen in Fig. IVa. The behavior of the curves with increasing vector couplings (R $\rightarrow$ B) are qualitatively similar to the results shown in Refs. Contrera et al. (2014); Wu and Shen (2017), although in those works only small vector couplings could be used because the counteracting effects of pairing were not taken into account. The overall softening of the QM equation of state and the resulting improvement of the stability along the line A $\rightarrow$ B of Fig. 1 continues further as we go on to the case C. In this case, the star mass at the onset of the phase transition is lower than in cases A and B, but the QM equation of state is stiff enough to support heavy hybrid stars with radii of about 11 km. Our results for all cases shown in Figs. IV follow qualitatively the “connected” topology of Fig. 2c in Ref. Alford et al. (2013), i.e., the pure neutron star and hybrid star configurations form a connected sequence. (As already mentioned, for case A this is not immediately evident on the scale of Fig. IVa.) A connected sequence of this type has also been obtained in the calculations of Ref. Abgaryan et al. (2018), see their Fig. 7 for the case of $\Delta_{P}>5\%$.

We finally remark that our scenario C is completely consistent with the PSR observation of massive neutron stars Demorest et al. (2010); Antoniadis et al. (2013), as shown in Fig. IVb. The data for the event observed by GW170817 Abbott et al. (2018, 2017); Evans et al. (2017) ($M_{\rm star}=1.17\sim 1.60\,M_{\odot}$, $R=11.9\pm 1.4\,$km) is also indicated in Figs. IV, and can be reproduced with a pure NM equation of state. Before drawing firm conclusions, however, it is most important to closely investigate the role of strangeness, both in the NM and the QM phase.

V SUMMARY

We have studied the equation of state for cold high density neutron star matter by using the two-flavor NJL model. Our principal aim was to use a model framework which is based on a reasonable description of hadron structure in vacuum and nuclear matter. We used the Gibbs conditions to construct the hadron-quark phase transition. Our emphasis was on the important counteracting roles played by the attractive pairing interaction (coupling constant $G_{S}$) and repulsive vector interaction (coupling constants $G_{\omega}$, $G_{\rho}$) in quark matter. For this purpose, we scaled the value of $G_{S}$ adjusted to the free nucleon mass by a factor $c_{s}$, and $G_{\omega}$, $G_{\rho}$ adjusted to the binding energy and symmetry energy at the saturation point of isospin symmetric nuclear matter by a common factor $c_{v}$, and investigated the dependence of the results on the parameters $(c_{s},c_{v})$. We found that there exists a narrow region of “allowed” values of these parameters, shown in yellow in Fig. 1, which give a reasonable description of the hadron-quark phase transition and the properties of neutron stars. Importantly, this region extends up to $c_{s}=0.98$, i.e., we found that practically the same strength of the scalar diquark interaction, which is required to reproduce the nucleon mass and other quark-quark correlation effects in baryons Cloët et al. (2014), can also be used to describe the phase transition to color superconducting quark matter and the stable massive hybrid stars. Concerning the vector interaction, we found that the coupling constants must be decreased in the quark matter sector by a factor of $0.45\leq c_{v}\leq 0.68$ to obtain a reasonable scenario for the phase transition and neutron stars. Nevertheless, we found that the inclusion of the vector repulsion in quark matter is very important, which is consistent with earlier reported results Baym et al. (2018); Whittenbury et al. (2016); Hell et al. (2013).

Our equation of state above the hadron-quark transition density is close to a phase where chiral symmetry is largely restored and color symmetry is strongly broken, which is consistent with the findings of many previous works based on the 2-flavor picture of quark pairing. We compared our resulting equation of state with the phenomenological parametrizations of Refs. Alford et al. (2013); Ranea-Sandoval et al. (2016), which are based on generic stability criteria for hybrid stars, and found a qualitative and quantitative consistency for the region of our allowed parameters. We found that the stability of massive hybrid stars favors the largest possible values of the pairing strength and the vector couplings in quark matter.101010Interestingly, a similar preference of large coupling constants has also been found in a different treatment of the phase transition based on hadron-quark continuityBaym et al. (2018). For the case of point C in Fig. 1, we found a phase transition to QM at about three times the normal nuclear matter density, and a connected sequence of pure neutron star and hybrid star configurations. The maximum star mass in this case exceeds two solar masses, which is in agreement with recent observations of a massive neutron star Demorest et al. (2010); Antoniadis et al. (2013). Our equation of state for quark matter, obtained by using the allowed region of parameters, is stiff enough to support a heavy hybrid star within a reasonably broad range of stability.

Although we used the Gibbs conditions to describe the hadron-quark phase transition, our results for the allowed range of parameters are very similar to the usual Maxwellian first-order phase transition: The variation of the pressure during the phase transition is small, and therefore the spatial extension of the mixed phase inside the star is small compared to its overall size, resulting in an almost sharp boundary between the quark matter core and the surrounding nuclear matter phase. This scenario suggests that the inclusion of finite size effects in the mixed phase, which are known to work against spatially extended mixed phases Endo (2011); Neumann et al. (2003); Voskresensky et al. (2003); Wu and Shen (2017), will not lead to qualitative changes of the overall physical picture. Quantitatively, however, finite size effects should be investigated in future applications of our model.

The calculations presented in this work should be extended to include the effects of strangeness in both the hadronic and the quark matter phases. For the case of single baryons, recent studies have shown that the properties of hyperons (masses, magnetic moments, and form factors) can be described in this NJL framework Carrillo-Serrano et al. (2016). In the quark matter phase different types of pairings, such as color-flavor locking, should be included, and the role of a chemical potential associated with color neutrality should be taken into account. In our view this work is an important step towards a unified description of single hadrons, nuclear systems, quark matter, and neutron stars in the framework of an effective quark theory of QCD.

Acknowledgements.

T.T. wishes to thank the staff and students of the Department of Physics at Tokai University for their discussions and advice. W.B. acknowledges the hospitality of Argonne National Laboratory where part of this work was performed. The work of I.C. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, contract no. DE-AC02-06CH11357.