# Effects of a Brueckner-Hartree-Fock-corrected effective mass on speed of sound, conformality, and observables of dark matter-admixed neutron stars

Arijit Das, Prashanth Jaikumar, Adarsh Karekkat, Tanumoy Mandal

► To cite this version:

Arijit Das, Prashanth Jaikumar, Adarsh Karekkat, Tanumoy Mandal. Effects of a Brueckner-Hartree-Fock-corrected effective mass on speed of sound, conformality, and observables of dark matter-admixed neutron stars. Physical Review C, 2025, 112 (5), pp.055803. 10.1103/9qlh-w64b. hal-05379157

HAL Id: hal-05379157

https://hal.science/hal-05379157v1

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Effects of a Brueckner-Hartree-Fock-corrected effective mass on speed of sound, conformality, and observables of dark matter-admixed neutron stars

Arijit Das<sup>1,1,*</sup> Prashanth Jaikumar<sup>2,2,†</sup> Adarsh Karekkat<sup>1,3,‡</sup> and Tanumoy Mandal<sup>1,§</sup>

<sup>1</sup>Indian Institute of Science Education and Research Thiruvananthapuram, Vithura, Kerala 695 551, India

<sup>2</sup>California State University Long Beach, Long Beach, California 90840, USA

<sup>3</sup>Université de Caen Normandie, ENSICAEN, CNRS/IN2P3, LPC Caen UMR6534, F-14000 Caen, France

<figure><img src="image_2.png" /></figure>

(Received 23 June 2025; accepted 12 September 2025; published 20 November 2025)

We construct an equation of state describing cold and dense matter in the core of neutron stars which includes an admixture of fermionic dark matter and incorporates nucleon effective masses derived from the relativistic Brueckner-Hartree-Fock (BHF) many-body approach within a relativistic mean-field model. Such a BHF-informed mixed-model approach increases stellar compactness, with mass-radius configurations that are consistent with smaller, lighter pulsars. The model displays the expected nonmonotonic behavior of sound speed hinted at by neutron-star data and is closer to the conformal bound at maximum mass. We find that the model displays tension with bounds on heavier pulsars, suggesting that the hypothesis of an aggregated dark component in neutron stars needs further critical study.

DOI: 10.1103/9qlh-w64b

## I. INTRODUCTION

Neutron stars (NSs) serve as natural astrophysical laboratories for studying strongly interacting matter at supernuclear densities. Observational data—such as that gathered from gravitational wave (GW) events, pulsars, and other ultracompact objects (refer to Table I for a list of such experimental data)—help constrain the NS equation of state (EoS). Neither perturbative QCD (pQCD) [1] nor lattice QCD [2] can be reliably applied to describe strongly interacting matter at the intermediate densities found in the inner core of a NS. In this density regime, a useful approach is to extrapolate relativistic mean-field (RMF) models—fitted to empirical properties of finite nuclei and nuclear matter at saturation—with subsequent refinements, such as cross-couplings, density-dependent couplings, short-range correlations, or constraints from chiral effective-field theory (ChEFT) [3].

Studies within the frameworks of density-dependent relativistic mean-field theory [27], density-dependent relativistic Hartree-Fock theory [28], and density-dependent quark mean-field theory [29] have demonstrated that the density dependence of the nucleon’s effective mass [30] impacts both the

static and dynamic properties of NSs. Hence, it is important to account for density-dependent nucleon masses when making precise predictions of NS observables. Within the RMF framework, such density-dependent mass ratios are typically obtained by self-consistently solving the equations of motion (EoMs) of the underlying model. However, in the poorly understood intermediate baryon density range of $2\rho_0-5\rho_0$ (where $\rho_0$ is the nuclear saturation density), many-body forces become increasingly important and may lead to structural changes in NS matter [31]. These effects have been shown to influence core-collapse supernovae and binary NS mergers [32], the cooling of NSs [33,34], and neutrino emissivity from NSs via the Urca process [35–37].

Given the pervasive importance of many-body effects, it is instructive to explore different ways of incorporating them into the RMF framework, albeit with a partial relaxation of the requirement for a fully self-consistent solution of the mean-field EoS. One way to implement this idea is to replace the (Dirac) effective nucleon mass in the RMF approach with an effective-mass ratio obtained from realistic microscopic considerations. This represents an attempt to reconcile the RMF framework with more realistic in-medium nucleon interactions. For cold and dense matter, we adopt values of the nucleon effective mass calculated within the relativistic Brueckner-Hartree-Fock (BHF) framework [38], which includes microscopic two-body and three-body forces [39,40]. These results have also been used in studies of mass-radius relationships for purely hadronic NSs [41–43]. Moreover, as noted in Ref. [44], a key limitation of many RMF models is their inability to support $\approx 2M_\odot$ NSs, mainly due to the absence of three-body interactions in standard RMF frameworks. Introducing realistic, density-dependent meson-nucleon couplings that capture the effects of these three-body forces has been proposed

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP³.

<footnote>*Contact author: arijit21@iisertvm.ac.in</footnote>

<footnote>†Contact author: prashanth.jaikumar@csulb.edu</footnote>

<footnote>‡Contact author: adarsh.karekkat@unicaen.fr</footnote>

<footnote>§Contact author: tanumoy@iisertvm.ac.in</footnote>
<table><caption>TABLE I. Summary of masses ($M/M_\odot$), radii ($R$), tidal deformabilities ($\Lambda$), and inferred moments of inertia ($I$) of various pulsars, GW events, and other compact objects. The $M-R$ constraints for 4U 1724-207, 4U 1608-52, and EXO 1745-248 are taken from Ref. [4]. Similar constraints for HESS J1731-347 and PSR J1824-2452I (M28I) are taken from Refs. [5] and [6] respectively. In our comparison with theoretical predictions, we have projected the moment of inertia (MoI) of PSR J0030+0451 [7], estimated at a mass of $M = 1.4 M_\odot$, across the full range of mass and radius estimates [8,9].</caption><thead><tr><th>Stellar object</th><th><i>M</i>/<i>M</i><sub>☉</sub></th><th><i>R</i> (km)</th><th>&Lambda;</th><th><i>I</i> (10<sup>45</sup> gm cm<sup>2</sup>)</th></tr></thead><tbody><tr><td>PSR J0030+0451</td><td>1.34<sup>+0.15</sup><sub>-0.16</sub> [8,9]</td><td>12.71<sup>+1.14</sup><sub>-1.19</sub> [8,9]</td><td>370<sup>+360</sup><sub>-130</sub> [7]</td><td>1.43<sup>+0.30</sup><sub>-0.13</sub> [7]</td></tr><tr><td>PSR J0348+0432</td><td>1.97–2.05 [10]</td><td>12.16–12.948 [11]</td><td>–</td><td>1.594–1.9073 [12]</td></tr><tr><td>PSR J0740+6620</td><td>2.08 &pm; 0.07 [13]</td><td>12.39<sup>+1.30</sup><sub>-0.98</sub> [14]</td><td>–</td><td>4.65<sup>+1.16</sup><sub>-0.82</sub> [15]</td></tr><tr><td>PSR J1614-2230</td><td>1.97 &pm; 0.04 [16]</td><td>8.6–13.8 [17]</td><td>–</td><td>–</td></tr><tr><td>PSR J0437–4715</td><td>1.418 &pm; 0.037 [18]</td><td>11.36<sup>+0.95</sup><sub>-0.63</sub> [18]</td><td>–</td><td>–</td></tr><tr><td>PSR J1231-1411</td><td>1.12 &pm; 0.07 [19]</td><td>9.91<sup>+0.88</sup><sub>-0.86</sub> [19]</td><td>–</td><td>–</td></tr><tr><td>PSR J0737-3039 A</td><td>1.377 &pm; 0.005 [20]</td><td>8.14–25.74 [21]</td><td>–</td><td>1.15<sup>+0.38</sup><sub>-0.24</sub> [22]</td></tr><tr><td>PSR J0737-3039 B</td><td>1.250 &pm; 0.005 [20]</td><td>–</td><td>–</td><td>–</td></tr><tr><td>GW170817</td><td>1.17–1.60 [23]</td><td>11.9 &pm; 1.4 [24]</td><td>70–580 [23]</td><td>1.43<sup>+0.30</sup><sub>-0.13</sub> [7]</td></tr><tr><td>GW190425</td><td>1.12–2.52 [25]</td><td>–</td><td>1.4<sup>+3.8</sup><sub>-1.2</sub> &times; 10<sup>3</sup> [26]</td><td>–</td></tr></tbody></table>

as a solution to this issue. In a similar manner, a valid motivation for studying the impact of incorporating BHF-informed, density-dependent nucleon effective masses into RMF frameworks is to investigate their effect on the maximum mass of NSs.

Many-body interactions can leave their imprints on the speed of sound ($C_S$) profile, which, in turn, influences the stiffness of the EoS and thereby influences the maximum mass and radius of NS. The sound speed also provides insight into the microscopic composition of matter at supernuclear densities. Unfortunately, the limitations of theoretical approaches based on first principles result in large uncertainties in the sound speed profile at intermediate densities. At very low baryon densities, $C_S^2 \ll 1$ [45], while at asymptotically high densities, pQCD calculations predict a gradual approach to the conformal limit, $C_S^2 \to 1/3$. Moreover, studies have shown that, contrary to the holographic hypothesis of $C_S^2 \le 1/3$ [46] at intermediate densities, this value is not a strict upper bound, and the profile may become nonmonotonic [31,47–50]. The holographic bound is in tension with observations of the most massive NSs [51], as well as with predictions from quarkyonic matter [52], high-density QCD models [53], and frameworks based on gauge-gravity duality [54]. The constraint $C_S^2 \le 1/3$ significantly reduces the predicted NS mass below $2M_\odot$ for a class of nonrelativistic and relativistic models, rendering them inconsistent with astrophysical observations [55]. Furthermore, in Ref. [56] it has been shown that imposing the upper bound $C_S^2 = 1/3$ significantly decreases the maximum mass, Kerr parameter, and moment of inertia (MoI) of NSs, preventing these quantities from reaching values derived using realistic EoSs or constrained by observations. Relativistic kinetic theory also provides additional density-dependent constraints on the speed of sound [57]. Given these uncertainties and constraints, it is instructive to ask whether the inclusion of many-body interactions in RMF models can shed further light on the behavior of $C_S$.

An additional motivation for pursuing a more microscopic approach lies in the study of conformality in NS matter. Threshold values of parameters such as the polytropic index [58], conformal distance [59], trace anomaly [60], and the

curvature of the energy per particle [61], which characterize the onset of conformality, have already been determined. In contrast with model-agnostic frameworks—which emphasize using state-of-the-art theoretical and multimessenger constraints to provide confidence intervals for these parameters [61]—we adopt a more microscopic approach by employing a BHF-informed model to study the emergence of conformality in NSs.

We also go beyond conventional NSs by examining the accumulation of dark matter (DM) in their cores. Large amounts of accreted DM may produce observable effects on NS properties. In Ref. [62], a study of dark matter-admixed neutron stars (DMANS) was conducted within the BHF framework, starting from two-body and three-body nuclear interactions derived from ChEFT. However, since ChEFT is applicable only within a limited density range ($1\rho_0-2\rho_0$) [45], it is reasonable to adopt an RMF model, informed by microscopic BHF results, to study the resulting properties of DMANS.

This paper is organized as follows: In Sec. II, we describe the model under consideration, present parameter sets and the EoS for hadronic and DMANS, compare the effective-mass ratios obtained from the RMF and BHF-informed models and discuss the crust-matching technique used to eliminate problematic regions of negative pressure in the EoS. Section III presents numerical results for static and slowly rotating NSs and compares the predictions of the RMF and BHF-informed models for both hadronic NSs and DMANS. These results are also compared with experimental observations. In Sec. IV, we examine the speed of sound profile and the relative degrees of conformality in the RMF and BHF-informed scenarios. Finally, Sec. V provides a summary and conclusions.

## II. THEORETICAL FRAMEWORK

Among the many models proposed to describe NSs, we adopt the DM admixed SU(2) chiral sigma model [63–65] for simplicity of illustration. Following Ref. [66], we present the relevant details of the model below.
A. Model Lagrangian

The hadronic part of the model is given as follows [64]:

$$ \mathcal{L}_H = \bar{\psi} \mathcal{D} \psi + \frac{1}{2}(\partial_\mu \bar{\pi} \partial^\mu \bar{\pi} + \partial_\mu \sigma \partial^\mu \sigma) - \frac{1}{4} \Omega_{\mu\nu} \Omega^{\mu\nu} \\ + \frac{1}{2} g_\omega^2 x^2 \omega_\mu \omega^\mu - \frac{1}{4} \bar{R}_{\mu\nu} \bar{R}^{\mu\nu} + \frac{1}{2} m_\rho^2 \bar{\rho}_\mu \bar{\rho}^\mu - U(\bar{x}). \quad (1) $$

Here, $\psi$ is the SU(2) nucleon isospin doublet. The operator $\mathcal{D}$ is defined as

$$ \mathcal{D} = i\gamma_{\mu}\partial^{\mu} - g_{\sigma}(\sigma + i\gamma_{5}\bar{\tau}\bar{\pi}) - g_{\omega}\gamma_{\mu}\omega^{\mu} - \frac{1}{2}g_{\rho}\gamma_{\mu}\bar{\tau}\bar{\rho}^{\mu}, \quad (2) $$

with $\pi$, $\sigma$, $\omega$, and $\rho$ are the mesonic fields, and $\Omega_{\mu\nu} = \partial_\mu \omega_\nu - \partial_\nu \omega_\mu$ and $\bar{R}_{\mu\nu} = \partial_\mu \bar{\rho}_\nu - \partial_\nu \bar{\rho}_\mu + g_\rho(\bar{\rho}_\mu \times \bar{\rho}_\nu)$ are the field strength tensors for $\omega$ and $\rho$ mesons, respectively. The mass of $\rho$ is denoted $m_\rho$, and $g_\sigma$, $g_\omega$, and $g_\rho$ are the coupling constants. The scalar potential $U(\bar{x})$ is defined as

$$ U(\bar{x}) = \frac{\lambda}{4}\bar{x}^4 + \frac{B}{6}\bar{x}^6 + \frac{C}{8}\bar{x}^8, \quad (3) $$

where $\bar{x}^2 = x^2 - x_0^2$ with $x^2 = \sigma^2 + \bar{\pi}^2$, and $x_0$ being the vacuum expectation value (vev) of $\sigma$. Here, $\lambda$, $B$, and $C$ are some constants in the potential.

The fermionic dark sector Lagrangian is given by

$$ \mathcal{L}_{DM} = \bar{\chi}(i\gamma_{\mu}\partial^{\mu} - m_{\chi})\chi. \quad (4) $$

Here, $\chi$ is the DM spinor and $m_\chi$ is its mass.

The new scalar mediator $\phi$ and the new vector mediator $\xi$ connect the hadronic sector and the dark sector via the following interactions,

$$ \mathcal{L}_I = -g_\phi \bar{\psi} \psi \phi - g_\xi \bar{\psi} \gamma_\mu \psi \xi^\mu - y_\phi \bar{\chi} \chi \phi - y_\xi \bar{\chi} \gamma_\mu \chi \xi^\mu. \quad (5) $$

Here, $g_\phi$, $g_\xi$, $y_\phi$, and $y_\xi$ are the interaction couplings. The kinetic and mass terms for the new mediators $\phi$ and $\xi$ are represented in the following Lagrangian:

$$ \mathcal{L}_M = -\frac{1}{4}\Xi_{\mu\nu}\Xi^{\mu\nu} + \frac{1}{2}m_\xi^2\xi_\mu\xi^\mu + \frac{1}{2}\partial_\mu\phi\partial^\mu\phi - \frac{1}{2}m_\phi^2\phi^2, \quad (6) $$

where the masses for $\xi$ and $\phi$ are denoted by $m_\xi$ and $m_\phi$, respectively, and $\Xi_{\mu\nu} = \partial_\mu \xi_\nu - \partial_\nu \xi_\mu$ is the field strength tensor for $\xi_\mu$.

The leptonic sector is completely decoupled from other sectors of the model and given by

$$ \mathcal{L}_{\ell} = \sum_{\ell} \bar{l}(i\gamma_{\mu}\partial^{\mu} - m_{\ell})\ell. \quad (7) $$

Here, $\ell$ (= $e, \mu$) is the spinor field for the charge lepton of flavor $\ell$ of mass $m_\ell$.

Collecting all the terms, the total model is given by

$$ \mathcal{L} = \mathcal{L}_H + \mathcal{L}_{DM} + \mathcal{L}_I + \mathcal{L}_M + \mathcal{L}_{\ell}. \quad (8) $$

The masses of the nucleons ($m_N$ where $N = p, n$), $\omega$ meson ($m_\omega$), and $\sigma$ meson ($m_\sigma$) arise from chiral symmetry breaking. If $\phi_0$ denotes the vev of the scalar field $\phi$, then the masses are given by

$$ m_N = g_\sigma x_0 + g_\phi \phi_0, \quad m_\omega = g_\omega x_0, \quad m_\sigma = \sqrt{2\lambda} x_0. \quad (9) $$

Due to the presence of the scalar mediator $\phi$, the nucleon mass now receives an additional contribution from $\phi_0$.

B. Equation of state

In the RMF approximation, one replaces the meson fields by their RMF expectation values and, hence, all the corresponding kinetic terms drop out of the RMF-approximated model. The effective masses of nucleons and the DM are then given by

$$ m_N^* = g_\sigma \sigma + g_\phi \phi, \quad m_\chi^* = m_\chi + y_\phi \phi, \quad (10) $$

where in the above two expressions, $\phi$ and $\sigma$ are now RMF values of those fields. In all the following discussions, except that for fermions, all other fields will be assumed to have their RMF values. Before presenting the EoS, we provide the reader with some necessary relations. The scalar density for a generic fermion ($f = N, \chi, \ell$) is given by

$$ \begin{aligned} \rho_f^S = \langle \bar{f} f \rangle &= \frac{\gamma_f}{2\pi^2} \int_0^{k_f} dk \frac{m_f^* k^2}{\sqrt{k^2 + m_f^{*2}}} \\ &= \frac{\gamma_f m_f^{*3}}{4\pi^2} \left[ x_f \sqrt{1 + x_f^2} - \sinh^{-1} x_f \right]. \end{aligned} \quad (11) $$

Here, $\gamma_f = 2$ is the spin degeneracy factor and $x_f = k_f/m_f^*$ where $k_f$ is the Fermi momentum. The vector densities [67] are given by

$$ \rho_f^V = \langle f^\dagger f \rangle = \langle \bar{f} \gamma^0 f \rangle = \frac{\gamma_f}{2\pi^2} \int_0^{k_f} dk k^2 = \frac{\gamma_f}{6\pi^2} k_f^3. \quad (12) $$

The general expressions (for a generic fermion $f = n, p, \ell, \chi$) for the pressure ($\mathcal{P}$) and energy density ($\mathcal{E}$) can be obtained as follows:

$$ \begin{aligned} \mathcal{P}_f &= \frac{\gamma_f}{6\pi^2} \int_0^{k_f} dk \frac{k^4}{\sqrt{k^2 + m_f^{*2}}} \\ &= \frac{\gamma_f m_f^{*4}}{48\pi^2} \left[ (-3x_f + 2x_f^3) \sqrt{1 + x_f^2} + 3 \sinh^{-1} x_f \right], \end{aligned} \quad (13) $$

$$ \begin{aligned} \mathcal{E}_f &= \frac{\gamma_f}{2\pi^2} \int_0^{k_f} dk k^2 \sqrt{k^2 + m_f^{*2}} \\ &= \frac{\gamma_f m_f^{*4}}{16\pi^2} \left[ (x_f + 2x_f^3) \sqrt{1 + x_f^2} - \sinh^{-1} x_f \right]. \end{aligned} \quad (14) $$

The expressions for the total pressure and energy density of DMANS matter are given as follows:

$$ \mathcal{P} = \sum_{f=N, \chi, \ell} \mathcal{P}_f - \Sigma + \Theta_{-}, \quad \mathcal{E} = \sum_{f=N, \chi, \ell} \mathcal{E}_f + \Sigma + \Theta_{+}. \quad (15) $$

The expressions for the quantities $\Sigma$ and $\Theta_\pm$ can be found from Eqs. (A1) and (A5) in Appendix A. In the above expressions, the density-dependent RMF fields $\sigma(\rho^V)$ and $\phi(\rho^V)$ ($\rho^V = \rho_n^V + \rho_p^V$) obtained from self-consistent solutions of Eqs. (A6) and (A10) (refer to Appendix A for discussions on the complete set of EoMs) should be used in order to numerically evaluate the EoS.
<table><caption>TABLE II. Parameter set for hadronic sector [66,68].</caption><thead><tr><th colspan="11">Hadronic sector</th></tr><tr><th>Parameter</th><th>C<sub>σ</sub></th><th>C<sub>ω</sub></th><th>C<sub>ρ</sub></th><th>B</th><th>C</th><th>m<sub>σ</sub></th><th>f<sub>π</sub> = x<sub>0</sub></th><th>K</th><th>E<sub>BE</sub>/A</th><th>J</th><th>L<sub>0</sub></th><th>ρ<sub>0</sub></th><th>Y(ρ<sub>0</sub>)</th></tr><tr><th>Units</th><th>(fm<sup>2</sup>)</th><th>(fm<sup>2</sup>)</th><th>(fm<sup>2</sup>)</th><th>(fm<sup>2</sup>)</th><th>(fm<sup>4</sup>)</th><th>(MeV)</th><th>(MeV)</th><th>(MeV)</th><th>(MeV)</th><th>(MeV)</th><th>(MeV)</th><th>(fm<sup>-3</sup>)</th><th></th></tr></thead><tbody><tr><td>Value</td><td>7.325</td><td>1.642</td><td>5.324</td><td>-6.586</td><td>0.571</td><td>444.614</td><td>153.984</td><td>231</td><td>-16.3</td><td>32</td><td>88</td><td>0.153</td><td>0.87</td></tr></tbody></table>

### C. Parameter sets and simplified EoS

There are five free parameters in the hadronic part of the model i.e., $C_σ = g_σ^2/m_σ^2$, $C_ω = g_ω^2/m_ω^2$, $C_ρ = g_ρ^2/m_ρ^2$, $B$, and $C$. These parameters of the hadronic sector are determined by fitting to symmetric nuclear matter (SNM) properties, as shown in Table II [66,68]. The tabulated values of SNM properties, e.g., saturation density ($ρ_0 = 0.153$ fm$^{-3}$), nucleon effective mass at saturation [$Y(ρ_0) = 0.87$], nuclear incompressibility ($K \approx 231$ MeV), binding energy per nucleon ($E_{BE}/A = -16.3$ MeV) and symmetry energy coefficient ($J = 32$ MeV), calculated using the specified parameter set, agree well with already-existing constraints. The value of the slope parameter at saturation ($L_0 = 88$) also falls well within already established constraints [69]. We use the parameter set for the DM sector presented in Table III. These values fall very well in agreement with the self-interaction cross-section constraints from the Bullet Cluster and the present thermal relic abundance of DM [70]. The couplings of the new mediators with the nucleons, $g_φ$ and $g_ξ$, are much smaller than the corresponding couplings with the DM, $y_φ$ and $y_ξ$. Therefore, the DM-nucleon scattering cross section agrees well with the available constraints [71]. In all our further considerations, we only consider NS composed solely of neutrons and DM [i.e., in Eq. (15), $f$ will run over $n$ and $χ$ only]. We also assume the following benchmark values to present our results: $M_χ = 5$ and 15 GeV and $k_χ = 80$ MeV. We will call these DM sector parameter sets BPI and BPII, respectively. The EoS $P = P(E)$ can be obtained by numerically eliminating $k_n$ from the expressions of $P$ and $E$.

### D. RMF vs BHF: Effective masses

As pointed out earlier, in the purely RMF-based approach to the DMANS model, Eqs. (A6) and (A10) must be solved self-consistently over a realistic range of $k_n$, with a fixed $k_χ$, to obtain $φ = φ(ρ^V)$ and $σ = σ(ρ^V)$ as functions of baryon density. Consequently, the effective neutron mass $m_n^*$ and the DM effective mass $m_χ^*$ both become functions of $ρ_n^V$. In all numerical calculations that follow in the RMF scenario, the functions $φ(ρ^V)$, $m_n^*(ρ^V)$, and $m_χ^*(ρ^V)$ must be used. In contrast, in the BHF scenario, we adopt the following approximate analytical fitting function for the effective-mass ratio at zero temperature

from Ref. [73]:

$$Y_{\text{BHF}} = a_1 + b_1 \beta + c_1 \beta^2 + (a_2 + b_2 \beta + c_2 \beta^2) \rho^V + \frac{d_1}{\rho^V} + d_4 \ln \rho^V. \quad (16)$$

The values of the fitting parameters $a_i$, $b_i$, $c_i$, and $d_i$ are presented in Table IV. The asymmetry parameter is $β = 1$ in our case, since it is defined as $β = (ρ_n^V - ρ_p^V)/(ρ_n^V + ρ_p^V)$. This fitting function is valid within the range $ρ_n^V \in [0.1, 0.8]$ fm$^{-3}$. In the BHF scenario, we use $m_n^* = Y_{\text{BHF}}m_n$ as the modified nucleon effective mass in Eq. (A6), and solve for $φ = φ_{\text{BHF}}(ρ^V)$. Consequently, the BHF-corrected effective DM mass becomes $m_χ^*(ρ^V) = m_χ + y_φφ_{\text{BHF}}(ρ^V)$. In Fig. 1(a), we compare the mass ratio $Y$ computed from the RMF approach with $Y_{\text{BHF}}$ for purely hadronic matter. One can clearly see that the BHF correction deviates significantly from the RMF mass ratio within the validity range of the BHF framework. Hence, notable differences in the EoS and macroscopic observables of NSs are expected. In Fig. 1(b), we show the radial density profiles $ρ(r)/ρ_0$ of $2M_⊙$ NSs in both the RMF and BHF scenarios. These have been obtained by solving the Tolman-Oppenheimer-Volkoff (TOV) equations as shown in Eq. (B1) in Appendix B. The density is noticeably higher in the BHF case, resulting in a smaller radius for NSs of the same mass compared with the RMF case.

Finally, we emphasize the choice of the BHF effective mass used in the following numerical analyses. Both the pressure $P$ and the energy density $E$ are evaluated at various values of the neutron Fermi momentum $k_n$, i.e., on different Fermi surfaces. On each Fermi surface, the Landau effective mass provides a more accurate characterization of quasiparticle dynamics in the presence of interactions, encapsulating modifications arising from two- and three-body forces. Therefore, to better account for the influence of three-body interactions in our numerical analysis for the BHF case, we employ the Landau effective mass rather than the Dirac effective mass.

In evaluating the BHF-corrected EoS, the energy density $E(r=R)$ at the stellar surface—defined by $P(r=R) = 0$— does not vanish. This may create an issue for the standard calculation of the tidal Love number $k_2$, which is required to determine the tidal deformability $Λ$. In such cases, a

<table><caption>TABLE III. Parameter sets for DM sector [66,72].</caption><thead><tr><th colspan="6">Dark matter sector</th></tr><tr><th>Parameters</th><th>m<sub>χ</sub> (MeV)</th><th>m<sub>φ</sub> (MeV)</th><th>m<sub>ξ</sub> (MeV)</th><th>y<sub>φ</sub> = y<sub>ξ</sub></th><th>g<sub>φ</sub> = g<sub>ξ</sub></th></tr></thead><tbody><tr><td>BP I</td><td>5</td><td>9</td><td>11</td><td>0.13</td><td>1.1 × 10<sup>-4</sup></td></tr><tr><td>BP II</td><td>15</td><td>20</td><td>34</td><td>0.21</td><td>1.1 × 10<sup>-4</sup></td></tr></tbody></table>
<table><caption>TABLE IV. BHF fitting parameters [73].</caption><thead><tr><th colspan="10">BHF fitting parameters</th></tr><tr><th>Parameter</th><th>a<sub>1</sub></th><th>b<sub>1</sub></th><th>c<sub>1</sub></th><th>a<sub>2</sub></th><th>b<sub>2</sub></th><th>c<sub>2</sub></th><th>d<sub>1</sub></th><th>d<sub>4</sub></th></tr></thead><tbody><tr><td>Value</td><td>0.102</td><td>-0.094</td><td>0.068</td><td>0.699</td><td>0.0354</td><td>0.0133</td><td>-0.0135</td><td>-0.367</td></tr></tbody></table>

crust-matching technique can be employed at low energy densities and pressures to reproduce crustal behavior using a reliable crust EoS [74]. Here, we adopt the well-known BSK19 crust EoS, matching it to the low-energy-density portion of our EoS. We locate the energy density at which the crust EoS intersects the hadronic (or DMANS) EoS—namely, the point where their pressures are equal—and then construct a piecewise-continuous EoS: to the left of the intersection, we replace the original EoS with the BSK19 crust branch, while to the right we retain the hadronic or DMANS branch, as appropriate.

### III. OBSERVABLES AND COMPARISON WITH DATA

In Fig. 2, we show a comparison of EoSs for the BPII parameter set obtained within the RMF and BHF approaches. The EoS lines are passing through a colored region which represents the constraint, $\Lambda_{1.4} < 400$ (where $\Lambda_{1.4}$ is the value of dimensionless tidal deformability $\Lambda$ at the benchmark mass of $M = 1.4M_{\odot}$), as obtained in Ref. [75]. In Ref. [75], an ensemble of EoS are generated with (i) ChEFT EoS for $\rho < 1.1\rho_0$, (ii) piecewise polytropes between the regions of validity of ChEFT and pQCD and (iii) NNLO pQCD results in the high density regime to describe deconfined quark matter. To obtain this constraint, the condition of existence of a $2M_{\odot}$ NS and the constraint on $\Lambda$ from GW170817 have been imposed. The trends of the slopes of EoSs in both scenarios have been found to be consistent with various RMF EoSs used in the literature [76]. The regions excluded by the

condition $P < 3(\rho^V m_n^* + \rho_\chi^V m_\chi^* + \mathcal{E}) = 3\mathcal{E}_{tot}$ from the relativistic kinetic theory considerations [57] always fall much above the EoS lines. The “matching” point on the EoS indicates the point where the crust EoS has been matched with the BHF-corrected EoS.

In Figs. 3(a) and 3(b), we compare the mass-radius ($M-R$) relations for purely hadronic and DMANS scenarios, in both RMF and BHF cases, respectively, with observational data from various pulsars listed in Table I. The $M-R$ relations are obtained by solving the TOV equations as shown in Eq. (B1) in Appendix B. We find that, at a fixed NS mass, the BHF corrections significantly reduce the radius compared with the RMF results, thereby increasing the compactness of the star. The RMF model yields maximum masses and radii that are broadly consistent with relatively heavier pulsars and GW observations. In contrast, the BHF-corrected results show tighter mass-radius curves, and while they provide better agreement with low-mass compact stars such as PSR J0437-4715, PSR J1231-1411, and the ultracompact object HESS J1731-347, they are not in good agreement with the masses of heavier pulsars, e.g., PSR J0740+6620.

Similarly, in Fig. 4, we compare the same $M-R$ relations with various GW events. The BHF cases are not in good agreement with the mass estimates of GW190425 M1(M2) [25], whereas the RMF-based results offer a better match. However, for GW170817 [24] event, the BHF predictions remain consistent, albeit with smaller radii compared with the RMF case. We show two special regions, as pointed out in Ref. [75], of tidal deformability $\Lambda_{1.4}$ at $M_{NS} = 1.4M_{\odot}$, i.e.,

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Comparison of RMF and BHF corrected nucleon effective-mass ratios for purely hadronic EoS. (b) The density profiles $\rho(r)/\rho_0$ for RMF and BHF hadronic $2.0M_{\odot}$ NS. The higher density in the BHF case leads to a more compact NS.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Comparison of EoSs for RMF and BHF-corrected cases for BP II parameter set. The constraint $\Lambda_{1.4} < 400$ is taken from Ref. [75].</figcaption></figure>

$\Lambda_{1.4} < 400$ and $400 < \Lambda_{1.4} < 800$. All the BHF curves lie within the region $\Lambda_{1.4} < 400$, which is consistent with the results shown later in Fig. 5(b).

Interestingly, while Ref. [77] required a large DM fraction (around 4.75%) to attain consistency with HESS J1731-347, our BHF-corrected DMANS model achieves good agreement with a significantly smaller DM fraction of $\approx 0.2\%$ (or even no DM at all). The purely RMF-based hadronic model fails to fall within the radius estimate of PSR J1824-2452I (M28I),

unlike its BHF counterpart. As the DM mass increases (for fixed $k_\chi$), both the maximum mass and corresponding radius of the NS decrease, consistent with previous studies [78–80]. However, we observed that for a fixed NS mass, the radius is less sensitive to $M_\chi$ in the BHF model compared with the RMF case. We note that hadronic configurations in the BHF scenario show better agreement with the transient low-mass X-ray binaries such as EXO 1745-248 and 4U 1608-52 [4], a feature not captured by the RMF-based hadronic model.

Although within our BHF approach, more compact objects can be better explained compared with the RMF approach, the ultracompact object XTE J1814-338 [81] cannot be explained in our framework. In the context of RMF models, studies using bosonic dark matter models [82], Higgs-mediated fermionic dark matter models [83], mirror dark matter ad-mixed strange star models [84] and hybrid star models [85] have all aimed to explain the ultracompact nature of XTE J1814-338. Given the lesser degree of sensitivity of stellar radius in the BHF case upon inclusion of DM, achieving consistency with XTE J1814-338 perhaps requires significantly larger DM masses compared with the RMF scenario.

In Fig. 5(a), we present the mass–compactness (M–C) relations for NSs within both the RMF and BHF frameworks, considering purely hadronic and DMANS. We find that the compactness is significantly less sensitive to variations in the DM mass in the BHF case compared with the RMF case. Using the derived compactness of PSR J0030+0451 [8], we find that in the BHF scenario, the inferred mass range for this pulsar lies between $(1.03 – 1.19)M_\odot$, which is lower than the observationally estimated mass range for PSR J0030+0451. In contrast, the RMF-based result yields a more consistent mass range of $(1.19 – 1.49)M_\odot$. Conversely, if one uses the known mass estimate of PSR J0030+0451 to infer

<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Mass-radius relationships for hadronic, BPI, and BPII parameters in the RMF scenario. Theoretical results are compared with relevant experimental results presented for pulsars in Table I. (b) Mass-radius relationships for hadronic, BPI, and BPII parameters in the BHF scenario. The M-R curves in the BHF scenario are completely within the $\Lambda(1.4M_\odot) < 400$ constraint obtained from Ref. [75], which is consistent with the result in Fig. 5(b).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 4. Mass-radius relationships for hadronic, BPI, and BPII parameters in both RMF and BHF scenarios are compared with various GW wave data.</figcaption></figure>

compactness, we obtain $C \in [0.13, 0.19]$ for the RMF case and $C \in [0.16, 0.20]$ for the BHF case. These results indicate that the BHF model is not fully compatible with both the mass and compactness constraints of PSR J0030+0451, while the RMF case provides a better fit. As for a particular BHF $M-R$ curve, the radius of the star is constrained in a very narrow range around $R = 10.9$ km, the compactness increases almost linearly with mass.

In Fig. 5(b), we present the variation of the dimensionless tidal deformability $\Lambda$ as a function of NS mass for both RMF and BHF models, considering pure hadronic as well as DMANS scenarios. While the BHF results do not agree well with the benchmark value $\Lambda_{1.4} \approx 370^{+360}_{-130}$ inferred from the combined analysis of GW170817 and NICER observations of PSR J0030+0451 [7], they remain well within the broader range $\Lambda_{1.4} \approx 70-580$ derived from the GW170817 data alone [23,24]. We also compare theoretical results with 68% confidence limits on GW190425 M1(M2) [86]. At $M_{NS} = 1.4M_{\odot}$, we observe that, for a fixed $k_{\chi}$, the tidal deformability $\Lambda$ in the BHF scenario exhibits reduced sensitivity to the DM mass compared with the RMF case. This reduced sensitivity is a consequence of the increased stiffness of the NS matter resulting from BHF corrections, as discussed earlier. For completeness, Table V provides a summary of the

key macroscopic observables—including the maximum mass $M_{max}$, corresponding radius $R_{peak}$ and $\Lambda_{1.4}$—for the benchmark DM parameter sets BPI and BPII, in both RMF and BHF frameworks.

The quantities $C$ and $\Lambda$ are connected via Eq. (B12), where another quantity called tidal Love number $k_2$ is also present (see Appendix B for its definition and related discussions). The variation of $k_2$ with DM mass in the RMF and BHF frameworks has been found to be qualitatively opposite. At a fixed stellar mass, $k_2$ mildly increases with DM mass in the RMF case, whereas it decreases mildly in the BHF case. Understanding the physical origin of this small variation of $k_2$ with respect to (w.r.t.) the DM mass is difficult, given its complex dependence on $C$ and the metric perturbation parameter $y_R$. However, such mild variations of $k_2$ also have marginal effects on $\Lambda$ compared with the effect of compactness $C$. The variation in $\Lambda$ is heavily dependent on the variation in $C$ as $\Lambda \sim C^{-5}$. In the RMF scenario, $C$ increases upon inclusion of DM (e.g., at $M = 1.4M_{\odot}$, $C$ is larger by 13% for the BPII parameter set compared with the hadronic scenario). Given the strong dependence of $\Lambda$ on $C$, this leads to $\approx 45\%$ decrease in $\Lambda$ for the BPII set compared with the hadronic case. A similar trend is also observed within the BHF results. Additionally, for a fixed $M$, considerable increase of $C$ in the BHF case (by $\approx 17\%$ in the BHF hadronic case compared with RMF) leads to a large decrease in $\Lambda$ ($\approx 82\%$) in the BHF case.

We present mass-MoI ($M-I$) relation in Figs. 5(c) in the slow-rotation approximation, where the rotation frequency of the star is always much smaller than the Kepler frequency and the MoI becomes largely independent of the rotation frequency of the star. Relevant expressions for the calculations of MoI are given in Appendix B. Although for a fixed $k_{\chi}$, the MoI decreases for increasing $M_{\chi}$ in both scenarios, variation of MoI is less sensitive to the DM mass in the BHF case. Similar behavior has been observed w.r.t. radius. Both the RMF and BHF cases are consistent with mass-MoI estimates of PSR J0030+0451 [7]. The mass-MoI estimate of light pulsar PSR J0737-3039A [22,87] is more consistent with the BHF scenario. Figure. 5(d) presents radius-MoI (R-I) relations. The BHF case has found not to be in good agreement with radius-MoI estimate of PSR J00340+0451. The inferred MoI of PSR J0737-3039A constrains the radius in the BHF star radii smaller ($\approx 10.6-11$ kms) than the RMF star ($\approx 11.4-13.3$ kms). The BHF results are also consistent with the radius measurement of light pulsar PSR J1231-1411. We observe that for stars with identical mass, $I_{BHF}$ turns out to be smaller than $I_{RMF}$. For instance at $M = 1.4M_{\odot}$ (for pure hadronic matter), the ratios,

<table><caption>TABLE V. Maximum masses $M_{max}$ of NSs, corresponding radii $R_{peak}$, and tidal deformability at $M_{NS} = 1.4M_{\odot}$ (denoted by $\Lambda_{1.4}$) for purely hadronic and DMANS (BPI and BPII) in RMF and BHF frameworks.</caption><thead><tr><th rowspan="2">Observables</th><th colspan="3">RMF</th><th colspan="3">BHF</th></tr><tr><th>Hadronic</th><th>BPI</th><th>BPII</th><th>Hadronic</th><th>BPI</th><th>BPII</th></tr></thead><tbody><tr><td><i>M</i><sub>max</sub> (<i>M</i><sub>☉</sub>)</td><td>2.13</td><td>2.08</td><td>2.00</td><td>2.03</td><td>1.99</td><td>1.93</td></tr><tr><td><i>R</i><sub>peak</sub> (km)</td><td>11.81</td><td>11.34</td><td>10.69</td><td>10.41</td><td>10.29</td><td>10.04</td></tr><tr><td>Λ<sub>1.4</sub></td><td>532</td><td>493</td><td>336</td><td>134</td><td>125</td><td>109</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIG. 5. Compactness (C), tidal deformability (Λ) and MoI (I) of NSs as functions of mass (M in M⊙ units) are presented in panels (a), (b), and (c), respectively. Panel (d) shows the variation of I with stellar radius (R). Results obtained for purely hadronic and DMANS (for BPI and BPII) cases in RMF and BHF scenarios are compared with the experimental findings presented in Table I.</figcaption></figure>

$R_{BHF}/R_{RMF} \approx 0.8$, $\Omega_{BHF}/\Omega_{RMF} \approx 1.1$ and $\tilde{\omega}'_{BHF}/\tilde{\omega}'_{RMF} \approx 1.8$ (here $\tilde{\omega}' = d\tilde{\omega}/dr$). These ratios give $I_{BHF}/I_{RMF} \approx 0.68$ which is consistent with the values obtained from Fig. 5(c).

IV. SPEED OF SOUND AND CONFORMALITY

The speed of sound ($C_s$) inside an NS is a direct measure of the stiffness of the EoS. This is defined as the slope of the EoS curve via

$$C_s^2 = \frac{\partial P}{\partial \mathcal{E}}. \qquad (17)$$

In the previous section, we observed that incorporating the BHF-corrected effective nucleon mass leads to significant differences in the M-R relations, tidal deformability, compactness, and MoI. Although $C_s$ cannot be directly measured, it can be inferred indirectly from EoS constraints and other NS observables obtained through pulsar and GW data. Bayesian

methods are commonly used to infer $C_s$ from observational priors by sampling various EoS parametrizations. These studies suggest that to support the existence of $\approx 2M_\odot$ NSs, $C_s$ must exceed the conformal limit $C_s \le c/\sqrt{3}$ (where c is the speed of light in vacuum), typically peaking around $\rho \approx 3\rho_0-4\rho_0$, attaining values $C_s^2 \approx 0.6$ [88]. These results strongly indicate possible stiffening of the EoS in the NS interior. We present a comparative analysis of the $C_s$ and associated measures of conformality in the RMF and BHF frameworks. We note that the inclusion of a small fraction of DM in the NS does not significantly alter $C_s$ [83].

In Fig. 6(a), we show the variation of $C_s^2$ as functions of $\rho^V$, along with the slope ($\alpha$) and curvature ($\beta$) of the energy per particle (see Appendix C for a brief discussion), in both the RMF and BHF frameworks. These theoretical predictions are compared with the confidence intervals derived in Ref. [61], in which a piecewise-linear speed-of-sound parametrization was adopted for the EoS in the density range [0.5, 1.1]$\rho_0$ [58],
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 6. (a) Comparison of $C_S^2$, $\alpha$, and $\beta$ [89] in RMF and BHF cases as functions of $\rho^V$. The green vertical shaded region bounded by black dotted lines encloses the $1\sigma$ credible region for the peak of $C_S^2$ [90]. The shaded regions bounded by purple and green vertical lines enclose $1\sigma$ credible intervals for central densities of $1.4 M_\odot$ and $2.0 M_\odot$ NSs [59]. (b) EoSs on the $C_S^2-(C_S^2)$ plane for hadronic and DMANS matter, with RMF and BHF curves. Colored dots mark $\rho^V \approx 1.05 \text{ fm}^{-3}$ ($\approx 7\rho_0$). The solid black line shows the mean EoS from a statistical ensemble consistent with constraints in Ref. [89]. (c) Variation of $\beta$ as a function of $\rho^V$ for hadronic and DMANS cases. (d) Variation of conformal distance $d_c$ with $\rho^V$.</figcaption></figure>

while using the Baym-Pethick-Sutherland EoS at subnuclear densities and pQCD EoS at asymptotically high densities. The observational inputs used in their analysis include the tidal deformability $\Lambda_{1.4}$ from GW170817 [24] and the mass of PSR J0952-0607 [91], leading to constraints expressed as functions of energy density $\mathcal{E}$. For comparison with our results expressed in terms of $\rho^V$, we have rescaled the energy density values in their constraints by a factor of $\rho_0/\mathcal{E}_0$, where $\mathcal{E}_0 \approx 150 \text{ MeV/fm}^3$ denotes the nuclear saturation energy density.

The nonmonotonic behavior of $C_S^2$ in the BHF scenario is consistent with earlier studies [31,47–49]. At sufficiently high baryon densities, the breach of the conformal limit is also supported by existing literature [48,88]. However, we

note a tension between the radii of maximally massive stars ($M \approx 2M_\odot$) in the BHF case ($R = 10.4$ km in the purely hadronic case) and the corresponding values reported in Ref. [48]. In Ref. [58], it was shown that for neutron stars with mass $M \approx 2M_\odot$, if the conformal bound is not strongly violated, the stars are predicted to contain quark-matter cores. In our BHF-based hadronic EoS, the maximum mass configuration occurs at central energy density $\mathcal{E}_c \approx 1412 \text{ MeV/fm}^3$, corresponding to $\rho^V \approx 1.05 \text{ fm}^{-3}$. This lies slightly outside the validity range of the BHF fitting function. Nonetheless, we extrapolate the BHF parametrization up to this density and find that at high $\rho^V$, the BHF corrections bring $C_S^2$ close to the conformal limit. This behavior could be indicative of a
transition toward weakly coupled quark matter [92,93]. Such a trend can be attributed to the large negative values of $\beta_{BHF}$ at high $\rho^V$, as illustrated in Fig. 6(a). Furthermore, the values of $C_S^2$ in the BHF framework lie within the $1\sigma$ confidence band presented in Ref. [61] over the range $\rho^V \approx [0.6-1.05] \text{ fm}^{-3}$.

We also compare our results with the $1\sigma$ credible region for the location of the peak of $C_S^2$ obtained in Ref. [90], as well as with the central densities corresponding to $1.4M_\odot$ and $2.0M_\odot$ NSs, as reported in Ref. [59]. To ensure that the $C_S$ in the BHF scenario remains causal, apart from $C_S^2 < 1$, we also verify that it satisfies the density-dependent upper bound derived from relativistic kinetic theory [57], (shown by RKT + RMF/BHF in the plot where “RKT” is used as an abbreviation for “relativistic kinetic theory”). The constraint is given by

$$C_S^2 \le \frac{\mathcal{E} + \rho^V m_n^* + \rho^V \chi m^* \chi - \mathcal{P}/3}{\mathcal{E} + \rho^V m_n^* + \rho^V \chi m^* \chi + \mathcal{P}}. \quad (18)$$

We have used the appropriate density-dependent effective masses for all fermionic species to evaluate this inequality.

From Eqs. (17) and (C2), the EoSs can be obtained on the $C_S^2 - \langle C_S^2 \rangle$ plane as shown in Fig. 6(b) in both RMF and BHF scenarios for both hadronic and DMANS matter. The contour $d_c = 0.2$ [89] (refer to Appendix C for a brief discussion on various measures of conformality) separates nonconformal ($d_c > 0.2$) and nearly conformal ($d_c \le 0.2$) behavior. Near $\rho^V \approx 7\rho_0 \approx 1.05 \text{ fm}^{-3}$, the BHF case is much closer to full conformality (the point of intersection of the lines $C_S^2 = 1/3$ and $\langle C_S^2 \rangle = 1/3$) than the RMF case. Similarly, one can get the EoS on the $\Delta - C_S^2$ plane to arrive at the same conclusion. Moreover, unlike the RMF scenario, the BHF EoS in the latter case has been found to be consistent with the $1\sigma$ interval of $\Delta$ as a function of $C_S^2$ at the center of maximally stable heavy NSs [90].

In Fig. 6(c), we show the variation of $\beta$ as functions of $\rho^V$. Conformality sets in for both RMF and BHF around the same density of $\rho^V \approx 0.8 \text{ fm}^{-3}$. At even higher densities, both are consistent with the $2\sigma$ confidence interval. Moreover, unlike RMF, within the intermediate density range of [0.3,0.6] $\text{fm}^{-3}$, the BHF scenario is consistent with the $1\sigma$ confidence interval. An enhanced values of $\beta$ in the BHF case in the range $\rho^V \in [0.3-0.8] \text{ fm}^{-3}$ implies an enhanced value of compression modulus [94].

In Fig. 6(d), we show the variation of the conformal distance $d_c$ [59] as functions of $\rho^V$. In the density range of [0.3-0.95] $\text{fm}^{-3}$, the BHF scenario is more consistent with the $1\sigma$ limits than RMF. However, for densities beyond 0.95 $\text{fm}^{-3}$, this trend is reversed.

Besides these results, we have also studied the variation of the polytropic index ($\gamma$) with $\rho^V$. In the range of 0.4-0.7 $\text{fm}^{-3}$, the BHF case has been found to be much more consistent with the $1\sigma$ limit [61] than the RMF case. At even higher densities, both cases are consistent with $2\sigma$ confidence limit. For the BHF case, the onset of conformality occurs ($\gamma \approx 1.75$, a necessary but not a sufficient condition [95]) at a baryon density slightly higher than the RMF case. We have also studied the variation of the conformal distance $d_c$ with $C_S^2$. The RMF EoS crosses the $d_c = 0.2$ conformality

threshold at lower values of $C_S^2$ compared with the BHF case. In the high-density limit, BHF cases are much closer to full conformality than the RMF cases, but at a lower $C_S^2$ compared with the RMF.

## V. CONCLUSION AND DISCUSSION

In this work, we have studied the properties of ultradense matter, speed of sound, and various aspects of conformality in NSs with and without an admixture of DM, using a relativistic mean-field model improved by the inclusion of realistic medium effects of many-body interactions. Specifically, we implemented BHF-informed effective nucleon masses within the RMF framework [73] to model interacting fermionic DM admixed with normal nuclear matter described by the SU(2) chiral sigma model [66]. We have shown the effects of such many-body interactions on macroscopic observables, speed of sound and conformality of purely hadronic and DMANS.

The BHF-corrected EoS for static neutron stars lowers the maximum mass ($1.9M_\odot-2.0M_\odot$) compared with the purely RMF case ($2M_\odot-2.2M_\odot$). The inclusion of BHF corrections increases the compactness of the star, allowing for improved consistency with smaller compact objects, e.g., HESS J1731-347 and PSR J1231-1411. In contrast with Ref. [77], our results are in agreement with HESS J1731-347 with a much smaller density of DM. In the BHF scenario, as the DM mass increases, the maximum masses and the corresponding radii decrease; also, the radii are less sensitive to DM mass compared with the RMF case. Large DM masses (e.g., 15 GeV or more) lead to tension with the mass of very heavy pulsars, but are still consistent with relatively low mass compact objects, e.g., HESS J1731-347. Dimensionless tidal deformabilities ($\Lambda$) are consistent with measurements on GW170817 [23,24] but in tension with joint analysis of PSR J0030+0451 and GW170817 [7]. These effects demonstrate that current DM admixed RMF models, which do not include many-body effects between nucleons, need to be improved if they are to achieve consistency with data on heavy NSs. The M-R relations obtained within the BHF framework are well within the constraint $\Lambda(1.4M_\odot) < 400$ [75], a result confirmed by the corresponding M-$\Lambda$ relations.

The BHF correction also changes the behavior of the speed of sound by making its profile nonmonotonic [31,47-49] and improving agreement in the high density regime with the available confidence limits from Ref. [61]. Such behavior of $C_S^2$ could be indicative of the onset of a weakly coupled quark-matter phase [92,93]. For both RMF and BHF scenarios, the $C_S^2$ profile is consistent with constraints from relativistic kinetic theory [57]. The BHF-corrected EoS approaches conformality more closely than the RMF case near $\rho_n^V \approx 7\rho_0 \approx 1.05 \text{ fm}^{-3}$, as illustrated on the $C_S^2-\langle C_S^2 \rangle$ plane. BHF results better match the $1\sigma$ confidence limits for conformality parameters, e.g., $\beta$ [59] and $d_c$ [61].

A brief discussion of some possible extensions of the present work is in order. The construction of the SU(2) chiral sigma model used in this work employs nucleons that are degenerate in mass. A more realistic model with isospin-breaking effect by introducing an isovector scalar field could be considered. Different effective-mass profiles for n and p
can be used [73] in this scenario to obtain more refined results. Effects of couplings between mesons, rapid rotation, and finite temperatures would further influence the present results, as would the introduction of conjectured additional degrees of freedom such as hyperons along with nucleon-hyperon three-body interactions. Effects of external parameters such as temperature and strong magnetic field on β-equilibrated, charge-neutral nuclear and quark matter [96–99] in the context of NS observables might be an interesting avenue to investigate in the future. Moreover, because of differences in the speed of sound profile of the BHF scenario compared with the RMF, it would be interesting to investigate the effects of such differences in the nonradial modes of oscillations of NS, such as r modes [100] and g modes [101].

## ACKNOWLEDGMENTS

T.M. acknowledges partial support from the SERB/ANRF, Government of India, through the Core Research Grant (CRG) No. CRG/2023/007031. P.J. is supported by a grant from the U.S. National Science Foundation Grant No. PHY-2310003. A.D. thanks Dr. Sanchari Goswami for carefully reading the manuscripts and sharing her valuable comments on the paper.

## DATA AVAILABILITY

No data were created or analyzed in this study.

## APPENDIX A: EOS AND EOM

The quantity $\Sigma$ introduced in Eq. (15) is defined as follows:

$$ \Sigma = \frac{g_{\sigma}^2 x_0^2}{C_{\sigma}} \left[ \frac{\bar{Y}^4}{8} - \frac{B \bar{Y}^6}{12 C_{\sigma}} + \frac{C \bar{Y}^8}{16 C_{\sigma}^2} \right], \quad (A1) $$

where $\bar{Y}^2 = 1 - Y^2$. This captures the contribution to the EoS coming from the scalar potential $U(\bar{x})$. Using the relations, $C_{\sigma} = g_{\sigma}^2/m_{\sigma}^2$ and $C_{\omega} = g_{\omega}^2/m_{\omega}^2 = x_0^{-2} = f_{\pi}^{-2}$, the above expression for $\Sigma$ simplifies to

$$ \Sigma = m_{\sigma}^2 f_{\pi}^2 \left[ \frac{\bar{Y}^4}{8} - \frac{f_{\pi}^2 B \bar{Y}^6}{12} + \frac{f_{\pi}^4 C \bar{Y}^8}{16} \right]. \quad (A2) $$

In the above, we have used the nucleon effective-mass ratio $Y$ defined as

$$ Y = \frac{g_{\sigma} \sigma + g_{\phi} \phi}{m_N}. \quad (A3) $$

The effective nucleon mass is given by $m_N^* = Ym_N$. Similarly, the effective DM mass ratio is given by

$$ Y_{DM} = \frac{m_{\chi} + y_{\phi} \phi}{m_{\chi}}. \quad (A4) $$

The quantity $\Theta_{\pm}$, introduced in Eq. (15), is defined as follows:

$$ \Theta_{\pm} = \frac{C_{\omega}(\rho^{\mathrm{V}})^2}{2Y^2} + \frac{1}{2}m_{\rho}^2(\rho^{03})^2 + \frac{1}{2}m_{\xi}^2(\xi^0)^2 \pm \frac{1}{2}m_{\phi}^2\phi^2, \quad (A5) $$

where the first (second) index on $\rho$ denotes the Lorentz [SU(2)] index. All other components of $\bar{\rho}_{\mu}$ are zero. This captures contributions coming from the mass terms of various mesons. In this expression $\rho^{\mathrm{V}} = \rho_{n}^{\mathrm{V}} + \rho_{p}^{\mathrm{V}}$. However, in the case of a hadronic sector composed solely of neutrons, $\rho^{\mathrm{V}} = \rho_{n}^{\mathrm{V}}$ (since $\rho_{p}^{\mathrm{V}} = 0$). The EoM for $\phi$ is

$$ \phi = - \frac{g_{\phi}}{m_{\phi}^2} \rho^{\mathrm{S}} - \frac{y_{\phi}}{m_{\phi}^2} \rho_{\chi}^{\mathrm{S}}, \quad (A6) $$

where $\rho^{\mathrm{S}} = \rho_{n}^{\mathrm{S}} + \rho_{p}^{\mathrm{S}}$. Using the familiar results of $\langle \bar{\psi} \gamma^i \psi \rangle = 0 = \langle \bar{\chi} \gamma^i \chi \rangle$ ($i = 1, 2, 3$) at zero temperature [67], and hence, $\xi^i = 0$, we get the EoM for $\xi^0$ as

$$ \begin{aligned} \xi^0 &= \frac{g_{\xi}}{m_{\xi}^2} \rho^{\mathrm{V}} + \frac{y_{\xi}}{m_{\xi}^2} \rho_{\chi}^{\mathrm{V}} \\ &= \frac{g_{\xi}}{m_{\xi}^2} \left( \frac{\gamma_n}{6\pi^2} k_n^3 + \frac{\gamma_p}{6\pi^2} k_p^3 \right) + \frac{y_{\xi}}{m_{\xi}^2} \frac{\gamma_{\chi}}{6\pi^2} k_{\chi}^3. \end{aligned} \quad (A7) $$

Here, $\rho^{\mathrm{V}} = \rho_{n}^{\mathrm{V}} + \rho_{p}^{\mathrm{V}}$. Similar reasoning leads to $\omega^i = 0$ and we get the EoM of $\omega^0$ as

$$ \omega^0 = \frac{\rho^{\mathrm{V}}}{g_{\omega} \sigma^2} = \frac{1}{g_{\omega} \sigma^2} \left( \frac{\gamma_p}{6\pi^2} k_p^3 + \frac{\gamma_n}{6\pi^2} k_n^3 \right). \quad (A8) $$

Similar assumptions lead to

$$ \rho^{03} = \frac{g_{\rho}}{2m_{\rho}^2} (\rho_{n}^{\mathrm{V}} - \rho_{p}^{\mathrm{V}}) = \frac{g_{\rho}}{2m_{\rho}^2} \left( \frac{\gamma_n}{6\pi^2} k_n^3 - \frac{\gamma_p}{6\pi^2} k_p^3 \right). \quad (A9) $$

The EoM for $\sigma$ is as follows:

$$ \frac{2C_{\sigma} \rho^{\mathrm{S}}}{Y(m_N - g_{\phi} \phi_0)} - \frac{2C_{\sigma} C_{\omega} \rho^{\mathrm{V}}^2}{Y^4(m_N - g_{\phi} \phi_0)^2} - \bar{Y}^2 + \frac{B}{C_{\omega}} \bar{Y}^4 \\ - \frac{C}{C_{\omega}^2} \bar{Y}^6 = 0. \quad (A10) $$

Using the expression for $\rho_{n}^{\mathrm{V}}$ from Eq. (12) and using the EoMs for $\rho^{03}$ and $\xi^0$ (shown above) in the pure neutron matter case, the expression of $\Theta_{\pm}$ simplifies to

$$ \begin{aligned} \Theta_{\pm} = \frac{\gamma_N^2 C_{\omega} k_n^6}{72\pi^4 Y^2} + \frac{\gamma_N^2 C_{\rho} k_n^6}{288\pi^4} + \frac{(\gamma_N g_{\xi} k_n^3 + \gamma_{\chi} y_{\xi} k_{\chi}^3)^2}{72\pi^4 m_{\xi}^2} \\ \pm \frac{1}{2} m_{\phi}^2 \phi^2. \end{aligned} \quad (A11) $$

## APPENDIX B: MACROSCOPIC PROPERTIES OF NEUTRON STARS

The mass-radius ($M-R$) relation of neutron stars is obtained by numerically solving the TOV equations, expressed here in gravitational units:

$$ \begin{aligned} \frac{d\mathcal{P}}{dr} &= -\frac{(\mathcal{E} + \mathcal{P})(M + 4\pi r^3 \mathcal{P})}{r(r - 2M)}, \\ \frac{dM}{dr} &= 4\pi r^2 \mathcal{E}, \end{aligned} \quad (B1) $$

where $M(r)$ is the enclosed mass within radius $r$, $\mathcal{P}(r)$ is the pressure, and $\mathcal{E}(r)$ is the energy density. These equations are integrated simultaneously with the following boundary conditions: at the center of the star ($r = 0$), the pressure and mass
satisfy $\mathcal{P}(0) = \mathcal{P}_c$ and $\mathcal{M}(0) = 0$, where $\mathcal{P}_c$ is the central pressure (or equivalently, $\mathcal{E}_c$ is the central energy density). The surface of the star is defined by the condition $\mathcal{P}(R) = 0$, which determines the stellar radius $R$. The stellar mass is obtained from $\mathcal{M}(r=R) = M$.

The compactness of a neutron star (NS) is a dimensionless quantity, defined as

$$C = \frac{GM}{c^2 R}, \qquad (B2)$$

where $G$ is Newton's gravitational constant, $c$ is the speed of light in vacuum, $M$ is the mass of the star, and $R$ is its radius. Since this is a dimensionless quantity, its numerical value remains the same regardless of the unit system. To obtain the compactness at an arbitrary radius $r$ within the star, one replaces $M \to \mathcal{M}(r)$ and $R \to r$ in the above expression.

To calculate the MoI of a NS, one must take into account its rotation. In a rotating NS, the centrifugal force experienced by a fluid element depends not only on the angular velocity $\Omega$ of the star but also on the local frame-dragging angular velocity $\omega(r)$ at a radial coordinate $r$. The relative angular velocity is defined as

$$\bar{\omega}(r) = \Omega - \omega(r). \qquad (B3)$$

In the slow-rotation approximation, $\bar{\omega}(r)$ satisfies the following second-order ordinary differential equation (ODE) [102]:

$$\frac{1}{r^4} \frac{d}{dr} \left( r^4 j(r) \frac{d\bar{\omega}}{dr} \right) + \frac{4}{r} \frac{d j(r)}{dr} \bar{\omega} = 0, \qquad (B4)$$

where $j(r) = e^{-\nu} \sqrt{1 - 2\mathcal{M}(r)/r}$ for $r < R$ and $j(r) = 1$ for $r \ge R$. Also, $\nu(r)$ satisfies the following ODE:

$$\frac{dv}{dr} = \frac{\mathcal{M}(r) + 4\pi r^3 \mathcal{P}(r)}{r[r - 2\mathcal{M}(r)]}. \qquad (B5)$$

The rotation frequency $\Omega(R)$ at the surface of the star is related to $\bar{\omega}(r)$ by the matching condition:

$$\Omega(R) = \bar{\omega}(R) + \frac{R}{3} \left( \frac{d\bar{\omega}}{dr} \right)_{r=R}. \qquad (B6)$$

To solve Eq. (B4), one imposes regular boundary conditions: $(d\bar{\omega}/dr)_{r=0} = 0$ and $\bar{\omega}(0) = \bar{\omega}_c$, where $\bar{\omega}_c$ is an arbitrarily chosen central value. Since Eq. (B4) is linear in $\bar{\omega}$, the solution can be scaled to match any desired $\Omega(R)$. The MoI can be calculated in the slow rotation approximation using the following expression [67]:

$$I = \frac{c^2}{6G} \left[ \frac{R^4}{\Omega(R)} \left( \frac{d\bar{\omega}}{dr} \right)_{r=R} \right]. \qquad (B7)$$

Although there is arbitrariness in the choice of the initial value $\bar{\omega}_c$, this does not affect the final value of $I$ in the slow rotation approximation, since both $\Omega(R)$ and $(d\bar{\omega}/dr)_{r=R}$ scale linearly with $\bar{\omega}_c$, and their ratio remains invariant.

Starting from the following unperturbed spherically symmetric, static metric,

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

one can perform a perturbative analysis on it to calculate $\Lambda$. The perturbed metric component $y(r)$ satisfies the following

first-order nonlinear ODE:

$$r \frac{dy}{dr} + y^2 + ye^{\lambda}[1 + 4\pi r^2(\mathcal{P} - \mathcal{E})] + r^2 \mathcal{Q} = 0, \qquad (B9)$$

where

$$\mathcal{Q}(r) = 4\pi e^{\lambda} \left( 5\mathcal{E} + 9\mathcal{P} + \frac{\mathcal{E} + \mathcal{P}}{C_S^2} \right) - \frac{6e^{\lambda}}{r^2} - \left( \frac{dv}{dr} \right)^2, \qquad (B10)$$

Using the boundary condition $y(0) = 2$, we calculate the tidal Love number $k_2$ [103,104] as

$$k_2 = \frac{8}{5} C^5 (1 - 2C)^2 [2 - y_R + 2C(y_R - 1)] \\ \times [2C\{6 - 3y_R + 3C \\ \times (5y_R - 8)\} + 4C^3\{13 - 11y_R + C(3y_R - 2) \\ + 2C^2(1 + y_R)\} \\ + 3(1 - 2C)^2\{2 - y_R + 2C(y_R - 1)\} \ln(1 - 2C)]^{-1}. \qquad (B11)$$

Here, $y_R = y(r=R)$ is the value of the perturbed metric element at the surface of the star and $\Lambda$ is given by [105]

$$\Lambda = \frac{2}{3} k_2 C^{-5} = \frac{2}{3} k_2 \left( \frac{Rc^2}{GM} \right)^5. \qquad (B12)$$

For the static NS, one solves the set of equations (B1), (B5), and (B9). However, in the slowly rotating NS case, one solves Eqs. (B1), (B4), (B5), and (B9) simultaneously to calculate all the relevant observables.

## APPENDIX C: MEASURES OF CONFORMALITY

Classical massless QCD turns out to exhibit conformal symmetry. Such a symmetry is characterized by the matter part of the trace anomaly [60] $\langle\Theta\rangle = \langle T^\mu_\mu \rangle = \mathcal{E} - 3\mathcal{P} = 0$, where $\mathcal{T}$ is the stress energy tensor of the model. So, fully conformal matter is characterized by the EoS $\mathcal{P} = \mathcal{E}/3$. However, quantizing the theory breaks this symmetry at the loop level due to finite masses of quarks. This deviation from conformality is measured by the deviation of $\langle\Theta\rangle$ from zero. As proposed in Ref. [60], however, we use

$$\Delta = \frac{1}{3} - \frac{\mathcal{P}}{\mathcal{E}} = \frac{\langle\Theta\rangle}{3\mathcal{E}}, \qquad (C1)$$

as a measure of trace anomaly. Thermodynamic stability ($\mathcal{P} \ge 0$) and causality ($0 \le C_S^2 < 1$) require that $-2/3 < \Delta < 1/3$. The average $\langle C_S^2 \rangle$, in an energy density interval $[0, \mathcal{E}]$ is defined as

$$\langle C_S^2 \rangle = \frac{1}{\mathcal{E}} \int_0^\mathcal{E} d\mathcal{E}' \frac{d\mathcal{P}}{d\mathcal{E}'} = \frac{\mathcal{P}(\mathcal{E})}{\mathcal{E}}. \qquad (C2)$$

This assumes a uniform probability distribution of $\mathcal{E}$. This leads to an alternate expression of $\langle C_S^2 \rangle$ as follows $\langle C_S^2 \rangle = 1/3 - \Delta$. There are two different representations of $C_S^2$ that one can use in the present context. One can decompose $C_S^2$ as a sum of slope ($\alpha$) and curvature ($\beta$) of energy per particle
$\mathcal{E}/\rho$ as follows [61]

$$C_S^2 = \frac{1}{\mu} \frac{dP}{d\rho} = \frac{2\rho}{\mu} \frac{d(\mathcal{E}/\rho)}{d\rho} + \frac{\rho^2}{\mu} \frac{d^2(\mathcal{E}/\rho)}{d\rho^2} = \alpha + \beta, \quad (C3)$$

where $\mu$ is the baryon chemical potential. The quantities $\alpha$ and $\beta$ are defined as

$$\alpha = \frac{2\rho}{\mu} \frac{d(\mathcal{E}/\rho)}{d\rho} = \frac{2C_S^2}{C_S^2 + \gamma} = 2\frac{1/3 - \Delta}{4/3 - \Delta}, \quad (C4)$$

and

$$\beta = \frac{\rho^2}{\mu} \frac{d^2(\mathcal{E}/\rho)}{d\rho^2} = C_S^2 - \alpha. \quad (C5)$$

The polytropic index $\gamma$ is defined as

$$\gamma = \frac{\mathcal{E}}{P} C_S^2 = \frac{C_S^2}{\langle C_S^2 \rangle} = \frac{C_S^2}{1/3 - \Delta}. \quad (C6)$$

An alternate representation of $C_S^2$ is as follows:

$$C_S^2 = \frac{1}{3} - \Delta - \Delta', \quad (C7)$$

[1] R. Somasundaram, I. Tews, and J. Margueron, Perturbative QCD and the neutron star equation of state, Phys. Rev. C 107, L052801 (2023).

[2] O. Philipsen, The QCD equation of state from the lattice, Prog. Part. Nucl. Phys. 70, 55 (2013).

[3] C. Drischler, J. W. Holt, and C. Wellenhofer, Chiral effective field theory and the high-density nuclear equation of state, Annu. Rev. Nucl. Part. Sci. 71, 403 (2021).

[4] F. Özel, D. Psaltis, T. Güver, G. Baym, C. Heinke, and S. Guillot, The dense matter equation of state from neutron star radius and mass measurements, Astrophys. J. 820, 28 (2016).

[5] V. Doroshenko, V. Suleimanov, G. Pühlhofer, and A. Santangelo, A strangely light neutron star within a supernova remnant, Nat. Astron. 6, 1444 (2022).

[6] E. Vurgun, M. Linares, S. Ransom, A. Papitto, S. Bogdanov, E. Bozzo, N. Rea, D. Garcia-Senz, P. Freire, and I. Stairs, The neutron star population in M28: A joint Chandra/GBT look at pulsar paradise, Astrophys. J. 941, 76 (2022).

[7] J.-L. Jiang, S.-P. Tang, Y.-Z. Wang, Y.-Z. Fan, and D.-M. Wei, PSR J0030+0451, GW170817 and the nuclear data: Joint constraints on equation of state and bulk properties of neutron stars, Astrophys. J. 892, 55 (2020).

[8] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink, and T. E. Strohmayer, A NICER view of PSR J0030+0451: Millisecond pulsar parameter estimation, Astrophys. J. Lett. 887, L21 (2019).

[9] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W.

where

$$\Delta' = \mathcal{E} \frac{d\Delta}{d\mathcal{E}} = \langle C_S^2 \rangle - C_S^2. \quad (C8)$$

The quantities $\Delta$ and $\Delta'$ are called nonderivative and derivative contributions to the speed of sound, respectively. The conditions $-2/3 < \Delta < 1/3$ and $C_S^2 \in [0, 1]$ lead to $\alpha \in [0, 1]$ and $\beta \in [-1, 1]$. At very low values of $\rho^V$, conformal symmetry is broken, i.e., $C_S^2 \to 0$ and $\Delta \to 1/3$. This implies that $\alpha \to 0$ and $\beta \to C_S^2$. On the other hand, at asymptotically high $\rho^V$, conformal symmetry is restored, i.e., $\Delta \to 0$, and hence, $C_S^2 \to 1/3$ and $\langle C_S^2 \rangle \to 1/3$. This leads to $\alpha \to 1/2$ and $\beta \to -1/6$. Hence, in the asymptotically large baryon density, $\Delta \to 0$ and $\Delta' \to 0$ in the fully conformal limit. This implies that the conformal distance [59]

$$d_c = \sqrt{\Delta^2 + \Delta'^2} = \sqrt{(1/3 - \langle C_S^2 \rangle)^2 + (C_S^2 - \langle C_S^2 \rangle)^2} \quad (C9)$$

vanishes in the fully conformal limit. However, it is useful to ask for values of various parameters described above (and the corresponding baryon density around which such values are reached), which can characterize the onset of conformality. This could prove to be useful as intermediate densities where conformality sets in could be close to densities reached in NSs. It has been shown in the literature that $\Delta = 0$ [60], $\beta = 0$ [61], $\gamma \approx 1.75$ [58], and $d_c \approx 0.2$ [59] serve as fairly good threshold conditions for the onset of conformality.

C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin et al., PSR J0030+0451 Mass and radius from NICER data and implications for the properties of neutron star matter, Astrophys. J. Lett. 887, L24 (2019).

[10] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen et al., A massive pulsar in a compact relativistic binary, Science 340, 1233232 (2013).

[11] J.-L. Huo and X.-F. Zhao, The moment of inertia of the proto neutron star PSR J0348+0432 under neutrino trapped, Chin. J. Phys. 56, 292 (2018).

[12] X.-F. Zhao, On the moment of inertia of PSR J0348+0432, Chin. J. Phys. 54, 839 (2016).

[13] E. Fonseca, H. T. Cromartie, T. T. Pennucci, P. S. Ray, A. Yu. Kirichenko, S. M. Ransom, P. B. Demorest, I. H. Stairs, Z. Arzoumanian, L. Guillemot, A. Parthasarathy, M. Kerr, I. Cognard, P. T. Baker, H. Blumer, P. R. Brook, M. DeCesar, T. Dolch, F. A. Dong, E. C. Ferrara et al., Refined mass and geometric measurements of the high-mass PSR J0740+6620, Astrophys. J. Lett. 915, L12 (2021).

[14] T. E. Riley, A. L. Watts, P. S. Ray, S. Bogdanov, S. Guillot, S. M. Morsink, A. V. Bilous, Z. Arzoumanian, D. Choudhury, J. S. Deneva, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, M. Loewenstein, R. M. Ludlam, C. B. Markwardt, T. Okajima, C. Prescod-Weinstein, R. A. Remillard et al., A NICER view of the massive pulsar PSR J0740+6620 informed
by radio timing and XMM-Newton spectroscopy, Astrophys. J. Lett. 918, L27 (2021).

[15] Y. Li, J. Wang, Z. Wu, and D. Wen, Inferring the gravitational binding energy and moment of inertia of PSR J0030 + 0451 and PSR J0740 + 6620 from new universal relations, Class. Quantum Grav. 39, 035014 (2022).

[16] P. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Shapiro delay measurement of a two solar mass neutron star, Nature (London) 467, 1081 (2010).

[17] A. Majid and M. Sharif, Quark stars in massive Brans-Dicke gravity with Tolman-Kuchowicz spacetime, Universe 6, 124 (2020).

[18] D. Choudhury, T. Salmi, S. Vinciguerra, T. E. Riley, Y. Kini, A. L. Watts, B. Dorsman, S. Bogdanov, S. Guillot, P. S. Ray, D. J. Reardon, R. A. Remillard, A. V. Bilous, D. Huppenkothen, J. M. Lattimer, N. Rutherford, Z. Arzoumanian, K. C. Gendreau, S. M. Morsink, and W. C. G. Ho, A NICER view of the nearest and brightest millisecond pulsar: PSR J0437-4715, Astrophys. J. Lett. 971, L20 (2024).

[19] L. Qi, S. Zheng, J. Zhang, M. Ge, Ang Li, S.-N. Zhang, F. Lu, H. Peng, Liang Zhang, H. Feng, Z. Zhang, Y. Xu, Z. Li, L. Song, S. Zhang, L. Tao, and W. Ye, PSR J1231-1411 revisited: Pulse profile analysis of x-ray observation, ApJ 981, 99 (2025).

[20] A. G. Lynn, A review of the double pulsar - PSR J0737-3039, Chin. J. Astron. Astrophys. 6, 162 (2006).

[21] Y.-Y. Yang, L. Chen, R.-F. Linghu, L.-Y. Zhang, and A. Taani, Constraints on estimation of radius of double pulsar PSR J0737-3039A and its neutron star nuclear matter composition, Chin. Phys. Lett. 34, 129701 (2017).

[22] P. Landry and B. Kumar, Constraints on the moment of inertia of PSR J0737-3039A from GW170817, Astrophys. J. Lett. 868, L22 (2018).

[23] B. P. Abbott et al. (LIGO Scientific, Virgo), GW170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119, 161101 (2017).

[24] B. P. Abbott et al. (LIGO Scientific, Virgo), GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121, 161101 (2018).

[25] B. P. Abbott et al. (LIGO Scientific, Virgo), GW190425: Observation of a compact binary coalescence with total mass ~3.4M☉, Astrophys. J. Lett. 892, L3 (2020).

[26] M.-Z. Han, S.-P. Tang, Y.-M. Hu, Y.-J. Li, J.-L. Jiang, Z.-P. Jin, Y.-Z. Fan, and D.-M. Wei, Is GW190425 consistent with being a neutron star-black hole merger? Astrophys. J. Lett. 891, L5 (2020).

[27] F. Hofmann, C. M. Keil, and H. Lenske, Application of the density dependent Hadron field theory to neutron star matter, Phys. Rev. C 64, 025804 (2001).

[28] B. Y. Sun, W. H. Long, J. Meng, and U. Lombardo, Neutron star properties in density-dependent relativistic Hartree-Fock theory, Phys. Rev. C 78, 065805 (2008).

[29] K. Huang, H. Shen, J. Hu, and Y. Zhang, Density-dependent quark mean-field model for nuclear matter and neutron stars, Phys. Rev. C 109, 045804 (2024).

[30] B.-A. Li, B.-J. Cai, L.-W. Chen, and J. Xu, Nucleon effective masses in neutron-rich matter, Prog. Part. Nucl. Phys. 99, 29 (2018).

[31] T. Kojo, QCD equations of state and speed of sound in neutron stars, AAPPS Bull. 31, 11 (2021).

[32] S. Li, J. Pang, H. Shen, J. Hu, and K. Sumiyoshi, Influence of effective nucleon mass on equation of state for super-nova simulations and neutron stars, Astrophys. J. 980, 54 (2025).

[33] A. Y. Potekhin and G. Chabrier, Magnetic neutron star cooling and microphysics, Astron. Astrophys. 609, A74 (2018).

[34] A. Y. Potekhin, A. I. Chugunov, and G. Chabrier, Thermal evolution and quiescent emission of transiently accreting neutron stars, Astron. Astrophys. 629, A88 (2019).

[35] M. Baldo, G. F. Burgio, H.-J. Schulze, and G. Taranto, Nucleon effective masses within the Brückner-Hartree-Fock theory: Impact on stellar neutrino emission, Phys. Rev. C 89, 048801 (2014).

[36] A. Dehghan Niri, H. R. Moshfegh, and P. Haensel, Nuclear correlations and neutrino emissivity from the neutron branch of the modified Urca process, Phys. Rev. C 93, 045806 (2016).

[37] P. S. Shternin, M. Baldo, and P. Haensel, In-medium enhancement of the modified Urca neutrino reaction rates, Phys. Lett. B 786, 28 (2018).

[38] F. Sammarruca, Temperature-dependence of single-particle properties in isospin-symmetric and asymmetric matter within the Dirac-Brückner-Hartree-Fock model, J. Phys. G 37, 085105 (2010).

[39] W. Zuo, A. Lejeune, U. Lombardo, and J. F. Mathiot, Interplay of three-body interactions in the EOS of nuclear matter, Nucl. Phys. A 706, 418 (2002).

[40] W. Zuo, A. Lejeune, U. Lombardo, and J. F. Mathiot, Microscopic three-body force for asymmetric nuclear matter, Eur. Phys. J. A 14, 469 (2002).

[41] A. W. Steiner and S. Gandolfi, Connecting neutron star observations to three-body forces in neutron matter and to the nuclear symmetry energy, Phys. Rev. Lett. 108, 081102 (2012).

[42] S. Gandolfi, J. Carlson, S. Reddy, A. W. Steiner, and R. B. Wiringa, The equation of state of neutron matter, symmetry energy, and neutron star structure, Eur. Phys. J. A 50, 10 (2014).

[43] Y. Yamamoto, H. Togashi, T. Tamagawa, T. Furumoto, N. Yasutake, and Th. A. Rijken, Neutron-star radii based on realistic nuclear interactions, Phys. Rev. C 96, 065804 (2017).

[44] B. T. Reed, M. Heinz, P. Arthuis, A. Schwenk, and I. Tews, Connecting relativistic density functional theory to microscopic calculations, Phys. Rev. C 112, 034331 (2025).

[45] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations, Astrophys. J. 860, 149 (2018).

[46] A. Cherman, T. D. Cohen, and A. Nellore, A bound on the speed of sound from holography, Phys. Rev. D 80, 066003 (2009).

[47] C. Ecker, C. Hoyos, N. Jokela, D. Rodríguez Fernández, and A. Vuorinen, Stiff phases in strongly coupled gauge theories with holographic duals, J. High Energy Phys. 11 (2017) 031.

[48] S. Altiparmak, C. Ecker, and L. Rezzolla, On the sound speed in neutron stars, Astrophys. J. Lett. 939, L34 (2022).

[49] L. Brandes, W. Weise, and N. Kaiser, Inference of the sound speed and related properties of neutron stars, Phys. Rev. D 107, 014011 (2023).
[50] S. Roy and T. Suyama, On the sound velocity bound in neutron stars, *Results Phys.* **61**, 107757 (2024).

[51] P. Bedaque and A. W. Steiner, Sound velocity bound and neutron stars, *Phys. Rev. Lett.* **114**, 031103 (2015).

[52] L. McLerran and S. Reddy, Quarkyonic matter and neutron stars, *Phys. Rev. Lett.* **122**, 122701 (2019).

[53] M. Leonhardt, M. Pospiech, B. Schallmo, J. Braun, C. Drischler, K. Hebeler, and A. Schwenk, Symmetric nuclear matter from the strong interaction, *Phys. Rev. Lett.* **125**, 142502 (2020).

[54] N. Kovensky, A. Poole, and A. Schmitt, Building a realistic neutron star from holography, *Phys. Rev. D* **105**, 034022 (2022).

[55] Ch. C. Moustakidis, T. Gaitanos, Ch. Margaritis, and G. A. Lalazissis, Bounds on the speed of sound in dense matter, and neutron star structure, *Phys. Rev. C* **95**, 045801 (2017); **95**, 059904(E) (2017).

[56] C. Margaritis, P. K. Koutmiridis, and C. Moustakidis, Speed of sound constraints on maximally-rotating neutron stars, *Phys. Rev. C* **27**, 155 (2020).

[57] T. S. Olson, Maximally incompressible neutron star matter, *Phys. Rev. C* **63**, 015802 (2000).

[58] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Evidence for quark-matter cores in massive neutron stars, *Nat. Phys.* **16**, 907 (2020).

[59] E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, J. Nättilä, and A. Vuorinen, Strongly interacting matter exhibits deconfined behavior in massive neutron stars, *Nat. Commun.* **14**, 8451 (2023).

[60] Y. Fujimoto, K. Fukushima, L. D. McLerran, and M. Praszalowicz, Trace anomaly as signature of conformality in neutron stars, *Phys. Rev. Lett.* **129**, 252702 (2022).

[61] Michał Marczenko, K. Redlich, and C. Sasaki, Curvature of the energy per particle in neutron stars, *Phys. Rev. D* **109**, L041302 (2024).

[62] D. Scordino and I. Bombaci, Dark matter admixed neutron stars with a realistic nuclear equation of state from chiral nuclear interactions, *JHEAp* **45**, 371 (2025).

[63] P. K. Sahu, R. Basu, and B. Datta, High density matter in the chiral sigma model, *Astrophys. J.* **416**, 267 (1993).

[64] P. K. Sahu and A. Ohnishi, SU(2) chiral sigma model and the properties of neutron stars, *Prog. Theor. Phys.* **104**, 1163 (2000).

[65] J. D. Walecka, A theory of highly condensed matter, *Ann. Phys. (NY)* **83**, 491 (1974).

[66] A. Guha and D. Sen, Feeble DM-SM interaction via new scalar and vector mediators in rotating neutron stars, *J. Cosmol. Astropart. Phys.* 09 (2021) 027.

[67] N. K. Glendenning, *Compact Stars: Nuclear Physics, Particle Physics, and General Relativity* (Springer, New York, NY, 1997).

[68] T. K. Jha and H. Mishra, Constraints on nuclear matter parameters of an effective chiral model, *Phys. Rev. C* **78**, 065802 (2008).

[69] M. Dutra, O. Lourenço, S. S. Avancini, B. V. Carlson, A. Delfino, D. P. Menezes, C. Providência, S. Typel, and J. R. Stone, Relativistic mean-field hadronic models under nuclear matter constraints, *Phys. Rev. C* **90**, 055203 (2014).

[70] S. Tulin, H.-B. Yu, and K. M. Zurek, Beyond collisionless dark matter: Particle physics dynamics for dark matter halo structure, *Phys. Rev. D* **87**, 115007 (2013).

[71] S. W. Randall, M. Markevitch, D. Clowe, A. H. Gonzalez, and M. Bradac, Constraints on the self-interaction cross-section of dark matter from numerical simulations of the merging galaxy cluster 1E 0657-56, *Astrophys. J.* **679**, 1173 (2008).

[72] D. Sen and A. Guha, Implications of feebly interacting dark sector on neutron star properties and constraints from GW170817, *Mon. Not. R. Astron. Soc.* **504**, 3354 (2021).

[73] X. L. Shang, A. Li, Z. Q. Miao, G. F. Burgio, and H. J. Schulze, Nucleon effective mass in hot dense matter, *Phys. Rev. C* **101**, 065801 (2020).

[74] A. G. Abac, C. C. Bernido, and Jose Perico H. Esguerra, Stability of neutron stars with dark matter core using three crustal types and the impact on mass-radius relations, *Phys. Dark Universe* **40**, 101185 (2023).

[75] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Gravitational-wave constraints on the neutron-star-matter equation of state, *Phys. Rev. Lett.* **120**, 172703 (2018).

[76] C.-J. Xia, T. Maruyama, A. Li, B. Y. Sun, W.-H. Long, and Y.-X. Zhang, Unified neutron star EOSs and neutron star structures in RMF models, *Commun. Theor. Phys.* **74**, 095303 (2022).

[77] V. Sagun, E. Giangrandi, T. Dietrich, O. Ivanytskyi, R. Negreiros, and C. Providência, What is the nature of the HESS J1731-347 compact object? *Astrophys. J.* **958**, 49 (2023).

[78] Q.-F. Xiang, W.-Z. Jiang, D.-R. Zhang, and R.-Y. Yang, Effects of fermionic dark matter on properties of neutron stars, *Phys. Rev. C* **89**, 025803 (2014).

[79] H. C. Das, A. Kumar, B. Kumar, and S. K. Patra, Dark matter effects on the compact star properties, *Galaxies (Basel)* **10**, 14 (2022).

[80] César H. Lenzi, M. Dutra, O. Lourenço, L. L. Lopes, and Débora P. Menezes, Dark matter effects on hybrid star properties, *Eur. Phys. J. C* **83**, 266 (2023).

[81] Central Bureau for Astronomical Telegrams, IAUC 8144: XTE J1814-338; 2003ed, http://www.cbat.eps.harvard.edu/iauc/08100/08144.html#Item1.

[82] S. L. Pitz and J. Schaffner-Bielich, Generating ultracompact neutron stars with bosonic dark matter, *Phys. Rev. D* **111**, 043050 (2025).

[83] L. L. Lopes and A. Issifu, XTE J1814-338 as a dark matter admixed neutron star, *Phys. Dark Universe* **48**, 101922 (2025).

[84] S.-H. Yang, C.-M. Pi, and F. Weber, Strange stars admixed with mirror dark matter: Confronting observations of XTE J1814-338, *Phys. Rev. D* **111**, 043037 (2025).

[85] P. Laskos-Patkos and Ch. C. Moustakidis, XTE J1814-338: A potential hybrid star candidate, *Phys. Rev. D* **111**, 063058 (2025).

[86] S. Traversi, P. Char, G. Pagliara, and A. Drago, Speed of sound in dense matter and two families of compact stars, *Astron. Astrophys.* **660**, A62 (2022).

[87] Z. Miao, A. Li, and Z.-G. Dai, On the moment of inertia of PSR J0737-3039 A from LIGO/Virgo and NICER, *Mon. Not. R. Astron. Soc.* **515**, 5071 (2022).

[88] S. Chatterjee, H. Sudhakaran, and R. Mallick, Analyzing the speed of sound in neutron star with machine learning, *Eur. Phys. J. C* **84**, 1291 (2024).

[89] Michał Marczenko, Conformality thresholds in neutron stars, *J. Subat. Part. Cosmol.* **3**, 100043 (2025).

[90] Michał Marczenko, L. McLerran, K. Redlich, and C. Sasaki, Reaching percolation and conformal limits in neutron stars, *Phys. Rev. C* **107**, 025802 (2023).
[91] R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, PSR J0952-0607: The fastest and heaviest known galactic neutron star, Astrophys. J. Lett. 934, L17 (2022).

[92] T. Gorda, A. Kurkela, P. Romatschke, S. Säppi, and A. Vuorinen, Next-to-next-to-next-to-leading order pressure of cold quark matter: Leading logarithm, Phys. Rev. Lett. 121, 202701 (2018).

[93] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Soft interactions in cold quark matter, Phys. Rev. Lett. 127, 162003 (2021).

[94] O. Ivanytskyi and D. B. Blaschke, Recovering the conformal limit of color superconducting quark matter within a confining density functional approach, Particles 5, 514 (2022).

[95] R. Somasundaram, I. Tews, and Jérôme Margueron, Investigating signatures of phase transitions in neutron-star cores, Phys. Rev. C 107, 025801 (2023).

[96] T. Mandal, P. Jaikumar, and S. Digal, Chiral and diquark condensates at large magnetic field in two-flavor superconducting quark matter, arXiv:0912.1413.

[97] T. Mandal and P. Jaikumar, Neutrality of a magnetized two-flavor quark superconductor, Phys. Rev. C 87, 045208 (2013).

[98] T. Mandal and P. Jaikumar, Effect of temperature and magnetic field on two-flavor superconducting quark matter, Phys. Rev. D 94, 074016 (2016).

[99] T. Mandal and P. Jaikumar, Effect of strong magnetic field on competing order parameters in two-flavor dense quark matter, Adv. High Energy Phys. 2017, 1 (2017).

[100] P. Jaikumar, G. Rupak, and A. W. Steiner, Viscous damping of $r$-mode oscillations in compact stars with quark matter, Phys. Rev. D 78, 123007 (2008).

[101] P. Jaikumar, A. Semposki, M. Prakash, and C. Constantinou, $g$-mode oscillations in hybrid stars: A tale of two sounds, Phys. Rev. D 103, 123009 (2021).

[102] J. B. Hartle, Slowly rotating relativistic stars. I. Equations of structure, Astrophys. J. 150, 1005 (1967).

[103] T. Hinderer, Tidal love numbers of neutron stars, Astrophys. J. 677, 1216 (2008); 697, 964 (2009).

[104] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral, Phys. Rev. D 81, 123016 (2010).

[105] K.-L. Leung, M.-c. Chu, and L.-M. Lin, Tidal deformability of dark matter admixed neutron stars, Phys. Rev. D 105, 123010 (2022).