Radial oscillations in neutron stars from QCD

We study the stability against infinitesimal radial oscillations of neutron stars generated by a set of equations of state obtained from first-principle calculations in cold and dense QCD and constrained by observational data. We consider mild and large violations of the conformal bound, $c_{s} = 1/\sqrt{3}$, in stars that can possibly contain a quark matter core. Some neutron star families in the mass-radius diagram become dynamically unstable due to large oscillation amplitudes near the core.


I. INTRODUCTION
The equation of state (EoS) for neutron star matter, including the constraints from electric charge neutrality and beta equilibrium, is needed as the input to determine the structural properties of neutron stars (NS). It has to cover a wide range of densities, from below the nuclear saturation density, n 0 = 0.16 fm −3 , up to about 10n 0 [1], a real challenge. At such high densities, the question whether quark matter (QM) would be present in the core of NS arises naturally (see, e.g., Refs. [2,3]). Some possibilities for its formation mechanism being core-collapse supernova explosions [4] and the merger of relatively light NS [5]. Ultimately, this question must be settled by astrophysical observations. For instance, the NICER mission has recently shown the radii of NS with masses M ∼ 1.4M to be R ∼ 13km with a very good precision [6,7]. Still, one expects to find QM only in the interior of very massive NS, 2M M 2.2M [8][9][10][11][12], where central densities could be high enough.
In order to build the EoS for neutron star matter starting from first principles and perform systematic, controlled approximations, one can combine chiral effective field theory for nuclear matter at low baryon densities (n B n 0 ) [13,14] with perturbative QCD (pQCD) for deconfined QM at very high densities (n B n 0 ) [15,16]. The non-trivial sector of intermediate densities can be parametrized by a polytropic interpolation taking into account the two-solar mass constraint [17] and also the gravitational wave event GW170817 [18]. The maximum mass of NS can also be constrained to be in the range 2.1M M 2.3M [19,20]. The EoS for dilute nuclear matter, n B n 0 , was obtained using the ab initio approach of chiral effective field theory for light nucleons up to next-to-next-to-nextto-leading order in chiral perturbation theory (χPT) by Tews et al. [13]. To describe neutron star physics, β equilibrium and electric charge neutrality were implemented by Hebeler et al. [14]. The authors, then, extrapolate to densities above n B = 1.1n 0 by performing successive matchings onto polytropic equations of state. Following Hebeler et al., Kurkela et al. [17] added the perturbative QCD description at very high densities, being the polytropes used now as interpolations. In doing so, the extremes of low and high densities could be described from first principles, constraining the behavior of the neutron star equation of state. The pressure at high density had been computed to three loops in Ref. [15], including renormalization group effects for the strange quark mass and strong coupling α s up to O(α 2 s ). By adding the external constraint from the measurements of NS with masses around 2M [8,9], the relative uncertainty of the EoS band is reduced to less than 30% at all relevant densities. Limits on the tidal deformability, Λ, provided by the binary merger GW170817 [21] restrict even more the EoS [18].
Recently, Annala et al. [22] allowed large deviations of the adiabatic squared speed of sound, c 2 s = dP/d , from the conformal bound of 1/3 going up to the causality limit of 1, when performing the matching of EoSs through their new speed-of-sound interpolation method (see also Ref. [23]). Their uncertainty bands display a rapid softening from stiff nuclear matter to soft quark matter, i.e. their modified polytropic index 1 , γ, changes from γ nucl 2.5 to almost conformal QCD matter, γ pQCD ∼ 1, at crit ≈ 400-700 MeV/fm 3 . These are typical energy densities for the deconfinement transition at high temperatures found in ultrarelativistic heavy ion collisions [24].
By means of a careful analysis of the sudden bending of their polytropic index, Annala et al. [22] suggested that cold QM might exist in the core of the most massive NS, their presence being seizeable when c 2 s,max < 0.5 and independent of the nature of the transition, whereas one would obtain pure hadronic stars only if c 2 s,max > 0.7 and one has a first-order transition (or a rapid crossover). The authors argued that the presence of QM at the core of heavy NS should be considered the standard scenario and not an exotic alternative. In particular, they find that, if the conformal limit is not strongly broken, i.e. c 2 s < 0.4, a significant amount of QM might exist in massive NS.
In this paper we investigate the dynamical stability against infinitesimal radial pulsations of NS generated by a set of EoSs extracted from the constrained bands considered in Ref. [22]. We use the framework of Gondek et al. [25] to compute the fundamental mode frequencies (f n=0 ) for each EoS. We use these results to show novel features of stable NS with (without) QM cores having sub-conformal (extreme) values of c 2 s , in terms of the adiabatic index Γ. We adopt three representative (tabulated) EoS from Ref. [22], producing heavy NS satisfying multi-messenger bounds. They correspond to a subconformal EoS with a crossover transition leading to a seizeable QM core (Case I), an intermediate-c 2 s EoS with a crossover transition leading also to a seizeable QM core (Case II), and a high-c 2 s EoS with a strong first-order phase transition without a QM core (Case III).
We also study three other EoSs from constrained bands [26] satisfying c 2 max,s > 0.7 with first-order transitions (Cases 1, 3, 4) and one having a crossover transition with c 2 s < 0.5 (Case 2). The large c 2 s in the former set implies a decreasing number of degrees of freedom at increasing densities [27], something which is in sharp contrast with hot QCD where quarks and gluons are liberated at high temperatures [28]. A partial resolution of this issue was given in Ref. [22], where the authors suggested that NS with c 2 s > 0.4 could be dismissed if one includes further astrophysical constraints, e.g. the electromagnetic counterpart of GW170817 [29]. However, the unusual lumps at intermediate NS masses on the right-hand side of their corresponding mass-radius diagram [22] cannot be discarded using only observational constraints.
Here, we use Cases 3 and 4 (without loss of generality) to show that those prominent lumps are dynamically unstable, although not in the usual sense of hav-ing the (squared) fundamental-mode frequency negative, i.e. f 2 n=0 < 0. Instead, their respective infinitesimal Lagrangian radial displacement, ξ n=0 = (∆r/r) n=0 , produces large values (∼ 30), in contrast to the condition of being small, i.e. 0 < ξ 1, thereby allowing us to classify them, strictly speaking, as metastable. In fact, since radial oscillations with finite amplitudes induce non-linear effects that dynamically destabilize the star [30], one can exclude such families of NS from the interpolation scheme.
For completeness, we also compare our results to the ones obtained from a pair of standard pure nuclear matter EoSs: APR [31], calculated using many-body techniques with phenomenological potentials, and TM1 [32], obtained within relativistic mean field theory. The form of these EoSs is depicted in Fig. 3.
The paper is organized as follows. In Sec. II we summarize the main features of the first-order radial stability equations used in this paper. Section III shows our results for the behavior of the fundamental mode frequencies as a function of different stellar properties. Section IV presents our summary.

II. STABILITY OF THE FUNDAMENTAL RADIAL MODE
The relation between the adiabatic index Γ (via c 2 s ) and the existence of QM in the core of NS brings the question of their dynamical stability against linear (small) radial perturbations, where Γ plays a fundamental role. Imposing the condition ∂M/∂ c ≥ 0 plus astrophysical constraints helps reducing the band of EoS possibilities and, therefore, also the allowed region in the mass-radius diagram. However, to ensure full dynamical stability one needs to analyze radial pulsations through the behavior of the fundamental mode frequencies, f n=0 , which may be seen as a further constraint imposed by full general relativity.
Radial pulsation in NS with QM cores were studied using simple models for the deconfined phase, e.g. the MIT bag model, and sometimes assuming a (dis)continuous transition. For instance, Gupta et al. [33] considered a Glendenning construction [1] for the mixed phase, i.e. a continuous transition (see Ref. [34] for similar studies). This allows them to use the Sturm-Liouville approach originally derived for non-hybrid NS by Chandrashekar [35]. They found a reduction in the maximum mass (∆M max ∼ 0.3M ) and a new "kink" in the ω n vs c plane at the start of the mixed phase. On the other hand, first-order transitions might also occur in the interior of NS, leading to a more complicated analysis that requires the solution of supplementary equations [36,37]. The constrained neutron star matter EoSs of Annala et al. [22], on the other hand, exhibit first-order transitions that behave numerically as continuous rapid crossovers at the borders of the approximately flat mixed phase region. So, we can adopt the formalism of Gondek et al. [25], which is mathematically equivalent to the one developed by Chandrashekar and has the advantage of being more amenable to numerical calculations: one can easily impose the boundary conditions at the surface of the star and there is no need for derivatives of the adiabatic index Γ (see Appendix).
The formalism consists of a pair of coupled first-order differential equations for the relative radial displacement ∆r/r ≡ ξ and Lagrangian perturbed pressure ∆P . They can be written in matrix form as [25]: where ω is the oscillation frequency and Z, Q, R, S are complicated functions of the radial profiles for the pres- sure P (r), energy densitie (r), Γ(r), metric functions ν(r) and λ(r), after solving the Tolman-Oppenheimer-Volkov (TOV) equations for each EoS (for more details, see Appendix).
Dynamically stable NS are usually only characterized by the behavior of the lowest eigenvalue (fundamental mode), ω 2 n=0 = (2πf n=0 ) 2 , f n=0 being the (zero-mode) linear frequency which we use for numerical convenience. If f 2 n=0 > 0, then all f 2 n > 0 and the star is stable. If f 2 n=0 < 0, then there is at least one unstable mode and the star becomes unstable. The transition between these two stellar states occurs when f 0 → 0 [38].

III. RESULTS
By solving the oscillation equations for the representative tabulated EoSs discussed in Section I, one can obtain the frequencies for the oscillation modes n = 0, 1 . . . ., f 0 being the most relevant for the stability analysis.
Here, we compute the frequency dependence on the total gravitational mass, M , and stellar central energy density, c . Besides the cases presented in Ref. [22], we also consider, for comparison, the three representative bitropic EoSs (Cases I, II, III) of Ref. [17]. As shown in Fig. 1, they surpass the bound c 2 s = 1/3 (without a 1storder transition) between monotropes rapidly in Cases I and III and slowly in Case II. As expected, all these features are manifested in their stellar structure as well as in the planes f 0 vs M and f 0 vs c , as is clear in Fig. 2. In this figure, and analogous figures in what follows, one can see that static, ∂M/∂ c = 0, and dynamical, f 0 → 0, stability criteria coincide at M max .
From Fig. 2 one can see that Case III produces larger values of f 0 (∼ 4.5 kHz) compared to the other cases due to the fast change from being less deformable (more compact) to more deformable (less compact), which is related to a rapid stiffening/softening of c 2 s,max = 0.95 when changing polytropes at intermediate densities (see Fig.  1). This non-standard behavior might be considered unphysical since there is no actual transition that supports this abrupt behavior of the EoS 2 . It also goes up with c much faster than the standard pure nuclear-matter EoSs even at low densities. On the other hand, Case I matches the behavior of the APR even when both EoSs take very distinct values of c 2 s at high densities, i.e. being always causal and becoming rapidly acausal, respectively. So, it is hard to distinguish Case I from a hadronic star, at least in the f 0 vs c plane. Case II stands between these extreme behaviors, reaching c 2 s,max = 0.76 slowly. Notice that the jumps in c 2 s do not affect notoriously the associated frequencies.
A final point we can make from these EoSs is that the behavior of f 0 as a function of M suggests a qualitative distinction between pure and hybrid NS [39]: nearly flat f 0 (M ) characterizes nuclear matter cores (such as in the case of APR and TM1 EoSs), whereas QM cores produce very narrow flat regions. For example, Fig. 2 shows that Cases I and II behave as hadronic NS and Case III as a hybrid. Nevertheless, as we have already discussed, the last case cannot be considered a realistic situation.
We now explore the sets of selected EoSs from Ref. [22]. Recall that, besides the addition of the bounds yielded by the GW170817 event, the authors considered the dependence of the NS matter EoS on c 2 s (reaching values up to the causality limit) through their speed-ofsound interpolation method. They were motivated by the observation that NS with very high masses and small radii require the EoS to suffer a rapid stiffening somewhere in the core, perhaps due to a hadron-quark transition. In fact, Ref. [40] concluded, after a detailed analysis, that c 2 s should behave non-monotonically at intermediate densities, surpassing the conformal bound and then coming down to approach it again only at very high densities.
In Fig. 3 we show the non-monotonic behavior of the speed of sound as a function of the central energy den- sity. In the first set, Cases I and II display one bump exhibiting a mild breaking of the conformal bound, so containing QM cores 3 . Case III, on the other hand, has a second bump at very high densities with much higher c 2 s,max , so yielding purely hadronic stars. Besides, only Case III displays a somewhat notorious kink associated to a fast decrease of c 2 s before reaching max c through a discontinuous transition (see Fig. 3), after which it becomes unstable and the second bump in c 2 s begins. This hydrostatic and dynamical equilibrium structure is shown in Fig. 4.
As in the case of Kurkela et al. [17], changes in compacteness (visible in the M-R diagram) produce a hierarchy of maximal values for the zero-mode frequencies at intermediate densities, from the highest (Case III) to the lowest (Case II) in the f 0 vs c plane. Thus, NS in Cases I and II, having sizeable QM cores, exhibit lower values of max(f 0 ). This can be understood by noticing that a large amount of QM dissipates oscillations more rapidly due to its high-density nature 4 .
The second set of equations of state (Cases 1 to 4) behaves very differently with respect to c 2 s : they exhibit two bumps at relevant NS densities, as shown in Fig. 3. This is contrast with the analysis of Ref. [27] in which the authors suggest that the speed of sound could only behave non-monotonically once (see also Ref. [37] on the effects from phase transitions). It is also remarkable, as shown in Fig. 5, that this set of EoSs (especially Cases 3 and 4) shows a peculiar feature in the mass-radius diagram, which is also manifest in the M vs c plot. The effect resembles that of a first-order phase transition, even though the densities are too low for a transition to QM 5 . This feature is not in line with X-ray observations [42] and brings the question whether the solutions obtained from Cases 1 to 4 would be dynamically stable.
Using the criteria of Ref. [39] in Fig. 5, we see that Case 1 behaves as a nucleonic NS, whereas Case 2 as a hybrid NS. This agrees with the criteria of Ref. [22], where large c 2 s favors nuclear matter cores and mild violations of the conformal bound ensure seizeable QM cores. Cases 3 and 4 deserve a more detailed discussion.
The simple approach of verifying the stability by looking at dM/dR = 0 [43] is not adequate to study the non-monotonic Cases 3 and 4. So, we use the full solutions of Gondek's equations, in particular ξ n and f n . In this sense, dynamical stability will play the role of a further consistency constraint on the mass-radius diagram. Again, in Fig. 5 we see that these cases display, not surprisingly, more pronounced oscillations not found in previous works on hydrid stars, e.g. Refs. [33,37]. This lead us to ask whether the associated radial amplitudes, ξ = ∆r/r, behave consistently.
A careful analysis of the behavior of ξ n (with n = 0, 1, 2) in Cases 3 and 4 indicates that they undergo large (continuous/discontinuous) increments along the interior of NS, i.e. in the region from the crust's outer layers down to near the core (r R/2), reaching values around 20 − 40. In other words, one has radial oscillations with amplitudes much larger than the size of the star. It has been shown that large amplitudes near the core rapidly destabilize the NS due to non-linear general-relativistic effects via, for instance, the breaking of weak equilibrium [30]. These families of stars will most likely undergo a large expansion before collapsing to a black hole. Thus, the corresponding region in the EoS band should be discarded for being dynamically unstable.
The case of pathological amplitudes near the surface, i.e. within the outer layers of NS, which occurs in Cases 1 and 2, has already been reported in Ref. [44] and is known to be harmless due to its dilute nature. Our calculations show that the same happens for the EoSs of Kurkela et al. [17] and the first set of Annala et al. [22] with amplitudes around ξ n 30.

IV. SUMMARY AND CONCLUSIONS
Using the first-order formalism of Gondek et al. [25], we have investigated the dynamical stability of hadronic and hybrid NS obtained from a set of constrained equations of state derived from first principles [17,22]. Using the method of Ref. [39], we showed that it is possible to distinguish between hadronic and hybrid stars from their behavior in f 0 as a function of M in the set from Ref. [17] and the first set from Ref. [22] (Cases I, II, III).
The second set of EoSs from Ref. [22], especially Cases 3 and 4, behave very differently, though. The NS produced are quasi-stable according to our linear analysis of small oscillations, in the sense that their zero-mode frequency remains positive and real. However, the amplitude of their ξ undergoes a resonance-like behavior in the region near the transition point at the core of NS, reaching values around ξ ∼ 30, many times larger than the size of the star. Using the results from the full (non-linear) general-relativistic study performed in Ref. [30], we conclude that these large amplitudes evolve dynamically and make these NS strictly unstable. So, the corresponding region in the EoS band should be discarded.
The linear stability analysis procedure could also be incorporated systematically in interpolations of the equation of state for NS matter besides requiring that ∂M/∂ > 0. In future astronomical observations, it might be possible to compare tidal deformabilities of NS with the zero-mode frequencies by mapping their dependence for different NS gravitational masses [45,46].

ACKNOWLEDGMENTS
The authors thank Tyler Gorda for kindly sharing his constrained equations of state for cold QCD matter and for his useful comments and suggestions. This Apart from giving us the P (r), (r), and Γ(r) profiles, solving the TOV equations gives us the metric profiles for ν(r) through dν dr = − 2 P + dP dr , ν(R) = log 1 − 2M R and for λ(r) (plus the mass profile M(r) as input) λ(r) = − log 1 − 2M(r) r , λ(R) = −ν(R), both entering Gondek's equations, where we have added their associated boundary conditions [1]. The appropriate boundary conditions for Gondek's variables are deduced as follows. For ∆P it is required that (∆P ) surface = 0 since the pressure also goes to zero. At the center it is found that ξ center = (∆r/r) center = 1 by construction since ∆r → 0 as r → 0. We note that for self-consistency of the equations, the solutions for ξ should also be small in order to obtain displace-ments ∆r smaller compared to the star's size R. Furthermore, the requirement of physical smoothness at the NS center and finiteness everywhere in the star gives us (∆P ) center = −3(ξP Γ) center . By following this procedure, we have verified that our code reproduces satisfactorily the mode frequencies of Kokkotas and Ruoff [47] who use a different numerical procedure.