Radial Oscillations of Quark Stars Admixed with Dark Matter

Abstract

We investigated compact stars consisting of cold quark matter and fermionic dark matter treated as two admixed fluids. We computed the stellar structures and fundamental radial oscillation frequencies of different masses of the dark fermion in the cases of weak and strong self-interacting dark matter. We found that the fundamental frequency can be dramatically modified and, in some cases, stable dark strange planets and dark strangelets with very low masses and radii can be formed.

1. Introduction

Compact stars offer a variety of possibilities for probing the inner structure of matter through astronomical observations. In particular, matter at extremely high densities can only be probed, so far, by investigating unique objects that represent one of the possible final stages of stellar evolution. The structure of compact stars can be determined by solving the Tolman–Oppenheimer–Volkov (TOV) equations, given the equation of state (EoS) for the matter under consideration [1,2].

For high enough central energy densities, one expects to find either hybrid stars, i.e., neutron stars with a quark matter (QM) core, or even more exotic objects, such as quark stars. Quark stars [3] and their structure [4] have been considered for more than half a century, even before the elaboration of quantum chromodynamics (QCD) in the 1970s. Later, after a seminal work by Witten [5], a rich phenomenology of self-bound strange stars [6,7] and quark (hybrid) stars emerged using the MIT bag model [8] as a framework for the EoS at high densities. For a review on quark matter in neutron stars, see Ref. [9].

On the other hand, dark matter (DM) represents about a quarter of the total mass–energy density content of the universe or, equivalently, ∼85% of its matter content. Apart from this, it is needed to explain structure formations without modifying general relativity in the current cosmological standard model [10,11,12]. Nevertheless, there is still no experimental evidence of DM-constituent particles, and its nature remains one of the greatest mysteries of particle physics. Over the years, many candidates have been proposed as being DM-constituent particles, with masses ranging from ${10}^{-33}$ GeV to ${10}^{15}$ GeV, including weakly interacting massive particles (WIMPs), axions and axion-like particles (ALPs), sterile neutrinos, neutralinos, and so on [13]. In spite of their nature, if DM-constituent particles do not self annihilate and are non-relativistic at freeze-out (cold dark matter (CDM)), the probability of their interaction with ordinary baryonic matter will increase within the extreme densities found in compact stars. In this case, DM can accumulate and thermalize in a small radius. So, if quark stars are to be found in the universe, they have most likely accumulated some amount of dark matter over the course of their lives.

In this paper, we investigate strange quark stars consisting of cold QM and non self-annihilating fermionic cold DM treated as two admixed fluids, attracted only gravitationally. As our main goal, we aim to compute the fundamental radial oscillation frequency of admixed QM and DM two-fluid stars for a relevant range of masses of the dark fermion. We consider the cases of weak and strong self-interacting DM and also study how the total mass and radius of quark and dark stars are modified by their mutual presence in the admixed star.

We describe the QM component on the framework of the MIT bag model, which represents a choice for simplicity, which was mainly motivated by the possibility of direct comparison to previous work. Moreover, for analyses that depend on a range of values for the dark fermion mass and the intensity of DM self-interaction, it is convenient to avoid other bands in parameters that would come about naturally in more realistic descriptions of the EoS, such as those relying on perturbative QCD [14,15,16,17,18,19,20,21,22,23,24].

Neutron stars and quark stars admixed with DM have previously been considered. The effects of fermionic and bosonic DM on the equilibrium features and radial oscillations of neutron stars (NS) have been discussed, e.g., in Refs. [25,26,27,28,29,30,31] (see Ref. [32] for a more complete list of references). Reference [33] considers hybrid NSs with an EoS for neutron star matter that uses perturbative QCD and effective field theory as high and low-density descriptions, respectively, and polytropes as interpolating functions, as discussed in Ref. [18]1, in addition to taking into account inner and outer crusts. The authors also consider white dwarfs admixed with asymmetric DM and find dark, compact (Jupiter-like) planets and limits on the DM content of the stars in order to satisfy the two-solar mass observational constraint [34,35]. A stability analysis is also performed by solving the usual Sturm–Liouville problem for one-fluid stars [36,37]. Reference [13] extends these results to a wider range of dark fermion masses, from 1 GeV to 500 GeV, assuming different amounts of DM at the stellar center. From this analysis, the authors inferred that the total mass decreases with $m_D$, putting constraints on $m_D$ and on DM capture.

Quark stars admixed with dark matter have been discussed in Ref. [38], where the authors use the MIT bag model to describe the EoS for QM admixed with DM made of dark fermions of mass $m_D=100$ GeV (on the typical WIMP mass scale). They considered two cases: free and strongly self-interacting DM. Solving the TOV equations with two fluids that interact only gravitationally, they found minor modifications to the maximum mass and radius, of the order of a few percent, though with higher values of the central energy density due to a greater gravitational pull.

So far, the stability against radial oscillations of quark stars admixed with dark matter has been studied using one-fluid formalism, usually with simplified interaction terms and unphysical dark-fermion masses [39,40]. A complete treatment of the stability of two-fluid stars requires a non-trivial extension of the Chandrasekhar second-order differential equation, where a coupled system of equations should be solved for the corresponding Lagrangian displacements associated with each fluid, but which also depend on the displacement of the other fluid [26,27,41,42], as will be discussed in the sequel.

The paper is organized as follows. In Section 2, we briefly describe the two-fluid hydrostatic equilibrium equations together with the general-relativistic formalism used to study radial pulsations of the admixed stars. Section 3 contains our main results and discussion. Section 4 presents our summary and perspectives. We adopt natural units, i.e., $\hslash =c=1$.

2. Framework

In this section we summarize the main features of the TOV equations for the admixture of two fluids that interact only gravitationally. The one-fluid radial oscillation equations are conveniently partitioned to analyze either the stability of the quark mater core or dark matter core of the whole compact star. For the bag constant we use $B^{1/4}=145$ MeV, which implements the Bodmer–Witten–Terazawa hypothesis for strange quark matter [2]. This choice yields a maximum mass of $2.01M_{\odot }$ and a radius of $R=11$ km. As mentioned previously, fermionic dark matter is considered as being either weakly ($y=0.1$) or strongly ($y={10}^3$) self interacting [43]. Here, $y\equiv m_D/m_{int}$, where $m_D$ is the dark fermion mass, and $m_{int}$ is the scale of interaction. One can consider that $m_{int}\sim 100$ MeV for strong interactions and $m_{int}\sim 300$ GeV for weak interactions. The values $y=0.1$ and $y={10}^3$ are commonly adopted as typical illustrations.

2.1. Two-Fluid Hydrostatic Equilibrium Equations

Since the structure of spherically symmetric, static one-fluid compact stars is determined from the usual TOV equations, they can be separated in order to deal with two fluids that only interact gravitationally. This can be performed as follows. The (perfect) one-fluid energy-momentum tensor is divided into two parts, i.e., $T^{\mu \nu }=T_1^{\mu \nu }+T_2^{\mu \nu }$. This, in turn, induces a separation of the total pressure and energy density in fluid components, as $p=p_1+p_2$ and $ϵ={ϵ}_1+{ϵ}_2$, respectively. Given that this separation does not affect the temporal $e^{\nu \left(r\right)}$ or radial $e^{\lambda \left(r\right)}$ metric functions, the corresponding equations keep their original forms, depending on the total pressure and energy densities. On the other hand, one has a set of coupled TOV equations for each of the fluids.

In our case, i.e., quark and dark matter fluids, we have these two-fluid TOV equations in its dimensionless form given by [33] (see Ref. [29] for a detailed variational derivation):

$$\begin{array}{lll}& & \frac{d{p}_{\mathrm{QM}}^{\prime }}{d{r}^{\prime }}=-\frac{(p_{\mathrm{QM}}^{\prime }+{ϵ}_{\mathrm{QM}}^{\prime })}{2}\frac{d\nu }{d{r}^{\prime }},\\ & & \frac{d{m}_{\mathrm{QM}}^{\prime }}{d{r}^{\prime }}=4\pi r^{\prime 2}{ϵ}_{\mathrm{QM}}^{\prime },\\ & & \frac{d{p}_{\mathrm{DM}}^{\prime }}{d{r}^{\prime }}=-\frac{(p_{\mathrm{DM}}^{\prime }+{ϵ}_{\mathrm{DM}}^{\prime })}{2}\frac{d\nu }{d{r}^{\prime }},\\ & & \frac{d{m}_{\mathrm{DM}}^{\prime }}{d{r}^{\prime }}=4\pi r^{\prime 2}{ϵ}_{\mathrm{DM}}^{\prime },\\ & & \frac{d\nu }{d{r}^{\prime }}=2\frac{(m_{\mathrm{QM}}^{\prime }+m_{\mathrm{DM}}^{\prime })+4\pi r^{\prime 3}(p_{\mathrm{QM}}^{\prime }+p_{\mathrm{DM}}^{\prime })}{r^{\prime }(r^{\prime }-2(m_{\mathrm{QM}}^{\prime }+m_{\mathrm{DM}}^{\prime }))},\end{array}$$

(1)

where $p^{\prime }$ and ${ϵ}^{\prime }$ are the dimensionless pressure and energy density, respectively, and $m_{\mathrm{QM},\mathrm{DM}}^{\prime }$ are the dimensionless gravitational masses enclosed inside the dimensionless radial coordinate $r^{\prime }$.

The set of equations above is solved simultaneously by specifying the (dimensionless) EoSs for QM, i.e., $p_{\mathrm{QM}}^{\prime }=p_{\mathrm{QM}}^{\prime }\left({ϵ}_{\mathrm{QM}}^{\prime }\right)$, and DM, i.e., $p_{\mathrm{DM}}^{\prime }=p_{\mathrm{DM}}^{\prime }\left({ϵ}_{\mathrm{DM}}^{\prime }\right)$. As usual, the conditions at the center should be given for QM and DM in the admixed star. The numerical integration stops when one of the pressures reaches zero, i.e., $p_{\mathrm{QM}/\mathrm{DM}}^{\prime }\left(R_{\mathrm{QM}/\mathrm{DM}}^{\prime }\right)=0$, characterizing the QM or DM core surface, allowing us to obtain the corresponding gravitational mass $m_{\mathrm{QM}/\mathrm{DM}}^{\prime }\left(R_{\mathrm{QM}/\mathrm{DM}}^{\prime }\right)=M_{\mathrm{QM}/\mathrm{DM}}^{\prime }$, where, in general, $R_{\mathrm{QM}}^{\prime }\ne R_{\mathrm{DM}}^{\prime }$. If $R_{\mathrm{QM}}^{\prime }>R_{\mathrm{DM}}^{\prime }$, the admixed star has a DM core and, if $R_{\mathrm{QM}}^{\prime }<R_{\mathrm{DM}}^{\prime }$, it has a DM halo surrounding a QM core.2

The boundary conditions for the metric function $\nu \left(r\right)$ come from ensuring that it matches the Schwarzschild metric outside the QM or DM core in the admixed star, i.e.,:

$$\nu \left(R_{\mathrm{QM}}^{\prime }\right)=ln\left(1-\frac{2(M_{\mathrm{QM}}^{\prime }+m_{\mathrm{DM}}^{\prime }\left(R_{\mathrm{QM}}^{\prime }\right))}{R_{\mathrm{QM}}^{\prime }}\right)$$

(2)

or

$$\nu \left(R_{\mathrm{DM}}^{\prime }\right)=ln\left(1-\frac{2(m_{\mathrm{QM}}^{\prime }\left(R_{\mathrm{DM}}^{\prime }\right)+M_{\mathrm{DM}}^{\prime })}{R_{\mathrm{DM}}^{\prime }}\right),$$

(3)

respectively.

2.2. Pulsations of Quark and Dark Matter Cores

The equations that describe the radial pulsations of one-fluid compact stars were obtained for the first time by S. Chandrasekhar [36]. He found that these equations could be arranged as a Sturm–Liouville problem where the eigenvalues are the oscillation frequencies squared, ${\omega }^2$ (the eigenfunctions being the radial Lagrangian displacements). For numerical purposes, these equations can be conveniently modified to a pair of first-order differential equations for each of the Lagrangian variables with more intuitive boundary conditions [44,45,46,47,48].

Strictly speaking, the dynamic stability of admixed stars must be studied using the full two-fluid formalism3 of Refs. [27,41], which would produce unified oscillation frequencies for the entire admixed star. However, since this calculation is computationally very expensive and time consuming, we decided to solve the equivalent of the two-fluid TOV equations, realized instead in the form of oscillation equations. In order to perform that function, we used the formalism of Ref. [47], which deals with the relative radial displacement $\Delta r^{\prime }/r^{\prime }\equiv {\xi }^{\prime }=\xi$ and the Lagrangian perturbation pressure $\Delta p^{\prime }$, both dimensionless.

Inspired by the previous separation for the total pressure and energy density, we separated the total Lagrangian variables as $\xi ={\xi }_{\mathrm{QM}}+{\xi }_{\mathrm{DM}}$ and $\Delta p^{\prime }=\Delta p_{\mathrm{QM}}^{\prime }+\Delta p_{\mathrm{DM}}^{\prime }$ (omitting the term $e^{i\omega t}$ in both variables), obtaining the following system of equations:

$$\frac{d{\xi }_{\mathrm{QM}/\mathrm{DM}}}{d{r}^{\prime }}\equiv -\frac{1}{r^{\prime }}\left(3{\xi }_{\mathrm{QM}/\mathrm{DM}}+\frac{\Delta p_{\mathrm{QM}}^{\prime }}{\mathsf{\Gamma }{p}^{\prime }}\right)-\frac{d{p}^{\prime }}{d{r}^{\prime }}\frac{{\xi }_{\mathrm{QM}/\mathrm{DM}}}{(p^{\prime }+{ϵ}^{\prime })},$$

(4)

$$\begin{array}{c}\frac{d\Delta p_{\mathrm{QM}/\mathrm{DM}}^{\prime }}{d{r}^{\prime }}\equiv {\xi }_{\mathrm{QM}/\mathrm{DM}}\left\{{\omega }^{\prime 2}{e}^{\lambda -\nu }(p^{\prime }+{ϵ}^{\prime })r^{\prime }-4\frac{d{p}^{\prime }}{d{r}^{\prime }}\right\}+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\xi }_{\mathrm{QM}/\mathrm{DM}}\left\{{\left(\frac{d{p}^{\prime }}{d{r}^{\prime }}\right)}^2\frac{r^{\prime }}{(p^{\prime }+{ϵ}^{\prime })}-8\pi e^{\lambda }(p^{\prime }+{ϵ}^{\prime })p^{\prime }{r}^{\prime }\right\}+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\Delta p_{\mathrm{QM}/\mathrm{DM}}^{\prime }\left\{\frac{d{p}^{\prime }}{d{r}^{\prime }}\frac{1}{p^{\prime }+{ϵ}^{\prime }}-4\pi (p^{\prime }+{ϵ}^{\prime })r^{\prime }{e}^{\lambda }\right\},\end{array}$$

(5)

where ${\omega }^{\prime }$ is the dimensionless oscillation frequency and $\mathsf{\Gamma }$ is the adiabatic index4$\mathsf{\Gamma }=(1+{ϵ}^{\prime }/p^{\prime })(\partial p^{\prime }/\partial {ϵ}^{\prime })$. The metric function $\lambda \left(r^{\prime }\right)$ is obtained from $\lambda \left(r^{\prime }\right)=-ln(1-2(m_{\mathrm{QM}}^{\prime }\left(r^{\prime }\right)+m_{\mathrm{DM}}^{\prime }\left(r^{\prime }\right))/r^{\prime })$ with boundary conditions given by Equations (2) and (3), i.e., $\lambda \left(R_{\mathrm{QM}}^{\prime }\right)=-\nu \left(R_{\mathrm{QM}}^{\prime }\right)$ and $\lambda \left(R_{\mathrm{DM}}^{\prime }\right)=-\nu \left(R_{\mathrm{DM}}^{\prime }\right)$.

So far we have not mentioned whether ${\omega }^{\prime }$ corresponds to the pulsation of a QM core or a DM core. Recall that these equations represent a Sturm–Liouville problem, which defines its eigenvalues in terms of the associated boundary conditions. In this case, they are

$${(\Delta p_{\mathrm{QM}/\mathrm{DM}}^{\prime })}_{\mathrm{center}}\equiv -3{\left({\xi }_{\mathrm{QM}/\mathrm{DM}}\mathsf{\Gamma }{p}_{\mathrm{QM}/\mathrm{DM}}^{\prime }\right)}_{\mathrm{center}},$$

(6)

demanding smoothness at the QM or DM stellar center, and

$${(\Delta p_{\mathrm{QM}/\mathrm{DM}}^{\prime })}_{\mathrm{surface}}\equiv 0,$$

(7)

since $p_{\mathrm{QM}/\mathrm{DM}}^{\prime }\left(R_{\mathrm{QM}/\mathrm{DM}}^{\prime }\right)=0$, with eigenfunctions normalized to ${\xi }_{\mathrm{QM}/\mathrm{DM}}\left(0\right)=1$ as usual. Thus, Equations (6) and (7) lead us to define ${\omega }^{\prime 2}\to {\omega }_{\mathrm{QM}/\mathrm{DM}}^{\prime 2}$ if we are dealing with a QM/DM oscillating core in the admixed star. In other words, only one of the cores oscillates depending on the boundary conditions. The other fluid only affects its oscillation indirectly, by the coupling of the total pressure and energy density.

A word of caution should be added at this point. Usually, two ways of dealing with the radial oscillations of two-fluid compact stars have been explored. In the simplest one, only the radial oscillation of the whole admixed star, i.e., treated as one fluid, is studied without explicitly considering the gravitational coupling between the QM and DM cores (see e.g., Refs. [33,39,40]). In the second, a consistent general-relativistic formalism to deal with the couplings between oscillation amplitudes and Lagrangian perturbations for each fluid is developed [41,42]. Unfortunately, dealing with a system of highly coupled and non-linear differential equations requires very time-consuming numerical calculations which we consider unnecessary when independently solving the oscillation equations for each DM or QM core while keeping the other fluid at rest but still coupled through the coefficients entering in the equations. In this sense, the formalism built in this work occurs more in the line of Ref. [27], which considers the independent oscillations of each fluid, thus forming an Sturm–Liouville-like problem. Notice that our reasoning agrees with the fact that each of the two-fluid TOV equations can be considered an independent ’one-fluid’ star only, coupled through $\nu \left(r^{\prime }\right)$ to the other ’one-fluid’ star. Thus, radial oscillations of each ’one-fluid’ star can be associated with a set of one-fluid5 oscillation equations coupled now by the total pressures, energy densities, and polytropic indices and metric functions, $\nu \left(r^{\prime }\right)$ and $\lambda \left(r^{\prime }\right)$. For consistency, we have verified that our formalism agrees with the results of Ref. [38] when a delay of the maximal central density is reached at higher densities for increasing amounts of DM when the zero frequencies are reached. We stress that, in this case, $f_{n=0}\to 0$ coincides with $\partial M/\partial {ϵ}_c\to 0$, since few amounts of DM were considered, whereas in this work we explore all the available DM densities which notoriously modifies the stability of the admixed stars, so that $\partial M/\partial {ϵ}_c$ must be used with caution.

In following sections, we focus on the fundamental mode frequency, ${\omega }_{n=0}$. It vanishes at the maximal stable QM or DM mass configuration, marking the onset of the instability of the corresponding oscillating core which, in turn, induces the gravitational collapse of the whole admixed star.

3. Results and Discussion

The parameter space for quark stars admixed with weakly or strongly self-interacting DM is large. In this section we show only results where the effects on observables are relevant. As mentioned before, we considered dark fermion masses $m_D=1,10,50,100,200,500$ GeV in order to include all possible dark fermion candidates.

Regarding the numerical values we chose for ${ϵ}_c^{\mathrm{QM}/\mathrm{DM}}$ in our calculations and showed in our plots: (i) for QM, the three values of ${ϵ}_c^{\mathrm{QM}}$ correspond to somewhat above, twice, and nearly twice the value of the maximal central energy density of pure quark stars with $B={\left(145\mathrm{MeV}\right)}^4$, i.e., $\sim 1\mathrm{GeV}/{\mathrm{fm}}^3$. The reason for this is that higher values of ${ϵ}_c^{\mathrm{QM}}$ are required when DM is present in the admixed star; (ii) the three values of ${ϵ}_c^{\mathrm{DM}}$ (for strongly or weakly interacting DM) correspond to near the minimum, intermediate, and near the maximal-mass central densities for corresponding pure DM stars. In

Table 1. Maximum masses $M^{\max }\left(M_{\odot }\right)$ and their corresponding minimum radii $R^{\min }$ (km) and maximum central energy densities ${ϵ}_c^{\max }$ (GeV/fm${}^3$) obtained for weakly interacting (wDM) pure dark matter stars.

$\mathbf{EoSs}$	${\mathit{ϵ}}_{\mathit{c}}^{\mathbf{max}}\left(\mathbf{wDM}\right)$	${\mathit{M}}^{\mathbf{max}}\left(\mathbf{wDM}\right)$	${\mathit{R}}^{\mathbf{min}}\left(\mathbf{wDM}\right)$
$\mathrm{DM}(m_D/\mathrm{GeV}=1)$	∼3	∼0.63	∼8.1
10	∼$3\times {10}^4$	∼$6.27\times {10}^{-3}$	∼$7.8\times {10}^{-2}$
50	∼$1.97\times {10}^7$	∼$2.50\times {10}^{-4}$	∼$3.2\times {10}^{-3}$
100	∼$3\times {10}^8$	∼$6.27\times {10}^{-5}$	∼$8.1\times {10}^{-4}$
200	∼$4.99\times {10}^9$	∼$1.56\times {10}^{-5}$	∼$2\times {10}^{-4}$
500	∼$1.97\times {10}^{11}$	∼$2.50\times {10}^{-6}$	∼$3\times {10}^{-5}$

and

Table 2. Same notation as in Table 1 but now for strongly interacting (sDM) pure dark matter stars.

$\mathbf{EoSs}$	${\mathit{ϵ}}_{\mathit{c}}^{\mathbf{max}}\left(\mathbf{sDM}\right)$	${\mathit{M}}^{\mathbf{max}}\left(\mathbf{sDM}\right)$	${\mathit{R}}^{\mathbf{min}}\left(\mathbf{sDM}\right)$
$\mathrm{DM}(m_D/\mathrm{GeV}=1)$	∼$4.9\times {10}^{-5}$	∼$2.67\times {10}^2$	∼$1.87\times {10}^3$
10	∼0.6	∼2.67	∼18.5
50	∼$4\times {10}^2$	∼$1.07\times {10}^{-1}$	$\sim 7.4\times {10}^{-1}$
100	∼$4.9\times {10}^3$	∼$2.67\times {10}^{-2}$	∼$18.7\times {10}^{-2}$
200	∼$8\times {10}^4$	∼$6.68\times {10}^{-3}$	∼$4.7\times {10}^{-2}$
500	∼$4\times {10}^6$	∼$1.07\times {10}^{-3}$	∼$7.5\times {10}^{-3}$

the maximal-mass values of central energy density for each $m_D$ are listed. This choice was made to quantify the full dependence of the stellar structure on the amounts of DM. We will show that, in some cases with a huge amount of DM, only very small objects with strangelet-like and planet-like masses are allowed. This was expected from the results of Ref. [38].

Although the usual criterion for static stability, $\partial M/\partial {ϵ}_c\ge 0$, consistently works for one-fluid stars, it should not be taken for granted in two-fluid stars; only the frequency analysis can decide on their stability. Our results for the oscillation frequency of the fundamental mode for QM and DM cores in the admixed stars are written in terms of the linear frequency $f_{n=0}^{\mathrm{QM}/\mathrm{DM}}={\omega }_{n=0}^{\mathrm{QM}/\mathrm{DM}}/\left(2\pi \right)$.

3.1. Admixtures of Quark Matter and Weakly $(y=0.1)$ Interacting Dark Matter

3.1.1. Solving the Two-Fluid TOV Equations

In Figure 1, we display the results obtained from solving the two-fluid TOV Equation (2) with the condition $p_{\mathrm{QM}}\left(R_{\mathrm{QM}}\right)=0$ for different central energy densities of weakly interacting DM. One can easily see that only the solutions for $m_D=1500$ GeV display sizable modifications on the QM stellar masses and radii. In particular, the case of $m_D=500$ GeV suffers a marked reduction of $1.2M_{\odot }$ due to the very high DM central energy densities ($\sim {10}^{11}\mathrm{GeV}/{\mathrm{fm}}^3$). Additionally, the central QM densities increased considerably, by a factor ∼7. Normally such QM densities would generate unstable pure quark stars with central energy densities at most ∼1 $\mathrm{GeV}/{\mathrm{fm}}^3$ without the DM component.

On the other hand, solutions for other values of the dark fermion mass show negligible effects. Again, these QM stars require higher central energy densities in order to compensate for the extra gravitational pull from the DM. Furthermore, solutions for all the dark fermion masses, except $m_D=1$ and 500 GeV, develop a plateau at low QM stellar masses in the mass vs. central energy density plots, which became wider for higher DM central energy densities. Our calculations show that these QM cores have masses between ${10}^{-18}$ to ${10}^{-4}{M}_{\odot }$ with radii between ${10}^{-4}$ and ${10}^{-2}$ km, depending on the value of the dark fermion mass. For example, $m_D\sim 10$ GeV mostly commonly produces stellar masses around ${10}^{-5}{M}_{\odot }$ with radii of ${10}^{-3}$ km. As we increase the mass of the dark fermion, the values of $M_{\mathrm{QM}}$ and $R_{\mathrm{QM}}$ are reduced by many orders of magnitude. We note that all these stars satisfy the criterion $\partial M_{\mathrm{QM}}/\partial {ϵ}_c^{\mathrm{QM}}>0$ and can be tentatively considered stable objects, “dark strange planets” in analogy to the results of Ref. [33], and “dark strangelets”.

Figure 2 shows our results from the two-fluid TOV Equation (2) with the condition $p_{\mathrm{DM}}\left(R_{\mathrm{DM}}\right)=0$ for different values of QM central energy densities. The DM stars that are most affected by the presence of QM are the ones with $m_D=1$ and 10 GeV. For stars with $m_D=1$ GeV, the central QM energy densities are high enough to convert the usual behavior of pure DM stars in the mass–radius diagram into a self-bound-like behavior, making them more compact. To see this more quantitatively, Table 1 shows values of masses and radii for pure $y=0.1$ DM stars for the whole range of dark fermion masses considered. One can see that the stellar masses and radii are slightly affected by the presence of QM near the maximum mass, but the radii of less massive DM stars are significantly modified. The same is true in the case of $m_D=10$ GeV. For $m_D=50,100,200,500$ GeV, the DM high central energy densities completely dominate the QM contribution.

3.1.2. Solving the Coupled Radial Oscillation Equations

The solutions to the coupled radial pulsation Equations (4) and (5), assuming an oscillating QM core with fixed DM central energy density with boundary conditions (6) and (7), are shown in Figure 3. We show the zero-mode frequency as a function of central energy density and stellar mass.

One can see in Figure 3 that the increments in the DM central energy density tend to delay the onset of radial instability (except in the case of $m_D=10$ GeV), which happens when $f_{n=0}^{\mathrm{QM}}=0$. At the same time, this results in the maximum QM stellar masses in the admixed star becoming smaller (in some cases by a factor of 10). This opens a new stability window of ultra-low QM masses (when surrounded by DM) in the range between ${10}^{-18}$ and ${10}^{-4}$$M_{\odot }$, depending on $m_D$, which correspond to the dark strange planets and strangelets discussed above.

In the same way, Figure 4 shows the results for the coupled radial pulsation Equations (4) and (5) assuming an oscillating DM core with boundary conditions (6) and (7) for different fixed central energy densities of QM. Clearly, the general behavior is qualitatively different, resembling the behavior of nucleonic stars. However, frequencies are very large, reaching ∼$3\times {10}^5$ kHz, in contrast with the few and tens of kHz for hadronic and quark stars, respectively [50]. In almost all cases, the DM core is essentially unaffected by QM due to its very large central density, the exception being the case with $m_D=1$ GeV. On the other hand, non-trivial effects show up in the $f_{n=0}^{\mathrm{DM}}$ vs. M diagrams: the maximum stable masses do not correspond to the the ones in the mass–radius diagram of Figure 2; they are smaller.

Furthermore, in Figure 4 one sees that, for the case of ${ϵ}_c^{\mathrm{QM}}=1.5\mathrm{GeV}/{\mathrm{fm}}^3$ with $m_D=1$ GeV, zero frequencies are not reached in the corresponding DM stars. After some point, the solutions become mechanically unstable in the sense of having negative QM pressure profiles inside the DM star. A similar phenomenon occurs for larger $m_D$, though it is much less visible. Systems exhibiting negative pressures, e.g., dark energy inducing an accelerated expansion of the universe, are not strictly prohibited but should be carefully interpreted. For instance, fluids develop negative-pressure states when stresses are applied for long periods. In the case of hybrid neutron stars with a first-order hadron–quark transition, oscillations to negative-pressure states may accelerate the nucleation of bubbles around the transition region, which, in the limit of large amplitudes, induce mechanical instabilities [51]. In our case of oscillating QM and DM cores with small amplitudes having negative-pressure profiles, the two-fluid TOV equations allows for their existence as hydrostatically-equilibrated configurations that are potentially unstable when disturbed by radial perturbations leading to the automatic collapse of the whole admixed star. In other words, negative-pressure interiors lead immediately to complex oscillation frequencies. Only by increasing ${ϵ}_c^{\mathrm{QM}/\mathrm{DM}}$ do the instabilities disappear and one is able to find only real frequencies.

3.2. Admixtures of Quark Matter and Strongly $(y={10}^3)$ Interacting Dark Matter

3.2.1. Solving the Two-Fluid TOV Equations

Similar to the case of weakly self-interacting DM, we solved the two-fluid TOV Equation (2) with the condition $p_{\mathrm{QM}}\left(R_{\mathrm{QM}}\right)=0$ for different central energy densities of strongly self-interacting DM. We present our results in Figure 5. As in the weak limit for DM, in most of the cases stellar masses, radii and central energy densities of the QM core are not appreciably affected. However, for increasing DM central energy densities, some relevant variations occur. In particular, when $m_D=1,100,200$ GeV, the maximum QM central energy densities are increased by a factor of ∼20. The cases with $m_D=10,50$ GeV show sizable variations of the masses and radii, especially near the maximum mass, where the presence of DM reduces the QM core masses down to ∼$0.4M_{\odot }$. Interestingly, in the case of $m_D=10$ GeV, the central QM energy density is almost unaffected by DM, whereas, for $m_D=50$ GeV, it is dramatically increased. In the case with $m_D=500$ GeV, the QM masses and radii increase by any amount of DM, and the QM central energy density is augmented by a factor of 10. Analogously to what we have seen before, in the cases with $m_D=50$ to 500 GeV, there are plateaus whose widths increase with $m_D$.

In order to study the structure of the opposite case, we solved the two-fluid TOV Equations (2) with the condition $p_{\mathrm{DM}}\left(R_{\mathrm{DM}}\right)=0$ for different central energy densities of QM, as displayed in Figure 6. Pure $y={10}^3$ DM stars display the same qualitative behavior as in the case of $m_D=500$ GeV, since, in this case, their masses and radii are almost unaffected by any amount of QM due to very large DM central energy densities. The same is true for $m_D=200$ GeV. See also Table 1.

Although not noticeable in Figure 6, when $m_D=1$ GeV, our calculations show that the QM in the DM core yields higher masses (not shown in the figure) that are in contradiction with the negative gradient of pressure required by the TOV equations. This happens, because we are considering unstable QM central energy densities for the DM star, as can be seen in Figure 5, which are manifested by producing increasing profiles of pressure then leading to mechanical instabilities associated with complex frequencies, thus destabilizing the whole admixed star. When $m_D=10$ GeV, we find a self-bound-like behavior for the DM star. This occurs since the DM and QM central energy densities are almost equal in the admixed star, and, in this case, the QM component dominates, modifying the behavior in the mass–radius diagram. As before, the central DM densities are increased by a factor of 10. The cases with $m_D=50,100$ GeV display a behavior that is a mixture of dark and quark matter, where QM mainly affects the sector of lower DM stellar masses, and the structure remains almost the same near the maximum mass. There, the central DM energy densities are almost the same since the DM energy densities are enormous compared the QM ones.

3.2.2. Solving the Coupled Radial Oscillation Equations

Finally, we solved the coupled radial pulsation Equations (4) and (5) assuming an oscillating QM core in the admixed star with the boundary conditions (6) and (7) for different fixed central energy densities of strongly interacting DM. The results are displayed in Figure 7. We found that only the case of $m_D=1$ GeV was unaffected by strongly self-interacting DM. As we increased $m_D$, the fundamental frequency was strongly affected. In fact, as occurred in the QM cases with weakly self-interacting DM, only low-mass QM stars survived radial oscillations and behaved as strange quark planets and strangelets. The oscillation frequencies of these objects can reach ∼${10}^5$ kHz for $m_D=500$ GeV.

In Figure 8, we show our results6 after solving the coupled radial pulsation Equations (4) and (5), assuming an oscillating DM core with boundary conditions (6) and (7) for different central energy densities of QM. In correspondence with the results of Figure 6 for the cases $m_D=1,10$ GeV, only a small family of DM stars survived the radial oscillation analysis for low mass stars. These DM stars increase their stability as long as one increases the QM component. The qualitative behavior resembles that of a strange star. This occurs due to the high QM central energy densities compared to the DM ones, with the QM component dominating the stability of the admixed star. On the other hand, the cases with $m_D=50,100,200,500$ GeV display the standard behavior of pure $y={10}^3$ DM stars due to the very high DM central energy densities. Interestingly, the same phenomenon of increasing stability for higher central QM energy densities occurs in all these cases. The physical picture indicates that low QM central energy densities support a small subset of DM stars against gravitational collapse. As we increase the central energy densities, the admixed star supports higher and higher central DM energy densities.

4. Summary and Outlook

We have studied the stability and main global features of strange quark stars admixed with fermionic dark matter in a large parameter space, allowing dark fermion masses from 1 to 500 GeV and considering weakly and strongly self-interacting dark matter. For simplicity, cold quark matter was described within the MIT bag model with $B^{1/4}=145$ MeV which, in the one-fluid case, produces strange quark stars.

After solving the two-fluid TOV equations, we computed the associated stellar structure. We found that, depending on $m_D$ and the interaction parameter y, some of the obtained QM and DM stars display significant modifications of their stellar masses and radii, whereas others show no change at all. Furthermore, some of the DM stars display a self-bound-like behavior in the mass–radius diagram. In most situations, the central QM and DM densities are increased by the presence of the other component in the admixed star.

For the radial pulsation analysis, a full general-relativistic two-fluid Sturm–Liouville-like problem should, in principle, be solved. Instead, inspired by the way one usually solves the two-fluid TOV equations by separating the total pressure and energy density into QM and DM components, we developed a framework where we separate the total Lagrangian variables entering the oscillation equations into QM and DM contributions. This method allowed us to solve the problem, assuming that we disturb only one component and the other is affected only indirectly.

Our calculations indicate that the static stability criterion $\partial M_{\mathrm{QM}/\mathrm{DM}}/\partial {ϵ}_c^{\mathrm{QM}/\mathrm{DM}}\ge 0$ alone might produce misleading and incomplete results when applied to two-fluid stars. We found that, in the case of QM stars admixed with DM, predominantly in the case of $y=0.1$, only very small QM stellar masses are dynamically stable leading to dark strange planets and dark strangelets. On the other hand, DM stars are mainly affected when small values of $m_D$ are considered, since larger dark fermion masses induce ultra-dense cores for which the QM contribution is almost negligible.

Although our results are still very sensitive, both to the dark fermion mass $m_D$ and the kind of self-interaction involved in the dark sector, whose scale is encoded in the dimensionless variable y, there is hope that the parameter space can be dramatically constrained by gravitational wave events, as discussed recently in Ref. [52].

In a future publications, we plan to refine our description using an equation of state obtained from perturbative QCD [17]. It would also be interesting to explore the effects of adding a nuclear mantle to our dark strange planets.