Stability boundary of neutron-dark matter mixed stars

Evaluating the stability of compact stars is important in astrophysics, particularly in the multi-messenger era marked by the gravitational wave observation of neutron stars (NSs) [1]. For single fluid models, the stability criterion has two equivalent formulations: one is based on the sign of fundamental radial oscillation frequencies [2] and the other on the extrema in mass-radius (M-R) curves (known as the BTM criterion [3] proposed by Bardeen, Thorne, and Meltzer). The former is time consuming but also yields non-zero oscillation modes and is applicable to multi-fluid cases. The latter is intuitive and convenient; however, it may not be extended to multi-fluid systems.

Two main scenarios are observed for multi-fluid stellar structures: hybrid stars [4–7], where the interactions among different fluids are strong, and mixed stars [8–17], where only gravitational interactions are present. The hybrid star case is relatively simple, and its stability rule can be found in [7]. As a byproduct in another study [17], other collaborators and one of the current authors conjectured necessary conditions from the BTM criterion for mixed stars, which requires further modifications. In fact, for the mixed star case, two independent approaches for determining stability are known. One is to compute the radial oscillation frequencies [2, 18], as mentioned previously, and the other is to locate critical curves in a certain parameter space [8]. The equivalence of these two criteria was suggested in [8] and applied with specific equations of state (EoSs) in [19]. In this paper, we provide a rigorous proof that the critical curve criterion is equivalent to the existence of zero-frequency radial oscillation modes, and further explore its application to dark matter (DM) – nuclear matter (NM) mixed stars.

NSs, with their extremely high compactness, serve as natural laboratories for testing theoretical physics. As relativistic objects, they should be described within the framework of general relativity. Their macroscopic properties (such as mass and radius) are determined by the Tolman-Oppenheimer-Volkoff (TOV) equations [20, 21], which describe gravitational equilibrium, and by the EoS, which governs microscopic nuclear interactions [22]. In our discussion, these objects are treated at zero temperature, as is standard in this context. Consequently, the pressure is considered as depending solely on the energy density.

Evidence for the existence of DM appears in multiple sources, including galactic rotation curves, large-scale structure formation, and gravitational lensing [15]. Precise cosmological observations and measurements further indicate that DM constitutes approximately 25% of the total energy budget of the universe [23].

At the particle physics level, one of the DM candidates is the weakly interacting massive particle (WIMP), such as the neutralino in supersymmetric theories, the thermal production mechanism of which can naturally explain the observed DM abundance [24, 25]. DM may possess other properties beyond the traditional WIMP paradigm. For example, asymmetric DM models, in which the DM particle is not self-conjugate [26, 27], provide another promising theoretical framework.

Additionally, DM could consist of bosonic particles with self-interactions, which would affect its dynamical evolution within NSs [25]. This so-called self-interacting DM (SIDM) extends the standard cold DM model by incorporating self-interactions, which can effectively alleviate small-scale structure issues [28, 29]. A key constraint on this model requires the ratio of its self-scattering cross-section to the mass to lie within a narrow range [28, 30, 31]. SIDM provides a viable theoretical pathway for forming compact stellar structures [32].

DM detection strategies can be broadly classified into direct and indirect detections. Direct detection experiments (e.g., DAMA, CoGeNT, CDMS, and XENON [33−36]) aim to detect signals from DM-nucleus scattering within terrestrial detectors [15]. By contrast, indirect detection is achieved by studying the effects of DM on astrophysical objects. For instance, the extremely strong gravitational field and high density of NS interiors may make them possible DM trappers. Subsequent physical changes, such as heating, accumulation, or collapse, in NSs can then be used to constrain the DM properties [25]. Specifically, the mechanisms by which DM may be captured into NSs have been studied extensively [37−43]. If the captured DM is self-annihilating, its subsequent annihilation can significantly heat the stellar interior, leading to observable increases in surface temperature or luminosity [38−40]. By contrast, non-annihilating DM would gradually accumulate within the star, potentially giving rise to stable mixed configurations [44], or alternatively, collapsing into a black hole once the accumulated DM reaches its self-gravitating limit [25]. Beyond capture, an alternative scenario is that the collapse leading to NS formation could be triggered within a pre-existing DM halo [45−47]; that is, the mixing of DM and NM may occur during the formation of the compact object itself.

Extensive research has been conducted on the influence of DM on the properties of compact stars (see, e.g., [15, 48−50]). The conditions for the stable existence of mixed stars composed of these components must be discussed; relevant studies include Refs. [17, 19, 46, 51−54]. Our study focuses primarily on mixed stars consisting of NM and SIDM or fermionic DM, and we systematically investigate their stability and macroscopic properties using different methods.

The remainder of this paper is organized as follows. In Section II, we establish the theoretical framework for multi-fluid stellar models, presenting the general-relativistic equilibrium equations and outlining the microphysical inputs, namely EoSs, used in our numerical examples. Section III presents the central theoretical result of this work: a rigorous proof, based on a variational principle, of the formal equivalence between the existence of the zero-frequency radial oscillation mode and static critical curve condition. Section IV provides a numerical validation of our theorem, showing that the stability boundaries computed via direct radial oscillation analysis precisely match the theoretically derived critical curves for various EoS combinations. Finally, Section V summarizes our main results and discusses their astrophysical implications. For completeness, we provide a background review of the static stability criteria and detailed framework of the standard radial pulsation equations in Appendices A and B, respectively.

To model the coexistence of DM and NM in compact astrophysical objects and their corresponding structure, we begin with a basic description of their possible internal composition, followed by the TOV equations for the spherical equilibrium structure within the framework of general relativity. The formulation is completed by specifying the EoS that relates pressure to energy density.

To model the coexistence of DM and NM in compact astrophysical objects and their corresponding structure, we begin with a basic description of their possible internal composition, followed by the TOV equations for the spherical equilibrium structure within the framework of general relativity. The formulation is completed by specifying the EoS that relates pressure to energy density.

To model the coexistence of DM and NM in compact astrophysical objects and their corresponding structure, we begin with a basic description of their possible internal composition, followed by the TOV equations for the spherical equilibrium structure within the framework of general relativity. The formulation is completed by specifying the EoS that relates pressure to energy density.

From the perspective of geometric configuration, two-fluid star models can be preliminarily classified into two categories. Hybrid star: a NM (DM) core enveloped by a DM (NM) crust; mixed star: a mixed core enveloped by a single-component shell [17, 51] ¹, as shown in Fig. 1.

Figure 1.  (color online) Schematic of two different possible structures: hybrid stars (I) and mixed stars (II).

The formation mechanisms of the layered and mixed-core configurations are fundamentally different from a microscopic interaction perspective. The emergence of a layered structure requires the existence of non-gravitational interactions [47] between DM and NM, such that the total energy of the two-component system in a phase-separated state is lower than that in the mixed state, whereas a mixed core arises when only gravitational interactions are present, lacking the driving force for a sharp interface.

This structural difference profoundly influences their macroscopic manifestations. In M-R relations, layered (hybrid) models depend on a single parameter variation, with a clear evolutionary path of a one-dimensional curve. By contrast, mixed-core models exhibit more complex parameter dependence corresponding to the central parameters of both DM and NM, which causes their M-R relation to form a two-dimensional continuous surface topologically.

In terms of stability, both types can have their stable regions rigorously determined via radial perturbation analysis. However, as shown in Fig. 2, unlike that for pure NS, the stability boundary characterization for layered stars exhibits a slight offset from the traditional BTM criterion [3] in the M-R relation, according to the radial mode results calculated in [7]. The stability analysis is more complex for mixed-core models, and their manifestation in macroscopic quantities must satisfy a generalization of the BTM criterion, as discussed in subsequent sections. In the following section, we focus on this mixed star case, which is scenario II in Fig. 1.

Figure 2.  (color online) Schematic of the stability turning point of the hybrid star deviates from the mass vertex.

From the perspective of geometric configuration, two-fluid star models can be preliminarily classified into two categories. Hybrid star: a NM (DM) core enveloped by a DM (NM) crust; mixed star: a mixed core enveloped by a single-component shell [17, 51] ¹, as shown in Fig. 1.

Figure 1.  (color online) Schematic of two different possible structures: hybrid stars (I) and mixed stars (II).

The formation mechanisms of the layered and mixed-core configurations are fundamentally different from a microscopic interaction perspective. The emergence of a layered structure requires the existence of non-gravitational interactions [47] between DM and NM, such that the total energy of the two-component system in a phase-separated state is lower than that in the mixed state, whereas a mixed core arises when only gravitational interactions are present, lacking the driving force for a sharp interface.

This structural difference profoundly influences their macroscopic manifestations. In M-R relations, layered (hybrid) models depend on a single parameter variation, with a clear evolutionary path of a one-dimensional curve. By contrast, mixed-core models exhibit more complex parameter dependence corresponding to the central parameters of both DM and NM, which causes their M-R relation to form a two-dimensional continuous surface topologically.

In terms of stability, both types can have their stable regions rigorously determined via radial perturbation analysis. However, as shown in Fig. 2, unlike that for pure NS, the stability boundary characterization for layered stars exhibits a slight offset from the traditional BTM criterion [3] in the M-R relation, according to the radial mode results calculated in [7]. The stability analysis is more complex for mixed-core models, and their manifestation in macroscopic quantities must satisfy a generalization of the BTM criterion, as discussed in subsequent sections. In the following section, we focus on this mixed star case, which is scenario II in Fig. 1.

Figure 2.  (color online) Schematic of the stability turning point of the hybrid star deviates from the mass vertex.

From the perspective of geometric configuration, two-fluid star models can be preliminarily classified into two categories. Hybrid star: a NM (DM) core enveloped by a DM (NM) crust; mixed star: a mixed core enveloped by a single-component shell [17, 51] ¹, as shown in Fig. 1.

Figure 1.  (color online) Schematic of two different possible structures: hybrid stars (I) and mixed stars (II).

The formation mechanisms of the layered and mixed-core configurations are fundamentally different from a microscopic interaction perspective. The emergence of a layered structure requires the existence of non-gravitational interactions [47] between DM and NM, such that the total energy of the two-component system in a phase-separated state is lower than that in the mixed state, whereas a mixed core arises when only gravitational interactions are present, lacking the driving force for a sharp interface.

This structural difference profoundly influences their macroscopic manifestations. In M-R relations, layered (hybrid) models depend on a single parameter variation, with a clear evolutionary path of a one-dimensional curve. By contrast, mixed-core models exhibit more complex parameter dependence corresponding to the central parameters of both DM and NM, which causes their M-R relation to form a two-dimensional continuous surface topologically.

In terms of stability, both types can have their stable regions rigorously determined via radial perturbation analysis. However, as shown in Fig. 2, unlike that for pure NS, the stability boundary characterization for layered stars exhibits a slight offset from the traditional BTM criterion [3] in the M-R relation, according to the radial mode results calculated in [7]. The stability analysis is more complex for mixed-core models, and their manifestation in macroscopic quantities must satisfy a generalization of the BTM criterion, as discussed in subsequent sections. In the following section, we focus on this mixed star case, which is scenario II in Fig. 1.

Figure 2.  (color online) Schematic of the stability turning point of the hybrid star deviates from the mass vertex.

The commonly used model for describing the compact star equilibrium configuration is provided by the TOV formalism. This framework comprises a set of differential equations for a static, non-rotating star modeled as an ideal fluid in spherical symmetry.

We consider a spherically symmetric spacetime with Schwarzschild-like coordinates $ (t, r, \theta, \phi) $. The line element is

where ${\rm d}\Omega ^2={\rm d}\theta ^2+\sin ^2\theta {\rm d}\phi ^2$, and the metric functions satisfy ${\rm e}^{2\nu(r)} = \left(1 - \dfrac{2G m(r)}{r c^{2}}\right)^{-1}$. The energy-momentum tensor for an ideal fluid is

where p is the fluid pressure, ζ is the specific energy density, ρ is the rest mass density, and $ u^{\mu} $ is the four-velocity vector. The TOV equations can be derived by substituting them into the Einstein equation

together with the energy-momentum conservation

The explicit TOV equations are expressed as follows:

where $ m\left(r \right) $ is the gravitational mass within radius r.

By defining the total energy density as $ \varepsilon = \rho c^2 + \rho\zeta $ and adopting the natural unit system, we simplify (5) to

We adopt the equations in (6) in our subsequent single-fluid numerical calculations.

The boundary conditions are prescribed as follows: at the star center, the mass within an infinitesimal radius must approach zero; at the surface (where $ r = R $, with $ R $ denoting the star radius), the metric potential $\phi$ assumes its vacuum form, $\phi = \dfrac{1}{2} \ln \left(1 - 2GM / (R c^2)\right)$, which matches the Schwarzschild exterior solution, where $ M \equiv m(R) $ represents the total mass enclosed within the star radius. By contrast, the central pressure $ p_c \equiv p(0) $ serves as a free parameter that can be specified independently.

Although the aforementioned TOV equations typically model the interior of a compact star as an ideal fluid composed of a single species of matter, the physical scenario under consideration here is more complex: a DM/NM mixed star. Generally, models with multiple matter components can be constructed using either a single-fluid approximation [55, 56] (where the system is still treated as one effective fluid, but with an EoS incorporating contributions from multiple components) or a multi-fluid formalism (where distinct matter components are described by separate EoSs). In the case of a weakly coupled multi-fluid model with only gravitational interactions—which means that the energy-momentum tensors of the individual components are decoupled—a multi-fluid generalization of the TOV equations can be derived [11, 57]:

where the subscript I indicates the components of the substance. $ m_{t}=\displaystyle\sum\nolimits_I{m_I} $, $ p_{t}=\displaystyle\sum\nolimits_I{p_I} $ , and $ \varepsilon _{t}=\displaystyle\sum\nolimits_I{\varepsilon_I} $ are the sums of the respective quantities from the contributions of each fluid component. When $ p_I $ reaches zero, the corresponding radius and mass for that component are $ R_I $ and $ M_I $. When $ p_{t} $ reaches zero, we obtain $ M_I $, the total star radius $ R_t $, and total star mass $ M=\displaystyle\sum\nolimits_I{M_I} $.

In addition to the nature units, we further apply astrophysical units for later convenience. The units of mass, pressure, energy density, and radius are $ M_{\odot} $, $ p_{\odot}= \varepsilon_{\odot} = 1/ M_{\odot} ^2 $, and $ {\rm{km}} $, respectively, where $ M_{\odot} $ denotes the solar mass.

The commonly used model for describing the compact star equilibrium configuration is provided by the TOV formalism. This framework comprises a set of differential equations for a static, non-rotating star modeled as an ideal fluid in spherical symmetry.

We consider a spherically symmetric spacetime with Schwarzschild-like coordinates $ (t, r, \theta, \phi) $. The line element is

where ${\rm d}\Omega ^2={\rm d}\theta ^2+\sin ^2\theta {\rm d}\phi ^2$, and the metric functions satisfy ${\rm e}^{2\nu(r)} = \left(1 - \dfrac{2G m(r)}{r c^{2}}\right)^{-1}$. The energy-momentum tensor for an ideal fluid is

where p is the fluid pressure, ζ is the specific energy density, ρ is the rest mass density, and $ u^{\mu} $ is the four-velocity vector. The TOV equations can be derived by substituting them into the Einstein equation

together with the energy-momentum conservation

The explicit TOV equations are expressed as follows:

where $ m\left(r \right) $ is the gravitational mass within radius r.

By defining the total energy density as $ \varepsilon = \rho c^2 + \rho\zeta $ and adopting the natural unit system, we simplify (5) to

We adopt the equations in (6) in our subsequent single-fluid numerical calculations.

The boundary conditions are prescribed as follows: at the star center, the mass within an infinitesimal radius must approach zero; at the surface (where $ r = R $, with $ R $ denoting the star radius), the metric potential $\phi$ assumes its vacuum form, $\phi = \dfrac{1}{2} \ln \left(1 - 2GM / (R c^2)\right)$, which matches the Schwarzschild exterior solution, where $ M \equiv m(R) $ represents the total mass enclosed within the star radius. By contrast, the central pressure $ p_c \equiv p(0) $ serves as a free parameter that can be specified independently.

Although the aforementioned TOV equations typically model the interior of a compact star as an ideal fluid composed of a single species of matter, the physical scenario under consideration here is more complex: a DM/NM mixed star. Generally, models with multiple matter components can be constructed using either a single-fluid approximation [55, 56] (where the system is still treated as one effective fluid, but with an EoS incorporating contributions from multiple components) or a multi-fluid formalism (where distinct matter components are described by separate EoSs). In the case of a weakly coupled multi-fluid model with only gravitational interactions—which means that the energy-momentum tensors of the individual components are decoupled—a multi-fluid generalization of the TOV equations can be derived [11, 57]:

where the subscript I indicates the components of the substance. $ m_{t}=\displaystyle\sum\nolimits_I{m_I} $, $ p_{t}=\displaystyle\sum\nolimits_I{p_I} $ , and $ \varepsilon _{t}=\displaystyle\sum\nolimits_I{\varepsilon_I} $ are the sums of the respective quantities from the contributions of each fluid component. When $ p_I $ reaches zero, the corresponding radius and mass for that component are $ R_I $ and $ M_I $. When $ p_{t} $ reaches zero, we obtain $ M_I $, the total star radius $ R_t $, and total star mass $ M=\displaystyle\sum\nolimits_I{M_I} $.

In addition to the nature units, we further apply astrophysical units for later convenience. The units of mass, pressure, energy density, and radius are $ M_{\odot} $, $ p_{\odot}= \varepsilon_{\odot} = 1/ M_{\odot} ^2 $, and $ {\rm{km}} $, respectively, where $ M_{\odot} $ denotes the solar mass.

The commonly used model for describing the compact star equilibrium configuration is provided by the TOV formalism. This framework comprises a set of differential equations for a static, non-rotating star modeled as an ideal fluid in spherical symmetry.

We consider a spherically symmetric spacetime with Schwarzschild-like coordinates $ (t, r, \theta, \phi) $. The line element is

where ${\rm d}\Omega ^2={\rm d}\theta ^2+\sin ^2\theta {\rm d}\phi ^2$, and the metric functions satisfy ${\rm e}^{2\nu(r)} = \left(1 - \dfrac{2G m(r)}{r c^{2}}\right)^{-1}$. The energy-momentum tensor for an ideal fluid is

where p is the fluid pressure, ζ is the specific energy density, ρ is the rest mass density, and $ u^{\mu} $ is the four-velocity vector. The TOV equations can be derived by substituting them into the Einstein equation

together with the energy-momentum conservation

The explicit TOV equations are expressed as follows:

where $ m\left(r \right) $ is the gravitational mass within radius r.

By defining the total energy density as $ \varepsilon = \rho c^2 + \rho\zeta $ and adopting the natural unit system, we simplify (5) to

We adopt the equations in (6) in our subsequent single-fluid numerical calculations.

The boundary conditions are prescribed as follows: at the star center, the mass within an infinitesimal radius must approach zero; at the surface (where $ r = R $, with $ R $ denoting the star radius), the metric potential $\phi$ assumes its vacuum form, $\phi = \dfrac{1}{2} \ln \left(1 - 2GM / (R c^2)\right)$, which matches the Schwarzschild exterior solution, where $ M \equiv m(R) $ represents the total mass enclosed within the star radius. By contrast, the central pressure $ p_c \equiv p(0) $ serves as a free parameter that can be specified independently.

Although the aforementioned TOV equations typically model the interior of a compact star as an ideal fluid composed of a single species of matter, the physical scenario under consideration here is more complex: a DM/NM mixed star. Generally, models with multiple matter components can be constructed using either a single-fluid approximation [55, 56] (where the system is still treated as one effective fluid, but with an EoS incorporating contributions from multiple components) or a multi-fluid formalism (where distinct matter components are described by separate EoSs). In the case of a weakly coupled multi-fluid model with only gravitational interactions—which means that the energy-momentum tensors of the individual components are decoupled—a multi-fluid generalization of the TOV equations can be derived [11, 57]:

where the subscript I indicates the components of the substance. $ m_{t}=\displaystyle\sum\nolimits_I{m_I} $, $ p_{t}=\displaystyle\sum\nolimits_I{p_I} $ , and $ \varepsilon _{t}=\displaystyle\sum\nolimits_I{\varepsilon_I} $ are the sums of the respective quantities from the contributions of each fluid component. When $ p_I $ reaches zero, the corresponding radius and mass for that component are $ R_I $ and $ M_I $. When $ p_{t} $ reaches zero, we obtain $ M_I $, the total star radius $ R_t $, and total star mass $ M=\displaystyle\sum\nolimits_I{M_I} $.

In addition to the nature units, we further apply astrophysical units for later convenience. The units of mass, pressure, energy density, and radius are $ M_{\odot} $, $ p_{\odot}= \varepsilon_{\odot} = 1/ M_{\odot} ^2 $, and $ {\rm{km}} $, respectively, where $ M_{\odot} $ denotes the solar mass.

Many EoS models exist for NS and quark stars, such as FPS [58], APR (ALF) [59], WFF [60], MPA1 [61], MSL [62], H4, [63], and the BSK family of models [64].

To verify the stability, we can use various types of numerical or analytical EoS. In our study, we use one numerical EoS and three analytical EoSs as examples.

For the NS EoS, we first consider the famous SLy4 [65−67], which survives from the GW170817 [1] and NICER [68] observation constraints. It is a Skyrme-type effective nuclear interaction parametrization that is tailored for describing asymmetric NM. Its parameters are derived by fitting to microscopic NM calculations and experimental data of doubly magic nuclei, thereby achieving self-consistency between nuclear and NS physics. However, it has many adjustable parameters and does not have an analytic form.

The second EoS is a double-polytropic model [69] motivated by AdS/QCD and the Witten-Sakai-Sugimoto framework [70, 71], in which an instanton-gas structure is introduced in the D-brane background and the corresponding free-energy variational equations are solved. This analytic EoS should be regarded as a simplified reference model, rather than a fully calibrated microphysical description. For the particular parameter choice used in this work, it supports a maximum mass of approximately $ 1.85M_{\odot} $, and therefore is not intended to represent a complete phenomenological fit to the most restrictive massive-pulsar constraints. Nevertheless, the double-polytropic form is flexible: with suitable choices of parameters, it can be made compatible with current observational constraints and the causality condition. In the present study, we use this model mainly because its closed analytic form makes the numerical comparison transparent and easily reproducible. The focus of this work is the proof and numerical verification of the equivalence between the dynamical and static stability criteria. The EoS for the double-polytropic model is

where $ {\cal{A}} =1.8\times 10^{-5}\times \ell ^{-7} $, and $ \ell $ (brane asymptotic separation) is the only adjustable parameter. We refer to this EoS as the Holographic EoS later.

For the DM component in mixed stars, the EoS is determined by the selected DM model. Studies have employed free Fermi gas DM combined with the SLy4 NS model to investigate mixed star stability [19]. In this work, we include this EoS for comparison. The EoS for such free Fermionic DM is expressed as [72, 73]

where $ m_{f} $ is the fermion mass and $ k_F $ is the Fermi momentum. Various interactions can be incorporated into this free Fermi gas model to construct more sophisticated EoSs.

For Bosonic DM, a stable mixed star model cannot be formed with free bosons owing to the absence of degeneracy pressure. Given the weak interaction between DM and NM, it is natural to consider self-interacting Bosonic DM to construct the DM EoS that balances the self-gravitational effect. In this study, we utilize the scalar field theory with the potential $ V\left(\phi \right) =\dfrac{{m_b}^2}{2}\left| \phi \right|^2+\dfrac{\lambda _4}{4}\left| \phi \right|^4 $, as proposed in [74], where $ m_b $ is the boson mass and $ \lambda_4 $ is the coupling. In the isotropic limit, this theory yields the following EoS [32]:

where the free parameter is $ B=\dfrac{0.08}{\sqrt{\lambda _4}}\left(\dfrac{m}{{\rm{GeV}}} \right) ^2 $.

Many EoS models exist for NS and quark stars, such as FPS [58], APR (ALF) [59], WFF [60], MPA1 [61], MSL [62], H4, [63], and the BSK family of models [64].

To verify the stability, we can use various types of numerical or analytical EoS. In our study, we use one numerical EoS and three analytical EoSs as examples.

For the NS EoS, we first consider the famous SLy4 [65−67], which survives from the GW170817 [1] and NICER [68] observation constraints. It is a Skyrme-type effective nuclear interaction parametrization that is tailored for describing asymmetric NM. Its parameters are derived by fitting to microscopic NM calculations and experimental data of doubly magic nuclei, thereby achieving self-consistency between nuclear and NS physics. However, it has many adjustable parameters and does not have an analytic form.

The second EoS is a double-polytropic model [69] motivated by AdS/QCD and the Witten-Sakai-Sugimoto framework [70, 71], in which an instanton-gas structure is introduced in the D-brane background and the corresponding free-energy variational equations are solved. This analytic EoS should be regarded as a simplified reference model, rather than a fully calibrated microphysical description. For the particular parameter choice used in this work, it supports a maximum mass of approximately $ 1.85M_{\odot} $, and therefore is not intended to represent a complete phenomenological fit to the most restrictive massive-pulsar constraints. Nevertheless, the double-polytropic form is flexible: with suitable choices of parameters, it can be made compatible with current observational constraints and the causality condition. In the present study, we use this model mainly because its closed analytic form makes the numerical comparison transparent and easily reproducible. The focus of this work is the proof and numerical verification of the equivalence between the dynamical and static stability criteria. The EoS for the double-polytropic model is

where $ {\cal{A}} =1.8\times 10^{-5}\times \ell ^{-7} $, and $ \ell $ (brane asymptotic separation) is the only adjustable parameter. We refer to this EoS as the Holographic EoS later.

For the DM component in mixed stars, the EoS is determined by the selected DM model. Studies have employed free Fermi gas DM combined with the SLy4 NS model to investigate mixed star stability [19]. In this work, we include this EoS for comparison. The EoS for such free Fermionic DM is expressed as [72, 73]

where $ m_{f} $ is the fermion mass and $ k_F $ is the Fermi momentum. Various interactions can be incorporated into this free Fermi gas model to construct more sophisticated EoSs.

For Bosonic DM, a stable mixed star model cannot be formed with free bosons owing to the absence of degeneracy pressure. Given the weak interaction between DM and NM, it is natural to consider self-interacting Bosonic DM to construct the DM EoS that balances the self-gravitational effect. In this study, we utilize the scalar field theory with the potential $ V\left(\phi \right) =\dfrac{{m_b}^2}{2}\left| \phi \right|^2+\dfrac{\lambda _4}{4}\left| \phi \right|^4 $, as proposed in [74], where $ m_b $ is the boson mass and $ \lambda_4 $ is the coupling. In the isotropic limit, this theory yields the following EoS [32]:

where the free parameter is $ B=\dfrac{0.08}{\sqrt{\lambda _4}}\left(\dfrac{m}{{\rm{GeV}}} \right) ^2 $.

Many EoS models exist for NS and quark stars, such as FPS [58], APR (ALF) [59], WFF [60], MPA1 [61], MSL [62], H4, [63], and the BSK family of models [64].

To verify the stability, we can use various types of numerical or analytical EoS. In our study, we use one numerical EoS and three analytical EoSs as examples.

For the NS EoS, we first consider the famous SLy4 [65−67], which survives from the GW170817 [1] and NICER [68] observation constraints. It is a Skyrme-type effective nuclear interaction parametrization that is tailored for describing asymmetric NM. Its parameters are derived by fitting to microscopic NM calculations and experimental data of doubly magic nuclei, thereby achieving self-consistency between nuclear and NS physics. However, it has many adjustable parameters and does not have an analytic form.

The second EoS is a double-polytropic model [69] motivated by AdS/QCD and the Witten-Sakai-Sugimoto framework [70, 71], in which an instanton-gas structure is introduced in the D-brane background and the corresponding free-energy variational equations are solved. This analytic EoS should be regarded as a simplified reference model, rather than a fully calibrated microphysical description. For the particular parameter choice used in this work, it supports a maximum mass of approximately $ 1.85M_{\odot} $, and therefore is not intended to represent a complete phenomenological fit to the most restrictive massive-pulsar constraints. Nevertheless, the double-polytropic form is flexible: with suitable choices of parameters, it can be made compatible with current observational constraints and the causality condition. In the present study, we use this model mainly because its closed analytic form makes the numerical comparison transparent and easily reproducible. The focus of this work is the proof and numerical verification of the equivalence between the dynamical and static stability criteria. The EoS for the double-polytropic model is

where $ {\cal{A}} =1.8\times 10^{-5}\times \ell ^{-7} $, and $ \ell $ (brane asymptotic separation) is the only adjustable parameter. We refer to this EoS as the Holographic EoS later.

For the DM component in mixed stars, the EoS is determined by the selected DM model. Studies have employed free Fermi gas DM combined with the SLy4 NS model to investigate mixed star stability [19]. In this work, we include this EoS for comparison. The EoS for such free Fermionic DM is expressed as [72, 73]

where $ m_{f} $ is the fermion mass and $ k_F $ is the Fermi momentum. Various interactions can be incorporated into this free Fermi gas model to construct more sophisticated EoSs.

For Bosonic DM, a stable mixed star model cannot be formed with free bosons owing to the absence of degeneracy pressure. Given the weak interaction between DM and NM, it is natural to consider self-interacting Bosonic DM to construct the DM EoS that balances the self-gravitational effect. In this study, we utilize the scalar field theory with the potential $ V\left(\phi \right) =\dfrac{{m_b}^2}{2}\left| \phi \right|^2+\dfrac{\lambda _4}{4}\left| \phi \right|^4 $, as proposed in [74], where $ m_b $ is the boson mass and $ \lambda_4 $ is the coupling. In the isotropic limit, this theory yields the following EoS [32]:

where the free parameter is $ B=\dfrac{0.08}{\sqrt{\lambda _4}}\left(\dfrac{m}{{\rm{GeV}}} \right) ^2 $.

An equilibrium configuration of a compact star determined by solving the TOV equations can be either stable or unstable against perturbations. The stability can be assessed through two main approaches: a dynamical method, which involves solving the equations for radial oscillations [2], and a static method, which analyzes the properties of sequences of equilibrium solutions.

Although the dynamical method provides a definitive answer, it is often computationally intensive. For single-fluid stars, an efficient static alternative exists in the form of the well-known BTM criterion [3]. For multi-fluid stars, where the BTM criterion fails, a generalization known as the "critical curve criterion" has been developed by Henriques et al. [8] and utilized in subsequent studies on DM admixed stars [19]. A detailed description of these static methods is provided in Appendix A. For completeness and formula consistency, the standard radial pulsation equations (the dynamical method), which was first derived by Kain in [18], is outlined in Appendix B.

Although the "critical curve criterion" has been numerically verified in several contexts, such as those in [16, 19], a formal proof founded upon first principles, such as variational methods and covariant conservation laws, has been lacking. In this section, we present the core contribution of our work, namely a rigorous demonstration of the formal equivalence between the existence of the zero-frequency oscillation mode and static critical curve criterion, and explain the connection with the dynamical stability boundary.

An equilibrium configuration of a compact star determined by solving the TOV equations can be either stable or unstable against perturbations. The stability can be assessed through two main approaches: a dynamical method, which involves solving the equations for radial oscillations [2], and a static method, which analyzes the properties of sequences of equilibrium solutions.

Although the dynamical method provides a definitive answer, it is often computationally intensive. For single-fluid stars, an efficient static alternative exists in the form of the well-known BTM criterion [3]. For multi-fluid stars, where the BTM criterion fails, a generalization known as the "critical curve criterion" has been developed by Henriques et al. [8] and utilized in subsequent studies on DM admixed stars [19]. A detailed description of these static methods is provided in Appendix A. For completeness and formula consistency, the standard radial pulsation equations (the dynamical method), which was first derived by Kain in [18], is outlined in Appendix B.

Although the "critical curve criterion" has been numerically verified in several contexts, such as those in [16, 19], a formal proof founded upon first principles, such as variational methods and covariant conservation laws, has been lacking. In this section, we present the core contribution of our work, namely a rigorous demonstration of the formal equivalence between the existence of the zero-frequency oscillation mode and static critical curve criterion, and explain the connection with the dynamical stability boundary.

An equilibrium configuration of a compact star determined by solving the TOV equations can be either stable or unstable against perturbations. The stability can be assessed through two main approaches: a dynamical method, which involves solving the equations for radial oscillations [2], and a static method, which analyzes the properties of sequences of equilibrium solutions.

Although the dynamical method provides a definitive answer, it is often computationally intensive. For single-fluid stars, an efficient static alternative exists in the form of the well-known BTM criterion [3]. For multi-fluid stars, where the BTM criterion fails, a generalization known as the "critical curve criterion" has been developed by Henriques et al. [8] and utilized in subsequent studies on DM admixed stars [19]. A detailed description of these static methods is provided in Appendix A. For completeness and formula consistency, the standard radial pulsation equations (the dynamical method), which was first derived by Kain in [18], is outlined in Appendix B.

Although the "critical curve criterion" has been numerically verified in several contexts, such as those in [16, 19], a formal proof founded upon first principles, such as variational methods and covariant conservation laws, has been lacking. In this section, we present the core contribution of our work, namely a rigorous demonstration of the formal equivalence between the existence of the zero-frequency oscillation mode and static critical curve criterion, and explain the connection with the dynamical stability boundary.

Preliminary Remarks

Before proceeding to the equivalence proof, several important considerations must be emphasized:

1. Stability vs. zero-mode condition: The fundamental mode vanishing ($ \omega_0^2=0 $) represents a stronger condition than the existence of any zero-mode ($ \exists \,n\in \mathbb{N} \text{ s.t. } \omega _{n}^{2}=0 $). The original critical curve method strictly identifies configurations in which some mode has zero frequency. Additional verification (which we discuss in the final part of the proof section) is required to determine the curve that corresponds exactly to the fundamental zero mode.

2. Mathematical foundations: The existence of fundamental modes, the fact that the eigenvalues are real, and the completeness of mode expansions were rigorously established in [75]. We adopt these conclusions without further proof.

3. Model assumptions and mode continuity: Our proof of the equivalence between the zero-frequency mode and critical curve criterion relies on several physical assumptions that should be made explicit. First, we consider stars composed of ideal fluids, for which the energy-momentum tensor takes the perfect-fluid form and no viscous or heat-conduction effects are present. Second, we assume that the fluids are non-dissipative, so that the dynamical equations are time-reversal invariant and no mode damping or growth arises from entropy production. Third, we consider smooth EoSs, ensuring that the relevant thermodynamic derivatives are well-defined and that the pulsation equations are regular throughout the star. Fourth, we fix the chemical composition and specific entropy of each fluid, so that the only dynamical variables are the fluid displacements and metric perturbations; in particular, no additional degrees of freedom (such as composition gradients or thermal diffusion) enter the problem. Under these assumptions, the dynamical equations are time-reversal invariant. Consequently, the eigenvalues $ \omega_n^2 $ are real and vary continuously as functions of the central pressures in the two-dimensional parameter space [75, 76]. This continuity guarantees that the critical curve, defined by the locus at which the fundamental mode frequency vanishes, properly separates stable and unstable regions, and that the stability boundary cannot be crossed without passing through a zero-frequency mode.

The equivalence proof proceeds by establishing that the critical curve condition is equivalent to the existence of zero-frequency modes. We begin by reformulating the problem and then proceed through two main theorems.

In practical analysis, we identify the zero-mode locus as the set of points in the two-dimensional parameter space spanned by the central pressures of the two fluids,

where the gradient operator is defined as

The quantities M, $ N_{{\rm{NM}}} $, and $ N_{{\rm{DM}}} $ are scalar functions on this parameter space. Then, the equivalence that we aim to prove can be formulated as

that is, the existence of a zero mode is equivalent to any two of the gradients $ \nabla_{{\boldsymbol{\sigma}}}M $, $ \nabla_{{\boldsymbol{\sigma}}}N_{{\rm{NM}}} $, and $ \nabla_{{\boldsymbol{\sigma}}}N_{{\rm{DM}}} $ becoming parallel.

Preliminary Remarks

Before proceeding to the equivalence proof, several important considerations must be emphasized:

1. Stability vs. zero-mode condition: The fundamental mode vanishing ($ \omega_0^2=0 $) represents a stronger condition than the existence of any zero-mode ($ \exists \,n\in \mathbb{N} \text{ s.t. } \omega _{n}^{2}=0 $). The original critical curve method strictly identifies configurations in which some mode has zero frequency. Additional verification (which we discuss in the final part of the proof section) is required to determine the curve that corresponds exactly to the fundamental zero mode.

2. Mathematical foundations: The existence of fundamental modes, the fact that the eigenvalues are real, and the completeness of mode expansions were rigorously established in [75]. We adopt these conclusions without further proof.

3. Model assumptions and mode continuity: Our proof of the equivalence between the zero-frequency mode and critical curve criterion relies on several physical assumptions that should be made explicit. First, we consider stars composed of ideal fluids, for which the energy-momentum tensor takes the perfect-fluid form and no viscous or heat-conduction effects are present. Second, we assume that the fluids are non-dissipative, so that the dynamical equations are time-reversal invariant and no mode damping or growth arises from entropy production. Third, we consider smooth EoSs, ensuring that the relevant thermodynamic derivatives are well-defined and that the pulsation equations are regular throughout the star. Fourth, we fix the chemical composition and specific entropy of each fluid, so that the only dynamical variables are the fluid displacements and metric perturbations; in particular, no additional degrees of freedom (such as composition gradients or thermal diffusion) enter the problem. Under these assumptions, the dynamical equations are time-reversal invariant. Consequently, the eigenvalues $ \omega_n^2 $ are real and vary continuously as functions of the central pressures in the two-dimensional parameter space [75, 76]. This continuity guarantees that the critical curve, defined by the locus at which the fundamental mode frequency vanishes, properly separates stable and unstable regions, and that the stability boundary cannot be crossed without passing through a zero-frequency mode.

The equivalence proof proceeds by establishing that the critical curve condition is equivalent to the existence of zero-frequency modes. We begin by reformulating the problem and then proceed through two main theorems.

In practical analysis, we identify the zero-mode locus as the set of points in the two-dimensional parameter space spanned by the central pressures of the two fluids,

where the gradient operator is defined as

The quantities M, $ N_{{\rm{NM}}} $, and $ N_{{\rm{DM}}} $ are scalar functions on this parameter space. Then, the equivalence that we aim to prove can be formulated as

that is, the existence of a zero mode is equivalent to any two of the gradients $ \nabla_{{\boldsymbol{\sigma}}}M $, $ \nabla_{{\boldsymbol{\sigma}}}N_{{\rm{NM}}} $, and $ \nabla_{{\boldsymbol{\sigma}}}N_{{\rm{DM}}} $ becoming parallel.

Preliminary Remarks

Before proceeding to the equivalence proof, several important considerations must be emphasized:

1. Stability vs. zero-mode condition: The fundamental mode vanishing ($ \omega_0^2=0 $) represents a stronger condition than the existence of any zero-mode ($ \exists \,n\in \mathbb{N} \text{ s.t. } \omega _{n}^{2}=0 $). The original critical curve method strictly identifies configurations in which some mode has zero frequency. Additional verification (which we discuss in the final part of the proof section) is required to determine the curve that corresponds exactly to the fundamental zero mode.

2. Mathematical foundations: The existence of fundamental modes, the fact that the eigenvalues are real, and the completeness of mode expansions were rigorously established in [75]. We adopt these conclusions without further proof.

3. Model assumptions and mode continuity: Our proof of the equivalence between the zero-frequency mode and critical curve criterion relies on several physical assumptions that should be made explicit. First, we consider stars composed of ideal fluids, for which the energy-momentum tensor takes the perfect-fluid form and no viscous or heat-conduction effects are present. Second, we assume that the fluids are non-dissipative, so that the dynamical equations are time-reversal invariant and no mode damping or growth arises from entropy production. Third, we consider smooth EoSs, ensuring that the relevant thermodynamic derivatives are well-defined and that the pulsation equations are regular throughout the star. Fourth, we fix the chemical composition and specific entropy of each fluid, so that the only dynamical variables are the fluid displacements and metric perturbations; in particular, no additional degrees of freedom (such as composition gradients or thermal diffusion) enter the problem. Under these assumptions, the dynamical equations are time-reversal invariant. Consequently, the eigenvalues $ \omega_n^2 $ are real and vary continuously as functions of the central pressures in the two-dimensional parameter space [75, 76]. This continuity guarantees that the critical curve, defined by the locus at which the fundamental mode frequency vanishes, properly separates stable and unstable regions, and that the stability boundary cannot be crossed without passing through a zero-frequency mode.

The equivalence proof proceeds by establishing that the critical curve condition is equivalent to the existence of zero-frequency modes. We begin by reformulating the problem and then proceed through two main theorems.

In practical analysis, we identify the zero-mode locus as the set of points in the two-dimensional parameter space spanned by the central pressures of the two fluids,

where the gradient operator is defined as

The quantities M, $ N_{{\rm{NM}}} $, and $ N_{{\rm{DM}}} $ are scalar functions on this parameter space. Then, the equivalence that we aim to prove can be formulated as

that is, the existence of a zero mode is equivalent to any two of the gradients $ \nabla_{{\boldsymbol{\sigma}}}M $, $ \nabla_{{\boldsymbol{\sigma}}}N_{{\rm{NM}}} $, and $ \nabla_{{\boldsymbol{\sigma}}}N_{{\rm{DM}}} $ becoming parallel.

The first theorem establishes a rigorous connection between the hydrostatic equilibrium and variational properties of the system, forming the foundation for Theorem III.2 and for simplifying the condition to the actual critical curve method. A similar proof of the variational principle for two-fluid systems was previously outlined in [77]. In the present work, we provide a self-contained and explicit derivation, introducing a unified auxiliary function $F(r) $ that consolidates the nonlocal integral contributions, and defining a common proportional function $\Gamma(r)=\lambda_I n_I/(p_I+\epsilon_I) $ to eliminate the Lagrange multipliers efficiently. This construction leads directly to the standard two-fluid TOV equations. Moreover, we establish the necessity and sufficiency of the relation between the variational extremum condition and the hydrostatic equilibrium, providing a complete and transparent foundation for the subsequent stability analysis.

Theorem III.1. A two-fluid isentropic mixed star configuration satisfies the coupled general-relativistic hydrostatic equilibrium equations (TOV equations) if and only if the total mass M takes its extremal value with respect to all time-independent variations that conserve the particle numbers of both fluids.

Proof. We provide the rigorous proof for Theorem III.1. The star configuration satisfies the equilibrium equations if and only if the quantity $ M $ defined by

is extremal under all variations of $ \varepsilon_I(r) $ that maintain the quantities

constant for $ I=1,2 $. This is equivalent to the existence of Lagrange multipliers $ \lambda_1 $ and $ \lambda_2 $ such that $ \delta H = \delta M - \lambda_1 \delta N_1 - \lambda_2 \delta N_2 $ is zero for all variations.

For given variations $ \delta \varepsilon_I(r) $, we obtain

We expand each term as follows (because the integrand vanishes above $ R_I+\delta R_I $, we write the upper limit of integration as positive infinity):

where $ m_t(r) = m_1(r) + m_2(r) $. Thus,

Using the isentropic conditions

and the mass variation

we perform substitution by steps. First, we substitute $ \delta n_I $:

Next, we substitute $ \delta m_J(r) $ into the final term:

By defining

and exchanging the order of integration in the double integral, the final term becomes

Combining all terms in $ \delta H $ yields

For $ \delta H $ to be zero for all variations $ \left[\left\{ \delta \varepsilon _I\left(r \right) \right\} \right] $, the integrand must vanish for each fluid component $ I $:

where $ A(r) = 1 - \dfrac{2 m_t(r)}{r} $. This is the extremization condition for the two-fluid system.

We define

Then, the condition becomes

From these two equations, we obtain a key constraint that reflects the effect of Lagrange multipliers in this two-fluid system:

which implies that

This reflects that the chemical potential between the two fluids satisfies a proportional relationship, which is somewhat different from the case of a single fluid [76].

We define a common function:

We differentiate the condition with respect to r:

where

We take the derivative in (32) and substitute the isentropic condition $ n_I' = \dfrac{n_I \varepsilon_I'}{p_I + \varepsilon_I} $:

We multiply by $ A^{3/2} $ and $ (p_I + \varepsilon_I)^2 $, and then rearrange to solve for $ p_I' $:

From the constraint definition $ \Gamma(r) = \lambda_I \dfrac{n_I}{p_I + \varepsilon_I} $, we obtain

Thus,

Substituting back into (35) yields

Now, we compute $ A' $: using $ m_t' = 4\pi r^2 \varepsilon_t $, we obtain

We substitute into the expression for $ p_I' $:

This is exactly the two-fluid equilibrium equation:

We have explicitly shown that the condition $ \delta H = 0 $ leads directly to the TOV equations. Conversely, since every step in the derivation above (substitution, differentiation, algebraic rearrangement) is reversible, if one starts from the TOV equations and chooses the constants $ \lambda_1,\lambda_2 $ such that $ \lambda_1/\lambda_2 = \mu_1(r)/\mu_2(r) $, one can integrate back to show that $ \delta H=0 $ (here $ \mu_I(r) \equiv (p_I+\varepsilon_I)/n_I $ is the local chemical potential of zero temperature fluid I). This completes the proof of necessity and sufficiency. □

The first theorem establishes a rigorous connection between the hydrostatic equilibrium and variational properties of the system, forming the foundation for Theorem III.2 and for simplifying the condition to the actual critical curve method. A similar proof of the variational principle for two-fluid systems was previously outlined in [77]. In the present work, we provide a self-contained and explicit derivation, introducing a unified auxiliary function $F(r) $ that consolidates the nonlocal integral contributions, and defining a common proportional function $\Gamma(r)=\lambda_I n_I/(p_I+\epsilon_I) $ to eliminate the Lagrange multipliers efficiently. This construction leads directly to the standard two-fluid TOV equations. Moreover, we establish the necessity and sufficiency of the relation between the variational extremum condition and the hydrostatic equilibrium, providing a complete and transparent foundation for the subsequent stability analysis.

Theorem III.1. A two-fluid isentropic mixed star configuration satisfies the coupled general-relativistic hydrostatic equilibrium equations (TOV equations) if and only if the total mass M takes its extremal value with respect to all time-independent variations that conserve the particle numbers of both fluids.

Proof. We provide the rigorous proof for Theorem III.1. The star configuration satisfies the equilibrium equations if and only if the quantity $ M $ defined by

is extremal under all variations of $ \varepsilon_I(r) $ that maintain the quantities

constant for $ I=1,2 $. This is equivalent to the existence of Lagrange multipliers $ \lambda_1 $ and $ \lambda_2 $ such that $ \delta H = \delta M - \lambda_1 \delta N_1 - \lambda_2 \delta N_2 $ is zero for all variations.

For given variations $ \delta \varepsilon_I(r) $, we obtain

We expand each term as follows (because the integrand vanishes above $ R_I+\delta R_I $, we write the upper limit of integration as positive infinity):

where $ m_t(r) = m_1(r) + m_2(r) $. Thus,

Using the isentropic conditions

and the mass variation

we perform substitution by steps. First, we substitute $ \delta n_I $:

Next, we substitute $ \delta m_J(r) $ into the final term:

By defining

and exchanging the order of integration in the double integral, the final term becomes

Combining all terms in $ \delta H $ yields

For $ \delta H $ to be zero for all variations $ \left[\left\{ \delta \varepsilon _I\left(r \right) \right\} \right] $, the integrand must vanish for each fluid component $ I $:

where $ A(r) = 1 - \dfrac{2 m_t(r)}{r} $. This is the extremization condition for the two-fluid system.

We define

Then, the condition becomes

From these two equations, we obtain a key constraint that reflects the effect of Lagrange multipliers in this two-fluid system:

which implies that

This reflects that the chemical potential between the two fluids satisfies a proportional relationship, which is somewhat different from the case of a single fluid [76].

We define a common function:

We differentiate the condition with respect to r:

where

We take the derivative in (32) and substitute the isentropic condition $ n_I' = \dfrac{n_I \varepsilon_I'}{p_I + \varepsilon_I} $:

We multiply by $ A^{3/2} $ and $ (p_I + \varepsilon_I)^2 $, and then rearrange to solve for $ p_I' $:

From the constraint definition $ \Gamma(r) = \lambda_I \dfrac{n_I}{p_I + \varepsilon_I} $, we obtain

Thus,

Substituting back into (35) yields

Now, we compute $ A' $: using $ m_t' = 4\pi r^2 \varepsilon_t $, we obtain

We substitute into the expression for $ p_I' $:

This is exactly the two-fluid equilibrium equation:

We have explicitly shown that the condition $ \delta H = 0 $ leads directly to the TOV equations. Conversely, since every step in the derivation above (substitution, differentiation, algebraic rearrangement) is reversible, if one starts from the TOV equations and chooses the constants $ \lambda_1,\lambda_2 $ such that $ \lambda_1/\lambda_2 = \mu_1(r)/\mu_2(r) $, one can integrate back to show that $ \delta H=0 $ (here $ \mu_I(r) \equiv (p_I+\varepsilon_I)/n_I $ is the local chemical potential of zero temperature fluid I). This completes the proof of necessity and sufficiency. □

The first theorem establishes a rigorous connection between the hydrostatic equilibrium and variational properties of the system, forming the foundation for Theorem III.2 and for simplifying the condition to the actual critical curve method. A similar proof of the variational principle for two-fluid systems was previously outlined in [77]. In the present work, we provide a self-contained and explicit derivation, introducing a unified auxiliary function $F(r) $ that consolidates the nonlocal integral contributions, and defining a common proportional function $\Gamma(r)=\lambda_I n_I/(p_I+\epsilon_I) $ to eliminate the Lagrange multipliers efficiently. This construction leads directly to the standard two-fluid TOV equations. Moreover, we establish the necessity and sufficiency of the relation between the variational extremum condition and the hydrostatic equilibrium, providing a complete and transparent foundation for the subsequent stability analysis.

Theorem III.1. A two-fluid isentropic mixed star configuration satisfies the coupled general-relativistic hydrostatic equilibrium equations (TOV equations) if and only if the total mass M takes its extremal value with respect to all time-independent variations that conserve the particle numbers of both fluids.

Proof. We provide the rigorous proof for Theorem III.1. The star configuration satisfies the equilibrium equations if and only if the quantity $ M $ defined by

is extremal under all variations of $ \varepsilon_I(r) $ that maintain the quantities

constant for $ I=1,2 $. This is equivalent to the existence of Lagrange multipliers $ \lambda_1 $ and $ \lambda_2 $ such that $ \delta H = \delta M - \lambda_1 \delta N_1 - \lambda_2 \delta N_2 $ is zero for all variations.

For given variations $ \delta \varepsilon_I(r) $, we obtain

We expand each term as follows (because the integrand vanishes above $ R_I+\delta R_I $, we write the upper limit of integration as positive infinity):

where $ m_t(r) = m_1(r) + m_2(r) $. Thus,

Using the isentropic conditions

and the mass variation

we perform substitution by steps. First, we substitute $ \delta n_I $:

Next, we substitute $ \delta m_J(r) $ into the final term:

By defining

and exchanging the order of integration in the double integral, the final term becomes

Combining all terms in $ \delta H $ yields

For $ \delta H $ to be zero for all variations $ \left[\left\{ \delta \varepsilon _I\left(r \right) \right\} \right] $, the integrand must vanish for each fluid component $ I $:

where $ A(r) = 1 - \dfrac{2 m_t(r)}{r} $. This is the extremization condition for the two-fluid system.

We define

Then, the condition becomes

From these two equations, we obtain a key constraint that reflects the effect of Lagrange multipliers in this two-fluid system:

which implies that

This reflects that the chemical potential between the two fluids satisfies a proportional relationship, which is somewhat different from the case of a single fluid [76].

We define a common function:

We differentiate the condition with respect to r:

where

We take the derivative in (32) and substitute the isentropic condition $ n_I' = \dfrac{n_I \varepsilon_I'}{p_I + \varepsilon_I} $:

We multiply by $ A^{3/2} $ and $ (p_I + \varepsilon_I)^2 $, and then rearrange to solve for $ p_I' $:

From the constraint definition $ \Gamma(r) = \lambda_I \dfrac{n_I}{p_I + \varepsilon_I} $, we obtain

Thus,

Substituting back into (35) yields

Now, we compute $ A' $: using $ m_t' = 4\pi r^2 \varepsilon_t $, we obtain

We substitute into the expression for $ p_I' $:

This is exactly the two-fluid equilibrium equation:

We have explicitly shown that the condition $ \delta H = 0 $ leads directly to the TOV equations. Conversely, since every step in the derivation above (substitution, differentiation, algebraic rearrangement) is reversible, if one starts from the TOV equations and chooses the constants $ \lambda_1,\lambda_2 $ such that $ \lambda_1/\lambda_2 = \mu_1(r)/\mu_2(r) $, one can integrate back to show that $ \delta H=0 $ (here $ \mu_I(r) \equiv (p_I+\varepsilon_I)/n_I $ is the local chemical potential of zero temperature fluid I). This completes the proof of necessity and sufficiency. □

This conjecture was initially proposed in the work of Henriques [8]. We demonstrate it by using conservation laws to extend a theorem that holds in a single-fluid spherically symmetric static compact star (see [76]) to the case of two fluids. A complete proof of the converse direction, together with a more comprehensive and mathematically rigorous analysis, are presented in a separate work [78].

Theorem III.2. For a two-fluid equilibrium configuration, the three gradients with respect to the central-pressure parameters, $ \nabla_{{\boldsymbol{\sigma}}} M $, $ \nabla_{{\boldsymbol{\sigma}}} N_{{\rm{NM}}} $, and $ \nabla_{{\boldsymbol{\sigma}}} N_{{\rm{DM}}} $, are mutually parallel in the central-pressure parameter space if there exists at least a zero mode with respect to radial perturbation.

Proof. If, for a certain static equilibrium configuration, $ n \in \mathbb{N} $ exists such that $ \omega_{n}^{2} = 0 $ under radial perturbation, the corresponding eigenmode is strictly time-independent:

This static solution satisfies the linearized field equations and corresponds to an infinitesimal variation between the equilibrium configuration $ {\boldsymbol{\sigma}} $ and a neighboring equilibrium configuration $ {\boldsymbol{\sigma}} + \delta{\boldsymbol{\sigma}} $, where $ \delta{\boldsymbol{\sigma}} = (\delta p_2^c, \delta p_1^c) $ is the direction in the central-pressure parameter space associated with the displacement field $ \xi_{In}(r) $.

The connection between the existence of static solutions and the critical curve relies on conserved quantities. Recalling the derivation of the two-fluid equilibrium and oscillation equations, the definitions of the physical quantities and their perturbations in terms of Lagrangian displacements originate from the Einstein field equations and covariant conservation of the particle number density for each fluid $ \nabla _{\mu}J_{I}^{\mu}=\nabla _{\mu}\left(n_Iu_{I}^{\mu} \right) =0 $, together with their linearizations. For a perturbation solution such as the zero-frequency mode, the linearized particle number conservation directly implies that the first-order variation of the total particle number vanishes:

We now explicitly verify this natural conclusion of particle number conservation:

● Symbol definition: Conserved charges and currents.

We define the conserved charges and currents as

where $\Sigma$ is a spacelike hypersurface of constant coordinate time $t$, $a_{\mu} = (-{\rm e}^{\phi},0,0,0)$ is the future-directed unit normal vector, $h$ is the determinant of the induced 3-metric, and $\sqrt{h}\,{\rm d}^{3}x = {\rm e}^{\nu} r^{2}\sin\theta\, {\rm d}r\,{\rm d}\theta\,{\rm d}\phi$.

● Explicit form of conserved charges:

We define $B = {\rm e}^{2\phi}, A = {\rm e}^{-2\nu}$. For $ i = N_I $, the conserved charge is explicitly expressed as

● Covariant conservation law:

The covariant conservation law is expressed as

● Proof of conservation:

Consider two adjacent spacelike hypersurfaces $\Sigma_1$ and $\Sigma_2$, which, together with a timelike boundary surface, form a closed hypersurface. As no current flows through the boundary, the flux of $ J_{N_I}^{\mu} $ through the boundary vanishes.

The change in the conserved charge from $\Sigma_1$ to $\Sigma_2$ is equal to the total flux of $ J_{N_I}^{\mu} $ through the closed hypersurface.

When applying the Stokes theorem, this flux equals the volume integral of the covariant divergence of $ J_{N_I}^{\mu} $ over the enclosed 4-volume:

From the covariant conservation law $ \nabla_{\mu} J_{N_I}^{\mu} = 0 $, we obtain

which implies that $ \Delta N_I = 0 $, thereby proving that the covariant conservation law leads to the global conservation of the total particle numbers $ N_1 $ and $ N_2 $.

Because the zero-frequency mode is a static solution of the radial perturbation equations and the equations incorporate the covariant conservation of the particle number density for each fluid, the linearized system directly forces the first-order variations of the particle numbers to vanish:

A similar conserved-current construction is not available for the mass variation: the Komar energy current built from a timelike Killing vector field does not yield the total gravitational mass M defined in our framework (the two differ by a metric-function factor), and its conservation relies on a global timelike Killing symmetry that is absent once radial perturbations are excited. Nevertheless, we can determine the first-order behaviour of M through the variational principle for the static zero-frequency mode. Theorem III.1 established that, for an equilibrium configuration, constants $ \lambda_1,\lambda_2 $ exist such that, for any static variation, $ \delta M - \lambda_1\,\delta N_1 - \lambda_2\,\delta N_2 = 0 $. For the zero-frequency static mode, the conservation $ \delta N_1= \delta N_2=0 $ reduces the variational relation to $ \delta M = 0 $ (to the first order).

Based on this, we can conclude that if an equilibrium configuration possesses a radial normal mode with a zero eigenvalue, this mode enters the tangent space of this configuration point in the equilibrium space, and the first order variations of the mass and particle number functionals vanish along this tangent vector, that is, $ \delta M[{\boldsymbol{\sigma}}; \delta{\boldsymbol{\sigma}}]=\delta N[{\boldsymbol{\sigma}}; \delta{\boldsymbol{\sigma}}]=0 $. Consequently, the directional derivatives of M and $ N_I $ along the same vector $ \delta{\boldsymbol{\sigma}} $ in parameter space are zero:

This is equivalent to the parallelism of the gradients of M and $ N_I $:

where the symbol $ \parallel $ denotes the parallel relation. □

Analysis and conclusion: As the Lagrange multipliers can be rescaled arbitrarily (e.g., by dividing $ \nabla M = \lambda_1 \nabla N_{NM} + \lambda_2 \nabla N_{DM} $ by $ \lambda_1 $ or $ \lambda_2 $), the statement of Theorem III.1 remains valid under any permutation of M, $ N_{NM} $, and $ N_{DM} $. This permutation symmetry implies that if any two of the three gradients ($ \nabla M $, $ \nabla N_{NM} $, and $ \nabla N_{DM} $) are parallel (indicating simultaneous extrema), the third must also be parallel to them. Hence, the critical curve can be determined by the parallelism of any two gradients.

Thus far, we have established that the stability boundary (the locus where the zero fundamental modes arise) constitutes a sufficient condition for the one-dimensional sequence derived via the critical curve method. Nevertheless, supplementary tests are necessary to verify whether the resulting curve indeed corresponds to the stability boundary. A straightforward approach involves conducting radial oscillation analyses on two sample configurations situated on opposite sides of the curve. Alternatively, by leveraging astrophysical insights (such as data from white dwarfs and NSs), we can often directly estimate the appropriate range of the central parameters, thereby enabling the direct determination of the stability boundary. Therefore, we can consider the static method as reliable and efficient for rigorous stability analysis.

This conjecture was initially proposed in the work of Henriques [8]. We demonstrate it by using conservation laws to extend a theorem that holds in a single-fluid spherically symmetric static compact star (see [76]) to the case of two fluids. A complete proof of the converse direction, together with a more comprehensive and mathematically rigorous analysis, are presented in a separate work [78].

Theorem III.2. For a two-fluid equilibrium configuration, the three gradients with respect to the central-pressure parameters, $ \nabla_{{\boldsymbol{\sigma}}} M $, $ \nabla_{{\boldsymbol{\sigma}}} N_{{\rm{NM}}} $, and $ \nabla_{{\boldsymbol{\sigma}}} N_{{\rm{DM}}} $, are mutually parallel in the central-pressure parameter space if there exists at least a zero mode with respect to radial perturbation.

Proof. If, for a certain static equilibrium configuration, $ n \in \mathbb{N} $ exists such that $ \omega_{n}^{2} = 0 $ under radial perturbation, the corresponding eigenmode is strictly time-independent:

This static solution satisfies the linearized field equations and corresponds to an infinitesimal variation between the equilibrium configuration $ {\boldsymbol{\sigma}} $ and a neighboring equilibrium configuration $ {\boldsymbol{\sigma}} + \delta{\boldsymbol{\sigma}} $, where $ \delta{\boldsymbol{\sigma}} = (\delta p_2^c, \delta p_1^c) $ is the direction in the central-pressure parameter space associated with the displacement field $ \xi_{In}(r) $.

The connection between the existence of static solutions and the critical curve relies on conserved quantities. Recalling the derivation of the two-fluid equilibrium and oscillation equations, the definitions of the physical quantities and their perturbations in terms of Lagrangian displacements originate from the Einstein field equations and covariant conservation of the particle number density for each fluid $ \nabla _{\mu}J_{I}^{\mu}=\nabla _{\mu}\left(n_Iu_{I}^{\mu} \right) =0 $, together with their linearizations. For a perturbation solution such as the zero-frequency mode, the linearized particle number conservation directly implies that the first-order variation of the total particle number vanishes:

We now explicitly verify this natural conclusion of particle number conservation:

● Symbol definition: Conserved charges and currents.

We define the conserved charges and currents as

where $\Sigma$ is a spacelike hypersurface of constant coordinate time $t$, $a_{\mu} = (-{\rm e}^{\phi},0,0,0)$ is the future-directed unit normal vector, $h$ is the determinant of the induced 3-metric, and $\sqrt{h}\,{\rm d}^{3}x = {\rm e}^{\nu} r^{2}\sin\theta\, {\rm d}r\,{\rm d}\theta\,{\rm d}\phi$.

● Explicit form of conserved charges:

We define $B = {\rm e}^{2\phi}, A = {\rm e}^{-2\nu}$. For $ i = N_I $, the conserved charge is explicitly expressed as

● Covariant conservation law:

The covariant conservation law is expressed as

● Proof of conservation:

Consider two adjacent spacelike hypersurfaces $\Sigma_1$ and $\Sigma_2$, which, together with a timelike boundary surface, form a closed hypersurface. As no current flows through the boundary, the flux of $ J_{N_I}^{\mu} $ through the boundary vanishes.

The change in the conserved charge from $\Sigma_1$ to $\Sigma_2$ is equal to the total flux of $ J_{N_I}^{\mu} $ through the closed hypersurface.

When applying the Stokes theorem, this flux equals the volume integral of the covariant divergence of $ J_{N_I}^{\mu} $ over the enclosed 4-volume:

From the covariant conservation law $ \nabla_{\mu} J_{N_I}^{\mu} = 0 $, we obtain

which implies that $ \Delta N_I = 0 $, thereby proving that the covariant conservation law leads to the global conservation of the total particle numbers $ N_1 $ and $ N_2 $.

Because the zero-frequency mode is a static solution of the radial perturbation equations and the equations incorporate the covariant conservation of the particle number density for each fluid, the linearized system directly forces the first-order variations of the particle numbers to vanish:

A similar conserved-current construction is not available for the mass variation: the Komar energy current built from a timelike Killing vector field does not yield the total gravitational mass M defined in our framework (the two differ by a metric-function factor), and its conservation relies on a global timelike Killing symmetry that is absent once radial perturbations are excited. Nevertheless, we can determine the first-order behaviour of M through the variational principle for the static zero-frequency mode. Theorem III.1 established that, for an equilibrium configuration, constants $ \lambda_1,\lambda_2 $ exist such that, for any static variation, $ \delta M - \lambda_1\,\delta N_1 - \lambda_2\,\delta N_2 = 0 $. For the zero-frequency static mode, the conservation $ \delta N_1= \delta N_2=0 $ reduces the variational relation to $ \delta M = 0 $ (to the first order).

Based on this, we can conclude that if an equilibrium configuration possesses a radial normal mode with a zero eigenvalue, this mode enters the tangent space of this configuration point in the equilibrium space, and the first order variations of the mass and particle number functionals vanish along this tangent vector, that is, $ \delta M[{\boldsymbol{\sigma}}; \delta{\boldsymbol{\sigma}}]=\delta N[{\boldsymbol{\sigma}}; \delta{\boldsymbol{\sigma}}]=0 $. Consequently, the directional derivatives of M and $ N_I $ along the same vector $ \delta{\boldsymbol{\sigma}} $ in parameter space are zero:

This is equivalent to the parallelism of the gradients of M and $ N_I $:

where the symbol $ \parallel $ denotes the parallel relation. □

Analysis and conclusion: As the Lagrange multipliers can be rescaled arbitrarily (e.g., by dividing $ \nabla M = \lambda_1 \nabla N_{NM} + \lambda_2 \nabla N_{DM} $ by $ \lambda_1 $ or $ \lambda_2 $), the statement of Theorem III.1 remains valid under any permutation of M, $ N_{NM} $, and $ N_{DM} $. This permutation symmetry implies that if any two of the three gradients ($ \nabla M $, $ \nabla N_{NM} $, and $ \nabla N_{DM} $) are parallel (indicating simultaneous extrema), the third must also be parallel to them. Hence, the critical curve can be determined by the parallelism of any two gradients.

Thus far, we have established that the stability boundary (the locus where the zero fundamental modes arise) constitutes a sufficient condition for the one-dimensional sequence derived via the critical curve method. Nevertheless, supplementary tests are necessary to verify whether the resulting curve indeed corresponds to the stability boundary. A straightforward approach involves conducting radial oscillation analyses on two sample configurations situated on opposite sides of the curve. Alternatively, by leveraging astrophysical insights (such as data from white dwarfs and NSs), we can often directly estimate the appropriate range of the central parameters, thereby enabling the direct determination of the stability boundary. Therefore, we can consider the static method as reliable and efficient for rigorous stability analysis.

This conjecture was initially proposed in the work of Henriques [8]. We demonstrate it by using conservation laws to extend a theorem that holds in a single-fluid spherically symmetric static compact star (see [76]) to the case of two fluids. A complete proof of the converse direction, together with a more comprehensive and mathematically rigorous analysis, are presented in a separate work [78].

Theorem III.2. For a two-fluid equilibrium configuration, the three gradients with respect to the central-pressure parameters, $ \nabla_{{\boldsymbol{\sigma}}} M $, $ \nabla_{{\boldsymbol{\sigma}}} N_{{\rm{NM}}} $, and $ \nabla_{{\boldsymbol{\sigma}}} N_{{\rm{DM}}} $, are mutually parallel in the central-pressure parameter space if there exists at least a zero mode with respect to radial perturbation.

Proof. If, for a certain static equilibrium configuration, $ n \in \mathbb{N} $ exists such that $ \omega_{n}^{2} = 0 $ under radial perturbation, the corresponding eigenmode is strictly time-independent:

This static solution satisfies the linearized field equations and corresponds to an infinitesimal variation between the equilibrium configuration $ {\boldsymbol{\sigma}} $ and a neighboring equilibrium configuration $ {\boldsymbol{\sigma}} + \delta{\boldsymbol{\sigma}} $, where $ \delta{\boldsymbol{\sigma}} = (\delta p_2^c, \delta p_1^c) $ is the direction in the central-pressure parameter space associated with the displacement field $ \xi_{In}(r) $.

The connection between the existence of static solutions and the critical curve relies on conserved quantities. Recalling the derivation of the two-fluid equilibrium and oscillation equations, the definitions of the physical quantities and their perturbations in terms of Lagrangian displacements originate from the Einstein field equations and covariant conservation of the particle number density for each fluid $ \nabla _{\mu}J_{I}^{\mu}=\nabla _{\mu}\left(n_Iu_{I}^{\mu} \right) =0 $, together with their linearizations. For a perturbation solution such as the zero-frequency mode, the linearized particle number conservation directly implies that the first-order variation of the total particle number vanishes:

We now explicitly verify this natural conclusion of particle number conservation:

● Symbol definition: Conserved charges and currents.

We define the conserved charges and currents as

where $\Sigma$ is a spacelike hypersurface of constant coordinate time $t$, $a_{\mu} = (-{\rm e}^{\phi},0,0,0)$ is the future-directed unit normal vector, $h$ is the determinant of the induced 3-metric, and $\sqrt{h}\,{\rm d}^{3}x = {\rm e}^{\nu} r^{2}\sin\theta\, {\rm d}r\,{\rm d}\theta\,{\rm d}\phi$.

● Explicit form of conserved charges:

We define $B = {\rm e}^{2\phi}, A = {\rm e}^{-2\nu}$. For $ i = N_I $, the conserved charge is explicitly expressed as

● Covariant conservation law:

The covariant conservation law is expressed as

● Proof of conservation:

Consider two adjacent spacelike hypersurfaces $\Sigma_1$ and $\Sigma_2$, which, together with a timelike boundary surface, form a closed hypersurface. As no current flows through the boundary, the flux of $ J_{N_I}^{\mu} $ through the boundary vanishes.

The change in the conserved charge from $\Sigma_1$ to $\Sigma_2$ is equal to the total flux of $ J_{N_I}^{\mu} $ through the closed hypersurface.

When applying the Stokes theorem, this flux equals the volume integral of the covariant divergence of $ J_{N_I}^{\mu} $ over the enclosed 4-volume:

From the covariant conservation law $ \nabla_{\mu} J_{N_I}^{\mu} = 0 $, we obtain

which implies that $ \Delta N_I = 0 $, thereby proving that the covariant conservation law leads to the global conservation of the total particle numbers $ N_1 $ and $ N_2 $.

Because the zero-frequency mode is a static solution of the radial perturbation equations and the equations incorporate the covariant conservation of the particle number density for each fluid, the linearized system directly forces the first-order variations of the particle numbers to vanish:

A similar conserved-current construction is not available for the mass variation: the Komar energy current built from a timelike Killing vector field does not yield the total gravitational mass M defined in our framework (the two differ by a metric-function factor), and its conservation relies on a global timelike Killing symmetry that is absent once radial perturbations are excited. Nevertheless, we can determine the first-order behaviour of M through the variational principle for the static zero-frequency mode. Theorem III.1 established that, for an equilibrium configuration, constants $ \lambda_1,\lambda_2 $ exist such that, for any static variation, $ \delta M - \lambda_1\,\delta N_1 - \lambda_2\,\delta N_2 = 0 $. For the zero-frequency static mode, the conservation $ \delta N_1= \delta N_2=0 $ reduces the variational relation to $ \delta M = 0 $ (to the first order).

Based on this, we can conclude that if an equilibrium configuration possesses a radial normal mode with a zero eigenvalue, this mode enters the tangent space of this configuration point in the equilibrium space, and the first order variations of the mass and particle number functionals vanish along this tangent vector, that is, $ \delta M[{\boldsymbol{\sigma}}; \delta{\boldsymbol{\sigma}}]=\delta N[{\boldsymbol{\sigma}}; \delta{\boldsymbol{\sigma}}]=0 $. Consequently, the directional derivatives of M and $ N_I $ along the same vector $ \delta{\boldsymbol{\sigma}} $ in parameter space are zero:

This is equivalent to the parallelism of the gradients of M and $ N_I $:

where the symbol $ \parallel $ denotes the parallel relation. □

Analysis and conclusion: As the Lagrange multipliers can be rescaled arbitrarily (e.g., by dividing $ \nabla M = \lambda_1 \nabla N_{NM} + \lambda_2 \nabla N_{DM} $ by $ \lambda_1 $ or $ \lambda_2 $), the statement of Theorem III.1 remains valid under any permutation of M, $ N_{NM} $, and $ N_{DM} $. This permutation symmetry implies that if any two of the three gradients ($ \nabla M $, $ \nabla N_{NM} $, and $ \nabla N_{DM} $) are parallel (indicating simultaneous extrema), the third must also be parallel to them. Hence, the critical curve can be determined by the parallelism of any two gradients.

Thus far, we have established that the stability boundary (the locus where the zero fundamental modes arise) constitutes a sufficient condition for the one-dimensional sequence derived via the critical curve method. Nevertheless, supplementary tests are necessary to verify whether the resulting curve indeed corresponds to the stability boundary. A straightforward approach involves conducting radial oscillation analyses on two sample configurations situated on opposite sides of the curve. Alternatively, by leveraging astrophysical insights (such as data from white dwarfs and NSs), we can often directly estimate the appropriate range of the central parameters, thereby enabling the direct determination of the stability boundary. Therefore, we can consider the static method as reliable and efficient for rigorous stability analysis.

The stability of an equilibrium configuration, which is an extremum of M, is determined by the second-order variation $ \delta^2 M $. The connection to dynamical stability is provided by the following thermodynamic relation [77, 79]:

where $ {\boldsymbol{\xi}} $ is the perturbation mode. The stability conditions are

● $ \delta^{2} M \gt 0 \Leftrightarrow \omega_0^2 \gt 0 $: Stable configuration (energy minimum).

● $ \delta^{2} M = 0 \Leftrightarrow \omega_0^2 = 0 $: Stability boundary (marginal stability).

● $ \delta^{2} M \lt 0 \Leftrightarrow \omega_0^2 \lt 0 $: Unstable configuration (saddle point).

This framework provides the physical intuition: a stable star resides at the bottom of an energy well. The stability boundary is reached when the well becomes flat in a certain direction, allowing for a transition to another state without an energy cost, which is precisely the correspondence of a zero-frequency fundamental mode.

The stability of an equilibrium configuration, which is an extremum of M, is determined by the second-order variation $ \delta^2 M $. The connection to dynamical stability is provided by the following thermodynamic relation [77, 79]:

where $ {\boldsymbol{\xi}} $ is the perturbation mode. The stability conditions are

● $ \delta^{2} M \gt 0 \Leftrightarrow \omega_0^2 \gt 0 $: Stable configuration (energy minimum).

● $ \delta^{2} M = 0 \Leftrightarrow \omega_0^2 = 0 $: Stability boundary (marginal stability).

● $ \delta^{2} M \lt 0 \Leftrightarrow \omega_0^2 \lt 0 $: Unstable configuration (saddle point).

This framework provides the physical intuition: a stable star resides at the bottom of an energy well. The stability boundary is reached when the well becomes flat in a certain direction, allowing for a transition to another state without an energy cost, which is precisely the correspondence of a zero-frequency fundamental mode.

The stability of an equilibrium configuration, which is an extremum of M, is determined by the second-order variation $ \delta^2 M $. The connection to dynamical stability is provided by the following thermodynamic relation [77, 79]:

where $ {\boldsymbol{\xi}} $ is the perturbation mode. The stability conditions are

● $ \delta^{2} M \gt 0 \Leftrightarrow \omega_0^2 \gt 0 $: Stable configuration (energy minimum).

● $ \delta^{2} M = 0 \Leftrightarrow \omega_0^2 = 0 $: Stability boundary (marginal stability).

● $ \delta^{2} M \lt 0 \Leftrightarrow \omega_0^2 \lt 0 $: Unstable configuration (saddle point).

This framework provides the physical intuition: a stable star resides at the bottom of an energy well. The stability boundary is reached when the well becomes flat in a certain direction, allowing for a transition to another state without an energy cost, which is precisely the correspondence of a zero-frequency fundamental mode.

This section presents the numerical results of our investigation into the stability properties and macroscopic characteristics of mixed compact stars. We start with a detailed stability analysis comparing single- and multi-fluid models, and then examine the resulting M-R relationships and their implications for astrophysical observations.

This section presents the numerical results of our investigation into the stability properties and macroscopic characteristics of mixed compact stars. We start with a detailed stability analysis comparing single- and multi-fluid models, and then examine the resulting M-R relationships and their implications for astrophysical observations.

This section presents the numerical results of our investigation into the stability properties and macroscopic characteristics of mixed compact stars. We start with a detailed stability analysis comparing single- and multi-fluid models, and then examine the resulting M-R relationships and their implications for astrophysical observations.

As shown in Fig. 3, we provide the eigenvalue-central pressure relations for single-fluid models with three different EoSs: SLy4, Bosonic DM (with $ B_4=0.1 $), and Holographic NM (with $ \ell^{-7}=10300 $). The stability distribution and M-R relations satisfy the BTM criterion. These numerical results serve as comparisons for the situation of mixed components.

Figure 3.  (color online) Single-fluid stability analysis. The left panel (Fig. 3(a)) shows the fundamental eigenvalue curves as functions of the central configuration for different EoSs: Bosonic DM ($ B_4 = 0.1 $), SLy4, and Holographic NM ($ \ell ^{-7}=10300 $). The right panel (Fig. 3(b)) displays the corresponding M-R relation curves.

As shown in Fig. 3, we provide the eigenvalue-central pressure relations for single-fluid models with three different EoSs: SLy4, Bosonic DM (with $ B_4=0.1 $), and Holographic NM (with $ \ell^{-7}=10300 $). The stability distribution and M-R relations satisfy the BTM criterion. These numerical results serve as comparisons for the situation of mixed components.

Figure 3.  (color online) Single-fluid stability analysis. The left panel (Fig. 3(a)) shows the fundamental eigenvalue curves as functions of the central configuration for different EoSs: Bosonic DM ($ B_4 = 0.1 $), SLy4, and Holographic NM ($ \ell ^{-7}=10300 $). The right panel (Fig. 3(b)) displays the corresponding M-R relation curves.

As shown in Fig. 3, we provide the eigenvalue-central pressure relations for single-fluid models with three different EoSs: SLy4, Bosonic DM (with $ B_4=0.1 $), and Holographic NM (with $ \ell^{-7}=10300 $). The stability distribution and M-R relations satisfy the BTM criterion. These numerical results serve as comparisons for the situation of mixed components.

Figure 3.  (color online) Single-fluid stability analysis. The left panel (Fig. 3(a)) shows the fundamental eigenvalue curves as functions of the central configuration for different EoSs: Bosonic DM ($ B_4 = 0.1 $), SLy4, and Holographic NM ($ \ell ^{-7}=10300 $). The right panel (Fig. 3(b)) displays the corresponding M-R relation curves.

The numerical results in Fig. 4(a) provide precise squared frequencies of the fundamental radial oscillation using Bosonic DM ($ B_4 = 0.1 $) and Holographic NM ($ \ell ^{-7}=10300 $). Thus, in addition to identifying stable regions, the results reveal the trend of fundamental frequency variation, showing non-monotonic behavior.

Figure 4.  (color online) Radial oscillation analysis results for the mixed-star model with EoS combinations of Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4 = 0.1 $). The left panel (Fig. 4(a)) shows the continuous distribution of the eigenvalues. The right panel (Fig. 4(b)) displays the distribution of stable and unstable configuration points in the central-pressure parameter space, where the black curve represents configurations in which the fundamental radial-mode frequency vanishes and its interior indicates the stable region.

To identify the stability boundary in the static approach numerically, we need to locate the points where the relevant gradients become parallel. Throughout this subsection, the gradients are defined in the two-dimensional central-pressure parameter space in Eq. (12), rather than with respect to the radial coordinate inside the star. Thus, the gravitational mass M and particle numbers $ N_I $ are treated as scalar functions on this parameter space. Previous work has usually relied on connecting the tangent points of several contour lines of mass and particle number [53, 80]. We adopt a more efficient numerical method by computing the normalized signed cross product between gradients. Specifically, for two scalar quantities $ X(\sigma) $ and $ Y(\sigma) $, selected from M, $ N_1 $, and $ N_2 $, we define

where $ \hat{u}_{\nabla_\sigma X}=\nabla_\sigma X/|\nabla_\sigma X| $ is the normalized gradient direction of X in the central-pressure parameter space. The magnitude $ |{\cal{C}}_{X,Y}| $ ranges from 0 to 1. Points with $ |{\cal{C}}_{X,Y}|=0 $ correspond to parallel gradients and form the critical curve, whereas points with $ |{\cal{C}}_{X,Y}|=1 $ correspond to orthogonal gradients. Figure 5 shows the resulting critical curves for different EoS combinations using the method introduced in Subsection A. In these plots, the black curves mark the zero contours $ |{\cal{C}}_{X,Y}|=0 $, namely the locations of the critical curves, whereas the white curves mark the contours $ |{\cal{C}}_{X,Y}|=1 $.

Figure 5.  (color online) Critical curves for different EoS combinations. The black curves indicate multiple stability boundaries. The units of both the horizontal and vertical axes are astrophysical units ($ p_\odot $), representing the magnitudes of the central pressure for the two components. The black and white heatmap in the background displays the absolute value of the normalized cross product of the gradient directions at each point of the total mass field and the particle number field of the second component; the darker the color, the more parallel the gradient directions of the two fields at that point tend to be. $ N_1\equiv N_{DM} $ represents the DM particle number.

Figure 5(a) shows Holographic NM with $ \ell ^{-7}=10300 $ and Bosonic DM with $ B_4 = 0.1 $. Figure 5(b) depics SLy4 and Bosonic DM with $ B _4=0.03 $. Figure 5(c) shows Holographic NM with $ \ell ^{-7}=10300 $ and Fermionic DM with $ m_f=0.5 $. Figure 5(d) presents SLy4 and Fermionic DM with $ m_f=0.5 $. Figures 4(b) and 5(a) clearly confirm that, for these EoS models, the critical curve obtained from the static parallel-gradient condition coincides with the boundary where the fundamental radial mode frequency vanishes. This numerically demonstrates that the zero mode identified by the critical curve method is indeed the fundamental mode, and thus, the curve constitutes the true stability boundary for these cases.

Figure 6 further illustrates the stability analysis using the signed cross product $ {\cal{C}}_{M,N_1} $. By replacing the absolute value used in Fig. 5(a) with the signed value of the cross product for the same EoS and parameters, we obtain the heatmap shown in Fig. 6(a). The green points where $ {\cal{C}}_{M,N_1}=0 $ mark the static critical curve. In Fig. 6(b), this static boundary is directly compared with the dynamical boundary obtained from the vanishing of the fundamental radial-mode frequency. The agreement between the two curves provides a numerical verification of the equivalence between the static criterion and dynamical stability condition for this EoS combination.

Figure 6.  (color online) Signed-cross-product critical curve and its comparison with the dynamical stability boundary. EoS combination: Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4 = 0.1 $).

Figure 7 extends the analysis of Fig. 5(a) to a larger parameter space. In this enlarged domain, the static gradient-alignment diagnostic reveals additional nested curve-like structures beside the innermost critical curve. These outer curves are interpreted as possible higher-order critical curves, as suggested by the macroscopic gradient-alignment structure of the equilibrium-configuration surface.

Figure 7.  (color online) Static gradient-alignment structures in an enlarged central-pressure parameter space. The units of both the horizontal and vertical axes are astrophysical pressure units ($ p_\odot $). This figure extends Fig. 5(a) using the same EoS combination. The dark regions indicate small values of $ |{\cal{C}}_{M,N_{{\rm{DM}}}}| $, and the black curves denote the zero contours $ |{\cal{C}}_{M,N_{{\rm{DM}}}}|=0 $. The innermost black curve corresponds to the stability boundary associated with the fundamental radial mode. The outer nested black curves are interpreted as possible higher-order critical curves, as suggested by the static gradient-alignment structure.

Figure 8(a) shows the M-R curves for mixed stars with a fixed Bosonic DM central pressure $ p_1^c $ and varying Holographic NM ($ \ell ^{-7}=10300 $) central pressure $ p_2^c $. These relations, obtained through single-parameter variations in mixed star models, no longer satisfy the BTM criterion. Furthermore, the deviation from the BTM criterion can be quantitatively analyzed, as demonstrated in Fig. 8(b). This panel plots the relative separation between the central pressure of the last stable point, $ p_{\text{last-stable}} $, and that of the maximum mass configuration, $ p_{\text{max-mass}} $, as a function of the DM properties, showing the trend of the systematic separation as it varies with the central pressure and properties of DM.

Figure 8.  (color online) Deviation from the BTM criterion in mixed stars.

Crucially, the magnitude of this deviation is highly sensitive to the EoS properties of the admixed DM. As shown, different DM models (e.g., Bosonic vs. Fermionic DM) and variations in their microphysical parameters (such as the interaction strength $ B_4 $ or fermion mass $ m_f $) lead to distinctly different deviation curves. This suggests that the degree to which the BTM criterion is violated is not universal, but rather, an indirect probe of the DM EoS and its gravitational coupling with NM. Deviations of the neutral stability point from the standard turning-point criterion under various additional physical effects have also been discussed in earlier literature [80−82].

The deviation of the stability boundary in the multi-fluid configuration space has direct implications for gravitational-wave observables. In particular, the stability boundary directly determines the maximum achievable compactness $ M/R $ and characteristic density scale $ M/R^{3} $ of stable configurations. As the fundamental f-mode frequency scales as $ f_f \sim \sqrt{M/R^{3}} $ (see, e.g., Refs. [83−85]), precise measurements of f-mode frequencies from gravitational-wave observations provide a direct investigation of this characteristic scale. Consequently, the deviation of the mixed-star stability boundary from the single-fluid maximum-mass point translates into a distinct upper-limit shift in the predicted $ M/R^{3} $ scale, which can be tested against future gravitational-wave data. That is, the configuration that maximizes the f-mode frequency is no longer necessarily aligned with the maximum-mass point of single-fluid sequences, but instead, follows from the full multi-dimensional stability condition in the central-pressure parameter space of the mixed system.

Therefore, in principle, a precise future measurement of the mass, radius, and stability threshold of a mixed star could provide a tool for constraining the nature of DM within the star.

The numerical results in Fig. 4(a) provide precise squared frequencies of the fundamental radial oscillation using Bosonic DM ($ B_4 = 0.1 $) and Holographic NM ($ \ell ^{-7}=10300 $). Thus, in addition to identifying stable regions, the results reveal the trend of fundamental frequency variation, showing non-monotonic behavior.

Figure 4.  (color online) Radial oscillation analysis results for the mixed-star model with EoS combinations of Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4 = 0.1 $). The left panel (Fig. 4(a)) shows the continuous distribution of the eigenvalues. The right panel (Fig. 4(b)) displays the distribution of stable and unstable configuration points in the central-pressure parameter space, where the black curve represents configurations in which the fundamental radial-mode frequency vanishes and its interior indicates the stable region.

To identify the stability boundary in the static approach numerically, we need to locate the points where the relevant gradients become parallel. Throughout this subsection, the gradients are defined in the two-dimensional central-pressure parameter space in Eq. (12), rather than with respect to the radial coordinate inside the star. Thus, the gravitational mass M and particle numbers $ N_I $ are treated as scalar functions on this parameter space. Previous work has usually relied on connecting the tangent points of several contour lines of mass and particle number [53, 80]. We adopt a more efficient numerical method by computing the normalized signed cross product between gradients. Specifically, for two scalar quantities $ X(\sigma) $ and $ Y(\sigma) $, selected from M, $ N_1 $, and $ N_2 $, we define

where $ \hat{u}_{\nabla_\sigma X}=\nabla_\sigma X/|\nabla_\sigma X| $ is the normalized gradient direction of X in the central-pressure parameter space. The magnitude $ |{\cal{C}}_{X,Y}| $ ranges from 0 to 1. Points with $ |{\cal{C}}_{X,Y}|=0 $ correspond to parallel gradients and form the critical curve, whereas points with $ |{\cal{C}}_{X,Y}|=1 $ correspond to orthogonal gradients. Figure 5 shows the resulting critical curves for different EoS combinations using the method introduced in Subsection A. In these plots, the black curves mark the zero contours $ |{\cal{C}}_{X,Y}|=0 $, namely the locations of the critical curves, whereas the white curves mark the contours $ |{\cal{C}}_{X,Y}|=1 $.

Figure 5.  (color online) Critical curves for different EoS combinations. The black curves indicate multiple stability boundaries. The units of both the horizontal and vertical axes are astrophysical units ($ p_\odot $), representing the magnitudes of the central pressure for the two components. The black and white heatmap in the background displays the absolute value of the normalized cross product of the gradient directions at each point of the total mass field and the particle number field of the second component; the darker the color, the more parallel the gradient directions of the two fields at that point tend to be. $ N_1\equiv N_{DM} $ represents the DM particle number.

Figure 5(a) shows Holographic NM with $ \ell ^{-7}=10300 $ and Bosonic DM with $ B_4 = 0.1 $. Figure 5(b) depics SLy4 and Bosonic DM with $ B _4=0.03 $. Figure 5(c) shows Holographic NM with $ \ell ^{-7}=10300 $ and Fermionic DM with $ m_f=0.5 $. Figure 5(d) presents SLy4 and Fermionic DM with $ m_f=0.5 $. Figures 4(b) and 5(a) clearly confirm that, for these EoS models, the critical curve obtained from the static parallel-gradient condition coincides with the boundary where the fundamental radial mode frequency vanishes. This numerically demonstrates that the zero mode identified by the critical curve method is indeed the fundamental mode, and thus, the curve constitutes the true stability boundary for these cases.

Figure 6 further illustrates the stability analysis using the signed cross product $ {\cal{C}}_{M,N_1} $. By replacing the absolute value used in Fig. 5(a) with the signed value of the cross product for the same EoS and parameters, we obtain the heatmap shown in Fig. 6(a). The green points where $ {\cal{C}}_{M,N_1}=0 $ mark the static critical curve. In Fig. 6(b), this static boundary is directly compared with the dynamical boundary obtained from the vanishing of the fundamental radial-mode frequency. The agreement between the two curves provides a numerical verification of the equivalence between the static criterion and dynamical stability condition for this EoS combination.

Figure 6.  (color online) Signed-cross-product critical curve and its comparison with the dynamical stability boundary. EoS combination: Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4 = 0.1 $).

Figure 7 extends the analysis of Fig. 5(a) to a larger parameter space. In this enlarged domain, the static gradient-alignment diagnostic reveals additional nested curve-like structures beside the innermost critical curve. These outer curves are interpreted as possible higher-order critical curves, as suggested by the macroscopic gradient-alignment structure of the equilibrium-configuration surface.

Figure 7.  (color online) Static gradient-alignment structures in an enlarged central-pressure parameter space. The units of both the horizontal and vertical axes are astrophysical pressure units ($ p_\odot $). This figure extends Fig. 5(a) using the same EoS combination. The dark regions indicate small values of $ |{\cal{C}}_{M,N_{{\rm{DM}}}}| $, and the black curves denote the zero contours $ |{\cal{C}}_{M,N_{{\rm{DM}}}}|=0 $. The innermost black curve corresponds to the stability boundary associated with the fundamental radial mode. The outer nested black curves are interpreted as possible higher-order critical curves, as suggested by the static gradient-alignment structure.

Figure 8(a) shows the M-R curves for mixed stars with a fixed Bosonic DM central pressure $ p_1^c $ and varying Holographic NM ($ \ell ^{-7}=10300 $) central pressure $ p_2^c $. These relations, obtained through single-parameter variations in mixed star models, no longer satisfy the BTM criterion. Furthermore, the deviation from the BTM criterion can be quantitatively analyzed, as demonstrated in Fig. 8(b). This panel plots the relative separation between the central pressure of the last stable point, $ p_{\text{last-stable}} $, and that of the maximum mass configuration, $ p_{\text{max-mass}} $, as a function of the DM properties, showing the trend of the systematic separation as it varies with the central pressure and properties of DM.

Figure 8.  (color online) Deviation from the BTM criterion in mixed stars.

Crucially, the magnitude of this deviation is highly sensitive to the EoS properties of the admixed DM. As shown, different DM models (e.g., Bosonic vs. Fermionic DM) and variations in their microphysical parameters (such as the interaction strength $ B_4 $ or fermion mass $ m_f $) lead to distinctly different deviation curves. This suggests that the degree to which the BTM criterion is violated is not universal, but rather, an indirect probe of the DM EoS and its gravitational coupling with NM. Deviations of the neutral stability point from the standard turning-point criterion under various additional physical effects have also been discussed in earlier literature [80−82].

The deviation of the stability boundary in the multi-fluid configuration space has direct implications for gravitational-wave observables. In particular, the stability boundary directly determines the maximum achievable compactness $ M/R $ and characteristic density scale $ M/R^{3} $ of stable configurations. As the fundamental f-mode frequency scales as $ f_f \sim \sqrt{M/R^{3}} $ (see, e.g., Refs. [83−85]), precise measurements of f-mode frequencies from gravitational-wave observations provide a direct investigation of this characteristic scale. Consequently, the deviation of the mixed-star stability boundary from the single-fluid maximum-mass point translates into a distinct upper-limit shift in the predicted $ M/R^{3} $ scale, which can be tested against future gravitational-wave data. That is, the configuration that maximizes the f-mode frequency is no longer necessarily aligned with the maximum-mass point of single-fluid sequences, but instead, follows from the full multi-dimensional stability condition in the central-pressure parameter space of the mixed system.

Therefore, in principle, a precise future measurement of the mass, radius, and stability threshold of a mixed star could provide a tool for constraining the nature of DM within the star.

The numerical results in Fig. 4(a) provide precise squared frequencies of the fundamental radial oscillation using Bosonic DM ($ B_4 = 0.1 $) and Holographic NM ($ \ell ^{-7}=10300 $). Thus, in addition to identifying stable regions, the results reveal the trend of fundamental frequency variation, showing non-monotonic behavior.

Figure 4.  (color online) Radial oscillation analysis results for the mixed-star model with EoS combinations of Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4 = 0.1 $). The left panel (Fig. 4(a)) shows the continuous distribution of the eigenvalues. The right panel (Fig. 4(b)) displays the distribution of stable and unstable configuration points in the central-pressure parameter space, where the black curve represents configurations in which the fundamental radial-mode frequency vanishes and its interior indicates the stable region.

To identify the stability boundary in the static approach numerically, we need to locate the points where the relevant gradients become parallel. Throughout this subsection, the gradients are defined in the two-dimensional central-pressure parameter space in Eq. (12), rather than with respect to the radial coordinate inside the star. Thus, the gravitational mass M and particle numbers $ N_I $ are treated as scalar functions on this parameter space. Previous work has usually relied on connecting the tangent points of several contour lines of mass and particle number [53, 80]. We adopt a more efficient numerical method by computing the normalized signed cross product between gradients. Specifically, for two scalar quantities $ X(\sigma) $ and $ Y(\sigma) $, selected from M, $ N_1 $, and $ N_2 $, we define

where $ \hat{u}_{\nabla_\sigma X}=\nabla_\sigma X/|\nabla_\sigma X| $ is the normalized gradient direction of X in the central-pressure parameter space. The magnitude $ |{\cal{C}}_{X,Y}| $ ranges from 0 to 1. Points with $ |{\cal{C}}_{X,Y}|=0 $ correspond to parallel gradients and form the critical curve, whereas points with $ |{\cal{C}}_{X,Y}|=1 $ correspond to orthogonal gradients. Figure 5 shows the resulting critical curves for different EoS combinations using the method introduced in Subsection A. In these plots, the black curves mark the zero contours $ |{\cal{C}}_{X,Y}|=0 $, namely the locations of the critical curves, whereas the white curves mark the contours $ |{\cal{C}}_{X,Y}|=1 $.

Figure 5.  (color online) Critical curves for different EoS combinations. The black curves indicate multiple stability boundaries. The units of both the horizontal and vertical axes are astrophysical units ($ p_\odot $), representing the magnitudes of the central pressure for the two components. The black and white heatmap in the background displays the absolute value of the normalized cross product of the gradient directions at each point of the total mass field and the particle number field of the second component; the darker the color, the more parallel the gradient directions of the two fields at that point tend to be. $ N_1\equiv N_{DM} $ represents the DM particle number.

Figure 5(a) shows Holographic NM with $ \ell ^{-7}=10300 $ and Bosonic DM with $ B_4 = 0.1 $. Figure 5(b) depics SLy4 and Bosonic DM with $ B _4=0.03 $. Figure 5(c) shows Holographic NM with $ \ell ^{-7}=10300 $ and Fermionic DM with $ m_f=0.5 $. Figure 5(d) presents SLy4 and Fermionic DM with $ m_f=0.5 $. Figures 4(b) and 5(a) clearly confirm that, for these EoS models, the critical curve obtained from the static parallel-gradient condition coincides with the boundary where the fundamental radial mode frequency vanishes. This numerically demonstrates that the zero mode identified by the critical curve method is indeed the fundamental mode, and thus, the curve constitutes the true stability boundary for these cases.

Figure 6 further illustrates the stability analysis using the signed cross product $ {\cal{C}}_{M,N_1} $. By replacing the absolute value used in Fig. 5(a) with the signed value of the cross product for the same EoS and parameters, we obtain the heatmap shown in Fig. 6(a). The green points where $ {\cal{C}}_{M,N_1}=0 $ mark the static critical curve. In Fig. 6(b), this static boundary is directly compared with the dynamical boundary obtained from the vanishing of the fundamental radial-mode frequency. The agreement between the two curves provides a numerical verification of the equivalence between the static criterion and dynamical stability condition for this EoS combination.

Figure 6.  (color online) Signed-cross-product critical curve and its comparison with the dynamical stability boundary. EoS combination: Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4 = 0.1 $).

Figure 7 extends the analysis of Fig. 5(a) to a larger parameter space. In this enlarged domain, the static gradient-alignment diagnostic reveals additional nested curve-like structures beside the innermost critical curve. These outer curves are interpreted as possible higher-order critical curves, as suggested by the macroscopic gradient-alignment structure of the equilibrium-configuration surface.

Figure 7.  (color online) Static gradient-alignment structures in an enlarged central-pressure parameter space. The units of both the horizontal and vertical axes are astrophysical pressure units ($ p_\odot $). This figure extends Fig. 5(a) using the same EoS combination. The dark regions indicate small values of $ |{\cal{C}}_{M,N_{{\rm{DM}}}}| $, and the black curves denote the zero contours $ |{\cal{C}}_{M,N_{{\rm{DM}}}}|=0 $. The innermost black curve corresponds to the stability boundary associated with the fundamental radial mode. The outer nested black curves are interpreted as possible higher-order critical curves, as suggested by the static gradient-alignment structure.

Figure 8(a) shows the M-R curves for mixed stars with a fixed Bosonic DM central pressure $ p_1^c $ and varying Holographic NM ($ \ell ^{-7}=10300 $) central pressure $ p_2^c $. These relations, obtained through single-parameter variations in mixed star models, no longer satisfy the BTM criterion. Furthermore, the deviation from the BTM criterion can be quantitatively analyzed, as demonstrated in Fig. 8(b). This panel plots the relative separation between the central pressure of the last stable point, $ p_{\text{last-stable}} $, and that of the maximum mass configuration, $ p_{\text{max-mass}} $, as a function of the DM properties, showing the trend of the systematic separation as it varies with the central pressure and properties of DM.

Figure 8.  (color online) Deviation from the BTM criterion in mixed stars.

Crucially, the magnitude of this deviation is highly sensitive to the EoS properties of the admixed DM. As shown, different DM models (e.g., Bosonic vs. Fermionic DM) and variations in their microphysical parameters (such as the interaction strength $ B_4 $ or fermion mass $ m_f $) lead to distinctly different deviation curves. This suggests that the degree to which the BTM criterion is violated is not universal, but rather, an indirect probe of the DM EoS and its gravitational coupling with NM. Deviations of the neutral stability point from the standard turning-point criterion under various additional physical effects have also been discussed in earlier literature [80−82].

The deviation of the stability boundary in the multi-fluid configuration space has direct implications for gravitational-wave observables. In particular, the stability boundary directly determines the maximum achievable compactness $ M/R $ and characteristic density scale $ M/R^{3} $ of stable configurations. As the fundamental f-mode frequency scales as $ f_f \sim \sqrt{M/R^{3}} $ (see, e.g., Refs. [83−85]), precise measurements of f-mode frequencies from gravitational-wave observations provide a direct investigation of this characteristic scale. Consequently, the deviation of the mixed-star stability boundary from the single-fluid maximum-mass point translates into a distinct upper-limit shift in the predicted $ M/R^{3} $ scale, which can be tested against future gravitational-wave data. That is, the configuration that maximizes the f-mode frequency is no longer necessarily aligned with the maximum-mass point of single-fluid sequences, but instead, follows from the full multi-dimensional stability condition in the central-pressure parameter space of the mixed system.

Therefore, in principle, a precise future measurement of the mass, radius, and stability threshold of a mixed star could provide a tool for constraining the nature of DM within the star.

To verify that the critical curve obtained from the static gradient-alignment algorithm is independent of numerical grid resolution, we recompute the curve using four different grid sizes in the central-pressure parameter space: $ 40\times40 $, $ 60\times60 $, $ 80\times80 $, and $ 100\times100 $. The EoS combination and parameters are the same as those in Fig. 5(a), that is, Holographic NM ($ \ell^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $). As shown in Fig. 9, the resulting critical curves are indistinguishable, confirming that the algorithm is resolution-free.

Figure 9.  (color online) Critical curves computed with four different grid resolutions ($ 40\times40 $, $ 60\times60 $, $ 80\times80 $, $ 100\times100 $) for the same EoS combination and parameters as those in Fig. 5(a) (Holographic NM $ \ell^{-7}=10300 $, Bosonic DM $ B_4=0.1 $). The curves coincide, demonstrating that the static method is resolution-free.

To verify that the critical curve obtained from the static gradient-alignment algorithm is independent of numerical grid resolution, we recompute the curve using four different grid sizes in the central-pressure parameter space: $ 40\times40 $, $ 60\times60 $, $ 80\times80 $, and $ 100\times100 $. The EoS combination and parameters are the same as those in Fig. 5(a), that is, Holographic NM ($ \ell^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $). As shown in Fig. 9, the resulting critical curves are indistinguishable, confirming that the algorithm is resolution-free.

Figure 9.  (color online) Critical curves computed with four different grid resolutions ($ 40\times40 $, $ 60\times60 $, $ 80\times80 $, $ 100\times100 $) for the same EoS combination and parameters as those in Fig. 5(a) (Holographic NM $ \ell^{-7}=10300 $, Bosonic DM $ B_4=0.1 $). The curves coincide, demonstrating that the static method is resolution-free.

To verify that the critical curve obtained from the static gradient-alignment algorithm is independent of numerical grid resolution, we recompute the curve using four different grid sizes in the central-pressure parameter space: $ 40\times40 $, $ 60\times60 $, $ 80\times80 $, and $ 100\times100 $. The EoS combination and parameters are the same as those in Fig. 5(a), that is, Holographic NM ($ \ell^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $). As shown in Fig. 9, the resulting critical curves are indistinguishable, confirming that the algorithm is resolution-free.

Figure 9.  (color online) Critical curves computed with four different grid resolutions ($ 40\times40 $, $ 60\times60 $, $ 80\times80 $, $ 100\times100 $) for the same EoS combination and parameters as those in Fig. 5(a) (Holographic NM $ \ell^{-7}=10300 $, Bosonic DM $ B_4=0.1 $). The curves coincide, demonstrating that the static method is resolution-free.

The stable region of equilibrium configurations maps to a continuous topological correspondence region in the M-$ R_t $-$ p_{I}^{c} $ space: a two-dimensional surface in three-dimensional space, as shown in Fig. 10. The EoS combination in Fig. 10 is Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $), with vertical axis $ p_{1}^{c} $ in Fig. 10(a) and $ p_{2}^{c} $ in Fig. 10(b). $ R_t $ is the total radius: $ R_t= \max(R_N,R_D) $. Here, $ R_N $ and $ R_D $ represent the radii of NM and DM, respectively. The red and blue dots indicate the configurations where $ R_N>R_D $ and $ R_N<R_D $, respectively. Any slice along the $ p_{I}^{c} $ axis corresponds to an M-R curve governed by a single parameter variation.

Figure 10.  (color online) The two-dimensional stable region from the central pressure parameter space is mapped onto a three-dimensional space, with the axes representing the total mass (in $ M_\odot $ units), total radius $R_t = \max\,(R_N, R_D)$ (in $\rm km$ units), and one of the central pressures (e.g., $ p_{1}^{c} $ and $ p_{2}^{c} $ in $ p_\odot $ units). This reveals the degeneracies of the two-dimensional M-R projection. The red and blue dots distinguish configurations in which the NM radius is larger than the DM radius ($ R_N \gt R_D $, blue) from those in which it is smaller ($ R_N \lt R_D $, red). (a) and (b) show the two-dimensional surfaces of the DM and NM pressures for the EoS combination in Fig. 5(a), namely Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $).

Figure 11 shows the projections of stable configurations onto both the mass-visible radius plane M-$ R_N $ and mass-total radius plane M-$ R_t $ with the EoS combination of Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $). Points within the stability boundary form a "stable group" in the projection plane, serving as a comparison tool for observational data. Overlapping regions of two colors imply the possible existence of "twin configurations" with identical masses and radii but distinct halo compositions. The dynamical stability of such twin configurations under non-linear perturbations has recently been investigated in Ref. [87].

Figure 11.  (color online) In both panels, $ R_N $ and $ R_D $ are the radii (in $ km $ units) where the NM and DM pressures vanish, respectively. The colored regions represent dynamically stable configurations, whereas the gray shaded areas are unstable. The color coding within the stable group distinguishes dark core configurations ($ R_N \gt R_D $, red) from dark halo configurations ($ R_N \lt R_D $, blue). The overlapping area between the red and blue regions, highlighted in green, indicates the parameter space in which twin stars [86] can exist: configurations with identical masses and radii but different internal structures (one being a dark core and the other a dark halo). The dashed line represents the boundary of the stable region (the critical curve).

This degeneracy between twin configurations can be resolved by combining gravitational-wave and electromagnetic observations. Gravitational-wave measurements constrain the stellar radius through tidal deformability [88], whereas electromagnetic observations such as pulse-profile modeling provide independent constraints on the surface emission geometry [89]. A systematic discrepancy between the inferred radii, particularly $ R_{{\rm{GW}}} \gt R_{{\rm{EM}}} $, may indicate the presence of an extended low-density component or a multi-fluid internal structure, offering a potential observational signature of twin stars [90].

The left boundary line of the "stable group" in Fig. 11(b) coincides with the curve where the radius size relationship between the two fluid components reverses, which means that, for mixed star models with the same mass, the stable configuration that minimizes the total star radius is always the one in which the two fluid components have equal radii.

The stable region of equilibrium configurations maps to a continuous topological correspondence region in the M-$ R_t $-$ p_{I}^{c} $ space: a two-dimensional surface in three-dimensional space, as shown in Fig. 10. The EoS combination in Fig. 10 is Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $), with vertical axis $ p_{1}^{c} $ in Fig. 10(a) and $ p_{2}^{c} $ in Fig. 10(b). $ R_t $ is the total radius: $ R_t= \max(R_N,R_D) $. Here, $ R_N $ and $ R_D $ represent the radii of NM and DM, respectively. The red and blue dots indicate the configurations where $ R_N>R_D $ and $ R_N<R_D $, respectively. Any slice along the $ p_{I}^{c} $ axis corresponds to an M-R curve governed by a single parameter variation.

Figure 10.  (color online) The two-dimensional stable region from the central pressure parameter space is mapped onto a three-dimensional space, with the axes representing the total mass (in $ M_\odot $ units), total radius $R_t = \max\,(R_N, R_D)$ (in $\rm km$ units), and one of the central pressures (e.g., $ p_{1}^{c} $ and $ p_{2}^{c} $ in $ p_\odot $ units). This reveals the degeneracies of the two-dimensional M-R projection. The red and blue dots distinguish configurations in which the NM radius is larger than the DM radius ($ R_N \gt R_D $, blue) from those in which it is smaller ($ R_N \lt R_D $, red). (a) and (b) show the two-dimensional surfaces of the DM and NM pressures for the EoS combination in Fig. 5(a), namely Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $).

Figure 11 shows the projections of stable configurations onto both the mass-visible radius plane M-$ R_N $ and mass-total radius plane M-$ R_t $ with the EoS combination of Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $). Points within the stability boundary form a "stable group" in the projection plane, serving as a comparison tool for observational data. Overlapping regions of two colors imply the possible existence of "twin configurations" with identical masses and radii but distinct halo compositions. The dynamical stability of such twin configurations under non-linear perturbations has recently been investigated in Ref. [87].

Figure 11.  (color online) In both panels, $ R_N $ and $ R_D $ are the radii (in $ km $ units) where the NM and DM pressures vanish, respectively. The colored regions represent dynamically stable configurations, whereas the gray shaded areas are unstable. The color coding within the stable group distinguishes dark core configurations ($ R_N \gt R_D $, red) from dark halo configurations ($ R_N \lt R_D $, blue). The overlapping area between the red and blue regions, highlighted in green, indicates the parameter space in which twin stars [86] can exist: configurations with identical masses and radii but different internal structures (one being a dark core and the other a dark halo). The dashed line represents the boundary of the stable region (the critical curve).

This degeneracy between twin configurations can be resolved by combining gravitational-wave and electromagnetic observations. Gravitational-wave measurements constrain the stellar radius through tidal deformability [88], whereas electromagnetic observations such as pulse-profile modeling provide independent constraints on the surface emission geometry [89]. A systematic discrepancy between the inferred radii, particularly $ R_{{\rm{GW}}} \gt R_{{\rm{EM}}} $, may indicate the presence of an extended low-density component or a multi-fluid internal structure, offering a potential observational signature of twin stars [90].

The left boundary line of the "stable group" in Fig. 11(b) coincides with the curve where the radius size relationship between the two fluid components reverses, which means that, for mixed star models with the same mass, the stable configuration that minimizes the total star radius is always the one in which the two fluid components have equal radii.

The stable region of equilibrium configurations maps to a continuous topological correspondence region in the M-$ R_t $-$ p_{I}^{c} $ space: a two-dimensional surface in three-dimensional space, as shown in Fig. 10. The EoS combination in Fig. 10 is Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $), with vertical axis $ p_{1}^{c} $ in Fig. 10(a) and $ p_{2}^{c} $ in Fig. 10(b). $ R_t $ is the total radius: $ R_t= \max(R_N,R_D) $. Here, $ R_N $ and $ R_D $ represent the radii of NM and DM, respectively. The red and blue dots indicate the configurations where $ R_N>R_D $ and $ R_N<R_D $, respectively. Any slice along the $ p_{I}^{c} $ axis corresponds to an M-R curve governed by a single parameter variation.

Figure 10.  (color online) The two-dimensional stable region from the central pressure parameter space is mapped onto a three-dimensional space, with the axes representing the total mass (in $ M_\odot $ units), total radius $R_t = \max\,(R_N, R_D)$ (in $\rm km$ units), and one of the central pressures (e.g., $ p_{1}^{c} $ and $ p_{2}^{c} $ in $ p_\odot $ units). This reveals the degeneracies of the two-dimensional M-R projection. The red and blue dots distinguish configurations in which the NM radius is larger than the DM radius ($ R_N \gt R_D $, blue) from those in which it is smaller ($ R_N \lt R_D $, red). (a) and (b) show the two-dimensional surfaces of the DM and NM pressures for the EoS combination in Fig. 5(a), namely Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $).

Figure 11 shows the projections of stable configurations onto both the mass-visible radius plane M-$ R_N $ and mass-total radius plane M-$ R_t $ with the EoS combination of Holographic NM ($ \ell ^{-7}=10300 $) and Bosonic DM ($ B_4=0.1 $). Points within the stability boundary form a "stable group" in the projection plane, serving as a comparison tool for observational data. Overlapping regions of two colors imply the possible existence of "twin configurations" with identical masses and radii but distinct halo compositions. The dynamical stability of such twin configurations under non-linear perturbations has recently been investigated in Ref. [87].

Figure 11.  (color online) In both panels, $ R_N $ and $ R_D $ are the radii (in $ km $ units) where the NM and DM pressures vanish, respectively. The colored regions represent dynamically stable configurations, whereas the gray shaded areas are unstable. The color coding within the stable group distinguishes dark core configurations ($ R_N \gt R_D $, red) from dark halo configurations ($ R_N \lt R_D $, blue). The overlapping area between the red and blue regions, highlighted in green, indicates the parameter space in which twin stars [86] can exist: configurations with identical masses and radii but different internal structures (one being a dark core and the other a dark halo). The dashed line represents the boundary of the stable region (the critical curve).

This degeneracy between twin configurations can be resolved by combining gravitational-wave and electromagnetic observations. Gravitational-wave measurements constrain the stellar radius through tidal deformability [88], whereas electromagnetic observations such as pulse-profile modeling provide independent constraints on the surface emission geometry [89]. A systematic discrepancy between the inferred radii, particularly $ R_{{\rm{GW}}} \gt R_{{\rm{EM}}} $, may indicate the presence of an extended low-density component or a multi-fluid internal structure, offering a potential observational signature of twin stars [90].

The left boundary line of the "stable group" in Fig. 11(b) coincides with the curve where the radius size relationship between the two fluid components reverses, which means that, for mixed star models with the same mass, the stable configuration that minimizes the total star radius is always the one in which the two fluid components have equal radii.

We employ the strongly self-repulsive Bosonic DM EoS introduced in Sec. II.C to discuss the observational implications of mixed stars. As shown by Colpi, Shapiro, and Wasserman [74], a complex scalar field with a quartic self-interaction potential $ V(\phi) = \dfrac{1}{4}\lambda |\phi|^4 $ can support compact equilibrium configurations with a maximum mass scaling $ M_{\max} \simeq 0.10\ \text{GeV}^2\, M_\odot\, \lambda^{1/2} / m^2 $. For $ \lambda \gg 1 $ (strong self-repulsion), this maximum mass becomes comparable to the Chandrasekhar mass of fermionic stars, making it possible to reach the extreme mass ranges discussed below. Therefore, strong self-repulsion (encoded in a small value of B in Eq. (10)) is essential for the mixed-star model to access the large masses shown in Fig. 12. Based on this justification, we now examine the predicted M-R relations.

Figure 12.  (color online) Prediction prospects for admixed star models displayed through the M-R relationship, with the EoS combination of Bosonic DM (10) and SLy4.

Collecting the latest M-R credible window from observational data of compact stars on the M-R stable groups can impose constraints on the DM EoS. Our mixed star model demonstrates the ability to cover a wide range of the M-R area by adjusting the EoS parameters. For instance, Fig. 11(a) covers the mass and radius values of the PSR J0437-4715 data from the latest NICER observations ($ M = 1.418 \pm 0.037 M_\odot $, $ R = 11.36^{+0.95}_{-0.63} $ km) [68].

Figure 12 and Table 1 show the predicted range of the M-R relations and compactness analysis of mixed stars under different EoS parameter selections and DM mixing ratios ($ ratio\equiv p^c_{\rm DM}/p^c_{\rm NM} $, the central pressure ratio of DM and NM), indicating that by appropriately adjusting the hardness of the EoS and ratios of the admixed DM, the mixed star models have the potential to become predictive for NSs and certain gap events, containing ranges of both 3−5 $ M_{\odot} $ and 50−100 $ M_{\odot} $ [32, 91]. Although such high-mass objects have conventionally been interpreted as black holes formed via hierarchical mergers or primordial mechanisms [92, 93], our model suggests an alternative: they could be DM-dominated exotic compact stars [94, 95]. The key question is whether these objects might exhibit non-black-hole signatures in gravitational-wave observations. Future third-generation detectors (e.g., Einstein Telescope and Cosmic Explorer) might be able to test these predictions by measuring the tidal Love numbers and quasi-normal mode frequencies of candidate events in the upper mass gap [96, 97]. Figure 12(b) shows the situations of massive stars with the observable radius of conventional compact stars and accompanied by a large dark halo. Additionally, the similar values of compactness $ M/R $ for the same EoS combinations but with different parameters and mixing ratios exhibit a slight violation of the scaling symmetry for a certain type of single EoS [32, 98].

We employ the strongly self-repulsive Bosonic DM EoS introduced in Sec. II.C to discuss the observational implications of mixed stars. As shown by Colpi, Shapiro, and Wasserman [74], a complex scalar field with a quartic self-interaction potential $ V(\phi) = \dfrac{1}{4}\lambda |\phi|^4 $ can support compact equilibrium configurations with a maximum mass scaling $ M_{\max} \simeq 0.10\ \text{GeV}^2\, M_\odot\, \lambda^{1/2} / m^2 $. For $ \lambda \gg 1 $ (strong self-repulsion), this maximum mass becomes comparable to the Chandrasekhar mass of fermionic stars, making it possible to reach the extreme mass ranges discussed below. Therefore, strong self-repulsion (encoded in a small value of B in Eq. (10)) is essential for the mixed-star model to access the large masses shown in Fig. 12. Based on this justification, we now examine the predicted M-R relations.

Figure 12.  (color online) Prediction prospects for admixed star models displayed through the M-R relationship, with the EoS combination of Bosonic DM (10) and SLy4.

Collecting the latest M-R credible window from observational data of compact stars on the M-R stable groups can impose constraints on the DM EoS. Our mixed star model demonstrates the ability to cover a wide range of the M-R area by adjusting the EoS parameters. For instance, Fig. 11(a) covers the mass and radius values of the PSR J0437-4715 data from the latest NICER observations ($ M = 1.418 \pm 0.037 M_\odot $, $ R = 11.36^{+0.95}_{-0.63} $ km) [68].

Figure 12 and Table 1 show the predicted range of the M-R relations and compactness analysis of mixed stars under different EoS parameter selections and DM mixing ratios ($ ratio\equiv p^c_{\rm DM}/p^c_{\rm NM} $, the central pressure ratio of DM and NM), indicating that by appropriately adjusting the hardness of the EoS and ratios of the admixed DM, the mixed star models have the potential to become predictive for NSs and certain gap events, containing ranges of both 3−5 $ M_{\odot} $ and 50−100 $ M_{\odot} $ [32, 91]. Although such high-mass objects have conventionally been interpreted as black holes formed via hierarchical mergers or primordial mechanisms [92, 93], our model suggests an alternative: they could be DM-dominated exotic compact stars [94, 95]. The key question is whether these objects might exhibit non-black-hole signatures in gravitational-wave observations. Future third-generation detectors (e.g., Einstein Telescope and Cosmic Explorer) might be able to test these predictions by measuring the tidal Love numbers and quasi-normal mode frequencies of candidate events in the upper mass gap [96, 97]. Figure 12(b) shows the situations of massive stars with the observable radius of conventional compact stars and accompanied by a large dark halo. Additionally, the similar values of compactness $ M/R $ for the same EoS combinations but with different parameters and mixing ratios exhibit a slight violation of the scaling symmetry for a certain type of single EoS [32, 98].

We employ the strongly self-repulsive Bosonic DM EoS introduced in Sec. II.C to discuss the observational implications of mixed stars. As shown by Colpi, Shapiro, and Wasserman [74], a complex scalar field with a quartic self-interaction potential $ V(\phi) = \dfrac{1}{4}\lambda |\phi|^4 $ can support compact equilibrium configurations with a maximum mass scaling $ M_{\max} \simeq 0.10\ \text{GeV}^2\, M_\odot\, \lambda^{1/2} / m^2 $. For $ \lambda \gg 1 $ (strong self-repulsion), this maximum mass becomes comparable to the Chandrasekhar mass of fermionic stars, making it possible to reach the extreme mass ranges discussed below. Therefore, strong self-repulsion (encoded in a small value of B in Eq. (10)) is essential for the mixed-star model to access the large masses shown in Fig. 12. Based on this justification, we now examine the predicted M-R relations.

Figure 12.  (color online) Prediction prospects for admixed star models displayed through the M-R relationship, with the EoS combination of Bosonic DM (10) and SLy4.

Collecting the latest M-R credible window from observational data of compact stars on the M-R stable groups can impose constraints on the DM EoS. Our mixed star model demonstrates the ability to cover a wide range of the M-R area by adjusting the EoS parameters. For instance, Fig. 11(a) covers the mass and radius values of the PSR J0437-4715 data from the latest NICER observations ($ M = 1.418 \pm 0.037 M_\odot $, $ R = 11.36^{+0.95}_{-0.63} $ km) [68].

Figure 12 and Table 1 show the predicted range of the M-R relations and compactness analysis of mixed stars under different EoS parameter selections and DM mixing ratios ($ ratio\equiv p^c_{\rm DM}/p^c_{\rm NM} $, the central pressure ratio of DM and NM), indicating that by appropriately adjusting the hardness of the EoS and ratios of the admixed DM, the mixed star models have the potential to become predictive for NSs and certain gap events, containing ranges of both 3−5 $ M_{\odot} $ and 50−100 $ M_{\odot} $ [32, 91]. Although such high-mass objects have conventionally been interpreted as black holes formed via hierarchical mergers or primordial mechanisms [92, 93], our model suggests an alternative: they could be DM-dominated exotic compact stars [94, 95]. The key question is whether these objects might exhibit non-black-hole signatures in gravitational-wave observations. Future third-generation detectors (e.g., Einstein Telescope and Cosmic Explorer) might be able to test these predictions by measuring the tidal Love numbers and quasi-normal mode frequencies of candidate events in the upper mass gap [96, 97]. Figure 12(b) shows the situations of massive stars with the observable radius of conventional compact stars and accompanied by a large dark halo. Additionally, the similar values of compactness $ M/R $ for the same EoS combinations but with different parameters and mixing ratios exhibit a slight violation of the scaling symmetry for a certain type of single EoS [32, 98].

In this study, we have comprehensively investigated the stability properties and macroscopic characteristics of DM-NM mixed stars within the framework of general relativity.

We have established a rigorous equivalence between the critical curve method and the emergence of a zero-frequency radial oscillation mode. This foundational link provides a formal justification for using the static critical curve method as a powerful proxy for dynamical stability analysis.

Our numerical calculations, employing four typical EoS combinations—SLy4, Holographic NM, Fermionic DM, and Bosonic DM—determined the stability regions of mixed stars and validated the equivalence between the two methods numerically while revealing complex stability distributions that differ fundamentally from those of single-fluid and hybrid stars. Unlike pure NSs that obey the classical BTM criterion, stable configurations of mixed stars form continuous two-dimensional surfaces in the three-dimensional parameter space spanned by (one of) the central pressure, mass, and radius. The projection of these stable configurations onto the M-R plane reveals a "stable group" that delineates natural regions in which mixed stars should be observationally found. These groups exhibit several noteworthy features, including the existence of twin configurations with identical M-R measurements but different internal compositions, and the characteristic property that the most compact stable configurations are those in which the two fluid components have equal radii. Moreover, by adjusting the EoS hardness and DM mixing ratio, our mixed star model can accommodate regular NS events and account for certain mass-gap events. It also has the potential to predict massive configurations featuring a compact star radius core accompanied by a large DM halo.

In the future, our theoretical framework can be extended to incorporate more sophisticated DM models, such as self-interacting fermionic DM and axion-like particles. The methodology also provides a foundation for investigating rotating mixed stars, magnetic field effects, and dynamical formation scenarios. As observational capabilities continue to advance, the stable groups and characteristic M-R relations identified through our analysis may serve as references for identifying potential mixed star candidates and constraining the properties of DM through astrophysical observations.

In this study, we have comprehensively investigated the stability properties and macroscopic characteristics of DM-NM mixed stars within the framework of general relativity.

We have established a rigorous equivalence between the critical curve method and the emergence of a zero-frequency radial oscillation mode. This foundational link provides a formal justification for using the static critical curve method as a powerful proxy for dynamical stability analysis.

Our numerical calculations, employing four typical EoS combinations—SLy4, Holographic NM, Fermionic DM, and Bosonic DM—determined the stability regions of mixed stars and validated the equivalence between the two methods numerically while revealing complex stability distributions that differ fundamentally from those of single-fluid and hybrid stars. Unlike pure NSs that obey the classical BTM criterion, stable configurations of mixed stars form continuous two-dimensional surfaces in the three-dimensional parameter space spanned by (one of) the central pressure, mass, and radius. The projection of these stable configurations onto the M-R plane reveals a "stable group" that delineates natural regions in which mixed stars should be observationally found. These groups exhibit several noteworthy features, including the existence of twin configurations with identical M-R measurements but different internal compositions, and the characteristic property that the most compact stable configurations are those in which the two fluid components have equal radii. Moreover, by adjusting the EoS hardness and DM mixing ratio, our mixed star model can accommodate regular NS events and account for certain mass-gap events. It also has the potential to predict massive configurations featuring a compact star radius core accompanied by a large DM halo.

In the future, our theoretical framework can be extended to incorporate more sophisticated DM models, such as self-interacting fermionic DM and axion-like particles. The methodology also provides a foundation for investigating rotating mixed stars, magnetic field effects, and dynamical formation scenarios. As observational capabilities continue to advance, the stable groups and characteristic M-R relations identified through our analysis may serve as references for identifying potential mixed star candidates and constraining the properties of DM through astrophysical observations.

In this study, we have comprehensively investigated the stability properties and macroscopic characteristics of DM-NM mixed stars within the framework of general relativity.

We have established a rigorous equivalence between the critical curve method and the emergence of a zero-frequency radial oscillation mode. This foundational link provides a formal justification for using the static critical curve method as a powerful proxy for dynamical stability analysis.

Our numerical calculations, employing four typical EoS combinations—SLy4, Holographic NM, Fermionic DM, and Bosonic DM—determined the stability regions of mixed stars and validated the equivalence between the two methods numerically while revealing complex stability distributions that differ fundamentally from those of single-fluid and hybrid stars. Unlike pure NSs that obey the classical BTM criterion, stable configurations of mixed stars form continuous two-dimensional surfaces in the three-dimensional parameter space spanned by (one of) the central pressure, mass, and radius. The projection of these stable configurations onto the M-R plane reveals a "stable group" that delineates natural regions in which mixed stars should be observationally found. These groups exhibit several noteworthy features, including the existence of twin configurations with identical M-R measurements but different internal compositions, and the characteristic property that the most compact stable configurations are those in which the two fluid components have equal radii. Moreover, by adjusting the EoS hardness and DM mixing ratio, our mixed star model can accommodate regular NS events and account for certain mass-gap events. It also has the potential to predict massive configurations featuring a compact star radius core accompanied by a large DM halo.

In the future, our theoretical framework can be extended to incorporate more sophisticated DM models, such as self-interacting fermionic DM and axion-like particles. The methodology also provides a foundation for investigating rotating mixed stars, magnetic field effects, and dynamical formation scenarios. As observational capabilities continue to advance, the stable groups and characteristic M-R relations identified through our analysis may serve as references for identifying potential mixed star candidates and constraining the properties of DM through astrophysical observations.

The static stability criteria determine the stability of a star by examining the properties of the entire equilibrium sequence, rather than analyzing perturbations of a single configuration. This approach is computationally efficient as it reuses the data for calculating the equilibrium models.

The static stability criteria determine the stability of a star by examining the properties of the entire equilibrium sequence, rather than analyzing perturbations of a single configuration. This approach is computationally efficient as it reuses the data for calculating the equilibrium models.

The static stability criteria determine the stability of a star by examining the properties of the entire equilibrium sequence, rather than analyzing perturbations of a single configuration. This approach is computationally efficient as it reuses the data for calculating the equilibrium models.

The BTM criterion provides a method for assessing the radial stability of general-relativistic star models. The key principle of this criterion (we apply the generalized version in [51]) is that when passing through each extremum of the M-R curve, the direction can either increase or decrease the core pressure:

● Stability changes occur exclusively at the extrema of the M-R curve.

● A counterclockwise bend at an extremum corresponds to gaining one unstable mode.

● A clockwise bend at an extremum corresponds to losing one unstable mode.

● The sequence is typically parameterized by the central pressure.

The original BTM criterion was formulated specifically for sequences of zero-temperature stars, representing a "family" of different stars composed of the same cold matter. This framework is directly applicable to the cold NS models considered in this study. By contrast, for sequences of hot isentropic stars with a fixed total baryon number, the stability change at an extremum follows the opposite rule [3]: a counterclockwise bend leads to the loss of an unstable mode, whereas a clockwise bend leads to the gain of one.

Owing to its intuitive nature, the BTM criterion has been widely used for decades as a primary tool for stability assessment in relativistic astrophysics. The equivalence between the oscillation method and BTM criterion has been proven in [99, 100].

Notably, the BTM criterion is specifically applicable to single-fluid models and becomes invalid for multi-fluid models containing interface or mixed phases. We provide a more detailed discussion of this fact and its manifestation in subsequent sections.

The BTM criterion provides a method for assessing the radial stability of general-relativistic star models. The key principle of this criterion (we apply the generalized version in [51]) is that when passing through each extremum of the M-R curve, the direction can either increase or decrease the core pressure:

● Stability changes occur exclusively at the extrema of the M-R curve.

● A counterclockwise bend at an extremum corresponds to gaining one unstable mode.

● A clockwise bend at an extremum corresponds to losing one unstable mode.

● The sequence is typically parameterized by the central pressure.

The original BTM criterion was formulated specifically for sequences of zero-temperature stars, representing a "family" of different stars composed of the same cold matter. This framework is directly applicable to the cold NS models considered in this study. By contrast, for sequences of hot isentropic stars with a fixed total baryon number, the stability change at an extremum follows the opposite rule [3]: a counterclockwise bend leads to the loss of an unstable mode, whereas a clockwise bend leads to the gain of one.

Owing to its intuitive nature, the BTM criterion has been widely used for decades as a primary tool for stability assessment in relativistic astrophysics. The equivalence between the oscillation method and BTM criterion has been proven in [99, 100].

Notably, the BTM criterion is specifically applicable to single-fluid models and becomes invalid for multi-fluid models containing interface or mixed phases. We provide a more detailed discussion of this fact and its manifestation in subsequent sections.

The BTM criterion provides a method for assessing the radial stability of general-relativistic star models. The key principle of this criterion (we apply the generalized version in [51]) is that when passing through each extremum of the M-R curve, the direction can either increase or decrease the core pressure:

● Stability changes occur exclusively at the extrema of the M-R curve.

● A counterclockwise bend at an extremum corresponds to gaining one unstable mode.

● A clockwise bend at an extremum corresponds to losing one unstable mode.

● The sequence is typically parameterized by the central pressure.

The original BTM criterion was formulated specifically for sequences of zero-temperature stars, representing a "family" of different stars composed of the same cold matter. This framework is directly applicable to the cold NS models considered in this study. By contrast, for sequences of hot isentropic stars with a fixed total baryon number, the stability change at an extremum follows the opposite rule [3]: a counterclockwise bend leads to the loss of an unstable mode, whereas a clockwise bend leads to the gain of one.

Owing to its intuitive nature, the BTM criterion has been widely used for decades as a primary tool for stability assessment in relativistic astrophysics. The equivalence between the oscillation method and BTM criterion has been proven in [99, 100].

Notably, the BTM criterion is specifically applicable to single-fluid models and becomes invalid for multi-fluid models containing interface or mixed phases. We provide a more detailed discussion of this fact and its manifestation in subsequent sections.

The stability analysis of multi-fluid compact stars can be significantly simplified through the application of the critical curve method. This approach, which was originally developed in [8] for determining the stability of boson-fermion stars, is in fact applicable to any spherically symmetric ideal multi-fluid star model in general relativity, including the DM-NM mixed stars considered in this study. The method provides an efficient criterion for determining star stability by relying solely on equilibrium solutions of the multi-fluid TOV equations, thereby circumventing the computationally intensive task of solving the coupled pulsation equations to identify the radial normal modes.

The fundamental premise of this method lies in characterizing a stability boundary in the parameter space of the central pressures for the constituent fluids. For a two-fluid system, the critical curve demarcating stable and unstable configurations is formed by a group of points characterizing the central conditions. At each point $ {\boldsymbol{\sigma}}=\left(p_{2}^{c}, p_{1}^{c} \right) $, there exists a direction $ \delta{\boldsymbol{\sigma}} $ that fits the condition [8]:

where $M(\infty)$ denotes the total gravitational mass of the star (equivalent to the mass within radius $R$, hereafter denoted simply as $M$), and $ N_{F} $ and $ N_{F} $ are the total particle numbers for fermions and bosons, respectively. For the specific case of DM-NM mixed stars considered in our study, this criterion takes the following form:

where $N_{\rm NM}$ and $N_{\rm DM}$ represent the particle numbers of DM and NM, respectively.

The stability analysis of multi-fluid compact stars can be significantly simplified through the application of the critical curve method. This approach, which was originally developed in [8] for determining the stability of boson-fermion stars, is in fact applicable to any spherically symmetric ideal multi-fluid star model in general relativity, including the DM-NM mixed stars considered in this study. The method provides an efficient criterion for determining star stability by relying solely on equilibrium solutions of the multi-fluid TOV equations, thereby circumventing the computationally intensive task of solving the coupled pulsation equations to identify the radial normal modes.

The fundamental premise of this method lies in characterizing a stability boundary in the parameter space of the central pressures for the constituent fluids. For a two-fluid system, the critical curve demarcating stable and unstable configurations is formed by a group of points characterizing the central conditions. At each point $ {\boldsymbol{\sigma}}=\left(p_{2}^{c}, p_{1}^{c} \right) $, there exists a direction $ \delta{\boldsymbol{\sigma}} $ that fits the condition [8]:

where $M(\infty)$ denotes the total gravitational mass of the star (equivalent to the mass within radius $R$, hereafter denoted simply as $M$), and $ N_{F} $ and $ N_{F} $ are the total particle numbers for fermions and bosons, respectively. For the specific case of DM-NM mixed stars considered in our study, this criterion takes the following form:

where $N_{\rm NM}$ and $N_{\rm DM}$ represent the particle numbers of DM and NM, respectively.

The stability analysis of multi-fluid compact stars can be significantly simplified through the application of the critical curve method. This approach, which was originally developed in [8] for determining the stability of boson-fermion stars, is in fact applicable to any spherically symmetric ideal multi-fluid star model in general relativity, including the DM-NM mixed stars considered in this study. The method provides an efficient criterion for determining star stability by relying solely on equilibrium solutions of the multi-fluid TOV equations, thereby circumventing the computationally intensive task of solving the coupled pulsation equations to identify the radial normal modes.

The fundamental premise of this method lies in characterizing a stability boundary in the parameter space of the central pressures for the constituent fluids. For a two-fluid system, the critical curve demarcating stable and unstable configurations is formed by a group of points characterizing the central conditions. At each point $ {\boldsymbol{\sigma}}=\left(p_{2}^{c}, p_{1}^{c} \right) $, there exists a direction $ \delta{\boldsymbol{\sigma}} $ that fits the condition [8]:

where $M(\infty)$ denotes the total gravitational mass of the star (equivalent to the mass within radius $R$, hereafter denoted simply as $M$), and $ N_{F} $ and $ N_{F} $ are the total particle numbers for fermions and bosons, respectively. For the specific case of DM-NM mixed stars considered in our study, this criterion takes the following form:

where $N_{\rm NM}$ and $N_{\rm DM}$ represent the particle numbers of DM and NM, respectively.

In his 1964 paper, Chandrasekhar rigorously established the pulsation equation for radial oscillations in single-fluid stars based on radial perturbations and variational principles [2], and it has been widely applied in many works, such as [101]. For the necessity of unifying notation, we briefly review the main steps and key equations in the derivation of this equation:

Linearization of Field Equations

An unperturbed, static, spherically symmetric equilibrium configuration satisfying the TOV equations is assumed. Its behavior after being subjected to infinitesimal radial perturbations can be used to determine its stability. We define ${\rm e}^{2\nu \left(r \right)}\equiv \left(1-\dfrac{2Gm\left(r \right)}{rc^2} \right) ^{-1}$. Then, the metric and fluid variables are expressed as

The Einstein field equations and Einstein-Euler equations are perturbed separately. Neglecting higher-order terms, we obtain the corresponding system of linearized radial perturbation equations:

Lagrangian Displacement and Reformulating Perturbations

The Lagrangian displacement ξ is defined such that the radial velocity satisfies $ v=\dfrac{{\rm d}r}{{\rm d}t}=\dfrac{\partial \xi}{\partial t} $, that is, $u^r={\rm e}^{-\phi_0}\dfrac{{\rm d}r}{{\rm d}t}= {\rm e}^{-\phi_0}\dfrac{\partial \xi}{\partial t}$. Then, using the first three equations of the linearized perturbation equations (B1), (B2), and (B3), all perturbations of the metric and fluid variables except for $ \delta p $ are expressed in terms of ξ and equilibrium quantities:

Notably, the perturbation quantities and Lagrangian displacement are all functions of $ r $ and $ t $ (e.g., $ \xi(r, t) $, $ \delta p(r, t) $). We now assume that all perturbation quantities $ F $ possess a time-harmonic dependence of the form $ F(r, t) = F(r) e^{i\omega t} $. Based on this separation of variables and by combining (B4) with (B5), (B6), and (B7), we obtain

where the perturbation equations are thus reduced from complex partial differential equations to a relatively simpler eigenvalue problem governed by an ordinary differential equation, with the perturbation terms now consisting of only the pressure variation.

Baryon Number Conservation

An additional relation is required to connect the pressure perturbation $ \delta p $ to the other variables. This is provided by imposing adiabaticity and baryon number conservation, which requires $ \nabla_\mu (n u^\mu) = 0 $. Defining the perturbation of $ n $ as $ n = n_0(r) + \delta n(r, t) $ and combining it with ξ defined previously leads to the form of $ \delta n $ constituted by the Lagrangian displacement and unperturbed terms:

Moreover, the first law of thermodynamics yields ${\rm d}\varepsilon = \dfrac{\varepsilon + p}{n} {\rm d}n$; thus, the baryon number density is now related to $ p $ or ε. This finally allows the pressure perturbation to be written solely in terms of the Lagrangian displacement:

Here, γ is defined as $\gamma =-\dfrac{n_0}{p_0}\dfrac{{\rm d}p_0}{{\rm d}n_0}$, which is the adiabatic index. This step is crucial as it incorporates the microphysics (EoS) into the dynamical description.

Pulsation Equation and Boundary Conditions

By substituting the expression for the pressure perturbation $ \delta p $ from (B10) into the perturbed equation of motion (B8) and after involving the use of the background TOV equations (6) to eliminate derivatives of the metric potentials, one arrives at the final pulsation equation:

Two key simplifications are required to cast it into the standard Sturm-Liouville form: multiplying the entire equation (B11) by ${\rm e}^{\nu_0 + 2\phi_0}$ to simplify the differential operator and introducing the new variable $\chi(r) \equiv r^2 {\rm e}^{-\phi_0(r)} \xi(r)$ to simplify the expression for the Lagrangian perturbation. Following manipulation, the equation takes the canonical Sturm-Liouville form:

where the coefficient functions are defined as

This formulation demonstrates the self-adjoint nature of the eigenvalue problem, which guarantees real eigenvalues $ \omega^2 $ as well as orthogonal and complete eigenfunctions.

Solutions to the pulsation equation must satisfy the boundary conditions imposed by physical regularity and spherical symmetry:

● At the center of the star ($ r = 0 $), all physical quantities must be finite. The radial displacement ξ must vanish at the origin to avoid a directional singularity, that is, $ \xi(0) = 0 $. A power series analysis shows that the leading-order of the displacement is linear: $ \xi \propto r $. Thus, owing to the factor of $ r^2 $, the scaling variable $\chi(r) = r^2 {\rm e}^{-\phi_0(r)} \xi(r)$ exhibits an asymptotic behavior $ \chi \propto r^3 $ near the origin. Therefore, the inner boundary condition is

In practice, numerical integration is started at a small but finite radius $ r = r_0 $ with the asymptotic form $ \chi(r_0) = \chi_0 r_0^3 $, where $ \chi_0 $ is a constant to be determined.

● At the surface of the star, the appropriate physical boundary condition is that the Lagrangian perturbation of the pressure must vanish:

Its relation to the Eulerian perturbation is

Therefore, in numerical calculations, we require the right side of B to be zero on the outer boundary.

The pulsation equation in Sturm-Liouville form, together with the boundary conditions (B13) and (B14), defines a characteristic value problem for the squared oscillation frequency $ \omega^2 $.

Stability Condition and Variational Principle

The stability of an equilibrium configuration is determined by the temporal evolution of the radial perturbations. From the ansatz $\xi(r, t) = \xi(r) {\rm e}^{{\rm i}\omega t}$, real values of ω (i.e., $ \omega^2 \gt 0 $) correspond to stable oscillations, whereas an imaginary ω ($ \omega^2 \lt 0 $) leads to exponentially growing or damping modes, indicating dynamical instability.

The self-adjoint nature of the pulsation equation (B11) permits the formulation of a variational principle for the characteristic frequencies $ \omega^2 $. Multiplication of (B11) by $ r^2 \xi {\rm e}^{(-\phi_0)} $ and integration over the star volume, followed by integration by parts and application of the boundary conditions, yield the integral expression

where the right-hand side "Numerator" is determined by the integral:

Equation (B16) expresses $ \omega^2 $ as the ratio of two integrals, which is known as the Rayleigh quotient:

Here, $ \omega^2[\xi] $ is a functional, the value of which depends on the choice of the trial function ξ.

The orthogonality of the eigenfunctions, as guaranteed by the self-adjoint Sturm-Liouville form, plays a crucial role in the variational principle. Any admissible trial function $ \widetilde{\xi} $ can be expanded as a series in the complete set of orthogonal eigenfunctions $ \xi^{(n)} $:

where $ c_n $ denotes the weight coefficients. The functional $ \widetilde{\omega }^2\left[\widetilde{\xi } \right] $ corresponding to the trial function is a weighted average of the true eigenvalues $ \omega_n^2 $:

where $ \langle \cdot, \cdot \rangle $ denotes the inner product with the corresponding weight functions. The Rayleigh quotient is stationary when $ \widetilde{\xi} $ is an eigenfunction (see [75] for further details), and it attains its minimum value $ \omega_0^2 $ when $ \widetilde{\xi} $ is proportional to the fundamental mode $ \xi^{(0)} $.

Therefore, for any trial function $ \widetilde{\xi} $ satisfying the boundary conditions, the computed value $ \widetilde{\omega}^2 $ is an upper bound of the true fundamental mode eigenvalue:

This formulation provides a powerful method for determining dynamical stability. A sufficient condition for dynamical instability is that any trial function $ \widetilde{\xi} $ exists for which the "Numerator" $ \leq 0 $. As the "Denominator" is always positive, this implies that $ \widetilde{\omega}^2 \leq 0 $, and consequently, $ \omega_0^2 \leq 0 $.

In his 1964 paper, Chandrasekhar rigorously established the pulsation equation for radial oscillations in single-fluid stars based on radial perturbations and variational principles [2], and it has been widely applied in many works, such as [101]. For the necessity of unifying notation, we briefly review the main steps and key equations in the derivation of this equation:

Linearization of Field Equations

An unperturbed, static, spherically symmetric equilibrium configuration satisfying the TOV equations is assumed. Its behavior after being subjected to infinitesimal radial perturbations can be used to determine its stability. We define ${\rm e}^{2\nu \left(r \right)}\equiv \left(1-\dfrac{2Gm\left(r \right)}{rc^2} \right) ^{-1}$. Then, the metric and fluid variables are expressed as

The Einstein field equations and Einstein-Euler equations are perturbed separately. Neglecting higher-order terms, we obtain the corresponding system of linearized radial perturbation equations:

Lagrangian Displacement and Reformulating Perturbations

The Lagrangian displacement ξ is defined such that the radial velocity satisfies $ v=\dfrac{{\rm d}r}{{\rm d}t}=\dfrac{\partial \xi}{\partial t} $, that is, $u^r={\rm e}^{-\phi_0}\dfrac{{\rm d}r}{{\rm d}t}= {\rm e}^{-\phi_0}\dfrac{\partial \xi}{\partial t}$. Then, using the first three equations of the linearized perturbation equations (B1), (B2), and (B3), all perturbations of the metric and fluid variables except for $ \delta p $ are expressed in terms of ξ and equilibrium quantities:

Notably, the perturbation quantities and Lagrangian displacement are all functions of $ r $ and $ t $ (e.g., $ \xi(r, t) $, $ \delta p(r, t) $). We now assume that all perturbation quantities $ F $ possess a time-harmonic dependence of the form $ F(r, t) = F(r) e^{i\omega t} $. Based on this separation of variables and by combining (B4) with (B5), (B6), and (B7), we obtain

where the perturbation equations are thus reduced from complex partial differential equations to a relatively simpler eigenvalue problem governed by an ordinary differential equation, with the perturbation terms now consisting of only the pressure variation.

Baryon Number Conservation

An additional relation is required to connect the pressure perturbation $ \delta p $ to the other variables. This is provided by imposing adiabaticity and baryon number conservation, which requires $ \nabla_\mu (n u^\mu) = 0 $. Defining the perturbation of $ n $ as $ n = n_0(r) + \delta n(r, t) $ and combining it with ξ defined previously leads to the form of $ \delta n $ constituted by the Lagrangian displacement and unperturbed terms:

Moreover, the first law of thermodynamics yields ${\rm d}\varepsilon = \dfrac{\varepsilon + p}{n} {\rm d}n$; thus, the baryon number density is now related to $ p $ or ε. This finally allows the pressure perturbation to be written solely in terms of the Lagrangian displacement:

Here, γ is defined as $\gamma =-\dfrac{n_0}{p_0}\dfrac{{\rm d}p_0}{{\rm d}n_0}$, which is the adiabatic index. This step is crucial as it incorporates the microphysics (EoS) into the dynamical description.

Pulsation Equation and Boundary Conditions

By substituting the expression for the pressure perturbation $ \delta p $ from (B10) into the perturbed equation of motion (B8) and after involving the use of the background TOV equations (6) to eliminate derivatives of the metric potentials, one arrives at the final pulsation equation:

Two key simplifications are required to cast it into the standard Sturm-Liouville form: multiplying the entire equation (B11) by ${\rm e}^{\nu_0 + 2\phi_0}$ to simplify the differential operator and introducing the new variable $\chi(r) \equiv r^2 {\rm e}^{-\phi_0(r)} \xi(r)$ to simplify the expression for the Lagrangian perturbation. Following manipulation, the equation takes the canonical Sturm-Liouville form:

where the coefficient functions are defined as

This formulation demonstrates the self-adjoint nature of the eigenvalue problem, which guarantees real eigenvalues $ \omega^2 $ as well as orthogonal and complete eigenfunctions.

Solutions to the pulsation equation must satisfy the boundary conditions imposed by physical regularity and spherical symmetry:

● At the center of the star ($ r = 0 $), all physical quantities must be finite. The radial displacement ξ must vanish at the origin to avoid a directional singularity, that is, $ \xi(0) = 0 $. A power series analysis shows that the leading-order of the displacement is linear: $ \xi \propto r $. Thus, owing to the factor of $ r^2 $, the scaling variable $\chi(r) = r^2 {\rm e}^{-\phi_0(r)} \xi(r)$ exhibits an asymptotic behavior $ \chi \propto r^3 $ near the origin. Therefore, the inner boundary condition is

In practice, numerical integration is started at a small but finite radius $ r = r_0 $ with the asymptotic form $ \chi(r_0) = \chi_0 r_0^3 $, where $ \chi_0 $ is a constant to be determined.

● At the surface of the star, the appropriate physical boundary condition is that the Lagrangian perturbation of the pressure must vanish:

Its relation to the Eulerian perturbation is

Therefore, in numerical calculations, we require the right side of B to be zero on the outer boundary.

The pulsation equation in Sturm-Liouville form, together with the boundary conditions (B13) and (B14), defines a characteristic value problem for the squared oscillation frequency $ \omega^2 $.

Stability Condition and Variational Principle

The stability of an equilibrium configuration is determined by the temporal evolution of the radial perturbations. From the ansatz $\xi(r, t) = \xi(r) {\rm e}^{{\rm i}\omega t}$, real values of ω (i.e., $ \omega^2 \gt 0 $) correspond to stable oscillations, whereas an imaginary ω ($ \omega^2 \lt 0 $) leads to exponentially growing or damping modes, indicating dynamical instability.

The self-adjoint nature of the pulsation equation (B11) permits the formulation of a variational principle for the characteristic frequencies $ \omega^2 $. Multiplication of (B11) by $ r^2 \xi {\rm e}^{(-\phi_0)} $ and integration over the star volume, followed by integration by parts and application of the boundary conditions, yield the integral expression

where the right-hand side "Numerator" is determined by the integral:

Equation (B16) expresses $ \omega^2 $ as the ratio of two integrals, which is known as the Rayleigh quotient:

Here, $ \omega^2[\xi] $ is a functional, the value of which depends on the choice of the trial function ξ.

The orthogonality of the eigenfunctions, as guaranteed by the self-adjoint Sturm-Liouville form, plays a crucial role in the variational principle. Any admissible trial function $ \widetilde{\xi} $ can be expanded as a series in the complete set of orthogonal eigenfunctions $ \xi^{(n)} $:

where $ c_n $ denotes the weight coefficients. The functional $ \widetilde{\omega }^2\left[\widetilde{\xi } \right] $ corresponding to the trial function is a weighted average of the true eigenvalues $ \omega_n^2 $:

where $ \langle \cdot, \cdot \rangle $ denotes the inner product with the corresponding weight functions. The Rayleigh quotient is stationary when $ \widetilde{\xi} $ is an eigenfunction (see [75] for further details), and it attains its minimum value $ \omega_0^2 $ when $ \widetilde{\xi} $ is proportional to the fundamental mode $ \xi^{(0)} $.

Therefore, for any trial function $ \widetilde{\xi} $ satisfying the boundary conditions, the computed value $ \widetilde{\omega}^2 $ is an upper bound of the true fundamental mode eigenvalue:

This formulation provides a powerful method for determining dynamical stability. A sufficient condition for dynamical instability is that any trial function $ \widetilde{\xi} $ exists for which the "Numerator" $ \leq 0 $. As the "Denominator" is always positive, this implies that $ \widetilde{\omega}^2 \leq 0 $, and consequently, $ \omega_0^2 \leq 0 $.

In his 1964 paper, Chandrasekhar rigorously established the pulsation equation for radial oscillations in single-fluid stars based on radial perturbations and variational principles [2], and it has been widely applied in many works, such as [101]. For the necessity of unifying notation, we briefly review the main steps and key equations in the derivation of this equation:

Linearization of Field Equations

An unperturbed, static, spherically symmetric equilibrium configuration satisfying the TOV equations is assumed. Its behavior after being subjected to infinitesimal radial perturbations can be used to determine its stability. We define ${\rm e}^{2\nu \left(r \right)}\equiv \left(1-\dfrac{2Gm\left(r \right)}{rc^2} \right) ^{-1}$. Then, the metric and fluid variables are expressed as

The Einstein field equations and Einstein-Euler equations are perturbed separately. Neglecting higher-order terms, we obtain the corresponding system of linearized radial perturbation equations:

Lagrangian Displacement and Reformulating Perturbations

The Lagrangian displacement ξ is defined such that the radial velocity satisfies $ v=\dfrac{{\rm d}r}{{\rm d}t}=\dfrac{\partial \xi}{\partial t} $, that is, $u^r={\rm e}^{-\phi_0}\dfrac{{\rm d}r}{{\rm d}t}= {\rm e}^{-\phi_0}\dfrac{\partial \xi}{\partial t}$. Then, using the first three equations of the linearized perturbation equations (B1), (B2), and (B3), all perturbations of the metric and fluid variables except for $ \delta p $ are expressed in terms of ξ and equilibrium quantities:

Notably, the perturbation quantities and Lagrangian displacement are all functions of $ r $ and $ t $ (e.g., $ \xi(r, t) $, $ \delta p(r, t) $). We now assume that all perturbation quantities $ F $ possess a time-harmonic dependence of the form $ F(r, t) = F(r) e^{i\omega t} $. Based on this separation of variables and by combining (B4) with (B5), (B6), and (B7), we obtain

where the perturbation equations are thus reduced from complex partial differential equations to a relatively simpler eigenvalue problem governed by an ordinary differential equation, with the perturbation terms now consisting of only the pressure variation.

Baryon Number Conservation

An additional relation is required to connect the pressure perturbation $ \delta p $ to the other variables. This is provided by imposing adiabaticity and baryon number conservation, which requires $ \nabla_\mu (n u^\mu) = 0 $. Defining the perturbation of $ n $ as $ n = n_0(r) + \delta n(r, t) $ and combining it with ξ defined previously leads to the form of $ \delta n $ constituted by the Lagrangian displacement and unperturbed terms:

Moreover, the first law of thermodynamics yields ${\rm d}\varepsilon = \dfrac{\varepsilon + p}{n} {\rm d}n$; thus, the baryon number density is now related to $ p $ or ε. This finally allows the pressure perturbation to be written solely in terms of the Lagrangian displacement:

Here, γ is defined as $\gamma =-\dfrac{n_0}{p_0}\dfrac{{\rm d}p_0}{{\rm d}n_0}$, which is the adiabatic index. This step is crucial as it incorporates the microphysics (EoS) into the dynamical description.

Pulsation Equation and Boundary Conditions

By substituting the expression for the pressure perturbation $ \delta p $ from (B10) into the perturbed equation of motion (B8) and after involving the use of the background TOV equations (6) to eliminate derivatives of the metric potentials, one arrives at the final pulsation equation:

Two key simplifications are required to cast it into the standard Sturm-Liouville form: multiplying the entire equation (B11) by ${\rm e}^{\nu_0 + 2\phi_0}$ to simplify the differential operator and introducing the new variable $\chi(r) \equiv r^2 {\rm e}^{-\phi_0(r)} \xi(r)$ to simplify the expression for the Lagrangian perturbation. Following manipulation, the equation takes the canonical Sturm-Liouville form:

where the coefficient functions are defined as

This formulation demonstrates the self-adjoint nature of the eigenvalue problem, which guarantees real eigenvalues $ \omega^2 $ as well as orthogonal and complete eigenfunctions.

Solutions to the pulsation equation must satisfy the boundary conditions imposed by physical regularity and spherical symmetry:

● At the center of the star ($ r = 0 $), all physical quantities must be finite. The radial displacement ξ must vanish at the origin to avoid a directional singularity, that is, $ \xi(0) = 0 $. A power series analysis shows that the leading-order of the displacement is linear: $ \xi \propto r $. Thus, owing to the factor of $ r^2 $, the scaling variable $\chi(r) = r^2 {\rm e}^{-\phi_0(r)} \xi(r)$ exhibits an asymptotic behavior $ \chi \propto r^3 $ near the origin. Therefore, the inner boundary condition is

In practice, numerical integration is started at a small but finite radius $ r = r_0 $ with the asymptotic form $ \chi(r_0) = \chi_0 r_0^3 $, where $ \chi_0 $ is a constant to be determined.

● At the surface of the star, the appropriate physical boundary condition is that the Lagrangian perturbation of the pressure must vanish:

Its relation to the Eulerian perturbation is

Therefore, in numerical calculations, we require the right side of B to be zero on the outer boundary.

The pulsation equation in Sturm-Liouville form, together with the boundary conditions (B13) and (B14), defines a characteristic value problem for the squared oscillation frequency $ \omega^2 $.

Stability Condition and Variational Principle

The stability of an equilibrium configuration is determined by the temporal evolution of the radial perturbations. From the ansatz $\xi(r, t) = \xi(r) {\rm e}^{{\rm i}\omega t}$, real values of ω (i.e., $ \omega^2 \gt 0 $) correspond to stable oscillations, whereas an imaginary ω ($ \omega^2 \lt 0 $) leads to exponentially growing or damping modes, indicating dynamical instability.

The self-adjoint nature of the pulsation equation (B11) permits the formulation of a variational principle for the characteristic frequencies $ \omega^2 $. Multiplication of (B11) by $ r^2 \xi {\rm e}^{(-\phi_0)} $ and integration over the star volume, followed by integration by parts and application of the boundary conditions, yield the integral expression

where the right-hand side "Numerator" is determined by the integral:

Equation (B16) expresses $ \omega^2 $ as the ratio of two integrals, which is known as the Rayleigh quotient:

Here, $ \omega^2[\xi] $ is a functional, the value of which depends on the choice of the trial function ξ.

The orthogonality of the eigenfunctions, as guaranteed by the self-adjoint Sturm-Liouville form, plays a crucial role in the variational principle. Any admissible trial function $ \widetilde{\xi} $ can be expanded as a series in the complete set of orthogonal eigenfunctions $ \xi^{(n)} $:

where $ c_n $ denotes the weight coefficients. The functional $ \widetilde{\omega }^2\left[\widetilde{\xi } \right] $ corresponding to the trial function is a weighted average of the true eigenvalues $ \omega_n^2 $:

where $ \langle \cdot, \cdot \rangle $ denotes the inner product with the corresponding weight functions. The Rayleigh quotient is stationary when $ \widetilde{\xi} $ is an eigenfunction (see [75] for further details), and it attains its minimum value $ \omega_0^2 $ when $ \widetilde{\xi} $ is proportional to the fundamental mode $ \xi^{(0)} $.

Therefore, for any trial function $ \widetilde{\xi} $ satisfying the boundary conditions, the computed value $ \widetilde{\omega}^2 $ is an upper bound of the true fundamental mode eigenvalue:

This formulation provides a powerful method for determining dynamical stability. A sufficient condition for dynamical instability is that any trial function $ \widetilde{\xi} $ exists for which the "Numerator" $ \leq 0 $. As the "Denominator" is always positive, this implies that $ \widetilde{\omega}^2 \leq 0 $, and consequently, $ \omega_0^2 \leq 0 $.

For a compact star composed of multiple ideal fluids in the absence of any non-gravitational coupling, the pulsation equations for a single fluid can be generalized to the multi-fluid case. The oscillation equations for such a system have previously been derived in the literature and applied to other mixed models (see, e.g., [18, 19, 84]). This section briefly outlines the main steps in deriving the equations using the same process as that in the references.

Linearized Field Equations and Energy-Momentum Tensor

Consider a spherically symmetric system composed of ${\cal{N}}$ ideal fluids. The total energy-momentum tensor is the sum of the decoupled contributions from each fluid:

For a spherically symmetric background, assuming that the EoS for each fluid depends only on its own energy density (i.e., the fluids only interact through gravity): $p_I = p_I(\varepsilon_I)$, the radial perturbation variables are set as

and the Lagrangian displacement $\xi_I$ of the $I$-th fluid satisfies $ v_I = \dot{\xi}_I $. Similar to the case of a single fluid, linearizing the Einstein field equations and energy-momentum conservation equations would yield a set of coupled perturbation equations.

Lagrangian Displacement and Expression of Perturbations

We introduce the harmonic time dependence assumption $\xi_I(t,r) = \xi_I(r) {\rm e}^{{\rm i}\omega t}$, and similarly for other perturbation variables. Similar to the case in the previous section, the metric perturbations $\delta\nu$, $\delta\phi$ can be expressed in terms of the time-independent $\xi_I$:

The energy density perturbation of each fluid can be expressed as

Particle Number Conservation and Pressure Perturbation

Similar to the case of a single fluid, particle number perturbations are constructed through covariant conservation, as follows:

For each fluid, the adiabatic index $\gamma_I =-\dfrac{n_{I0}}{p_{I0}}\dfrac{{\rm d}p_{I0}}{{\rm d}n_{I0}}=\left(\dfrac{p_{I0}+\varepsilon _{I0}}{p_{I0}} \right) \dfrac{{\rm d}p_{I0}}{{\rm d}\varepsilon _{I0}}$ is introduced. The pressure perturbation can then be expressed as

These perturbation expressions couple the fluid variables and displacements $\xi_I$ of all fluids, reflecting the gravitational interaction between the multiple fluids.

Multi-fluid Pulsation Equations and Boundary Conditions

Substituting the aforementioned expressions into the linearized equations and introducing the variable $\chi _I(r)=r^2{\rm e}^{-\phi _0}\xi _I(r)$ (where $\phi _c=\phi _0\left(r=0 \right) $) yield a set of coupled second-order ordinary differential equations:

where the coefficients are

The coupling terms are explicitly expressed as

where

The boundary conditions are as follows:

● At the star center ($r=0$), $\chi_I \propto r^3$ for each fluid.

● At the boundary of each fluid $r = R_I$, the Lagrangian pressure perturbation must vanish, that is,

Here, $ p_{I0} $ tends to 0 at the corresponding outer boundary of the fluid. Therefore, as the other quantities ($ \xi _I $ and unperturbed quantities) in the final term must be finite, the contribution of this term to the boundary constraint is ignored, and thus, the boundary condition reduces to

This outer boundary condition is also consistent with those in [53, 75]. In [18], the corresponding term of ${\rm e}^{\phi _0}$ is ${\rm e}^{2\phi _0}$ after unifying the convention, which appears to be a typo. However, this slight discrepancy merely introduces a negligible scaling factor in the numerical results and does not affect the stability boundary.

Stability Condition and Variational Principle

The stability criterion from the single-fluid analysis carries over to the two-fluid case: for perturbations of the form $\xi_I(r, t) = \xi_I(r) {\rm e}^{{\rm i}\omega t}$, the configuration is stable for real ω and dynamically unstable for imaginary ω.

Although the multi-fluid pulsation equations cannot be written in a standard Sturm-Liouville form, the coupled Rayleigh quotient can still be provided:

which could be constructed by the integration of the pulsation equation, including various complex coupling terms. The extreme value of its trial functions are the eigenvalues, and the minimum value is the square of the fundamental frequency. If a set of trial functions $\{\widetilde{\xi}_I\}$ exists such that $\omega^2 \lt 0$, the system is unstable against that perturbation mode.

For a compact star composed of multiple ideal fluids in the absence of any non-gravitational coupling, the pulsation equations for a single fluid can be generalized to the multi-fluid case. The oscillation equations for such a system have previously been derived in the literature and applied to other mixed models (see, e.g., [18, 19, 84]). This section briefly outlines the main steps in deriving the equations using the same process as that in the references.

Linearized Field Equations and Energy-Momentum Tensor

Consider a spherically symmetric system composed of ${\cal{N}}$ ideal fluids. The total energy-momentum tensor is the sum of the decoupled contributions from each fluid:

For a spherically symmetric background, assuming that the EoS for each fluid depends only on its own energy density (i.e., the fluids only interact through gravity): $p_I = p_I(\varepsilon_I)$, the radial perturbation variables are set as

and the Lagrangian displacement $\xi_I$ of the $I$-th fluid satisfies $ v_I = \dot{\xi}_I $. Similar to the case of a single fluid, linearizing the Einstein field equations and energy-momentum conservation equations would yield a set of coupled perturbation equations.

Lagrangian Displacement and Expression of Perturbations

We introduce the harmonic time dependence assumption $\xi_I(t,r) = \xi_I(r) {\rm e}^{{\rm i}\omega t}$, and similarly for other perturbation variables. Similar to the case in the previous section, the metric perturbations $\delta\nu$, $\delta\phi$ can be expressed in terms of the time-independent $\xi_I$:

The energy density perturbation of each fluid can be expressed as

Particle Number Conservation and Pressure Perturbation

Similar to the case of a single fluid, particle number perturbations are constructed through covariant conservation, as follows:

For each fluid, the adiabatic index $\gamma_I =-\dfrac{n_{I0}}{p_{I0}}\dfrac{{\rm d}p_{I0}}{{\rm d}n_{I0}}=\left(\dfrac{p_{I0}+\varepsilon _{I0}}{p_{I0}} \right) \dfrac{{\rm d}p_{I0}}{{\rm d}\varepsilon _{I0}}$ is introduced. The pressure perturbation can then be expressed as

These perturbation expressions couple the fluid variables and displacements $\xi_I$ of all fluids, reflecting the gravitational interaction between the multiple fluids.

Multi-fluid Pulsation Equations and Boundary Conditions

Substituting the aforementioned expressions into the linearized equations and introducing the variable $\chi _I(r)=r^2{\rm e}^{-\phi _0}\xi _I(r)$ (where $\phi _c=\phi _0\left(r=0 \right) $) yield a set of coupled second-order ordinary differential equations:

where the coefficients are

The coupling terms are explicitly expressed as

where

The boundary conditions are as follows:

● At the star center ($r=0$), $\chi_I \propto r^3$ for each fluid.

● At the boundary of each fluid $r = R_I$, the Lagrangian pressure perturbation must vanish, that is,

Here, $ p_{I0} $ tends to 0 at the corresponding outer boundary of the fluid. Therefore, as the other quantities ($ \xi _I $ and unperturbed quantities) in the final term must be finite, the contribution of this term to the boundary constraint is ignored, and thus, the boundary condition reduces to

This outer boundary condition is also consistent with those in [53, 75]. In [18], the corresponding term of ${\rm e}^{\phi _0}$ is ${\rm e}^{2\phi _0}$ after unifying the convention, which appears to be a typo. However, this slight discrepancy merely introduces a negligible scaling factor in the numerical results and does not affect the stability boundary.

Stability Condition and Variational Principle

The stability criterion from the single-fluid analysis carries over to the two-fluid case: for perturbations of the form $\xi_I(r, t) = \xi_I(r) {\rm e}^{{\rm i}\omega t}$, the configuration is stable for real ω and dynamically unstable for imaginary ω.

Although the multi-fluid pulsation equations cannot be written in a standard Sturm-Liouville form, the coupled Rayleigh quotient can still be provided:

which could be constructed by the integration of the pulsation equation, including various complex coupling terms. The extreme value of its trial functions are the eigenvalues, and the minimum value is the square of the fundamental frequency. If a set of trial functions $\{\widetilde{\xi}_I\}$ exists such that $\omega^2 \lt 0$, the system is unstable against that perturbation mode.

For a compact star composed of multiple ideal fluids in the absence of any non-gravitational coupling, the pulsation equations for a single fluid can be generalized to the multi-fluid case. The oscillation equations for such a system have previously been derived in the literature and applied to other mixed models (see, e.g., [18, 19, 84]). This section briefly outlines the main steps in deriving the equations using the same process as that in the references.

Linearized Field Equations and Energy-Momentum Tensor

Consider a spherically symmetric system composed of ${\cal{N}}$ ideal fluids. The total energy-momentum tensor is the sum of the decoupled contributions from each fluid:

For a spherically symmetric background, assuming that the EoS for each fluid depends only on its own energy density (i.e., the fluids only interact through gravity): $p_I = p_I(\varepsilon_I)$, the radial perturbation variables are set as

and the Lagrangian displacement $\xi_I$ of the $I$-th fluid satisfies $ v_I = \dot{\xi}_I $. Similar to the case of a single fluid, linearizing the Einstein field equations and energy-momentum conservation equations would yield a set of coupled perturbation equations.

Lagrangian Displacement and Expression of Perturbations

We introduce the harmonic time dependence assumption $\xi_I(t,r) = \xi_I(r) {\rm e}^{{\rm i}\omega t}$, and similarly for other perturbation variables. Similar to the case in the previous section, the metric perturbations $\delta\nu$, $\delta\phi$ can be expressed in terms of the time-independent $\xi_I$:

The energy density perturbation of each fluid can be expressed as

Particle Number Conservation and Pressure Perturbation

Similar to the case of a single fluid, particle number perturbations are constructed through covariant conservation, as follows:

For each fluid, the adiabatic index $\gamma_I =-\dfrac{n_{I0}}{p_{I0}}\dfrac{{\rm d}p_{I0}}{{\rm d}n_{I0}}=\left(\dfrac{p_{I0}+\varepsilon _{I0}}{p_{I0}} \right) \dfrac{{\rm d}p_{I0}}{{\rm d}\varepsilon _{I0}}$ is introduced. The pressure perturbation can then be expressed as

These perturbation expressions couple the fluid variables and displacements $\xi_I$ of all fluids, reflecting the gravitational interaction between the multiple fluids.

Multi-fluid Pulsation Equations and Boundary Conditions

Substituting the aforementioned expressions into the linearized equations and introducing the variable $\chi _I(r)=r^2{\rm e}^{-\phi _0}\xi _I(r)$ (where $\phi _c=\phi _0\left(r=0 \right) $) yield a set of coupled second-order ordinary differential equations:

where the coefficients are

The coupling terms are explicitly expressed as

where

The boundary conditions are as follows:

● At the star center ($r=0$), $\chi_I \propto r^3$ for each fluid.

● At the boundary of each fluid $r = R_I$, the Lagrangian pressure perturbation must vanish, that is,

Here, $ p_{I0} $ tends to 0 at the corresponding outer boundary of the fluid. Therefore, as the other quantities ($ \xi _I $ and unperturbed quantities) in the final term must be finite, the contribution of this term to the boundary constraint is ignored, and thus, the boundary condition reduces to

This outer boundary condition is also consistent with those in [53, 75]. In [18], the corresponding term of ${\rm e}^{\phi _0}$ is ${\rm e}^{2\phi _0}$ after unifying the convention, which appears to be a typo. However, this slight discrepancy merely introduces a negligible scaling factor in the numerical results and does not affect the stability boundary.

Stability Condition and Variational Principle

The stability criterion from the single-fluid analysis carries over to the two-fluid case: for perturbations of the form $\xi_I(r, t) = \xi_I(r) {\rm e}^{{\rm i}\omega t}$, the configuration is stable for real ω and dynamically unstable for imaginary ω.

Although the multi-fluid pulsation equations cannot be written in a standard Sturm-Liouville form, the coupled Rayleigh quotient can still be provided:

which could be constructed by the integration of the pulsation equation, including various complex coupling terms. The extreme value of its trial functions are the eigenvalues, and the minimum value is the square of the fundamental frequency. If a set of trial functions $\{\widetilde{\xi}_I\}$ exists such that $\omega^2 \lt 0$, the system is unstable against that perturbation mode.