Neutron stars serve as natural laboratories for investigating the properties of matter under extreme densities and strong gravitational fields (Lattimer and Prakash, 2004). Understanding the properties of dense nuclear matter is essential for describing the structure and evolution of neutron stars (Burgio et al., 2021; Sedrakian et al., 2023). At the core of this pursuit lies the nuclear equation of state (EoS), which connects microscopic nuclear interactions to macroscopic observables such as neutron star masses and radii (Lattimer and Prakash, 2000; Lattimer and Prakash, 2007; Oertel et al., 2017; Huth et al., 2022). The EoS essentially encapsulates the relationship between pressure and density in nuclear matter, determining how matter behaves under the extreme conditions found in neutron star interiors. In particular, the recent detection of neutron stars with masses exceeding 2 solar masses and the advent of multi-messenger astronomy have placed stringent constraints on the EoS, highlighting the necessity of developing and refining theoretical models that are not only consistent with laboratory nuclear physics data but also aligned with the latest astrophysical observations (Demorest et al., 2010; Antoniadis et al., 2013; Fonseca et al., 2016; Arzoumanian et al., 2018; Cromartie et al., 2020; Fonseca et al., 2021; Abbott et al., 2017; Abbott et al., 2018; Tong et al., 2020; Han et al., 2023).

Over the years, considerable theoretical efforts have been devoted to determining the EoS of neutron star matter using various nuclear many-body approaches. In general, these approaches can be categorized into two classes: density functional theories (DFTs) employing effective nucleon-nucleon ( N N ) interactions, and ab initio methods based on realistic interactions. The effective N N interactions in DFTs, either non-relativistic or relativistic, are tuned to reproduce the properties of finite nuclei and nuclear matter around saturation density, as in the Skyrme (1956), Dechargé and Gogny (1980), and relativistic mean-field (RMF) models (Ring, 1996; Meng et al., 2006; Meng, 2016). However, due to weak constraints on the isovector channels, their predictions for nuclear matter properties, such as the nuclear symmetry energy at higher densities, remain uncertain (Li et al., 2008). In contrast, ab initio methods based on realistic N N interactions stand out for their predictive power, free from uncertainties associated with adjustable parameters. In recent years, a growing variety of ab initio methods have been developed for nuclear many-body studies, including quantum Monte Carlo (Carlson et al., 2015), coupled-cluster (Hagen et al., 2014), no-core shell model (Barrett et al., 2013), self-consistent Green’s function (Dickhoff and Barbieri, 2004), lattice effective field theory (Lee, 2009; Lähde and Meißner, 2019; Tong et al., 2025a; Tong et al., 2025b; Tong et al., 2025c), in-medium similarity renormalization group (Hergert et al., 2016), Monte Carlo shell model (Otsuka et al., 2001; Liu et al., 2012), and Brueckner-Hartree-Fock (BHF) theory (Shang et al., 2021). Among these, the relativistic Brueckner-Hartree-Fock (RBHF) theory stands out as one of the most successful ab initio methods based solely on bare two-body forces. Benefiting from the relativistic framework–which is essential at high densities due to the crucial role of Lorentz covariance, the RBHF theory has been successfully applied to both finite nuclei (Shen et al., 2016; Shen et al., 2019; Wang et al., 2019) and dense matter systems such as nuclear matter and neutron stars (Brockmann and Machleidt, 1990; Sehn et al., 1997; de Jong and Lenske, 1998; van Dalen et al., 2005; Katayama and Saito, 2013; Tong et al., 2018; Wang et al., 2020).

The RBHF theory provides a self-consistent framework to study the nuclear many-body problem by combining Dirac phenomenology with the in-medium scattering equation. In this approach, the interaction between two nucleons in the nuclear medium is described by the in-medium scattering matrix G , obtained by summing ladder diagrams with a realistic N N potential. The effective single-particle potential is derived from G matrix, which in turn modifies the nucleon spinors via the Dirac equation, thereby closing the self-consistent loop. To simplify RBHF calculations, earlier studies adopted the average center of mass (c. m.) momentum approximation for computing the binding energy (Brueckner et al., 1968; Alonso and Sammarruca, 2003; Sammarruca et al., 2012; Sammarruca, 2014). With modern computational capabilities, this approximation can be avoided. Recent work derived exact analytic expressions for the angular integration over the c. m. momentum, with a focus on asymmetric nuclear matter (Tong et al., 2018). A significant contribution to the saturation properties was found when treating the total momentum exactly, underscoring its impact on higher-order quantities in both the energy of symmetric matter and the symmetry energy. Another key challenge in RBHF theory is the self-consistent extraction of the nucleon single-particle potential from the in-medium G matrix, where symmetry arguments dictate its decomposition into scalar and vector components (Serot and Walecka, 1986). Traditional approaches include the momentum-independence approximation (Brockmann and Machleidt, 1990), which neglects momentum dependence and fails to capture the correct isospin dependence of the single-particle potential in asymmetric nuclear matter (Ulrych and Müther, 1997; Schiller and Müther, 2001), and the projection method (Horowitz and Serot, 1987; Nuppenau et al., 1989; Gross-Boelting et al., 1999), which retains momentum dependence but is limited to positive-energy states (PESs). Notably, these methods yield contradictory predictions for the isospin dependence of the single-particle potential (Ulrych and Müther, 1997). Recently, a fully self-consistent RBHF framework in the full Dirac space has been developed (Wang et al., 2021), where the Lorentz structure and momentum dependence are determined without approximations. This advance resolves the long-standing discrepancy and provides a unique description of isospin effects in nuclear matter. As a result, the RBHF theory in the full Dirac space has been successfully and systematically applied to diverse nuclear systems, including the nuclear matter (Wang et al., 2022a; Wang et al., 2022b; Qu et al., 2023; Wang et al., 2023; Wang et al., 2024; Qin et al., 2025; Huang et al., 2025), the properties of ²⁰⁸Pb with a liquid droplet model (Tong et al., 2023), neutron star properties (Tong et al., 2022; Wang et al., 2022c; Qu et al., 2025; Laskos-Patkos et al., 2025), optical potential for proton-nucleus scattering (Qin et al., 2024), and in-medium nucleon-nucleon cross sections (Wang et al., 2025).

In this review, we summarize these recent advances in RBHF theory formulated in the full Dirac space and their implications for the physics of dense matter and neutron stars.

In the RBHF theory, nucleons within the nuclear medium are treated as dressed particles due to their interactions with surrounding nucleons. The single-particle motion of these nucleons is described by the Dirac equation α ⋅ p + β M + U u p , s = E p u p , s , (1) where α and β are the Dirac matrices, M is the nucleon mass, p and E p are the momentum and the single-particle energy, and s denotes the spin. According to the translational and rotational invariance, time-reversal invariance, hermiticity, and parity conservation, the single-particle potential U can be decomposed in its Lorentz form (Serot and Walecka, 1986) U p = U S p + γ 0 U 0 p + γ ⋅ p ̂ U V p . (2)

The quantities U S ( p ) , U 0 ( p ) , and U V ( p ) are the scalar potential, timelike, and spacelike parts of the vector potential respectively with p = | p | the magnitude of nucleon momentum. p ̂ = p / | p | is the unit vector. By using the following effective quantities in Equations 3a–3c: p * = p + p ̂ U V p , (3a) M p * = M + U S p , (3b) E p * = E p − U 0 p , (3c) the solution of Equation 1 leads to the in-medium positive-energy spinor u ( p , s ) and negative-energy spinor v ( p , s ) u p , s = E p * + M p * 2 M p * 1 σ ⋅ p * E p * + M p * χ s , (4a) v p , s = γ 5 u p , s = E p * + M p * 2 M p * σ ⋅ p * E p * + M p * 1 χ s , (4b) where χ s is the spin wave function.

The Dirac equation can be solved exactly once the single-particle potentials are determined. To this end, three matrix elements of U ( p ) in the full Dirac space are introduced, Σ + + p = u ̄ p , 1 / 2 U p u p , 1 / 2 = U S p + E p * M p * U 0 p + p * M p * U V p , (5a) Σ − + p = v ̄ p , 1 / 2 U p u p , 1 / 2 = p * M p * U 0 p + E p * M p * U V p , (5b) Σ − − p = v ̄ p , 1 / 2 U p v p , 1 / 2 = − U S p + E p * M p * U 0 p + p * M p * U V p . (5c)

After obtaining Σ + + ( p ) , Σ − + ( p ) , and Σ − − ( p ) , single-particle potentials in Equation 2 can be determined uniquely through U S p = Σ + + p − Σ − − p 2 , (6a) U 0 p = E p * M p * Σ + + p + Σ − − p 2 − p * M p * Σ − + p , (6b) U V p = − p * M p * Σ + + p + Σ − − p 2 + E p * M p * Σ − + p . (6c)

This approach avoids approximations in the Dirac space with PESs only. The matrix elements Σ + + ( p ) , Σ − + ( p ) , and Σ − − ( p ) can be calculated alternatively by summing up the effective two-body interaction G matrix with all the nucleons inside the Fermi sea in the Hartree-Fock approximation Σ + + p = ∑ s ′ ∫ 0 k F d 3 p ′ 2 π 3 M p ′ * E p ′ * ⟨ u ̄ p , 1 / 2 u ̄ p ′ , s ′ | G ̄ + + + + W | u p , 1 / 2 u p ′ , s ′ ⟩ , (7a) Σ − + p = ∑ s ′ ∫ 0 k F d 3 p ′ 2 π 3 M p ′ * E p ′ * ⟨ v ̄ p , 1 / 2 u ̄ p ′ , s ′ | G ̄ − + + + W | u p , 1 / 2 u p ′ , s ′ ⟩ , (7b) Σ − − p = ∑ s ′ ∫ 0 k F d 3 p ′ 2 π 3 M p ′ * E p ′ * ⟨ v ̄ p , 1 / 2 u ̄ p ′ , s ′ | G ̄ − + − + W | v p , 1 / 2 u p ′ , s ′ ⟩ . (7c)

In Equations 7a–7c, the anti-symmetrized G matrix is expressed with G ̄ , where the ± -signs in the superscript denote the positive- or negative-energy states. W is the starting energy which equals to the total single-particle energies in the initial states.

The G matrix is obtained by solving the in-medium Thompson equation (Brockmann and Machleidt, 1990) G q ′ , q | P , W = V q ′ , q | P + ∫ d 3 k 2 π 3 V q ′ , k | P × M P + k * M P − k * E P + k * E P − k * Q k , P W − E P + k − E P − k + i ϵ G k , q | P , W , (8) where P = 1 2 ( k 1 + k 2 ) is half the total momentum and k = 1 2 ( k 1 − k 2 ) is the relative momentum of the two interacting nucleons with momenta k 1 and k 2 . The initial, intermediate, and final relative momenta of the two nucleons scattering in nuclear matter are denoted by q , k , and q ′ , respectively. V is derived from a realistic N N interaction. Here we introduced the one-boson-exchange potential (OBEP) as an example (Machleidt, 1989), which has been well constrained by experimental data on both N N bound states and scattering observables. The N N scattering in the nuclear medium is restricted with the Pauli operator in Equation 9: Q k , P = 1 , | P + k | , | P − k | > k F , 0 , otherwise . (9)

The inclusion of an infinitesimal i ϵ term in the denominator is necessary when the starting energy W = E P + k + E P − k approaches a pole within the continuous choice. The RBHF equations were solved self-consistently for symmetric nuclear matter in the full Dirac space within the continuous choice for the single-particle potential (Wang et al., 2022a). The resulting single-particle energy and Dirac mass exhibit smooth continuity across the Fermi surface. Equation 1, Equations 6a–6c, Equations 7a–7c, Equation 8 constitute a coupled system that has to be solved in a self-consistent way.

After the solution of G matrix and the calculation of single-particle potentials converge, the binding energy per nucleon in nuclear matter can be calculated using E / A = 1 ρ ∑ s ∫ 0 k F d 3 p 2 π 3 M p * E p * ⟨ u ̄ p , s | γ ⋅ p + M | u p , s ⟩ − M + 1 2 ρ ∑ s , s ′ ∫ 0 k F d 3 p 2 π 3 ∫ 0 k F ′ d 3 p ′ 2 π 3 M p * E p * M p ′ * E p ′ * × ⟨ u ̄ p , s u ̄ p ′ , s ′ | G ̄ W | u p , s u p ′ , s ′ ⟩ , (10) where ρ denotes the total density of nucleons.

It should also be noted that the calculation of the binding energy yields a three dimension integrals over the c. m. momentum P and relative momentum q . The three-dimensional integrals over the c. m. momentum P are numerically challenging. Normally, the widely used averaged c. m. momentum approximation in Equation 11 is adopted (Brueckner et al., 1968; Alonso and Sammarruca, 2003), P av 2 = ∫ 0 k F d 3 k 1 ∫ 0 k F d 3 k 2 P 2 δ q − 1 2 | k 1 − k 2 | ∫ 0 k F d 3 k 1 ∫ 0 k F d 3 k 2 δ q − 1 2 | k 1 − k 2 | . (11)

It does not depend on the direction and this value is usually applied in the G -matrix in Equation 10. By relieving this approximation and deriving the exact expressions of the angular integrations of the momentum P within RBHF theory, the exact results has been calculated in Ref. Tong et al. (2018), especially for the case of asymmetric nuclear matter. For the calculations in the full Dirac space, the exact evaluation of P and its integral during the iteration are also performed.

One of the motivations for developing a microscopic and fully relativistic theory of dense nuclear matter is its application to neutron star. The neutron star matter here is assumed to be composed of nucleons and leptons (mainly electrons and muons), while neglecting possible phase transitions or the appearance of exotic degrees of freedom at densities above nuclear saturation. The matter is considered to be in beta equilibrium and charge neutrality, leading to the following equilibrium conditions for the chemical potentials of the nucleons and leptons in Equation 12: μ p = μ n − μ e , μ μ = μ e , (12) where μ e , μ μ , μ p , and μ n denote the chemical potentials of electrons, muons, protons, and neutrons, respectively. Charge neutrality is maintained in Equation 13: ρ p = ρ e + ρ μ , (13) where ρ p , ρ e , and ρ μ are the number densities of protons, electrons, and muons, respectively. The energy density of the beta equilibrium nuclear matter is then obtained as ε = ρ E ρ , α / A + Y p M p + 1 − Y p M n + ε e + ε μ , (14) where Y i = ρ i / ρ   ( i = e ,   μ ,   p ,   n ) are the equilibrium particle fractions. The chemical potential for each particle i is given in Equation 15: μ i = ∂ ε / ρ ∂ Y i . (15)

For a given density ρ , the particle fractions Y i are determined by solving the equilibrium conditions and charge neutrality, allowing the calculation of the energy density ε using Equation 14. The pressure P is then derived from Equation 16: P = − ∂ ε / ρ ∂ 1 / ρ = ρ ∂ ε ∂ ρ − ε . (16)

This yields the EoS of beta equilibrium nuclear matter in the form of P ( ε ) .

Once the EoS in the form P ( ε ) is obtained, the mass and radius of a cold, spherically symmetric, static, and relativistic star can be described by the Tolman-Oppenheimer-Volkov (TOV) equations (Oppenheimer and Volkoff, 1939; Tolman, 1939), d P r d r = − P r + ε r M r + 4 π r 3 P r r r − 2 M r , (17a) d M r d r = 4 π r 2 ε r , (17b) where P ( r ) is the pressure at neutron star radius r , M ( r ) is the total neutron star mass inside a sphere of radius r . Besides the masses and radii, another key property of neutron stars is their dimensionless tidal deformability (Damour et al., 1992; Hinderer, 2008; Flanagan and Hinderer, 2008), denoted by Λ . This quantity characterizes the star’s response to an external tidal field and is defined in Equation 18: Λ = 2 3 k 2 C − 5 . (18) C = M / R is the compactness parameter, where M is the neutron star mass and R is the radius and they are determined from the following two conditions: P ( R ) = 0 and M = M ( R ) . k 2 is the second love number quantifying the tidal response of the star in Equation 19, k 2 = 8 C 5 5 1 − 2 C 2 2 − y R + 2 C y R − 1 × 6 C 2 − y R + C 5 y R − 8 + 4 C 3 13 − 11 y R + C 3 y R − 2 + 2 C 2 1 + y R + 3 1 − 2 C 2 2 − y R + 2 C y R − 1 ln 1 − 2 C − 1 , (19) where y R = y ( R ) characterizes the response of the metric perturbation to the external tidal field at the stellar surface, and it can be calculated by solving the following differential equation, r d y r d r + y 2 r + y r F r + r 2 Q r = 0 , (20) with F r = 1 − 2 M r r − 1 1 − 4 π r 2 ε r − P r , (21a) Q r = 4 π 5 ε r + 9 P r + ε r + P r ∂ P ∂ ε r − 6 r 2 × 1 − 2 M r r − 1    − 2 M r r 2 + 2 × 4 π r P r 2 × 1 − 2 M r r − 2 . (21b)

The differential Equation 20 can be integrated together with the TOV equations with the boundary condition y ( 0 ) = 2 . In addition to tidal deformability, the rotational properties of neutron stars also provide crucial insights into their internal structure. The moment of inertia is calculated within the slow-rotation approximation (Hartle, 1967; Hartle and Thorne, 1968), where the frequency Ω of a uniformly rotating neutron star is significantly lower than the Kepler frequency at the equator. In this approximation, the moment of inertia I of a uniformly rotating, axially symmetric neutron star is given in Equation 22; Fattoyev and Piekarewicz (2010). I = 8 π 3 ∫ 0 R r 4 e − ν r ω ̄ r Ω ϵ r + P r 1 − 2 M r / r d r . (22)

Here, ν ( r ) is a radially-dependent metric function defined in Equation 23: ν r = 1 2 ln 1 − 2 M R − ∫ r R M x + 4 π x 3 P x x 2 1 − 2 M x / x d x . (23)

The frame-dragging angular velocity ω ̄ represents the angular velocity of the fluid as measured in a local inertial reference frame, which is usually expressed through the dimensionless relative frequency ω ̃ ≡ ω ̄ / Ω , which satisfies the following second-order differential equation in Equation 24: d d r r 4 j r d ω ̃ r d r + 4 r 3 d j r d r ω ̃ r = 0 , (24) where j ( r ) = e − ν ( r ) 1 − 2 M ( r ) / r for r ≤ R . The relative frequency ω ̃ ( r ) is subject to the boundary conditions ω ̃ ′ 0 = 0 ,    ω ̃ R + R 3 ω ̃ ′ R = 1 . (25)

It should be noted that, under the slow-rotation approximation, the moment of inertia is independent of the stellar frequency Ω .

The quadrupole moment characterizes the degree of rotational deformation of the neutron star away from spherical symmetry (Yagi and Yunes, 2013). It can be computed by numerically solving for the interior and exterior gravitational field of a neutron star in a slow-rotation (Hartle, 1967; Hartle and Thorne, 1968) and a small-tidal-deformation approximation (Hinderer, 2008; Hinderer et al., 2010). To explore the universal dimensionless moment of inertia-tidal deformability-quadrupole moment ( I –Love– Q ) relations, which are nearly independent of the EoS, we introduce the following quantities in Equation 26: I ̄ ≡ I M 3 , Q ̄ ≡ − Q M I Ω 2 . (26)

In addition, to describe the rapidly rotating and axisymmetric neutron star configurations in general relativity, the stellar matter is treated as a perfect fluid, characterized by the energy-momentum tensor in Equation 27: T μ ν = ε + P u μ u ν − g μ ν P , (27) where ε , P , and u μ are the energy density, pressure, and fluid’s four-velocity, respectively. The Einstein field equations are solved assuming an axisymmetric and stationary spacetime with the metric in Equation 28: d s 2 = − e γ + ρ d t 2 + e 2 α d r 2 + r 2 d θ 2 + e γ − ρ r 2 ⁡ sin 2 ⁡ θ d ϕ − ω d t 2 , (28) where the metric potentials γ , ρ , α , and ω are functions of the radial coordinates r and the polar angle θ . For numerical calculations, we utilize the RNS code (Stergioulas and Friedman, 1995; Paschalidis and Stergioulas, 2017) for rapidly rotating neutron stars.

In this review, we have focused on recent advances in the study of neutron star properties based on RBHF theory formulated in the full Dirac space.

Figure 1 from Ref. Qu et al. (2025) illustrates the gravitational mass of both static and rotating neutron stars as a function of their equatorial radius. The left panel presents results for fixed spin ratios, χ ≡ f / f K = 0.0 , 0.2, 0.4, 0.6, 0.8, and 1.0, while the right panel displays cases with fixed spin frequencies, f = 0 , 400, 600, 800, 1,000 Hz, and the Keplerian limit f = f K . For a given spin ratio or frequency, the gravitational mass decreases with increasing equatorial radius, but rises as the rotation rate increases, reflecting the additional centrifugal support provided by rotation. These results highlight the capability of RBHF theory to consistently describe both static and rapidly rotating neutron stars within the same microscopic framework. For the Bonn A potential, the static case yields a radius of 10.93 km at the maximum mass. In contrast, at the Keplerian frequency, the radius expands to 13.84 km, representing a 26.2% increase. In addition to the EoS obtained with the Bonn A potential, results based on the Bonn B and C potentials (Brockmann and Machleidt, 1990) are also presented. Overall, the mass–radius relations derived from these three parameterizations of realistic N N interactions exhibit very similar patterns, regardless of whether the stars are static or rotating. This robust consistency across different rotation rates suggests that the influence of rotational dynamics on the mass–radius relation is relatively insensitive to the specific details of the underlying interactions governing the EoS. Specifically, the maximum masses for χ = 0.0 , 0.2 , 0.4 , 0.6 , and 0.8 are 2.43, 2.46, 2.52, 2.63, and 2.80 M ⊙ , and the corresponding central energy densities are ε c = 2.26 , 2.25 , 2.21 , 2.12 , and 2.00 × 1 0 15 g / c m 3 , respectively. In particular, the maximum mass for rotating configurations M max , can reach up to 2.93 M ⊙ , which is 20.6% higher than the static result M TOV = 2.43 M ⊙ .

Specifically, the radii of a canonical neutron star with mass 1.4 M ⊙ in the static case are calculated to be R 1.4 = 11.98 , 12.17, and 12.32 km for the Bonn A, B, and C potential, respectively. The smallest radius predicted by the Bonn A potential implies that the RBHF calculations with this potential yields the softest EoS. This softness is attributed to the weakest tensor force in the Bonn A potential, which leads to the strongest attraction between nucleons. These differences underscore the sensitivity of neutron star properties to the underlying N N interactions and emphasize the importance of accurately modeling these interactions to predict astrophysical observables. Further discussions on the tensor force effects in nuclear matter, derived from realistic N N interactions, can be found in Ref. Wang et al. (2024). Moreover, the other results for neutron star radii shown in Ref. Tong et al. (2022) are also consistent with various empirical and observational constraints.

In the multimessenger era, the tidal deformability Λ of neutron stars has emerged as a crucial astrophysical constraint. The neutron star tidal deformabilties at 1.4 M ⊙ from the RBHF theory in full Dirac space are given as Λ 1.4 M ⊙ = 376,405,433 for Bonn A, B, C, respectively (Tong et al., 2022). Notably, the Bonn A potential predicts smaller Λ 1.4 M ⊙ values compared to Bonn B and C. This trend can be understood through the stiffness of the EoS: for a given neutron star mass, a softer symmetry energy results in more compact stellar configurations, leading to both smaller radii and reduced tidal deformabilities. These theoretical predictions can be contextualized with current observational constraints. The initial estimation for tidal deformability Λ 1.4 M ⊙ has an upper bound Λ 1.4 M ⊙ < 800 (Abbott et al., 2017) from the observation of binary neutron star (BNS) merger event GW170817. Subsequent revised analysis from LIGO and Virgo collaborations narrowed this to Λ 1.4 M ⊙ = 19 0 − 120 + 390 (Abbott et al., 2018). Importantly, the results from three potentials fall within these observational bounds, with the Bonn A results exhibiting the closest agreement with the central values inferred from GW170817.

Beyond tidal deformability, universal relations among neutron star observables offer an additional, largely EoS-independent avenue for cross-checking theoretical models against observations. Figure 2 examines the EoSs derived from RBHF theory in the full Dirac space, with the projection method, and the momentum-independence approximation with Bonn potentials–in light of the universal I -Love- Q relations (Wang et al., 2022c). The I -Love and Q -Love are illustrated in the top panels of Figure 2. The I - Q relations can also be found in Ref. Wang et al. (2022c). Along each curve, the mass or compactness serves as the single varying parameter, increasing towards the left in the plots. The universal relations are found to hold with high accuracy across different EoSs. Owing to their weak dependence on the specific EoS, a single empirical fit (black solid curves) can be applied, given by the functional form (Yagi and Yunes, 2017): ln y i = a i + b i ⁡ ln x i + c i ln x i 2 + d i ln x i 3 + e i ln x i 4 , (29) where the fitting coefficients are summarized in Table 4 of Wang et al. (2022c). These coefficients closely agree with those obtained in Ref. Yagi and Yunes (2017), based on a broad ensemble of EoSs. The bottom panels of Figure 2 present the absolute fractional deviations between the data and the fit, which remain below 1% over the entire mass range examined. The universal relation between I ̄ and the tidal deformability Λ enables the inference of the moment of inertia for a 1.4 M ⊙ neutron star, I ̄ 1.4 M ⊙ , from the tidal deformability Λ 1.4 M ⊙ measured in the GW170817 event. The updated analysis from the LIGO and Virgo Collaborations reports Λ 1.4 M ⊙ = 19 0 − 120 + 390 (Abbott et al., 2018), corresponding to I ̄ 1.4 M ⊙ = 10.3 0 − 2.10 + 3.39 , as shown in the left panel of Figure 2. Using the relation I ̄ = I / M 3 , this yields I 1.4 M ⊙ = 1.2 2 − 0.25 + 0.40 × 1 0 45 g c m 2 .

We have reviewed recent developments in RBHF theory within the full Dirac space, with particular emphasis on their implications for the properties of dense nuclear matter and neutron stars. This relativistic ab initio calculations enhance the internal consistency of relativistic many-body calculations and represent a significant advancement in the microscopic description of dense matter under extreme conditions. Further progress in the RBHF theory is anticipated through the inclusion of higher-order many-body correlations, in particular by extending beyond the two-hole-line expansion currently employed in standard RBHF theory. The incorporation of three-hole-line contributions and other higher-order terms is essential for achieving a more complete and quantitatively accurate description of in-medium nuclear interactions at supranuclear densities. In parallel, while a leading order and next-to-leading order covariant chiral nuclear forces have recently been applied within RBHF calculations under the momentum-independence approximation (Zou et al., 2024; Zou et al., 2025b; Zou et al., 2025a; Zheng et al., 2025; Shen et al., 2025), a natural next step is to implement the high-fidelity chiral nuclear forces (Ren et al., 2018; Lu et al., 2022) in the full Dirac space. Such an extension would enable a more consistent and comprehensive treatment of relativistic effects, thereby improving the predictive power of relativistic ab initio calculations for the EoS and neutron star properties.