DDark matter admixed neutron stars
Ben Kain
Department of Physics, College of the Holy Cross, Worcester, Massachusetts 01610 USA
Neutron stars could contain a mixture of ordinary nuclear matter and dark matter, such that darkmatter could influence observable properties of the star, such as its mass and radius. We study thesedark matter admixed neutron stars for two choices of dark matter: a free Fermi gas and mirror darkmatter. In addition to solving the multi-fluid Tolmon–Oppenheimer–Volkoff equations for staticsolutions and presenting mass-radius diagrams, we focus on two computations that are lacking inthe literature. The first is a rigorous determination of stability over the whole of parameter space,which we do using two different methods. The first method is based on harmonic time-dependentperturbations to the static solutions and on solving for the radial oscillation frequency. The secondmethod, which is less well-known, conveniently makes use of unperturbed, static solutions only. Thesecond computation is of the radial oscillation frequency, for fundamental modes, over large swathsof parameter space.
I. INTRODUCTION
If sufficient amounts of dark matter mix with theordinary matter in a neutron star, then dark mat-ter could influence measurable properties of the star.This suggests the exciting possibility that neutronstars could act as laboratories for indirectly measur-ing dark matter properties.An important question is how dark matter couldbecome mixed with ordinary matter in a neutronstar. One well-studied possibility is through cap-ture [1–7]. If, in addition to gravitational interac-tions, dark matter has non-gravitational interactionswith the ordinary matter in neutron stars, the ex-treme densities of neutron stars make them excel-lent targets. Captured self-annihilating dark mat-ter could potentially heat the star [2–4], while non-self-annihilating dark matter would accumulate. Forbosonic dark matter, this accumulation could lead tothe formation of a small black hole which destroysthe parent neutron star [1, 3–6], while for fermionicdark matter, degeneracy pressure is able to stabi-lize the star unless the dark matter particle massexceeds 10 GeV [8]. Over the lifetime of a neutronstar, the mass of the accumulated dark matter doesnot typically exceed 10 − M (cid:12) [1, 2, 9, 10], whichhas a negligibly small effect on the mass of the star.Another possibility for mixing is that the accumu-lation of dark matter occurs during stellar formation.In [10], it is argued that a natural possibility is thatthe formation of a dark matter minihalo causes aneutron star to form from collapse. (For other pos-sibilities, see [9, 11].) Detailed studies of this processare lacking and would be interesting. In this paper,we focus on bulk properties of the star, such as itsmass and radius, and have this latter possibility formixing in mind.An appealing model for non-self-annihilating darkmatter is asymmetric dark matter [12, 13], in which a conserved charge keeps dark matter from self-annihilating and an imbalance in the early universebetween dark matter and anti-dark matter leads tothe abundance of dark matter observed today. Acommon description is a gas of Dirac fermions, pos-sibly with self-interactions (see, for example, [8]).Another description is as mirror dark matter [14],which is motivated by the assumption that the Uni-verse is parity symmetric. The addition of new par-ticles to restore parity to the Standard Model leadsto mirror baryons as viable candidates for dark mat-ter (see [15–22] and references therein).Dark matter admixed neutron stars are two-fluidsystems, in which the first fluid describes ordi-nary nuclear matter through an equation of statefor a neutron star without dark matter and thesecond fluid describes dark matter. Propertiessuch as the mass and radius of the star are foundby solving the multi-fluid Tolmon–Oppenheimer–Volkoff (TOV) equations [19, 23, 24]. Null resultsfrom dark matter direct detection experiments [25–27] have placed stringent constraints on the darkmatter-nucleon coupling strength. From the per-spective of the TOV equations, this is generallytaken to mean that the dark matter-nucleon cou-pling strength is negligibly small [8, 11] and thatdark matter admixed neutron stars are two-fluid sys-tems in which the only inter-fluid interactions aregravitational.Early work on dark matter admixed neutron starswas undertaken by Henriques, Liddle, and Moor-house [24, 28, 29] in their study of boson-fermionstars. In these papers, the fermions were a freeFermi gas of neutrons and the scalar bosons couldbe interpreted as dark matter, though this was notexplicitly stated. Reference [29] presented an un-derappreciated method for determining stability intwo-fluid systems, which we make use of in Sec. III.Ciarcelluti and Sandin [19, 20] used mirror baryonsas dark matter. Subsequent studies with mirror dark a r X i v : . [ g r- q c ] F e b atter by Goldman et al. [21, 22] allowed for a mir-ror baryon mass smaller than the ordinary baryonmass. A series of papers by Leung et al. [30–32]used a free Fermi gas as dark matter and studied sit-uations in which dark matter forms either the coreor the halo of the star. As far as we are aware,Refs. [30–32] are the only papers that have computedradial oscillation frequencies for such systems usingtwo-fluid methods. A large number of studies havesince followed [8–11, 33–42], studying such things asself-interacting dark matter [33, 34], ordinary mat-ter that includes hyperons [33, 41] or strange quarkmatter [37], and a computation of the tidal deforma-bility [9, 11, 42].In this paper, we study dark matter admixed neu-tron stars using two different models for asymmetricdark matter. For the first model, we use a free Fermigas. Although self-interactions have been consideredin a number of works and shown to lead to interest-ing effects (see, for example, [8, 10, 11, 33, 34]), forsimplicity we do not include them. For the secondmodel, we use mirror dark matter, which was one ofthe first considerations in the study of dark matteradmixed neutron stars [19].In addition to solving the TOV equations for staticsolutions, which gives the mass and radius of thestar, we make a careful study of the stability ofthese solutions. Rigorous determinations of stabil-ity with respect to small perturbations over largeswaths of parameter space is lacking in the litera-ture. We present two different methods for deter-mining stability. The first method is to perturb thestatic solutions with harmonic perturbations and tosolve for the squared radial oscillation frequency. Wedo this using an approach developed in [43], whichderived a system of pulsation equations for an ar-bitrary number of perfect fluids with only gravi-tational inter-fluid interactions and whose solutiongives the squared radial oscillation frequency. Thesecond method we use was developed in [29] and con-veniently makes use of only unperturbed, static solu-tions. Interestingly, we find regions of stable param-eter space for which a naive analysis of the single-fluid equations of state would not have deemed sta-ble.Using the pulsation equations of [43], we alsomake a systematic determination of the radial os-cillation frequencies for large swaths of stable pa-rameter space. This too appears to be lacking inthe literature. Although radial oscillation modes donot couple to gravitational waves, they are, in princi-ple, observable [44] and the hope is that their studyand detection can reveal details of the inner struc-ture of the star. We find interesting results here aswell, in that the oscillation frequencies of dark mat-ter admixed stars can be larger than the maximum possible frequencies of single-fluid stars made fromthe same equations of state.In the next section, we review the multi-fluid TOVequations and the equations of state that we will beusing. In Sec. III, we study stability, with some ofthe details given in the Appendix. In Sec. IV, wepresent mass-radius diagrams. In Sec. V, we com-pute radial oscillation frequencies. We conclude inSec. VI. II. EQUATIONS AND EQUATIONS OFSTATE
Dark matter admixed neutron stars are solutionsto the multi-fluid TOV equations [19, 23, 24]. TheTOV equations follow from the Einstein field equa-tions and the equations of motion in a sphericallysymmetric spacetime with static matter. In writingequations, we use units such that c = G = (cid:126) = 1.The spherically symmetric metric may then be writ-ten as ds = − e ν ( r ) dt + dr − m ( r ) /r + r d Ω , (1)where d Ω = dθ +sin θ dφ and the metric function m ( r ) gives the total mass inside a radius r . Forstatic solutions, the metric function ν ( r ) decouplesand is not needed. It is needed when determiningradial oscillation frequencies and is discussed in theAppendix.Null results from direct detection experiments fordark matter have placed stringent constraints onthe dark matter-nucleon coupling strength [25–27].From the perspective of the TOV equations, thisis generally taken to mean that any interaction be-tween dark matter and the ordinary matter of theneutron star is negligibly small [8, 11] and that darkmatter and ordinary matter can be modeled as sep-arate fluids with only gravitational inter-fluid inter-actions. The energy-momentum tensor, then, sep-arates, T µν = (cid:80) i T µνi , where the subscripted i in-dicates the fluid (either ordinary or dark matter),and each energy-momentum tensor takes the perfectfluid form, T µνi = diag( (cid:15) i , p i , p i , p i ) , (2)where (cid:15) i ( r ) and p i ( r ) are the fluid’s energy densityand pressure. Each energy-momentum tensor is alsoconserved, ∇ µ T µνi = 0, which gives the equations ofmotion. The TOV equations are then dm i dr = 4 πr (cid:15) i dp i dr = 4 πr p + mr (1 − m/r ) ( (cid:15) i + p i ) , (3)2here the first equation follows from the Einsteinfield equations, the second from the equations of mo-tion, m = (cid:80) i m i , and p = (cid:80) i p i . In addition to theabove equations, we shall need an equation for thenumber of particles inside a radius r , N i ( r ), whichis d N i dr = 4 πr n i (cid:112) − m/r , (4)where n i ( r ) is the number density.With only gravitational inter-fluid interactions,the equations of state also separate, (cid:15) i = (cid:15) i ( p i ),where the energy density only depends on the pres-sure of the same fluid. For ordinary matter, we usethe analytical fit [45] to the SLy equation of state[46]. SLy is a unified equation of state, obtainedfrom a single effective nuclear Hamiltonian, allow-ing for a smooth transition between core and crustsof a neutron star. The analytical fit further smoothsthe equation of state. This level of smoothness isunnecessary for one of the methods we use to de-termine stability in the next section, but is helpfulfor the other method, which is also used to computeradial oscillation frequencies, because it requires tak-ing derivatives of the equation of state.It is worth noting that many equations of state inthe literature for the ordinary matter of a neutronstar, including [45], list the baryonic number den-sity and not the number density for the fluid (i.e.the density of fluid elements). But, as we will see, itis the fluid’s number density that is needed for deter-mining stability. From the thermodynamic identityat zero temperature, (cid:15) + p − µn = 0, where µ = d(cid:15)/dn is the chemical potential for the fluid, one finds thatthe number density for the fluid can be computedfrom the energy density and pressure, n ∝ exp (cid:18)(cid:90) d(cid:15)(cid:15) + p (cid:19) , (5)where the proportionality constant does not affectthe determination of stability and therefore doeshave to be known.As mentioned in the Introduction, we consider twopossibilities for asymmetric dark matter. For thefirst possibility, we use a simple and common de-scription, modeling dark matter as a free Fermi gas.The well-known equation of state and number den-sity for a free Fermi gas are [47, 48] (cid:15) = 12 π (cid:90) k F dk k (cid:113) k + m f = 18 π (cid:34) k F (cid:113) k F + m f (2 k F + m f ) − m f ln k F + (cid:113) k F + m f m f (cid:35) p = 16 π (cid:90) k F dk k (cid:113) k + m f = 124 π (cid:34) k F (cid:113) k F + m f (2 k F − m f )+ 3 m f ln k F + (cid:113) k F + m f m f (cid:35) n = k F π , (6)where m f is the fermion mass and k F is the Fermimomentum. The Fermi momentum is eliminatedwhen forming (cid:15) = (cid:15) ( p ) and n = n ( p ), making thefermion mass the only free parameter.For the second possibility, we consider mirror darkmatter. For this case, following [19–22], we use the same equation state for dark matter that we use forordinary matter.It is straightforward to show that the inner bound-ary conditions for the TOV equations in (3) and theparticle number equation in (4) are m i ( r ) = O ( r ), p i ( r ) = p ci + O ( r ), where p ci is the central pressurefor fluid i , and N i = O ( r ). The central pressuresuniquely identify solutions. Upon specifying centralpressures, the TOV and particle number equationsmay be integrated outward from some small r . Atsome point during the integration, the pressure ofone of the fluids will hit zero, p i ( R i ) = 0, markingthe edge of fluid i at r = R i . At this point, the inte-gration is broken and restarted using the single -fluidequations and the equation of state of the remainingfluid. When the pressure of the remaining fluid hitszero, p j ( R j ) = 0, the edge of fluid j , as well as theedge of the star, is at r = R j . We then have for thetotal mass of the star, M = (cid:80) i m i ( R i ), and for thetotal number of particles for fluid i , N i = N i ( R i ).Solutions to the TOV equations, when using theequations of state presented in this section, are ex-amples of dark matter admixed neutron stars. If R dm < R om , the star has a dark matter core, whileif R dm > R om , it has a dark matter halo. In the fol-lowing, we display the parameter space of solutionsusing the central pressures ( p c om , p c dm ), since theyuniquely identify solutions. III. STABILITY: CRITICAL CURVES
Once a solution is found, an important question iswhether it is stable with respect to small perturba-3ions. In systems with only a single fluid, this ques-tion is straightforward to answer, since it is well-known that the transition from stable to unstableoccurs at the solution with the largest mass [47–49].Since each solution is uniquely identified by a sin-gle quantity (the central pressure of the fluid), thestatic solution with the largest mass constitutes asingle point in the parameter space of solutions andis called the critical point. The critical points forsingle-fluid stars constructed with the SLy and freeFermi gas equations of state from the previous sec-tion are ( p c SLy ) crit = 860 .
24 MeVfm ( p c Fermi ) crit = 291 . (cid:16) m f (cid:17) MeVfm . (7)The situation is more complicated in two-fluid sys-tems, such as dark matter admixed neutron stars.Solutions are identified by two quantities (the cen-tral pressure of each fluid). The transition from sta-ble to unstable is then marked by a critical curve in parameter space. Determining the critical curveis not as straightforward as determining the criticalpoint in the single-fluid case. We offer two methodsfor its determination.The first method is to perturb the static solutionswith time-dependent harmonic perturbations, whichdepend on the radial oscillation frequency. Thisleads to a system of pulsation equations, whose so-lutions give the squared radial oscillation frequency.If the squared radial oscillation frequency is posi-tive for the fundamental solution, the correspond-ing static solution is stable; otherwise it is unstable.Such a system of pulsation equations was derived in[43] for an arbitrary number of perfect fluids withonly gravitational inter-fluid interactions and is re-viewed in the Appendix (see also [50]).The second method was developed by Henriques,Liddle, and Moorhouse in their study of boson-fermion stars [29]. The details are reviewed in theAppendix. The conclusion of their analysis is thatthe critical curve is defined by dMd p = dN om d p = dN dm d p = 0 , (8)where M and N i are the total mass and fluid num-ber of a static solution and p is a vector in parame-ter space that is simultaneously tangent to the levelcurves of M and N i . It can be shown that if two ofthe quantities in (8) are zero, then the third is also[29, 51]. We stress that N om is the total number offluid elements and not the baryonic number. Thisis the reason why the number density for the fluidmust be known, which can be computed from theenergy density and pressure using Eq. (5). m f = 1 GeV p c om [MeV/fm ] p c d m = ( m f = G e V ) [ M e V / f m ] FIG. 1. The parameter space of static solutions is shownas a function of the central pressures p c om and p c dm , wheredark matter is taken to be a free Fermi gas with fermionmass m f = 1 GeV. Each point represents a static solu-tion to the multi-fluid TOV equations. The thick blackline is the critical curve, separating stable static solu-tions from unstable ones. Stable solutions are colored,with green indicating a static solution with a dark matterhalo and red indicating a dark more core. In the original paper [29], the critical curve wasfound by plotting contour lines for N om and N dm and determining those points where the contour linesmeet, but do not cross, so that their tangents areequal. Such points give the critical curve. An alter-native procedure [52], which is the one we use here,is to first compute contour lines of either N om or N dm in the two-fluid system. Moving along a sin-gle contour line, we determine the point where M isan extremum (in practice, we find that it is a maxi-mum). These points give the critical curve.A benefit of the first method is that it can domore than just determine stability, since it can com-pute the radial oscillation frequency for an arbitrarystatic solution. Its disadvantages are that it is timeconsuming to find a solution and it requires takingderivatives of the equation of state, which may bedifficult to do if the equation of state is insufficientlysmooth. The second method does not suffer from ei-ther of these disadvantages, but can only determinestability. We have confirmed numerically that bothmethods give the same answer, which gives confi-dence that the code we are using is working properly.The figures presented in this section were made us-ing the second method. As far as we are aware, thisis the first time that the second method has beenapplied to dark matter admixed neutron stars whena realistic equation of state is used for the ordinarymatter.In Fig. 1, we show results for the free Fermi gaswith fermion mass m f = 1 GeV. The thick black line4 a) m f = 0 : p c om [MeV/fm ] p c d m = ( m f = G e V ) [ M e V / f m ] (b) m f = 0 : p c om [MeV/fm ] (c) m f = 2 GeV p c om [MeV/fm ] (d) m f = 5 GeV p c om [MeV/fm ] FIG. 2. The same as Fig. 1, except with fermion masses, m f , as indicated above each plot. is the critical curve. The colored parameter space in-dicates stable static solutions, with green indicatingstatic solutions with a dark matter halo and red adark matter core. For sufficiently small p c om or p c dm ,the critical curve is seen to agree with the single-fluidcritical points in (7). This is expected, since if one of p c om or p c dm is small while the other is large, the fluidwith the large central pressure dominates and weeffectively have a single-fluid system. Interestingly,there is a region of stable parameter space, in theupper-right corner, where p c om and p c dm are greaterthan their single-fluid critical values ( p c om ) crit and( p c dm ) crit in (7). m f = 1 GeV can be taken to approximate thetransitional mass, where masses above and belowthis value lead to qualitatively different results. Thetransitional mass is expected to be somewhere nearthe baryon mass of 938 MeV. This is evident inFig. 2, where we show critical curves for fermionmasses above and below 1 GeV. First, for m f be-low 1 GeV, we see from Figs. 2(a) and 2(b) thatthe critical curve moves toward extending beyondthe single-fluid ( p c dm ) crit , but not beyond the single-fluid ( p c om ) crit . This flips for m f above 1 GeV,where we see from Figs. 2(c) and 2(d) that now thecritical curve moves toward extending beyond thesingle-fluid ( p c om ) crit , but not beyond the single-fluid( p c dm ) crit . We do not compute the precise mass wherethis transition occurs, but simply take m f = 1 GeVto approximate its value.In Fig. 3 we show the critical curve for mirror darkmatter. Since the same equation of state is usedfor both ordinary and dark matter, both the criticalcurve and the line separating a dark matter corefrom a dark matter halo are symmetric in parameterspace. We find again that there is a stable region ofparameter space where p c om and p c dm are greater thantheir single-fluid critical values, although the regionis smaller than with the free Fermi gas. p c om [MeV/fm ] p c d m [ M e V / f m ] FIG. 3. The same as Fig. 1, except for mirror dark mat-ter, in which dark matter has the same equation of stateas ordinary matter.
In Sec. V, when we look at radial oscillation fre-quencies, we will gain some insight as to why thecritical curves extend past their single-fluid criticalvalues and how this depends on fermion mass.
IV. MASS-RADIUS DIAGRAMS
Two of the most important properties to computeare the mass and radius of the star, since these prop-erties are observable. They are typically presentedin mass-radius (MR) diagrams by plotting the massas a function of radius for a given solution to theTOV equations. In systems with only a single fluid,the mass-radius parameter space is a curve.Things are not as simple in two-fluid systems, suchas dark matter admixed neutron stars, since the pa-rameter space is larger. In the literature, a singleMR curve or a sequence of MR curves is often pre-sented. Unfortunately, this can give the false im-5 f = 1 GeV R om [km] M [ M - ] FIG. 4. The mass as a function of the visible radius,which is the radius of ordinary matter, for the stablestatic solutions shown in Fig. 1 (dark matter is a freeFermi gas with fermion mass m f = 1 GeV). The thickblack curve is for the single-fluid star with only ordinarymatter, i.e. for a neutron star without dark matter. Thecolor scheme is the same as in Fig. 1, with green indicat-ing a dark matter halo and red a dark matter core. pression that there is an important property thatconnects two points on the same curve that does notconnect two points on different curves. While thisis sometimes the case, commonly the property thatconnects two points on the same curve is arbitrary.Another way to understand this is that any curvethat slices through the parameter space in Figs. 1–3 gives an MR curve, but there are infinitely manyways to slice through the parameter space.As an alternative, we present MR diagrams thatshow mass-radius relations for the entire stable pa-rameter space shown in Figs. 1–3. In Fig. 4, wedisplay the MR diagram for the free Fermi gas withfermion mass m f = 1 GeV. We plot the total massof the system versus the visible radius, which is theradius of ordinary matter. The thick black line givesthe MR curve for the single-fluid star with the SLyequation of state, i.e. for a neutron star without darkmatter. The color scheme is the same as in Figs. 1–3, with green indicating a dark matter halo and redindicating a dark matter core. For m f = 1 GeV,we can see that much of the MR parameter spaceis taken up by a dark matter core. Further, the in-clusion of dark matter does not allow the mass orradius to extend past the values for a neutron starwithout dark matter. In this sense, the inclusionof dark matter leads to a decrease in the mass andradius of the star [19].We mention again that all points plotted in Fig. 4are stable, as determined using the methods of theprevious section. Those points underneath the thickblack line correspond to static solutions in which or- dinary matter is dominating and we effectively havea single-fluid system. As expected in such a case,stable static solutions do not extend past the peakof the thick black line.We previously mentioned that m f = 1 GeV is atransitional mass, in that larger and smaller massesgive qualitatively different results. This was evidentin the previous section with stability and is also evi-dent with MR diagrams. In Fig. 5, we show MR di-agrams for fermion masses above and below 1 GeV.Beginning with Fig. 5(b), we see that as the fermionmass is lowered, the parameter space with a darkmatter halo expands and the total mass of the sys-tem is able to increase past the values for a neutronstar without dark matter. This continues in Fig. 5(a)with an even smaller fermion mass. We concludethat by decreasing the fermion mass, the total massof the system can increase [21] and dark matter gen-erally encompasses ordinary matter, forming a halo.What is happening is that as the fermion mass is de-creased, which requires decreasing the dark mattercentral pressure to retain stability, dark matter hasan increasingly weaker effect on ordinary matter. Inother words, as the fermion mass is decreased, or-dinary matter begins acting as if it is a single fluidgiven by the black curve. Not shown is that whendark matter forms a halo, the smaller fermion massallows dark matter to extend farther out. With darkmatter extending farther out, more dark matter par-ticles can stably exist, which can compensate for thesmaller fermion mass and increase the total mass ofthe system. On the other hand, when dark matterforms a core, there are not enough dark matter par-ticles to compensate for the smaller fermion mass,which is why the red region is squeezing close to theblack curve.Now consider when the fermion mass is increasedabove 1 GeV. We can see in Figs. 5(c) and (d) that asthe fermion mass is increased, the parameter spacebecomes increasingly dominated by a dark mattercore. Not shown is that the radius of the dark mat-ter core tends to decrease with increasing fermionmass. If ordinary matter has a relatively large cen-tral pressure, then the small dark matter core haslittle effect, as seen by the red region squeezing closeto the black curve in 5(d). On the other hand, if or-dinary matter has a relatively small central pressure(which, in the absence of dark matter, would be lo-cated on the tail end of the black curve with largeradii), the dark matter core pulls this ordinary mat-ter in to smaller radii, giving the low mass solutionsthat can be seen at the very bottom of Fig. 5(d). Asthe fermion mass is increased further, one continuesto find these low mass solutions with ever smallermasses and radii around 10 km. Eventually solu-tions with planet-like masses are found, which have6een dubbed “dark compact planets” in [36].In Fig. 6 we show the MR diagram for mirror darkmatter. Compared to the free Fermi gas with m f =1 GeV in Fig. 4, the region with a dark matter core issimilar. A notable difference is an expanded regionwith a dark matter halo, in which we find 1–2 M (cid:12) solutions, but with visible radii in the 2–8 km range. V. RADIAL OSCILLATIONS
Chandrasekhar initiated the study of stellar oscil-lations of neutron stars when he derived a pulsationequation whose solution gives the squared radial os-cillation frequency [53]. His pulsation equation wassubsequently rewritten in various ways [54–57], insome cases to facilitate numerical solutions. As men-tioned in Sec. III, the squared radial oscillation fre-quency can be used to determine the stability of astatic solution. Although radial oscillation modesdo not couple to gravitational waves, they are, inprinciple, observable by the emission of electromag-netic radiation from the surface of the star (see, forexample, [44]). The hope is that their study canreveal details of the inner structure of the star. Ra-dial oscillation frequencies have been computed for alarge number of equations of state (see, for example,[44, 56–61]).Chandrasekhar derived his pulsation equation fora single-fluid system. Recently, a system of pulsa-tion equations was derived for an arbitrary numberof perfect fluids with only gravitational inter-fluid in-teractions [43] (see also [50]). In [43], these equationswere used to compute radial oscillation frequenciesin one-, two-, and three-fluid systems, where all flu-ids were taken to be a free Fermi gas. In this sectionwe use these equations to study radial oscillationsof dark matter admixed neutron stars. This is thefirst time the pulsation equations of [43] have beenapplied to a system where one of the fluids is de-scribed with a realistic equation of state. The pul-sation equations and how they are solved is reviewedin the Appendix.Radial oscillations of dark matter admixed neu-tron stars is understudied. As far as we are aware,the only works that have computed such oscillationfrequencies using two-fluid methods are [30–32], us-ing the equations of [50], and only [31] presentedresults beyond those used in determining stability.(References [62, 63] computed radial oscillation fre-quencies of two-fluid systems, but did so using Chan-drasekhar’s single-fluid pulsation equation.) Fre-quencies have also been computed by Fourier trans-forming results from simulations using full numericalrelativity [52]. Our aim in this section is to make a systematic computation of oscillation frequencies forthe fundamental radial mode.We begin first by solving Chandrasekhar’s single -fluid pulsation equation for the radial oscillation fre-quency, ω , of the fundamental solution for single-fluid stars using the SLy and free Fermi gas equa-tions of state. The results are shown in Fig. 7 asa function of the central pressure, which uniquelyidentifies the static solution. The solid black line isfor SLy and the dashed blue line is for the free Fermigas. Note that the frequencies hit zero at the criti-cal central pressures given in (7). This is expected,since the critical central pressures mark the point atwhich the static solutions transition from stable tounstable, which occurs when the (squared) radial os-cillation frequency (of the fundamental mode) tran-sitions from positive to negative. Note also that thefrequencies of a free Fermi gas with fermion mass m f = 1 GeV are smaller than those for the neutronstar (with the SLy equation of state) over much ofthe parameter space. For m f = 0 . m f = 2GeV, the dashed blue curve grows by a factor of 4,making it larger than the neutron star over much ofthe parameter space. It is useful to keep these factsin mind in the following.Unfortunately, the computation of the radial os-cillation frequency in a two-fluid system is time con-suming. For this reason, we do not present resultsfor all cases, nor over the entirety of the parame-ter space, that we considered previously. In Fig. 8we show results for a free Fermi gas with fermionmasses m f = 0 .
5, 1, and 2 GeV. In each plot, thethick black line is the critical curve first shown inFigs. 1 and 2. We can see that near the criticalcurves, the oscillation frequencies head toward zero,as expected. This offers some visual evidence thatthe two methods for determining stability [29, 43]discussed in Sec. III agree, though we note that wehave confirmed that they agree to much higher pre-cision than that shown in Fig. 8.Consider now the bottom of the plot in Fig. 8(a)for m f = 0 . p c om . This tells us that only when we have largervalues of p c dm and smaller values of p c om is dark mat-ter able to affect substantially the frequency. Oncewe compare this to the other plots in Fig. 8, weconclude that this is because of the smaller fermionmass of m f = 0 . a) m f = 0 : R om [km] M [ M - ] -1 (b) m f = 0 : R om [km] (c) m f = 2 GeV R om [km] (d) m f = 5 GeV R om [km] FIG. 5. The same as Fig. 4, except for the stable static solutions shown in Fig. 2 (dark matter is a free Fermi gaswith fermion masses, m f , as indicated above each plot). In (a), note that the vertical axis has a log scale and thatthe red region is plotted above the black curve so that it is visible. In (b) note that the entire red region for a darkmatter core overlaps green for a dark matter halo. R om [km] M [ M - ] FIG. 6. The same as Fig. 4, except for the stable staticsolutions shown in Fig. 3 (dark matter is mirror darkmatter). why the critical curve extends past the single-fluidvalue of ( p c dm ) crit , but not ( p c om ) crit , first noticed inSec. III. With the small fermion mass of m f = 0 . p c dm substantially isdark matter able to affect the system and cause itto be unstable. Even then, for large p c om , the or-dinary matter is always dominating and we do notfind stability past the single-fluid value of ( p c om ) crit .Now consider Fig. 8(b), where we see some,though not significant, changes compared to 8(a).That the changes are not significant is expected, be-cause for both m f = 0 . p c SLy ; p c Fermi = ( m f = [MeV/fm ] ! S L y ; ! F e r m i = ( m f = G e V ) [ k H z ] FIG. 7. The fundamental radial oscillation frequency, ω ,is plotted as a function of the central pressure for single-fluid stars using the SLy (solid black curve) and freeFermi gas (dashed blue curve) equations of state. Thecurves are seen to hit zero at the critical central pressuresgiven in Eq. (7). For SLy, the maximum frequency is ω =19 .
96 kHz, which occurs for a central pressure of p c =36 .
02 MeV/fm . For the free Fermi gas, the maximumfrequency is ω = 7 . m f / kHz, which occurs fora central pressure of p c = 61 . m f / MeV/fm . mass, dark matter has a bigger influence over thefrequency. This connects with the fact that p c dm doesnot have to be raised as high before it makes the sys-tem go unstable and that for large p c dm , we can push p c om past its single-fluid critical value and still havestable solutions.We do see significant changes in Fig. 8(c). Thisis expected, since for m f = 2 GeV the single-fluidfrequency for a free Fermi gas can be larger than8 a) m f = 0 : p c om [MeV/fm ] ) l og ( p c d m = ( m f = G e V ) [ M e V / f m ] ) (b) m f = 1 GeVlog( p c om [MeV/fm ] ) (c) m f = 2 GeVlog( p c om [MeV/fm ] ) ! [ k H z ] FIG. 8. The fundamental radial oscillation frequency, ω , is plotted as a function of the central pressures p c om and p c dm , where dark matter is taken to be a free Fermi gas with fermion masses, m f , as indicated above each plot. Thethick black lines are critical curves, first shown in Figs. 1 and 2. the single-fluid frequency for ordinary matter in Fig.7. The larger frequencies in Fig. 8(c) coming fromdark matter, in a sense, “collide” with the smallerfrequencies coming from normal matter in the centerof the figure. Interestingly, this causes the frequencyto increase, as can be seen by the darker region nearthe center, to values larger than the maximum pos-sible single-fluid frequencies.In Fig. 9 we show results for mirror dark mat-ter. The thick black line is the critical curve firstshown in Fig. 3. The figure is symmetric in param-eter space, since the same equation of state is usedfor dark matter as is used for ordinary matter. Sim-ilar to Fig. 8(c), we find a “collision” in the centerof the plot with a frequency that is larger than themaximum single-fluid frequency given by the solidblack line in Fig. 7. VI. CONCLUSION
We studied dark matter admixed neutron stars,which are two-fluid systems with only gravitationalinter-fluid interactions. The first fluid describes or-dinary nuclear matter and the second fluid describesdark matter. We considered two possibilities fordark matter: a free Fermi gas and mirror dark mat-ter. Static solutions were found by solving the two-fluid TOV equations.Our study focused on three computations. Thefirst was the stability of static solutions with respectto small perturbations. Rigorous determinations ofstability over large swaths of parameter space werelacking in the literature. We presented two differentways to determine stability and computed critical log( p c om [MeV/fm ] ) l og ( p c d m [ M e V / f m ] ) ! [ k H z ] FIG. 9. The same as Fig. 8, except for mirror dark mat-ter, in which dark matter has the same equation of stateas ordinary matter. The thick black line is the criticalcurve, first shown in Fig. 3. curves, which separate stable solutions from unsta-ble ones in parameter space. Interestingly, we foundstable regions of parameter space for which a naiveanalysis of the individual equations of state wouldnot have deemed stable.The second computation was for mass-radius re-lations from static solutions. As an alternative towhat is commonly presented in the literature, wegave mass-radius diagrams over the whole of the sta-ble static parameter space, highlighting when darkmatter acts as a core or as a halo in the star.The third computation was the radial oscillationfrequency. As with stability, computations of theradial oscillation frequency over large swaths of pa-rameter space were lacking in the literature. Inter-9stingly, our results showed that the frequencies ofdark matter admixed neutron stars could be largerthan the maximum possible frequencies of single-fluid stars made with the individual equations ofstate.
Appendix A: Radial oscillations and stability
In Sec. III, we discussed two methods for comput-ing the critical curve, which separates stable staticsolutions from unstable ones in parameter space. InSec. V, we computed radial oscillation frequencies.In this appendix we review the equations used inthese computations.
1. Radial oscillations
The method we use to solve for the radial oscil-lation frequencies was derived in [43]. We brieflyreview the equations and how they are solved hereand refer the reader to [43] for details. To make atime-dependent perturbation around a static solu-tion, we must first write the energy-momentum ten-sor, T µν = (cid:80) i T µνi , in terms of the full perfect fluidform, T µνi = ( (cid:15) i + p i ) u µi u νi + p i g µν , (A1)where u µi is the four-velocity of the fluid, and notin terms of the static form, as was done in Eq. (2).Spherical symmetry sets u θ = u φ = 0 and we define v i ≡ e ν/ u ri , where ν is the metric function in (1)(but now with a time dependence). We can thenwrite the metric functions, energy densities, andpressures as perturbations about their static solu-tions and then write the Einstein field equations andequations of motion to first order in the perturba- tions. Note that v i is at the level of a perturbation,since it vanishes in the static limit.Defining the quantity ξ i through ∂ t ξ i ≡ v i , theperturbations are all taken to be of harmonic form, ξ i ( t, r ) = ξ i ( r ) e iωt , (A2)which defines the radial oscillation frequency, ω . Wefurther define ζ i ( r ) ≡ r e − ν ( r ) / ξ i ( r ) , (A3)where a subscripted 0 in this section refers to a staticsolution. The idea is to combine the Einstein fieldequations and equations of motion such that we ob-tain a system of pulsation equations which dependon ζ i and its derivatives, and not on any other per-turbations.Before presenting the pulsation equations, werewrite the metric in (1) as ds = − H ( t, r ) σ ( t, r ) dt + dr H ( t, r ) + r d Ω , (A4)where H ( t, r ) ≡ − m ( t, r ) r , σ ( t, r ) ≡ e ν ( t,r ) / (cid:112) H ( t, r ) , (A5)which is better suited for numerical solutions. Theequation for the static metric function σ ( r ) is dσ dr = 4 πrσ H ( (cid:15) + p ) , (A6)which follows from the Einstein field equations andis one of the TOV equations, but was not listed in(3) because, for static solutions, it decouples and isnot needed.The system of pulsation equations is [43] ∂ r ( (cid:98) Π ζ (cid:48) i ) + ( (cid:98) Q i + ˆ ω W i )ˆ ζ i + (cid:98) R (cid:18) (cid:15) i + p i r − p (cid:48) i (cid:19) (cid:88) j ( (cid:15) j + p j )ˆ ζ j + r ( (cid:15) i + p i )ˆ σ H (cid:88) j η j = (cid:98) S i (cid:88) j ( (cid:15) j + p j ) (cid:16) ˆ ζ j − ˆ ζ i (cid:17) + r ˆ σ H (cid:98) R ( (cid:15) i + p i ) (cid:88) j (cid:88) k p j γ j ( (cid:15) k + p k ) (cid:16) ˆ ζ k − ˆ ζ j (cid:17) + (cid:98) Rγ i p i (cid:88) j (cid:104) ( (cid:15) (cid:48) j + p (cid:48) j ) (cid:16) ˆ ζ j − ˆ ζ i (cid:17) + ( (cid:15) j + p j ) (cid:16) ˆ ζ (cid:48) j − ˆ ζ (cid:48) i (cid:17)(cid:105) , (A7)where a prime denotes an r derivative, where p (cid:48) i is given by the TOV equation in (3), and where (cid:98) Π i = 1 r p i γ i ˆ σ H W i = 1 r H ( (cid:15) i + p i ) 10 Q i = − ˆ σ H r (cid:26) r p (cid:48) i + (cid:20) πH p ( (cid:15) i + p i ) + (cid:18) πrH (cid:15) − m r H (cid:19) (cid:18) (cid:15) i + p i r − p (cid:48) i (cid:19)(cid:21)(cid:27)(cid:98) R = 4 π ˆ σ r (cid:98) S i = (cid:98) R (cid:26) ( γ i − p (cid:48) i + γ (cid:48) i p i + γ i p i (cid:20) πrH ( (cid:15) + p ) − r (cid:21)(cid:27) γ i = (cid:18) (cid:15) i p i (cid:19) ∂p i ∂(cid:15) i . (A8)Those quantities with a hat have been scaled by pow-ers of σ c , the central value of σ . This has the ef-fect of changing the boundary conditions and mak-ing the equations easier to solve. Note that the righthand side of Eq. (A7) vanishes for a single fluid, inwhich case the left hand side is equivalent to Chan-drasekhar’s pulsation equation [53]. Though equiv-alent, the left hand side of Eq. (A7) is not writtenin an identical form to Chandrasekhar’s because, inthe presence of multiple fluids, terms cannot canceland combine in the same way.To actually solve the pulsation equations, we de-fine η i ≡ (cid:98) Π i ˆ ζ (cid:48) i and solve the system of first order dif-ferential equations made up of Eq. (A7), but writtenin terms of η i , and ˆ ζ (cid:48) i = η i (cid:98) Π i . (A9)In the two-fluid case, the inner boundary conditionsare m = ζ i = 0 and ˆ σ = η = 1. The outerboundary conditions for the static variables are asdiscussed in Sec. II and for the perturbations are η i ( R i ) = 0. The undetermined parameters are theinner boundary condition for η and the value ofscaled squared radial oscillation frequency ˆ ω . Theseare determined using the shooting method. Once asolution to the pulsation equations is found, so tooare σ c , which follows from the outer value of ˆ σ , andˆ ω . The squared radial oscillation frequency is thengiven by ω = ( σ c ˆ ω ) .
2. Stability
In Sec. III, we discussed two methods for deter-mining whether a solution to the TOV equations isstable with respect to small perturbations. The firstmethod is to compute the squared radial oscillationfrequency using the methods of the previous sub-section. Importantly, one must compute oscillation frequencies for the fundamental solution, which hasthe smallest frequency, and not for excited solutions.In practice, this is easily done by making sure thatthe η i do not have nodes, i.e. values for which theyequal zero before the edge of their respective fluid.If ω < ω >
0, the corresponding static solution is stable.It is possible for the squared frequency of the fun-damental solution to be negative while the excitedfrequencies are positive, which is why it is the fun-damental frequency that must be computed. In thispaper, including in Sec. V, we only compute radialoscillation frequencies for fundamental solutions.The transition from stable to unstable for staticsolutions occurs at those points in parameter spacewhere ω = 0. Such points map out the criticalcurves shown in Figs. 1–3. In this way, the com-putation of the squared radial oscillation frequencyusing the methods of the previous subsection can beused to compute critical curves.The second method we use for computing the crit-ical curve is from [29]. The idea is straightforward.We still wish to find those points in parameter spacewhere ω = 0. From Eq. (A2) we find that, in thiscase, the perturbation is time independent and thusthe perturbed static solution is itself a static solu-tion. This is convenient, since it means we needonly deal with static solutions. Further, the pertur-bations (including those with time-dependence) pre-serve the total fluid number, N i , and total mass, M ,of the system. Putting these facts together, if thereis a point in the parameter space of static solutionswhere ω = 0, then there must be some direction inparameter space, given by vector p , that preserves ω = 0 as well as the total fluid numbers and mass ofthe system. This immediately gives Eq. (8), which iswhat we used to compute the critical curves in Sec.III. [1] I. Goldman and S. Nussinov, Weakly InteractingMassive Particles and Neutron Stars, Phys. Rev. D , 3221 (1989).
2] C. Kouvaris, WIMP Annihilation and Cooling ofNeutron Stars, Phys. Rev. D , 023006 (2008),arXiv:0708.2362 [astro-ph].[3] G. Bertone and M. Fairbairn, Compact Stars asDark Matter Probes, Phys. Rev. D , 043515(2008), arXiv:0709.1485 [astro-ph].[4] A. de Lavallaz and M. Fairbairn, Neutron Starsas Dark Matter Probes, Phys. Rev. D , 123521(2010), arXiv:1004.0629 [astro-ph.GA].[5] C. Kouvaris and P. Tinyakov, Can Neutron starsconstrain Dark Matter?, Phys. Rev. D , 063531(2010), arXiv:1004.0586 [astro-ph.GA].[6] R. Brito, V. Cardoso, and H. Okawa, Accretion ofdark matter by stars, Phys. Rev. Lett. , 111301(2015), arXiv:1508.04773 [gr-qc].[7] M. Cerme˜no, M. ´A. P´erez-Garc´ıa and J. Silk,Fermionic Light Dark Matter Particles and theNew Physics of Neutron Stars, Publ. Astron. Soc.Austral. , e043 (2017), arXiv:1710.06866 [astro-ph.HE].[8] M. I. Gresham and K. M. Zurek, Asymmetric DarkStars and Neutron Star Stability, Phys. Rev. D ,083008 (2019), arXiv:1809.08254 [astro-ph.CO].[9] J. Ellis, G. H¨utsi, K. Kannike, L. Marzola,M. Raidal, and V. Vaskonen, Dark Matter EffectsOn Neutron Star Properties, Phys. Rev. D ,123007 (2018), arXiv:1804.01418 [astro-ph.CO].[10] M. Deliyergiyev, A. Del Popolo, L. Tolos, M. Le Del-liou, X. Lee, and F. Burgio, Dark compact objects:an extensive overview, Phys. Rev. D , 063015(2019), arXiv:1903.01183 [gr-qc].[11] A. Nelson, S. Reddy, and D. Zhou, Dark ha-los around neutron stars and gravitational waves,JCAP , 012, arXiv:1803.03266 [hep-ph].[12] D. E. Kaplan, M. A. Luty, and K. M. Zurek,Asymmetric Dark Matter, Phys. Rev. D , 115016(2009), arXiv:0901.4117 [hep-ph].[13] K. M. Zurek, Asymmetric Dark Matter: Theories,Signatures, and Constraints, Phys. Rept. , 91(2014), arXiv:1308.0338 [hep-ph].[14] L. Okun, Mirror particles and mirror matter: 50years of speculations and search, Phys. Usp. , 380(2007), arXiv:hep-ph/0606202.[15] S. I. Blinnikov and M. Yu. Khlopov, “On the possi-ble signatures of mirror particles,” Yadernaya Fizika , 809 (1982) [Sov. J .Nucl. Phys.
472 (1982)].[16] S. I. Blinnikov and M. Yu. Khlopov, “On the pos-sible astronomical effects of ‘mirror’ particles,” As-tron. Zh.
632 (1983) [Sov. Astron.
371 (1983)].[17] S. I. Blinnikov and M. Yu. Khlopov, “Excitation ofthe Solar oscillations by the objects consisting ofy-matter,” Solar physics , 383 (1983).[18] M. Yu. Khlopov, G. M. Beskin, N. G. Bochkarev,L. A. Pustilnik and S. A. Pustilnik, “Observationalphysics of mirror world,” Astron. Zh.
121 (1991)[Sov. Astron.
21 (1991)][19] F. Sandin and P. Ciarcelluti, Effects of mirror darkmatter on neutron stars, Astropart. Phys. , 278(2009), arXiv:0809.2942 [astro-ph].[20] P. Ciarcelluti and F. Sandin, Have neutron stars adark matter core?, Phys. Lett. B , 19 (2011), arXiv:1005.0857 [astro-ph.HE].[21] I. Goldman, Implications of Mirror Dark Matter onNeutron Stars, Acta Phys. Polon. B , 2203 (2011),arXiv:1112.1505 [astro-ph.CO].[22] I. Goldman, R. Mohapatra, S. Nussinov, D. Rosen-baum, and V. Teplitz, Possible Implications ofAsymmetric Fermionic Dark Matter for Neu-tron Stars, Phys. Lett. B , 200 (2013),arXiv:1305.6908 [astro-ph.CO].[23] T. Kodama and M. Yamada, Theory of superdensestars, Prog. Theor. Phys. , 444 (1972).[24] A. B. Henriques, A. R. Liddle, and R. G. Moor-house, Combined Boson-Fermion Stars, Phys. Lett.B , 99 (1989).[25] D. Akerib et al. (LUX), Limits on spin-dependentWIMP-nucleon cross section obtained from the com-plete LUX exposure, Phys. Rev. Lett. , 251302(2017), arXiv:1705.03380 [astro-ph.CO].[26] X. Cui et al. (PandaX-II), Dark Matter ResultsFrom 54-Ton-Day Exposure of PandaX-II Ex-periment, Phys. Rev. Lett. , 181302 (2017),arXiv:1708.06917 [astro-ph.CO].[27] E. Aprile et al. (XENON), Dark Matter SearchResults from a One Ton-Year Exposure ofXENON1T, Phys. Rev. Lett. , 111302 (2018),arXiv:1805.12562 [astro-ph.CO].[28] A. B. Henriques, A. R. Liddle, and R. G. Moor-house, Combined Boson-Fermion Stars: Configura-tions and Stability, Nucl. Phys. B , 737 (1990).[29] A. B. Henriques, A. R. Liddle, and R. G. Moor-house, Stability of boson-fermion stars, Phys. Lett.B , 511 (1990).[30] S. Leung, M. Chu, and L. Lin, Dark-matter admixedneutron stars, Phys. Rev. D , 107301 (2011),arXiv:1111.1787 [astro-ph.CO].[31] S. Leung, M. Chu, and L. Lin, Equilibrium Struc-ture and Radial Oscillations of Dark Matter Ad-mixed Neutron Stars, Phys. Rev. D , 103528(2012), arXiv:1205.1909 [astro-ph.CO].[32] S.-C. Leung, M.-C. Chu, L.-M. Lin, and K.-W.Wong, Dark-matter admixed white dwarfs, Phys.Rev. D , 123506 (2013), arXiv:1305.6142 [astro-ph.CO].[33] A. Li, F. Huang, and R.-X. Xu, Too massive neutronstars: The role of dark matter?, Astropart. Phys. , 70 (2012), arXiv:1208.3722 [astro-ph.SR].[34] Q.-F. Xiang, W.-Z. Jiang, D.-R. Zhang, and R.-Y.Yang, Effects of fermionic dark matter on propertiesof neutron stars, Phys. Rev. C , 025803 (2014),arXiv:1305.7354 [astro-ph.SR].[35] X. Y. Li, F. Y. Wang, and K. S. Cheng, Gravita-tional effects of condensate dark matter on com-pact stellar objects, JCAP , 031, arXiv:1210.1748[astro-ph].[36] L. Tolos and J. Schaffner-Bielich, Dark Com-pact Planets, Phys. Rev. D , 123002 (2015),arXiv:1507.08197 [astro-ph.HE].[37] P. Mukhopadhyay and J. Schaffner-Bielich, Quarkstars admixed with dark matter, Phys. Rev. D ,083009 (2016), arXiv:1511.00238 [astro-ph.HE].[38] G. Panotopoulos and I. Lopes, Gravitational effects f condensed dark matter on strange stars, Phys.Rev. D , 023002 (2017), arXiv:1706.07272 [gr-qc].[39] G. Panotopoulos and I. Lopes, Dark matter effecton realistic equation of state in neutron stars, Phys.Rev. D , 083004 (2017), arXiv:1709.06312 [hep-ph].[40] S. A. Bhat and A. Paul, “Cooling of Dark-MatterAdmixed Neutron Stars with density-dependentEquation of State,” Eur. Phys. J. C , 544 (2020),arXiv:1905.12483 [hep-ph].[41] A. Del Popolo, M. Deliyergiyev, and M. Le Del-liou, Solution to the hyperon puzzle using darkmatter, Phys. Dark Univ. , 100622 (2020),arXiv:2011.00984 [gr-qc].[42] K. Zhang and F.-L. Lin, Constraint on hybrid starswith gravitational wave events, Universe , 231(2020), arXiv:2011.05104 [astro-ph.HE].[43] B. Kain, Radial oscillations and stability ofmultiple-fluid compact stars, Phys. Rev. D ,023001 (2020), arXiv:2007.04311 [gr-qc].[44] A. Brillante and I. N. Mishustin, Radial oscillationsof neutral and charged hybrid stars, EPL , 39001(2014), arXiv:1401.7915 [astro-ph.SR].[45] P. Haensel and A. Y. Potekhin, Analytical repre-sentations of unified equations of state of neutron-star matter, Astron. Astrophys. , 191 (2004),arXiv:astro-ph/0408324.[46] F. Douchin and P. Haensel, A unified equation ofstate of dense matter and neutron star structure,Astron. Astrophys. , 151 (2001), arXiv:astro-ph/0111092.[47] S. L. Shapiro and S. A. Teukolsky, Black Holes,White Dwarfs and Neutron Stars: The Physics ofCompact Objects (John Wiley & Sons, 1983).[48] N. K. Glendenning,
Compact Stars: NuclearPhysics, Particle Physics, and General Relativity,2nd ed. (Springer, New York, 2000).[49] B. K. Harrison, K. S. Thorne, M. Wakano, and J. A.Wheeler,
Gravitation Theory and Gravitational Col-lapse (University of Chicago Press, Chicago, 1965).[50] G. Comer, D. Langlois, and L. M. Lin, Quasinor-mal modes of general relativistic superfluid neutronstars, Phys. Rev. D , 104025 (1999), arXiv:gr-qc/9908040.[51] P. Jetzer, Stability of Combined Boson - Fermion Stars, Phys. Lett. B , 36 (1990).[52] S. Valdez-Alvarado, C. Palenzuela, D. Alic, andL. A. Ure˜na-L´opez, Dynamical evolution of fermion-boson stars, Phys. Rev. D , 084040 (2013),arXiv:1210.2299 [gr-qc].[53] S. Chandrasekhar, The Dynamical Instability ofGaseous Masses Approaching the SchwarzschildLimit in General Relativity, Astrophys. J. , 417(1964), [Erratum: Astrophys.J. 140, 1342 (1964)].[54] C. W. Misner, K. Thorne, and J. Wheeler, Gravita-tion (W. H. Freeman, San Francisco, 1973).[55] G. Chanmugam, Radial Oscillations of Zero-Temperature White Dwarfs and Neutron Stars Be-low Nuclear Densities, Astrophys. J. , 799(1977).[56] D. Gondek, P. Haensel, and J. Zdunik, Radialpulsations and stability of protoneutron stars, As-tron. Astrophys. , 217 (1997), arXiv:astro-ph/9705157.[57] K. Kokkotas and J. Ruoff, Radial oscillations of rel-ativistic stars, Astron. Astrophys. , 565 (2001),arXiv:gr-qc/0011093.[58] E. N. Glass and L. Lindblom, The Radial Oscilla-tions of Neutron Stars, Astrophys. J. Supp , 93(1983).[59] O. G. Benvenuto and J. E. Horvath, Radial pulsa-tions of strange stars and the internal compositionof pulsars, Mon. Not. Roy. Astron. Soc. , 679(1991).[60] C. Vasquez Flores and G. Lugones, Radialoscillations of color superconducting self-boundquark stars, Phys. Rev. D , 063006 (2010),arXiv:1008.4882 [astro-ph.HE].[61] F. Di Clemente, M. Mannarelli, and F. Tonelli,Reliable description of the radial oscillations ofcompact stars, Phys. Rev. D , 103003 (2020),arXiv:2002.09483 [gr-qc].[62] G. Panotopoulos and I. Lopes, Radial oscillationsof strange quark stars admixed with condenseddark matter, Phys. Rev. D , 083013 (2017),arXiv:1709.06643 [gr-qc].[63] G. Panotopoulos and I. Lopes, Radial oscillationsof strange quark stars admixed with fermionic darkmatter, Phys. Rev. D , 083001 (2018)., 083001 (2018).