Cold quark matter with heavy quarks and the stability of charm stars

We study the effects of heavy quarks on the equation of state for cold and dense quark matter obtained from perturbative QCD, yielding observables parametrized only by the renormalization scale. We investigate the behavior of charm quark matter under the constraints of $\beta$-equilibrium and electric charge neutrality in a region of densities where perturbative QCD is in principle much more reliable. We also analyze the stability of charm stars, and find that such quark stars are unstable under radial oscillations.


I. INTRODUCTION
Heavy quark matter, i.e. quark matter including heavy flavors, could play a relevant role in extreme situations in the primordial quark-hadron transition [1][2][3]. Furthermore, cold quark matter brings about the possibility of charm stars. Those were investigated in the past within the crudest version of the MIT bag model, being ruled out due to instabilities under radial pulsations [4][5][6] (see also Refs. [7,8]). Experimentally, it is expected that FAIR with its Compressed Baryonic Matter (CBM) experiment will be able to produce charm quarks immediately after heavy-ion collisions with energies close or above the charm threshold [9].
In this work we investigate cold quark matter with heavy quarks using in-medium perturbative quantum chromodynamics (pQCD) [10]. We build the framework for the case with two massive flavors, strange and charm quarks, and determine the equation of state (EoS), which is the main ingredient for astrophysical applications. In practice, we extend the formalism developed for N f = N l + 1 flavors in Ref. [11], N l being the number of massless quarks and "1" the massive flavor, to the case with any number of massive (heavy) flavors, and explore their effects on the EoS of quark matter 1 . Given the equation of state, we consider β−equilibrated and electrically neutral charm quark matter and the possibility of charm stars. Using the method of first-order coupled oscillation equations of Gondek et al. [20], we find that such quark stars are unstable under radial oscillations.
This work is organized as follows. In Sec. II we summarize the main aspects of the perturbative QCD formalism for N f = N l + 1 flavors and presents our systematic extension to include heavy quarks in the framework. In Sec. III we build the EoS for charm quark matter. We solve the structure equations for charm stars, and study their 1 Quark mass effects on the equation of state for cold quark matter were first considered several decades ago [12,13], but only after two decades of being mostly ignored were computed within the modern MS renormalization scheme [11,[14][15][16][17][18][19].
stability under radial acoustic perturbations in Section IV. Section V presents our summary and outlook.

II. FRAMEWORK
A. Cold N f = N l + 1 quark matter The equation of state for quark matter at high densities and zero temperature was first obtained in perturbative QCD by Freedman and McLerran [12,21]; and Baluni [22] in a modified momentum subtraction scheme over four decades ago (cf. also Ref. [13,23]). Later, it was computed in the modern MS renormalization scheme for massless quarks in Refs. [24][25][26][27][28]. These results were then extended to include the role of a massive quark at two loops by Fraga and Romatschke [14] and three loops by Kurkela et al. [11].
The latter framework was designed to deal with N f = N l + 1 quark flavors, i.e. N l massless quarks plus 1 massive quark. Originally, the massive flavor was chosen to be the strange quark in order to study its influence on the stellar structure of quark stars. The perturbative QCD thermodynamic potential up next-to-next-toleading-order (NNLO) in the strong coupling α s , including the massless contribution plus a massive term, together with the mixed vacuum-matter (VM) diagrams and the corresponding ring terms, can be written as whereμ corresponds to the massive quark chemical potential, µ ≡ (µ 1 , ..., µ N l ) represents the vector chemical potential for the massless quarks, and m is the physical (renormalized up to the same order in the strong coupling) mass 2 associated to the massive quark flavor.
Although not explicit in Eq. 1, the free energy depends on a renormalization scale parameterΛ, which is an additional scale generated by the perturbative expansion. 2 One must be careful when defining these quark masses since they only make sense in the UV regime where asymptotic freedom takes place.
Then, after solving the renormalization group equations in the MS scheme for the strong coupling at this order, one obtains [14] α s (Λ) = 4π where β 0 = 11 − 2N f /3, β 1 = 51 − 19N f /3, L = 2 ln Λ /Λ MS . Since α s depends on N f , fixing the massive quark at some energy scale also depends on the number of flavors 3 . One defines the renormalization scale parameterΛ in terms of the natural scale at ultra-high densities, where quarks are massless, as beingΛ = 2µ s for any N f , and considers an uncertainty band defined by variations by a factor of 2. It is convenient to write this parameter in the formΛ = X i µ i /N f , where the sum runs over all quark flavors that are present in the system, and the dimensionless parameter X sits between 1 and 4 [11].
It is expected that at high densities not only the light quarks will present in a system of quark matter, but also some heavy flavors. In QCD, heavy quarks are meant to be the ones satisfying m Λ MS , i.e. quarks with masses that are very large compared the QCD natural scale. Usually, their influence is neglected in most calculations by invoking the heavy-quark decoupling theorem [30][31][32], which states that these quarks do not affect sizeably the observables calculated for light quarks, since their masses are of the order of the QCD energy scale.
We start by writing conveniently the total number of flavors in the form where N m is the number of massive quarks present in the system, and respecting the constraint N l + N m = N f . So, we add at least one massless quark for each massive flavor included. For example, for charm quark matter it will be convenient to write this sum over flavors as The usefulness of this way of writing will become clear when summing the massless and massive contributions to the total free energy. Of course, this represents only a convenient way of writing the degrees of freedom at the level of the formalism. Additional physical conditions are needed in order to control when a heavy partner appears actively. Such conditions can be introduced by choosing appropriate values of the renormalization scaleΛ, depending on the chosen heavy flavor to be introduced in the system. 3 For the massive quark running m(Λ), see Ref. [29].
With these points in mind, we write the thermodynamic potential (up to NNLO) for N l massless and N m massive quarks as where one must choose the number of massless flavors first when adding a massive one, so that is the massless contribution and the mixed massive contribution, where µ (i) = (µ 1 , ..., µ i ) is the massless vector chemical potential,μ (i) the massive (heavy) quark chemical potentials, and m (i) their corresponding masses. Here, Ω[...] indicates just the implicit parameter dependence (e.g. on N f ), whereas Ω(...) represents an explicit function dependence.
In the following, we apply these results to the case of charm quark matter and charm stars. We will also show in the next section that including heavy quarks makes the QCD thermodynamic potential less sensitive to the renormalization scaleΛ. In practice, its range of values is reduced in order to obtain a consistent thermodynamic transition between quark flavors, similar to results obtained in hot QCD [33,34].

III. CHARM MATTER
In this section we consider the simplest case of heavy quark matter, charm quark matter, which is composed of light quark matter plus charm quarks. Of course, it can only be realized above a given critical charm chemical potential. As we go to higher values of quark mass, asymptotic freedom makes the perturbative QCD formalism more reliable. To build the EoS for charm matter we first need to establish the constraints this system must respect below and above the charm threshold, in order to generate matter configurations stable under electroweak interactions.
A. Below the charm threshold: The condition of electric charge neutrality for a system with N f = 2 + 1 quarks (plus electrons) is given by where n i (µ i ) are the associated particle number densities for quarks and electrons in the system. The electron number density is approximated, as usual, by that of a free Fermi gas, i.e. n e = µ 3 e /(3π 2 ).
Weak reactions among light quark flavors are given by and yield the following relations between chemical potentials: We neglect the neutrino chemical potential since their mean free path are large compared to the size of a typical compact star. Solving simultaneously Eqs. (7) and (10) one is able to write all the quark and electron chemical potentials in terms of only the strange chemical potential, µ s .
When µ s crosses the charm quark threshold, the following weak equilibrium reaction is allowed to take place: yielding the condition 4 The electric charge neutrality condition turns into where we have included free muons, with n µ = (µ 2 µ − m 2 µ ) 3/2 /(3π 2 ), which appear when µ µ > m µ = 105.7 MeV. Lepton number conservation gives us µ µ = µ e .
Again, by solving simultaneously Eqs. (10), (12) and (13), we can express the quark and lepton chemical potentials only in terms of µ s . In the notation of Sec. II, the charm matter free energy corresponds to the case Now we need to fix the parameters entering the thermodynamic potential, i.e. running quark masses and strong coupling at some specific energy scale. Solving the renormalization group equations for the quark mass parameters up to second order in the strong coupling α s , one obtains the following results for the strange and charm quarks [29] m s (Λ) =m s α s π 4 This is in contrast to the high temperature case of heavy ion collisions, where only thermal equilibrium of charm quarks is reached with the surrounding medium [35].
with {m q } being the renormalization group invariant quark masses, i.e.Λ independent 5 . Since Eq. (2)  We assume that charm quarks are allowed in the system when where m medium c is the (unknown) in-medium charm mass 6 . Then, the renormalization scale parameter below and above the charm threshold are given by 7 where the approximations in the inequalities of Eqs.
(17a) and (17b) represent that fact that just before the threshold point the electron chemical potential takes its lowest value compared to the strange one, thus allowing us to make the approximation µ c ≈ µ s . The only way to go from Eq. (17a) to Eq. (17b) continuously through the N f transition is by requiring the factors X * and X to have the same range of possible values. To have them greater than 1, one needs values greater than 4/3 = 1.3333... for both, which implies a reduction in the renormalization scale band of the EoS when heavy quarks are included, something already found in thermal perturbative QCD with charm quarks [33]. 5 Expressing the quark masses in this way, one can see that their invariant masses can be fixed at independent energy scales, which is not obvious when using the quark mass function defined in Ref. [11]. 6 An exact value for the in-medium charm mass at finite density is still not known, whereas its vacuum mass at some fixed energy scale, m 0 c , can be extracted from lattice calculations.  (2) ] , so that the flavors are counted as N f = (1 + 1) (1) + (1 + 1) (2) = (u + c) (1) + (d + s) (2) . From this, one can identify the pressure as P = −Ω[N f = 2 + 1 + 1], and compute quark number densities using the standard thermodynamic relation n f = (dP/dµ f ) 8 . We define the total quark number density for charm matter, for a given X, as the baryon density as n B = n q /3, the total particle density as n = n q + n L , where n L = n e + n µ . In Fig. 1 we show the behavior of the relative particle populations for our β-equilibrated and electrically neutral charm quark matter system in the case of X = 4. Only above the charm threshold, charm quarks begin to contribute to the total number density, n. The location of the threshold depends on the value we choose for X and is within µ s ≈ 1.2−1.4 GeV for the band we consider.
To build the total pressure, one should be careful with the fact that the derivatives of the thermodynamic potential give rise to terms which have the form of ∂α s (Λ)/∂µ f sinceΛ ∝ µ f , as can be seen in Eqs. (17a)-(17b) (see also Ref. [11]). To keep thermodynamic consistency, one can take the quark and lepton number densities as the fundamental ingredients and build the other thermodynamic observables (e.g., pressure and energy density) imposing consistency on the number densities.
Thus, we define the total pressure of the system as where we have separated the contributions coming from N f = (u + c) (1) + (d + s) (2) , defining each term as and the lepton contribution as including strange quarks even at zero pressure, from which we start the integration of the particle densities, and adding the charm and leptons when crossing their respective thresholds.
We define the energy density as where the quark and lepton contributions are (2) Following this recipe, we can build the EoS, P = P ( ), by combining Eqs. (19) and (23) for a given X.
In Fig. 2 we plot the total pressure for charm matter, normalized by a Stefan-Boltzmann gas of quarks with N f = 4. From this one can see the usual behavior of the pressure for N f = 2+1 at intermediate densities, followed by a kink representing the charm threshold which softens the total (normalized) pressure. The charm quark contribution reduces the renormalization-scale uncertainty band for X at high densities, which also affects the behavior of the EoS a lower densities. An additional kink appears due to the muons. So, the charm EoS is largely softened, generating an apparent instability which could have astrophysical effects. In particular, it suggests the possibility of another kind of ultradense compact star: charm stars.

IV. CAN CHARM STARS EXIST IN NATURE?
The first quantitative study of the possibility of the existence of charm stars, i.e. strange stars satisfying the Bodmer-Witten hypothesis and having finite charm quark fractions at their cores, was carried out more than two decades ago Ref. [4] (see also Ref. [7]). The star bulk was described using the simplest version of the MIT bag model. After performing the stability analysis, the conclusion was that charm stars would be unstable.
We revisit this question by using our first-principle perturbative QCD for EoS for charm quark matter. We choose the parameter space to be in the range X ≥ 3, which satisfies the Bodmer-Witten hypothesis, as shown in Ref. [11], and perform the stability analysis as follows. The Tolman-Oppenheimer-Volkov (TOV) equations ensure the relativistic hydrostatic equilibrium of stellar configurations [5]. However, these configurations must also satisfy the thermodynamic condition ∂M/∂ c ≥ 0 [5]. The maximum mass configuration for a given stellar family is identified with the point where ∂M/∂ c = 0. In Fig. 3 we show our results for the mass as a function of the central energy density. It can be seen that the necessary condition for thermodynamic stability is satisfied in the two branches, one at relative low and another at much higher energy densities. However, we note that for the case X = 2 this condition is not satisfied when charm quarks appear, which is indicated by the black dots in Fig. 3. This is somewhat expected since it is difficult to have heavy quarks present in low-mass strange stars. In Table I we show the values of these observables at the charm threshold. On the other hand, for X > 3 the thermodynamic condition is satisfied when charm quarks are present, which would correspond to charm stars. In Fig.  4 we show the mass-radius diagram for quark stars made of N f = 2 + 1 + 1 quarks plus electrons and muons for different values of X.   The previous analysis provides a necessary but insufficient condition for stability of star configurations. One must still test for dynamical stability under radial pulsations. For that we use the method of Gondek et al. [20], solving a pair of first-order differential equations, one for the relative-displacement variable, ξ ≡ ∆r/r, and another for the Lagrangian perturbation, ∆P , with appropriate boundary conditions (see Appendix A). In the first-order radial pulsation formalism, amplitudes oscillate harmonically when the frequencies are such that Re(ω n ) > 0 and Im(ω n ) = 0, or increase exponentially if Re(ω n ) = 0 and Im(ω n ) > 0 9 . Since the radial oscillation equations represent a Sturm-Liouville problem [20], their eigenvalues (the eigenfrequencies) satisfy the ordering ω 2 0 < ω 2 1 < ω 2 2 < · · · < ω 2 n . So, if Im(ω 0 ) > 0 (and Re(ω 0 ) = 0) from some value of central energy density c , then all the higher modes will become complex too, representing the onset of the instability. For convenience we express our results in the terms of the linear frequency defined by f n ≡ ω n /(2π), So, if Im(f 0 ) > 0 continues finite and increasing for higher densities, all the configurations becomes unstable. In Fig. 5 we show that for densities above the maximum-mass strange star configuration (for X = 3), the stellar configurations increase their amplitudes even in the region where charm stars are expected, thus making them dynamically unstable. Since the same behavior was obtained for larger values of X, one can conclude from a perturbative QCD analysis that charm stars are unstable. One could ask if higher-order perturbative terms could in some way stabilize charm stars. However, a recent N 3 LO weak coupling expansion, which also includes nonperturbative terms, yielded minor modifications to the EoS [19].

V. SUMMARY AND OUTLOOK
In this paper we have extended the perturbative QCD N f = N l + 1 formalism in order to allow for the inclusion of heavy quark flavors in the EoS for cold and dense quark matter. In particular, we have investigated the effects of charm quarks in the equation of state in the case of βequilibrium and electric charge neutrality, and the possibility of charm stars, spanning a range in quark chemical potentials where pQCD is in principle much more reliable. After performing a radial stability analysis, it was concluded these stars would be unstable.
Although charm stars apparently do not exist, it is possible to have small amounts of charm quark matter in the core of the heaviest observed neutron stars (or, rather, hybrid stars), where a matching between a nuclear and a quark phases could be possible via the Glendenning construction for first-order phase transitions. Recently, a related possibility was investigated under the consideration of strange quark matter contaminated by charm quark impurities (in the sense of condensed matter physics), producing a QCD Kondo effect [41,42]. Moreover, a non-negligible amount of charm quarks could contribute to the EoS at the early stage of neutron star mergers, when very high densities are reached [43,44].
Our extended framework is appropriate to study the heavy sector of the QCD phase diagram (see Ref. [45] for related studies) which could exhibit new features, although it was shown in Refs. [46,47] that heavy quarks affect negligibly the chiral and deconfinement transitions at finite temperature.
Physical smoothness at the star's center requires that the coefficient of the 1/r term vanishes for r → 0. So, we impose that (∆P ) center = −3(ξΓP ) center , and normalize the eigenfunctions at the origin to be ξ(0) = 1. Since P (r → R) → 0, the Lagrangian perturbation of the pressure at the surface vanishes, i.e. (∆P ) surface = 0. Our code to study stellar charm stars reproduces the pulsation frequencies for the EoSs listed in Ref. [48].