Saturation of nuclear matter in the relativistic Brueckner Hatree-Fock approach with a leading order covariant chiral nuclear force

Nuclear saturation is a crucial feature in nuclear physics that plays a fundamental role in understanding various nuclear phenomena, ranging from properties of finite nuclei to those of neutron stars. However, a proper description of nuclear saturation is highly nontrivial in modern nonrelativistic~\textit{ab initio}~studies because of the elusive three-body forces. In this letter, we calculate the equation of state for nuclear matter in the relativistic Brueckner-Hartree-Fock (RBHF) framework with the leading order covariant chiral nuclear force. We show that a simultaneous description of the nucleon-nucleon scattering data and the saturation of the symmetric nuclear matter can be achieved. In this regard, the relativistic effects nicely explain the saturation of nuclear matter. As a result, the present study based on the covariant chiral nuclear force shows that in the RBHF framework, one can achieve saturation with a leading order covariant chiral nuclear force with only two-body forces, in contrast to the vast majorities of studies in the non-relativistic framework, where the next-to-next-to-leading order two-body and three-body chiral forces are needed. This study sets the foundation for studying nuclear saturation with the covariant chiral force in the RBHF framework, which allows for a systematic understanding of one of the key features of nuclear physics more microscopically.


I. INTRODUCTION
Nuclear matter is an ideal nuclear system of uniform density consisting of infinite neutrons and protons.The saturation of symmetric nuclear matter (SNM), i.e., the energy per nucleon reaches a minimum of about −16 MeV around a density of 0.16 fm −3 neglecting the Coulomb interaction, plays a vital role in our understanding of nuclear physics [1][2][3][4].The historical success of some of the early phenomenological models, such as the semi-empirical mass formula [5,6] and the liquid drop model [7], lies in the fact that a nucleus can be treated as an incompressible quantum liquid drop with a constant density ρ 0 at first approximation.This fact alone tells that any microscopic theory of atomic nuclei should be able to describe nuclear matter satisfactorily [8].
Given that the volume term is the most important one in the liquid drop model of nuclear binding energies, one may argue that saturation is a bulk property and theoretically should be independent of fine details of the nucleon-nucleon interaction.Nevertheless, a proper description of the empirical saturation properties of SNM remains a highly nontrivial task in modern non-relativistic ab initio approaches.That is, one needs to supplement the two-body NN forces (2NF) with the threebody NNN forces (3NF) to reasonably describe the saturation of SNM [9][10][11][12][13][14][15].In non-relativistic theories, the repulsion required for nuclear matter saturation is provided by many-body forces such as 3NF, which can be understood as the leadingorder (in powers of ρ B ) term of the velocity dependence inherent in the Lorentz scalar interaction [16].We note that the development of consistent microscopic 3NF for the two extreme systems, finite nuclei and nuclear matter, remains a matter of current research [17].More recently, many ab initio calculations of nuclear structure employing chiral two-and three-nucleon interactions (at N 2 LO and N 3 LO) found that the bulk properties of medium-mass nuclei and nuclear matter saturation cannot be consistently described [8,18,19].This prompted the development of interactions that use nuclear structure data of medium-mass nuclei or the saturation properties of SNM to constrain the nucleon-nucleon interaction.One such example is the N 2 LO SAT interaction [12] that simultaneously fitted the ground-state energies and radii of oxygen isotopes to obtain an NN+3N interaction at N 2 LO.Though very successful phenomenologically, such efforts seem to deviate from the conventional definition of ab initio studies (using the nuclear force calibrated with the nucleon-nucleon scattering and few-body data, and uniquely and correctly predicting the properties of nuclear matter and the structure of finite nuclei) [8,20,21].We note that how to define an ab initio theory is still a topic under debate; see, e.g., Ref. [22].
A different line of research is to treat atomic nuclei and nuclear matter as relativistic systems.In Ref. [23], with 2NF only, it was shown that the RBHF approach could yield saturation properties closer to the empirical value, in contrast with the non-relativistic BHF theory.Such studies hinted that nuclear matter saturation is a relativistic effect, namely, the Lorentz nature of the scalar and vector interactions leads to density-dependent repulsion, as pointed out in the previ-ous studies [24][25][26][27].In recent years, this line of research has gained much more momentum following the overcoming of several long-standing technical problems [28][29][30][31].However, the two-body bare NN interaction used in these studies is based on the phenomenological meson exchange picture, i.e., Bonn A, B, and C [32], which does not allow a systematic improvement for the description of SNM.From this perspective, a relativistic nuclear force more closely related to the theory of the strong interaction, i.e., Quantum Chromodynamics (QCD), is needed to gain further insight into the role of relativity in the saturation of nuclear matter.The present work aims to fill this gap.
Chiral effective field theory (χEFT), as the low-energy effective theory of QCD, can provide microscopic and QCDbased nucleon-nucleon interactions [33,34] needed in ab initio studies.In the past, most studies focused on nonrelativistic nucleon-nucleon interactions [35,36].Recently, it was proposed that one can derive a relativistic chiral nucleonnucleon interaction based on the covariant χEFT [37][38][39], which provides the much-needed microscopic NN forces for the RBHF studies [40].The leading order covariant chiral nuclear force can describe the nucleon-proton scattering phaseshits reasonably well [38,41], and its renormalizability has been studied, which indeed shows better performance compared to the Weinberg force.For instance, the much-discussed 3 P 0 channel is cutoff independent, different from its counterpart in the Weinberg formulation [38,42].
In this letter, we study the properties of SNM in the RBHF theory.We employ the relativistic chiral nuclear force to demonstrate how relativity allows one to achieve saturation already with the leading order (LO) two-body force, indicating that saturation can indeed be seen as originating from relativistic effects and is undoubtedly a bulk property less affected by the fine details of the nucleon-nucleon interaction.This letter is organized as follows.In Sec.II, we briefly introduce the theoretical framework of the RBHF theory based on the LO covariant chiral nuclear force.We analyze the results in Sec.III, followed by a summary and outlook in Sec.IV.

II. THEORETICAL FRAMEWORK
In the RBHF theory, one employs the Dirac equation to describe the single-particle motion of the nucleon in nuclear matter: where α and β are the Dirac matrices, u(p, λ) is the Dirac spinor with momentum p, single-particle energy E p and he-licity λ, M is the mass of the free nucleon, and U is the singleparticle potential operator, providing the primary medium effects.Due to time-reversal invariance, one can neglect the spacelike component of the vector fields, and as a result, the single-particle potential operator can be expressed as [23,43]: The momentum dependence of the scalar field (U S ) and the timelike component of the vector fields (U 0 ) is weak [31,44] and therefore neglected.The in-medium Dirac equation can be expressed in the form of the free Dirac equation: by introducing the following effective quantities: whose solution reads where χ λ is the Pauli spinor helicity basis.The covariant normalization is ū(p, λ)u(p, λ) = 1.
Once U S and U 0 of the single-particle potential operator are determined, the in-medium Dirac equation can be solved analytically.To achieve this, the matrix element of U is constructed following Refs.[23,45] as where the direction of p is taken along the z axis.Once Σ is obtained, U S and U 0 can be determined via the following relations, where p 1 and p 2 describe two momenta that are in the Fermi sea.
On the other hand, the matrix elements of U in Eq. ( 6) can be calculated as the integrals of the effective G matrix where Ḡ is the antisymmetrized G matrix.In the no-sea ap-proximation [27], the integral is only performed for the single-particle states in the Fermi sea.
In the nonrelativistic BHF theory, the Brueckner-Goldstone expansion has been shown to converge in the number of hole lines [46] (see also Ref. [47]).The lowest order of the holeline expansion, the two-hole-line, is the same as the first order G-matrix expansion, i.e., the RBHF framework, where the Gmatrix is obtained by solving the in-medium relativistic scattering equation.One of the most widely used scattering equations in the RBHF theory is the Thompson equation [48], a relativistic three-dimensional reduction of the Bethe-Salpeter equation [49].To include medium effects, the Thompson equation in the rest frame of nuclear matter reads, where is the relative momentum of the two interacting nucleons with momenta k 1 and k 2 , and q, q ′ , and k are the initial, final, and intermediate relative momenta of the two nucleons scattering in nuclear matter, respectively.W = E P +q + E P −q is used to describe the starting energy.M * and E * P ±k are effective masses and energies.The Pauli operator Q(k, P ) only allows the scattering of nucleons to un-occupied states, i.e., where k F is the Fermi momentum.
One key feature of the RBHF theory is that the V in Eq. ( 9) should be a bare NN interaction of covariant form.This work adopts the LO covariant chiral nuclear force [38,50].In this order, the relativistic potential V is the sum of the contact term and the one-pion-exchange diagram, where the contact potential (CTP) is where C S,V,AV,T are low-energy constants (LECs) to be determined by fitting to the nucleon-nucleon scattering phase shifts.The one-pion-exchange potential (OPEP) reads where the pion decay constant f π = 92.4MeV, the axial vector coupling g A = 1.29 [51].m π is the pion mass, q = (E p ′ − E p , p ′ − p) represents the four momentum transferred, and τ 1 , τ 2 are the isospin Pauli matrix.Furthermore, the potential has to be regularized to avoid ultraviolet divergences and facilitate numerical calculations.Here, we choose the commonly used separable cutoff function [41,52]: where n = 2 and Λ is the cutoff momentum.Eqs. ( 3),( 7), (8) , and (9) constitute a set of coupled equations that need to be solved self-consistently.Starting from chosen initial values of 0 , one solves the in-medium Dirac equation ( 3) to obtain the Dirac spinors.Next, one solves the Thompson equation ( 9) to obtain the G matrix and uses the integrals in Eq. ( 8) to get Σ.With Eq. ( 7), one obtains a new set of fields, U

S , U
(1) 0 , to be used in the next iteration.
When U S and U 0 of the single-particle potential converge, the binding energy per nucleon in SNM can be calculated as where the isospin indexes are suppressed.The starting energy is W = E p + E p ′ .The density ρ is related to the Fermi momentum k F through ρ = 2k 3 F /3π 2 .

III. RESULTS AND DISCUSSIONS
The unknown LECs are determined in Refs.[38,39] by performing a simultaneous fit to the J ≤ 1 Nijmegen partial wave phase shifts of the np channel up to the laboratory kinetic energy (E lab ) of 100 MeV at six energies [53].In the fitting process, the χ2 = i (δ i −δ i PWA93 ) 2 [38] is minimized.In Refs.[38,39], a sharp cutoff is used.In the present work, we adopt the exponential regulator in Eq. ( 14), which is more amenable for RBHF studies.In addition, to study the impact of the cutoff on the results, we vary Λ from 450 MeV to 600 MeV.The corresponding LECs C S,V,AV,T are listed in Table I  With the four sets of LECs, the descriptions of the Nijmegen multi-energy np phase shifts up to E lab = 300 MeV are shown in Fig. 8.The red bands cover the variations from the best-fit results obtained with the cutoff ranging from 450 MeV to 600 MeV.The LO covariant NN force can provide a reasonable description of the np phase shifts of 1 S 0 , 3 P 0 , 1 P 1 and 3 S 1 , but the description of the phase shifts of 3 P 1 , 3 D 1 , ε 1 is not very satisfactory at higher energies.Furthermore, the cutoff variation from 450 MeV to 600 MeV does not qualitatively change the overall picture.
In Fig. 10, we show the corresponding equations of state (EOS) for SNM.We note that the EOSs become more and more repulsive as the cutoff Λ decreases.At low densities, we observe a weak cutoff dependence.In contrast, at higher densities, the energy per particle of SNM strongly depends on the momentum cutoff Λ.We realized in retrospect that the cutoff dependence of the potential could be exacerbated in nuclear matter compared to that in vacuum, similar to the nonrelativistic case [13].The origin of such an exacerbation is not clear currently.In the future, we will investigate whether full cutoff independence can be achieved in the relativistic framework.We further note that for Λ of 567 MeV, the binding energy per nucleon obtained in the RBHF theory is −15.82MeV, in good agreement with the empirical value of −16±1 MeV, and the saturation density ρ 0 is 0.16 fm −3 which is also in good agreement with the empirical one of 0.16 ± 0.01 fm −3 .For Λ = 567 MeV, the compression modulus K ∞ at the saturation density is 252 MeV, again in agreement with the empirical value of 240±20 MeV [56].We expect that the dependence of the results on the cutoff will be alleviated when higher-order relativistic chiral forces are employed.On the other hand, compared with the results obtained with the Bonn A potential, at densities ρ > ρ 0 , the EOSs of SNM are much stiffer, although they agree with each other at lower densities.This might be traced to the fact that the LO covariant chiral nuclear force provides a relatively poor description of the phase shifts of P partial waves.
To better understand the EOSs of SNM, we show the contributions of different partial waves to the potential energy and those of the scalar and vector potentials to the potential energy in Fig. 5.Because the leading-order covariant chiral nuclear force is only fitted to partial waves with J ≤ 1, we only show the contribution of partial waves with J ≤ 1 to the nucleon potential energy.The contribution of higher partial waves with J ≥ 2 to the potential energy is small compared to that of the partial waves with J ≤ 1.However, the contributions of all the partial waves with J ≤ 7 are considered in obtaining the equation of state.As can be seen from Fig. 6, for both the chiral force and the Bonn A potential, the 1 S 0 and 3 S 1 -3 D 1 channels provide attraction and the P -waves provide repulsion.Although their magnitude and density dependence differ, the sums of the attraction and repulsion to the potential energy are similar below the saturation density.One can conclude that the saturation mechanism of the LO covariant chiral nuclear force is similar to but different from that of the Bonn A potential.This may explain why the EOSs of SNM obtained with the LO covariant chiral nuclear force are much stiffer at densities ρ > ρ 0 .
In the RBHF theory, it is known that the larger the scalar potential U S , the larger the relativistic effect [40].In other words, the small component of the Dirac spinor can become non-negligible even in the low-velocity limit, as M * can be significantly smaller than the bare mass M .In Fig. 7, we compare the full RBHF results and the one obtained by replacing the effective mass M * with the bare nucleon mass M in iteratively solving the set of coupled equations in the RBHF theory.The latter can be viewed as the non-relativistic result [23].The results show that the relativistic effect provides repulsion, which is consistent with the results obtained with the Bonn A potential.On the other hand, the density dependencies are different.However, due to the exploratory nature of the LO covariant chiral force employed, we do not discuss further the differences between the results obtained with the Bonn A potential and those obtained with the LO covariant chiral force.
Next, we compare in Fig. 5 the RBHF results with the BHF results obtained with the non-relativistic LO, NLO, and N 2 LO chiral nuclear forces [14].The solid line shows the EOS obtained with the LO covariant chiral nuclear force and a Λ = 567 MeV cutoff.There is no saturation in the BHF results up to NLO.In Refs.[13,14,38], it is shown that the NLO non-relativistic chiral nuclear force can provide a reasonable description of the np phase shifts of partial waves of J ≤ 1 similar to the LO covariant chiral nuclear force.On the other hand, the NLO non-relativistic chiral nuclear force cannot provide nuclear matter saturation in the BHF framework and all other non-relativistic frameworks.Only at N 2 LO, when 3NF is considered, a reasonable description of nuclear matter saturation can be achieved.On the other hand, the RBHF theory with the LO covariant chiral nuclear force is sufficient to provide a decent description of the saturation properties of SNM, which indicates that nuclear matter saturation can indeed be understood as a relativistic effect [23,25] and is independent of the fine details of the nucleon-nucleon interaction.Although the P wave contributions to the EOS

IV. SUMMARY AND OUTLOOK
In summary, we have studied the equation of state in the relativistic Brueckner-Hartree-Fock theory with the LO covariant chiral nuclear force.We found that the saturation properties of symmetric nuclear matter can be reasonably described with a cutoff of Λ = 567 MeV.This contrasts with the non-relativistic Brueckner-Hartree-Fock theory, where reasonable saturation can only be achieved at N 2 LO (with 3NF).On the other hand, in the RBHF theory, the saturation properties of symmetric nuclear matter can be reasonably described already at leading order, indicating that nuclear matter saturation can indeed be viewed as a relativistic effect, reinforcing the longheld belief [23,25,[28][29][30][31].In addition, we can achieve a simultaneous description of the np phase shifts and the saturation properties of symmetric nuclear matter at leading order without explicitly introducing three-body forces.In this respect, the results of the present work are consistent with the recent works seeking the essential elements of nuclear binding [57,58], but in a relativistic framework and highlights the unique role played by consistently treating relativistic degrees of freedom.We stress that this is the first RBHF study employing (covariant) chiral nuclear forces, narrowing the gap between RBHF studies and their non-relativistic counterparts.On the other hand, it should be noted that the dependence of the RBHF results on the momentum cutoff is much stronger than the bare relativistic chiral nuclear force.It is, therefore, desirable to check whether such a dependence becomes milder at higher chiral orders.Indeed, our preliminary results support such a trend.
In the future, we would like to perform studies with higherorder covariant chiral nuclear forces and check the convergence of our results.In addition, we would like to extend such studies to pure neutron matter, whose EOS plays an essential role in heavy-ion physics as well as in models of neutron stars, gravitational collapse supernovae, and neutron star mergers.

SUPPLEMENTARY MATERIAL FOR "SATURATION OF NUCLEAR MATTER IN THE RELATIVISTIC BRUECKNER HATREE-FOCK APPROACH WITH A LEADING ORDER COVARIANT CHIRAL NUCLEAR FORCE"
This supplemental material provides further details useful for understanding the results presented in the main text.

THE POTENTIALS OF SYMMETRIC NUCLEAR MATTER
The scalar and vector potentials are two important quantities in the RBHF model to connect the Dirac equation and G matrices through the nucleon single-particle potential, denoting the attraction and repulsion parts of the NN interaction at different ranges.We present in Fig. 6 the momentum cutoff dependence of the single-particle potentials U S and U 0 .Both attractive U S and repulsive U 0 , as large as several hundred MeV, are generated for all cutoffs.It can be seen from Fig. 6 that the sensitivity on the momentum cutoff mainly comes from the Lorentz vector potential.The results hinted that nuclear matter saturation is a relativistic effect, namely, the Lorentz nature of the scalar and vector interactions leads to density-dependent repulsion, as pointed out in the previous studies [24][25][26][27].

PARTIAL WAVE CONTRIBUTIONS TO THE POTENTIAL ENERGY
It is instructive to decompose the partial wave contributions to the potential energy and compare them with the results obtained with the Bonn A potential.As can be seen from Fig. 7, for both the chiral force and the Bonn A potential, the 1 S 0 and 3 S 1 -3 D 1 channels provide attraction, and the P -waves provide repulsive contributions.Although their magnitude and dependence on density differ, the sums of the attraction and repulsion to the potential energy are similar below the saturation density.One can conclude that the saturation mechanism for the leading-order (LO) covariant chiral nuclear force is similar but different from that for the Bonn A potential.One should not over-interpret the difference because the chiral force is only at leading order while the Bonn A potential is more accurate.
As seen from Fig. 7, the contribution of the partial waves J ≥ 2 to the potential energy is very small compared to the contribution of the partial waves J ≤ 1.
We show the contribution of 3 P 2 , 1 D 2 , and 3 D 2 in Fig. 8.We note that the contributions of 3 P 2 , 1 D 2 , and 3 D 2 depend weakly on the momentum cutoff.In addition, they are all smaller than those of the corresponding Bonn A potential.Given that the contributions of P -waves are more repulsive than those of the Bonn A, in the covariant chiral force we do not have the cancellations observed in the Bonn A potential, as shown in Fig. 9.This is unsurprising because the leading order covariant chiral force differs significantly from the Bonn A potential.Nonetheless, the sum of all the partial waves is similar to that of the Bonn A potential below the saturation point, as shown in Fig. 7.
We also show the Lorentz scalar and vector contributions to the potential energy in Fig.

1 FIG. 1 .
FIG. 1. (Color online) Neutron-proton phase shifts for partial waves of J ≤ 1.The red bands are those of the LO covariant chiral nuclear force, with the cutoff ranging from 450 MeV to 600 MeV.The solid dotted lines represent the np phase shift analyses of Nijmegen [53].The gray backgrounds denote the energy regions where the theoretical results are predictions.

FIG. 2 .
FIG. 2. (Color online) Energy per nucleon (E/A) in SNM as a function of the density ρ in the RBHF theory obtained with the leadingorder covariant chiral nuclear in comparison with the results obtained with the Bonn A potential (magenta dash line)[23].The pentagrams denote the saturation point.The shaded area indicates the empirical values[54,55].

FIG. 3 .
FIG. 3. (Color online) Partial wave and the scalar and vector potentials' contributions to the potential energy in SNM obtained with the LO chiral nuclear force, in comparison with the results obtained with the Bonn A potential (magenta dash line).(a): shows the contributions of the 1 S0 and 3 S1-3 D1 channels; (b): shows the contributions of P-waves; (c): shows the sum of the contributions of the partial waves with J ≤ 1; (d): shows the sum of the contributions of the partial waves with J ≥ 2; (e): shows the contributions of the scalar and vector potentials; (f): shows the energy per nucleon in SNM.

FIG. 4 .ΛFIG. 5 .
FIG. 4.Repulsive relativistic effect on the potential energy of SNM in the RBHF obtained with the LO covariant chiral force, in comparison with the results obtained with the Bonn A potential.

FIG. 6 .
FIG. 6. (Color online) Momentum cutoff dependence of the single-particle scalar US and vector U0 potentials.

FIG. 8 .FIG. 9 .
FIG. 8. Contributions of 3 P2, 1 D2, and 3 D2 to the potential energy in SNM with the LO chiral nuclear force, in comparison with the results obtained with the Bonn A potential.
. In the following, we employ the four LO covariant chiral nuclear forces listed in TableIto calculate phase shifts and properties of SNM.

TABLE I .
Values of the LO LECs (in units of GeV −2 ) with the cutoff ranging from 450 MeV to 600 MeV.