Dark Matter Admixed Neutron Stars

• Dark matter admixed neutron stars are two-fluid systems where baryonic matter and dark matter interact solely via gravity, allowing configurations from dark cores to extended halos.
• They exhibit a range of equations of state and structural properties that modify mass-radius relations, stability thresholds, and electromagnetic as well as gravitational observables.
• The two-fluid formalism enables independent analyses of rotation and oscillations, offering insights into multimessenger signals and constraints on dark matter microphysics in compact stars.

Searching arXiv for recent and foundational papers on dark matter admixed neutron stars to ground the article. Dark matter admixed neutron stars (DANSs) are compact stars whose total mass-energy contains both ordinary baryonic or nuclear matter and a non-negligible dark matter component. In the literature summarized here, they are modeled primarily as two-fluid relativistic stars in which the visible neutron-star matter and the dark component are treated as separate perfect or degenerate fluids that interact only through gravity, sharing the same spacetime but generally retaining distinct pressures, energy densities, radii, and in some treatments distinct rotational states. This framework permits configurations ranging from modest dark cores embedded within otherwise ordinary neutron stars to dark-matter-dominated objects with a small normal-matter core inside an extended dark halo, and it has been used to study equilibrium structure, stability, oscillations, cooling thresholds, rotation, tidal response, binary initial data, X-ray observables, and observational bounds on dark matter content (Leung et al., 2011).

1. Foundational definition and two-fluid formalism

The foundational relativistic description of DANSs treats them as two gravitationally coupled fluids: ordinary nuclear matter (NM), or more generally baryonic matter (BM), and dark matter (DM). In the early general-relativistic treatment, the system is formulated through a master function

$\Lambda(n^2,p^2,x^2),$

with

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$

where the NM and DM currents are separately conserved. Under the assumption that NM and DM interact only through gravity, the master function is taken to be separable,

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$

with no dependence on $x^2$, so entrainment-type coupling vanishes in that model (Leung et al., 2011).

In static spherical symmetry, the spacetime is written as

$$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$

and the coupled equilibrium equations generalize the Tolman–Oppenheimer–Volkoff system. In one formulation, the structure equations are

$$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$

$$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$

$$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$

$$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$

while many later treatments adopt the equivalent two-fluid TOV form in which each fluid satisfies its own hydrostatic equation in the common gravitational potential (Leung et al., 2011). A representative form is

$$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$

with

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$0

so that the fluids are dynamically distinct but gravitationally inseparable (Xiang et al., 2013).

A persistent feature of this formalism is the existence of two characteristic radii. The baryonic or nuclear surface is defined where the visible fluid pressure vanishes, while the dark surface is defined by the vanishing of the dark pressure. Thus a DANS may be core-dominated, with $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$1, or halo-dominated, with $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$2. Because electromagnetic emission originates at the NM or BM surface, the observed radius may be the visible core radius rather than the full gravitational extent of the object (Xiang et al., 2013).

2. Equations of state and dark-sector microphysics

The DANS literature is distinguished less by the two-fluid gravitational framework, which is broadly shared, than by the diversity of equations of state assumed for the baryonic and dark sectors. For ordinary matter, the cited studies use realistic neutron-star EOSs including APR (Leung et al., 2011), relativistic mean-field models such as SLC, SLCd, and IU-FSU (Xiang et al., 2013), microscopic chiral-interaction matter through the BL EOS derived in Brueckner-Hartree-Fock theory (Scordino et al., 2024), SLy (Rüter et al., 2023, Kain, 2021), NL3$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$3L55 (Shawqi et al., 25 Aug 2025), DDME2 (Zou et al., 11 Apr 2026), unified NL3 from CompOSE (Zhang et al., 25 Jun 2026), and covariant density functional families spanning PKA1, PKO1–3, PKDD, DD-ME2, DD-MEX, DD-LZ1, NL3, GM1, TM1e, and NL3$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$4 (Hu et al., 29 Dec 2025).

The dark sector is treated in several distinct ways. A minimal and influential choice is a non-self-annihilating ideal relativistic Fermi gas of fermions with mass $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$5 GeV (Leung et al., 2011). Closely related studies vary the fermion mass and introduce scalar attraction and vector repulsion through a relativistic mean-field-like Lagrangian,

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$6

with effective potential

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$7

so vector repulsion stiffens the DM EOS and scalar attraction softens it (Xiang et al., 2013).

Other studies use repulsively self-interacting fermionic dark matter in the compact-star parametrization

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$8

with weakly interacting and strongly interacting regimes specified by $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$9 and $$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$0, respectively (Zou et al., 11 Apr 2026). A separate first-principles direction models DM as a confining QCD-like sector based on lattice $$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$1-QCD, with the dark matter candidate taken to be the stable three-dark-quark baryon and the EOS interpolated and extrapolated from lattice data (Dengler et al., 25 Mar 2025).

Bosonic dark matter has also been explored. One rapidly rotating study treats DM as a self-interacting bosonic fluid derived from the complex-scalar Lagrangian

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$2

leading in the zero-temperature Bose-condensed limit to

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$3

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$4

(Cipriani et al., 25 Feb 2025). Another study compares self-interacting bosonic and ideal-fermion DM and shows that both can generate halo-induced departures in the mass–radius–tidal response relation at fixed stellar mass (Sun et al., 2023).

This breadth of microphysical input is central to the field. It implies that DANS phenomenology is not controlled by a single dark parameter but by particle mass, self-interaction strength, annihilation properties, and, in practice, the central abundance or enthalpy ratio of the two fluids.

3. Internal morphology: dark cores, dark halos, and composition parameters

A defining structural distinction in DANS research is the difference between dark-core and dark-halo configurations. Early general-relativistic calculations established that sufficiently DM-rich equilibria can consist of a small NM core of radius of order a few km surrounded by a DM halo of total size $$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$5 km (Leung et al., 2011). In the strongly DM-dominated example with composition parameter

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$6

set to $$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$7, the density profiles show NM concentrated at the center and DM forming the extended halo. For one representative configuration, the visible radius of the NM core is

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$8

while the mass enclosed within the NM core is

$$\Lambda(n^2,p^2)=\Lambda_{\rm NM}(n^2)+\Lambda_{\rm DM}(p^2),$$9

(Leung et al., 2011).

Later work sharpened the core–halo distinction. A DM halo forms when the DM distribution extends beyond the NM distribution, $x^2$0, and such halos are found to be more likely for low-mass DM, in the presence of repulsive DM self-interactions, and at lower central energy density (Xiang et al., 2013). In contrast, attractive self-interactions concentrate DM and suppress halo formation (Xiang et al., 2013). More recent analysis showed that the relative distribution is essentially determined by the ratio of the central enthalpy of the DM component to that of baryonic matter: $x^2$1 (Miao et al., 2022). At the critical configuration $x^2$2, the critical DM fraction satisfies the analytical estimate

$x^2$3

implying $x^2$4 at fixed self-interaction parameter and $x^2$5 in the strong self-interaction regime (Miao et al., 2022).

Several papers provide useful classification schemes. One crust study defines the halo mass fraction

$x^2$6

and labels configurations as dark core for $x^2$7, compact halo for $x^2$8, intermediate halo for $x^2$9, and diffuse halo for $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$0. In its fixed-$$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$1 sequences, the transition proceeds from diffuse halo for $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$2 MeV through compact halo for $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$3 MeV to dark core for $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$4 MeV (Zhang et al., 25 Jun 2026).

A consistent trend across the literature is that heavier DM tends to be more centrally concentrated, while lighter or more repulsive DM tends to spread outward and form halos. This suggests that the core–halo dichotomy is not a secondary detail but one of the primary observables linking dark-sector microphysics to neutron-star phenomenology.

4. Mass–radius relations, compactness, and maximum-mass systematics

In ordinary neutron-star theory, a fixed EOS produces a one-parameter mass–radius sequence. A central result of DANS research is that this property is lost once dark matter is admitted. Because the equilibrium depends on two EOSs and on parameters such as DM fraction and DM self-interaction strength, DANSs populate a family of mass–radius curves rather than a unique sequence (Xiang et al., 2013). This “spread of mass-radius relationships” implies that a measured mass and radius cannot be mapped to a unique nuclear EOS without additional assumptions about DM content (Xiang et al., 2013).

The earliest detailed DANS study found that increasing the DM fraction shifts the mass–radius relation away from that of ordinary neutron stars, decreases the maximum stable mass, and decreases the stellar radius. For the maximum stable-mass configuration, increasing $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$5 from 0 to 0.2 reduces $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$6 by about 35% and $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$7 by about 9% (Leung et al., 2011). A later two-fluid TOV analysis confirmed that for fixed DM particle mass, increasing the DM fraction lowers the maximum stable mass, reduces the stellar radius, and increases the spread of possible mass–radius curves (Xiang et al., 2013).

The effect of DM mass is more nuanced. In one study, heavy DM candidates with mass $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$8 GeV have little effect on the mass–radius relation, whereas light DM of order $$ds^2 = -e^{\nu(r)}dt^2 + e^{\lambda(r)}dr^2 + r^2(d\theta^2+\sin^2\theta\,d\phi^2),$$9 MeV produces major changes in both maximum mass and radius (Xiang et al., 2013). Another study using a realistic chiral nuclear EOS finds opposite trends depending on DM mass: for $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$0, DM reduces the maximum gravitational mass relative to the ordinary-star case, but for $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$1, DANSs form with a DM halo and the total gravitational maximum mass increases relative to the ordinary-star case (Scordino et al., 2024). The same work identifies a transition at $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$2: for $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$3, the maximum mass remains compatible with observed $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$4 neutron stars for any DM fraction, whereas for larger masses only restricted DM fractions are permitted (Scordino et al., 2024).

A related scaling appears in the self-gravitating pure-DM limit. For fermionic dark stars, one study quotes

$$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$5

(Xiang et al., 2013). A 2026 extension to extremely light DM emphasizes the corresponding general scaling

$$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$6

or, including self-interaction,

$$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$7

and argues that for $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$8 GeV the DANS becomes DM-dominated, with the NM core reduced to an embedded seed inside an enormous DM halo (Zou et al., 11 Apr 2026). For $$A^0_0 p' + B^0_0 n' + \frac{1}{2}(Bn + Ap)\nu' = 0,$$9, the paper reports

$$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$0

while the NM core remains near

$$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$1

(Zou et al., 11 Apr 2026). This suggests that the DANS concept can interpolate continuously between modestly admixed neutron stars and supermassive dark objects with embedded neutron-star cores.

Not all approaches use a strict two-fluid separation. Some mean-field studies incorporate DM directly into the core EOS through Higgs exchange or via fixed DM Fermi momentum. In those models, larger DM masses and larger DM Fermi momenta generally produce more compact stars with reduced radii, lower maximum masses, and smaller tidal deformabilities (Kumar et al., 2024, Lopes et al., 2024). A specific case study interprets XTE J1814-338, with inferred

$$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$2

as a DANS requiring a DM Fermi momentum $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$3 GeV, yielding $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$4 km and $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$5 in the nucleonic version of the model (Lopes et al., 2024). Because that work does not present a dedicated radial-stability analysis, its ultracompact configurations are best interpreted as phenomenological TOV solutions rather than fully established stable equilibria (Lopes et al., 2024).

5. Stability, oscillations, and crustal seismology

The stability of DANSs has been examined through several complementary methods. The original relativistic two-fluid study solved the radial oscillation eigenvalue problem and found that for the DM-dominated sequence with $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$6, the fundamental mode satisfies $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$7 along the stable branch, crosses zero at the maximum-mass configuration, and becomes negative beyond it. Thus the onset of instability coincides with the maximum-mass point, preserving the standard one-fluid intuition in the two-fluid setting (Leung et al., 2011).

A later comprehensive analysis emphasized that in a two-fluid star the stability boundary is not a single point but a critical curve in the central-pressure plane $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$8. It established stability in two ways: by solving the harmonic radial perturbation problem with $$C^0_0 p' + A^0_0 n' + \frac{1}{2}(An + Cp)\nu' = 0,$$9, and by applying a turning-point criterion along directions tangent to contours of conserved particle numbers. Along the stability boundary there exists a direction $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$0 such that

$$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$1

and the two methods agree numerically (Kain, 2021). A central conclusion is that there exist stable regions in parameter space where either the ordinary-matter or dark-matter central density exceeds its corresponding one-fluid critical value, so one-fluid stability criteria are insufficient for DANSs (Kain, 2021, Scordino et al., 2024).

Nonlinear stability has also been verified dynamically. In the first dynamical evolution of DANSs, the equations of motion for multiple perfect fluids were cast in conservation form and evolved using finite volume and high-resolution shock-capturing methods. Stable static DANS remained stable nonlinearly, while unstable ones collapsed to black holes. The same simulations recovered radial oscillation frequencies by Fourier transforming central-pressure time series and showed excellent agreement with linear pulsation theory (Gleason et al., 2022). That work also demonstrated a dynamical formation scenario in which a preexisting neutron star surrounded by a Gaussian dark-matter cloud evolves into a stable mixed star, including a case where a dark-matter core forms by $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$2 (Gleason et al., 2022).

The crust has recently been incorporated explicitly. In DANS models with unified NL3 baryonic matter and repulsively self-interacting fermionic DM, the crust thickness

$$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$3

is always smaller than in the corresponding pure neutron star at fixed total mass and fixed DM fraction $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$4. For $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$5, the maximum reduction is about 12% near $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$6 GeV; for $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$7, the maximum reduction is about 16% (Zhang et al., 25 Jun 2026). The effect is negligible for diffuse halos, where little DM lies inside the baryonic surface, and strongest for compact halos and dark cores (Zhang et al., 25 Jun 2026). The same paper derives the approximation

$$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$8

with $$\lambda' = \frac{1-e^\lambda}{r} - 8\pi r e^\lambda \Lambda,$$9 for NL3, and finds that the rational form reproduces numerical crust thicknesses typically below 0.5% error and at worst about 1% (Zhang et al., 25 Jun 2026).

This crustal compression raises torsional mode frequencies. In the relativistic Cowling approximation, the torsional eigenvalue problem is

$$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$0

with boundary conditions at the crust-core interface and surface (Zhang et al., 25 Jun 2026). For $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$1, the maximum frequency shifts relative to pure neutron stars are about 6.4% for the quadrupole fundamental $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$2 and about 16.4% for the first overtone $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$3 as $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$4 GeV (Zhang et al., 25 Jun 2026). Because electron screening softens the crust and lowers frequencies, while DM thins the crust and raises them, the two effects can be degenerate; the paper finds that this degeneracy can be broken in some regions of parameter space, especially for massive stars and the first overtone (Zhang et al., 25 Jun 2026).

6. Rotation, tidal response, binaries, and multimessenger observables

Rotation exposes the multi-fluid nature of DANSs particularly clearly. In the slow-rotation Hartle formalism generalized to multiple fluids, each fluid has its own angular velocity $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$5, four-velocity

$$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$6

and moment of inertia $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$7, while the total angular momentum is $$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$8 (Cronin et al., 2023). The first-order frame-dragging equation is sourced by the sum of all fluids, but each fluid contributes separately to the angular-momentum budget (Cronin et al., 2023). A key observational point is that electromagnetic measurements probe ordinary matter alone, whereas gravitational observables such as the tidal deformability depend on the total mass distribution. Accordingly, the observationally relevant dimensionless quantities are

$$\nu' = -\frac{1-e^\lambda}{r} + 8\pi r e^\lambda \Psi,$$9

not their total-fluid analogues (Cronin et al., 2023). Using these definitions, the paper shows that DANS can deviate significantly from standard single-fluid I–Love–Q relations (Cronin et al., 2023).

Rapid rotation has been studied in two distinct regimes. One axisymmetric two-fluid calculation extends the RNS code to rapidly rotating DANS in which BM is spun up by accretion and rotates rigidly, while DM remains torque-free and differentially rotates through frame dragging: $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$0 (Shawqi et al., 25 Aug 2025). For three representative halo types with non-rotating DM fraction $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$1, corresponding at $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$2 to $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$3 km, $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$4 km, and $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$5 km, rapid BM rotation reduces DM halo sizes when central densities are held fixed. At $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$6, the baryonic equatorial radius increases by about 5% and the polar radius decreases by about 4%, while the diffuse halo radius can shrink by about 80% and the cloud mass outside the baryonic surface by about 97%; the corresponding reductions are about 20% and 39% for the intermediate halo and about 5% and 37% for the compact halo (Shawqi et al., 25 Aug 2025). The same study emphasizes that local TOV-like component masses and global Komar-type component masses differ and should not be interpreted as separately measurable gravitational masses (Shawqi et al., 25 Aug 2025).

A complementary rapid-rotation study allows BM and DM to rotate independently, including co-rotating and counter-rotating states, in models with self-interacting bosonic DM (Cipriani et al., 25 Feb 2025). There the defining feature is independent fluid angular velocities,

$$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$7

with even opposite signs permitted (Cipriani et al., 25 Feb 2025). In the baseline case where BM rotates rapidly but DM is non-rotating, pure baryonic stars show the expected rotational increase of maximum mass by about 20% and radius by about 40%, while DM core configurations reduce both maximum mass and moment of inertia relative to no-DM stars (Cipriani et al., 25 Feb 2025). Turning on DM rotation reveals strong qualitative differences between co-rotation and counter-rotation. In one example the ratio of DM polar to equatorial radius drops from about 0.51 in co-rotation to about 0.34 in counter-rotation, and counter-rotating sequences become numerically difficult before the last stable configuration because of extreme deformation (Cipriani et al., 25 Feb 2025).

Tidal response has likewise emerged as a major diagnostic. One study proposes a new criterion: at fixed stellar mass, DANS can exhibit a negative correlation between normal-matter radius $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$8 and tidal deformability $$\frac{dp_N}{dr}= -\big[\mathcal{E}_N(r)+p_N(r)\big]\frac{d\nu}{dr},\qquad \frac{dp_D}{dr}= -\big[\mathcal{E}_D(r)+p_D(r)\big]\frac{d\nu}{dr},$$9, so that a star with smaller visible radius has larger $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$00. This does not occur for ordinary neutron stars in that analysis and arises because $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$01 depends on the larger of $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$02 and $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$03, while the observed radius tracks only $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$04 (Sun et al., 2023). For APR3 with self-interacting bosonic DM of mass $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$05 MeV, the paper reports that $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$06 can reach about 8000 for large DM fraction $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$07, compared with $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$08 in ordinary APR3 stars (Sun et al., 2023).

Binary DANS systems have also begun to be modeled self-consistently. Using an adapted SGRID code, the first constraint-solved quasi-equilibrium initial data for equal-mass binary neutron stars each containing dark matter were constructed in both dark-core and dark-halo scenarios (Rüter et al., 2023). The stars are treated as two non-interacting perfect fluids minimally coupled to gravity, with irrotational flow $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$09, and the data solve the coupled XCTS equations plus fluid equilibrium conditions (Rüter et al., 2023). DM modifies tidal response in a morphology-dependent way: halo configurations are much more easily distorted, and the DM component in the halo case reaches a mass-shedding diagnostic $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$10 already at relatively large separation; mass shedding occurs at separations between about $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$11 and $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$12, whereas in the dark-core case the baryonic fluid reaches mass shedding only around $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$13 in the representative sequence (Rüter et al., 2023).

7. Observational signatures, constraints, and open interpretations

The principal observational motivation for DANSs is that dark matter can decouple the gravitational and electromagnetic appearance of a compact star. A DANS may have a total mass in the neutron-star range and a total radius of order 10 km, yet only a smaller NM-emitting core contributes to thermal radiation. The original DANS paper therefore emphasized that such objects may be observed as extraordinarily small neutron stars incompatible with realistic nuclear matter models (Leung et al., 2011). A later two-fluid study illustrated this with a comparison at fixed total mass $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$14: the DANS had a smaller NM core radius of about 9 km inside a total radius $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$15 km and a larger redshift, roughly $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$16 versus $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$17 for the ordinary neutron star (Xiang et al., 2013).

Moment of inertia offers an additional discriminator. In the non-interacting two-fluid limit of the original relativistic model, one can write

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$18

and define the scaled quantity

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$19

(Leung et al., 2011). Ordinary neutron stars lie near the Bejger–Haensel universal fit in terms of compactness $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$20, but DANS do not; for DM-dominated stars with $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$21, $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$22 is significantly smaller than expected for an ordinary neutron star of the same compactness (Leung et al., 2011). This line of argument is strengthened in the slow-rotation I–Love–Q analysis, which shows that once observationally relevant BM-only moments of inertia are used, the deviations from single-fluid relations are larger still (Cronin et al., 2023).

X-ray pulse-profile modeling is particularly sensitive to dark halos because the metric exterior to the baryonic surface is not vacuum if DM extends beyond it. A focused study of halo/core formation and pulse profiles found that the peak flux deviation is controlled by the compactness-like ratio $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$23, with

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$24

and $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$25 (Miao et al., 2022). For a compact halo with $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$26, the peak flux deviation can reach $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$27, while nearly halo-free systems with $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$28 should distort the flux by less than $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$29 (Miao et al., 2022). Bayesian analysis using NICER posteriors for PSR J0030+0451 and PSR J0740+6620 then yielded a 90% upper limit $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$30 in the nearly halo-free case, while in the pure dark-core case the posterior peaks around $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$31; the self-interaction strength $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$32 remained unconstrained (Miao et al., 2022).

Several papers convert recent multimessenger observations into bounds on dark content. With twelve covariant density functional EOSs and self-interacting fermionic DM fixed at $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$33 MeV and $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$34, one study finds a strong linear correlation

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$35

between the maximum allowed DM fraction and the maximum mass of the corresponding pure neutron star, with Pearson coefficient $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$36 (Hu et al., 29 Dec 2025). Propagating the observational maximum-mass distribution then yields

$$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$37

at 68% confidence (Hu et al., 29 Dec 2025). A more targeted study of PSR J0740+6620 argues that the simultaneous mass–radius measurement reduces the uncertainty in the DM central density ratio $$n^2=-n_\alpha n^\alpha,\qquad p^2=-p_\alpha p^\alpha,\qquad x^2=-n_\alpha p^\alpha,$$38 by more than 50% compared with using mass or radius separately, and concludes that the total DM fraction should be smaller than 2% when constrained by the observed neutron-star maximum mass alone and smaller than 0.3% with the simultaneous mass–radius measurement (Zhang et al., 9 Sep 2025).

Observational interpretation remains model dependent. Some studies invoke DANSs to explain unusually compact objects such as EXO 1745-248, 4U 1608-52, 4U 1820-30 (Rezaei, 2016), HESS J1731-347 (Lopes et al., 2024), or XTE J1814-338 (Lopes et al., 2024). In each case the mechanism is the same: dark matter softens or compresses the visible star, allowing smaller apparent radii than standard hadronic EOSs typically admit. A plausible implication is that DANSs are most difficult to distinguish from ordinary neutron stars when the dark fraction is small and core-like, and easiest to distinguish when the star carries a sizable halo, an anomalously small visible radius, or a moment of inertia or tidal response incompatible with one-fluid expectations.

The main unresolved issue across this literature is not whether mathematically consistent DANS equilibria exist—they do in many formulations—but how much dark matter real neutron stars can accumulate and which microphysical model is realized in nature. The review literature therefore treats DANSs not as a single astrophysical object class with settled parameters, but as a broad two-fluid framework connecting dark-sector microphysics to neutron-star structure, stability, and multimessenger phenomenology (Grippa et al., 2024).