THE QUARKYONIC STAR

Quantum chromodynamics (QCD) is one fundamental theory in the Standard Model in which non-Abelian gauge fields (gluons) and fermions (quarks) interact nonperturbatively and QCD-related phenomena are characterized by an energy scale: ${{\rm{\Lambda }}}_{{\rm{QCD}}}\simeq 0.2\;{\rm{GeV}}$. The QCD vacuum has a rich structure with a variety of condensates, and it dynamically confines any colored excitation. At finite temperature of $T\sim {{\rm{\Lambda }}}_{{\rm{QCD}}}$ and low baryon density of ${n}_{{\rm{B}}}\ll {{\rm{\Lambda }}}_{{\rm{QCD}}}^{3}$, ultrarelativistic nucleus–nucleus collision at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC), supported by the lattice-QCD Monte Carlo calculations, has revealed quantitative properties of a novel strongly correlated state of QCD matter, hadronic gas at $T\lesssim {{\rm{\Lambda }}}_{{\rm{QCD}}}$ and the quark–gluon plasma (QGP) at $T\gtrsim {{\rm{\Lambda }}}_{{\rm{QCD}}}$. In contrast, when $T\ll {{\rm{\Lambda }}}_{{\rm{QCD}}}$ and ${n}_{{\rm{B}}}\sim {{\rm{\Lambda }}}_{{\rm{QCD}}}^{3}$, our knowledge from QCD is severely limited because the first-principle numerical method based on the importance sampling fails at finite baryon density. Therefore, it is important to impose constraints on theoretical uncertainties from the experimental side such as the neutron star observation and the beam-energy scan program at RHIC and future heavy-ion facilities. The purpose of this paper is to import the state-of-the-art idea of a color deconfinement phenomenon at high T and high ${n}_{{\rm{B}}}$, especially a special feature of dense and large-${N}_{{\rm{c}}}$ QCD, into the neutron star phenomenology, which in turn provides us with useful constraints.

In some limiting cases at low and high baryon densities we have a reasonable understanding of cold QCD matter. At ${n}_{{\rm{B}}}=(1\sim 2){n}_{0}$ (where ${n}_{0}\simeq 0.16\;{{\rm{fm}}}^{-3}$ is the nuclear saturation density), on the one hand, one can utilize empirical knowledge from nuclear physics and well-developed theoretical methods to analyze nuclear matter (NM) properties (Weinberg 1990; Epelbaum et al. 2009). Recent developments include not only the chiral perturbation theory (Fiorilla et al. 2012; Krueger et al. 2013) and the chiral effective model (Holt et al. 2014; Lastowiecki et al. 2015) but also the functional renormalization group (Drews & Weise 2015). On the other hand, at asymptotically high density ${n}_{{\rm{B}}}\gtrsim 100{n}_{0}$ (or even smaller densities hopefully), perturbative QCD (pQCD) calculations are validated for the bulk quantities of quark matter (QM) as seen from the convergence of three-loop perturbative results (Freedman & McLerran 1977; Kurkela et al. 2010; Fraga et al. 2015). Besides, pQCD works well to describe the color-superconducting (CSC) states. In particular, the ground state of QCD in the high-density limit has been identified as the color-flavor-locked (CFL) state (Alford et al. 1999b; see also a review by Fukushima & Hatsuda 2011 and references therein).

The most problematic is the QCD matter study in the intermediate-density region: $2{n}_{0}\lt {n}_{{\rm{B}}}\lt 100{n}_{0}$ or $1.1\;{\rm{GeV}}\lt {\mu }_{{\rm{B}}}\lt 3\;{\rm{GeV}}$ in terms of the baryon chemical potential. Neutron stars are unique cosmic laboratories to access such domains experimentally (Buballa et al. 2014). Useful observations include the mass–radius (M–R) relation, the cooling curve (Blaschke et al. 2000, 2004; Tsuruta et al. 2002; Page et al. 2004, 2006, 2011; Shternin et al. 2011), the surface and toroidal magnetic fields (Duncan & Thompson 1992; Cardall et al. 2001; Kitamoto et al. 2014; Olausen & Kaspi 2014), the gravitational waves from the merger of binary neutron stars (Abbott et al. 2009; Hotokezaka et al. 2011, 2013), and so on. In the present work we will specifically pay attention to the mass–radius relation.

The interesting point of the M–R relation is that it has a one-to-one correspondence to the QCD equation of state (EOS) through the Tolman–Oppenheimer–Volkoff (TOV) equation. The overall size of neutron star radii is determined by the EOS around ${n}_{{\rm{B}}}=(1\sim 2){n}_{0}$ (Lattimer & Prakash 2001, 2007), while the typical neutron star masses $M\gtrsim {M}_{\odot }$ (where ${M}_{\odot }$ denotes the solar mass) are largely correlated with the pressure at the central cores with ${n}_{{\rm{B}}}=(2\sim 10){n}_{0}$. According to the discoveries of PSR J1614-2230 (Demorest et al. 2010) and PSR J0348+0432 (Antoniadis et al. 2013), there should be neutron stars whose mass exceeds $2{M}_{\odot }$, indicating that the QCD EOS must be very stiff at ${n}_{{\rm{B}}}\gt 2{n}_{0}$ as compared to what was naïvely considered.

The requirement of stiffness challenges conventional EOSs beyond the nuclear regime. In typical hadronic models the strangeness appears at ${n}_{{\rm{B}}}=(2\sim 3){n}_{0}$ (Glendenning & Schaffner-Bielich 1998, 1999; Tsubakihara et al. 2010), which significantly softens the EOS (i.e., hyperon puzzle). Seminal works (Weissenborn et al. 1990; Nishizaki et al. 2002; Vidana et al. 2011; Yamamoto et al. 2014; Lonardoni et al. 2015) introduced repulsive interactions of strangeness to circumvent the hyperon puzzle, and the results imply that the repulsion not only in two-body force of baryons but also in three- and more-body forces is necessary to pass the $2{M}_{\odot }$ constraint. While many-body forces including hyperons have been poorly constrained, new data on the nucleon–hyperon scattering and hypernuclei at J-PARC (Nagae 2010; Nakazawa & Takahashi 2010; Tamura 2010), as well as the lattice simulations at the physical quark masses (Inoue et al. 2012; Beane et al. 2013), will provide important clues to resolve the hyperon puzzle.

In general, contributions from many-body force of baryons start growing rapidly at ${n}_{{\rm{B}}}\sim 2{n}_{0}$ (Hebeler et al. 2010). From the microscopic point of view strong many-body correlations imply that baryons should exchange many mesons and quarks at ${n}_{{\rm{B}}}\gtrsim 2{n}_{0}$. With more and more exchanged quarks, the identity of isolated baryons should be diminished, and quark degrees of freedom gradually take over physical degrees of freedom. Eventually, baryons with radii 0.5–0.8 fm begin to overlap with each other at ${n}_{{\rm{B}}}\gtrsim (5\sim 10){n}_{0}$, leading to the percolated QM (Baym & Chin 1976). This physics picture contrasts with other crossover scenarios and smoothly fills in a gap between NM and QM with exchanged mesons and quarks. The matter in the confinement/deconfinement crossover domain should inherit properties from both NM and QM. The idea of quarkyonic matter, which was first recognized by McLerran & Pisarski (2007), is an important step toward a correct understanding of the crossover, as we will closely discuss later. We will further push this idea and consider the concrete realization through the construction of an EOS.

In this paper we delineate the crossover scenario bridged by quarkyonic matter. As we emphasize later, the most essential point in this scenario is that there exists an overlapping region that can be described in terms of strongly interacting baryons or strongly interacting quarks equally. For a concrete realization we utilize a Nambu–Jona–Lasinio (NJL) type model including vector and diquark interactions (for a review, see Buballa 2005), which are generally density dependent. Applying the quark–hadron duality picture, we determine the low-density behavior of parameters by fitting the Akmal–Pandaripande–Ravenhall (APR) nuclear EOS (Akmal et al. 1998). The high-density behavior is constrained by confronting the EOS with the $2{M}_{\odot }$ constraint. In this way we cast the M–R relation onto the quark model parameters, and we analyze the transitional behavior with increasing density.

One robust conclusion from our analysis is that the interactions should remain large even at ${n}_{{\rm{B}}}\sim 10{n}_{0}$ to account for the existence of massive neutron stars. One might consider that the asymptotic freedom in QCD would rather validate a weakly interacting picture just in percolated QM, since the typical distance among quarks is small. Such intuition, however, does not always work unless soft gluon contributions are cut off. For instance, in high-T QCD matter magnetic gluons remain unscreened and nonperturbative contributions survive even for asymptotically high temperature (Hietanen et al. 2009). Also, QCD in two dimensions as an asymptotic free theory ('t Hooft 1974) has dense QM in which nonperturbative gluons survive to asymptotically high density, since in spatially one dimension screening effects are not enhanced at finite density (Schon & Thies 2000; Bringoltz 2009; Kojo 2012). In dense QCD, while the hard-dense-loop (HDL) type calculations at weak coupling predict the electric screening mass of the QCD coupling constant times quark chemical potential, $\sim {g}_{s}{\mu }_{{\rm{q}}}$, actually the soft region in the gluon polarization function is highly dependent on the nonperturbative physics near the quark Fermi surface and can deviate from the HDL results at a qualitative level (Rischke et al. 2001; Huang & Shovkovy 2004b; Fukushima 2005; Kojo & Baym 2014). Therefore, the significance of the screening effects is still an open question until first-principle calculation determines the phase structure. In this situation it should be useful to construct an argument in which nonperturbative gluons survive in percolated QM, as discussed in the quarkyonic scenario. In this paper we try to address this issue from the neutron star phenomenology and use the $2{M}_{\odot }$ constraint to discuss why nonperturbative gluons should survive to ${n}_{{\rm{B}}}\sim 10{n}_{0}$.

The treatment of the crossover EOS proposed in the present work has some overlap with preceding works. Attempts to directly interpolate pQCD and nuclear EOS have been made by Freedman & McLerran (1978) and Fraga et al. (2015). Also, Alford et al. (2005) discussed that an EOS for QM with CSC can mimic the APR EOS. In these studies the descriptions of QM are based on the extrapolation of the pQCD picture. On the other hand, recent studies on the crossover scenario have used more model dependent but direct descriptions to conceptualize strongly correlated QM for ${n}_{{\rm{B}}}\gtrsim (3\sim 5){n}_{0}$ and then interpolated the resulting quark EOS and nuclear EOS at ${n}_{{\rm{B}}}\lesssim 2{n}_{0}$ (Masuda et al. 2013a, 2013b, 2015; Alvarez-Castillo et al. 2014; Hell & Weise 2014; Kojo et al. 2015; Kojo 2015). These studies, however, do not manifestly treat microscopic dynamics in the interpolated domain, so one cannot address, for instance, how the strangeness changes in the crossover domain. In the present work we will take a more direct and concrete path to the microscopic description with the interplay among quarks, diquarks, and baryons. We particularly put our emphasis on the role played by diquarks, with which there emerges a natural classification scheme of crossover scenarios, as we see below.

This paper is organized as follows. In Section 2 we classify the crossover scenarios and discuss the underlying physics. In Section 3 we define our model and elucidate how to treat the model parameters. In Section 4 we construct an EOS and extract the density dependence of the model parameters. In Section 5 we discuss the resulting M–R relations, and finally we make concluding remarks in Section 6.

There is no rigorous order parameter for quark deconfinement as soon as dynamical quarks are included in the theory, irrespective of various theoretical efforts (Fukushima 2003). The Polyakov loop serves for an only approximate order parameter at high T, while we have no clue about even an approximate order parameter if the baryon density is high (see also an attempt by Dexheimer & Schramm (2010) to introduce a parameter analogous to the Polyakov loop). Here we discuss the general features of deconfinement crossover not relying on any models.

2.1. Three Scenarios

2.2. Quark–Hadron Continuity

In the low-density region we adopt the EOS of NM, and in this section let us explain how we compute the EOS of QM at high baryon density. We will treat the NJL model within the mean-field approximation (for details see, e.g., Buballa 2005; Klahn et al. 2007; Kojo 2015). We have three gap energies, ${{\rm{\Delta }}}_{{ud}}$, ${{\rm{\Delta }}}_{{ds}}$, and ${{\rm{\Delta }}}_{{su}}$, associated with diquark condensates. Moreover, we should introduce three chiral condensates, $\langle \bar{u}u\rangle$, $\langle \bar{d}d\rangle$, and $\langle \bar{s}s\rangle$, which are all dynamically determined from the gap equations. As we discuss in detail, the vector coupling is an essential ingredient for our prescription, so we need to treat $n=\langle \bar{q}{\gamma }^{0}q\rangle$ as a mean-field variable, as well as three chemical potentials ${\mu }_{e}$, ${\mu }_{3}$, and ${\mu }_{8}$ to impose electric and color charge neutrality. Therefore, in total, there are 10 mean-field variables, and we must solve 10 gap equations simultaneously.

With these mean-field variables (with a notation ${\phi }_{u}=\langle \bar{u}u\rangle$, ${\phi }_{d}=\langle \bar{d}d\rangle$, and ${\phi }_{s}=\langle \bar{s}s\rangle$) we can write the pressure in the following way:

$$\begin{eqnarray}P & = & \displaystyle \frac{1}{8{\pi }^{2}}{\displaystyle \int }_{0}^{{\rm{\Lambda }}}{dp}\;{p}^{2}\displaystyle \sum _{i=1}^{72}{\epsilon }_{i}(p)\\ & & -{G}_{s}({\phi }_{u}^{2}+{\phi }_{d}^{2}+{\phi }_{s}^{2})-4{G}_{d}{\phi }_{u}{\phi }_{d}{\phi }_{s}+{G}_{{\rm{v}}}{n}^{2}\\ & & -\displaystyle \frac{1}{H}({{\rm{\Delta }}}_{{ud}}^{2}+{{\rm{\Delta }}}_{{ds}}^{2}+{{\rm{\Delta }}}_{{su}}^{2})+{P}_{e}+{P}_{\mu },\end{eqnarray} \tag{ 1 }$$

where ${\epsilon }_{i}(p)$ represents the energy dispersion relations obtained from the Dirac Hamiltonian in the Nambu–Gorkov doubled basis. The last terms, P_e and P_μ, represent the pressure contributions from electrons and muons.

In the presence of the mean-field n, the quark chemical potential is renormalized as ${\mu }_{r}=\mu -2{G}_{{\rm{v}}}n$. Then, the quark chemical potential takes the form of a 9 × 9 matrix given by

$$\begin{eqnarray}&&{\boldsymbol{\mu }}={\mu }_{r}{1}_{{\rm{c}}}\otimes {1}_{{\rm{f}}}-{\mu }_{e}{1}_{{\rm{c}}}\otimes {Q}_{e}+({\mu }_{3}{T}_{{\rm{c}}}^{3}+{\mu }_{8}{T}_{{\rm{c}}}^{8})\otimes {1}_{{\rm{f}}},\end{eqnarray} \tag{ 2 }$$

where ${Q}_{e}={\rm{diag}}(2/3,-1/3,-1/3)$ is the charge matrix in flavor space. The explicit form of the gap matrix with ${{\rm{\Delta }}}_{{ud}}$, ${{\rm{\Delta }}}_{{ds}}$, and ${{\rm{\Delta }}}_{{su}}$ can be found in the literature (Fukushima et al. 2005).

For the model parameters, we use the set of Hatsuda & Kunihiro (1994):

$$\begin{eqnarray}{\rm{\Lambda }} & = & 631.4\;{\rm{MeV}},\quad \quad {G}_{s}{{\rm{\Lambda }}}^{2}=3.67,\quad \quad {G}_{d}{{\rm{\Lambda }}}^{5}=-9.29,\\ {m}_{u,d} & = & 5.5\;{\rm{MeV}},\quad \quad {m}_{s}=135.7\;{\rm{MeV}},\end{eqnarray} \tag{ 3 }$$

which are fixed to reproduce hadron phenomenology in the vacuum. The other parameters H and ${G}_{{\rm{v}}}$ are left unfixed, and we will use this freedom to adjust H and ${G}_{{\rm{v}}}$ to realize the quark–hadron duality interpolating the APR and CFL EOSs.

Because the CFL solution does not exist in the density region near NM, which indicates a first-order phase transition as discussed later, a natural candidate for the QM ground state there is the 2SC state. As elucidated comprehensively in Section 2, the diquark condensate $\langle {ud}\rangle$ in the 2SC phase breaks no new symmetry, and so it can coexist in the hadronic phase. This means that the 2SC phase could be smoothly connected to NM and, furthermore, the quarkyonic scenario requires the presence of a duality region where NM and 2SC can represent the same physical system.

It is still nontrivial how to formulate such duality in a practical way. We use the APR hadronic EOS for the baseline, while the 2SC phase has several unconstrained model parameters such as H and ${G}_{{\rm{v}}}$. Our strategy is to make use of these uncertainties positively to reproduce APR in the low-density region in a spirit of the quarkyonic scenario. To this end we treat ${G}_{{\rm{v}}}$ in the 2SC phase as a ${\mu }_{{\rm{B}}}$-dependent control parameter and fit the resulting EOS with APR. Besides, we emphasize that this choice of ${G}_{{\rm{v}}}$ as a control parameter would provide us with convenient and transparent intuition about modified QM EOS and that, at the same time, ${G}_{{\rm{v}}}$ should be a representative of all unknown effects such as confinement near the normal nuclear density. The important point is that in this way, once ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ is determined, we can smoothly extend APR toward QM without any artificial manipulations.

The pressure P in the 2SC phase depends monotonically on ${G}_{{\rm{v}}}$, so that we can easily find ${G}_{{\rm{v}}}$ that reproduces APR at each ${\mu }_{{\rm{B}}}$. One might think that a ${\mu }_{{\rm{B}}}$-dependent ${G}_{{\rm{v}}}$ may change one of the gap equations, $\partial P/\partial n=0$. For technical simplicity we do not change the gap equation but treat n as an internal mean-field variable that takes a different value from ${n}_{{\rm{B}}}$. Thus, to obtain ${n}_{{\rm{B}}}$ correctly, we have to compute $\partial P/\partial {\mu }_{{\rm{B}}}$ including the derivative on ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$.

In order to find a suitable running coupling, we found it necessary to take a diquark coupling larger than 1. Otherwise, the NJL pressure at any positive ${G}_{{\rm{v}}}$ undershoots APR at typical ${\mu }_{{\rm{B}}}$, because the NJL effective quark mass, ${M}_{{\rm{q}}}\simeq 336\;{\rm{MeV}}$, allows the pressure to develop only when ${\mu }_{{\rm{B}}}\gt 3{M}_{{\rm{q}}}\simeq 1010\;{\rm{MeV}}(\gt {m}_{N})$ (m_N: nucleon mass). This situation is cured by including the sufficiently large diquark coupling that effectively reduces the effective quark mass (Kojo et al. 2015). Similar effects can be seen in constituent quark models where the color hyperfine interaction makes the nucleon mass smaller than three times the constituent quark mass. We will be interested in ${G}_{{\rm{v}}}\gtrsim 0.5{G}_{s}$ to achieve sufficient stiffness, and then $H\gtrsim 1.4{G}_{s}$ is required to pass the $2{M}_{\odot }$ and causality constraints. Therefore, we will examine $H=1.5{G}_{s}$ and $H=1.6{G}_{s}$ in the following.

In our fitting procedures one of the questions is how to choose the fitting range. Let us first consider a relatively wide fitting range, ${\mu }_{{\rm{B}}}=(1.02\sim 1.20)\;{\rm{GeV}}$ or ${n}_{{\rm{B}}}\simeq (1.5\sim 3.5){n}_{0}$. (We omitted a very dilute region ${\mu }_{{\rm{B}}}\lesssim 1.02\;{\rm{GeV}}$ because it requires too much fine-tuning in our fit functions.) We show our numerical results with running ${G}_{{\rm{v}}}$ for $H=1.5{G}_{s}$ (lower solid curve) and $H=1.6{G}_{s}$ (upper solid curve) in Figure 1. We emphasize that we did not assume any functional form a priori and ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ shown in Figure 1 results solely from the fit to APR once we make a choice of the diquark coupling H. We vary H to check the sensitivity and will see that this choice is near the upper limit not to violate the causality. It is important to note that the chiral phase transition is a very smooth crossover in the presence of large ${G}_{{\rm{v}}}$, so that this 2SC phase can accommodate both diquark and chiral condensates for any ${\mu }_{{\rm{B}}}$. This is clear in Figure 2, where the constituent quark masses, M_u, M_d, M_s, and the gap energy ${{\rm{\Delta }}}_{{ud}}$ are given as functions of ${\mu }_{{\rm{B}}}$. There is a small discrepancy between M_u and M_d because of the electric charge neutrality condition that breaks isospin symmetry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our fit to APR tells us that ${G}_{{\rm{v}}}$ must be very large at low ${\mu }_{{\rm{B}}}$ and relax to ${G}_{{\rm{v}}}\sim 0.5{G}_{s}$ at large ${\mu }_{{\rm{B}}}$. The importance of ${G}_{{\rm{v}}}$ can be inferred from nuclear physics (e.g., ω-meson exchange), and typical magnitude ${G}_{{\rm{v}}}\sim 0.5{G}_{s}$ in our fit is understandable within the conventional context. On the other hand, the very large value at low ${\mu }_{{\rm{B}}}$ represents something else and requires explanations. Let us recall that using the diquark coupling, we were able to make the onset chemical potential of the NJL pressure close to APR, as required from the quark–hadron duality. But there remains another problem in our fitting procedures: pressure and number density grow much faster in the NJL model than in APR. The obvious reason is the lack of confinement in our model; a quark excites individually, not as a part of a baryon, producing artificial excess of pressure. But large ${G}_{{\rm{v}}}$ can be used to imitate confining effects by tempering the growth of quark number density. This effect should be stronger at lower density; as baryons overlap, ${G}_{{\rm{v}}}$ becomes less responsible for confining effects and should relax to the value expected from the repulsive nuclear forces.

Interestingly, we have found that the best-fit form of ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ is an inverse logarithm for both $H=1.5{G}_{s}$ and $H=1.6{G}_{s}$. Such an inverse logarithmic is quite suggestive because it is consistent with the common form of the running coupling constant at one-loop level. However, the validity of this fitting should be lost at some point of the baryon density. In fact, at sufficiently high baryon density there is a phase transition from the 2SC to CFL phase that softens the EOS. If we use the running ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ determined above, the EOS including the CFL phase does not support the massive neutron star with $M\gtrsim 2{M}_{\odot };$ ${G}_{{\rm{v}}}$ should be $\sim 0.5{G}_{s}$ or greater in the CFL phase. Thus, we need stiffer 2SC pressure with which the EOS can remain stiff enough even after the transition to the CFL phase occurs.

The above determination used fitting over the range ${n}_{{\rm{B}}}\simeq (1.5\sim 3.5){n}_{0}$, but below we shall relax the condition and use only the range ${n}_{{\rm{B}}}\simeq (1.5\sim 2){n}_{0}$, which corresponds to ${\mu }_{{\rm{q}}}=(340\sim 345)\;{\rm{MeV}}$ (i.e., ${\mu }_{{\rm{B}}}=(1.02\sim 1.035)\;{\rm{GeV}}$). Then the 2SC pressure can be stiffer than before in the range of ${n}_{{\rm{B}}}\simeq (2\sim 3.5){n}_{0}$. It is reasonable to stop fitting at ${\mu }_{{\rm{B}}}=1.035\;{\rm{GeV}}$, because the baryon density is ${n}_{{\rm{B}}}\simeq 2{n}_{0}$, beyond which APR is no longer reliable.

To satisfy the boundary conditions, i.e., the smooth connection to APR on the lower-density side and to the CFL phase with ${G}_{{\rm{v}}}\gtrsim 0.5{G}_{s}$ on the higher-density side, we modify ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ from an inverse logarithm to the following form:

$$\begin{eqnarray}&&{G}_{{\rm{v}}}({\mu }_{{\rm{B}}})/{G}_{s}=\displaystyle \frac{a}{\mathrm{log}[({\mu }_{{\rm{B}}}-b)/c]}+d\end{eqnarray} \tag{ 4 }$$

with an offset by d. Once we fix d, we can determine the other three parameters, a, b, and c, using the smooth connection to APR. We changed d to find that the massive neutron star with $M\gt 2{M}_{\odot }$ is impossible with $d\lesssim 0.3$. We shall therefore choose $d=0.4$ throughout this work. The parameters fixed in such a way, for $H=1.5{G}_{s}$ and $1.6{G}_{s}$, respectively, are listed in Table 1, and the corresponding ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ that interpolates between APR and the CFL phase is overlaid by dashed curves in Figure 1.

Table 1.  Parameters for Interpolating ${G}_{{\rm{v}}}({\mu }_{{\rm{B}}})$ between APR and the CFL Phase

H/${G}_{s}$	d	a	b (GeV)	c (GeV)
1.5	0.4	0.05283	0.4049	0.5735
1.6	0.4	0.1127	0.2942	0.6804

Download table as:  ASCIITypeset image

With this running ${G}_{{\rm{v}}}$ we can find the CFL solution and the 2SC phase, and then we can locate a first-order phase transition between them by comparing the pressure. In Figure 3 we show an example for $H=1.5{G}_{{\rm{v}}}$ to find a first-order phase transition at ${\mu }_{{\rm{B}}}=1.31\;{\rm{GeV}}$, where the pressure of the interpolating 2SC phase with running ${G}_{{\rm{v}}}$ and the CFL phase crosses. We make a remark on the connection between APR and the 2SC phase around ${\mu }_{{\rm{B}}}\sim 1\;{\rm{GeV}}$. From Figure 3 one might think that APR has a slightly larger pressure above the fitting region, and so APR would be rather favored. To resolve such confusion, we here again emphasize our picture of the quarkyonic scenario. The change from NM to QM is not any phase transition, but what we assume is a dual regime around ${\mu }_{{\rm{B}}}\sim 1\;{\rm{GeV}}$ in which NM is gradually taken over by QM. In contrast to this smooth crossover from NM to QM, the change from the 2SC phase to the CFL phase is a genuine physical phase transition with different symmetry-breaking patterns. In many model studies including the present work, this phase transition turns out to be of first order.

Zoom In Zoom Out Reset image size

Now that we have the EOS for the whole range of ${\mu }_{{\rm{B}}}$ from NM to CSC, we can compute not only the pressure P but the energy $\varepsilon ={\mu }_{{\rm{B}}}{n}_{{\rm{B}}}-P$ as well. Actually, the relation of P versus ε is essential for the estimation of the neutron star mass. Because ε involves a first derivative in ${n}_{{\rm{B}}}$, its value jumps discontinuously at the first-order phase transition. We can see this behavior in our numerical results shown in Figure 4. It is also clear in Figure 4 that the 2SC part hardly changes with different choices of H. We can explain this minor dependence from the fact that we impose the same boundary condition of APR at lower density. The other boundary condition of the CFL phase side is, on the other hand, loosely constrained by APR, and so there remains H dependence in the CFL part, as is the case in Figure 4. But as we will see below, the value of H is constrained by the causality condition.

Zoom In Zoom Out Reset image size

One can deduce the sound velocity (squared) ${c}_{s}^{2}=\partial P/\partial \varepsilon$ from the slope of curves in Figure 4. We can numerically take the derivative of P with respect to ε and present the obtained ${c}_{s}^{2}$ in Figure 5. It is important to check the causal condition that ${c}_{s}^{2}$ should not exceed unity (i.e., the speed of light in the natural unit). In the present setup, as understood from Figure 5, the causality is not violated.

Zoom In Zoom Out Reset image size

The structure of ${c}_{s}^{2}$ is almost the same irrespective of the choice of H; with increasing ε it monotonically increases in the APR region and then has a peak in the 2SC phase. Once a first-order phase transition to the CFL phase occurs, ${c}_{s}^{2}$ is pushed down to ∼0.5. In the limit of large ε or high baryon density, ${c}_{s}^{2}$ asymptotically approaches ∼0.6, which slightly depends on H. An interesting observation in Figure 5 is that the peak height strongly depends on H, and if H is greater than $\sim 1.6{G}_{s}$, it would go beyond unity, which would violate the causality. Therefore, such a large H is not allowed, and $H=1.6{G}_{s}$ is close to the upper limit for our choice of $d=0.4$ in Equation (4). In this way we find that there is not much uncertainty in the EOS determination after all.

Finally, we comment on the connection of our result to the pQCD EOS. If we extrapolate our $P({\mu }_{{\rm{B}}})$ curve to ${\mu }_{{\rm{B}}}\sim 3\;{\rm{GeV}}$, it undershoots the pQCD pressure. This is not quite surprising: in the pQCD domain, one should expect smaller ${G}_{s}$, ${G}_{{\rm{v}}}$, and H than those used in our analyses, since nonperturbative effects (gluons) are irrelevant by definition of the pQCD domain. By simply implementing the reduction of NJL parameters at large ${\mu }_{{\rm{B}}}$, one can connect our EOS to the pQCD EOS at ${\mu }_{{\rm{B}}}\gtrsim 3\;{\rm{GeV}}$.

To solve the M–R relation from the TOV equation, what we need is the EOS shown in Figure 4. Plugging our EOS into the TOV equation and changing the initial condition that is the central pressure at r = 0, we can get a curve in the plane of the mass and the radius of the neutron star as is just the standard procedure.

We summarize our results in Figure 6. The M–R curves below ${M}_{\odot }$ are essentially determined by the APR EOS up to ${n}_{{\rm{B}}}\sim 2{n}_{0}$. (The tail at large R is determined by the crust EOS, for which we adopt the SLy model by Douchin & Haensel (2001). Then, the 2SC curves start to take off from the APR by degrees. The 2SC curve's having larger radii than the APR originates from a stiffer EOS in the 2SC phase at ${n}_{{\rm{B}}}\gt 2{n}_{0}$. In the vicinity of $M\sim 2{M}_{\odot }$ the curves have a turning point, which signals the phase transition from the 2SC phase to the CFL phase. The maximal mass in our parameterization reaches $\simeq 2.2{M}_{\odot }$.

Zoom In Zoom Out Reset image size

The neutron star radii at the canonical mass $1.4{M}_{\odot }$ typically spread over $R=(9\sim 16)$ km, depending on the EOS for asymmetric NM. The radius of $R\simeq 12$ km in our results belongs to a group of small radii, reflecting that the APR at low density is softer than other typical hadronic EOSs such as the mean-field EOS. The small radii are consistent with recent observation that favors $R=(10\sim 13)$ km (Ozel et al. 2010, 2015; Steiner et al. 2010, 2013), though those analyses still contain several assumptions to be verified. On the other hand, there is also a group insisting that $R\gt 14$ km (Suleimanov et al. 2011), which favors a stiff hadronic EOS at low density (for a hybrid EOS with a stiff hadronic EOS, see, e.g., Bonanno & Sedrakian 2012; Benic et al. 2015). Another hint for the small radii comes from recent quantum Monte Carlo calculations, which indicate that the nuclear EOS at low density is relatively soft, or even slightly softer than the APR (Gandolfi et al. 2012). The heavy-ion data by Danielewicz et al. (2002) also favor a soft EOS at low density.

Our crossover EOS yields relatively simple topology for the M–R curves, similar to those in the purely hadronic EOS. In case of the hybrid EOS with a first-order phase transition, there are more varieties of the M–R curves. The classification was done by Alford et al. (2013) and Alford & Han (2015), assuming the constant speed of sound.

Quarkyonic matter is a likely candidate to bridge a gap between NM and QM. It is an intermediate state of matter with a duality region that is describable by baryons or quarks. Thus, the EOS construction is not a patchwork of distinct EOSs of NM and QM, as is a conventional hybrid description having first-order phase transitions.

We discussed the interplay among quarks, diquarks, and baryons. Through processes of exchanging mesons or quarks between baryons, the diquark correlation develops in the intermediate state. While the density is low, the diquark is combined with other colored excitations to form local and color-singlet objects. For this reason diquarks should become manifest increasingly as baryons involve more and more quark exchange. We classified scenarios of deconfinement crossover according to the diquark characters and addressed what quarkyonic matter implies.

To interpret the crossover EOS in the language of microscopic degrees of freedom, we introduced a running vector coupling for quark models, with which the quark–hadron duality is implemented in a practically feasible manner. The running coupling is constrained by the conventional nuclear EOS on the low-density side and by massive neutron stars with $M\gt 2{M}_{\odot }$ on the high-density side. Because the idea is quite generic, one can employ any quark models, and we adopted an NJL-type model for a concrete demonstration in this work.

The important finding is that the model parameters around ${n}_{{\rm{B}}}\sim 10{n}_{0}$ must be comparable to those in the vacuum (i.e., matter at $\sim 10{n}_{0}$ still belongs to a nonperturbative regime); otherwise, the $2{M}_{\odot }$ constraint is not satisfied. The upper bound of the diquark interaction comes from the causality condition. After all, there is not much uncertainty remaining in the parameter determination. In principle, all the interactions in the quark model are mediated by gluons, and so substantially large couplings require nonperturbative gluon dynamics. We need to clarify the origin of nonperturbative dynamics in the high-density limit in a similar sense to the high-temperature limit at which the magnetic sector is still confining.

In this paper we interpret the enhancement of many-quark (meson) exchanges at ${n}_{{\rm{B}}}\gtrsim 2{n}_{0}$ as an indicator for a new state of matter. While we consider diquarks in this paper, another branch of our picture would be the condensation of mesonic objects since many mesons are available among baryons. Although the previous studies in hadronic models typically predicted strong first-order phase transitions associated with meson condensed phases, recent studies in quark models indicate that the appearance of inhomogeneous mesonic condensates does not necessarily accompany strong first-order phase transitions and thus can avoid the significant softening (Buballa & Carignano 2015). The Skyrme crystal descriptions at high density (Paeng et al. 2015) have some overlap with such quark descriptions since both of them can address several effects beyond purely hadronic descriptions, such as the structural change of hadrons.

In this work we focused on the bulk properties of matter, i.e., the EOS only. To explore the refined characters of matter and to promote the advantage of our proposed method, we should study transport processes such as the cooling problem involving the strangeness. We shall leave this for future works.

We thank G. Baym, K. Masuda, L. McLerran, and J. Wambach for insightful conversations during Quark Matter 2015 in Kobe. This work was partially supported by JSPS KAKENHI grant Nos. 15H03652 and 15K13479 (K.F.) and by NSF grants PHY09-69790 and PHY13-05891 (T.K.).