# Phenomenological QCD equation of state for massive neutron stars

###### Abstract

We construct an equation of state for massive neutron stars based on quantum chromodynamics phenomenology. Our primary purpose is to delineate the relevant ingredients of equations of state that simultaneously have the required stiffness and satisfy constraints from thermodynamics and causality. These ingredients are: (i) a repulsive density-density interaction, universal for all flavors; (ii) the color-magnetic interaction active from low to high densities; (iii) confining effects, which become increasingly important as the baryon density decreases; (iv) nonperturbative gluons, which are not very sensitive to changes of the quark density. We use the following ”3-window” description: At baryon densities below about twice normal nuclear density, , we use the Akmal-Pandharipande-Ravenhall (APR) equation of state, and at high densities, , we use the three-flavor Nambu-Jona-Lasinio (NJL) model supplemented by vector and diquark interactions. In the transition density region, we smoothly interpolate the hadronic and quark equations of state in the chemical potential-pressure plane. Requiring that the equation of state approach APR at low densities, we find that the quark pressure in nonconfining models can be larger than the hadronic pressure, unlike in conventional equations of state. We show that consistent equations of state of stiffness sufficient to allow massive neutron stars are reasonably tightly constrained, suggesting that gluon dynamics remains nonperturbative even at baryon densities .

###### pacs:

26.60.Kp, 21.65.Mn, 21.65.Qr, 21.65.Cd, 26.60.-c, 05.70.Ce## I Introduction

As the Relativistic Heavy Ion Collider and the Large Hadron Collider continue to push the limits of our experimental knowledge of hot dense quantum chromodynamics (QCD), neutron stars are the only cosmic laboratories in which we can study the structure of cold dense QCD BC1976 (); Lattimer:2006xb (); Fukushima:2010bq (); Alford:2007xm (); Buballa:2014tba (). The recent discoveries of neutron stars with masses ( is the solar mass), including the binary millisecond pulsar J1614-2230 with mass Demorest () and the pulsar J0348+0432 with mass Antoniadis2013 () (also PSR J1311-3430 Romani2012 ()), together with recent simultaneous determinations of neutron star masses and radii Ozel2010 (); Steiner2012 () pose particular challenges to the theoretical construction of the neutron star equation of state.

On the one hand, the existence of these massive stars suggests that the equation of state must be stiffer than conventional hadronic descriptions of matter including hyperons. Furthermore, the central baryon density in neutron stars with masses well exceeds twice nuclear matter density , and may reach as high as . To understand why such high mass stars are stable requires a knowledge of the equation of state at baryon densities over a range . However we cannot at present reliably calculate the equation of state over such a range; the densities are too high to apply reliably conventional hadronic equations of state, and too low to apply perturbative QCD.

This situation motivates us to investigate the properties of strongly correlated quark matter, intermediate between the hadronic and perturbative QCD phases, and ask how the properties of such matter is constrained by neutron star observations. Using a schematic quark model, we manifestly take into account quark degrees of freedom, while including interaction effects such as vector repulsion between quarks, known from hadron spectroscopy, color-magnetic diquark interactions, and six-quark interactions arising from the axial anomaly. We examine the roles of these interactions and find that it is possible, within a reasonable parameter range, to construct an equation of state that (i) is sufficiently stiff to include stable stars with ; (ii) satisfies the thermodynamic constraint that the baryon number density be an increasing function of the baryon chemical potential, ; and (iii) is consistent with the (suggestive) causality constraint that the speed of sound (at zero frequency) not exceed the speed of light ruderman-bludman (); Ellis2007 (). While these conditions provide relatively tight constraints on the quark matter equation of state, it is nonetheless possible to construct the desired equation of state using quark model parameters compatible with the hadron spectroscopy.

To further motivate the picture of strongly correlated quark matter, we briefly review the domain of applicability of hadronic and perturbative QCD equations of state. Conventional hadronic equations of state are constrained by experimental data at low energy and density, e.g., two-body hadronic scattering below the pion production threshold, the masses and level structure of light nuclei, and nuclear matter around nuclear saturation density . While hadronic equations of state include the relevant physics in the low density regime in which their parameters are fit, with increasing density multiple meson exchanges, many-baryon interactions, and virtual baryonic excitations become increasingly important. (The systematics can be most clearly seen in the chiral effective theory approach Weinberg:1978kz (); Epelbaum:2008ga (); Bogner:2009bt ().) In nucleonic potential models Akmal:1998cf (), the three-body nucleon interaction is crucial to reproducing nuclear matter properties at , and its contribution to the energy density can be even comparable (and of opposite sign) to that of the two-body force at . Beyond baryon densities a well defined expansion in terms of static two-, three-, or more, body forces no longer exists.

The equation of state of perturbative QCD Freedman:1976ub (); Fraga:2001id (); Kurkela:2009gj () relies on a picture of weakly coupled quarks and gluons. A current state-of-the-art calculation in this regime, to second order in the strong interaction fine structure constant , with strange quark mass corrections Kurkela:2009gj (), finds a relatively strong dependence of the QCD equation of state on the renormalization scale below the quark chemical potential GeV, which corresponds to a baryon density . Such dependence indicates that nonperturbative effects remain quite important in the lower density range relevant to neutron stars.

In constructing a phenomenological QCD equation of state here, we follow the spirit of the “3-window” approach of Masuda, Hatsuda, and Takatsuka Masuda:2012kf () that interpolates between a nuclear equation of state at low density and a quark equation of state at high density. At densities below we adopt the hadronic Akmal-Pandaripande-Ravenhall (APR) equation of state Akmal:1998cf () (denoted in their paper as A18++UIX*). At densities above 4-7 , where a gas of baryons of radius 0.4-0.5 fm would begin to percolate perc (), we employ a three-flavor Nambu–Jona-Lasinio (NJL) quark model including vector and diquark interactions. While we use the NJL model to be specific, our discussions are more general. In the intermediate region, where purely hadronic or purely quark descriptions are not appropriate, we construct an equation of state using a smooth polynomial interpolation in the baryon chemical potential–pressure (,) plane. In this plane the pressure must be a continuous and monotonically increasing function of . One cannot rule out the possibility of a first order transition as the baryon density increases; such a transition would appear as a discontinuity in the first derivative, .

The 3-window approach is quite different from the conventional hybrid one in which one regards the quark and hadronic phases as distinct. In the latter, the quark pressure at given must, with increasing , intersect the hadronic pressure from below, and moreover must remain larger than the hadronic pressure at larger . By regarding the hadronic equation of state at density larger than as a valid description of matter, such a hybrid construction implicitly selects out possible forms of quark equations of state; in order to have an intersection, the quark pressure must be larger than that of the hadronic phase as large . Such quark equations of state are typically soft, and as a consequence hybrid stars with larger quark cores tend to have smaller masses. In contrast, in the 3-window approach, we construct high density quark equations of state independently of assuming a trustworthy high density hadronic pressure; at high density, the resulting quark pressure at given does not have to grow fast and may remain smaller than the pressure extrapolated from the hadronic phase. Such quark equations of state tend to be stiff, and a star with large quark core can have a large mass.

Within this schematic 3-window description, we aim to incorporate the following effects known from observed hadronic spectroscopy: (i) A repulsive flavor-independent density-density interaction Kunihiro:1991qu (), which stiffens the equation of state (Fig.1a). (ii) The attractive color-magnetic interaction, relevant at all densities, which reduces the average single quark energy (Fig.1.(b)). This effect is similar to that observed in the constituent quark model De Rujula:1975ge (), in which the average quark energy in a nucleon is reduced from the constituent quark mass, MeV, to one-third of the nucleon mass, MeV. As we show, this effect plays an important role in ensuring that the interpolated equation of state satisfies thermodynamic constraints. (iii) Confinement, which, at low densities, traps quarks into baryons and forbids quarks to contribute significantly to the pressure. As we describe, the NJL model, which does not include confinement, has a higher pressure at low density than nuclear models. The requirement that the interpolated pressure merges smoothly into APR at low densities (Fig. 1c) effectively suppresses such excess pressure.

The present approach of interpolating between a hadronic and an NJL based quark picture makes the tacit assumption that the behavior of the gluon sector does not change appreciably over the range ; and in particular, that the gluons do not add a bag constant, , to the energy (and subtracted from the pressure) when the nonperturbative gluons become perturbative. The bag constant measures the energy difference between the trivial (perturbative vacuum) and the nonperturbative vacuum. With the zero-point of the energy set to make the QCD (nonperturbative) vacuum energy zero, the perturbative vacuum has positive energy. Thus, whenever we consider the extreme conditions under which nonperturbative effects disappear, we must include the bag constant in addition to the contributions of the perturbative effects. However, the stability of massive neutron stars does not permit gluon condensation at the QCD scale GeV to produce a gluonic bag constant, , since such a term would, as we argue, too greatly soften the equation of state; we thus exclude the possibility of such a term. On the other hand, a quark bag constant, of order associated with restoration of chiral symmetry, is unavoidable in the NJL model Asakawa:1989bq ().

We extrapolate NJL parameters obtained via hadron phenomenology at to high density quark matter () Hatsuda:1994pi (); Buballa:2003qv (), an approach that is consistent with the observation from analyses for a large number of colors, , that gluon dynamics is insensitive to quark loop effects McLerran:2007qj ().

This paper is organized as follows. In Sec. II, we briefly describe the hadronic and quark models adopted in this study. In Sec. III.1, we examine interaction effects on the quark equation of state. In Sec. IV, we construct the interpolated equation of state. In particular, we explain the difference between the present 3-window description and conventional equations of state which introduce a first order phase transition between hadronic and quark matter in BC1976 (); Alford:2002rj (); Alford:2004pf (); Klahn:2006iw (); Bonanno:2011ch (). As we see the constraints from thermodynamics and causality are quite important. In Sec. V, we solve the Tolman-Oppenheimer-Volkoff (TOV) equation and examine the resulting mass-radius (-) relation of neutron stars. Section VI is devoted to a summary and outlook.

We use the following conventions: , , and the charge conjugation operator is . The flavor and color generator matrices and are composed, respectively, of the identity and Gell-Mann matrices, and are normalized as and . We work in units in which and .

## Ii Models

In this section we briefly summarize the features of the hadronic (APR) and quark (NJL) matter equations of state employed in this paper. APR will be used to describe low density matter, , while the NJL model will be used at high densities, . The precise density beyond which we adopt a fully quark description of matter will depend upon details of the interpolation, as discussed in section IV.

### ii.1 The APR equation of state

In this work we adopt the A18++UIX version of the APR equation of state to describe low density hadronic matter Akmal:1998cf (). This equation of state, based on the Argonne two-body potential, which fits hadronic scattering data very well, and the Urbana IX three-body interaction, which is important to explain nuclear saturation properties, includes charge neutrality and -equilibrium. The indicates the inclusion of relativistic corrections. For simplicity, we adopt the APR version excluding neutral pion condensate, which emerges at . The APR model includes only nucleonic degrees of freedom, and does not take into account hyperons, whose interactions with nucleons and among themselves are not well determined. Typical models of nucleon-hyperon interactions predict hyperon onset at a density . We restrict our application of APR to .

### ii.2 The NJL equation of state

#### ii.2.1 The Lagrangian

In descriptions of quark matter, we adopt a three flavor Nambu–Jona-Lasinio model with Lagrangian density

(1) |

where is a quark field with color, flavor, and Dirac indices, is the quark current mass matrix, and and are four- and six-quark interaction terms, respectively, chosen to reflect the symmetries of QCD. The four-quark interactions possess symmetry for flavors, while the six-quark interactions reflect the anomaly.

The first of the four-quark interactions describes spontaneous chiral symmetry breaking:

(2) |

where (attractive), and is the chiral operator with flavor indices (summed over the color index ).

The second of the four-quark terms Kitazawa:2002bc (),

(3) |

describes the repulsive density-density interaction, analogous to -meson exchange in nuclear matter.

The third of the four-quark terms,

(4) | |||||

describes attractive diquark pairing, where and () are the antisymmetric generators of U(3) flavor and SU(3) color, respectively. The structure of the interaction can be understood as the color-magnetic interaction in the scatterings of quarks in -wave, spin-singlet, flavor- and color- anti-triplet channel. The operators are diquark operators of left- and right-handed chirality.

Next we discuss the six-quark interactions responsible for the anomaly Kobayashi:1970ji (). The first term involves the product of the chiral condensates of different flavors:

(5) |

where denotes the determinant with respect to flavor indices. The second term couples the chiral and diquark condensates Hatsuda:2006ps (),

(6) |

At tree level, these two interactions may be related via a Fierz transformation, which leads to the conclusion . However, renormalization effects will, in general, destroy this equality so that at the mean field level we may treat and as independent parameters, but with .

#### ii.2.2 Electric and color charge neutrality constraints

In order to avoid energetically expensive static long-range electric Coulomb interactions and color flux tube configurations in stable homogeneous quark matter, we impose the local electric and color charge neutrality constraints

(7) |

where is the local electric charge density, and the are the local color densities. These conditions are enforced via standard Lagrange multipliers – with the appropriate chemical potentials coupled to the electric and color charge densities, respectively Iida:2000ha ().

Introducing the charge chemical potential , we add to the Lagrangian the terms

(8) |

where is the quark charge operator for quarks, in units of the proton charge . The are lepton fields and the lepton masses. We may safely omit contributions from - and - leptons since their populations are vanishingly small in the density range of interest Masuda:2012kf ().

Colors and flavors in dense matter are coupled through the diquark interactions of . Thus, an asymmetry in quark flavor densities (e.g., a 2SC phase) leads to a corresponding net quark color density. Most generally case, we should introduce eight independent color chemical potentials Buballa:2005bv (). However, for the diquark pairing structures considered in this paper, all color densities except and automatically vanish Abuki:2005ms (). Thus, we only need to add the terms

(9) |

to constrain the system. The values of will be tuned to satisfy the neutrality conditions.

#### ii.2.3 Mean field equation of state

The mean fields for the chiral condensate and quark densities are

(10) |

Below we write for later convenience. For the diquark mean fields, we write

(11) |

where

(12) |

With these definitions, the diquark condensates correspond to quark pairings, respectively.

The thermodynamic potential may be computed from the mean field particle propagators in terms of these mean fields; the inverse of the propagator , can be read off from the mean field Lagrangian Hatsuda:2006ps (),

(13) |

where the effective mass matrix has diagonal elements

(14) |

while the three diquark pairing amplitudes,

(15) |

and the effective chemical potential matrix,

(16) |

are color- and flavor-dependent.

For each momentum, the inverse propagator is a matrix. There is spin degeneracy, so maximally there are -independent eigenstates in the presence of the flavor and color asymmetry. In the Nambu-Gor’kov formalism, the eigenenergies appear as pairs, . The single particle contribution to the thermodynamic potential is

where ; here is an ultraviolet cutoff. The -dependence is hidden in the eigenvalue . Because Eq. (13) cannot be inverted analytically, the eigenvalues must be computed numerically for each momentum Blaschke:2005uj ().

In order to remove the double-counting of interactions typical of mean-field treatments, we must also include in the thermodynamic potential the terms

(18) | |||||

These terms are positive (), except for the final term.

The quark matter thermodynamic potential is . However, there still remains the nontrivial choice of the “zero” of the thermodyanmic potential. For discussions of neutron star masses, this procedure is extremely important because in general relativitity, the absolute energy density, as in the TOV equation, and not simply its deviation from the QCD vacuum, is physically relevant. Given that the cosmological constant is extremely small compared to the QCD scale, we set the origin of the thermodynamic potential to zero at zero quark density and temperature. Thus in constructing the quark matter equation of state, we will use the renormalized thermodynamic potential

(19) |

which vanishes at .

Finally, the electron contribution to the thermodynamic potential is the standard

(20) |

with , and where we recall that the electron chemical potential is .

Writing the total thermodynamic potential as , the thermodynamic state of the system is determined minimizing the free energy with respect to the seven condensates {,,} under the neutrality conditions

(21) |

which yields the “gap equations.”

(22) |

Below we solve these self-consistent equations using the method outlined in Abuki:2005ms (). Whenever we encounter regions in the solution suggestive of first order phase transitions, we explicitly compare in the relevant phases to determine which local minima is gives the lower free energy. For the ground state we calculate at a nonzero but very small temperature , which makes the numerical calculations faster and more stable.

#### ii.2.4 The NJL parameters

For the model outlined above, we identify two distinct sets of parameters: and . The first set is fixed by matching to QCD vacuum phenomenology. In this work we will use the set by Hatsuda and Kunihiro (HK)Hatsuda:1994pi () (Table 1), which gives the vacuum effective masses for light flavors, MeV, and strange quark, MeV. The second set of parameters does not manifestly affect the quantities in QCD vacuum at the mean field level; we therefore treat them as free parameters, but of the same order of magnitude as the first set, based upon Fierz transformations connecting the associated interaction vertices in the absence of any known anomalous suppression. Briefly, we will investigate , , and and . The choice of these values will be explained in Sec.III.

HK Hatsuda:1994pi () | 631.4 | 5.5 | 135.7 | 1.835 | 9.29 |

RHK Rehberg () | 602.3 | 5.5 | 140.7 | 1.835 | 12.36 |

LKW Lutz1992 () | 750.0 | 3.6 | 87.0 | 1.820 | 8.90 |

## Iii Qualitative effects on the quark equation of state

In this section we examine a number of qualitative effects related to the quark matter equation of state.

### iii.1 When do the quark equations of state become stiff?

We begin our discussion of stiffening of the equation of state by considering a schematic expression for the quark sector equation of state (see also Alford:2004pf ()). For simplicity, we presently consider only matter in a single phase, ignoring any complications arising from phase transitions. In this context, the energy density may be parameterized in terms of the quark density as

(23) |

where with the quark Fermi momentum. The first term is the kinetic energy contribution. The second term contains contributions from both diquark pairing on the Fermi surface () and mass corrections (), where we assume that . (The limit is applicable for all quarks, even strange, at high density, where chiral restoration occurs, and is applicable for u and d quarks at intermediate density.) The third term represents the density-density interaction. The last term is the bag constant . We neglect, for large , the density dependence of the pairing gaps, as well as that of the bag constant.

Differentiating (23) yields the chemical potential , from which the pressure is:

(24) |

The first term is the kinetic pressure, while the remaining terms correspond to corrections arising from the mechanisms discussed above. For given , the pressure becomes large when and . The former condition is met when the quarks interact attractively near the Fermi surface. The latter condition simply expresses the requirement that the density-density interaction should be repulsive in order to stiffen the equation of state. More generally, for stiff equations of state, the coefficients should be negative, while should be positive. Finally, a smaller value of the bag constant also tends to stiffen the equation of state.

### iii.2 Quark and gluon bag constants: and

To consider the impact of the quark and gluon bag constants, we begin by supposing that both the quark and gluon sectors are weakly interacting and that all quark and gluon condensates are vanishingly small. In the absence of perturbative corrections, the equation of state is then

(25) |

where is a function of the number of quark colors and flavors , and the net bag constant is sum of quark and gluon contributions, .

The existence of the bag constant changes the energy-pressure relation from to . Therefore, a smaller enhances at given , stiffening the equation of state. In fact, for a three flavor ideal quark gas with a bag constant, the maximum neutron star mass scales as Book (); Witten:1984rs ()

(26) |

while the corresponding radius scales as

(27) |

Thus, smaller values of give rise to more massive, larger neutron stars.

In the NJL model, the quark bag constant at large density appears automatically when the gap in the Dirac sea is closed through chiral restoration. As a result, its value be computed explicitly as

(28) |

where is the effective mass in the quark energy; becomes the current quark mass () in the perturbative vacuum and the dynamically generated mass () in the chiral symmetry-broken vacuum. In the HK parameter set with vacuum effective masses and , the bag constant is

(29) |

Naively substituting this value into Eq. (26), we obtain a maximum neutron star mass , which less than 1/2 the mass of observed massive stars. This low maximum mass indicates the importance of interaction effects in order to sustain massive neutron stars.

Typical values of the bag constant used previously Book (); Akmal:1998cf () are in the range MeV. The bag constant used in Book () to construct strange quark stars, , is one-fourth of the NJL bag constant, . For a three-flavor free quark gas with , the star mass is relatively large, , but still does not reach ; to do so requires including interactions. Since the value of the bag constant is not precisely known, we employ the NJL value (29) for consistency. As noted, we should also consider a gluon bag constant, when gluons become perturbative at some large quark density. However, because the quark bag constant in the NJL model alone already provides considerable softening of the equation of state, a significant contribution from the gluonic bag constant is unlikely in our equation of state, as we show later in Sec. IV C.

One might argue that considering a gluonic bag constant is unnecessary because the gluons are integrated out in determining the interactions in the NJL model, and thus the quark bag constant already contains the gluonic contributions. This is not quite correct. In the NJL model, the long-range components of the gluons such as those producing confinement are certainly not taken into account; integrating out the long-range components generally produces nonlocal interactions among quarks, which are not present in the NJL model. Therefore, we must consider the contributions from long-range gluons separately, and not simply ignore the gluonic bag constant.

### iii.3 Repulsive density-density interaction:

The repulsive quark vector interaction is inspired by the repulsive density-density interaction in nuclear matter, described, e.g., by omega meson exchange Walecka:1974qa (). Extrapolating the picture of nuclear matter to the strongly correlated quark matter domain, we anticipate that the quark vector coupling is of similar magnitude to the hadronic coupling scale .

Reference Bratovic:2012qs () demonstrated that the vector coupling should be in order to explain the lattice results on the curvature of the chiral restoration line near zero density Kaczmarek:2011zz (). (Note that the coupling constant here corresponds to half that used in Ref. Bratovic:2012qs ().) In the following, we focus mainly on the value .

The inclusion of a vector coupling smooths out chiral symmetry restoration Bratovic:2012qs () because the density-density repulsion forbids a rapid increase of the baryon density, and as a result the chiral transition also does not occur rapidly. Indeed, beyond a particular critical coupling, a first order chiral transition is turned into a smooth crossover.

Intuitively, an increasing repulsive vector force stiffens the equation of state, vs. , as shown in Fig. 2. While the NJL equation of state is considerably softer than APR for small vector couplings, when is sufficiently large the NJL equation of state can achieve a stiffness on par with APR across a wide range of densities. By increasing sufficiently, we can obtain an equation of state within the present framework stiff enough to support neutron star masses .

On the other hand, increasing makes it more difficult to interpolate between the APR and NJL regimes. This challenge is seen in the plots of and in Fig. 3, where for both and , the APR and NJL curves become more widely separated in as increases. One might imagine that the matching could be performed rather simply by allowing a first order phase transition in the interpolated region; however a first order transition, a sudden increase in is simply a kink in vs. , which does not help the interpolation.

A part of the difficulty of interpolating between APR and NJL is that the constituent quark mass for light flavors is MeV, larger than the one-third of the nucleon mass. Accordingly, the curve in the NJL model tends be below that of APR. We next discuss two-body correlations mediated by the color magnetic interaction, which tend to shift the curves to the left, rendering the interpolation procedure more feasible.

### iii.4 Two-body correlations: the color magnetic interaction

At high density, quarks undergo BCS pairing (diquark condensation) as a consequence of the color magnetic interaction. Pairing reduces the energy density by an amount , or equivalently, enhances the pressure by .

We expect correlation effects among the quarks to increase with decreasing density. Eventually three-quark correlations must be dominant in the hadronic phase. One path to three-quark correlations is increasing diquark correlations plus diquark-single quark correlations beyond that described by the standard choice of the diquark coupling = 0.75-1.0 , based on Fierz transformation of the one-gluon exchange type vertex. In this range, we do not find significant effects of on the equation of state in the density range of interest. A diquark mean field appears only when the Fermi surface becomes sufficiently large to overcome chiral symmetry breaking effects. However, diquark correlations, which can exist even without a large Fermi sea, reduce the energy of a pair to less than twice the effective quark mass. To simulate such effects within the present mean field approach, we allow the diquark mean field to reduce the single quark energy at all densities by exploring somewhat larger values = 1.0-1.5 than the standard.

As shown in Sec. III.1, pairing tends to stiffen the equation of state in the high density regime. Figure 4 shows the development of constituent quark masses and mean field pairing gaps; as increases, quark matter first appears in a 2SC phase in which only up and down quarks are paired (, ), and later evolves into a CFL phase in which all three quark flavors pair (). At the 2SC-CFL transition appears to be first order for all NJL parameter sets. However, given that this transition occurs at relatively low density (), the quark model results must be treated with caution.

Two-body correlations are also important at low densities. For example, in the constituent quark model, color magnetic interactions between quarks, in the presence of confinement, reduces the nucleon mass from three times the constituent quark mass, MeV by some 80 MeV to its physical value. Since confinement, by localizing the quarks into a spatial region , increases the quark kinetic energies, as well as adding the energy of color flux tubes – of typical length and energy , where is the string tension – the energy gain from the color magnetic interaction must exceed MeV.

The diquark interaction, , treated in mean field, qualitatively simulates the reduction in the average quark energy at low density that results from pairing effects. As the magnitude of the diquark interaction increases, the curves of the thermodynamic variables as functions of are shifted leftwards to lower chemical potential, as shown in Fig. 5. Thus, by including effects of pairing, one is able to maintain the stiff equation of state produced by a relatively large vector coupling, while at the same time enabling a smooth interpolation between the NJL and APR equations of state for all thermodynamic variables.

Figure 6 demonstrates the impact of pairing on the stiffness of the NJL equation of state. The discontinuous change of at fixed reflects the first order 2SC-CFL phase transition. While for , the equations of state exhibit softening immediately following the 2SC-CFL transition, as the quark density increases further, pairing effects stiffen the equation of state for all NJL parameter sets, relative to the no-pairing case. Moreover, when the pairing is sufficiently strong (), the equation of state is stiffer than without pairing () across the entire density range.

We note that for large , the quark pressure at given exceeds the APR pressure even at very low densities. Taken at face value, this would suggest that even at very low densities the ground state of QCD matter is quark rather than hadronic matter. However, as we discuss in Sec. IV, this high pressure is an unphysical consequence of the NJL model not being confining at low densities.

### iii.5 Chiral-diquark coupling:

We now turn to the axial anomaly-induced coupling between the chiral and diquark condensates. The importance of the anomalous coupling depends on the size of the chiral and diquark condensates. For , this coupling favors the coexistence of chiral and diquark condensates Hatsuda:2006ps (). Thus, as increases from zero the diquark condensate emerges at lower chemical potential and the chiral condensate persists to higher chemical potential.

Figure 7 shows the impact of on the NJL equation of state. We note that while increasing from to slightly stiffens the equation of state, its impact is much smaller than that of or . Since the term in the Lagrangian can be read as a diquark interaction with an interaction strength proportional to the chiral condensate, the effect of can be largely absorbed by the variation of and ; thus in the following we do not study the variation of in detail but take as a canonical value.

## Iv Interpolated EOS

We now discuss constructing an interpolated equation of state, that smoothly joins the equations of state of the low density APR model to the high density NJL model. The first step in defining an interpolation method for joining the hadronic and quark sectors is to determine the “overlap region” in which the two equations of state will be merged. Beginning in the low density hadronic regime described by APR, we expect that as density increases, corrections from many-body forces, hyperon degrees of freedom, multiple meson exchanges and the like, will become important above . Thus, we fix the lower boundary of the interpolation to .

As the density decreases in the quark regime , confining effects, which trap quarks into baryons, become increasingly important. Assuming that the radius of a typical baryon is fm, baryons begin to percolate at around . Thus, we set the upper boundary of the interpolation to , with the precise location determined by the details of the given NJL parameter set being used.

In the intermediate density regime at a given chemical potential, the pressure of the interpolated equation of state should be lower than the NJL pressure extrapolated to lower densities (Fig.8). This follows from the fact that nonconfining models yield excess quark populations at given chemical potential, due to the unphysically small energy cost of having quarks present. In other words, were the confining effects of QCD incorporated in the description of the quark phase, the pressure, especially in dilute region, would be significantly suppressed. This situation is quite analogous to the “semi”-QGP picture for finite temperature QCD Pisarski:2000eq () in which an “overpopulated” quark pressure is suppressed by Polyakov loop effects, until thermal quarks and gluons exhibit quasiparticle behavior at temperatures beyond , where is the pseudocritical temperature for deconfinement Andersen:2011sf ().

The present 3-window description is quite different from the conventional one involving a first order hadron-quark phase transition BC1976 (). In the latter case, the quark pressure at low density of nonconfining models must be smaller than the hadronic one, in order to ensure the intersection of the quark and hadronic pressure at reasonable density. This is achieved either by restricting the quark model parameters or by introducing a bag constant to lower the quark pressure. Such choices, however, generally affect the quark matter equation of state not only in the (presumably unreliable) low density limit, but also in the high density regime in which it should be reliable.

### iv.1 Thermodynamic constraints

Having discussed the qualitative aspects of a hadron-quark interpolation, we now briefly review the thermodynamic constraints imposed on this interpolation which are necessary to ensure that the interpolated equation of state is physical. These constraints are as follows: (i) the pressure must be continuous everywhere; (ii) must be a monotonically increasing function in order to ensure stability of the system with respect to phase separation: . (iii) In addition one physically expects that the speed of sound must be less than the speed of light: ruderman-bludman (); Ellis2007 (). These conditions tightly constrain possible interpolations of the equation of state.

### iv.2 Interpolation method

We now describe a particular method for constructing a phenomenological quark-hadron equation of state. To interpolate in the variables -, we employ a simple polynomial interpolation function, defining an th order polynomial interpolant for the pressure:

(30) |

where and are defined as the points where and . The coefficients are chosen to satisfy the matching conditions at the boundaries of the interpolating interval. At :

(31) |

and at :

(32) |

The number of derivatives that one matches at each boundary is a matter of choice. In general, matching more derivatives results in a smoother interpolation, but at the same time increases the probability of producing unphysical artifacts in the interpolating region (e.g., inflection points in which violate thermodynamic constraints). Here we match up to second order derivatives at each boundary, which ensures that the pressure, number density, and number susceptibility are continuous. Correspondingly, these six boundary conditions require that Eq. (30) has six terms (.

### iv.3 Interpolated EOS

Figure 9 shows the interpolated equation of state , with the interpolation boundaries . For illustration, we consider two NJL parameter sets:

(33) |

with in both cases. For , Set I satisfies all the conditions demanded by the thermodynamic constraints, as we verify shortly. One cannot, however, within the present polynomial interpolation, construct a sensible interpolation for Set II, because at the interpolation boundaries the APR and NJL pressures are rather widely separated in , a possibility noted in Sec. III.4. This wide separation in requires a small slope of the interpolated pressure, but at the same time the slope must be larger than the slope of the APR pressure at the lower boundary, because of the compressibility condition . The interpolated equation of state has a region of , as is clearly seen in the plot of in Fig. 10.

The result presented in Fig. 10 does not preclude constructing a sensible interpolated pressure for Set II. As one sees in Fig. 10, it is possible to join the low and high baryon density curves with a nondecreasing function of . The present exercise shows that the class of interpolating functions for Set II is much more restrictive than for Set I. For example, if there is a first order phase transition between the hadronic and quark regions, the possible density discontinuities are smaller for Set II than Set I. More detailed treatments of the interpolation region are beyond the scope of the present work and we henceforth restrict our consideration to the simple polynomial interpolation, rejecting NJL parameter sets, such as II, incapable of being joined with APR in a thermodynamically consistent manner.

Figure 11 shows the pressure vs. energy density for parameter set I. In this case, the high density equation of state is as stiff as APR extrapolated into the region beyond . From Fig. 12, we observe that the causality constraint is satisfied when the high density boundary, , is varied from to . If we take , however, we find that , a putative violation of causality.

Finally, we consider the possible impact of a gluonic bag constant, , which should be included when the gluon sector becomes perturbative. This contribution reduces the pressure in the quark matter region by %, which makes it extremely difficult to interpolate between the hadronic and quark regimes without violating the condition (cf. Fig. 9). At the same time, the bag constant increases the energy density by , so the resulting equation of state becomes significantly softer. Strictly speaking, even in this situation it would be possible to construct equations of state by increasing and significantly from our current choices; however the current choices for these couplings are already relatively large and it is difficult to identify a mechanism that would significantly increase either coupling in the dense regime. Thus, we conclude that should be very small in the quark matter equation of state, even at ; the gluons remain condensed and nonperturbative.

## V Mass-radius relation

By solving the TOV equation for a given value of the baryon density at the center, we construct a family of stars whose masses and radii are functions of . The - relation for the equation of state with NJL parameters is shown in Fig. 13, where we have taken as in the previous section. To examine effects of the quark bag constant associated with chiral restoration, we also show results for the three-flavor free quark gas with bag constants and . In the NJL model by itself, the quark bag constant is so large that neutron star masses are restricted to . However, we note that as the bag constant decreases, neutron star masses rise, so that models yielding smaller bag constants allow more massive stars.

We see in Fig. 13 that the curves for our interpolated equation of state and APR are quite similar, a not too surprising result given that our chosen NJL parameter set yields an equation of state quite similar to APR. However, while the thermodynamic properties of the two systems are similar, the underlying effective degrees of freedom are quite different. Indeed, APR is well known for its extreme stiffness at high densities, but the effect of hyperonic degrees of freedom in the hadronic sector are expected to reduce this stiffness. On the other hand, our NJL treatment of the quark sector includes strange quarks from the beginning, and is capable of producing a sufficiently stiff equation of state to support stars whose masses exceed Masuda:2012kf ().

Figure 14 shows the - relation for the interpolated equation of state with parameter set I. For most central densities the stellar radius is km, which is compatible with observational data Ozel2010 (); Steiner2012 (). However, the radius at which the mass starts to rise in the - plot is sensitive to the properties of the hadronic equation of state for , with different hadronic equations of state yielding radii differing by up to km Masuda:2012kf ().

## Vi Summary

In this paper we have constructed a phenomenological equation of state over the range of baryon densities by interpolating between the low density hadronic APR equation of state and the high density NJL quark model. In so doing, we explored a number of relevant ingredients of the equation of state needed to realize neutron stars of mass , while satisfying necessary thermodynamic and causality constraints. These requirements constrain the form of the interpolated equation of state and allow one to infer qualitative effects regarding the intermediate density region between the hadronic and quark regimes. The repulsive density-density interaction, color-magnetic interaction, and confining effects play a vital role in determining the structure of the equation of state. A crucial result is that the gluonic bag constant must be small, i.e., the gluons must remain strongly coupled throughout the density region of interest in massive neutron stars; the gluon sector remains nonperturbative even at . One reaches a similar conclusion from studies of quarkyonic matter McLerran:2007qj (); Kojo:2011fh (); Fukushima:2014pha (), in which nonperturbative gluons in quark matter play a crucial role.

We have also emphasized why the three window model of interpolation Masuda:2012kf () is capable of producing stiffer equations of state than conventional hybrid equations of state involving first order phase transitions. As opposed to conventional models, which require the intersection of quark and hadronic curves, we propose that the quark pressure based on nonconfining models need not (and even should not) intersect the hadronic equation of state before the inclusion of confining effects. While this observation is not directly applicable at low density (where we did not use the NJL model), it results in a much wider range of possible high density quark equations of state. As a result, we are able to explore a region of parameter space that has been omitted from prior studies, while still producing a stiff equation of state required to support massive neutron stars.

In this work we have not taken into account possible meson condensed phases, by which we mean condensates in which the order parameter has the quantum number of a mesonic field. Such condensates have been studied in a nuclear context Migdal:1978az (); Baym-picond (); Kaplan:1986yq () as well in quark matter, e.g., inhomogeneous diquark Alford:2007xm () and chiral Buballa:2014tba () condensates. If extant, such exotic phases would likely occur in the the neither purely hadronic nor purely quark density region in which we have interpolated. Thus one cannot directly take over previous results for meson condensates, including the strength, or density discontinuity, of the first order phase transition to the condensed phase. For a given hadronic equation of state at low density and quark equation of state at high density, the strength of such a phase transition is bounded. Although we have, for simplicity, considered only a smooth interpolation scheme, one should more generally allow for such exotic phases; then the smooth curve used in this work would be replaced by one with a small kink, keeping the positive curvature of in the interpolation. We anticipate that even with condensates at intermediate density, it will still be possible to find a reasonable parameter set for the color-magnetic and vector interactions that is compatible with the existence of massive stars. The issue of exotic phases remains open until we can reliably estimate the high density quark equations of state. Further studies are needed to understand better the impact of exotic phases on the in-medium NJL parameters and in turn, their implications for the description of massive neutron stars.

In order to improve the description of the intermediate density matter in neutron stars it is important to further refine our understanding of the hadronic equation of state near . In particular, a more careful assessment of the importance of many-body interactions and the emergence of hyperons is required. Further constraints may be obtained from heavy-ion collisions, including strangeness production Klahn:2006ir (), and lattice QCD calculations of hyperon-nucleon interactions Aoki:2012tk (). Studies of the density dependence of nuclear forces in terms of quarks and gluons play a crucial role in determining when (and why) quarks emerge as the proper degrees of freedom at high density. It would be desirable in the future to extend to finite temperature the present approach to the equation of state to enable us to address dynamical questions such as applications to heavy ion collisions, and neutron star cooling. It is also important to obtain an improved estimate of the quark bag constant, for were it much smaller than the NJL estimate used here, the softening associated with chiral restoration would be significantly reduced and large vector and diquark couplings would not be necessary to obtain a stiff quark matter equation of state within the present context.

## Vii Acknowledgements

Author T.K. thanks Y. Hidaka, A. Ohnishi, and R. Pisarski for enlightening discussions during the Fourth APS-JPS joint meeting HAWAII 2014, and M. Alford, S. Han, and K. Schwenzer for discussions during his visit to Washington University. He is also grateful to D. Blaschke for discussions and his lectures given at the University of Bielefeld in 2013. Authors G.B. and T.K. thank T. Hatsuda and K. Masuda for numerous discussions of the equation of state at intermediate densities. This research was supported in part by NSF Grants PHY09-69790 and PHY13-05891.

## References

- (1) G. Baym and S. A. Chin, Phys. Lett. 62B, 241 (1976).
- (2) Reviewed in J. M. Lattimer and M. Prakash, Phys. Rept. 442 (2007) 109 [astro-ph/0612440].
- (3) K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74 (2011) 014001 [arXiv:1005.4814 [hep-ph]]; K. Fukushima and C. Sasaki, Prog. Part. Nucl. Phys. 72, 99 (2013) [arXiv:1301.6377 [hep-ph]].
- (4) For a review of the color superconductivity, M. G. Alford, A. Schmitt, K. Rajagopal and T. SchÃ¤fer, Rev. Mod. Phys. 80 (2008) 1455 [arXiv:0709.4635 [hep-ph]].
- (5) For a review of the inhomogeneous chiral phases, M. Buballa and S. Carignano, arXiv:1406.1367 [hep-ph].
- (6) P. Demorest, T. Pennucci, S. M. Ransom, R. M. S. Roberts, and J. W. T. Hessels, Nature (London) 467, 1081 (2010).
- (7) J. Antoniadis, et al., Science 340, 1233232 (2013)
- (8) R. W. Romani, A. V. Filippenko, J. M. Silverman, S. B. Cenko, J. Greiner, A. Rau, J. Elliott, and H. J. Pletsch, Astrophys. J. Lett. 760, L36 (2012).
- (9) F. Özel, G. Baym, and T. Güver, Phys. Rev. D 82 010301 (2010).
- (10) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J, 722 33, (2010); Astrophys. J. Lett. 765, L5 (2013).
- (11) M. A. Ruderman and S. A. Bludman, Phys. Rev. 170, 1176 (1968).
- (12) G. F. R. Ellis, R. Maartens, and M. A. H. MacCallum, Gen. Relativ. Gravit. 39, 1651 (2007).
- (13) S. Weinberg, Physica A 96 (1979) 327; ibid. Nucl. Phys. B 363 (1991) 3; ibid. Phys. Lett. B 251 (1990) 288.
- (14) E. Epelbaum, H. W. Hammer and U. G. Meissner, Rev. Mod. Phys. 81 (2009) 1773 [arXiv:0811.1338 [nucl-th]]; E. Epelbaum, Prog. Part. Nucl. Phys. 57 (2006) 654 [nucl-th/0509032].
- (15) S. K. Bogner, R. J. Furnstahl and A. Schwenk, Prog. Part. Nucl. Phys. 65 (2010) 94 [arXiv:0912.3688 [nucl-th]].
- (16) A. Akmal, V. R. Pandharipande and D. G. Ravenhall, Phys. Rev. C 58 (1998) 1804 [nucl-th/9804027].
- (17) B. A. Freedman and L. D. McLerran, Phys. Rev. D 16 (1977) 1169; ibid. 17 (1978) 1109.
- (18) A. Kurkela, P. Romatschke and A. Vuorinen, Phys. Rev. D 81 (2010) 105021 [arXiv:0912.1856 [hep-ph]].
- (19) E. S. Fraga, R. D. Pisarski and J. Schaffner-Bielich, Phys. Rev. D 63 (2001) 121702 [hep-ph/0101143]; E. S. Fraga, A. Kurkela and A. Vuorinen, Astrophys. J. Lett. 781, 2, L25 (2014) [arXiv:1311.5154 [nucl-th]]; A. Kurkela, E. S. Fraga, J. Schaffner-Bielich and A. Vuorinen, Astrophys. J. 789 (2014) 127 [arXiv:1402.6618 [astro-ph.HE]].
- (20) K. Masuda, T. Hatsuda and T. Takatsuka, Astrophys. J. 764, 12 (2013) [arXiv:1205.3621 [nucl-th]]; Prog. Theor. Exp. Phys. (2013) 073D01. [arXiv:1212.6803 [nucl-th]].
- (21) G. Baym, Physica (Amsterdam) 96A, 131 (1979).
- (22) T. Kunihiro, Phys. Lett. B 271 (1991) 395.
- (23) A. De Rujula, H. Georgi and S. L. Glashow, Phys. Rev. D 12 (1975) 147.
- (24) M. Asakawa and K. Yazaki, Nucl. Phys. A 504 (1989) 668.
- (25) T. Hatsuda and T. Kunihiro, Phys. Rept. 247 (1994) 221 [hep-ph/9401310];
- (26) M. Buballa, Phys. Rept. 407 (2005) 205 [hep-ph/0402234].
- (27) L. McLerran and R. D. Pisarski, Nucl. Phys. A 796 (2007) 83 [arXiv:0706.2191 [hep-ph]].
- (28) M. Alford and S. Reddy, Phys. Rev. D 67 (2003) 074024 [nucl-th/0211046].
- (29) M. Alford, M. Braby, M. W. Paris and S. Reddy, Astrophys. J. 629 (2005) 969 [nucl-th/0411016].
- (30) T. Klahn, D. Blaschke, F. Sandin, C. Fuchs, A. Faessler, H. Grigorian, G. Ropke and J. Trumper, Phys. Lett. B 654 (2007) 170 [nucl-th/0609067].
- (31) L. Bonanno and A. Sedrakian, Astron. Astrophys. 539 (2012) A16 [arXiv:1108.0559 [astro-ph.SR]].
- (32) M. Kitazawa, T. Koide, T. Kunihiro and Y. Nemoto, Prog. Theor. Phys. 108 (2002) 929 [hep-ph/0207255, hep-ph/0307278].
- (33) M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 44 (1970) 1422; G. ’t Hooft, Phys. Rev. Lett. 37 (1976) 8; ibid. Phys. Rev. D 14 (1976) 3432 [Erratum-ibid. D 18 (1978) 2199].
- (34) T. Hatsuda, M. Tachibana, N. Yamamoto and G. Baym, Phys. Rev. Lett. 97 (2006) 122001 [hep-ph/0605018]; H. Abuki, G. Baym, T. Hatsuda and N. Yamamoto, Phys. Rev. D 81 (2010) 125010 [arXiv:1003.0408 [hep-ph]]; P. D. Powell and G. Baym, Phys. Rev. D 88 (2013) 1, 014012 [arXiv:1302.0416 [hep-ph]].
- (35) K. Iida and G. Baym, Phys. Rev. D 63 (2001) 074018 [Erratum-ibid. D 66 (2002) 059903] [hep-ph/0011229]; M. Alford and K. Rajagopal, JHEP 0206 (2002) 031 [hep-ph/0204001]; A. W. Steiner, S. Reddy and M. Prakash, Phys. Rev. D 66 (2002) 094007 [hep-ph/0205201].
- (36) M. Buballa and I. A. Shovkovy, Phys. Rev. D 72 (2005) 097501 [hep-ph/0508197].
- (37) H. Abuki and T. Kunihiro, Nucl. Phys. A 768 (2006) 118 [hep-ph/0509172].
- (38) D. Blaschke, S. Fredriksson, H. Grigorian, A. M. Oztas and F. Sandin, Phys. Rev. D 72 (2005) 065020 [hep-ph/0503194].
- (39) P. Rehberg, S. P. Klevansky, and J. Hefner, Phys. Rev. C 53, 410 (1996).
- (40) M. Lutz, S. Klimt, and W. Weise, Nucl. Phys. A 542, 521 (1992).
- (41) E. Witten, Phys. Rev. D 30 (1984) 272.
- (42) Kohsuke Yagi, Tetsuo Hatsuda, and Yasuo Miake, Quark-Gluon Plasma: From Big Bang to Little Bang, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (Cambridge University, Cambridge, England, 2008)
- (43) J. D. Walecka, Annals Phys. 83 (1974) 491; B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16 (1986) 1; ibid., Int. J. Mod. Phys. E 06 (1997) 515 [nucl-th/9701058].
- (44) N. M. Bratovic, T. Hatsuda and W. Weise, Phys. Lett. B 719 (2013) 131 [arXiv:1204.3788 [hep-ph]].
- (45) O. Kaczmarek, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt and W. Soeldner et al., Phys. Rev. D 83 (2011) 014504 [arXiv:1011.3130 [hep-lat]].
- (46) R. D. Pisarski, Phys. Rev. D 62 (2000) 111501 [hep-ph/0006205]; A. Dumitru, Y. Guo, Y. Hidaka, C. P. K. Altes and R. D. Pisarski, Phys. Rev. D 83 (2011) 034022 [arXiv:1011.3820 [hep-ph]]; K. Fukushima, Phys. Lett. B 591 (2004) 277 [hep-ph/0310121].
- (47) J. O. Andersen, L. E. Leganger, M. Strickland and N. Su, JHEP 1108 (2011) 053 [arXiv:1103.2528 [hep-ph]]; J. O. Andersen, S. Mogliacci, N. Su and A. Vuorinen, Phys. Rev. D 87 (2013) 074003 [arXiv:1210.0912 [hep-ph]]; N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, JHEP 1405 (2014) 027 [arXiv:1402.6907 [hep-ph]].
- (48) T. Kojo, Nucl. Phys. A 877 (2012) 70 [arXiv:1106.2187 [hep-ph]].
- (49) K. Fukushima, Nucl. Phys. A931, 257 (2014). arXiv:1408.0547 [hep-ph].
- (50) A. B. Migdal, Rev. Mod. Phys. 50 (1978) 107.
- (51) G. Baym, Phys. Rev. Lett. 30, 1340 (1973); G. Baym and D. Campbell, in Mesons in Nuclei edited by M. Rho and D. Wilkinson (North Holland, Amsterdam, 1978), Vol. III, pp. 1031.
- (52) D. B. Kaplan and A. E. Nelson, Phys. Lett. B 175 (1986) 57.
- (53) T. Klahn, D. Blaschke, S. Typel, E. N. E. van Dalen, A. Faessler, C. Fuchs, T. Gaitanos and H. Grigorian et al., Phys. Rev. C 74 (2006) 035802 [nucl-th/0602038].
- (54) S. Aoki et al. [HAL QCD Collaboration], PTEP (2012) 01A105 [arXiv:1206.5088 [hep-lat]]; T. Doi et al. [HAL QCD Collaboration], Prog. Theor. Phys. 127 (2012) 723 [arXiv:1106.2276 [hep-lat]].