Nuclear Matter and Neutron Stars from Relativistic Brueckner-Hartree-Fock Theory

The momentum and isospin dependence of the single-particle potential for the in-medium nucleon are the key quantities in the Relativistic Brueckner-Hartree-Fock (RBHF) theory. It depends on how to extract the scalar and the vector components of the single-particle potential inside nuclear matter. In contrast to the RBHF calculations in the Dirac space with the positive-energy states (PESs) only, the single-particle potential can be determined in a unique way by the RBHF theory together with the negative-energy states (NESs), i.e., the RBHF theory in the full Dirac space. The saturation properties of symmetric and asymmetric nuclear matter in the full Dirac space are systematically investigated based on the realistic Bonn nucleon-nucleon potentials. In order to further specify the importance of the calculations in the full Dirac space, the neutron star properties are investigated. The direct URCA process in neutron star cooling will happen at density $\rho_{\rm{DURCA}}=0.43,~0.48,~0.52$ fm$^{-3}$ with the proton fractions $Y_{p,\rm{DURCA}}=0.13$. The radii of a $1.4M_\odot$ neutron star are predicated as $R_{1.4M_\odot}=11.97,~12.13,~12.27$ km, and their tidal deformabilities are $\Lambda_{1.4M_\odot}=376,~405,~433$ for potential Bonn A, B, C. Comparing with the results obtained in the Dirac space with PESs only, full-Dirac-space RBHF calculation predicts the softest symmetry energy which would be more favored by the gravitational waves (GW) detection from GW170817. Furthermore, the results from full-Dirac-space RBHF theory are consistent with the recent astronomical observations of massive neutron stars and simultaneous mass-radius measurement.


INTRODUCTION
Compact matter in the universe, such as neutron stars, provides an intriguing interplay between nuclear processes and astrophysical observables (Lattimer & Prakash 2004;Burgio et al. 2021). The central density in massive neutron stars reaches several times the empirical nuclear matter saturation density (ρ 0 ≈ 0.16 fm −3 ), which is far from those encountered in terrestrial experiments. Thus, astrophysical observations of neutron star properties are crucial in constraining the equation of state (EOS) for dense nuclear matter, since the latter constitutes the basic input quantity for the theoretical reconstruction of a neutron star.
During the past decade, considerable progress has been achieved with precise measurements for the massive neutron stars PSR J1614-2230 (1.928 ± 0.017M ⊙ ) (Demorest et al. 2010;Fonseca et al. 2016), PSR J0348+0432 (2.01 ± 0.04M ⊙ ) (Antoniadis et al. 2013), and PSR J0740+6620 (2.14 +0.10 −0.09 M ⊙ ) (Cromartie et al. 2020). Recently, for the first time the mass and radius of PSR J0030+0451 were simultaneously measured by the Neutron star Interior Composition Explorer (NICER) collaboration (Raaijmakers et al. 2019). From NICER's data, two independent analysis groups have reported a mass of 1.44 +0.15 −0.14 M ⊙ with a radius of 13.02 +1.24 −1.06 km (Miller et al. 2019) and a mass of 1.34 +0.15 −0.16 M ⊙ with a radius of 12.71 +1.14 −1.19 km ). In addition, the advanced LIGO and Virgo Collaborations detected a gravitational wave (GW) signal from a binary neutron star (BNS) merger, i.e., GW170817 (Abbott et al. 2017). This has greatly advanced multimessenger astronomy and stimulated an intense research activity in the field of nuclear physics. Besides the masses and radii of BNS, a significant signature carried by GW is the tidal deformability which represents the mass quadrupole moment response of a neutron star to the strong gravitational field induced by its com-sbwang@cqu.edu.cn panion (Damour et al. 1992;Hinderer 2008;Flanagan & Hinderer 2008;Damour & Nagar 2009;Postnikov et al. 2010). In particular, the limits on the tidal deformability inferred from the GW signal have been widely used to constrain the neuron star radius Annala et al. 2018;Most et al. 2018;Tews et al. 2018), the neutron skin thickness of 208 Pb , and the nuclear matter EOS (Malik et al. 2018;Zhang et al. 2018;Tong et al. 2020). Two years after GW170817, the detection of a compact binary coalescence involving a 22.2 − 24.3M ⊙ black hole and a compact object with a mass of 2.50 − 2.67M ⊙ was reported by the LIGO and Virgo Collaborations, i.e., GW190814 (Abbott et al. 2020). The secondary object of GW190814 still challenges the physics community nowadays, since it may be either the heaviest neutron star or the lightest black hole ever discovered (Fattoyev et al. 2020;Vattis et al. 2020;Safarzadeh & Loeb 2020;Yang et al. 2020;Tsokaros et al. 2020;Huang et al. 2020;Bombaci et al. 2021).
Theoretically, many attempts have also been made to obtain the EOS of neutron star matter via different nuclear many-body theories. They can roughly be divided into two categories: density functional theories (DFTs) with effective nucleon-nucleon (N N ) interactions and ab-initio methods with realistic ones.
The effective N N interactions in DFTs, either non-relativistic or relativistic, are constructed to reproduce the properties of finite nuclei and nuclear matter around saturation density, such as Skyrme (Skyrme 1956;Vautherin & Brink 1972), Gogny (Dechargé & Gogny 1980), or relativistic mean-field (RMF) models (Boguta & Bodmer 1977;Reinhard 1989;Ring 1996;Meng et al. 2006;Meng 2016). However, since the isovector channels of the density functionals are loosely constrained in the fitting procedures, their predictions for the nuclear matter properties such as symmetry energy at higher densities still have large uncertainties (Li et al. 2008;Jiang et al. 2015).
In the RBHF calculations, one of the most important procedure is the self-consistent determination of the singleparticle potential of the nucleons with the effective G matrix inside the nuclear medium. According to the limitations of symmetries (Serot & Walecka 1986), the single-particle potential operator U is generally divided into scalar and vector components. However, due to the complexity of the procedure, there exists different approaches to extract the components of single-particle potential in the RBHF theory, i.e., the momentum-independence approximation (Brockmann & Machleidt 1990), the projection method (Horowitz & Serot 1987;Nuppenau et al. 1989;Gross-Boelting et al. 1999) and the RBHF theory in the full Dirac space (Wang et al. 2021).
In the momentum-independence approximation (Brockmann & Machleidt 1990), the momentum dependence of the scalar potential U S and the timelike component of vector potential U 0 are neglected by regarding them as constants at a given baryon density. When investigating the isospin asymmetric nuclear matter (ANM), this method can not determine the correct behavior of the isospin dependence of the single-particle potential (Ulrych & Müther 1997;Schiller & Müther 2001). In the projection method (Horowitz & Serot 1987;Nuppenau et al. 1989;Gross-Boelting et al. 1999), the effective G matrix is projected to Lorentz invariant amplitudes which are used to determine the scalar and vector components of the single-particle potential. This method is more accurate than the momentum-independence approximation, since it keeps the spacelike component of the vector potential and the momentum dependence of the single-particle potential. However, this method is still restricted to the positive-energy states (PESs) without considering the negative-energy states (NESs). In paticular, contradictory results for the isospin dependence of the single-particle potential are found between the two approximations (Ulrych & Müther 1997).
Recently, to avoid the approximations used in the momentum-independence approximation and the projection method in the Dirac space with PESs only, a new self-consistent RBHF calculation in the full Dirac space for both SNM and ANM has been developed (Wang et al. 2021(Wang et al. , 2022. The Lorentz structure and the momentum dependence of the single-particle potential are determined uniquely by decomposing the matrix elements of U in the full Dirac space. In particular, the long-standing controversy of the isospin dependence of the single-particle potential has been clarified for the opposite tendency predicted by the momentum-independence approximation and projection method.
In this work, it is timely and necessary to study the nuclear matter properties above the saturation density and neutron star properties with the latest RBHF theory in the full Dirac space using the Bonn potentials. In order to compare with the results obtained in the Dirac space with PESs only, we also perform the calculations with the momentum-independence approximation and the projection method. This paper is arranged as follows. The theoretical framework of the RBHF theory in the full Dirac space, the nuclear matter properties, and the neutron star properties are reviewed in Section 2. In Section 3, the properties of nuclear matter and neutron star are presented and discussed. The summary is given in Section 4. In the RBHF theory, the essential point is to describe the single-particle motion in nuclear matter by using the Dirac equation where α and β are the Dirac matrices, M is the nucleon mass, p and E p,τ are the momentum and single-particle energy, s and τ denote the spin and isospin. According to the translational and rotational invariance, time-reversal invariance, hermiticity, and parity conservation, the single-particle potential U τ in the infinite, uniform nuclear matter can be expressed as (Serot & Walecka 1986) where U S,τ is the scalar part of single-particle potential, U 0,τ and U V,τ are the timelike and spacelike parts of the vector potential.p = p/|p| is the unit vector. It should be emphasized that the single-particle potential in Eq.
(2) is the function of the nucleon momentum p.
With the definitions of the following effective quantities: the solutions of Eq. (1) in the full Dirac space are where u τ and v τ are in-medium baryon spinors with positive and negative energy, χ s and χ τ are the spin and isospin wave functions. The baryon spinor can be obtained once the single-particle potentials are determined. The singleparticle potentials in Eq.
(2) can be obtained uniquely from three matrix elements Σ ++ τ , Σ −+ τ , and Σ −− τ of the operator U τ in the full Dirac space (Anastasio et al. 1981;Poschenrieder & Weigel 1988), This indicates that one needs the matrix elements not only for the positive-energy solutions given in Eq. (4), but also the elements coupling positive-with negative-energy solutions given in Eq. (4) and those for the negative-energy solutions, i.e., Σ ++ τ , Σ −+ τ , and Σ −− τ . These three matrix elements can be evaluated in the the mean-field approximation through the effective N N interactions G matrix in the full Dirac space where k τ ′ F specifies the Fermi momentum for neutrons and protons, which has a relation to the total baryonic density ρ = ρ n + ρ p .Ḡ is the antisymmetrized G matrix, where ± in the superscript denotes PESs or NESs.
To derive the effective interactive G matrix in nuclear matter, one of the most widely used equations in the RBHF theory is the in-medium covariant Thompson equation (Brockmann & Machleidt 1990), where τ τ ′ = nn, pp or np, W is the starting energy, V τ τ ′ denotes a realistic bare N N interaction and it is constructed in terms of effective Dirac spinors in the full Dirac space as explained in Eq. (4). Three parametrizations of Bonn potentials, namely, Bonn A, B, and C (Machleidt 1989) are chosen here. P = 1 2 (k 1 + k 2 ) is the center-of-mass 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 . q, k and q ′ are the initial, intermediate, and final relative momenta of the two nucleons scattering in nuclear matter, respectively. Q τ τ ′ (k, P ) is the Pauli operator, which prevents N N scattering into occupied states in the nuclear medium.
The relativistic G matrix is self-consistently calculated with the single-particle potentials in the standard RBHF iterative procedure through equations (1), (5), (6), and (7). Once the iteration has converged, the binding energy per nucleon in nuclear matter can be calculated by

Nuclear matter properties
The binding energy per nucleon of ANM can be generally expressed as a power series in the asymmetry parameter where E/A(ρ, 0) is the binding energy per nucleon in SNM and E sym (ρ) is the nuclear symmetry energy, The binding energy per nucleon in SNM can be expanded around the saturation density ρ 0 , where E/A(ρ 0 , 0) denotes the binding energy per nucleon at ρ 0 . The second derivative of E/A(ρ, 0) with respect to ρ is given by the incompressibility K(ρ), and K ∞ is the value at ρ 0 . Around the saturation density ρ 0 , the symmetry energy can be similarly expanded as where E sym (ρ 0 ) is the value of the symmetry energy at saturation density, L and K sym are the slope parameter and curvature parameter of the nuclear symmetry energy at ρ 0 : Furthermore, K asy = K sym − 6L has been widely used to describe the isospin dependence of the incompressibility of ANM in Refs. (Baran et al. 2005;Chen et al. 2005;Li et al. 2008;Sun et al. 2008;Centelles et al. 2009).

Neutron star properties
The neutron star is described as the β-stable nuclear matter system in this work, which consists of nucleons and leptons (mainly e − and µ − ). The chemical potentials of nucleons and leptons satisfy the equilibrium conditions where µ e , µ µ , µ p , and µ n are the chemical potentials of leptons, proton and neutron. In addition, charge neutrality is required as With these constraints, the energy density of the β-stable nuclear matter is then obtained as where the leptons are treated as non-interacting Fermi gas (Bombaci & Lombardo 1991;Maieron et al. 2004) and Y i = ρ i /ρ (i = e, µ, n, p) are the equilibrium particle fractions. The various chemical potential for each particle i is Through solving Eqs. (16), (17), (18), and (19) simultaneously, one can obtain the particle fractions Y i at a given density ρ. Then the energy density ε can be easily obtained from Eq. (18). The pressure P for the β-stable nuclear matter is defined as Once the EOS of β-stable nuclear matter 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 & Volkoff 1939;Tolman 1939), 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. For a given EOS, the TOV equations have the unique solution that depends on the conditions of β-stable matter at the center, such as the central density or the central pressure. In addition, to solve the TOV equations, the EOS must cover full regions of the neutron star from the crust to the core. In the present work, we mainly concentrate on discussing the core region within the RBHF theory. For the crust, the EOS from the Baym-Pethick-Sutherland (BPS) (Baym et al. 1971b) and the Baym-Bethe-Pethick (BBP) (Baym et al. 1971a) model is used. Besides the masses and radii, another important property of neutron star, the tidal deformability Λ is defined as which represents the mass quadrupole moment response of a neutron star to the strong gravitational field induced by its companion (Damour et al. 1992;Hinderer 2008;Flanagan & Hinderer 2008;Damour & Nagar 2009;Postnikov et al. 2010). C = M/R is the compactness parameter, M and R are the neutron star mass and radius, and k 2 is the second love number where y R = y(R) can be calculated by solving the following differential equation: with The differential equation (24) can be integrated together with the TOV equations with the boundary condition y(0) = 2.
3. RESULTS AND DISCUSSIONS

Nuclear matter properties
In ANM, one of the most important questions is the isospin and density dependence of the single-particle potentials (van Dalen et al. 2005), which is related to a wildly-discussed topic concerning the proton-neutron mass splitting. In Fig. 1, the isospin and density dependence of single-particle potentials U S,τ , U 0,τ and U V,τ at the Fermi momentum k τ F for densities ρ = 0.16 fm −3 , 0.32 fm −3 and 0.48 fm −3 are obtained from RBHF calculation in the full Dirac space. At the density ρ = 0.16 fm −3 , the values of U V,τ for both neutron and proton are close to zero. For the other two components, the RBHF calculation in the full Dirac space gives U S,n < U S,p and U 0,n > U 0,p . It is also found that with increasing neutron excess, U S,n and U V,n decrease, while U 0,n increases. The case for proton shows an opposite behavior. In the RMF theory (Ring 1996), the isospin splittings of the scalar and the timelike part of the vector potential can be calculated as (Ulrych & Müther 1997) The parameters g i and m i (i = δ, ρ) refer to the meson-nucleon coupling constants and masses of the mesons. ρ s,τ is the scalar density. For ANM, Eqs. (26) lead to U S,n < U S,p and U 0,n > U 0,p , which is consistent with the results by the RBHF theory in the full Dirac space. This consistence underlines the crucial importance for the δ-meson exchange in the RMF theory (Kubis & Kutschera 1997;Hofmann et al. 2001;Liu et al. 2002;Roca-Maza et al. 2011). Further, since neutron star radii and tidal deformabilities are sensitive to the nuclear matter properties at suprasaturation density Zhang, Nai-Bo & Li, Bao-An 2019;Tong et al. 2020). We also study the isospin dependence of single-particle potentials at supra-saturation densities ρ = 0.32 fm −3 and 0.48 fm −3 . The differences between neutron and proton becomes larger when increasing nuclear matter density. In comparison to the density ρ = 0.16 fm −3 , the amplitudes of U V,τ for both neutron and proton are much larger than zero. However, U V,τ is neglected in the momentum-independence approximation even at supra-saturation densities. It means that determining the correct behavior of the isospin dependence of the single-particle potentials at supra-saturation densities is important to describe nuclear matter and neutron star properties. From single-particle potentials one can calculate the binding energies per nucleon E/A for ANM. In Fig. 2, E/A are depicted as functions of the density ρ with the asymmetry parameter α ranging from 0 to 1. The binding energy and the saturation density become progressively smaller when α increases. The EOS of ANM exhibits a minimum which disappears for α 0.8. For comparison, Fig. 2 also displays the EOS for SNM and PNM obtained with the variational calculations (Akmal et al. 1998) using the Argonne v 18 interaction  with relativistic boost corrections to the two-nucleon interaction, as well as three-nucleon interactions modeled with the Urbana force (Pudliner et al. 1995). It can be seen that the results obtained with two ab initio methods are comparable.
In Table 1, we summarize our results for the properties of nuclear matter from RBHF theory. In the second row, we show our results calculated by the RBHF theory in the full Dirac space. The results obtained by the projection method (Gross-Boelting et al. 1999) with the ps representation for the subtracted T matrix (van Dalen et al. 2004) are shown in the third row. The results obtained by the momentum-independence approximation are also shown in the fourth row, the scalar potential and the time component of the vector potential are extracted directly from the singleparticle potential energy at two casually selected momenta, 0.5k τ F and k τ F . The one-boson-exchange potential (OBEP) is used for the realistic N N interactions, which is defined as a sum of one-particle-exchange amplitudes of six non-strange bosons with given spin-parity properties (Machleidt 1989). The coupling strengths and form-factor parameters are determined by fitting to the N N scattering data and deuteron properties. Three different parametrizations of OBEP were constructed for the application to relativistic nuclear structure physics, which were denoted by Bonn A, B, and C. The three potentials differ essentially in the πN N form-factor parameter. Since the nuclear tensor force is essentially provided by the pion, the main difference among the three potentials is the strength of the tensor force which can be reflected in their predictions for D-state probability of the deuteron (P D = 4.5%, 5.1%, 5.5% for Bonn A, B, C,  (Akmal et al. 1998) are also shown with black circles and green squares, respectively. respectively) (Machleidt et al. 1987). As is well known, the strength of the tensor force determines the location of the nuclear matter saturation point (Coester et al. 1970;Brockmann & Machleidt 1990). Compared to potential Bonn B and C, all three methods by using potential Bonn A with the weakest tensor force can reproduce the empirical value (Bethe 1971;Sprung 1972) for saturation density ρ 0 and binding energy E/A. With the potential Bonn A, the incompressibility of SNM at saturation density is 258 MeV for the RBHF calculation in the full Dirac space and 254 MeV in the projection method, which agree well with the commonly accepted value of 240 ± 20 MeV (Garg et al. 2007;Garg & Colò 2018). As for ANM, the symmetry energy E sym and the slope parameter L in the full Dirac space with potential Bonn A are 33.1 MeV and 65.2 MeV, which are smaller than the values suggested in the projection method and the momentum-independence approximation. This implies that the full Dirac space predicts the softest symmetry energy at higher densities. Ksym (left panel) and Esym (right panel) as functions of L calculated by RBHF theory using the potential Bonn A in the full Dirac space (red star), in comparison with results obtained by the projection method (blue down-pointing triangle), the momentum-independence approximation (gray up-pointing triangle) and various density functionals (circles and squares) (Sagawa et al. 2007;Zhao et al. 2010). The shaded regions denote a recent constraint on Esym, L and Ksym by combining astrophysical data with PREX-II and chiral effective field theory (Essick et al. 2021). The solid orange line is a linear fit to the results of density functionals.
In several previous works, slope parameter L and curvature parameter K sym of symmetry energy E sym , as well as their correlations, were shown to have significant effects on the neutron star properties, e.g., the composition of the core (Lattimer et al. 1991), crust-core transition density and pressure (e.g., Vidaña et al. 2009;Zhang et al. 2018;Li & Magno 2020). At the saturation density, E sym has been well constrained by recent terrestrial nuclear experiments and astrophysical observations. However, L and K sym still remain relatively large uncertainties which characterize the E sym at supra-saturation densities. In Fig. 3, K sym , L and E sym are calculated by a series of DFTs and the ab initio RBHF theory. The DFTs are adopted from nonrelativistic models (e.g., SLy4 and those starting with S) (Dutra et al. 2012), to relativistic models (e.g., PK1, NL3, DD-ME) (Long et al. 2006;Dutra et al. 2014). Ab initio calculations have been carried out with the RBHF theory using the potential Bonn A in the full Dirac space, together with results obtained by the projection method and the momentum-independence approximation. E sym and K sym obtained from DFTs exhibit a positive correlation with L, and the results of RBHF theory are all in good agreement with this correlation. We also perform a linear fitting to this correlation, with the correlation coefficients r = 0.93 and r = 0.92. The results of fitting are depicted with the solid orange line. Our results are close to K sym = −307.862 + 3.292L with r = 0.972 given in Ref. (Dong et al. 2012) which used several different DFTs from the present work. It should be mentioned that K sym predicted by RBHF theory are negative values, i.e., -67.3 MeV in full Dirac space, -70.1 MeV in the projection method and -65.7 MeV in the momentum-independence approximation. In addition, the shaded regions in Fig. 3 denote a recent constraint on E sym , L, and K sym by combining astrophysical data with PREX-II and chiral effective field theory (Essick et al. 2021). It can be seen that the results in the full Dirac space located inside the region labeled by χEFT Astro+PREX-II posterior. Our results are also consistent with other recent works, such as K sym = −120 +80 −100 MeV at 68% confidence level from a Bayesian analysis of neutron star properties (Xie & Li 2020) and K sym = −209 +270 −182 MeV at 95% confidence level from combining the neutron skin data and the pure neutron matter EOS from chiral effective field theory calculations (Newton & Crocombe 2021  On the basis of the Eqs. (16), (17), (18), and (19), the particle fractions Y i (i = n, p, e, µ) for β-stable matter in Fig. 4 are calculated by RBHF theory using the potential Bonn A in the full Dirac space, together with the projection method, and the momentum-independence approximation. The neutron fractions Y n decrease gradually with increasing ρ in three methods. In contrast, the values of Y p , Y e , and Y µ grow steadily. According to the energetic constraints (Wiringa et al. 1988), muons only emerge at higher densities, and the onset densities for muons from three methods are 0.144 fm −3 , 0.132 fm −3 and 0.142 fm −3 . The muons appear latest in the full Dirac space, since it has the softest symmetry energy and generates the smallest proton fraction. We also observe that the behavior of particle fractions in the full Dirac space is qualitatively similar to the one seen in other two methods.
The direct URCA (DURCA) process n → p + e − +ν e , p → n + e + + ν e , is one of the scenario of neutron star cooling (Pethick 1992;Prakash et al. 1992;Prakash 1994). This process is allowed when the momenta for nucleons and leptons can be conserved, which can predict the threshold density of the DURCA process. It has been pointed out in Ref. (Lattimer & Prakash 2004) that the β-stable matter with a proton-to-neutron ratio in excess of 1/8, or the proton fraction Y p 1/9 can be cooled by the DURCA process. As for the results from RBHF theory, the DURCA process happens after muon occurs for three methods. The threshold densities are calculated as ρ = 0.432 fm −3 and Y p = 0.127 in the full Dirac space, ρ = 0.372 fm −3 and Y p = 0.126 in the projection method, and ρ = 0.464 fm −3 and Y p = 0.127 in the momentum-independence approximation. These results from three methods are consistent with the astronomical observations, which does not allow the DURCA process to be too early.
Based on the EOS from RBHF theory, the mass-radius relations of neutron stars can be calculated from the TOV equation. In Fig. 5, the neutron star mass-radius relations obtained are shown from the RBHF theory in the full Dirac space with potential Bonn A, together with the projection method, and the momentum-independence approximation. The radii R 1.4M⊙ of a 1.4M ⊙ neutron star from three methods are 11.97 km, 12.38 km, and 12.35 km, respectively. The prediction from the full Dirac space leads to the smallest radius for R 1.4M⊙ . The relatively small neutron star radius suggested by full Dirac space implies that the symmetry energy at higher densities is C a u s a l i t y Figure 5. Neutron star mass-radius relations obtained by RBHF theory using the potential Bonn A in the full Dirac space (red solid line), together with the projection method (blue dashed line), and the momentum-independence approximation (gray dotted line). The shaded areas represent the constraints from astronomical observables for the massive neutron star (Demorest et al. 2010;Fonseca et al. 2016;Antoniadis et al. 2013;Cromartie et al. 2020). The inner and outer contours indicate the allowed area of mass and radius of neutron stars by NICER's analysis of PSR J0030+451 (Miller et al. 2019). The excluded causality region is also shown (Lattimer & Prakash 2007).
soft. It should be noted that soft EOSs, which produce small values of radii for fixed M are favored by the GW170817 . It can be found that our results for neutron star radii are also consistent with other works, such as R 1.4M⊙ 13.76 km , R 1.4M⊙ 13.6 km (Annala et al. 2018), 12.00 km R 1.4M⊙ 13.45 km (Most et al. 2018) and 9.0 km R 1.4M⊙ 13.6 km (Tews et al. 2018) from the tidal deformability (Abbott et al. 2017(Abbott et al. , 2019, R 1.4M⊙ = 10.9 +1.9 −1.5 km from the Bayesian approach (Kumar & Landry 2019) and the recent constraints by NICER (Miller et al. 2019) for the mass and radius of PSR J0030+0451, i.e., mass 1.44 +0.15 −0.14 M ⊙ with radius 13.02 +1.24 −1.06 km. The 68% and 95% contours of the joint probability density distribution of mass and radius from the NICER analysis are also shown.
The neutron star matter is only considered to be composed of nucleons and leptons here, the central densities of massive neutron stars in this work can reach 5 to 7 times the nuclear saturation density ρ 0 = 0.16 fm −3 of neutrons and protons found in laboratory nuclei, which is far away from the region ρ 0.57 fm −3 where RBHF theory in the full Dirac space is applicable (Wang et al. 2022). Therefore, for RBHF theory in the full Dirac space, we also follow the strategy in Refs. (Rhoades & Ruffini 1974;Gandolfi et al. 2012;Wang et al. 2022). The maximally stiff or causal EOS given by p(ε) = c 2 ε − ε c is used to replace the EOS for neutron star matter above a critical density ρ c , where p is the pressure, ε is the energy density, c is the speed of light, and ε c is a constant. This EOS predicts the most rapid increase of p with ε without violating causality limit and the upper bound on the maximum mass of the neutron star. In our calculations, ρ c is given as 0.57 fm −3 and ε c is determined by ensuring that the energy density is continuous between the low-and high-density EOS. The upper bound maximum masses of neutron stars are predicted as 2.43M ⊙ in the full Dirac space, which is still smaller than the compact object with a mass of 2.50-2.67M ⊙ by LIGO Scientific and Virgo collaborations in GW190814 (Abbott et al. 2020). Moreover, the maximum masses are predicted as 2.31M ⊙ and 2.18M ⊙ from the projection method and the momentum-independence approximation, which are consistent with the measurements of massive neutron stars, such as PSR J1614-2230 (Demorest et al. 2010;Fonseca et al. 2016), PSR J0348+0432 (Antoniadis et al. 2013), and PSR J0740+6620 (Cromartie et al. 2020).
In the multimessenger era, another important constraint of the neutron star is the tidal deformability Λ. In Fig. 6, the tidal deformabilities Λ as functions of neutron star mass M from RBHF theory with potential Bonn A in the full Dirac space, the projection method, and the momentum-independence approximation are shown. These tidal deformabilities decrease with increasing the neutron star mass. For a given neutron star mass, the method which has a softer symmetry energy yield smaller neutron star radii and correspondingly smaller tidal deformabilities. The tidal deformabilities of a 1.4M ⊙ neutron star from three methods are predicted as Λ 1.4M⊙ = 376, 473, and 459, respectively. It can be found that the tidal deformabilities from the full Dirac space are significantly lower than those from the projection method and the momentum-independence approximation. The initial estimation for tidal deformability Λ 1.4M⊙ has an upper bound Λ 1.4M⊙ < 800 (Abbott et al. 2017) from the observation of BNS merger event GW170817. Then a revised analysis from LIGO and Virgo collaborations gives Λ 1.4M⊙ = 190 +390 −120 (Abbott et al. 2018). It is important to underscore that the results from three methods are located in these regions, especially that from the full Dirac space, which is the closest one to the central values of the GW170817 shown as the circle symbol. In Table 2, our results for neutron star properties, such as threshold of DURCA process, the radii and tidal deformability for 1.4M ⊙ , and central densities, are obtained by the RBHF theory in the full Dirac space, the projection method and the momentum-independence approximation with the potential Bonn A, B, C. It was pointed out that the strength of tensor force plays an important role in the investigation of neutron star properties by the RBHF theory in the momentum-independence approximation (Wang et al. 2020). Therefore, it is interesting to investigate the effect of tensor force on the neutron star properties in the RBHF theory, especially in the full Dirac space. The threshold densities for DURCA process are listed in the third column. As for the full Dirac space, ρ DURCA from potential Bonn A, B, C are 0.43 fm −3 , 0.48 fm −3 and 0.52 fm −3 , with the same DURCA proton fractions Y p = 0.13. An inverse correlation can be found between the DURCA threshold density and the symmetry energy, since Bonn A generates the largest symmetry energy at a given baryon density with the weakest tensor force. The critical mass for DURCA process are also listed in the fifth column, M DURCA from potential Bonn A, B, C are 1.28M ⊙ , 1.55M ⊙ , and 1.78M ⊙ in the full Dirac space. Recently, the analysis of a neutron star in the transient system MXB 1659-29 has founded that the critical mass is around 1.6M ⊙ (Brown et al. 2018), and we found that the results from Bonn B is the closest one to the current astrophysical observations. The neutron star radii and tidal deformabilties at 1.4M ⊙ from the full Dirac space are shown in the sixth and the eighth column, which are given as R 1.4M⊙ = 11.97 km and Λ 1.4M⊙ = 376 for Bonn A, R 1.4M⊙ = 12.13 km and Λ 1.4M⊙ = 405 for Bonn B, and R 1.4M⊙ = 12.27 km and Λ 1.4M⊙ = 433 for Bonn C. We found that R 1.4M⊙ and Λ 1.4M⊙ predicted by Bonn A are smaller than those from Bonn B and C. This indicates that Bonn A with the lowest strength of tensor force has the softest EOS and Bonn C has the stiffest EOS. In the ninth column, the upper bound on the maximum masses of neutron stars from the full Dirac space are shown as 2.43M ⊙ , 2.43M ⊙ , 2.44M ⊙ with Bonn A, B, C. The difference among the three kinds of results is smaller than 0.01M ⊙ . This is reasonable since neutrons dominate the nucleonic component of the maximum masses of neutron stars while the main difference between the three Bonn potentials is in their tensor force strength in the (T = 0) 3 S 1 -3 D 1 state that does not contribute to the (T = 1) neutron-neutron state.

SUMMARY
The single-particle potential can be used to describe the quasi-particle properties of a nucleon inside a strongly interacting nuclear matter system. However, because the incompleteness of the Dirac space, different approximations have been introduced to extract the momentum and isospin dependence of the single-particle potential in the RBHF theory. In contrast to the RBHF calculations in the Dirac space with PESs only, the single-particle potential can be determined in a unique way from RBHF theory in the full Dirac space. The scalar component for neutron is found to be smaller than that for proton at and above nuclear saturation density, i.e., U S,n < U S,p , which leads to the neutronproton Dirac mass splitting M * D,n < M * D,p . As for the timelike component of the vector potential, the RBHF theory in the full Dirac space gives U 0,n > U 0,p . At higher densities, the amplitudes of U V,τ for both neutron and proton are much larger than zero, it means U V,τ at higher densities can not be neglected. The saturation properties of symmetric and asymmetric nuclear matter are systematically investigated based on the Bonn potentials. The symmetry energy E sym and its slope parameter L at the saturation density are 33.1 and 65.2 MeV for potential Bonn A, which are consistent with the empirical and experimental values. Especially, the slope parameter L and the curvature parameter K sym predicted in the full Dirac space agree well with the recent results from χEFT Astro+PREX-II posterior.
To further study the importance of the calculations in the full Dirac space at supra-saturation densities, the neutron star properties are studied with the latest RBHF theory in the full Dirac space. The threshold densities of the DURCA process are calculated as ρ = 0.43 fm −3 and Y p = 0.13 in the full Dirac space. These results are consistent with the astronomical observations, which does not allow the DURCA process to be too early. The radii of a 1.4M ⊙ neutron star are calculated as R 1.4M⊙ = 11.97, 12.13, 12.27 km, and their tidal deformabilities are Λ 1.4M⊙ = 376, 405, 433 for potential Bonn A, B, C. Comparing the results obtained with the projection method and the momentum-independence approximation, the RBHF calculations in the full Dirac space predict the softest EOS which would be more favored by GW170817. Especially, the tidal deformability from the full Dirac space is the closest one to the central values of the GW170817 from a revised analysis by LIGO and Virgo collaborations. Among the potential Bonn A, B, C, it is also found that RBHF calculations with the potential Bonn A has the softest EOS and potential Bonn C has the stiffest EOS since the former one has the lowest strength of tensor force. Furthermore, the results from RBHF theory are consistent with the recent astronomical observations of massive neutron stars and simultaneous mass-radius measurement.
S.W. thanks Qiang Zhao for helpful discussions. This work was partly supported by the National Key R&D Program of China under Contracts No. 2017YFE0116700, the National Natural Science Foundation of China (NSFC) under Grants 12147102, the Fundamental Research Funds for the Central Universities under Grants No. 2020CDJQY-Z003 and No. 2021CDJZYJH-003, the MOST-RIKEN Joint Project "Ab initio investigation in nuclear physics". Part of this work was achieved by using the supercomputer OCTOPUS at the Cybermedia Center, Osaka University under the support of Research Center for Nuclear Physics of Osaka University.