Fully self-consistent relativistic Brueckner-Hartree-Fock theory for finite nuclei

TL;DR

This paper develops a fully self-consistent relativistic ab initio method using the Bonn potential to accurately compute ground state properties of finite nuclei without free parameters or three-body forces.

Contribution

It introduces a self-consistent relativistic Brueckner-Hartree-Fock approach for finite nuclei starting from the Bonn potential, improving upon previous fixed-basis and non-relativistic methods.

Findings

Accurately calculates binding energies, radii, and spin-orbit splittings for selected nuclei.

Demonstrates convergence properties of the RBHF equations for finite nuclei.

Provides results consistent with previous RBHF and non-relativistic calculations.

Abstract

Starting from the relativistic form of the Bonn potential as a bare nucleon-nucleon interaction, the full Relativistic Brueckner-Hartree-Fock (RBHF) equations are solved for finite nuclei in a fully self-consistent basis. This provides a relativistic ab initio calculation of the ground state properties of finite nuclei without any free parameters and without three-body forces. The convergence properties for the solutions of these coupled equations are discussed in detail at the example of the nucleus $^{16}$O. The binding energies, radii, and spin-orbit splittings of the doubly magic nuclei $^{4}$He, $^{16}$O, and $^{40}$Ca are calculated and compared with the earlier RBHF calculated results in a fixed Dirac Woods-Saxon basis and other non-relativistic ab initio calculated results based on pure two-body forces.

Tables

Table 1. Table 1: Total energy, charge radius, and c.m. correction for 16 O calculated by RBHF for PBV and PAV.

	$PBV$	$PAV$
$E$ (MeV)	$-110.1$	$-101.4$
$r_c$ (fm)	$2.566$	$2.577$
$E_{cm}$ (MeV)	$-11.83$	$-11.12$

Table 2. Table 2: Parameters for DWS potentials in Eqs. ( 27 ) and ( 28 ).

$V_0$ (MeV)	$R$ (fm)	$a$ (fm)	$W_0$ (MeV)	$R^{ls}$ (fm)	$a^{ls}$ (fm)
$-60$ to $-80$	3.10697	0.615	725.9136	2.8827	0.648

Table 3. Table 3: Binding energy per nucleon, charge radius, and proton 1 ​ p 1 𝑝 1p spin-orbit splitting of 16 O in the RBHF calculations with different choices for U a ​ b subscript 𝑈 𝑎 𝑏 U_{ab} in Eq. ( 18 ).

	${\epsilon }^{\prime }={\epsilon }_{\pi 1p1/2}$	${\epsilon }^{\prime }={\epsilon }_{\nu 1s1/2}$	Gap	Exp.
$E/A$ (MeV)	$-7.10$	$-6.88$	$-5.41$	$-7.98$
$r_c$ (fm)	2.56	2.57	2.64	2.70
$\Delta E_{\pi 1p}^{ls}$ (MeV)	5.4	5.4	5.4	6.3

Table 4. Table 4: Total energy, charge radius, matter radius, and π ​ 1 ​ p 𝜋 1 𝑝 \pi 1p spin-orbit splitting of 16 O calculated by RBHF theory with the interaction Bonn A Machleidt ( 1989 ) , in comparison with experimental data. The corresponding results from RBHF in fixed DWS basis Shen et al. ( 2016 ) , DDRHF with PKO1 Long et al. ( 2006 ) and PKA1 Long et al. ( 2007 ) , non-relativistic BHF  Hu et al. ( 2017 ) with V low − k subscript 𝑉 low 𝑘 V_{{\rm low}-k} derived from Argonne v 18 subscript 𝑣 18 v_{18} , Coupled-Cluster (CC) method Hagen et al. ( 2009 ) and No Core Shell Model (NCSM) Roth et al. ( 2011 ) with N 3 LO  Entem and Machleidt ( 2003 ) , and Nuclear Lattice Effective Field Theory (NLEFT)  Lähde et al. ( 2014 ) with N 2 LO  Epelbaum et al. ( 2009b ) are also included.

	$E$ (MeV)	$r_c$ (fm)	$r_m$ (fm)	$\Delta E_{\pi 1p}^{ls}$ (MeV)
Exp. Wang et al. (2012); Angeli and Marinova (2013); Ozawa et al. (2001); Coraggio et al. (2003)	$-127.6$	$2.70$	$2.54\left(2\right)$	$6.3$
RBHF, Bonn A	$-113$.$5$	$2.56$	$2.42$	$5$.$4$
RBHF (DWS) Shen et al. (2016)	$-120$.$7$	$2.52$	$2.38$	$6$.$0$
DDRHF, PKO1 Long et al. (2006)	$-128.3$	$2.68$	$2.54$	$6.4$
DDRHF, PKA1 Long et al. (2007)	$-127.0$	$2.80$	$2.67$	$6.0$
BHF Hu et al. (2017), AV18	$-134.2$		$1.95$	$13.0$
CC Hagen et al. (2009), N3LO	$-120.9$		2.30
NCSM Roth et al. (2011), N3LO	$-119.7\left(6\right)$
NLEFT Lähde et al. (2014), N2LO	$-121.4\left(5\right)$

Table 5. Table 5: Total energy, charge radius, single-particle energies, and π ​ 1 ​ p 𝜋 1 𝑝 \pi 1p spin-orbit splitting of 16 O calculated by RBHF theory with the interactions Bonn A, B, and C. All energies are in MeV and radius in fm.

	Exp.	Bonn A	Bonn B	Bonn C
$E$	$-127.6$	$-113.5$	$-102.7$	$-95.0$
$r_c$	$2.70$	$2.56$	$2.61$	$2.65$
${\epsilon }_{\nu 1s_{1/2}}$	$-47$	$-48.1$	$-45.2$	$-43.1$
${\epsilon }_{\nu 1p_{3/2}}$	$-21.8$	$-26.4$	$-24.9$	$-23.7$
${\epsilon }_{\nu 1p_{1/2}}$	$-15.7$	$-21.0$	$-20.2$	$-19.7$
${\epsilon }_{\pi 1s_{1/2}}$	$-44\pm 7$	$-43.9$	$-41.1$	$-39.1$
${\epsilon }_{\pi 1p_{3/2}}$	$-18.5$	$-22.5$	$-21.1$	$-20.0$
${\epsilon }_{\pi 1p_{1/2}}$	$-12.1$	$-17.1$	$-16.5$	$-16.0$
$\Delta E_{\pi 1p}^{ls}$	$6.3$	$5.4$	$4.6$	$4.0$

Table 6. Table 6: Total energy, charge radius, and proton radius of 4 He calculated by RBHF theory using the interaction Bonn A  Machleidt ( 1989 ) with PBV and PAV, in comparison with experimental data  Wang et al. ( 2012 ); Angeli and Marinova ( 2013 ); Lu et al. ( 2013 ) . The corresponding results from DDRHF with PKO1 Long et al. ( 2006 ) and PKA1 Long et al. ( 2007 ) , solution of the Faddeev-Yakubovsky (FY) equation  Nogga et al. ( 2000 ) with CD-Bonn, FY  Binder et al. ( 2016 ) with N 4 LO, NCSM  Navrátil ( 2007 ) with N 3 LO, NLEFT  Lähde et al. ( 2014 ) with N 2 LO, and BHF  Hu et al. ( 2017 ) with V low − k subscript 𝑉 low 𝑘 V_{{\rm low}-k} derived from Argonne v 18 subscript 𝑣 18 v_{18} are also included.

	$E$ (MeV)	$r_c$ (fm)	$r_p$ (fm)
Exp. Wang et al. (2012); Angeli and Marinova (2013); Lu et al. (2013)	$-28.30$	$1.68$	$1.46$
RBHF (PBV), Bonn A	$-35.05$	$1.83$	$1.64$
RBHF (PAV), Bonn A	$-26.31$	$1.90$	$1.73$
DDRHF, PKO1 Long et al. (2006)	$-28.45$	$1.90$	$1.72$
DDRHF, PKA1 Long et al. (2007)	$-28.28$	$2.06$	$1.90$
FY Nogga et al. (2000), CD-Bonn	$-26.26$
FY Binder et al. (2016), N4LO	$-24.27\left(6\right)$		1.547(2)
NCSM Navrátil (2007), N3LO	$-25.39\left(1\right)$		$1.515\left(2\right)$
NLEFT Lähde et al. (2014), N2LO	$-25.60\left(6\right)$
BHF Hu et al. (2017), AV18	$-25.90$

Table 7. Table 7: Total energy, charge radius, matter radius, and proton spin-orbit splitting for the 1 ​ d 1 𝑑 1d shell of 40 Ca calculated by RBHF with the interaction Bonn A Machleidt ( 1989 ) , in comparison with experimental data  Wang et al. ( 2012 ); Angeli and Marinova ( 2013 ); Coraggio et al. ( 2003 ) . The corresponding results from DDRHF with PKO1  Long et al. ( 2006 ) and PKA1  Long et al. ( 2007 ) , BHF  Hu et al. ( 2017 ) , NCSM Roth and Navrátil ( 2007 ) , and CC Hagen et al. ( 2007 ) with V low − k subscript 𝑉 low 𝑘 V_{{\rm low}-k} derived from Argonne v 18 subscript 𝑣 18 v_{18}   Wiringa et al. ( 1995 ) , and CC Hagen et al. ( 2010 ) with N 3 LO  Entem and Machleidt ( 2003 ) are also included.

	$E$ (MeV)	$r_c$ (fm)	$r_m$ (fm)	$\Delta E_{\pi 1d}^{ls}$ (MeV)
Exp. Wang et al. (2012); Angeli and Marinova (2013); Coraggio et al. (2003)	$-342.1$	$3.48$		$6.6\pm 2.5$
RBHF , Bonn A	$-290.8$	$3.23$	$3.11$	$5.8$
DDRHF, PKO1 Long et al. (2006)	$-343.3$	$3.44$	$3.33$	$6.6$
DDRHF, PKA1 Long et al. (2007)	$-341.7$	$3.53$	$3.41$	$7.2$
BHF Hu et al. (2017), AV18	$-552.1$		$2.20$	$24.9$
NCSM Roth and Navrátil (2007), AV18	$-461.8$		$2.27$
CC Hagen et al. (2007), AV18	$-502.9$
CC Hagen et al. (2010), N3LO	$-345.2$

Equations

$${\mathcal{L}}_{NNpv}$$

$${\mathcal{L}}_{NNs}$$

$${\mathcal{L}}_{NNv}$$

$$\mathcal{H}=\sum _i\frac{\partial \mathcal{L}}{\partial ({\partial }^0{\varphi }_i)}{\partial }^0{\varphi }_i-\mathcal{L},$$

$$H=\int {\mathrm{d}}^{3}r\mathcal{H}(\mathbf{r}).$$

$$H=\int {\mathrm{d}}^{3}r\bar{\psi }\left(-i\mathbf{\gamma }\cdot \nabla +M\right)\psi +\frac{1}{2}\sum _{\alpha }\int {\mathrm{d}}^3{r}_1{\mathrm{d}}^3{r}_2\bar{\psi }({\mathbf{r}}_1){\Gamma }_{\alpha }^{(1)}\psi ({\mathbf{r}}_1)D_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2)\bar{\psi }({\mathbf{r}}_2){\Gamma }_{\alpha }^{(2)}\psi ({\mathbf{r}}_2),$$

$${\Gamma }_s=$$

$${\Gamma }_{pv}=$$

$${\Gamma }_v^{\mu }=$$

$$\frac{{\Lambda }_{\alpha }^2-m_{\alpha }^2}{{\Lambda }_{\alpha }^2+{\mathbf{q}}^2},$$

$$D_{\alpha }(x_1,x_2)=\pm \int \frac{d^{4}q}{(2\pi {)}^4}\frac{1}{m_{\alpha }^2-q^2}{e}^{-iq\cdot (x_1-x_2)},$$

$$D_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2)$$

$${\psi }^{†}(\mathbf{r})=\sum _k{\psi }_k^{†}(\mathbf{r})b_k^{†},\psi (\mathbf{r})=\sum _k{\psi }_k(\mathbf{r})b_k,$$

$$H=\sum _{k{k}^{\prime }}⟨k|T|k^{\prime }⟩b_k^{†}{b}_{k^{\prime }}+\frac{1}{2}\sum _{kl{k}^{\prime }{l}^{\prime }}⟨kl|V|k^{\prime }{l}^{\prime }⟩b_k^{†}{b}_l^{†}{b}_{l^{\prime }}{b}_{k^{\prime }},$$

$$⟨k|T|k^{\prime }⟩$$

$$⟨kl|V_{\alpha }|k^{\prime }{l}^{\prime }⟩$$

$$\times D_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2){\bar{\psi }}_l({\mathbf{r}}_2){\Gamma }_{\alpha }^{(2)}{\psi }_{l^{\prime }}({\mathbf{r}}_2).$$

$$⟨ab|\bar{G}(W)|a^{\prime }{b}^{\prime }⟩=⟨ab|\bar{V}|a^{\prime }{b}^{\prime }⟩+\frac{1}{2}\sum _{cd}⟨ab|\bar{V}|cd⟩\frac{Q(c,d)}{W-{\epsilon }_c-{\epsilon }_d}⟨cd|\bar{G}(W)|a^{\prime }{b}^{\prime }⟩,$$

$$Q(c,d)=\{\begin{array}{ll}1,& \mathrm{for}{\epsilon }_{\mathrm{c}}>{\epsilon }_{\mathrm{F}}\mathrm{and}{\epsilon }_{\mathrm{d}}>{\epsilon }_{\mathrm{F}},\\ 0,& \mathrm{otherwise}.\end{array}$$

$$(T+U)|a⟩=e_a|a⟩,$$

$$U_{ab}=\frac{1}{2}\sum _{c=1}^A⟨ac|{\bar{V}}_{\mathrm{eff}}|bc⟩.$$

$$U_{ab}=\frac{1}{2}\sum _{c=1}^A⟨ac|\bar{G}({\epsilon }_a+{\epsilon }_c)+\bar{G}({\epsilon }_b+{\epsilon }_c)|bc⟩,$$

$$U_{ab}=\sum _{c=1}^A⟨ac|\bar{G}({\epsilon }_a+{\epsilon }_c)|bc⟩,$$

$$U_{ab}=\frac{1}{2}\sum _{c=1}^A⟨ac|\bar{G}({\epsilon }_a^{\prime }+{\epsilon }_c)+\bar{G}({\epsilon }_b^{\prime }+{\epsilon }_c)|bc⟩,$$

$$|a⟩=\sum _k{D}_{ka}|k⟩.$$

$$\sum _l\left(T_{kl}+U_{kl}\right)D_{la}=e_a{D}_{ka}.$$

$$\sum _l\left(T_{kl}+U_{kl}\right)D_{l{k}^{\prime }}=e_{k^{\prime }}{D}_{k{k}^{\prime }},$$

$$T_{k^{\prime }{l}^{\prime }}$$

$${\bar{V}}_{k^{\prime }{l}^{\prime }{m}^{\prime }{n}^{\prime }}$$

$$E_{cm}=\frac{{⟨\mathbf{P}⟩}^2}{2AM}$$

$$|a\rangle=\frac{1}{r}\left(\begin{array}[]{c}F_{n_{a}\kappa_{a}}(r)\Omega_{j_{a}m_{a}}^{l_{a}}(\theta,\varphi)\\ iG_{n_{a}\kappa_{a}}(r)\Omega_{j_{a}m_{a}}^{\tilde{l}_{a}}(\theta,\varphi)\end{array}\right),$$

$$\left(\begin{array}[]{cc}M+\Sigma(r)&-\frac{d}{dr}+\frac{\kappa}{r}\\ \frac{d}{dr}+\frac{\kappa}{r}&-M+\Delta(r)\end{array}\right)\left(\begin{array}[]{c}F_{a}(r)\\ G_{a}(r)\\ \end{array}\right)=e_{a}\left(\begin{array}[]{c}F_{a}(r)\\ G_{a}(r)\\ \end{array}\right),$$

$$\Sigma (r)=$$

$$\Delta (r)=$$

Full text

I Introduction

To understand the nuclear system from the underlying interaction between nucleons has been one of the central problems in nuclear physics. Because of the strong repulsive core at short distance Jastrow (1951), the realistic nucleon-nucleon ($NN$) interaction is notoriously difficult to be solved in the usual many-body framework. Many methods have been proposed in the past to treat this singular behavior, such as Brueckner theory Brueckner et al. (1954), variational method Jastrow (1955), Lee-Suzuki method Suzuki and Lee (1980), the unitary correlation operator method Feldmeier et al. (1998), the low momentum effective $NN$ interaction $V_{low-k}$ Bogner et al. (2002), and the similarity renormalization group Bogner et al. (2007). Recently, with the great progress of the high-precision $NN$ interactions, such as Reid93 Stoks et al. (1994), AV18 Wiringa et al. (1995), CD Bonn Machleidt (2001), or chiral potentials Epelbaum et al. (2009a); Machleidt and Entem (2011), and with the rapid increase of computational power, more and more ab initio methods have been developed to study the nuclear many-body system. Celebrated examples include the quantum Monte Carlo method Carlson et al. (2015), the coupled-cluster method Hagen et al. (2014), the no core shell model Barrett et al. (2013), the self-consistent Green’s function method Dickhoff and Barbieri (2004), the lattice chiral effective field theory Lee (2009), the in-medium similarity renormalization group Hergert et al. (2016), the Monte Carlo shell model Liu et al. (2012), or the Brueckner-Hartree-Fock (BHF) theory Day (1967).

Among all ab initio methods, the Brueckner-Hartree-Fock theory is one of the most promising theories for an extension to heavy nuclei. Historically, the Brueckner theory was introduced to deal with the hard core of the nuclear force in nuclear many-body calculations Brueckner et al. (1954). The basic idea is to describe nuclear structure in a mean-field approximation, but replacing the bare nuclear force by an effective interaction in the medium. This effective interaction is the reaction $G$-matrix, which takes into account the two-nucleon short range correlations by summing up all the ladder diagrams in the nuclear medium. In this way, the saturation property of nuclear matter can be obtained qualitatively Brueckner et al. (1954). A formal derivation of Brueckner theory was provided by Goldstone Goldstone (1957), and substantial progress has been made by Bethe and co-workers Bethe and Goldstone (1957); Bethe et al. (1963). This was one of the breakthroughs in microscopic nuclear many-body theory, and many developments have been done along this direction. The readers are referred to the review papers Day (1967); Rajaraman and Bethe (1967); Baranger (1969) for the basic idea of Brueckner theory, the three hole-line expansion beyond BHF, and BHF for finite nuclei, respectively.

However, in the 1970s it was realized that all non-relativistic potentials failed to reproduce the saturation properties of infinite nuclear matter in detail. The saturation points obtained with various forces are distributed along the so-called Coester line Coester et al. (1970), which systematically deviates from the empirical value. It is the general opinion that this discrepancy was caused by the missing of the three-body forces Fujita and Miyazawa (1957); Brown and Green (1969), which have been used in all the modern non-relativistic investigations phenomenologically. In this way, at the cost of additional phenomenological parameters, one was able to reproduce the saturation properties of nuclear matter Zuo et al. (2002) as well as the ground states and a few excited states of light nuclei Pieper and Wiringa (2001).

On the other hand, nuclear structure has also been investigated in a relativistic framework. Johnson and Teller showed already in 1955 that proper nuclear saturation properties can be obtained provided that the potential depends on the velocity of the nucleons Johnson and Teller (1955). Later this theory was reformulated by Duerr in a relativistically invariant way Duerr (1956). In this way, the collapse of the nucleus occurring in the non-relativistic theories for high kinetic energies was avoided, and at the same time an extremely strong spin-orbit coupling was found Duerr (1956), in agreement with experimental data on the magic numbers Goeppert-Mayer and Jensen (1955). The Hamiltonian proposed by Duerr was then applied to finite nuclei Rozsnyai (1961). After the one-boson-exchange potential was established gradually Green and Sawada (1967), Miller and Green developed a Relativistic Hartree-Fock (RHF) approach based on attractive scalar and repulsive vector meson-exchange forces Miller and Green (1972). The relativistic approach became popular when Walecka established the $\sigma+\omega$ model and applied it successfully to highly condensed matter Walecka (1974). See Ref. Meng (2016) for a recent review.

Inspired by the success of early phenomenological relativistic investigations for nuclear structure, several groups proposed the relativistic versions of BHF theory and developed the Relativistic Brueckner-Hartree-Fock (RBHF) method. In the pioneering works by the Brooklyn group  Anastasio et al. (1980, 1983), the wave function of nucleon was chosen as a Dirac spinor in free space and the relativistic effects were taken into account in the first-order perturbation theory. Further developments were advanced by Horowitz and Serot Horowitz and Serot (1984, 1987), Brockmann and Machleidt Brockmann and Machleidt (1984, 1990), and ter Haar and Malfliet ter Haar and Malfliet (1986); Haar and Malfliet (1987). In these studies, the relativistic effects produced a strong density-dependent repulsion and therefore the saturation point was shifted remarkably close towards the empirical value. Using perturbation theory, it was found that relativistic effects lead to three-body forces through virtual nucleon-antinucleon excitations in the intermediate states (the so-called Z diagrams) Brown et al. (1987). Later, using the newly developed Bonn A/B/C potentials Machleidt (1989), the nuclear matter saturation points obtained within the RBHF theory were located on a new Coester line improving significantly the old one Brockmann and Machleidt (1990). The study of RBHF theory in nuclear matter has also been extended to the investigations of optical potential Nuppenau et al. (1990); Xu et al. (2016), asymmetric nuclear matter Huber et al. (1995); van Dalen et al. (2004), or neutron stars Huber et al. (1996); Katayama and Saito (2013). See Ref. Müther et al. (2017) for a recent review.

Even though the RBHF theory has achieved great success in the study of nuclear matter properties, the corresponding progress in finite nuclei was rather slow. Due to its enormous computational requirement, for a long time, RBHF theory for finite nuclei has been available only with certain approximations, such as the effective density approximation (EDA) Müther et al. (1988), or the local density approximation (LDA) Celenza et al. (1990); Marcos et al. (1991); Gmuca (1991); Brockmann and Toki (1992); Fritz et al. (1993). In the LDA approach, the density-dependence of the effective interaction, i.e. of the $G$-matrix, in nuclear matter is mapped onto a density-dependent relativistic Hartree or Hartree-Fock (DDRH or DDRHF) model, which is easy to be solved for finite nuclei. However, this mapping is far from unique and therefore this method suffers from large ambiguities as discussed in detail in Ref. Van Giai et al. (2010). As a consequence, different LDA approaches lead to rather different results. Only very recently, for the first time the RBHF equations were solved directly for finite nuclei Shen et al. (2016). The adopted Dirac Woods-Saxon (DWS) basis, which is obtained by solving the Dirac equation with a Woods-Saxon potential Zhou et al. (2003), guarantees the full relativistic structure of the Dirac spinors. Furthermore, the angle averaging Schiller et al. (1999); Suzuki et al. (2000) is avoided by solving the Bethe-Goldstone (BG) equation in the rest frame. Taking the nucleus 16O as an example, convergence has been achieved with an energy cut-off close to 1.1 GeV and good descriptions of binding energy and radius have been obtained without any adjustable parameters.

In this work, we will adopt the self-consistent RHF basis in solving the BG equation and go beyond the previous work Shen et al. (2016), where the BG equation was solved in a fixed DWS basis. All the cut-offs in the calculation will be checked and presented in detail. We will discuss the self-consistent single-particle potential and its uncertainties. The center-of-mass motion will be treated in both projection before and after variation (PBV and PAV) methods. Finally, as examples, the doubly magic nuclei 4He, 16O, and 40Ca will be calculated in the RBHF framework. The ground state properties will be studied and compared with other ab initio calculations in the literature.

In Sec. II, the RBHF framework will be given. All the numerical details will be discussed at length using 16O as an example in Sec. III. Results for 4He, 16O, and 40Ca will be presented in Sec. IV. A summary and perspectives for future investigations will be given in Sec. V.

II Theoretical Framework

In this Section, we will outline the theoretical framework of the relativistic Brueckner-Hartree-Fock theory for finite nuclei. In particular, the relevant formulas will be shown explicitly with spherical symmetry.

The RBHF theory leads to a system of coupled non-linear equations. The relativistic Hartree-Fock equations produce an optimized self-consistent single-particle potential $U$ with the single-particle wave functions $|a\rangle$, and the single-particle energies $\varepsilon_{a}$. Filling them up to the Fermi surface leads to a product state $|\Phi\rangle$ and the corresponding mean-field energy $E_{\rm RHF}$. These results depend on the effective interaction $V_{\rm eff}$ used in the mean-field equations. As compared to conventional mean-field theory, where the effective interaction $V_{\rm eff}$ is phenomenological, in the RBHF theory this effective interaction is replaced by the $G$-matrix, derived in an ab initio calculation from the bare nucleon-nucleon interaction using the Bethe-Goldstone equation. Its solution depends in a self-consistent way on the mean-field potential $U$ and its single-particle wave functions and energies. Therefore, RBHF theory presents a coupled system of equations which has to be solved by iteration. Its starting point is a relativistic form of the bare nucleon-nucleon interaction.

II.1 Relativistic Brueckner-Hartree-Fock theory for finite nuclei

II.1.1 Realistic one-Boson-exchange Lagrangian

We start with a relativistic one-boson-exchange $NN$ interaction which describes the $NN$ scattering data Machleidt (1989):

[TABLE]

where $\psi$ denotes the nucleon field. The bosons to be exchanged are characterized by the index $\alpha$ and include the pseudoscalar mesons ($\eta,\pi$) with pseudovector ($pv$) coupling, the scalar ($s$) mesons ($\sigma,\delta$), and the vector ($v$) mesons ($\omega,\rho$). For each pair, e.g., ($\eta,\pi$), the first (second) meson has isoscalar (isovector) character. For the isovector mesons, the field operator $\varphi_{\alpha}$ will be replaced by $\vec{\varphi}_{\alpha}\cdot\vec{\tau}$ with $\vec{\tau}$ being the usual Pauli matrices.

II.1.2 Hamiltonian

The Hamiltonian density is obtained using the Legendre transformation,

[TABLE]

where $\phi_{i}$ represent the nucleon field $\psi$, the meson fields $\varphi_{\alpha}$, and the photon field.

In the stationary case the Hamiltonian is found as a three-dimensional integral over the Hamiltonian density:

[TABLE]

Eliminating the meson fields one finds the following many-body Hamiltonian for the nucleons Brockmann (1978),

[TABLE]

where $\Gamma_{\alpha}^{(1)},\Gamma_{\alpha}^{(2)}$ are the interaction vertices for particles 1 and 2, with the coordinates $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, respectively:

[TABLE]

For the Bonn interaction Machleidt (1989), a form factor of monopole-type is attached to each vertex. In momentum space it has the form:

[TABLE]

where $\Lambda_{\alpha}$ is the cut-off parameter for meson $\alpha$ and $\mathbf{q}$ is the momentum transfer.

In Minkowski space, the meson propagators $D_{\alpha}(x_{1},x_{2})$ are the retarded solutions of the Klein-Gordon equations,

[TABLE]

where $q$ is the four-momentum transfer between the two particles. The sign $-$ holds for the scalar (and pseudoscalar) mesons and the sign $+$ holds for the vector fields. The dependence on the zero-component momentum transfer $q_{0}$ (energy) reflects the retardation of the interaction between the two particles. In the Bonn interaction of Ref. Machleidt (1989), this effect was deemed to be small and was ignored from the beginning. In this way, the meson propagators are just Yukawa functions:

[TABLE]

Note, however, that with the form factor in Eq. (6) the meson propagators are no longer simple Yukawa functions. In practice, the relevant matrix elements are calculated numerically in the momentum space, see Appendix A for details.

Now we expand the nucleon-field operators $\psi(\mathbf{r}),$ $\psi^{{\dagger}}(\mathbf{r})$ in terms of a complete orthonormal static relativistic basis $|k\rangle$:

[TABLE]

where $b_{k}^{\dagger}$ and $b_{k}$ form a complete set of creation and annihilation operators for nucleons in the state $|k\rangle$, which can be of positive energy or of negative energy. Here $\psi_{k}(\mathbf{r})$ is the corresponding Dirac spinor. The quantum number $k$ characterizing the state $|k\rangle$ contains also the isospin $\tau=n,\,p$ for neutrons and protons. We then have the Hamiltonian for nuclear system in the second quantized form as:

[TABLE]

where the matrix elements are given by

[TABLE]

The two-body interaction $V$ contains contributions from the different mesons $\alpha$.

The indices $k,l$ run over an arbitrary complete basis of Dirac spinors with positive and negative energies, as, for instance, over plane wave states $u(\mathbf{k},s)$ and $v(\mathbf{k},s)$ in the momentum space Itzykson and Zuber (1980) or over the eigensolutions of a Dirac equation with potentials of Woods-Saxon shapes discussed in Refs. Koepf and Ring (1991); Zhou et al. (2003).

II.1.3 Bethe-Goldstone equation

As is well known, the matrix elements of the bare nucleon-nucleon interaction $\langle ab|V|cd\rangle$ are very large and difficult to be used directly in nuclear many-body theory. Within the Brueckner theory, one takes into account the fact that nucleons in the nuclear medium do not feel the same interaction as that in free space. All the states below the Fermi surface are occupied and therefore the Pauli principle allows only scattering processes into intermediate states above the Fermi surface. The $T$-matrix, which describes scattering processes in free space is therefore, in the nuclear medium, replaced by the $G$-matrix Brueckner et al. (1954). It sums up all the ladder diagrams with two particles in intermediate states above the Fermi surface. It is deduced from the Bethe-Goldstone equation Bethe and Goldstone (1957),

[TABLE]

where $\langle ab|\bar{V}|a^{\prime}b^{\prime}\rangle=\langle ab|V|a^{\prime}b^{\prime}-b^{\prime}a^{\prime}\rangle$ is the antisymmetrized two-body matrix element, $W$ is the starting energy, and $\varepsilon_{c}$, $\varepsilon_{d}$ are the single-particle energies of the two particles in the intermediate states. The Pauli operator $Q(c,d)$ allows the scattering only to states $c$ and $d$ above the Fermi surface. It is defined as

[TABLE]

In this paper, we use the convention that indices $a,b,c,\dots$ run over the single-particle states being the solutions of RBHF equation, whereas the indices $k,l,m,\dots$ run over an arbitrary complete single-particle basis. Of course, the single-particle energies $\varepsilon_{c}$ and $\varepsilon_{d}$ depend on the solution of the corresponding Hartree-Fock equation and therefore we are left with a coupled system of equations, which has to be solved by iteration.

II.1.4 Relativistic Hartree-Fock equation

The relativistic Hartree-Fock equation reads

[TABLE]

where $e_{a}=\varepsilon_{a}+M$ is the single-particle energy with the rest mass of nucleon $M$, and $U$ is the self-consistent single-particle potential. In conventional relativistic Hartree-Fock theory Brockmann (1978); Long et al. (2006) based on an effective interaction $V_{\rm eff}$, this potential is defined as

[TABLE]

According to the no-sea approximation, the index $c$ runs over all the occupied states in the Fermi sea. In the relativistic Brueckner-Hartree-Fock framework, the interaction $V_{\rm eff}$ in the definition of $U$ is replaced by the $G$-matrix. Because of the fact that the effective interaction in the medium $G(W)$ depends on the starting energy, this definition is not as straightforward as in the conventional Hartree-Fock case, where the effective interaction does not depend on energy. The connection between the matrix element $U_{ab}$ and $G(W)$ was first discussed in Ref. Bethe et al. (1963) in nuclear matter. In the framework of perturbation theory, it was shown that, according to the Bethe-Brandow-Petschek (BBP) theorem Bethe et al. (1963), a specific choice of the starting energy in terms of the single-particle energies causes a large set of diagrams beyond the Hartree-Fock level to vanish. The extension to finite nuclei gives the following results Baranger (1969); Davies et al. (1969):

[TABLE]

if $|a\rangle$ and $|b\rangle$ are both hole (i.e. occupied) states, and

[TABLE]

if $|a\rangle$ is a hole state and $|b\rangle$ is a particle (i.e. unoccupied) state, and

[TABLE]

if $|a\rangle$ and $|b\rangle$ are both particle states.

In the above expressions, $\varepsilon$ labels the self-consistent single-particle energy, while $\varepsilon^{\prime}$ is somewhat uncertain Davies et al. (1969). This means the matrix elements of the self-consistent potential $U_{ab}$ with both states $|a\rangle$ and $|b\rangle$ above the Fermi level are not well defined in the Brueckner-Hartree-Fock theory. Different choices have been proposed in the literature Rajaraman and Bethe (1967); Davies et al. (1969). One extreme is to set the potential $U_{ab}=0$ in Eq. (18). This is known as the gap choice in nuclear matter. Another extreme is to set $\varepsilon_{a}^{\prime}=\varepsilon_{a}$, which is known as the continuous choice in nuclear matter. The BHF theory can be viewed as the first-order approximation of the so-called hole-line expansion Day (1967), which orders the Bethe-Brueckner-Goldstone expansion diagrams according to the number of independent hole lines. In principle, the final result will not depend on different choices of the single-particle potential $U$ if this hole-line expansion is taken into account up to high orders. This has already been confirmed in a non-relativistic calculation for nuclear matter up to the three-hole-line level: as shown in Ref. Song et al. (1998) the resulting equation of state does not depend much on the choice of $U$. Meanwhile, it was found in Ref. Song et al. (1998) on the BHF level, i.e., on the two-hole-line level, that near the saturation density the continuous choice produces several MeV$/A$ more binding than the gap choice. The results in the three-hole-line level agree with each other for these two different choices, and lie in between the BHF results.

Following the discussions in Ref. Davies et al. (1969), we choose in the present investigation a prescription in between the above two extremes. Precisely, $\varepsilon^{\prime}_{a}=\varepsilon^{\prime}_{b}=\varepsilon^{\prime}$ is fixed as an energy among the occupied states and we discuss the difference of the results by fixing $\varepsilon^{\prime}$ as the highest and as the lowest energy of the occupied states in the Fermi sea.

II.1.5 Solution in the RHF basis

The coupled RBHF equations in finite nuclei are solved within a complete Dirac basis $\{|k\rangle\}$, i.e. in a basis in Dirac space with states of positive and negative energies. This means that the RBHF single-particle states $|a\rangle$ are expressed as linear combinations

[TABLE]

Because of the requirement of the completeness of the basis, $|k\rangle$ runs not only over the positive-energy states in and above the Fermi sea but also over the negative-energy states in the Dirac sea Zhou et al. (2003), even though the no-sea approximation is adopted for the RBHF or RHF calculations, which means that the sums in the evaluation of various densities run only over the occupied states in the Fermi sea. As a consequence, also the index $c$ in Eqs. (16–18) is restricted to the $A$ occupied states in the Fermi sea. The RHF equation in the basis $\{|k\rangle\}$ reads

[TABLE]

It is solved by diagonalization in the calculations of the present work.

It should be noticed that both the BG equation (12) and the single-particle potential (16–18) are defined in the RBHF single-particle basis $\{|a\rangle\}$. Therefore, it requires a double-iteration procedure.

We start with a set of trial single-particle states, a discrete DWS basis Zhou et al. (2003) with the single-particle wave functions $|k\rangle$ and the corresponding energies $\varepsilon_{k}$, and we evaluate the matrix elements of the kinetic energy $T_{kl}$ and the antisymmetrized two-body matrix elements of the bare interaction $\bar{V}_{klmn}$ in this basis. This is a rather lengthy process, but it is done only once. Details are given in Appendix A.

Then we start the iteration:

We solve the BG equation (12) in this basis by matrix inversion. This yields a first set of $G$-matrix elements in the DWS basis. 2. 2.

These $G$-matrix elements are used to evaluate the single-particle potentials $U_{kl}$ in Eqs. (16–18). 3. 3.

The RHF equation (14) is solved as

[TABLE]

and a new set of single-particle states $\{|k^{\prime}\rangle\}$ with single-particle energies $\varepsilon_{k^{\prime}}$ is obtained. The iteration has converged if $\{|k^{\prime}\rangle\}=\{|k\rangle\}\equiv\{|a\rangle\}$. Otherwise the set $\{|k^{\prime}\rangle\}$ forms a new RHF basis and we continue. 4. 4.

The single-particle matrix elements of the kinetic energy and the two-body matrix elements of the bare interaction are transformed to the new RHF basis

[TABLE]

Then we go back to Step 1 and solve the BG equation (12) in this new RHF basis. This yields a second set of $G$-matrix elements and we continue with Step 2.

In practice, we found that performing an RHF iteration in Step 3 in each step of the RBHF iteration can speed up the convergence, and avoid the very time consuming Step 1, i.e. the solution of the BG equation.

This complicated iteration scheme allows the fully self-consistent solution of the RBHF problem. In Ref. Shen et al. (2016) a simplified version has been applied. In that case Step 4 was omitted and for each step of the iteration the BG equation was solved in the original DWS basis. Only the changes in the single-particle energies $\varepsilon_{k^{\prime}}$ in the propagator were taken into account. In practice this means, that the self-consistent single-particle wave functions in the Pauli operator in Eq. (13) were replaced by the corresponding one DWS basis.

II.1.6 Center-of-mass motion

In the above formulation, the spurious center-of-mass (c.m.) motion is included in the total Hamiltonian in Eq. (9). It is not of interest and should be removed. Since the Hamiltonian is invariant against translations, the exact many-body eigenstates of the system should be eigenfunctions of the total momentum $\mathbf{P}=\sum_{i}^{A}\mathbf{p}_{i}$. The spurious center-of-mass motion is therefore removed by projection onto the eigenspaces with vanishing eigenvalues of this operator. It has been shown in Refs. Beck et al. (1970); Ring and Schuck (1980) that, for large values of the particle number $A$, the projected energy is obtained in a good approximation by removing the center of mass energy

[TABLE]

from the total energy $\langle H\rangle=\langle T+V\rangle$.

In most of the (R)HF calculations the variation is carried out without projection, i.e., the (R)HF equations are solved for the total Hamiltonian $H$ and the spurious center of mass energy in Eq. (24) is removed after the variation, as for instance in Ref. Long et al. (2004). This is a projection after variation (PAV). A more strict treatment Zeh (1965) would be to exclude this term also in the (R)HF equation, i.e., to carry out a projection before variation (PBV), as it has been done for instance in Ref. Chabanat et al. (1998). We will discuss these two different choices in the RBHF framework as it has been done in Ref. Becker et al. (1974) for the non-relativistic BHF framework. We also should mention that the two-body part of the center of mass correction $\langle{\bf{P}}\rangle^{2}/2AM$ as well as the Coulomb force (see Appendix C) has only been taken into account in the RHF equation (21) and not in the BG equation for the calculation of the $G$-matrix.

II.2 RBHF theory for spherical nuclei

II.2.1 Spherical DWS basis

The eigenfunctions of a Dirac equation with spherical symmetry can be written as (for simplicity we neglect here the isospin indices):

[TABLE]

where $\Omega_{jm}^{l}(\theta,\varphi)$ are the spinor spherical harmonics. The radial, orbital angular momentum, total angular momentum, and magnetic quantum numbers are denoted by $n,\,l,\,j,$ and $m$, respectively, while the quantum number $\kappa$ is defined as $\kappa=\pm(j+1/2)$ for $j=l\mp 1/2$. $\tilde{l}=2j-l$ is the orbital angular momentum for the lower component. $F(r),\,G(r)$ are the radial wave functions which satisfy the radial Dirac equation:

[TABLE]

where $\Sigma=V+S$ and $\Delta=V-S$ are the sum and the difference of vector and scalar potentials. For a DWS basis Zhou et al. (2003), $\Sigma(r)$ and $\Delta(r)$ are potentials with Woods-Saxon shape parameterized as

[TABLE]

II.2.2 Symmetries of the BG equation

The BG equation has to be solved in the space of particle pairs. The quantum numbers of each particle are labeled as $a=(n_{a},j_{a},\pi_{a},t_{a})$ with the parity $\pi=(-)^{l}$ and the isospin quantum number $t=n,\,p$. Symmetries including rotation, parity, and charge can be used to reduce the dimension of this space. We therefore consider particle pairs with angular momentum ${\bf{J}}={\bf{j}}_{a}+{\bf{j}}_{b}$, parity $P=\pi_{a}\pi_{b}$, and $z$-component of the isospin (charge) $T=t_{a}+t_{b}$:

[TABLE]

Introducing the particle-particle ($pp$) $jj$-coupled matrix elements in Eq. (43), the BG equation can be reduced into sub-systems (channels) with the channel quantum numbers $JPT$:

[TABLE]

where the $m$ degeneracy has been summed up by Eq. (29). The starting energy $W$ is chosen according to Eqs. (16–18).

In practice, the BG equation has to be solved independently for each channel with the quantum numbers $JPT$. For the isospin, there are three channels ($T=-1,0,+1$), which can be solved independently as

[TABLE]

The RHF equation is a single-particle equation of the structure $b^{\dagger}b$. Therefore we need the antisymmetrized particle-hole (ph) $jj$-coupled matrix elements as given in Eqs. (44) and (45). In spherical nuclei, only $I=0$ matrix elements are relevant Tarbutton and Davies (1968):

[TABLE]

where $\hat{j}=\sqrt{2j+1}$. The antisymmetrized $ph$-coupled matrix elements $\langle ac|\bar{G}(W)|bc\rangle_{ph}^{I=0}$ are derived from the $pp$-coupled matrix elements, i.e. from the solutions of the BG equation (30) using Eq. (52):

[TABLE]

The maximum $J$ in $pp$ coupling used in the BG equation is determined by the single-particle angular momentum cut-off $j_{\mathrm{cut}}$ of the basis. In practice, the high $J$ matrix elements give relatively small contributions to $\langle 12|\bar{G}|34\rangle_{ph}^{I=0}$, thus we introduce a cut-off $J_{\mathrm{cut}}$ and solve the $G$-matrix for $0\leq J\leq J_{\mathrm{cut}}$.

In practice, we first construct the $ph$ coupled matrix elements $\bar{V}_{ph}^{I}$ (for details see Appendix A). From those we obtain the $pp$ coupled matrix elements $\bar{V}_{pp}^{J}$ by an inverse recoupling in Eq. (51). They are used to solve the BG equation (30) and to obtain $\bar{G}_{pp}^{J}$ for all $J$ values satisfying $0\leq J\leq J_{\mathrm{cut}}$. Finally using Eq. (II.2.2) and recoupling again to the $ph$ channel, the RHF equation can be solved for the next step.

II.2.3 Observables

After getting the solution of the coupled system of RBHF equations, the total energy is expressed as

[TABLE]

with the starting energy $W=\varepsilon_{a}+\varepsilon_{b}$. In the present framework, the Coulomb energy $E_{C}$ is not included in the evaluation of the $G$-matrix and is calculated separately. $E_{cm}$ is the center of mass energy in Eq. (24). For details in the calculation of $E_{cm}$ and $E_{C}$ see Appendices B and C.

The charge density distribution is obtained from Campi and Sprung (1972); Long et al. (2004):

[TABLE]

where $\lambda^{2}=1/(a^{2}-B^{2})$, with $a=\sqrt{\frac{2}{3}}\times 0.8$ fm the correction due to the finite proton size, and $B^{2}=\frac{3}{2}\langle\mathbf{P}^{2}\rangle^{-2}$ the correction due to the center-of-mass motion. For PBV, $B^{2}=0$ since the center-of-mass correction has already been considered in the single-particle wave functions. The charge radius is calculated as

[TABLE]

III Numerical details

As an example we consider the nucleus 16O. We use the realistic $NN$ interaction Bonn A which has been adjusted to the $NN$ scattering data Machleidt (1989). The Woods-Saxon potential in Eqs. (27,28) is taken from Ref. Koepf and Ring (1991). The DWS basis is then obtained by solving the spherical Dirac equation (26) in a box with the box size $R$ and mesh size $dr=0.05$ fm.

The RHF equation (14) can be solved in either the DWS basis or the obtained RHF basis. The validity of this RHF code is confirmed by reproducing the results of other RHF codes in the oscillator basis Lalazissis et al. (2009) and in coordinate space Long et al. (2006).

The BG equation (30) is solved by matrix inversion Haftel and Tabakin (1970) in the space of pair states $|ab\rangle$ given in Eq. (29). These are pairs of Dirac spinors with the full relativistic structure coupled to good angular momentum $J$ ($pp$ coupling). The indices $a$ and $b$ run over all solutions of the Dirac equation (with positive and negative energies). The BG equations are solved for each of the channels characterized by the quantum numbers $J$, parity $\pi$, and $z$-component of the isospin $T=t_{a}+t_{b}$. This leads, for various $J$ values, to a set of $pp$-coupled matrix elements of the $G$-matrix. Because we work always in the RHF basis during the iteration, the Pauli operator in Eq. (13) and all its relativistic structure is here fully taken into account. In particular, there is no angle averaging involved as it is the case in most Brueckner calculations in nuclear matter Schiller et al. (1999); Suzuki et al. (2000). The fully self-consistent RHF basis is used in comparison with previous investigation with a fixed DWS basis Shen et al. (2016). The BG equation (30) is solved for four different values of the starting energy $W$, equally distributed between the lowest and the highest single-particle energies in the Fermi sea. The $G$-matrix with specific starting energy $W$ is obtained by a four point Lagrange polynomial interpolation Davies et al. (1969).

In the following convergence check, $\varepsilon^{\prime}$ in Eq. (18) is chosen as the minimum single-particle energy of the hole states ($\varepsilon^{\prime}=\varepsilon_{\nu 1s1/2}$). The center-of-mass correction is not considered. After the convergence check, discussion will be focused on the single-particle potential of particle states in Eq. (18) as well as PAV and PBV. Without explicitly stating, all the calculations are performed in the self-consistent RHF basis.

III.1 Convergence check

It is well known that the bare $NN$ interaction contains a repulsive core and a strong tensor part connecting the nucleons below the Fermi surface to the states with high momentum in the continuum. In order to take this coupling fully into account, one needs a relatively large basis space and it is crucial to investigate the corresponding convergence of the RBHF calculations. In our previous investigation with a fixed basis Shen et al. (2016), a reasonable convergence is achieved near an energy cut-off $\varepsilon_{\rm cut}=1.1$ GeV. In this work, we will carry out the same check for the fully self-consistent RBHF calculation. Moreover, we will give a convergence check for the other cut-offs, i.e., the single-particle orbital angular momentum cut-off $l_{\rm cut}$, the single-particle energy cut-off in the Dirac sea $\varepsilon_{\rm Dcut}$, and the total angular momentum cut-off $J_{\rm cut}$ in the derivation of the $ph$ matrix elements of the G-matrix in Eq. (II.2.2).

Except these single-particle cut-offs for the basis space, we have introduced pair cut-offs with $\varepsilon_{1}+\varepsilon_{2}\leq E_{\rm cut}$ and $l_{1}+l_{2}\leq L_{\rm cut}$. The difference between the single-particle cut-off and the pair cut-off is illustrated in Fig. 1. Since the matrix elements with small total angular momentum $J$ in the $pp$ coupling have a larger contribution to the matrix elements of total angular momentum $I=0$ in the $ph$ coupling, the pair cut-offs will be introduced for $J>J_{h}$ with $J_{h}$ defined as the largest total angular momentum that the hole states can couple to. For 16O, we have $J_{h}=3~{}\hbar$ coupled from two particles in the $1p_{3/2}$ orbit. In present work, we set $E_{\rm cut}=\varepsilon_{\rm cut},\,L_{\rm cut}=l_{\rm cut}$, and the convergence in the pair cut-off will be achieved automatically in the convergence of the single-particle cut-off.

First, we give the convergence check for single-particle angular momentum cut-off $l_{\rm cut}$ in Fig. 2. In this check, the other cut-offs are chosen as $J_{\rm cut}=6~{}\hbar,\,\varepsilon_{\rm cut}=1100$ MeV, and $\varepsilon_{\rm Dcut}=-1700$ MeV. With $\varepsilon_{\rm Dcut}=-1700$ MeV, single-particle states with high angular momentum ($l>7~{}\hbar$) in the Dirac sea are not included in the basis, this will be discussed later in the check of $\varepsilon_{\rm Dcut}$. At $l_{\rm cut}=20~{}\hbar$, the convergence is achieved. Increasing $l_{\rm cut}$ to $25~{}\hbar$ will change the energy by $0.6$ MeV and charge radius by $0.003$ fm.

As shown in the RBHF calculation with fixed DWS basis Shen et al. (2016), one needs a very large energy cut-off to take into account the short-range correlation of the strong repulsive core. The general feature of this convergence still holds in the self-consistent RBHF calculation as shown in Fig. 3. The other cut-offs are $J_{\rm cut}=6~{}\hbar,\,l_{\rm cut}=20~{}\hbar$, and $\varepsilon_{\rm Dcut}=-1700$ MeV. In comparison with RBHF calculation with fixed DWS basis Shen et al. (2016), we have included a pair cut-off for the energy, that is, only pair states $|ab\rangle$ with $\varepsilon_{a}+\varepsilon_{b}\leq\varepsilon_{\rm cut}$ are included in the intermediate states in the BG equation (30) for the channels with $J>J_{h}$. As in the case with the fixed DWS basis, we find good convergence at $\varepsilon_{\rm cut}=1.1$ GeV. Increasing $\varepsilon_{\rm cut}$ to 1.3 GeV will change the energy by $0.8$ MeV and charge radius by $0.005$ fm.

As has been pointed out in Ref. Zhou et al. (2003), the single-particle states with negative energy in the Dirac sea must be included in the basis space for completeness. Because of spherical symmetry, the Dirac equation is solved independently in different blocks with the quantum numbers $\kappa=(l,j)$. We include at least 2 states in the Dirac sea for each block, and then we add more states by gradually decreasing the energy cut-off in the Dirac sea $\varepsilon_{\rm Dcut}$ to check the convergence. The other cut-offs are $J_{\rm cut}=6~{}\hbar,\,l_{\rm cut}=20~{}\hbar$, and $\varepsilon_{\rm cut}=1100$ MeV. The result is shown in Fig. 4. Convergence is achieved at $\varepsilon_{\rm Dcut}=-1700$ MeV. Decreasing $\varepsilon_{\rm Dcut}$ to $-1800$ MeV will change the energy by $0.4$ MeV and charge radius by $0.002$ fm.

Up to now, convergence has been checked with respect to the single-particle basis space and the final choice is $l_{\rm cut}=20~{}\hbar,\,\varepsilon_{\rm cut}=1100$ MeV, and we include at least 2 states in the Dirac sea for each block, and then add more states by choosing the energy cut-off in the Dirac sea $\varepsilon_{\rm Dcut}=-1700$ MeV. For a given DWS potential with a box size $R=7$ fm, there will be 686 single-particle states distributed among 41 blocks (i.e. $0\leq l\leq 20$), where 92 are states with negative energy and 594 are states with positive energy. The dimension of pair states $|ab\rangle$ is different for different channels ($J$, parity $P$, and isospin $T$), for example, for ($0,-,1$) it is 5834, and for ($2,+,0$) it is 54654.

In Fig. 5, we show the convergence with total angular momentum cut-off $J_{\rm cut}$ introduced in Eq. (II.2.2). Increasing $J_{\rm cut}$ from $6$ to $7$ will change the energy by $0.3$ MeV and charge radius by $0.0007$ fm. Therefore $J_{\rm cut}=6$ will be used in the following calculations.

All the above convergence checks were performed at box size $R=7$ fm. In Fig. 6, we show that this is enough for 16O. Further increasing the box size to 8 fm will cause changes by $0.3$ MeV in the energy and by $0.006$ fm in the charge radius.

III.2 Center-of-mass motion

In the above convergence check, the center-of-mass motion in Eq. (24) has not been corrected. Usually there are two ways to treat this term:

Treat it as a first-order correction after the solution of the RHF equation. As discussed in Sec II.1.6, we will label this method as a projection after variation (PAV). 2. 2.

Exclude this term in the RHF equation similarly as in Ref. Chabanat et al. (1998), which we will label as a projection before variation (PBV).

In none of these cases the center-of-mass term is included in the solution of the BG equation (30).

The results are listed in Table 1. It can be seen that the total energy given by PBV is about $9$ MeV smaller than PAV, while the charge radii and center-of-mass correction energy $E_{cm}$ are almost the same for both cases. In order to understand it more clearly, we plot out the total energy at each iteration step in Fig. 7. As can be seen from the figure, there is not so much difference between PBV and PAV during the first RBHF iteration step, where the $G$-matrix $G_{1}$ is calculated from the initial DWS basis. Accordingly, the $G$-matrices $G_{2},\,G_{3},$ and $G_{4}$ are respectively calculated from the convergent RHF basis. One may say that with the same interaction, PAV can be viewed as a good approximation of PBV. But in the next RBHF iteration step, the energy of PBV suddenly becomes smaller than that of PAV.

The reason for this sudden change can be understood from Fig. 8, where the single-particle spectra in $s$ and $p$ blocks are given after the first ($G_{1}$) and last ($G_{4}$) RBHF iteration. It can be seen, the single-particle energies in the RBHF calculation with PBV are generally lower than those of PAV, especially for those high-lying states. In the conventional RHF calculations, the results for PAV and PBV are similar because only the occupied states are concerned and they are similar for PAV and PBV. But for RBHF calculations, the difference in single-particle spectra, in particular high-lying states, will lead to different $G$-matrices in next iteration. Generally, the lower the unoccupied states, the more attractive the $G$-matrix elements of occupied states. As a result, the $G$-matrix in PBV is more attractive and the total binding energy is larger.

PBV and PAV have already been discussed in the non-relativistic BHF calculation Becker et al. (1974). As the single-particle energies of BHF basis states are not fed back into the $G$-matrix, the conclusion in Ref. Becker et al. (1974) was that the two methods give almost the same results for 16O. This is the case in Fig. 7 calculated with the $G$-matrix $G_{1}$. In the fully self-consistent RBHF calculation, however, the total energy given by PBV is about 9 MeV smaller and the single-particle energies are generally lower than those by PAV.

III.3 Self-consistent basis and fixed basis

To take into account the relativistic structure of the Pauli operator in Eq. (30), we adopt a relativistic basis in our calculation. The relativistic DWS basis Zhou et al. (2003) has advantages in comparison with the harmonic oscillator basis Gambhir et al. (1990) like a proper asymptotic behavior of nuclear density distribution, which is crucial for describing, e.g., halo nuclei. More important here is that the nucleon single-particle potential is close to the DWS shape, which serves as a good approximation for the final converged RBHF single-particle states.

The DWS basis is obtained by solving the radial Dirac equation  (26) in potentials of Woods-Saxon shape given in Eqs. (27) and (28). The parameters are chosen with reference to Koepf and Ring (1991) for the neutron potential, with potential depths $V_{0}$ from $-60$ MeV to $-80$ MeV. They are listed in Table 2.

The RBHF calculation has been previously carried out for finite nuclei with a fixed DWS basis in Ref. Shen et al. (2016). We will compare the results with fixed basis and results with the self-consistent RBHF basis in the following.

The total energies at each iteration step are plotted in Fig. 9. Solid symbols stand for the self-consistent RBHF calculations, while open symbols stand for fixed DWS basis calculations. Different shapes represent different DWS potential depths $V_{0}$. Similar to Fig. 7, the first RBHF iteration is represented by $G_{1}$.

We mention that in Fig. 7 or Fig. 9, different calculations may have different number of iteration steps. We have adjusted them slightly in the figures to make them the same without losing too much precision.

For the first RBHF iteration the results are identical to those of the DWS basis calculation for each value of $V_{0}$, because in both cases the same basis is used. After the first RBHF iteration, there appears a difference between the self-consistent calculation and the fixed basis one. In the end, the results of the self-consistent calculation do not depend on the initial DWS basis, while the results of a fixed basis one do. As expected, self-consistency is very important to get an unambiguous result.

III.4 Single-particle potential $U$ of particle states

One uncertainty in (R)BHF theory is the choice of the single-particle potential for particle states in Eq. (18). Different choices have been discussed in BHF for finite nuclei and the readers are referred to Ref. Davies et al. (1969) for more details. Here we choose Eq. (18) for $U_{ab}$ where $\varepsilon^{\prime}_{a}$ is different from the single-particle energy $\varepsilon_{a}$, with $\varepsilon^{\prime}_{a}=\varepsilon^{\prime}_{b}=\varepsilon^{\prime}$ as the single-particle energy of hole states, and compare the results between the lowest value $\varepsilon^{\prime}=\varepsilon_{\nu 1s1/2}$ and the highest value $\varepsilon^{\prime}=\varepsilon_{\pi 1p1/2}$. We also give as a comparison the results for the gap choice: $U_{ab}=0$. These results are presented in Table 3.

It can be seen that the gap choice gives the least bound system. While fixing $\varepsilon^{\prime}$ as an energy among the occupied states gives $1.4$ to $1.7$ MeV per nucleon more binding.

In the framework of Bethe-Brueckner-Goldstone (BBG) (see for instance Ref. Day (1967)), the gap choice and the continuous choice have been discussed in nuclear matter in Ref. Song et al. (1998). The continuous choice gives more binding than the gap choice at the BHF level. By performing the BBG expansion up to the three hole-line level, the above two choices give similar results, and lie between the results of these two choices at the BHF level Song et al. (1998). From this point of view, the choice of $\varepsilon^{\prime}$ as an energy among the occupied states is a reasonable choice.

In Fig. 10 we present the single-particle spectrum of RBHF calculations with different choices and we find: First, in the gap choice there is a larger gap between the unoccupied states and occupied states as compared to other choices. This has been observed in earlier non-relativistic calculations and therefore this name was chosen Day (1967). Second, as compared to the gap choice, the choices of $\varepsilon^{\prime}$ fixed as an energy among the occupied states give more attraction to low-lying states. This can be understood because the single-particle energies $\varepsilon$ of low-lying particle states are close to the chosen $\varepsilon^{\prime}$. Therefore it can be viewed as the continuous choice and the $G$-matrix is attractive. Third, the single-particle energies of low-lying states with choices of $\varepsilon^{\prime}=\varepsilon_{\pi 1p1/2}$ and $\varepsilon^{\prime}=\varepsilon_{\nu 1s1/2}$ are close to each other.

Summing up the previous discussions, the numerical details for the solution of RBHF equations in the following applications include $l_{\rm cut}=20~{}\hbar$ (Fig. 2), $\varepsilon_{\rm cut}=-1.1$ GeV (Fig. 3), $\varepsilon_{\rm Dcut}=-1700$ MeV (Fig. 4), $J_{\rm cut}=6~{}\hbar$ (Fig. 5), and $R_{\rm box}=7$ fm (Fig. 6). The center-of-mass term is treated with PBV. For the $pp$-potential $U_{ab}$ in Eq. (18), we consider a continuous choice with $\varepsilon^{\prime}$ as the last occupied state in the Fermi sea. This means $\varepsilon^{\prime}=\varepsilon_{\pi 1p1/2}$ for 16O. For the nucleus 40Ca we use $\varepsilon^{\prime}=\varepsilon_{\pi 1d3/2}$, $l_{\rm cut}=25~{}\hbar,J_{\rm cut}=9~{}\hbar$, and the same values for the remaining quantities.

IV Results and discussion

IV.1 The nucleus 16O

The total energy, charge radius $r_{c}$, matter radius $r_{m}$, and proton spin-orbit splitting for the $1p$ shell of 16O calculated by RBHF with the interaction Bonn A Machleidt (1989) are listed in Table 4, in comparison with experimental data Wang et al. (2012); Angeli and Marinova (2013); Ozawa et al. (2001); Coraggio et al. (2003). The corresponding results from RBHF in fixed DWS basis Shen et al. (2016), Density-Dependent Relativistic Hartree-Fock (DDRHF) with PKO1 Long et al. (2006) and PKA1 Long et al. (2007), non-relativistic BHF Hu et al. (2017) with $V_{{\rm low}-k}$ derived from Argonne $v_{18}$ Wiringa et al. (1995), Coupled-Cluster (CC) method Hagen et al. (2009) and No Core Shell Model (NCSM) Roth et al. (2011) with N3LO Entem and Machleidt (2003), and Nuclear Lattice Effective Field Theory (NLEFT) Lähde et al. (2014) with N2LO Epelbaum et al. (2009b) are also included.

The phenomenological functional PKO1 includes tensor correlations induced by the $\pi$-exchange and PKA1 includes in addition tensor correlations by the $\rho$-exchange. As expected, the phenomenological results, which are both fitted to the experimental data Long et al. (2006, 2007), are in much better agreement with data than all the ab initio results.

The total energy of our self-consistent RBHF calculation is underbound by 14.1 MeV (or by 11%) and the charge radius is smaller by 0.14 fm (or by 5%) as compared to the experimental values. This is similar to the fixed basis calculation Shen et al. (2016) and consistent with the infinite nuclear matter results Brockmann and Machleidt (1990). In the self-consistent non-relativistic BHF calculations Hu et al. (2017), the interaction $V_{{\rm low}-k}$ derived from the Argonne interaction $v_{18}$ was used. The result shows an overbinding by $6.6$ MeV (or by 5%) and the rms radius is too small. In comparison, the values of $E=-120.9$ MeV is obtained with the CC Hagen et al. (2009) using the chiral $NN$ interaction N3LO, $E=-119.7(6)$ MeV obtained with the NCSM Roth et al. (2011) using the same interaction N3LO, and $E=-121.4(5)$ obtained with the NLEFT Lähde et al. (2014) using the interaction N2LO.

The experimental matter radius $r_{m}=2.54$ fm is larger than $r_{m}=2.30$ fm calculated with the CC method. It has been shown that this is a general feature of all the ab initio calculations with conventional forces based on chiral effective field theory Lapoux et al. (2016). Only modern chiral forces Ekström et al. (2015), which include also radii in the adjustment of the parameters for the bare forces, are able to cure this problem.

The spin-orbit splittings of $1p$ proton shell in RBHF theory $\Delta E_{\pi 1p}^{ls}=5.4$ MeV is smaller than the previous RBHF results with fixed DWS basis $\Delta E_{\pi 1p}^{ls}=6.0$ MeV, and the deviation with the data is respectively $14\%$ and $5\%$.

Figure 11 shows the energy per nucleon $E/A$ for 16O as a function of the charge radius $r_{c}$ calculated in RBHF using the interactions Bonn A, B, and C, in comparison with BHF and the relativistic effective density approximation (EDA) Müther et al. (1990). It can be seen in all cases that relativistic effects improve the results considerably. By comparing EDA with RBHF, one can see that self-consistency has important effect.

In Table 5, the total energy, charge radius, single-particle energies, and $\pi 1p$ spin-orbit splitting of 16O calculated by RBHF theory with the interactions Bonn A, B, and C are listed, and in comparison with the experimental data. The experimental single-particle energies are taken from Ref. Coraggio et al. (2003).

Figure 12 shows the charge density distributions of 16O calculated by RBHF theory with the interactions Bonn A, B, and C, in comparison with experimental data. The experimental data is from Ref. De Vries et al. (1987). It can be seen Bonn A, B, and C interactions give too large central distribution, as a result the charge radius is smaller than the experimental value.

IV.2 The nuclei 4He and 40Ca

The total energy, charge radius, and proton radius of 4He calculated by RBHF theory using the interaction Bonn A Machleidt (1989) with PBV and PAV are listed in Table 6, in comparison with experimental data Wang et al. (2012); Angeli and Marinova (2013); Lu et al. (2013). The corresponding results from DDRHF with PKO1 Long et al. (2006) and PKA1 Long et al. (2007), solution of the Faddeev-Yakubovsky (FY) equation Nogga et al. (2000) with CD-Bonn Machleidt et al. (1996), FY Binder et al. (2016) with N4LO Epelbaum et al. (2015), NCSM Navrátil (2007) with N3LO Entem and Machleidt (2003), NLEFT Lähde et al. (2014) with N2LO, and BHF Hu et al. (2017) with $V_{{\rm low}-k}$ derived from Argonne $v_{18}$ Wiringa et al. (1995) are also included.

It can be seen that by considering only the two-body interaction, all the non-relativistic calculations predicated underbinding for 4He. In contrast, RBHF with PBV gives too much binding and a too large radius. As expected, the center-of-mass correction plays an important role for such a light nucleus: PBV and PAV give very different results.

In Table 7, the total energy, charge radius, matter radius, and proton spin-orbit splitting for the $1d$ shell of 40Ca calculated by RBHF with the interaction Bonn A Machleidt (1989) are listed, in comparison with experimental data Wang et al. (2012); Angeli and Marinova (2013); Coraggio et al. (2003). The corresponding results from DDRHF with PKO1 Long et al. (2006) and PKA1 Long et al. (2007), BHF Hu et al. (2017), NCSM Roth and Navrátil (2007), and CC Hagen et al. (2007) with $V_{{\rm low}-k}$ derived from Argonne $v_{18}$ Wiringa et al. (1995), and CC Hagen et al. (2010) with N3LO Entem and Machleidt (2003) are also included.

Similar as for 16O, the total energy of 40Ca calculated by RBHF with the interaction Bonn A is underbound by 51.3 MeV (or by 15%) and the charge radius is smaller by 0.25 fm (or by 7%) as compared to the experimental values. For the proton $1d$ spin-orbit splitting, RBHF with Bonn A gives a very good description for the data. Most of the non-relativistic results, because of the missing three-body force, give too large binding energy and too small radius, except the CC method with N3LO which reproduces the experimental binding energy well.

V Summary

The relativistic Brueckner-Hartree-Fock equations have been solved for finite nuclei in a self-consistent relativistic Hartree-Fock basis. For this purpose a basis transformation is performed in each step of the iteration, and the $G$-matrix is obtained by solving the Bethe-Goldstone equation in this self-consistent RHF basis. The Pauli operator has been taken into account fully self-consistently without any approximation. Relativistic versions of the bare $NN$ interactions Bonn A, B, and C Machleidt (1989) have been used.

Taking 16O as an example, the relevant convergence properties in the RBHF calculation have been checked, including a single-particle angular momentum cut-off $l_{\rm cut}$, a single-particle energy cut-off $\varepsilon_{\text{cut}}$, a single-particle energy cut-off in the Dirac sea $\varepsilon_{\text{Dcut}}$, a total angular momentum cut-off $J_{\rm cut}$, and a box size $R$.

The binding energy calculated by the self-consistent RBHF does not depend on the initial basis, and is generally smaller than the results of RBHF calculations in the fixed DWS basis. This underbinding is, at a first glance, unexpected, because in conventional non-relativistic Hartree-Fock theory the variational principle yields the lowest energy for the fully self-consistent calculation. However, this can be understood by the well known fact that the Brueckner theory does not obey the variational principle: the matrix-elements of the $G$-matrix depend on the Pauli operator and on the single-particle energies, and the single-particle potential is not defined by the variational principle as the Hartree-Fock theory does.

Two different approaches for the treatment of the center-of-mass motion have been discussed: approximate projection before variation (PBV) and after variation (PAV). While the two approaches give similar results for RHF calculations, the total energy of 16O given by PBV is about 9 MeV smaller than that of PAV for RBHF calculations. This is due to the fact that the single-particle energies are different in the two approaches, especially for those high-lying states, which have a strong impact on the self-consistent $G$-matrix. Thus, the more consistent approach PBV should be used in the RBHF framework.

We have also discussed different choices for the single-particle potential of particle states, which brings a major ambiguity in the (R)BHF framework. This ambiguity, in the more general framework of the hole-line expansion, is connected to the limitation of two hole-lines used here. We have chosen Eq. (18) for the single-particle potential of particle states, and found that the choice of $\varepsilon^{\prime}$ in Eq. (18) as an energy among the occupied states is a reasonable choice.

We have performed RBHF calculations for the doubly closed shell nuclei 4He, 16O, and 40Ca, and the results are compared with experimental data, phenomenological DDRHF calculations, and other state-of-the-art ab initio calculations. For 4He, because of its very light nature, the treatment of the center-of-mass motion is essential. The total energy of 16O and 40Ca calculated by RBHF with the interaction Bonn A is underbound by 1 MeV per nucleon and the charge radius is smaller by 5% to 7% as compared to the experimental values. The spin-orbit splittings for both 16O and 40Ca are very well reproduced.

Except 4He, we find considerable underbinding by RBHF for the heavier nuclei. There remains the open question of the origin of this underbinding.

First, we cannot exclude, at present, the importance of additional three-body forces, even in a relativistic description. In the relativistic framework, three-body forces resulting from virtual nucleon-antinucleon excitations ($Z$-diagram) are included Brown et al. (1987). There are also other non-relativistic origins for three-body forces in nuclei, as for instance, the Fujita-Miyazawa force Fujita and Miyazawa (1957). Microscopic investigations of various contributions to the three-body force have shown that for lower densities, up to the saturation density in nuclear matter, the relativistic $Z$-diagrams, which are included in our RBHF calculations, play a major role Zuo et al. (2002). However, they produce too much repulsion, i.e., one should expect underbinding in finite nuclei, as we find in 16O and 40Ca. As shown in Ref. Zuo et al. (2002) additional three-body forces of non-relativistic nature produce attraction and their contributions are not negligible. For lower densities they play a less important role, but details and more quantitative conclusions are still open questions for the future.

Second, we stayed in these calculations completely on the Brueckner-Hartree-Fock level, i.e. we did not include higher-order diagrams with more than two hole-lines in the hole-line expansion. They are known to have certain contributions in non-relativistic nuclear matter investigations Rajaraman and Bethe (1967). Among them, the third-order saturation-potential diagrams (or rearrangement diagrams) can be taken into account by the so-called renormalized BHF approach Becker (1970) and it is definitely interesting to investigate its relativistic extension in the future.

The current work is intended to establish a firm ground for a relativistic ab initio framework in finite nuclei. With recent progress in covariant chiral nucleon-nucleon interaction Ren et al. (2016) and hyperon-nucleon interaction Li et al. (2016), we are looking forward to study nuclear many-body problem with RBHF based on chiral effective field theory. The ultimate goal is to extend it to heavy nuclei and help us to understand nuclear structure in a microscopic way. We hope to learn from such ab initio calculations and to settle down some of the open questions in modern nuclear energy density functional theory, such as the isospin dependence or the importance of the tensor terms Lalazissis et al. (2009); Long et al. (2007).

ACKNOWLEDGMENTS

We thank Jinniu Hu for discussions and Wenhui Long for providing his RHF code. This work was partly supported by the Major State 973 Program of China No. 2013CB834400, Natural Science Foundation of China under Grants No. 11335002, No. 11375015, and No. 11621131001, the Overseas Distinguished Professor Project from Ministry of Education No. MS2010BJDX001, the Research Fund for the Doctoral Program of Higher Education under Grant No. 20110001110087, and the DFG cluster of excellence “Origin and Structure of the Universe” (www.universe-cluster.de). SS would like to thank the short-term PhD student exchange program of Peking University and the RIKEN IPA project, and HL would like to thank the RIKEN iTHES project and iTHEMS program.

Appendix A Two-body matrix elements

Here we give the details for the calculation of relativistic two-body matrix elements for the Bonn interaction Machleidt (1989), which is defined in terms of meson-exchange terms:

[TABLE]

with $\alpha$ for scalar ($\sigma,\delta$), vector ($\omega,\rho$), and pseudoscalar ($\eta,\pi$) mesons. The interacting vertices are given in Eq. (5). Here we discuss the isoscalar mesons only since the isovector mesons are identical except for the isospin matrix elements, which can be considered separately.

Using the propagator in Eq. (8) in momentum space, we can rewrite the two-body matrix elements as

[TABLE]

The form factor in Eq. (6) depends only on the meson $\alpha$ and the absolute momentum transfer $|\mathbf{q}|$, so that it can be added at the last step.

The plane wave can be expanded as

[TABLE]

where $j_{L}(qr)$ is the spherical Bessel function, $Y_{LM}$ is the spherical harmonic function, and $\mathbf{\hat{q}}$ represents the angular part of the vector $\mathbf{q}$. The integration over $\hat{\mathbf{q}}$ yields

[TABLE]

For spherical nuclei we consider matrix elements coupled to good angular momentum. There are several ways to couple pairs of indices

• •

$pp$-coupled matrix elements:

[TABLE]

• •

$ph$-coupled matrix elements (direct term):

[TABLE]

• •

$ph$-coupled matrix elements (exchange term):

[TABLE]

These different coupling schemes are related to each other by the following recoupling rules de Shalit and Talmi (1963)

[TABLE]

and the inverse relation:

[TABLE]

The antisymmetrized $ph$-coupled matrix elements are obtained as

[TABLE]

Expressing the various vertices $\gamma^{0}\Gamma_{\alpha}Y_{LM}$ in Eq. (42) in terms of spherical tensor operators $\hat{O}_{\lambda\mu}$, we can use the Wigner-Eckart theorem Varshalovich et al. (1988) and express the $ph$-coupled matrix elements by the reduced matrix elements of the corresponding tensors

[TABLE]

where $\hat{J}=\sqrt{2J+1}$. Therefore we present in the following the direct terms of the $ph$-coupled matrix elements for the different mesons.

A.1 Scalar meson

For the scalar meson,

[TABLE]

where

[TABLE]

with $F_{a}(r)$ and $G_{a}(r)$ are the large and small components of the Dirac spinor in Eq. (25).

A.2 Pseudoscalar meson

For the pseudoscalar meson,

[TABLE]

where $[Y_{L}\sigma]_{I}$ is the reduced form of

[TABLE]

The matrix elements read

[TABLE]

A.3 Vector meson

The expression for the vector meson is more complicated. According to the Lorentz structure, we divide it into a time component ($v0$) and a space component ($v1$). Take the vector-vector coupling (vv) as an example,

[TABLE]

The results for the time component are

[TABLE]

for the vector-vector, tensor-tensor (tt), and vector-tensor (vt) components, respectively.

The results for the space component are (with $\bm{\alpha}=\gamma^{0}\bm{\gamma}$):

[TABLE]

with

[TABLE]

The reduced matrix element reads

[TABLE]

Notice that in the Bonn interaction Machleidt (1989) the tensor part of $\omega$ meson is deemed to be small and omitted, only the $\rho$ meson has a tensor part $f_{\rho}$.

A.4 Isospin matrix elements

The isospin operator has been separated in the vertex expression in Eq. (5). For the isovector mesons ($\delta,\pi,\rho$), the full interaction vertices are accompanied with the isovector operator as

[TABLE]

For the isoscalar operator $\mathbb{I}$, the non-zero matrix elements are

[TABLE]

In contrast, for the isovector operator $\vec{\tau}_{1}\cdot\vec{\tau}_{2}$, the non-zero matrix elements are

[TABLE]

with the conventions $\tau_{z}|n\rangle=|n\rangle$ and $\tau_{z}|p\rangle=-|p\rangle$ for neutron and proton states.

Appendix B Center-of-mass Motion

B.1 Matrix element for center-of-mass motion

In second quantization the operator for the center of mass correction is given by

[TABLE]

It contains a one-body operator $T_{cm}$ and a two-body operator $V_{cm}$.

The matrix element of the one-body term $T_{cm}$ is similar to the non-relativistic kinetic energy. Because of spherical symmetry and parity conservation, it is diagonal in the block index $\kappa=(l,j)$ and can be derived as

[TABLE]

The matrix element of the two-body term is

[TABLE]

In this equation, summation indices $\eta,\,\zeta$ sum over the upper ($+$) and lower ($-$) components, and $A,\,B$ are defined as

[TABLE]

Here we also give the antisymmetrized $jj$-coupled form of $V_{cm}$. In the spherical RHF equation, only $ph$ coupling $I=0$ matrix elements in Eq. (44) are needed,

[TABLE]

B.2 Center-of-mass correction for total energy

The energy of the center-of-mass motion is

[TABLE]

The first term can be obtained from Eq. (81),

[TABLE]

The second term can be obtained from Eq. (B.1),

[TABLE]

Indeed, the direct contribution in this expression vanishes. Note that the $\sum^{\prime}$ does not sum over magnetic quantum number $m$. This energy $E_{cm}$ should be subtracted in the total energy expression in Eq. (36).

Appendix C Coulomb Interaction

The Coulomb interaction can be derived from the interaction vertex,

[TABLE]

In the R(B)HF framework, its contribution to the single-particle potential is

[TABLE]

which is not included in the G-matrix and calculated separately.

C.1 Direct term

The direct term of Coulomb field in the spherical case reads

[TABLE]

with $\theta$ the step functions. The contribution to total energy is

[TABLE]

C.2 Exchange term

In principle, the exchange term of Coulomb field is non-local. For simplicity, in the present calculations, we use the relativistic local density approximation (RLDA) for the Coulomb exchange term Gu et al. (2013), which in the spherical case reads

[TABLE]

This is the well-known Slater approximation Slater (1951) plus a relativistic correction. The contribution to total energy is

[TABLE]

Appendix D Extrapolation with Respect to Box Size

For fixed $l_{\rm cut}$ and $\varepsilon_{\rm cut}$, the basis space increases dramatically as $R$ increases. As shown in the coupled-cluster calculations with the harmonic oscillator basis, when the ultraviolet condition is fulfilled, the total energy and radius have specific relations with the box size $R$ due to the infrared cutoffs induced by the box boundary condition Furnstahl et al. (2012)

[TABLE]

where $\beta=2k_{\infty}R$, and $E_{\infty}$ and $\langle r^{2}\rangle_{\infty}$ are the expectation values of total energy and radius when the box size $R$ is extrapolated to infinity. Together with $k_{\infty},a_{0},c_{0},c_{1}$, they are fitted to several calculated points. Now we can evaluate these relations within the RBHF framework.

We perform the RBHF calculations for 16O with several different box sizes from $R=4$ to $R=7$ fm. The results are shown in Fig. 13, where squares symbols represent the results calculated by RBHF theory with Bonn A interaction. The extrapolation of total energy and charge radius is carried out by using the calculated points from $R=4$ to $5.5$ fm indicated with the solid symbols. The corresponding uncertainties are evaluated by the Jackknife resampling method Efron (1982). It is clearly seen that the calculated points with open symbols, which are not included in the fit of Eq. (110), are well located on the lines of extrapolation.

The extrapolated values of total energy and charge radius of 16O by the RBHF theory with Bonn A interaction are

[TABLE]

where the uncertainty comes from Jackknife resampling.