### Unified Interacting Quark Matter and its Astrophysical Implications

#### Abstract

We investigate interacting quark matter (IQM), including the perturbative QCD correction and color superconductivity, for both up-down quark matter ($$ud$$​QM) and strange quark matter (SQM). We first derive an equation of state (EOS) unifying all cases by a simple reparametrization and rescaling, through which we manage to maximally reduce the number of degrees of freedom. We find, in contrast to the conventional EOS $$p=1/3(\rho-4B_{\rm eff})$$​ for non-interacting quark matter, that taking the extreme strongly interacting limit on the unified IQM EOS gives $$p=\rho-2B_{\rm eff}$$​, where $$B_{\rm eff}$$​ is the effective bag constant. We employ the unified EOS to explore the properties of pure interacting quark stars (IQSs) composed of IQM. We describe how recent astrophysical observations, such as the pulsar-mass measurements, the NICER analysis, and the binary merger gravitational-wave events GW170817, GW190425, and GW190814, further constrain the parameter space. An upper bound for the maximum allowed mass of IQSs is found to be $$M_{\rm TOV}\lesssim 3.23\, M_{\odot}$$​. Our analysis indicates a new possibility that the currently observed compact stars, including the recently reported GW190814’s secondary component ($$M=2.59^{+0.08}_{-0.09}\, M_{\odot}$$​), can be quark stars composed of interacting up-down quark matter.

# 1 Introduction

It has been long expected that quark matter, a state consisting purely of quark and gluon degrees of freedom without confining into individual nucleons, can form at high density or high temperature. On the other hand, Bodmer [1], Witten [2] and Terazawa [3] proposed that quark matter with comparable numbers of $$u, \,d, \,s$$​ quarks, also called strange quark matter (SQM), might be the ground state of baryonic matter at the zero temperature and pressure.

However, the original proposals were based on the bag model that failed to account for the flavor-dependent feedback of the quark gas on the QCD vacuum. With this being adequately included in a phenomenological quark-meson model, a recent study [4] demonstrated that $$u, d$$​ quark matter ($$ud$$​QM) is in general more stable than SQM, and it can be more stable than the ordinary nuclear matter at sufficiently large baryon number beyond the periodic table. This has been connected to a series of recent phenomenological explorations [5]–[13] and experimental searches [14]–[16]. Interacting quark matter (IQM) includes the interquark effects from perturbative QCD (pQCD) and color superconductivity. The pQCD corrections are due to the gluon-mediated interaction [17]–[19]. Color superconductivity is caused by the spin-0 Cooper-pair condensate that is antisymmetric in color-flavor space, as expected to form to lower the energy [20]–[22]. This can result in two-flavor color superconductivity, where $$u$$​ quarks pair with $$d$$​ quarks (conventionally termed “2SC" /”2SC+s" without/with strange quarks)1, or the color-flavor locking (CFL) phase in which $$u,d,s$$​ quarks pair with each other antisymmetrically.

All previous related studies of color superconductivity assumed an effective bag constant that is independent of the flavor composition, resulting in the conclusion that 2SC phases are absent in compact star physics [23]. However, improved NJL and quark-meson models [4], [24]–[27] suggest that the value of the effective bag constant is very likely dependent on the flavor composition, opening up new possibilities. A relatively smaller bag constant can make the two-flavor color superconductivity stable, and thus we should re-consider its possibility.

The binary merger of compact stars produces strong gravitational wave fields, the waveforms of which encode information about tidal deformation that is sensitive to the matter equation of state (EOS). In general, stars with stiff EOSs can easily be tidally deformed due to their large radii. The GW170817 event detected by LIGO [28], [29] is the first confirmed merger event of compact stars. In conjunction with the more recent GW190425 event [30], this has inspired many investigations into equation of state and the gravitational properties of nuclear matter [31]–[41], SQM [42], [43], and $$ud$$​QM [5], [11]. In particular, it has recently been shown [5], [11] that all compact stars observed are likely to be up-down quark stars ($$ud$$​QSs) composed of non-interacting $$ud$$​QM, taking into consideration considering GW170817, GW190425, and pulsar observations. Alternative explanations for observed glitches and quasi-periodic oscillations might exist for quark stars having possible complicated structures of quark matter.

More recently, a binary merger event GW190814 was reported [44], featuring a primary black hole with mass $$23.2^{+1.1}_{+1.0} M_{\odot}$$​, and a secondary companion of $$2.59^{+0.08}_{-0.09}\, M_{\odot}$$​, which is much larger than the upper bound $$M_{\rm TOV}\lesssim 2.3 M_{\odot}$$​ of the maximum mass of a non-rotating neutron star, set by various analyses of GW170817 [45]–[48]. Conventional non-interacting SQM and $$ud$$​QM have bag constant values not sufficiently small to account for this large star mass [5], [11], [42]. It is consequently of interest to see how strongly interacting quark stars (IQSs), which are composed of IQM, can fit all these constraints.

We begin by first providing a unified framework for all possible strongly interacting phases of SQM and $$ud$$​QM by a simple reparametrization, and derive a simple but universal EOS for IQM. We then explore the properties of IQSs, investigating how astrophysical constraints such as observed large pulsar masses [49]–[51], recent analysis of the NICER X-ray spectral-timing event data [52], [53], and the LIGO GW170817, GW190425, GW190814 events [28]–[30], [44] constrain the IQS parameter space. It turns out all these constraints are compatible with compact stars composed of strongly interacting up-down quark matter for a determined parameter space.

# 2 Properties of IQM

We first rewrite the free energy of the superconducting quark matter [23] into a general form with the pQCD correction included [18], [54], [55]:

$\begin{gather} \Omega=&-\frac{\xi_4}{4\pi^2}\mu^4+\frac{\xi_4(1-a_4)}{4\pi^2}\mu^4- \frac{ \xi_{2a} \Delta^2-\xi_{2b} m_s^2}{\pi^2} \mu^2 \\ &-\frac{\mu_{e}^4}{12 \pi^2}+B_{\rm eff} , \label{omega95mu} \end{gather}$ where $$\mu$$​ and $$\mu_e$$​ are the respective average quark and electron chemical potentials2. The first term represents the unpaired free quark gas contribution. The second term with $$(1-a_4)$$​ represents the pQCD contribution from one-gluon exchange for gluon interaction to $$O(\alpha_s^2)$$​ order. To phenomenologically account for higher-order contributions, we can vary $$a_4$$​ from $$a_4=1$$​, corresponding to a vanishing pQCD correction, to very small values where these corrections become large [18], [54], [55]. The term with $$m_s$$​ accounts for the correction from the finite strange quark mass if applicable, while the term with the gap parameter $$\Delta$$​ represents the contribution from color superconductivity. The constant coefficients are

$\begin{gather} (\xi_4,\xi_{2a}, \xi_{2b}) = \left\{ \begin{array} {ll} ( \left(\frac{1}{3}\right)^{\frac{4}{3}}+ \left(\frac{2}{3}\right)^{\frac{4}{3}})^{-3},1,0) & \textrm{2SC phase}\\ (3,1,3/4) & \textrm{2SC+s phase}\\ (3,3,3/4)& \textrm{CFL phase} \end{array} \nonumber \right.\end{gather}$ for the various types of quark matter. $$B_{\rm eff}$$​ is the effective bag constant that accounts for the non-perturbative contribution from QCD vacuum, the size of which can be flavor-dependent [4].

From the thermodynamic relations

$$$p=-\Omega, \,\, n_{q}=-\frac{\partial\Omega}{\partial \mu},\,\, n_{e}=-\frac{\partial\Omega}{\partial \mu_e},\,\, \rho=\Omega+n_q \mu+n_e \mu_e , \label{thermo}$$$ we obtain the relevant thermodynamic quantities such as pressure $$p$$​, quark and electron number density $$n_{q,e}$$​, and total energy density $$\rho$$​, in terms of the chemical potential, where $$\Omega$$​ is the free energy. To reduce the size of the parameter space, we define

$$$\lambda=\frac{\xi_{2a} \Delta^2-\xi_{2b} m_s^2}{\sqrt{\xi_4 a_4}} \label{lam}$$$ characterizing the strength of the related strong interaction. The relations ($$\ref{thermo}$$) then imply

$\begin{gather} n_q&=&\frac{\xi_4 a_4}{\pi^2}\mu^3 + \frac{\lambda \sqrt{\xi_4 a_4} }{\pi^2} 2\mu, \quad n_e=\frac{\mu_e^3}{3\pi^2}, \label{n95mu}\end{gather}$

and

$\begin{gather} \rho&=&\frac{3\xi_4 a_4}{4\pi^2}\mu^4+\frac{\mu_{e}^4}{4 \pi^2}+B_{\rm eff} + \frac{ \lambda\sqrt{\xi_4 a_4} }{\pi^2} \mu^2. \label{rho95mu}\end{gather}$

Using Eq. ($$\ref{omega95mu}$$) and Eq. ($$\ref{rho95mu}$$) to eliminate the dependence on $$\mu$$​ we obtain the following analytic expression

$$$p=\frac{1}{3}(\rho-4B_{\rm eff})+ \frac{4\lambda^2}{9\pi^2}\left(-1+\sqrt{1+3\pi^2 \frac{(\rho-B_{\rm eff})}{\lambda^2}}\right) \label{eos95tot}$$$ for the unified IQM EOS 3. We see that the general EOS expression above unifies the 2SC, "2SC+s", and CFL phases. It only has 2 independent parameters $$(B_{\rm eff},\lambda)$$​, whilst all other parameters $$(a_4, \Delta,m_s)$$​ and $$(\xi_4,\xi_{2a}, \xi_{2b})$$​ are subsumed in $$\lambda$$​ using Eq. ($$\ref{lam}$$)4 .

We can further remove the $$B_{\rm eff}$$​ parameter by doing the dimensionless rescaling:

$$$\bar{\rho}=\frac{\rho}{4\,B_{\rm eff}}, \,\, \bar{p}=\frac{p}{4\,B_{\rm eff}}, \,\, \label{scaling95prho}$$$ and

$$$\bar{\lambda}=\frac{\lambda^2}{4B_{\rm eff}}= \frac{(\xi_{2a} \Delta^2-\xi_{2b} m_s^2)^2}{4\,B_{\rm eff}\xi_4 a_4}, \label{scaling95lam}$$$

so that the equation of state Eq. ($$\ref{eos95tot}$$) reduces to the dimensionless form

$$$\bar{p}=\frac{1}{3}(\bar{\rho}-1)+ \frac{4}{9\pi^2}\bar{\lambda} \left(-1+\sqrt{1+\frac{3\pi^2}{\bar{\lambda}} {(\bar{\rho}-\frac{1}{4})}}\right), \label{eos95p}$$$ or conversely

$$$\bar{\rho}=3\bar{p}+1- \frac{4}{\pi^2}\bar{\lambda} \left(-1+\sqrt{1+\frac{\pi^2}{\bar{\lambda}} {(\bar{p}+\frac{1}{4})}}\right). \label{eos95rho}$$$ With Eq. ($$\ref{eos95p}$$), one can easily show that $$\partial P/\partial \bar{\lambda}$$​ is always positive, so a larger $$\bar{\lambda}$$​ (i.e., smaller $$B_{\rm eff}, a_4, m_s$$​ or larger $$\Delta$$​) makes the EOS stiffer, resulting in a larger star mass and radius, as we will subsequently demonstrate numerically.

As $$\bar{\lambda}\to0$$​, Eq. ($$\ref{eos95p}$$) reduces to the conventional non-interacting rescaled quark matter EOS $$\bar{p}=(\bar{\rho}-1)/3$$​. At the opposite limit where $$\bar{\lambda}$$​ goes extremely large, Eq. ($$\ref{eos95p}$$) approaches the special form

$$$\bar{p}\vert_{\bar{\lambda}\to \infty}=\bar{\rho}-\frac{1}{2}, \label{eos95infty}$$$ which is equivalent to $$p={\rho}-2B_{\rm eff}$$​ after scaling back with Eq. ($$\ref{scaling95prho}$$). We see that strong interaction effects can reduce the surface mass density of the quark star from $$\rho_0= 4B_{\rm eff}$$​ down to $$\rho_0=2B_{\rm eff}$$​, and increase the quark matter sound speed $$c_s^2=\partial p/\partial \rho$$​ from $$1/3$$​ up to $$1$$​ (the light speed) maximally.

Since the energy per baryon number $$E/A=3\mu\vert_{P=0}=3\mu\vert_{\Omega=0}$$​, we have

$$$\frac{E}{A}=\frac{3\sqrt{2} \pi}{(\xi_4 a_4)^{1/4}}\frac{ {B_{\rm eff}}^{1/4}}{\sqrt{\sqrt{4\bar{\lambda}+\pi^2}+2\sqrt{\bar{\lambda}}}} \label{EA}$$$ from ($$\ref{omega95mu}$$). We see that a smaller bag constant or a larger $$\bar{\lambda}$$​ yields a smaller $$E/A$$​. As $$\bar{\lambda} \to 0, \,a_4\to1$$​, we recover the results for non-interacting quark matter [4], [57].

# 3 properties of IQS$$s$$​

To study gravitational effects of interacting quark stars, we further rescale the mass and radius into dimensionless form [58], [59] in geometric units where $$c=G=1$$

$$$\bar{m}=m{\sqrt{4\,B_{\rm eff}}}, \quad \bar{r}={r}{\sqrt{4\,B_{\rm eff}}}, \,\, \label{scaling95mr}$$$ so that the Tolman-Oppenheimer-Volkov (TOV) equation [60], [61]

\begin{aligned} {dp(r)\over dr}&=-{\left[m(r)+4\pi r^3p(r)\right]\left[\rho(r)+p(r)\right]\over r(r-2m(r))}\,,\,\,\\ {dm(r)\over dr}&=4\pi\rho(r)r^2,\, \end{aligned} \label{eq:tov} can also be rescaled into dimensionless form (simply replace non-bar symbols with barred ones). The rescaled TOV solution on $$(\bar{M}, \bar{R}) =(M\sqrt{4B_{\rm eff}}, R\sqrt{4B_{\rm eff}})$$​ can thus be obtained with the rescaled EOS Eq. ($$\ref{eos95p}$$) with respect to any given value of $$\bar{\lambda}$$​. We depict the solutions to Eq. $$\ref{eq:tov}$$ in Figure 1 for the rescaled mass and radius for several different values of $$\bar{\lambda}$$​. The physical $$(M, R)$$​ with respect to any specific $$B_{\rm eff}$$​ value can then straightforwardly be obtained directly from rescaling $$(\bar{M}, \bar{R})$$​ back with Eq. ($$\ref{scaling95mr}$$).

From Fig. 1, we easily see that a larger $$\bar{\lambda}$$​ leads to a larger $$M_{\rm TOV}$$​ as expected, since larger $$\bar{\lambda}$$​ maps to a stiffer EOS. We can interpolate the $$\bar{M}_{\rm TOV}(\bar{\lambda})$$​ numerical results with the following sigmoid-type function

$$$\bar{M}_{\rm TOV}(\bar{\lambda})= \frac{\bar{M}_{\infty} }{1+c_1\, e^{-\bar{\lambda}^{c_2}}}+ \left(\bar{M}_0- \frac{\bar{M}_{\infty}}{1+c_1\, e^{-\bar{\lambda}^{c_3}}} \right)e^{-\bar{\lambda}^{c_4}} \label{Mbarmax}$$$ where the coefficients $$c_1\approx 0.8220 , c_2\approx0.4537, c_3\approx0.3313,c_4\approx0.2676$$​ are the best-fit values, with error only at $$0.1\%$$​ level. And $$\bar{M}_{0}=\bar{M}_{\rm TOV}(\bar{\lambda}\to 0)=0.05168$$​, $$\bar{M}_{\infty}=\bar{M}_{\rm TOV}(\bar{\lambda}\to\infty)=0.1204$$​, corresponding to

$\begin{gather} M_{\rm TOV}(\bar{\lambda}\to 0)&\approx&\frac{15.17\, M_{\odot}}{ (B_{\rm eff}/\rm MeV\, fm^{-3})^{1/2}},\\ M_{\rm TOV}(\bar{\lambda}\to \infty)&\approx&\frac{35.35\, M_{\odot}}{(B_{\rm eff}/\rm MeV\, fm^{-3})^{1/2}}\end{gather}$

referring to Eq. ($$\ref{scaling95mr}$$). We see that the strongly interacting limit has a maximum star mass 2.33 times larger than that of the non-interacting case. In general, we have the function $$M_{\rm TOV} (B_{\rm eff}, \bar{\lambda})$$​ from Eq. ($$\ref{Mbarmax}$$) after inserting Eq. ($$\ref{scaling95mr}$$). The largest measured pulsar mass therefore imposes a constraint on the $$(B_{\rm eff}, \bar{\lambda})$$​ space.

Recent NICER analysis of PSR J0030+0451 [52], [53] points to a star with a mass around $$1.4\,M_{\odot}$$​ with a radius around 13 km ($$90\%$$​ C.L.). The inferred contour of joint probability density distribution on the $$M-R$$​ plane can then be translated to a range for constraints on the $$(B_{\rm eff}, \bar{\lambda})$$​ space, utilizing the derived $$\bar{M}(\bar{R})$$​ results presented in Figure 1.

# 4 Tidal deformability

The response of compact stars to external disturbances is characterized by the Love number $$k_2$$​ [62]–[65]

\begin{aligned} k_2 &= \frac{8 C^5}{5} (1-2C)^2 [ 2 + 2 C (y_R-1) -y_R ] \\ &\times \{ 2 C [ 6- 3y_R + 3C (5y_R-8)]+4C^3[ 13-11 y_R \\&+ C (3y_R-2)+2C^2(1+y_R)] \\ &+ 3 (1-2C)^2 [ 2 - y_R + 2C (y_R-1)] \log (1-2C )\}^{-1} , \end{aligned}\label{eqn:k2} where $$C=M/R=C(\bar{M})$$​ for given $$\bar{\lambda}$$​. The quantity $$y_R$$​ is $$y(r)$$​ evaluated at the star surface, which can be obtained by solving the following equation [65]:

\begin{aligned} &ry^\prime(r)+y(r)^2+ r^2Q(r) \\ &+y(r)e^{\lambda(r)}\left[1+4\pi r^2(p(r)-\rho(r))\right]=0\,, \end{aligned}\label{eqn:y} with the boundary condition $$y(0)=2$$​. Here

\begin{aligned} Q(r)&=4\pi e^{\lambda(r)} \left(5\rho(r)+9p(r)+\frac{\rho(r)+p(r)}{c_s^2(r)} \right) \\ &-6\frac{e^{\lambda(r)}}{r^2}-\left(\nu^\prime(r)\right)^2, \end{aligned}\label{eq:Q}

$$$e^{\lambda(r)}=\left[1-{2m(r)\over r}\right]^{-1} \quad \nu^\prime(r)=2e^{\lambda(r)}{m(r)+4\pi p(r)r^3\over r^2} \label{eq:met}$$$ and $$c_s^2(r)\equiv dp/d\rho$$​ denotes the sound speed squared. For stars with a finite surface density like quark stars, a matching condition [66], [67] $$y_R^{\rm ext}=y_R^{\rm int}- 4\pi R^3\rho_s/M=y_R^{\rm int}- 4\pi \bar{R}^3\bar{\rho}_s/\bar{M}$$​ should be used at the boundary. Solving ($$\ref{eqn:y}$$) with $$\rho(r)$$​ and $$p(r)$$​ obtained from ($$\ref{eq:tov}$$), we obtain the function $$k_2(C)$$​ for a given $$\bar{\lambda}$$​. The dimensionless tidal deformability $$\Lambda=2k_2/(3C^5)$$​ as a function of $$\bar{M}$$​ and $$\bar{\lambda}$$​ is thus obtained accordingly.

We depict the results of $$\Lambda(\bar{M},\bar{\lambda})$$​ in Figure 2. Note that the lower end of each curve is determined by requiring the star mass not to exceed its maximum allowed mass. We can see that a larger $$\bar{\lambda}$$​ results in a larger tidal deformability for a given $$\bar{M}$$​, as expected since a larger $$\bar{\lambda}$$​ maps to a stiffer EOS.

Assuming the compact objects detected by the recent LIGO event are pure IQSs, we can use the LIGO constraint on $$\Lambda(1.4 M_\odot)$$​ to narrow down the parameter space $$(B_{\rm eff},\bar{\lambda})$$​ with the $$\Lambda(\bar{M},\bar{\lambda})$$​ calculated above.

The average tidal deformability of a binary system is defined as

$\begin{gather} \tilde{\Lambda} &=& \frac{16}{13} \frac{ (1+12q)}{(1+q)^5} {\Lambda} (M_1)+ \frac{16}{13} \frac{q^4 (12+q)}{(1+q)^5}{\Lambda}({M}_2), \label{LamLam}\end{gather}$ where $$M_1$$​ and $$M_2$$​ are the masses of the binary components, and $$q=M_2/M_1=\bar{M}_2/\bar{M}_1$$​, with $$M_2$$​ being the smaller mass so that $$0<q\lesssim1$$​. For any given chirp mass $$M_{c} = (M_1 M_2)^{3/5}/(M_1+M_2)^{1/5}$$​, one thus has

$$$M_2=(q^2(q+1))^{1/5} M_c \text{ and } M_1=((1+q)/q^3)^{1/5} M_c\; .$$$

Using the rescaled mass parameter $$\bar{M}=M\sqrt{4 B_{\rm eff}}$$​, we eventually obtain $$\tilde{\Lambda}=\tilde{\Lambda}(\bar{M}_1,\bar{M}_2,\bar{\lambda})=\tilde{\Lambda}(M_c,q, B_{\rm eff},\bar{\lambda})$$​. LIGO constraints on $$\tilde{\Lambda}$$​ can then be used to narrow down the parameter space $$(B_{\rm eff},\bar{\lambda})$$​ for given $$M_c$$​ and $$q$$​, assuming the objects detected by the recent LIGO events are pure IQSs.

We thus obtain the constrained parameter space shown in Figure 3, considering related astrophysical observations. We do not use the constraints that assumed hadronic EOSs, since we are studying quark stars. In our results, only the stability lines have an explicit dependence on the flavor composition and the size of $$a_4$$​.

As described previously, the solid black curves are determined from Eq. ($$\ref{Mbarmax}$$) for the measured pulsar masses. The blue and red-colored bands (bounded by dot-dash lines) represent the $$\tilde{\Lambda}$$​ constraint of GW170817 and GW190425 translated from the previously obtained $$\tilde{\Lambda}(M_c,q, B_{\rm eff},\bar{\lambda})$$​ result. For a given $$M_c$$​, the left edge of the band is determined by the upper bound of $$\tilde{\Lambda}$$​ with the smallest allowed $$q$$​ value, while the right edge is determined jointly by the lower bound of $$\tilde{\Lambda}$$​ with $$q=1$$​ and the requirement that $$\bar{M}\lesssim\bar{M}_{\rm TOV}$$​.

Considering the solid black curve always goes leftward for a larger mass, there is then a critical mass value $$M_c\approx3.23 M_{\odot}$$​ beyond which the associated solid black curve no longer intersects with the blue band (GW170817) for any $$\bar{\lambda}$$​ value. Therefore, the possibility of IQSs is only allowed for $$M_{\rm TOV}\lesssim 3.23\, M_{\odot}$$​.

One can easily observe that for any $$\bar{\lambda}$$​ value, the lower bound of $$B_{\rm eff}$$​ is determined by the LIGO constraint $$\bar{\Lambda}\lesssim720$$​ of GW170817. The upper bound of $$B_{\rm eff}$$​ is set by the constraint of $$M_{\rm TOV}$$​ (the solid black curves) for small $$\bar{\lambda}$$​, and by the NICER constraint (the right edge of the cyan band) for large $$\bar{\lambda}$$​. More explicitly, when $$M_{\rm TOV}\gtrsim 2.59 M_{\odot}$$​ (the left of the left solid black curve) as suggested by GW190814, the joint constraints tell that any $$B_{\rm eff}$$​ is excluded for $$\bar{\lambda}^{1/4} \lesssim 1.17$$​, while $$92.7 \lesssim B_{\rm eff}/ (\rm MeV/fm^3)\lesssim 111$$​ for $$\bar{\lambda}^{1/4} \lesssim 1.37$$​, and $$111 \lesssim B_{\rm eff}/ (\rm MeV/fm^3)\lesssim 130$$​ for $$\bar{\lambda}^{1/4} \gtrsim 1.37$$​. For comparison, relaxing the maximum mass to $$M_{\rm TOV}\gtrsim 2.14\, M_{\odot}$$​ (the left of the middle solid black curve), we have the constrained parameter space $$49.5 \lesssim B_{\rm eff}/ (\rm MeV/fm^3) \lesssim 76.2$$​ for $$\bar{\lambda}^{1/4} \lesssim 0.67$$​, and $$76.2 \lesssim B_{\rm eff}/ (\rm MeV/fm^3)\lesssim 130$$​ for $$\bar{\lambda}^{1/4} \gtrsim 0.67$$​. We see that the constrained region can be well above the two-flavor stability lines within the uncertainties of $$a_4$$​, so that interacting $$ud$$​QM can be more stable than ordinary nuclei in the allowed parameter space.

# 5 summary

We have explored IQM, i.e., quark matter with strong interaction effects, including one-gluon exchange corrections and color superconductivity in a general parameterization, and managed to reduce the dependent parameter freedom into one single parameter $$\bar{\lambda}$$​ for the rescaled dimensionless EOS, or two parameters $$(B_{\rm eff},\bar{\lambda})$$​ for the physical dimensionful result. The parameter $$\bar{\lambda}$$​, as defined in Eq. ($$\ref{scaling95lam}$$), characterizes the relative size of strong interaction effects. A larger $$\bar{\lambda}$$​ results in a stiffer EOS. At large $$\bar{\lambda}$$​ limit, the EOS become $$\bar{p}=\bar{\rho}-1/2$$​, or equivalently $$p=\rho-2B_{\rm eff}.$$

We then studied the properties of IQSs (compact stars composed of IQM) utilizing the rescaled EOS. It turns out a larger $$\bar{\lambda}$$​ results in a larger (rescaled) maximum mass, with an upper bound given by the large $$\bar{\lambda}$$​ limit in which $$M_{\rm TOV}(\bar{\lambda}\to \infty)=35.35\, M_{\odot} (B_{\rm eff}/\rm MeV\, fm^{-3})^{-1/2}$$​. We also obtained a general expression of $$\bar{M}_{\rm TOV}(\bar{\lambda})$$​ in Eq. ($$\ref{Mbarmax}$$) with error only at $$0.1\%$$​ level, and computed the deformability of IQSs.

We have translated recent astrophysical observations into various constraints on the $$(B_{\rm eff}, \bar{\lambda}^{1/4})$$​ parameter space, assuming the related compact objects are IQSs. We obtained an upper bound for the maximum allowed mass for IQSs: $$M_{\rm TOV}\lesssim 3.23 M_{\odot}$$​. We found that the parameter space is confined into a window of moderate $$B_{\rm eff}$$​ and large $$\bar{\lambda}$$​ above the two-flavor stability line, enclosed by the constraints from the upper bound of $$\tilde{\Lambda}$$​ of GW170817, the recent NICER analysis, and the mass of the most massive compact star identified, as shown in Figure  3. Therefore, compact stars identified by the recent astrophysical observations can indeed be consistently interpreted as quark stars composed of interacting $$ud$$​QM within the determined parameter space. This study paves the way for future astrophysical observations to either confirm this possibility, further constrain it, or rule it out entirely.

Acknowledgments. We thank Bob Holdom and Jing Ren for helpful discussions. This research is supported in part by the Natural Sciences and Engineering Research Council of Canada.

# References

[1] A. R. Bodmer, Phys. Rev. D 4, 1601 (1971).

[2] E. Witten, Phys. Rev. D 30, 272 (1984).

[3] H. Terazawa, INS-Report-336 (INS, University of Tokyo, Tokyo) May, 1979.

[4] B. Holdom, J. Ren and C. Zhang, Phys. Rev. Lett. 120, no.22, 222001 (2018) [arXiv:1707.06610 [hep-ph]].

[5] C. Zhang, Phys. Rev. D 101, no.4, 043003 (2020) [arXiv:1908.10355 [astro-ph.HE]].

[6] Q. Wang, C. Shi and H. S. Zong, Phys. Rev. D 100, no.12, 123003 (2019) [arXiv:1908.06558 [hep-ph]].

[7] Q. Wang, T. Zhao and H. Zong, [arXiv:1908.01325 [hep-ph]].

[8] T. Zhao, W. Zheng, F. Wang, C. M. Li, Y. Yan, Y. F. Huang and H. S. Zong, Phys. Rev. D 100, no.4, 043018 (2019) [arXiv:1904.09744 [nucl-th]].

[9] K. Iida and T. Fujie, JPS Conf. Proc. 31, 011057 (2020).

[10] C. J. Xia, S. S. Xue, R. X. Xu and S. G. Zhou, Phys. Rev. D 101, no.10, 103031 (2020) [arXiv:2001.03531 [nucl-th]].

[11] J. Ren and C. Zhang, [arXiv:2006.09604 [hep-ph]].

[12] Z. Q. Miao, C. J. Xia, X. Y. Lai, T. Maruyama, R. X. Xu and a. P. Zhou, [arXiv:2008.06932 [nucl-th]].

[13] Z. Cao, L. W. Chen, P. C. Chu and Y. Zhou, [arXiv:2009.00942 [astro-ph.HE]].

[14] G. Aad et al.[ATLAS Collaboration], arXiv:1905.10130 [hep-ex].

[15] B. Acharya et al.[MoEDAL], [arXiv:2002.00861 [hep-ex]].

[16] L. W. Piotrowski, K. Małek, L. Mankiewicz, M. Sokołowski, G. Wrochna, A. Zadrożny and A. F. Żarnecki, Phys. Rev. Lett. 125, no.9, 091101 (2020) [arXiv:2008.01285 [astro-ph.HE]].

[17] E. Farhi and R. L. Jaffe, Phys. Rev. D 30, 2379 (1984).

[18] E. S. Fraga, R. D. Pisarski and J. Schaffner-Bielich, Phys. Rev. D 63, 121702 (2001) [arXiv:hep-ph/0101143 [hep-ph]].

[19] E. S. Fraga, A. Kurkela and A. Vuorinen, Astrophys. J. Lett. 781, no.2, L25 (2014) [arXiv:1311.5154 [nucl-th]].

[20] M. G. Alford, K. Rajagopal and F. Wilczek, Nucl. Phys. B 537, 443-458 (1999) [arXiv:hep-ph/9804403 [hep-ph]].

[21] K. Rajagopal and F. Wilczek, Phys. Rev. Lett. 86, 3492-3495 (2001) [arXiv:hep-ph/0012039 [hep-ph]].

[22] G. Lugones and J. E. Horvath, Phys. Rev. D 66, 074017 (2002) [arXiv:hep-ph/0211070 [hep-ph]]. Buballa:1998pr,PWang01,Ratti:2002sh,OsipovSM,HRZ2017.

[23] M. Alford and K. Rajagopal, JHEP 0206, 031 (2002) [hep-ph/0204001].

[24] M. Buballa and M. Oertel, Phys. Lett. B 457, 261 (1999) [hep-ph/9810529].

[25] P. Wang, V. E. Lyubovitskij, T. Gutsche and A. Faessler, Phys. Rev. C 67 (2003) 015210 [hep-ph/0205251].

[26] C. Ratti, Europhys. Lett. 61, 314 (2003) [hep-ph/0210295].

[27] J. Moreira, J. Morais, B. Hiller, A. A. Osipov and A. H. Blin, Phys. Rev. D 91, 116003 (2015) [arXiv:1409.0336 [hep-ph]].

[28] B. Abbott et al.[LIGO Scientific and Virgo], Phys. Rev. Lett. 119, no.16, 161101 (2017) [arXiv:1710.05832 [gr-qc]].

[29] B. P. Abbott et al.[LIGO Scientific and Virgo Collaborations], Phys. Rev. X 9, no. 1, 011001 (2019) [arXiv:1805.11579 [gr-qc]].

[30] B. Abbott et al.[LIGO Scientific and Virgo], Astrophys. J. Lett. 892, L3 (2020) [arXiv:2001.01761 [astro-ph.HE]].

[31] B. P. Abbott et al.[LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 121, no. 16, 161101 (2018) [arXiv:1805.11581 [gr-qc]].

[32] D. Radice, A. Perego, F. Zappa and S. Bernuzzi, Astrophys. J. 852, no. 2, L29 (2018) [arXiv:1711.03647 [astro-ph.HE]].

[33] A. Bauswein et al., AIP Conf. Proc. 2127, no. 1, 020013 (2019) [arXiv:1904.01306 [astro-ph.HE]].

[34] E. Annala, T. Gorda, A. Kurkela and A. Vuorinen, Phys. Rev. Lett. 120, no. 17, 172703 (2018) [arXiv:1711.02644 [astro-ph.HE]].

[35] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger and C. M. Biwer, Phys. Rev. Lett. 121, no. 9, 091102 (2018) Erratum: [Phys. Rev. Lett. 121, no. 25, 259902 (2018)][arXiv:1804.08583 [astro-ph.HE]].

[36] E. R. Most, L. R. Weih, L. Rezzolla and J. Schaffner-Bielich, Phys. Rev. Lett. 120, no. 26, 261103 (2018) [arXiv:1803.00549 [gr-qc]].

[37] E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, and L. Rezzolla, Phys. Rev. Lett. 122, no. 6, 061101 (2019) [arXiv:1807.03684 [astro-ph.HE]].

[38] L. R. Weih, E. R. Most and L. Rezzolla, Astrophys. J. 881, 73 (2019) [arXiv:1905.04900 [astro-ph.HE]].

[39] G. Montana, L. Tolos, M. Hanauske and L. Rezzolla, Phys. Rev. D 99, no. 10, 103009 (2019) [arXiv:1811.10929 [astro-ph.HE]].

[40] V. Dexheimer, L. T. T. Soethe, J. Roark, R. O. Gomes, S. O. Kepler and S. Schramm, Int. J. Mod. Phys. E 27, no. 11, 1830008 (2018) [arXiv:1901.03252 [astro-ph.HE]].

[41] A. Drago and G. Pagliara, Astrophys. J. 852, no. 2, L32 (2018) [arXiv:1710.02003 [astro-ph.HE]].

[42] E. P. Zhou, X. Zhou and A. Li, Phys. Rev. D 97, no. 8, 083015 (2018) [arXiv:1711.04312 [astro-ph.HE]].

[43] G. F. Burgio, A. Drago, G. Pagliara, H. J. Schulze and J. B. Wei, Astrophys. J. 860, no. 2, 139 (2018) [arXiv:1803.09696 [astro-ph.HE]].

[44] R. Abbott et al.[LIGO Scientific and Virgo], Astrophys. J. Lett. 896, no.2, L44 (2020) [arXiv:2006.12611 [astro-ph.HE]].

[45] B. Margalit and B. D. Metzger, Astrophys. J. Lett. 850, no.2, L19 (2017) [arXiv:1710.05938 [astro-ph.HE]].

[46] L. Rezzolla, E. R. Most and L. R. Weih, Astrophys. J. Lett. 852, no.2, L25 (2018) [arXiv:1711.00314 [astro-ph.HE]].

[47] M. Ruiz, S. L. Shapiro and A. Tsokaros, Phys. Rev. D 97, no.2, 021501 (2018) [arXiv:1711.00473 [astro-ph.HE]].

[48] M. Shibata, E. Zhou, K. Kiuchi and S. Fujibayashi, Phys. Rev. D 100, no.2, 023015 (2019) [arXiv:1905.03656 [astro-ph.HE]].

[49] P. Demorest, T. Pennucci, S. Ransom, M. Roberts and J. Hessels, Nature 467, 1081 (2010) [arXiv:1010.5788 [astro-ph.HE]].

[50] J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. Verbiest and D. G. Whelan, Science 340, 6131 (2013) [arXiv:1304.6875 [astro-ph.HE]].

[51] H. T. Cromartie et al., arXiv:1904.06759 [astro-ph.HE].

[52] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink and T. E. Strohmayer, Astrophys. J. Lett. 887, no.1, L21 (2019) [arXiv:1912.05702 [astro-ph.HE]].

[53] M. Miller, F. Lamb, A. Dittmann, S. Bogdanov, Z. Arzoumanian, K. Gendreau, S. Guillot, A. Harding, W. Ho, J. Lattimer, R. Ludlam, S. Mahmoodifar, S. Morsink, P. Ray, T. Strohmayer, K. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin and Y. Soong, Astrophys. J. Lett. 887, no.1, L24 (2019) [arXiv:1912.05705 [astro-ph.HE]].

[54] M. Alford, M. Braby, M. W. Paris and S. Reddy, Astrophys. J. 629, 969-978 (2005) [arXiv:nucl-th/0411016 [nucl-th]].

[55] S. Weissenborn, I. Sagert, G. Pagliara, M. Hempel and J. Schaffner-Bielich, Astrophys. J. Lett. 740, L14 (2011) [arXiv:1102.2869 [astro-ph.HE]].

[56] J. P. Pereira, C. V. Flores and G. Lugones, Astrophys. J. 860, no.1, 12 (2018) [arXiv:1706.09371 [gr-qc]].

[57] F. Weber, Prog. Part. Nucl. Phys. 54, 193-288 (2005) [arXiv:astro-ph/0407155 [astro-ph]].

[58] J. L. Zdunik, Astron. Astrophys. 359, 311 (2000) [astro-ph/0004375].

[59] P. Haensel, A. Y. Potekhin and D. G. Yakovlev, Astrophys. Space Sci. Libr. 326, pp.1 (2007).

[60] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).

[61] R. C. Tolman, Phys. Rev. 55, 364 (1939).

[62] A. E. H. Love, Proc. R. Soc. A 82, 73 (1909).

[63] T. Hinderer, Astrophys. J. 677, 1216 (2008) [arXiv:0711.2420 [astro-ph]].

[64] T. Hinderer, B. D. Lackey, R. N. Lang and J. S. Read, Phys. Rev. D 81, 123016 (2010) [arXiv:0911.3535 [astro-ph.HE]].

[65] S. Postnikov, M. Prakash and J. M. Lattimer, Phys. Rev. D 82, 024016 (2010) [arXiv:1004.5098 [astro-ph.SR]].

[66] T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009) [arXiv:0906.0096 [gr-qc]].

[67] J. Takátsy and P. Kovács, Phys. Rev. D 102, no.2, 028501 (2020) [arXiv:2007.01139 [astro-ph.HE]].

1. Another variant of two-flavor color superconductivity is the 2SCus phase, where $$u$$​ pairs with $$s$$​. Ref. [23] showed that the 2SCus phase has the same free energy as the "2SC+s" phase to order $$m^4_s$$​, so we neglect the discussion of it in this paper.↩︎

2. For 2SC, $$\mu=(\mu_u+2\mu_d)/3$$​. For "2SC+s" and CFL phase, $$\mu=(\mu_u+\mu_d+\mu_s)/3$$​.↩︎

3. We removed the electron contribution in this derivation. A numerical check approves this approximation.↩︎

4. As a necessary check of our general formula Eq. ($$\ref{eos95tot}$$), inserting Eq. ($$\ref{lam}$$) into it with the CFL factor where $$(\xi_4,\xi_{2a}, \xi_{2b})=(3,3/4,3)$$​ can reproduce the CFL result in Ref. [22], [56].↩︎