Towards an ab initio covariant density functional for nuclear structure

TL;DR

This paper reviews covariant density functional theory (CDFT) in nuclear physics, discusses its challenges, and explores developing an ab initio relativistic density functional based on realistic nucleon-nucleon interactions for improved nuclear structure predictions.

Contribution

It introduces the concept of constructing an ab initio covariant density functional for nuclear structure from realistic nucleon-nucleon interactions, advancing beyond phenomenological models.

Findings

Review of covariant density functional theory success in nuclear property prediction

Discussion of challenges in predicting unstable nuclei and high-density matter

Overview of current relativistic Brueckner-Hartree-Fock calculations for finite nuclei

Abstract

Nuclear structure models built from phenomenological mean fields, the effective nucleon-nucleon interactions (or Lagrangians), and the realistic bare nucleon-nucleon interactions are reviewed. The success of covariant density functional theory (CDFT) to describe nuclear properties and its influence on Brueckner theory within the relativistic framework are focused upon. The challenges and ambiguities of predictions for unstable nuclei without data or for high-density nuclear matter, arising from relativistic density functionals, are discussed. The basic ideas in building an ab initio relativistic density functional for nuclear structure from ab initio calculations with realistic nucleon-nucleon interactions for both nuclear matter and finite nuclei are presented. The current status of fully self-consistent relativistic Brueckner-Hartree-Fock (RBHF) calculations for finite nuclei or neutron drops (ideal systems composed of a finite number of neutrons and confined within an external field) is reviewed. The guidance and perspectives towards an ab initio covariant density functional theory for nuclear structure derived from the RBHF results are provided.

Tables

Table 1. Table 1: Saturation properties of symmetric nuclear matter calculated by RBHF using the Bonn A, B, and C interactions [ 165 ] by different methods: Saturation density ρ 0 subscript 𝜌 0 \rho_{0} in fm -3 and energy per nucleon E / A 𝐸 𝐴 E/A in MeV.

Bonn A		Bonn B		Bonn C
${\rho }_0$	$E/A$	${\rho }_0$	$E/A$	${\rho }_0$	$E/A$	Description	Ref.
$0.174$	$-16.49$	$0.172$	$-15.73$	$0.170$	$-14.38$	Full basis, momentum-averaged s.p. potential	[189]
$0.174$	$-15.72$	$0.170$	$-14.81$	$0.162$	$-13.73$	Full basis, momentum-dependent s.p. potential	[189]
$0.149$	$-10.5$	$0.130$	$-7.3$	$0.112$	$-5.2$	Full basis	[169]
$0.181$	$-16.15$	$0.163$	$-14.59$	$0.145$	$-13.69$	Projection, $ps$ for $T_{\mathrm{Sub}}$	[285]
$0.181$	$-15.72$	$0.159$	$-13.99$	$0.138$	$-13.00$	Projection, complete $pv$ for $T_{\mathrm{Sub}}$	[285]
$0.189$	$-15.81$	$0.166$	$-13.70$	$0.148$	$-12.31$	Projection, conventional $pv$	[285]
$0.185$	$-15.59$	$0.174$	$-13.60$	$0.155$	$-12.26$	Momentum-dependence analysis	[160]

Table 2. Table 2: Bulk parameters (see text for definitions) of symmetric and asymmetric nuclear matter at saturation density ρ 0 subscript 𝜌 0 \rho_{0} obtained by RBHF theory using the Bonn A, B, and C interactions, with the exact and averaged center-of-mass (c.m.) momentum. The quantities Δ Δ \Delta are defined as the differences between the exact and the averaged treatments of the c.m. momentum. Results are compared with those obtained by non-relativistic BHF theory with and without three-body forces (TBF)  [ 160 , 292 , 293 ] and with empirical values  [ 295 ] . Table taken from Ref.  [ 290 ] .

Model Potential			${\rho }_0$	$E/A$	$K_{\infty }$	$Q_0$	$M_0$	$E_{\mathrm{sym}}$	$L$	$K_{\mathrm{sym}}$	$K_{\mathrm{asy}}$	$K_{\tau }$	$K_{\mathrm{Coul}}$
			(fm-3)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)
		exact	$0.180$	$-15.38$	$286$	$731$	$4163$	$33.7$	$75.8$	$-57.0$	$-512$	$-705$	$-8.30$
	A	average	$0.182$	$-15.04$	$289$	$650$	$4118$	$32.6$	$74.7$	$-53.1$	$-501$	$-669$	$-8.09$
		$\Delta$	$-0.002$	$-0.34$	$-3$	$81$	$45$	$1.1$	$1.1$	$-3.9$	$-11$	$-36$	$-0.21$
		exact	$0.164$	$-13.44$	$222$	$547$	$3211$	$29.9$	$63.0$	$-56.3$	$-434$	$-590$	$-7.98$
RBHF	B	average	$0.165$	$-13.08$	$220$	$791$	$3431$	$28.7$	$65.3$	$-47.5$	$-439$	$-674$	$-8.86$
		$\Delta$	$-0.001$	$-0.36$	$2$	$-244$	$-220$	$1.2$	$-2.3$	$-8.8$	$5$	$84$	$0.88$
		exact	$0.149$	$-12.12$	$176$	$260$	$2372$	$26.8$	$51.7$	$-55.6$	$-366$	$-442$	$-7.00$
	C	average	$0.150$	$-11.75$	$168$	$638$	$2654$	$25.6$	$58.8$	$-41.1$	$-394$	$-618$	$-8.74$
		$\Delta$	$-0.001$	$-0.37$	$8$	$-378$	$-282$	$1.2$	$-7.1$	$-14.5$	$28$	$176$	$1.74$
	A		$0.428$	$-23.55$	$204$			$32.1$
BHF	B		$0.309$	$-18.30$	$160$			$31.8$
	C		$0.247$	$-15.75$	$143$			$28.5$
	AV18	no TBF	$0.240$	$-17.30$	$214$	$-225$	$2343$	$35.8$	$63.1$	$-27.8$	$-406$	$-340$	$-6.01$
	AV18	TBFa	$0.187$	$-15.23$	$196$	$-281$	$2071$	$34.3$	$66.5$	$-31.3$	$-430$	$-335$	$-5.23$
	AV18	TBFb	$0.176$	$-14.62$	$186$	$-225$	$2007$	$33.6$	$66.9$	$-23.4$	$-425$	$-344$	$-5.30$
Empirical			$0.155$	$-15.8$	$230$	300		$32$	$60$	$-100$		$-400$
			$\pm 0.005$	$\pm 0.3$	$\pm 20$	$\pm 400$		$\pm 2$	$\pm 15$	$\pm 100$		$\pm 100$

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

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

Table 4. Table 4: Binding energy per nucleon ( B / A 𝐵 𝐴 B/A in MeV), charge radius ( r c subscript 𝑟 𝑐 r_{c} in fm), and proton 1 ​ p 1 𝑝 1p spin-orbit splitting ( Δ ​ E π ​ 1 ​ p l ​ s Δ superscript subscript 𝐸 𝜋 1 𝑝 𝑙 𝑠 \Delta E_{\pi 1p}^{ls} in MeV) of 16 O in the RBHF calculations with different choices for U a ​ b subscript 𝑈 𝑎 𝑏 U_{ab} . Formula I for Eq. ( 74 ), and formula II for Eq. ( 57 ). The experimental data are from [ 343 , 344 , 345 ] .

	I, Gap	I, $e^{\prime }=e_{\nu 1s1/2}$	I, $e^{\prime }=e_{\pi 1p1/2}$	II, $e^{\prime }=e_{\pi 1p1/2}$	Exp.
$B/A$	$5.41$	$6.88$	$7.10$	$7.51$	$7.98$
$r_c$	$2.64$	$2.57$	$2.56$	$2.53$	$2.70$
$\Delta E_{\pi 1p}^{ls}$	$5.4$	$5.4$	$5.4$	$5.3$	$6.3$

Table 5. Table 5: Total energy, charge radius, matter radius, and π ​ 1 ​ p 𝜋 1 𝑝 \pi 1p spin-orbit splitting of 16 O calculated by RBHF theory [ 37 , 232 ] with the interactions Bonn A, B, and C, in comparison with experimental data. The corresponding results from DDRHF with effective interactions PKO1 [ 212 ] and PKA1 [ 207 ] , non-relativistic BHF  [ 348 ] with V low − k subscript 𝑉 low 𝑘 V_{{\rm low}-k} derived from Argonne v 18 subscript 𝑣 18 v_{18} , CC method [ 349 ] , IM-SRG [ 350 ] , NCSM [ 351 ] , and SCGF method [ 352 ] with N 3 LO, NLEFT  [ 353 ] with N 2 LO, and QMC method [ 355 ] with local chiral force N 2 LO are also included.

	$E$ (MeV)	$r_c$ (fm)	$r_m$ (fm)	$\Delta E_{\pi 1p}^{ls}$ (MeV)	Ref.
Exp.	$-127.6$	$2.70$	$2.54\left(2\right)$	$6.3$	[343, 344, 347, 345]
RBHF, Bonn A	$-120.2$	$2.53$	$2.39$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5.3}$	[37]
RBHF, Bonn B	$-107.1$	$2.59$	$2.45$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}4.5}$	[37]
RBHF, Bonn C	$-98.0$	$2.64$	$2.50$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3.9}$	[37]
DDRHF, PKO1	$-128.3$	$2.68$	$2.54$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}6.4}$	[212]
DDRHF, PKA1	$-127.0$	$2.80$	$2.67$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}6.0}$	[207]
BHF, AV18	$-134.2$		$1.95$	$13.0$	[348]
CC, N3LO	$-120.9$		$2.30$		[349]
IM-SRG, N3LO	$-122.71\left(9\right)$				[350]
NCSM, N3LO	$-119.7\left(6\right)$				[351]
SCGF, N3LO	$-122$				[352]
NLEFT, N2LO	$-121.4\left(5\right)$				[353]
QMC, N2LO	$-87\left(3\right)\left(11\right)$	$2.76\left(5\right)\left(12\right)$			[355]

Table 6. Table 6: Total energy, charge radius, and proton radius of 4 He calculated by RBHF theory using the interaction Bonn A with PBV and PAV  [ 37 ] , in comparison with experimental data  [ 343 , 344 , 360 ] . The corresponding results from DDRHF with PKO1 [ 212 ] and PKA1 [ 207 ] , the FY equation  [ 361 ] with CD-Bonn, FY  [ 346 ] with N 4 LO, NCSM  [ 363 ] with N 3 LO, NLEFT  [ 353 ] with N 2 LO, and BHF  [ 348 ] 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)	Ref.
Exp.	$-28.30$	$1.68$	$1.46$	[343, 344, 360]
RBHF (PBV), Bonn A	$-35.05$	$1.83$	$1.64$	[37]
RBHF (PAV), Bonn A	$-26.31$	$1.90$	$1.73$	[37]
DDRHF, PKO1	$-28.45$	$1.90$	$1.72$	[212]
DDRHF, PKA1	$-28.28$	$2.06$	$1.90$	[207]
FY, CD-Bonn	$-26.26$			[361]
FY, N4LO	$-24.27\left(6\right)$		1.547(2)	[346]
NCSM, N3LO	$-25.39\left(1\right)$		$1.515\left(2\right)$	[363]
NLEFT, N2LO	$-25.60\left(6\right)$			[353]
BHF, AV18	$-25.90$			[348]

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 [ 37 , 232 ] with the interaction Bonn A, in comparison with experimental data  [ 343 , 344 , 345 ] . The corresponding results from DDRHF with PKO1  [ 212 ] and PKA1  [ 207 ] , BHF  [ 348 ] , NCSM [ 364 ] , and CC [ 365 ] with V low − k subscript 𝑉 low 𝑘 V_{{\rm low}-k} derived from Argonne v 18 subscript 𝑣 18 v_{18} , CC [ 366 ] and IM-SRG [ 350 ] with N 3 LO are also included.

	$E$ (MeV)	$r_c$ (fm)	$r_m$ (fm)	$\Delta E_{\pi 1d}^{ls}$ (MeV)	Ref.
Exp.	$-342.1$	$3.48$		$6.6\pm 2.5$	[343, 344, 345]
RBHF, Bonn A	$-306.1$	$3.22$	$3.10$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5.8}$	[37]
DDRHF, PKO1	$-343.3$	$3.44$	$3.33$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}6.6}$	[212]
DDRHF, PKA1	$-341.7$	$3.53$	$3.41$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}7.2}$	[207]
BHF, AV18	$-552.1$		$2.20$	$24.9$	[348]
NCSM, AV18	$-461.8$		$2.27$		[364]
CC, AV18	$-502.9$				[365]
CC, N3LO	$-345.2$				[366]
IM-SRG, N3LO	$-358.5$				[350]

Table 8. Table 8: Rms radius R N subscript 𝑅 𝑁 R_{N} of N = 50 𝑁 50 N=50 neutron drop in a HO trap ( ℏ ​ ω = 10 Planck-constant-over-2-pi 𝜔 10 \hbar\omega=10 MeV) calculated by RBHF theory using the interaction Bonn A. The symmetry energy a sym subscript 𝑎 sym a_{\rm sym} and slope parameter L 𝐿 L calculated in nuclear matter [ 414 , 184 , 415 ] are also listed, and they are compared with the results of relativistic functionals NL3  [ 104 ] , PK1  [ 107 ] , DD-ME2  [ 114 ] , PKDD  [ 107 ] , PC-PK1  [ 117 ] , and PKO1  [ 206 ] .

	$R_{N=50}$ (fm)	$a_{\mathrm{sym}}$ (MeV)	$L$ (MeV)
Bonn A	$3.61$	$34.8$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}71}$
NL3	$4.04$	$36.6$	$119$
PK1	$4.04$	$37.6$	$116$
DD-ME2	$3.72$	$32.3$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}51}$
PKDD	$3.99$	$36.8$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}90}$
PC-PK1	$4.08$	$35.6$	$113$
PKO1	$3.90$	$34.4$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}98}$

Equations

$$E_U[\rho ]=F_{\mathrm{HK}}[\rho ]+\int d^{3}r{U}_{\mathrm{ext}}(\mathbf{r})\rho (\mathbf{r}),$$

$${\rho }_{\mathrm{gs}}(\mathbf{r})=\sum _{i=1}^N|{\varphi }_i(\mathbf{r}){|}^2,$$

$$\left[-\frac{{\hslash }^2}{2M}{\nabla }^2+v_{\mathrm{KS}}(\mathbf{r})\right]{\varphi }_i(\mathbf{r})={\epsilon }_i{\varphi }_i(\mathbf{r}).$$

$$F_{\mathrm{HK}}[\rho ]=T_s[\rho ]+E_H[\rho ]+E_{xc}[\rho ].$$

$$E[\rho ]=⟨\Phi |T+V[\rho ]|\Phi ⟩,$$

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

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

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

$$|{\Phi }_0⟩=\prod _k^A{a}_k^{†}|0⟩,$$

$$E^{\mathrm{HF}}=⟨{\Phi }_0|H|{\Phi }_0⟩=\sum _k^A⟨k|T|k⟩+\frac{1}{2}\sum _{kl}^A⟨kl|\bar{V}|kl⟩,$$

$$(T+U^{\mathrm{HF}})|k⟩=e_k|k⟩,$$

$$⟨k^{\prime }|U^{\mathrm{HF}}|k⟩=\sum _l^A⟨k^{\prime }l|\bar{V}|kl⟩.$$

$$⟨k^{\prime }{l}^{\prime }|\bar{G}(W)|kl⟩=⟨k^{\prime }{l}^{\prime }|\bar{V}|kl⟩+\frac{1}{2}\sum _{mn}⟨k^{\prime }{l}^{\prime }|\bar{V}|mn⟩\frac{Q(m,n)}{W-e_m-e_n}⟨mn|\bar{G}(W)|kl⟩,$$

$$E^{\mathrm{BHF}}=\sum _k^A⟨k|T|k⟩+\frac{1}{2}\sum _{kl}^A⟨kl|\bar{G}(W=e_k+e_l)|kl⟩.$$

$$⟨k|U^{\mathrm{BHF}}|k⟩=\sum _l^A⟨kl|\bar{G}(W)|kl⟩,$$

$$⟨k^{\prime }|U^{\mathrm{BHF}}|k⟩=\{\begin{array}{ll}\frac{1}{2}{\sum }_{l=1}^A⟨k^{\prime }l|\bar{G}(e_{k^{\prime }}+e_l)+\bar{G}(e_k+e_l)|kl⟩,& e_{k^{\prime }},e_k\le e_F\\ {\sum }_{l=1}^A⟨k^{\prime }l|\bar{G}(e_k+e_l)|kl⟩,& e_{k^{\prime }}\le e_F,e_k>e_F\\ {\sum }_{l=1}^A⟨k^{\prime }l|\bar{G}(e^{\prime }+e_l)|kl⟩,& e_{k^{\prime }},e_k>e_F,\end{array}$$

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

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

$$⟨{\mathbf{q}}^{\prime }{\lambda }_1^{\prime }{\lambda }_2^{\prime }|V_s|\mathbf{q}{\lambda }_1{\lambda }_2⟩$$

$$⟨{\mathbf{q}}^{\prime }{\lambda }_1^{\prime }{\lambda }_2^{\prime }|V_{pv}|\mathbf{q}{\lambda }_1{\lambda }_2⟩$$

$$⟨{\mathbf{q}}^{\prime }{\lambda }_1^{\prime }{\lambda }_2^{\prime }|V_v|\mathbf{q}{\lambda }_1{\lambda }_2⟩$$

$$\times \left[g_v^2\bar{u}(-{\mathbf{q}}^{\prime },{\lambda }_2^{\prime }){\gamma }_{\mu }u(-\mathbf{q},{\lambda }_2)+\frac{f_v}{2M}\bar{u}(-{\mathbf{q}}^{\prime },{\lambda }_2^{\prime }){\sigma }_{\mu \nu }i(q^{\prime }-q{)}^{\nu }u(-\mathbf{q},{\lambda }_2)\right],$$

$$u(\mathbf{q},\lambda)=\sqrt{\frac{E_{\mathbf{q}}+M}{2M}}\left(\begin{array}[]{c}1\\ \frac{2\lambda|\mathbf{q}|}{E_{\mathbf{q}}+M}\end{array}\right)|\lambda\rangle,$$

$$u(\mathbf{q},\lambda)=\sqrt{\frac{E_{\mathbf{q}}+M}{2E}}\left(\begin{array}[]{c}1\\ \frac{2\lambda|\mathbf{q}|}{E_{\mathbf{q}}+M}\end{array}\right)|\lambda\rangle.$$

$$\frac{{\Lambda }_{\varphi }^2-m_{\varphi }^2}{{\Lambda }_{\varphi }^2+({\mathbf{q}}^{\prime }-\mathbf{q}{)}^2},$$

$$f_{ps}=g_{ps}\frac{m_{ps}}{2M}.$$

$$\mathcal{M}=\mathcal{V}+\mathcal{VGM},$$

$$\mathcal{M}=$$

$$\mathcal{W}=$$

$$\mathcal{M}=V+Vg\mathcal{M},$$

$$\mathcal{M}({\mathbf{q}}^{\prime },\mathbf{q}|\mathbf{P})=V({\mathbf{q}}^{\prime },\mathbf{q})+\int \frac{d^{3}k}{(2\pi {)}^3}V({\mathbf{q}}^{\prime },\mathbf{k})\frac{M^2}{E_{\mathbf{k}}{E}_{\mathbf{P}/2+\mathbf{k}}}\frac{1}{2E_{\mathbf{q}}-2E_{\mathbf{k}}+iϵ}\mathcal{M}(\mathbf{k},\mathbf{q}|\mathbf{P}),$$

$$\mathcal{H}=\sum _i\frac{\partial \mathcal{L}}{\partial ({\partial }^0{\varphi }_i)}{\partial }^0{\varphi }_i-\mathcal{L},H=\int d^{3}r\mathcal{H}(\mathbf{r}),$$

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

Full text

††journal: Progress in Particle and Nuclear Physics

1. 1 Introduction

2. 1.1 Brief introduction on nuclear theory

3. 1.2 The realistic nucleon-nucleon interaction

4. 1.3 Description of nucleus as a relativistic system

5. 1.4 Relativistic Brueckner-Hartree-Fock theory

6. 1.5 Connections between RBHF and CDFT

7. 1.6 Structure of the Review

8. 2 Basic Concepts

9. 2.1 Hohenberg-Kohn and Kohn-Sham density functional theory

10. 2.2 Challenges in nuclear DFT

11. 2.3 Current status in nuclear density functional

12. 2.4 Towards ab initio nuclear density functional

13. 3 Relativistic Brueckner-Hartree-Fock Theory and Applications

14. 3.1 BHF formalism

15. 3.2 RBHF formalism and its application for nuclear matter

16. 3.2.1 RBHF formalism

17. 3.2.2 Nuclear matter

18. 3.3 RBHF formalism for finite nuclei

19. 3.4 Local density approximation on RBHF and its promotion to CDFT

20. 3.5 Self-consistent RBHF calculations for finite nuclei

21. 3.5.1 Convergence check

22. 3.5.2 Center-of-mass motion

23. 3.5.3 Self-consistent basis

24. 3.5.4 Choice of the single-particle potential

25. 3.6 Ground-state properties for finite nuclei

26. 3.7 Pseudospin and spin symmetries

27. 4 RBHF Calculations for Neutron Drops

28. 4.1 Neutron drops

29. 4.2 General properties

30. 4.3 Spin-orbit splitting and tensor effects

31. 4.4 Connections to CDFT

32. 5 Summary and Perspectives

1 Introduction

1.1 Brief introduction on nuclear theory

The discoveries of radioactivity by Becquerel [1] and the Curies [2, 3] and the existence of a compact nucleus at the center of an atom by Rutherford et al. [4] opened the door of nuclear physics. During the hundred years of development in nuclear physics, there emerged several significant milestones, including the discovery of the neutron by Chadwick [5] which verified the composition of the nucleus as protons and neutrons, the meson-exchange theory for the strong interaction between nucleons by Yukawa [6], the independent-particle shell model of the nucleus by Goeppert-Mayer [7], Haxel, Jensen, and Suess [8], and the collective Hamiltonian for nuclear rotation and vibration by Rainwater [9], Bohr and Mottelson [10, 11], etc.

With the understanding of the composition of a nucleus as protons and neutrons [5] and the meson-exchange theory for the strong interaction between the nucleons [6], nuclear physicists hoped to describe the nucleus, a quantum many-body system, from the underlying nucleon-nucleon interaction. Euler, a student of Heisenberg, assumed the nuclear force as a two-body (2N) interaction with a Gaussian shape and calculated the infinite nuclear system, i.e., homogeneous nuclear matter, using second-order perturbation theory [12]. However, the strong repulsive core of the realistic nuclear force [13] prevents the application of perturbation theory.

On the other hand, the nuclear structure model with a phenomenological mean field achieved great success. Goeppert-Mayer [7], Haxel, Jensen, and Suess [8] introduced a strong spin-orbit potential and proposed the nuclear independent-particle shell model that successfully explained the conventional magic numbers in nuclei. Rainwater [9], Bohr and Mottelson [10, 11] explored the nuclear deformation and proposed the nuclear collective model that provided successful descriptions for nuclear collective excitations. Although the independent-particle shell model could describe the single-particle motion in a nucleus with a phenomenological mean-field potential, it could not provide even a qualitative description for nuclear bulk properties. On the contrary, a unified phenomenological description of nuclear vibrations and rotations could be achieved by the collective Hamiltonian, whereas it was helpless in understanding the motion of a single nucleon.

Since the 1970s, nuclear density functionals that start from effective nucleon-nucleon interactions (or Lagrangians) emerged, such as the Skyrme-Hartree-Fock theory [14], the relativistic mean-field (RMF) theory [15], and the Gogny-Hartree-Fock theory [16]. The parameters in these density functionals are determined by fitting to the properties of nuclear matter and the ground or even excited states of selected nuclei, and thus they are called the phenomenological density functionals or effective interactions. These effective interactions can describe the ground and excited states for almost all nuclei in the whole nuclear chart very well and became important tools to study heavy nuclei and exotic nuclei far away from the stability region [17, 18, 19].

The efforts to describe nuclear structure from realistic bare nucleon-nucleon interactions, later on known as nuclear ab initio calculations, have been started in the 1950s. Brueckner, Levinson, and Mahmoud [20] defined the reaction matrix in the nuclear medium, the so-called $G$-matrix, which took into account the two-body short-range correlations induced by the strong repulsive core of the nuclear force, and reproduced in the Hartree-Fock approximation qualitatively the saturation properties of nuclear matter. Different from that, Jastrow [21] introduced the correlation function to take into account correlation effects induced by the strong repulsive core, and to determine the nuclear ground state using a variational method. By starting from a Slater determinant, Coester [22] proposed the coupled-cluster method. Inspired by Brueckner’s $G$-matrix, Kuo and Brown [23] developed the configuration-interaction shell model. After the 1980s, with improvements of the nuclear force and with increasing computational resources, ab initio calculations starting from realistic nucleon-nucleon interactions have been largely promoted with more and more advanced many-body methods, such as the quantum Monte-Carlo method [24], the self-consistent Green’s function method [25], the no-core shell model [26], the Monte-Carlo shell model [27], the nuclear lattice effective field theory [28], or the in-medium similarity renormalization group [29]. For recent reviews, see Refs. [30, 31, 32, 33, 34, 35].

Recently, relativistic Brueckner-Hartree-Fock (BHF) theory, a fully self-consistent relativistic version of ab inito calculations has been successfully applied to finite nuclei [36, 37], and this will be reviewed in more detail in the following text. On the one hand, ab initio calculations are of fundamental significance by themselves. On the other hand, they can also be used to guide the development of nuclear density functionals, which was one of the original ideas of nuclear density functional theory [38]. See, in particular, Sections 1.3-1.5 for the motivations of RBHF and CDFT as well as the connections in between.

1.2 The realistic nucleon-nucleon interaction

The nuclear force models can be dated back to 1935 when Yukawa proposed the meson-exchange model [6]. The one-pion-exchange model achieved great successes in explaining the nucleon-nucleon scattering and deuteron properties. It provided a good description of the long-range part of the nuclear force. However, investigations of multi-pion-exchange processes in the 1950s encountered severe difficulties, and there existed considerable ambiguities in the theoretical models [39, 40, 41]. On the other hand, experiments on the nuclear force had made much progress at that period. By analyzing the proton-proton scattering data at $340$ MeV, Jastrow pointed out that there is a very strong repulsive core in the short-range part of the nuclear force and introduced the hard-core model [13].

In parallel, many efforts were made to construct phenomenological nuclear forces according to symmetries. In 1941, Eisenbud and Wigner [42] considered possible interaction forms which linearly depended on the relative momentum of the two nucleons, including the central, spin-spin, tensor, and spin-orbit terms. If a quadratic dependence of the momentum is allowed, there will be a quadratic spin-orbit term as well [43]. With such interaction forms, Gammel and Thaler obtained in 1957 a quantitative nuclear force model by fitting to the nucleon-nucleon scattering data for the first time [44]. In the Gammel-Thaler potential, a hard-core was used for the short-range part ($r\leq 0.4$ fm). During the 1960s, phenomenological nuclear forces were further developed. The Yale group [45] and Hamada and Johnston [46] independently developed their nuclear force models with a quadratic spin-orbit term. Reid not only improved further the hard-core potential but also proposed a soft-core potential [47]. Comparing with the hard-core potential, the short-range part of the nuclear force in the soft-core potential is no longer infinite repulsive but of finite values. One of the Reid soft-core potentials became one of the most widely used nuclear force later [47].

With more mesons discovered in the 1960s [48, 49], the nuclear force based on the meson-exchange picture was been revived again. The nuclear force based on the one-boson-exchange (OBE) model was proposed and achieved great success [50, 51]. It is still one of the most economical models in describing the nucleon-nucleon interaction quantitatively. It provides a simple parametrization and can describe the nucleon-nucleon scattering data well using only a few parameters.

However, the scalar isoscalar $\sigma$ boson, which provides the middle-range attraction, has a very large width ($400\sim 700$ MeV as summarized by the Particle Data Group [52]). It is considered to be related to the two-pion exchange, and many efforts were made to derive the middle-range attraction from the correlated two-pion exchange, including the one-boson-exchange model of the Stony Brook group [53], the Paris potential [54] derived from dispersion relations, the Partovi-Lomon model [55] and the Bonn potential [56] from field theory.

With the establishment of quantum chromodynamics (QCD) as the basic theory of strong interaction, the meson-exchange theory was no longer the basic starting point for the nuclear force. However, because of its non-perturbative nature at low energies, it is difficult to derive the nuclear force directly from QCD. The quark model which started from QCD began to be developed [57, 58]. It tried to describe the hadron structure and the hadron-hadron interaction in a unified way from the quark level, and many features of the nuclear force were explained successfully. Nevertheless, the quark model is still a phenomenological model with a set of parameters and therefore it cannot be viewed as a basic theory either.

One breakthrough appeared in 1990 when Weinberg applied the chiral effective field theory [59] to the low-energy region of QCD [60, 61]. The idea is that the effective degrees of freedom of the nuclear system in the low-energy region should be nucleons and pions instead of quarks and gluons. The pions are, as Goldstone bosons, connected with the violation of chiral symmetry in QCD. Then the nuclear force can be described by a Lagrangian of nucleons and pions satisfying the chiral symmetry of QCD. The chiral Lagrangian provides a perturbative expansion framework in deriving the nucleon-nucleon interaction, where the small quantity in the expansion is $Q/\Lambda_{\chi}$ with the low momentum $Q$ and the scale of chiral symmetry breaking $\Lambda_{\chi}\sim 1$ GeV. This low-energy perturbation theory is also called chiral perturbation theory, with the potential being iterated at least to the leading order. The perturbative order is characterized by the power of $Q$, which defines a power counting scheme [61]. Ordóñez and van Kolck [62] constructed the next-to-next-to-the-leading order (NNLO) chiral force in time-ordered perturbation theory reproducing most features of the nuclear force in a quantitative way. Using the unitary transformation method, Epelbaum, Glöckle, and Meissner (Bochum-Jülich group) solved the energy-dependent problem in time-ordered perturbation theory and constructed the NNLO chiral force [63, 64]. Afterward, the three-body force from NNLO was introduced [65]. Entem and Machleidt (Idaho group) [66] and the Bochum-Jülich group [67] constructed the chiral force up to N3LO, and then up to N4LO [68, 69]. By investigating the major terms in N5LO, it was found that the expansion of the chiral force achieved good convergence at N4LO [70]. Recently, a relativistic chiral force was developed, and here the description for the nucleon-nucleon scattering data at the leading order is comparable with that of the non-relativistic NLO chiral force [71, 72]. This indicates that the expansion of the chiral force may achieve faster convergence in the relativistic framework. For recent reviews on the nuclear chiral force, see Refs. [73, 74].

Since the 1990s, several versions of high-precision nuclear forces have been constructed. With the effort of the Nijmegen group, the analysis of nucleon-nucleon scattering phase shifts has been improved notably [75]. They made a multi-energy partial-wave analysis for all nucleon-nucleon scattering data below $T_{\rm lab}=350$ MeV, published between 1955 and 1992. The final database consists of $1787$ proton-proton and $2514$ neutron-proton scattering data after careful examination. The fitting result is $\chi^{2}/N_{\rm df}=1.08$, with the total number of degrees of freedom $N_{\rm df}=3945$. This result is remarkably better than the preceding phase-shift fitting accuracy [75] and forms the ground for a new generation of high-precision nuclear forces. Other well known high-precision nuclear forces include Nijmegen [76], Reid93 [76], and Argonne V18 [77] potentials based on phenomenological models, or CD-Bonn [78] based on the one-boson-exchange model. For the forces based on chiral perturbation theory, for example, the $\chi^{2}$/datum for the description of the Nijmegen phase shifts has achieved 0.3 (np) and 0.6 (pp) for data with an energy less than 200 MeV calculated by N4LO (Bochum group) [68], and it has achieved 1.15 for data with an energy less than 290 MeV by N4LO (Idaho group) [79]. These high-precision nuclear forces can reproduce experimental nucleon-nucleon scattering data accurately and provide the foundation for nuclear many-body calculations.

The lattice QCD approach solving the QCD on a lattice provides a promising way to derive the nuclear force directly from QCD. There are two major approaches to study hadron-hadron interactions with lattice QCD. One is based on Lüscher’s finite volume formula [80], which extracts the ground-state energy using the temporal correlation functions on a finite lattice space, see Ref. [81] for a review. Another one is the HAL QCD method [82], which defines and extracts the interaction between hadrons through the Nambu-Bethe-Salpeter wave function, see Ref. [83] for a review. Recent lattice simulations, extensively developed by the HAL QCD collaboration, have obtained successfully the baryon-baryon interactions with a realistic pion mass [84]. This provides a new possibility to describe the nuclear system from QCD.

1.3 Description of nucleus as a relativistic system

The relativistic description for nuclear systems goes in parallel with the non-relativistic one. In 1955, Johnson and Teller found that a good description of nuclear matter saturation properties can be obtained if a velocity-dependent nuclear force is introduced [85]. Duerr reformulated this theory in a covariant form, avoided the collapse of the nucleus occurring in the non-relativistic theory for high kinetic energies, obtained similar saturation properties, and predicted naturally a large spin-orbit splitting [86]. Rozsnyai extended Duerr’s Hamiltonian to study finite nuclei [87]. By readjusting the coupling strengths and the range of the force, he obtained a reasonable binding energy and density distribution. After the discovery of scalar and vector mesons, the one-boson-exchange potential was gradually established [50]. Miller and Green developed the relativistic Hartree-Fock (RHF) theory based on a model of scalar and vector mesons exchange, which produces strong attractive and repulsive interactions, respectively [88].

In 1974, Walecka proposed a renormalizable relativistic mean-field $\sigma$-$\omega$ model that gave good saturation properties of nuclear matter [15]. In the calculation of the densities, only the particles in the Fermi sea are considered. Within the no-sea approximation [89], the contribution from the Dirac sea is neglected, because they would lead to couplings to high energy nucleon-antinucleon pairs. From the point of view of effective field theory [60], it can be understood that the high-energy intermediate states, including nucleon-antinucleon pairs, can be integrated out using an appropriate renormalization scheme [90]. In practice, all successful relativistic calculations based on this model have used the no-sea approximation.

The success of the simple Walecka model promoted the development of relativistic many-body theory. Boguta and Bodmer [91] solved the problem of too large incompressibility in Walecka’s model by introducing nonlinear self-couplings of the $\sigma$ meson. Serot extended Walecka’s $\sigma$-$\omega$ model by including the isovector $\pi$ and $\rho$ mesons [92]. In nuclei with open valence shells the Munich group presented a successful relativistic description of axial [93] and triaxial [94] deformations, pairing correlations have been introduced in the framework of relativistic Hartree-Bogolubov (RHB) theory [95] and rotating nuclei are well described by relativistic cranking theory [96, 97]. For the development of relativistic nuclear many-body theory in this period, one can see the review papers [89, 98, 99, 100]. In these RMF models, the coupling strengths were determined by fitting to the saturation properties of nuclear matter and to the selected data of finite nuclei, similar to the zero-range Skyrme force [14] or the finite-range Gogny force [16] in the non-relativistic framework.

The phenomenological relativistic mean-field theory [15, 101], also called covariant density functional theory (CDFT) in recent years, is able to reproduce very successfully basic properties for most of the nuclei in the nuclear chart. Along the branch with the nonlinear self-couplings of the meson fields, also referred as NLRMF, the parameter set NL1 [102] gives a good description of stable nuclei, the sets NL-SH [103], NL3 [104], and NL3* [105] extended this for nuclei far from stablity, TM1 [106] further considered the nonlinear self-couplings of $\omega$ mesons, and PK1 [107] takes into account a microscopic center-of-mass (c.m.) correction and neutron-proton mass differences. The idea of nonlinear meson couplings goes back to the 1970s, a time when the relativistic models were considered as a first step to a fully fledged quantum field theory [89] for the description of nuclei, called Quantum Hadrodynamics (QHD) [108]. For a review see [109]. Such a theory had to be renormalizable [90] and therefore Boguta and Bodmer [91] introduced a $\phi^{4}$-coupling scheme in the meson sector. When it turned out that there were serious problems with renormalization in the second order [110], the RMF theory was interpreted as a relativistic density functional theory based on the ab initio model of Brockmann and Toki [111] with density-dependent coupling strengths adjusted to relativistic Brueckner-Hartree-Fock calculations. This model was too simple to compete with the very successful phenomenological nonlinear meson-coupling models, but a phenomenological density dependence (referred as DDRMF) with several parameters, adjusted in a similar way to experimental data as the nonlinear models, turned out to be very successful. Moreover, DDRMF is closer to the spirit of density functional theory. Therefore, it become more and more popular with extensive applications and developments of effective Lagrangians, such as TW99 [112], PKDD [107], DD-ME1 [113], DD-ME2 [114], etc.

In the RMF Lagrangians mentioned above, the nuclear force is propagated by massive mesons, thus keeping the feature of finite range. This needs some numerical efforts, in particular for excited states in the density dependent relativistic Random Phase Approximation and for methods going beyond mean field approach. Therefore, similar to the zero-range Skyrme force in the non-relativistic case, the relativistic point-coupling model, a new branch of the RMF approach, was developed by replacing the meson exchange with the local many-body contact interactions. Compared to NLRMF and DDRMF, the relativistic point-coupling model avoids solving the complicated meson fields and therefore has the advantages of a simple form, fast computing, and easy to be extended to methods beyond mean field approach. Commonly used point-coupling RMF sets include the nonlinear PC-LA [115], PC-F1 [116], PC-PK1 [117], and the density-dependent DD-PC1 [118].

With these effective interactions, the relativistic many-body theory has achieved great success in describing properties of nuclear ground-states and excited-states. The related reviews include, for example, the description of nuclear bulk properties, exotic nuclei, pairing correlations, high-spin states, and time-dependent problems [101, 18], relativistic continuum Hartree-Bogoliubov theory for exotic nuclei [119], relativistic quasiparticle random phase approximation for collective modes in exotic nuclei [120], relativistic density functional theory beyond mean field [121, 122], the study of nuclear magnetic and antimagnetic rotation with the tilted axis cranking covariant density functional theory [123], the relativistic description of pseudospin and spin symmetries [124], the investigation of halo phenomena in medium-heavy and heavy nuclei [125]. The current status in the relativistic many-body theory for nuclear structure has been summarized in Ref. [19].

Although the relativistic many-body theory achieved great success in describing many properties of nuclear structure, it is based on phenomenological interactions, which are determined by fitting to ground-state properties of stable nuclei. Substantial ambiguities exist in the prediction for unstable nuclei without data and for nuclear matter at high densities. On the other hand, some effects which do not show up clearly in stable nuclei, such as the tensor-force effect, are difficult to study by phenomenological interactions. Meanwhile, the tensor force is very important in describing the properties of neutron-rich exotic nuclei, in particular, the evolution of the shell-structure [126]. Particularly in systems with non-spin-saturated configurations it has a strong influence on the single-particle structure.

One has to distinguish between the tensor part of the bare nucleon-nucleon interaction, which is in principle determined by the scattering data, and the tensor term in the effective interaction in the nuclear medium. In microscopic nuclear matter calculations, the bare tensor force plays a crucial role. The second-order diagram containing a two-pion exchange coupled to $T=0$ is the origin of the major part of the middle-range attraction, which is, in phenomenological models, usually described by the exchange of the isoscalar $\sigma$ meson. Nuclear matter is spin-saturated. There is no spin-orbit splitting, and thus the so-called “first-order” tensor in the effective force, which acts between particles in open spin configurations, can be neglected. This is very different in finite heavy nuclei with open spin configurations, where the two spin-orbit partners are on different sides of the Fermi surface. Here the tensor force can lead to considerable shifts in the single-particle energies, and these in turn can lead to new magic numbers, to changes of the deformation, and to quantum phase transitions. In most of the successful relativistic and non-relativistic density functionals, effective tensor terms are not taken into account. The phenomenological fit of such functionals to bulk properties includes such effects in a somewhat averaged way, and as a consequence these functionals are relatively successful in the region of stable nuclei. Systematic investigations of relativistic density functionals [127, 128] have shown that the various types of functionals produce rather similar results in the region of stable nuclei, where their parameters have been adjusted. There are, however, essential differences far from stability, and they are often caused by differences in the single-particle energies in those regions [129]. There is hope that a proper treatment of the tensor force [130] can help to increase the predictive power of such future functionals.

In such cases, ab initio calculations become very important. By solving the nuclear system on the basis of realistic bare nucleon-nucleon forces, one can establish the link between the microscopic nucleon-nucleon interaction and nuclear structure properties. This will provide valuable information in describing exotic nuclei, nuclear matter under extreme conditions, and the effects of the tensor-force. In this way one will be able improve the nuclear density functional [131].

1.4 Relativistic Brueckner-Hartree-Fock theory

Brueckner theory [20] has been introduced in the 1950s to overcome the difficulties caused by the extremely strong short-range repulsion of the bare nuclear force [13] in many-body calculations. For the description of nucleon-nucleon scattering, the scattering amplitude calculated in first-order, i.e., in Born approximation, is totally wrong. The correct scattering amplitude can be obtained [132] only by the exact solution of the two-particle Lippmann-Schwinger equation [133]. In multiple scattering of hadrons on the nucleus, instead of using the realistic nuclear force, Watson et al. [134] constructed an optical potential via the scattering matrix ($T$ matrix) derived from the solution of the Lippmann-Schwinger equation and found rather good descriptions of the experimental results. By replacing the bare nuclear force with an effective interaction in the medium, later known as the $G$-matrix, Brueckner et al. [20] extended this idea to describe nuclear structure in a mean-field approximation. By summing all the ladder diagrams and taking into account the Pauli principle in nuclear medium, the obtained $G$-matrix incorporates the two-body short-range correlation and describes the saturation properties of nuclear matter within the Hartree-Fock approximation qualitatively.

Bethe elaborated the Brueckner theory by constructing a model of the nucleus, in which each nucleon moves in a self-consistent potential, which is obtained from the reaction matrix for two nucleons derived from scattering theory in the nuclear medium [135]. This provides the foundation for the extension of Brueckner theory from nuclear matter to realistic finite nuclei. A formal derivation of Brueckner theory was provided by Goldstone [136]. These are breakthroughs in microscopic nuclear many-body theory, and many developments have been made along this direction in the 1960s [132, 137, 138].

In Brueckner theory, the $G$-matrix is obtained by solving the scattering equation in the nuclear medium. In contrast to two-body scattering in free space, the Pauli principle must be taken into account, because all the states below the Fermi surface are occupied in the nuclear medium and therefore they cannot be used in the intermediate scattering processes. This effect of the Pauli principle was already pointed out by Brueckner very early [20], but, at the beginning, it was neglected as an approximation. Later, Brueckner and Wada [139] took this effect into account in perturbation theory. A full treatment of the Pauli principle was first given by Bethe and Goldstone [140], and therefore, the scattering equation in Brueckner theory is also called the Bethe-Goldstone (BG) equation.

Because of the Pauli operator and the energy denominator, the Bethe-Goldstone equation represents a nonlinear integral equation, which was very difficult to solve even for nuclear matter. By replacing the spectrum of the intermediate states with a ”reference spectrum”, Bethe, Brandow, and Petschek reduced the Brueckner integral equation to a differential one which can be solved easily [141]. Afterwards they took the Pauli operator into account order by order in a perturbative method. Separating the two-nucleon interaction into short-range and long-range parts, Moszkowski and Scott [142] proposed a method, which is very useful in solving the diagonal matrix elements of the $G$-matrix but it turns out to be relatively difficult for the non-diagonal parts. These methods are mainly applicable to the static nuclear force. For a velocity dependence, one usually needs to work in momentum space. Brown, Jackson, and Kuo [143] as well as Haftel and Tabakin [144] independently developed the matrix inversion method to solve the Bethe-Goldstone equation in momentum space. It accurately reproduces the results obtained in coordinate space.

Compared with the bare nuclear force, the $G$-matrix has included some two-body short-range correlations and can be viewed as an effective nuclear interaction suitable to be used in a perturbation expansion. For a hard-core potential, the matrix elements of the bare interaction are divergent, but the elements of the $G$-matrix are finite [145]. This indicates that one can use the $G$-matrix for a perturbative expansion. The result is called the Brueckner-Bethe-Goldstone expansion, or the hole-line expansion as the expansion order is indicated by the number of independent hole lines in each diagram [137]. The convergence of this expansion is roughly determined by the particle density times the wound integral of the two-body wave function [132]. The lowest order gives the Brueckner-Hartree-Fock framework. In the BHF or general Hartree-Fock framework, a good choice of the single-particle mean-field potential can cancel contributions of certain higher order diagrams. This is good for a rapid convergence of the expansion. For a review of the three-hole-line expansion beyond BHF theory, see Ref. [137], and for applications for finite nuclei, see Ref. [138].

Early theoretical works showed that BHF theory was a good approximation for the nuclear matter system near or below the saturation density. However, from a quantitative point of view, the saturation properties given by BHF were not good enough compared with the empirical data. Coester et al. [146] found that all the BHF results with various nuclear forces were located on or near the now-called Coester line, which deviated from the empirical saturation region systematically. Neither can BHF theory give a good description of binding energies and charge radii of finite nuclei at the same time. On the other hand, in the late 1970s, a totally different approach, the variational method based on a correlation operator acting on a Slater determinant was developed to deal with the short-range repulsion of the nuclear force, and the results obtained with this method were similar as those of the Brueckner method [147]. This leaded to the discussion on the missing of a three-body (3N) force, which has been proposed at earlier times [148, 149, 150, 151]. By introducing a three-body force, good descriptions of nuclear matter saturation properties as well as the ground state and a few excited states properties of light nuclei were obtained [152]. Including a three-body force, BHF has been extended to study neutron star and nuclear matter properties [153, 154, 155].

Inspired by the success of the relativistic phenomenological many-body theory in the 1970s, different research groups tried to extend Brueckner theory to the relativistic framework. In the pioneering work of the Brooklyn group, the single-particle wave function was chosen as the Dirac spinor in free space, and the relativistic effect was taken into account in first-order perturbation theory [156]. The results obtained in this way showed a significant improvement compared to non-relativistic investigations. Both binding energy and saturation density were much closer to the empirical values. Horowitz and Serot [157, 158], Brockmann and Machleidt [159, 160], and ter Haar and Malfliet [161, 162] independently solved the nuclear matter system using relativistic Brueckner-Hartree-Fock (RBHF) theory self-consistently, see also Refs. [163, 164, 98, 162]. With the one-boson-exchange potentials Bonn A, B, and C [165], the RBHF results improved the non-relativistic Coester line greatly [160]. RBHF theory was also extended to study asymmetric nuclear matter [166], the optical potential [167], neutron stars [168], neutron stars with hyperon degree of freedom [169], etc. Latest reviews can be found in Refs. [170, 171, 172, 173, 174].

The fact, that relativistic effects play such an important role, is somehow beyond expectations, as the Fermi momentum of a nuclear system near saturation density is small compared with the rest mass of the nucleon. It raises questions of where this relativistic effect comes from and why it is essential in improving the theoretical description. In the relativistic framework, the exchange of scalar mesons induces a strong attraction between nucleons. For the particle states with positive energy, it is to a large extent compensated by the repulsion caused by the strong zeroth component of the vector field induced by the omega exchange. For the negative-energy states, the effects of both fields add up. This has two essential consequences: (i) A very strong velocity-dependent part of the field in the form of a strong spin-orbit term. It cannot be neglected even for the relatively small velocities in the Fermi sea, because the factor in front of the velocity is very large. (ii) An additional repulsion between the nucleons increasing with the density. The scalar field is proportional to the Dirac effective mass of the nucleon. This effective mass decreases as the density increases. Therefore, the average attraction felt by the nucleons decreases too. On the other hand, the vector mesons producing repulsion have no such feature. The net effect is that the repulsion felt by nucleons increases with increasing density.

This density-dependent effect leads to the saturation of nuclear matter in the relativistic framework [175, 176, 177, 160]. In non-relativistic density functionals such as Skyrme [14] or Gogny [16], this effect does not exist and one needs a strongly repulsive phenomenological density-dependent term $t_{3}$ to find proper saturation. By expanding the Dirac spinor for plane waves with an effective mass $M^{*}/M$ in nuclear matter in terms of the free Dirac spinor with effective mass $M^{*}/M=1$, Brown et al. [176] pointed out that, in a relativistic description, one effectively takes into account nucleon-antinucleon excitations (or the so-called $Z$-diagram). Zuo et al. [155] showed in a non-relativistic BHF investigation of the possible contributions to three-body forces in nuclear matter that this $Z$-diagram contributes a large part to the total three-body forces, at least for densities not much higher than the saturation density. It has also been shown in nuclear matter that the results of RBHF with 2N interactions is very similar to BHF with chiral 2N + 3N interactions [178]. Furthermore, a good description of both nuclear matter and finite nuclei has been achieved in a non-relativistic framework with a 3N interaction that simulates the relativistic effects  [173]. This explains clearly the reason why, without three-body forces, the relativistic BHF calculations lead to a much better description of the saturation properties of nuclear matter than the non-relativistic ones.

Although RBHF theory with only two-body forces can already give rather good saturation properties of nuclear matter, there are still several open problems [170]. Because of the cluster effects, the computation at low density becomes unstable [179]. Since the solution of the Bethe-Goldstone equation in most of the RBHF calculations of nuclear matter do not include the antinucleon degrees of freedom, the scalar and vector channels of the nucleon self-energy cannot be decomposed uniquely  [180]. Such a decomposition is necessary for the solution of the RHF equations in the next step of the iteration. Different decompositions produce different results. The method proposed by Brockmann and Machleidt [160] gives reasonable results for symmetric nuclear matter, but not for asymmetric nuclear matter [181, 182]. Another possibility is the mapping of the $G$-matrix onto 5 (or 6 for asymmetric nuclear matter) complete Lorentz invariants [158], but such a mapping is not unique either because of the missing antinucleon degrees of freedom. Nowadays, a commonly used method introduced by Müther and his collaborators [181, 183] is to decompose the $G$-matrix in the form $G=V+\Delta G$ where $V$ is the bare nuclear force. This is similar to the optimal $T$-matrix representation developed by the Tübingen group [184, 185, 186]. Boersma et al. [187] parameterized the $G$-matrix in the form of a Yukawa function with form factors by fitting the masses and the coupling constants. De Jong and Lenske [188] and Huber et al. [189] tried to solve this problem by projecting in a full Dirac space where both the positive-energy and negative-energy solutions are included. However, because the bare nucleon-nucleon interaction was obtained by considering the positive-energy states only, such a full-space projection is not consistent [188]. More discussions on RBHF theory in nuclear matter can be found in the review [170].

With the success of RBHF theory in nuclear matter, it is natural to extend RBHF for finite nuclei. However, due to considerable numerical difficulties, such a project was extremely challenging and a full solution was for a long time not possible. Müther et al. [190, 191] started with a non-relativistic BHF code for finite nuclei [192] and took relativistic effects into account using the effective density approximation. More investigations were mainly based on the local density approximation (LDA) [38]. The idea of LDA is to find an effective interaction which can reproduce the RBHF results for nuclear matter at a fixed density. This lead to an effective density-dependent interaction which is used to study finite nuclei. However, because of the uncertainties of RBHF theory in nuclear matter mentioned above, the mapping from nuclear matter to finite nuclei is far from unique [111, 193, 194, 187, 195, 196, 181, 182, 197, 198, 199].

1.5 Connections between RBHF and CDFT

Based on the meson exchange picture of the nuclear force, the RMF approach, consistent with the spirit of density functional theory [200, 201], has achieved many successes in describing properties of nuclear ground and excited states. In the RMF approach, i.e. in the relativistic Hartree approach, some important ingredients are missing, such as the $\pi$- and $\rho$-tensor couplings which enter via the exchange (Fock) terms. Almost in parallel to the RMF approach, the relativistic Hartree-Fock (RHF) theory was initialized by Miller and Green [88] and applied systematically to nuclear matter [202] and finite nuclei [203]. Similar to the RMF approach, the problem of the too large incompressibility was solved by nonlinear self-couplings of the $\sigma$ meson [204] or, in the zero-range limit, by three- and four-body contact terms [205]. However, due to the numerical complexity of the Fock terms, it remained a long-standing problem for RHF theory to provide a quantitative description of nuclear properties comparable with RMF.

The situation changed when computer power was increasing and when an explicit density-dependence of the meson-nucleon coupling strengths was introduced into the RHF scheme, leading to the so-called density-dependent relativistic Hartree-Fock (DDRHF) theory [206, 207, 208]. In this case it became possible, at least in the spherical case, to carry out extensive calculations, such that reasonable new RHF Lagrangians have been obtained by the Peking-Orsay group, such as PKO1 [206], PKO2 [209], and PKO3 [209], and by the Peking-Aizu group, such as PKA1 [207]. Within these models, similar accuracy as in standard RMF has been achieved for a quantitative description of nuclear properties. Pairing correlations can been included in relativistic Hartree-Fock-Bogoliubov (RHFB) theory [95]. The corresponding equation can now be solved also for finite spherical nuclei [208].

The explicit treatment of the Fock terms brings significant improvements. The missing degrees of freedom associated are taken into account naturally with the $\pi$ and $\rho$-tensor couplings via the Fock terms, which improve the self-consistent description of nuclear shell evolution [209, 210, 211]. The inclusion of $\rho$-tensor couplings in PKA1, which corresponds to a strong attractive potential, leads to a better preserved pseudo-spin symmetry [212, 207, 213], and eliminates the spurious shell closures $N,Z=58$ and $92$ that commonly exist in the conventional covariant density functional calculations [214, 207, 210]. With the presence of Fock terms, the nuclear tensor force, an important ingredient of the nuclear force, can be naturally taken into account [209, 215, 216, 217, 218] in a phenomenological way. Substantial effects due to the Fock terms were also revealed in extending the RHF scheme to describe exotic nuclei [208, 213, 219], superheavy magic shells [220], the nuclear symmetry energy, and the equation of state [221, 222, 223].

Besides the nuclear ground-state properties, the Fock terms also lead to essential improvements in the investigations of spin- and isospin excitations in nuclei. Combined with the random phase approximation (RPA), the RHF+RPA approach provides a fully self-consistent description for nuclear Gamow-Teller and spin-dipole excitations. Here the Fock terms of the isoscalar $\sigma$ and $\omega$ meson fields were found to play an essential role [224, 225]. Based on RHFB theory, a fully self-consistent quasiparticle random phase approximation (QRPA) was developed [226, 227] for the study the $\beta$-decay half-lives of spherical neutron-rich nuclei. The available data are well reproduced by including an isospin-dependent proton-neutron pairing interaction in the isoscalar channel of the RHFB+QRPA approach [226]. On the other side, with all these successes of RHF theory including Fock terms, one has to emphasize that the solution of the RHF equations with finite meson masses in deformed systems is still connected with enormous numerical difficulties even on the most modern computers. There exist codes for axially deformed nuclei working in an oscillator basis [228], but the number of oscillator shells is limited and therefore such calculations have only been carried out for relatively light nuclei. The RHF applications for triaxial or rotating nuclei and for investigations beyond mean field are not yet feasible at the moment.

As discussed in Sections 1.3 and 1.4, the success of the early CDFT investigations in the 1970s that take into account Lorentz invariance for the nuclear system stimulated relativistic Brueckner-Hartree-Fock theory for nuclear matter which has achieved great success. However, for a long period, RBHF studies for finite nuclei have stayed on the level of the local density approximation. This situation has changed recently. In Ref. [36], the RBHF equations were solved directly for finite nuclei in the Dirac Woods-Saxon (DWS) basis using the Bonn A interaction [165]. The DWS basis is obtained by solving the Dirac equation with fixed scalar and vector potentials of Woods-Saxon type [229]. It keeps the full relativistic structure of the Dirac spinors. Therefoe the solution of the BG equation in this basis does not need approximations like angle averaging [230, 231]. Taking 16O as an example, convergent results are achieved for an energy cutoff of around $1.1$ GeV. The resulting binding energies, charge radii, and spin-orbit splittings are considerably improved in comparison with non-relativistic BHF results. In Ref. [37], the fully self-consistent RBHF theory for finite nuclei was realized using a self-consistent RHF basis. In particular, the Pauli operator in the BG equation is defined in this basis. A detailed discussion of the RBHF formalism and its application for finite nuclei is given in Ref. [37]. Spin symmetry in the Dirac sea (or equivalently pseudospin symmetry in the Fermi sea) was discussed in Ref. [232].

It is noticeable that the fully self-consistent RBHF theory, without any adjustable parameters, has achieved such good results even at the Hartree-Fock level, using only the two-body interaction fitted to the nucleon-nucleon scattering and deuteron properties. Of course, the accuracy of RBHF calculations (as any other ab initio methods) cannot be compared with that of phenomenological CDFTs. On the other hand RBHF calculations can provide essential information to guide the construction of DFTs. For example, the difficulty of nuclear DFTs to fix some terms of the functional by experimental data, such as the tensor term [130], can be solved by nuclear ab initio calculations. As a simple ideal system, neutron drops have been studied using RBHF theory and a specific pattern due to the tensor force was found in the evolution of spin-orbit splittings as a function of the neutron number [233, 234]. This forms an important guide for determining the tensor force and microscopic derivations of relativistic and non-relativistic DFTs. As an example, recently a new Skyrme functional with a tensor force, SAMi-T, has been developed along this line [235].

1.6 Structure of the Review

In this Review, we will focus on the progress of the relativistic Brueckner-Hartree-Fock theory for finite nuclei, and its guidance towards an ab initio covariant density functional theory. The paper is organized as follows. The basic concepts connecting relativistic Brueckner-Hartree-Fock theory and covariant density functional theory will be given in Section 2. The relativistic Brueckner-Hartree-Fock calculations for finite nuclei and neutron drops will be discussed in Sections 3 and 4. Finally, a summary and perspectives will be given in Section 5.

2 Basic Concepts

2.1 Hohenberg-Kohn and Kohn-Sham density functional theory

Density functional theory is one of the most important and successful theories in dealing with quantum many-body systems [200, 201]. It was developed first in atomic systems with the Hohenberg-Kohn (HK) theorem [200]: for an electromechanical many-body system in an external potential $U_{\rm ext}(\mathbf{r})$, there exists a density functional $E_{U}[\rho]$ of the density $\rho(\mathbf{r})$,

[TABLE]

where the functional $F_{\rm HK}[\rho]$ does not depend on the external potential $U_{\rm ext}(\mathbf{r})$ nor on the particle number of the system (the particle-number dependence is included implicitly in the normalization of the density). It only depends on the interaction between the particles. It is, therefore, identical for all Coulombic systems, such as atoms, molecules, and solids. In this sense, this density functional is supposed to be universal. The variational principle determines the exact ground-state density $\rho_{\rm gs}(\mathbf{r})$ and the corresponding exact ground-state energy by minimizing $E_{U}[\rho]$ with respect to the density $\rho(\mathbf{r})$.

Although this theorem is exact, it does not show explicit clues on how to construct this functional $F_{\rm HK}[\rho]$. Furthermore, it is highly non-trivial, if not impossible, to derive shell effects from the variation of a density functional $E_{U}[\rho]$ with respect to the local density $\rho(\mathbf{r})$. This problem of HK theory has been solved by Kohn and Sham (KS) [201] in a very elegant scheme, by mapping the exact local density $\rho(\mathbf{r})$ in a unique way to a local auxiliary single-particle potential $v_{\rm KS}(\mathbf{r})$ defined in such a way that the exact ground-state density $\rho_{\rm gs}(\mathbf{r})$ of the interacting system is the same as the ground-state density of an auxiliary non-interacting system

[TABLE]

which is expressed in terms of the $N$ lowest occupied single-particle orbitals of the Kohn-Sham equations

[TABLE]

To date, most practical applications of density functional theory use this Kohn-Sham scheme.

The exact density functional $F_{\rm HF}[\rho]$ can be then decomposed into three separate terms:

[TABLE]

Here $T_{s}$ is the kinetic energy of the auxiliary non-interacting system, which is supposed to capture the main part of the kinetic energy of the realistic interacting system as well as the shell effects; $E_{H}$ is the Hartree energy, and all the rest is left in the exchange-correlation energy $E_{xc}$ which contains, by definition, the Fock term as well as all the other many-body effects. It is clear that the robustness of the Kohn-Sham scheme crucially depends on how accurately the universal exchange-correlation energy functional can be built.

The extension of the Hohenberg-Kohn theorem [200] to relativistic systems was first formulated in Ref. [236] by utilizing a quantum electrodynamics (QED)-based Hamiltonian with a four-current. Within the no-sea approximation, which is also called no-pair approximation in atomic physics, i.e. neglecting vacuum polarization, i.e. all effects due to the creation of particle-antiparticle pairs, the total energy of the system can be expressed as a functional with respect to the four-current $j^{\mu}(\mathbf{r})=(\rho(\mathbf{r}),\mathbf{j}(\mathbf{r}))$ [237, 238] instead of only $\rho(\mathbf{r})$ in the non-relativistic case. In the early development, however, questions related to zero-point energies, radiative corrections, and UV-divergences were not discussed. The UV-divergences showed up in the derivation of gradient corrections to the kinetic energy [239]. Detailed discussions of the renormalization and related issues can be found in Ref. [240].

At the local density approximation level, it is found that the relativistic LDA gives bigger errors than the non-relativistic counterpart. This comes from the fact that, in the non-relativistic LDA, there is a cancellation of errors between the exchange and the correlation energy, which does not hold in the relativistic framework [241]. In Ref. [242], the optimized-potential-method has been extended to the relativistic framework on the longitudinal no-pair level. Later on, this method was applied to the study of atoms, it was found that the relativistic optimized-potential-method agrees well with the relativistic Hartree-Fock calculation and relativistic effects become drastic for high-$Z$ atoms [243]. Though it achieved an accurate description, the optimized-potential-method is limited by its computational cost, because, besides the kinetic energy, the exchange energy is also treated in orbital representations. Later, the relativistic generalized gradient approximation has been developed, and a good description could be achieved with rather small computational cost [241]. For more recent developments and investigations of relativistic density functional theory in atomic systems, see Ref. [244].

2.2 Challenges in nuclear DFT

DFT in nuclear physics was introduced in the early 1970s in a slightly different way. First of all, there is an essential difference between the DFTs of Coulombic and nuclear systems.. The Hohenberg-Kohn theorem was derived under the assumption that there is an external potential $U_{\rm ext}$, which, in atoms, is the Coulomb potential provided by the nucleus. However, the nucleus itself is a self-bound system, and such an external potential does not exist . It is still not clear whether the Hohenberg-Kohn theorem holds in nuclei as an exact statement, or only as a reasonable approximation [245]. Recent discussions on the Hohenberg-Kohn theorem for a self-bound wave-packet state can be found in the reviews [246, 247]. On the other hand, the shell structure in nuclei pointed out by Goeppert-Mayer and Jensen [7, 8] indicates that there should be an effective mean field, in which eventually the Hohenberg-Kohn theorem and the Kohn-Sham scheme could be applied at least in a very good approximation.

The second challenge comes from the fact that the microscopic origin of such an effective mean field is not apparent at all. Because of the existence of a strong repulsive core in the nucleon-nucleon interaction, or in general because of the non-perturbative nature of this interaction, the attempt of deriving the mean field using Hartree-Fock theory or using the decomposition (4) into Hartree and exchange-correlation terms fails in nuclear physics, whereas it is very successful in Coulombic systems.

It was only the work of Brueckner et al. [20] that gave a better understanding, why the concept of a mean field works so well in nuclei. They treated the scattering process between two nucleons in the nuclear medium similar to the case in the vacuum summing up the ladder diagrams in the Born series. Furthermore, they took into account the Pauli principle, that two nucleons in the nuclear medium can only be scattered to the unoccupied states above the Fermi levels. The obtained scattering matrix is the so-called $G$-matrix, which is treated as an effective interaction between two nucleons in the nuclear medium and gives a reasonable description for nuclear systems even at the Hartree-Fock level. This is the BHF theory mentioned in the introduction.

BHF calculations in nuclear matter show, that the $G$-matrix is strongly density dependent and that Hartree-Fock calculations with this effective interactions failed [38]. Based on this observation, Vautherin and Brink [14] adopted a simple phenomenological interaction introduced by Skyrme [248] containing a zero-range three-body term. At the mean-field level, this leads to an effective density-dependent two-body interaction. By an appropriate fit, they were able to reproduce not only the experimental binding energies but also the radii at the same time. This illustrated the importance of density dependence in effective nuclear interactions. Therefore, from the beginning, even in many cases today, nuclear DFT does not start with a density functional, but rather with a density-dependent interaction $V[\rho]$. The underlying relation is clear as the corresponding functional can be derived from the Hartree-Fock approximation,

[TABLE]

where $|\Phi\rangle$ is a single Slater determinant and $T$ is the kinetic operator.

Besides the above difference, the way to derive the functional is also quite different in atomic and nuclear systems. Because the interaction in atoms, the Coulomb force, is well known and has a simple form, one can derive the density functional microscopically without any additional parameters [249]. On the contrary it is very difficult for nuclear systems because of the complexity of nuclear force [73, 74]. Essential degrees of freedom in nuclear interaction include the spin, isospin, pairing, and relativistic degrees of freedom. Furthermore, the components of the nuclear interaction include not only the central force but also the spin-orbit force, tensor terms, and so on.

2.3 Current status in nuclear density functional

Because of the above challenges, all current successful nuclear density functionals are based on phenomenological parameters, which are determined by fitting to selected experimental ground-state (and excited-state) properties. The form of the functional, i.e. the structure of the various terms, is chosen based on the consideration of symmetries and simplicity, for example, the translational and rotational symmetry in the non-relativistic case, and the Poincare symmetry (including the Lorentz invariance) in the relativistic case.

The importance of some terms, such as the central force and the spin-orbit force, is easy to be seen from the experimental observables. This might not be the case for other terms, as for example, for the tensor force. It is known that the tensor force is an essential component in the bare nucleon-nucleon interaction from the scattering data and deuteron properties [165]. However, the role of the tensor force in nuclear density functional is much less clear [130]. In non-relativistic functionals, the zero-range tensor terms can be easily included [250], but their strengths are hard to be adjusted as they have little influences on the binding energies or charge radii. Even though the tensor terms have strong influences on the single-particle energies, according to the Kohn-Sham concept, in the framework of DFT these single-particle energies are purely artificial quantities and should not be used in a fit. Indeed, the single-particle energies in nuclei are usually fragmented, and thus beyond-mean-field effects are crucial, such as particle-vibrational coupling [251, 252]. Therefore, it is difficult to take into account such effects in the fitting both conceptually and technically.

2.4 Towards ab initio nuclear density functional

In such a case mentioned above, the derivation of the density functional from ab initio calculations is important although such a derivation for nuclear systems is much more difficult than its counterpart for atomic systems and it is still at the infancy stage.

Following the concepts in Coulombic systems, firstly the nuclear matter problem is solved, where the volume term can be studied. Then one could try to derive exact results for gradient terms, surface terms, and so on. Along with this direction, Fayans proposed to start with a density functional for nuclear matter from an ab initio calculation and to add on top it surface, Coulomb and spin-orbit terms, reducing in this way the total number of phenomenological parameters [253]. This idea has been realized in Ref. [254] together with further developments [255, 256, 257, 258, 259], where the functional is based on the results of microscopic nuclear and neutron matter calculations with the Brueckner $G$-matrix at various densities. With only four to five adjustable parameters, the authors have shown that it is enough to reproduce nuclear binding energies and radii with the same quality as obtained with other commonly used effective interactions. Progress has also been made for the relativistic density functional. In Ref. [260], the scalar-isovector term, which in most cases has been neglected for simplicity, is included in the functional DD-ME$\delta$ by fitting to the isovector effective mass difference derived from the RBHF theory in nuclear matter. Together with the constraints from relativistic and non-relativistic BHF results in nuclear matter, only four phenomenological parameters in DD-ME$\delta$ need to be adjusted to experimental data in finite nuclei. Of course, at this moment, it does not necessarily mean that this functional is superior to other relativistic density functionals. For example, the DD-ME$\delta$ functional does not predict octupole deformation in actinide nuclei which are known to be octupole deformed [261], and it gives significantly lower inner fission barriers in super-heavy nuclei [262].

These promising results demonstrate the feasibility in deriving a microscopic ab initio density functional for nuclei. However, the concept to derive functionals only from the properties of infinite nuclear matter cannot teach us much in cases, where one wants to describe effects which do not show up in nuclear matter calculations, e.g., the influence of tensor terms in spin non-saturated systems. For that, we need ab initio calculations for finite nuclei. In particular, for covariant density functional theory, we need the ab initio calculations for finite nuclei in a relativistic framework. Therefore, in the following Section 3, we will introduce in detail the progress of RBHF theory for both nuclear matter and finite nuclei, and in Section 4 we will discuss explicitly how the RBHF results are able to guide the developments of covariant DFT.

3 Relativistic Brueckner-Hartree-Fock Theory and Applications

3.1 BHF formalism

Reviews for non-relativistic Brueckner theory can be found in Refs. [132, 137, 263, 264, 265, 266]. For a nuclear system composed of $A$ nucleons, the Hamiltonian including a one-body kinetic-energy term $T$ and a two-body interaction term $V$ reads as,

[TABLE]

where $a_{k}^{\dagger}$ and $a_{k}$ form a complete set of creation and annihilation operators, referring to the state $|k\rangle$ with the single-particle wave function $\psi_{k}(\mathbf{r})$, and the one-body and two-body matrix elements are

[TABLE]

Within BHF theory, the ground state of the nuclear system is given by the Hartree-Fock approximation,

[TABLE]

where $|0\rangle$ is the vacuum and $a_{k}|0\rangle\equiv 0$. The total energy of the system can be derived as

[TABLE]

with the antisymmetrized two-body matrix element $\langle kl|\bar{V}|kl\rangle=\langle kl|V|kl\rangle-\langle kl|V|lk\rangle$. The single-particle states $|k\rangle$ satisfy the HF equation

[TABLE]

where $e_{k}$ is the eigenenergy of the state $|k\rangle$, and $U^{\rm HF}$ is the HF single-particle potential

[TABLE]

However, as the bare nuclear interaction has a very strong repulsive core [13], it cannot be used directly in the HF approximation. In modern ab initio calculations, the renormalization group method has been widely used to decouple the low-energy physics, which is of the interest, from the high-energy degrees of freedom, which is induced by the repulsive core. This greatly simplifies the numerical problem and benefits much the modern ab initio calculations [267].

In BHF theory, the bare interaction $V$ is replaced by an effective interaction, the $G$-matrix in the nuclear medium. It satisfies the Bethe-Goldstone equation,

[TABLE]

where $W$ is the starting energy, and $e_{m},e_{n}$ are the single-particle energies of the two particles in the intermediate states. The Pauli operator $Q(m,n)$ forbids the two particles being scattered to the already occupied states below the Fermi surface, i.e., $Q(m,n)=1$ for $e_{m}>e_{F}$ and $e_{n}>e_{F}$, otherwise $Q(m,n)=0$. Thus, in the BHF framework, the ground-state energy can be expressed as,

[TABLE]

The basis and correctness of the replacement of the bare interaction with an effective interaction of this form in the nuclear medium is rooted in many-body perturbation theory, or more specifically, in the Goldstone expansion [136]. Using Goldstone diagrams to represent the contributions of the full perturbation series, such a replacement is equivalent to consider the series of ladder diagrams in all orders of of $V$. Numerically, the matrix elements of $G$ are much smaller than those of $V$. Therefore $G$-matrix is more suitable to be applied in perturbation theory. Already at the first order the expression (14) gives most of the contributions [137].

Following the replacement of $V$ by $G$, the definition of the HF single-particle potential in Eq. (12) has to be modified. Unlike the expression of the total energy in Eq. (14) where the starting energy is well defined, the single-particle potential in Eq. (12) is an auxiliary quantity from the beginning and its definition depends on specific choices. Even if the form of the HF single-particle potential is used, i.e.,

[TABLE]

it still remains a problem how to choose the starting energy $W$. Fortunately, it was proven [141] that for occupied states $k$ a large amount of higher-order diagrams cancel with each other, if the starting energy is chosen on shell: $W=e_{k}+e_{l}$. This is good for a faster convergence and improves the BHF approximation.

However, there is no exact guideline for the choice of the starting energy when $k$ is an unoccupied state. At this point, there exists an uncertainty in the BHF scheme, and different recipes have been proposed, e.g., the gap choice ($U=0$) [268], the continuous choice (the same form as the one when $k$ is occupied) [269], or choices in between the above two. In principle, the auxiliary single-particle potential is introduced to make the convergence of the perturbation expansion faster. If the expansion is carried out to high enough order, the potential itself does not affect the final result. In Ref. [270] it is shown that in the next order beyond the BHF, i.e., the three hole-line expansion, the equation of state for nuclear matter near or below the saturation density calculated with the gap choice and the continuous choice agree with each other. Therefore, it seems reasonable to use a choice in between the gap choice and continuous choice.

Summarizing these discussions, the single-particle potential in the BHF framework can be chosen as

[TABLE]

where $e^{\prime}$ is uncertain and its choice for finite nuclei has been discussed both in the non-relativistic BHF [138, 271] and relativistic BHF [37, 232] theories. In nuclear matter, only the diagonal matrix elements of $U$ exist, that is $k^{\prime}=k$ in the above expression. Whereas the off-diagonal matrix elements appear in finite nuclei, and the related discussion can be seen in Ref. [138].

3.2 RBHF formalism and its application for nuclear matter

3.2.1 RBHF formalism

In the relativistic Brueckner-Hartree-Fock framework, the nucleon-nucleon (NN) interaction is of covariant form. Using the one-boson-exchange Bonn interaction as an example, the nuclear interaction is mediated by six bosons [165],

[TABLE]

including the scalar (s) $(\sigma,\delta)$, the pseudovector (pv) $(\eta,\pi)$, and the vector (v) $(\omega,\rho)$ meson-nucleon couplings, and each pair contains isoscalar or isovector character. For the isovector mesons, an additional isospin operator $\vec{\tau}$ is present. The nucleon field is denoted by $\psi$, and operator $\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$. Without specification, bold letters will be used for three-vectors (space only) such as $\mathbf{x},\mathbf{q},\bm{\gamma}$, and usual letters will be used for four-vectors (time and space) such as $x,q,\gamma$. By default, the Einstein summation convention is adopted. For the isovector couplings, the densities or currents in Eqs. (17) are replaced by the isovector ones, namely $\bar{\psi}\vec{\tau}\Gamma\psi$ with $\Gamma=1,\gamma^{\mu},\gamma^{5}\gamma^{\mu}$ and $\sigma^{\mu\nu}$.

In free space, the above Lagrangians lead to the following OBE two-body matrix elements (or OBE amplitudes) in the plane-wave helicity representation [165]

[TABLE]

where the reference frame is chosen at zero center-of-mass momentum. The momenta for two initial (final) states are $\mathbf{p}_{1}=\mathbf{q}$ and $\mathbf{p}_{2}=-\mathbf{q}$ ($\mathbf{p}_{1}^{\prime}=\mathbf{q}^{\prime}$ and $\mathbf{p}_{2}^{\prime}=-\mathbf{q}^{\prime}$), respectively. The Dirac spinor with positive energy is labelled as $u$, which in the helicity representation is

[TABLE]

with $E_{\mathbf{q}}^{2}=M^{2}+\mathbf{q}^{2}$. Notice that with this form, the Dirac spinor is normalized covariantly as $\bar{u}(\mathbf{q},\lambda)u(\mathbf{q},\lambda)=1$. This is convenient for solving the scattering equation as the two-body matrix elements in Eqs. (18) are Lorentz scalars. For calculating physical observables in many-body theory, we should have creation and annihilation operators obeying Fermi commutation relations and the normalization condition should be $u^{\dagger}u=1$, and the spinor becomes

[TABLE]

In the Bonn interaction, a form factor is applied to each meson-nucleon vertex [165]:

[TABLE]

where $m_{\phi}$ is the mass for meson $\phi$. The cut-off parameters $\Lambda_{\phi}$, together with the coupling strengths $g_{\phi}$, are determined by fitting to the NN scattering data and deuteron properties. Three examples of nuclear interactions, Bonn A, Bonn B, and Bonn C, have been obtained with the Thompson equation [165], in which the average nucleon mass $M=938.926$ MeV is used. For the pseudoscalar mesons, the coupling strength of the pseudovector channel ($f_{ps}$) is related to the pseudoscalar one ($g_{ps}$) by

[TABLE]

In the relativistic framework, the two-nucleon scattering equation is the Bethe-Salpeter equation [272], which can be written in the operator form as

[TABLE]

where $\mathscr{M}$ is the relativistic two-nucleon scattering amplitude and $\mathscr{G}$ is the relativistic free two-nucleon propagator. In this equation, $\mathscr{V}$ is the infinite sum of all irreducible diagrams and is usually replaced by the one-boson exchange contribution within the “ladder” approximation, $\mathscr{V}\approx V_{\rm OBE}$ (the subscript “OBE” will be omitted without causing ambiguity).

The Bethe-Salpeter equation (23) is a four-dimensional integral equation, and it is very difficult to solve. In practice three-dimensional reductions are used, and such reductions still remain covariant and the relativistic elastic unitarity is satisfied. Since unitarity does not uniquely determine the Green’s function, the three-dimensional reduction is not unique and in principle an infinite number of choices exist [273]. Typically a propagator $g$ is introduced,

[TABLE]

so that the scattering equation (24a) involves only a three-dimensional integration. With the approximation $\mathscr{W}\approx V$, Eq. (24a) becomes

[TABLE]

which can be easily solved in momentum space.

Take the Thompson choice as an example, one finds the Thompson equation [274] (the spin indices have been omitted) as,

[TABLE]

where $\mathbf{P}$ is the total momentum of two particles, and $\mathbf{q}^{\prime},\mathbf{k}$, and $\mathbf{q}$ are the relative momenta of the initial, intermediate, and final states respectively. For the OBE interaction, the expressions of the two-body interaction matrix element $V(\mathbf{q}^{\prime},\mathbf{q})$ are given in Eqs. (18). By solving the Thompson equation (26), the scattering amplitudes and phase shifts can be calculated and compared with experimental data. Three examples for the bare nuclear interaction, Bonn A, Bonn B, and Bonn C, have been provided by Ref. [165].

The Hamiltonian density is obtained by a Legendre transformation, and a three-dimensional integration leads to the total Hamiltonian $H$ in the stationary case,

[TABLE]

where $\phi_{i}$ represents the nucleon field $\psi$, the meson fields $\phi$, and the photon field $A$.

Eliminating the meson fields, the many-body Hamiltonian for nuclear systems is derived as [275],

[TABLE]

where the field operators $\hat{\psi}$ (and $\hat{\bar{\psi}}=\hat{\psi}^{\dagger}\gamma^{0}$) have been hatted in order to distinguish them from the single-particle wave functions, and $\Gamma_{\phi}^{(1)}$ and $\Gamma_{\phi}^{(2)}$ are the interaction vertices for particles at the coordinates $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, respectively,

[TABLE]

In Minkowski space, the meson propagators $D_{\phi}(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, and the sign $-$ holds for the scalar and pseudoscalar mesons and $+$ for the vector ones. The dependence on the zero-component of momentum transfer $q_{0}$ (energy) reflects the retardation of the interaction. For the Bonn interaction in Ref. [165], this effect was deemed to be small and was ignored from the beginning. In this way, the meson propagators are just Yukawa functions:

[TABLE]

Notice that, with the form factor in Eq. (21), the meson propagators are no longer simple Yukawa functions.

The nucleon-field operators $\hat{\psi}(\mathbf{r})$ and $\hat{\psi}^{{\dagger}}(\mathbf{r})$ can be expanded on a static relativistic basis $|k\rangle$,

[TABLE]

where $a_{k}^{\dagger}$ and $a_{k}$ form a complete set of creation and annihilation operators of nucleons, referred to the state $|k\rangle$ with positive or negative energies, and $\psi_{k}(\mathbf{r})$ is the corresponding Dirac spinor. The quantum number $k$ characterizing the state $|k\rangle$ contains also the isospin $t=n,\,p$ for neutrons and protons. Then the Hamiltonian for the nuclear system in second quantized form reads,

[TABLE]

where the matrix elements are given by

[TABLE]

The two-body interaction $V$ contains the contributions from various meson (photon) coupling channels $\phi$. The indices $k,l$ run over an complete basis of Dirac spinors with positive and negative energies, for instance, the plane wave states $u(\mathbf{k},s)$ and $v(\mathbf{k},s)$ in momentum space [276] or the eigensolutions of a Dirac equation with Woods-Saxon potentials [277, 229]. The two-body interaction matrix elements for the OBE model in the plane-wave helicity basis have been given in Eqs. (18), and the retardation is ignored in the Thompson choice so that $(q^{\prime}-q)_{0}=0$ in these equations.

3.2.2 Nuclear matter

In the relativistic scheme, there are solutions with both positive and negative single-particle energies. The no-sea approximation [101] is usually adopted, i.e. the ground state of nuclear system is taken as a simple Slater determinant composed of the single-particle states in the Fermi sea. The nucleon Dirac equation for nuclear matter reads

[TABLE]

The single-particle energy in momentum space is represented by a capital letter “$E_{\mathbf{p}}$”, which is equivalent to the lower case “$e_{k}$” in Eq. (11). Another convention is the single-particle potential $U$, which is usually referred as self-energy labelled by $\Sigma$. Although Eq. (35) does not distinguish solutions with positive or negative energy, usually the spinor $u(\mathbf{p},s)$ is used to represent a positive-energy solution and the spinor $v(\mathbf{p},s)$ a solution with negative energy.

In free space where $U=0$, the plane-wave solutions for Eq. (35) in helicity representation are given in Eq. (19). For nuclear matter with good parity, time-reversal invariance, and hermiticity, the single-particle potential has in general the form [89]

[TABLE]

Introducing the starred quantities,

[TABLE]

where ${E^{*}}^{2}={M^{*}}^{2}+{\mathbf{p}^{*}}^{2}$ and $\hat{\mathbf{p}}$ is the unit vector of $\mathbf{p}$, the positive-energy solution of the Dirac equation (35) is

[TABLE]

Unlike the covariantly normalized spinor in Eq. (19), the spinor in Eq. (38) is normalized in the standard way $u^{\dagger}u=1$. Different from the free space, the spinors in nuclear matter (38) depend on the single-particle potentials (self-energies). They are often called the “dressed” spinors. Changing the covariantly normalized spinor to the standardly normalized spinor will give an extra factor $\sqrt{E_{\mathbf{p}}^{*}/M_{\mathbf{p}}^{*}}$.

The RHF single-particle potential in nuclear matter can be obtained from Eq. (12) as

[TABLE]

where the integration is over the Fermi sphere with radius ${k_{F}}_{t_{2}}$, the Fermi momentum of particle $2$.

With the relativistic two-body matrix element (18), the in-medium spinor (38) and the components (36), the single-particle potentials for the isoscalar mesons are [203],

[TABLE]

In these expressions, the convention for the isospin is $t=0$ for proton and $t=1$ for neutron. The scalar density $\rho_{s}$ and baryon density $\rho_{b}$ are defined as

[TABLE]

For meson $i$, the definitions of $\Theta_{i}(p,q)$ and $\Phi_{i}(p,q)$ are

[TABLE]

For the isovector mesons, the single-particle potentials are almost identical to Eqs. (40) with the replacements: $\rho_{s}\to\rho_{s,t_{1}}-\rho_{s,1-t_{1}},\,\rho_{b}\to\rho_{b,t_{1}}-\rho_{b,1-t_{1}},\,\delta_{t_{1}t_{2}}\to 2-\delta_{t_{1}t_{2}}$. For the Bonn potentials [165], there is no isoscalar-tensor channel, i.e., $f_{\omega}=0$.

In the relativistic framework, the single-particle wave functions in nuclear matter must be solved iteratively. As shown in Eqs. (38) and (37), the single-particle wave functions depend on the potentials $U_{s}(p),U_{0}(p),$ and $U_{v}(p)$, which are calculated from Eqs. (40). After solving the RHF equation iteratively, the total energy (10) is

[TABLE]

In the RBHF calculation, the scattering equation in the nuclear medium has to be solved, and then the bare interaction can be replaced by the effective one, i.e., scattering matrix in the nuclear medium. Following the work by Brockmann and Machleidt [160], replacing the plane wave in free space (19) by the plane wave in nuclear medium (38) and adding the Pauli operator, the relativistic scattering Thompson equation in the nuclear medium can be obtained as

[TABLE]

Several points should be commented for obtaining Eq. (44) from the Thompson equation in the free space (26):

The scattering matrix $\mathscr{M}$ has been relabelled as $G$ to be consistent with the BHF convention. 2. 2.

The two-body matrix elements $V(\mathbf{q}^{\prime},\mathbf{q})$ have the same forms as free space (18) but with the dressed spinors (38). 3. 3.

The starting energy is denoted as $W$, for on-shell definition $W=2E_{\mathbf{P}/2+\mathbf{q}}^{*}$. 4. 4.

The Pauli operator $Q(\mathbf{k},\mathbf{P})$ ensures that the two intermediate particles are not scattered into the occupied states. 5. 5.

Same as in Ref. [160], the considered approximations include replacing $E_{\mathbf{k}}$ by $E_{\mathbf{P}/2+\mathbf{k}}$, $E_{\mathbf{q}}$ by $E_{\mathbf{P}/2+\mathbf{q}}$, and the angle-averaging approximation, i.e., the energy of the intermediate particle $E_{\mathbf{P}/2+\mathbf{k}}=E_{\mathbf{P}/2-\mathbf{k}}$ and the energy of the initial particle $E_{\mathbf{P}/2+\mathbf{q}}=E_{\mathbf{P}/2-\mathbf{q}}$. Then the energy denominator is the difference between two initial single-particle energies and two intermediate ones. 6. 6.

Because of the choices of three-dimensional reduction for the Bethe-Salpeter equation, and the approximations mentioned above, the relativistic Bethe-Goldstone (RBG) equation is not unique.

Equation (44) is formally identical to the non-relativistic BG equation in nuclear matter , and usually solved in the helicity basis. By decomposing the two-body matrix elements in different partial-wave channels $V^{JT}(q^{\prime},q)$, with $J$ the total angular momentum and $T$ the total isospin, the scattering equation is solved in different $J,T$ channels [51].

In Ref. [160], the RBG equation is solved in the rest frame. Alternatively, one can solve the $G$ matrix in the two-nucleon center-of-mass frame, as usually done in the calculations of two-nucleon scattering, and then transform to the rest frame to calculate the single-particle potential (16) and ground-state energy (14). This has been adopted, for example, by the Brooklyn group [156, 163], Horowitz and Serot [157, 158], and the Groningen group [161, 162].

In the relativistic framework, the single-particle wave functions are not known from the beginning and must be solved iteratively. The difficulty is that the single-particle wave functions (38), which are needed to calculate the two-body matrix elements $V(\mathbf{q}^{\prime},\mathbf{q})$, and the single-particle energies appearing in the denominator of RBG equation (44), depend explicitly on different components of the single-particle potential $U_{s},U_{0},U_{v}$ (37). On the other hand, the effective interaction $G$-matrix has mixed all the channels in solving the RBG equation (44) and the different components of the single-particle potential cannot be easily calculated using the form of Eqs. (40). In RBHF theory, several methods have been proposed to deal with this problem, including the solution in the full basis [156, 278, 163, 279, 189, 188, 280, 281, 282, 283], projection [157, 158, 161, 162, 187, 284, 285, 286, 183, 184], and a momentum-dependence analysis [159, 160, 287].

In Table 1, the saturation properties of symmetric nuclear matter calculated by RBHF using the interactions Bonn A, B, and C [160] are shown for different methods. For the results of the projection method, because of the uncertainty mentioned above, different schemes of projections have been used. Details are given in the corresponding references. Neglecting the very different results from Ref. [169], the energies per nucleon by different methods are receptively $-16.49\sim-15.59$ for Bonn A, $-15.73\sim-13.60$ for Bonn B, and $-14.38\sim-12.26$ for Bonn C, and the saturation densities $0.174\sim 0.189$ for Bonn A, $0.159\sim 0.174$ for Bonn B, and $0.138\sim 0.170$ for Bonn C. It is clear that the uncertainties with the different methods exist and an accurate calculation of the single-particle potential is still an open problem for RBHF theory in nuclear matter.

Besides the uncertainty in calculating different components of the single-particle potential, other uncertainties may exist because of the adopted approximations, such as the angle-average approximation of the Pauli operator [230, 288, 231, 289] or the average approximation of total momentum. Recently, an exact treatment for the total momentum in RBHF has been reported [290].

Table 2 shows properties of symmetric and asymmetric nuclear matter at saturation density obtained by RBHF theory using the Bonn A, B, and C interactions for the exact and for averaged center-of-mass momentum. In Ref. [290], the calculation of the self-energy is based on a momentum-dependence analysis similar to that of Ref. [160] used in Table 1, but there are uncertainties, such as which momenta are chosen to calculate the self-energy. The values given by the different calculations may also slightly depend on numerical details.

The incompressibility $K_{\infty}$ and skewness parameter $Q_{0}$ are defined by the expansion of the binding energy per nucleon $E/A$ in symmetric nuclear matter around the saturation density $\rho_{0}$ [291],

[TABLE]

The slope of the incompressibility is $M_{0}=3\rho\frac{\partial K(\rho)}{\partial\rho}|_{\rho=\rho_{0}}$. The slope parameter $L$ and curvature parameter $K_{\rm sym}$ are defined by the expansion of the symmetry energy $E_{\rm sym}$ around the saturation density [291],

[TABLE]

The incompressibility for finite nuclei is parameterized as,

[TABLE]

with $\alpha=(\rho_{n}-\rho_{p})/\rho,K_{\tau}=K_{\rm sym}-6L-\frac{Q_{0}}{K_{\infty}}L,$ and $K_{\rm Coul}=\frac{3}{5}\frac{e^{2}}{r_{0}}\left(-8-\frac{Q_{0}}{K_{\infty}}\right)$. When $Q_{0}\approx 0$, one defines parameter $K_{\tau}\approx K_{\rm asy}=K_{\rm sym}-6L$.

The results of RBHF in Table 2 are compared with those obtained by non-relativistic BHF theory with and without three-body force (TBF) [160, 292, 293]. Results of two sets of TBF, namely TBFa [294] and TBFb [293], are presented. The empirical values are also given for comparison [295]. It can be seen that the average approximation is acceptable for properties like $\rho_{0},E/A,K_{\infty},E_{\rm sym}$, and $L$, but not necessarily for the higher-order parameters such as $Q_{0},M_{0},K_{\rm sym},K_{\rm asy}$, and $K_{\rm Coul}$.

In Fig. 1, the equations of state for symmetric nuclear matter calculated with Bonn A, B, and C interactions by RBHF theory are shown in comparison with the results by non-relativistic BHF theory [160]. The empirical saturation properties $E/A=-15.8\pm 0.3$ MeV and $k_{F}=1.319\pm 0.014$ fm*-1* are also plotted as the black star [295]. The shadowed bands, which connect the three saturation points obtained by the interactions Bonn A, B, and C, are a schematic representation of the Coester lines. It was originally found in Ref. [146] that the saturation points obtained in BHF calculations with interactions of different tensor strength locate on a line that systematically deviates from the empirical region. As the three Bonn interactions are different in their tensor strengths, their nuclear matter saturation points are at different locations on the Coester line, as seen in Fig. 1. For the results of both RBHF and BHF, Bonn A with the weakest tensor force gives the largest binding energy and saturation density, and Bonn C with the strongest tensor force gives the smallest.

Looking at the non-relativistic results, the Coester line connecting the saturation points given by different interactions misses the empirical region systematically. Compared to the non-relativistic results, the relativistic framework leads to considerably improved results and the saturation points given by different interactions form a new “Coester” line, which is much closer to the empirical region. The relativistic effect produces a strong density-dependent repulsion, which can be explained by the following simple estimation in Ref. [176].

Considering the simplest case, the effective interaction $G$-matrix can be viewed as the effective interaction in the Walecka model with the $\sigma$ and $\omega$ mesons in the mean-field (Hartree) approximation. The single-particle potentials in Eqs. (40) become

[TABLE]

Typical values of the scalar and vector potentials are $U_{s}=-400(\rho/\rho_{0})$ MeV and $U_{0}=300(\rho/\rho_{0})$ MeV around the saturation density. From free space to nuclear matter, the positive-energy wave functions change as

[TABLE]

with $M$ being replaced by the effective mass $M^{*}=M+U_{s}$. Here $p_{0}=\sqrt{\mathbf{p}^{2}+M^{2}}$.

In free space, the relativistic effect reveals itself as a kinetimatic effect: the larger the velocity (or equivalently the momentum $\mathbf{p}$) the larger the relativistic effect (i.e. the small component becomes larger). For $\mathbf{p}=0$ the spinor becomes trivial unity and the wave function is the same as the non-relativistic plane wave. As $\mathbf{p}$ increases, the small component becomes larger and the relativistic wave function becomes more different from the non-relativistic plane wave.

In nuclear matter, situation gets more complicated as there is an additional origin of the relativistic effect, the dynamic effect: the larger the scalar potential $U_{s}$, the larger the relativistic effect. In other words, the small component can become non-negligible even in the small velocity limit, as $M^{*}$ can be significantly smaller than the bare mass $M$. As discussed in Section 1.4 , this has the consequence that in nuclear physics where the average velocity of a nucleon is generally small, effects of relativity can still be important.

To evaluate the size of the relativistic effect caused by replacing $M$ with $M^{*}$, one can expand the positive single-particle energy as

[TABLE]

In the small velocity limit $|\mathbf{p}|<<M$, the leading relativistic correction to the energy is

[TABLE]

Given the average kinetic energy per nucleon is $\langle\mathbf{p}^{2}/(2M)\rangle\approx 23(\rho/\rho_{0})^{2/3}$ MeV, the average correction can be evaluated as [176]

[TABLE]

Taking the results of the Bonn B interaction in Fig. 1 as an example, the relativistic correction to the energy per nucleon is well fitted by this ansatz with a different factor $\overline{\delta E}\approx 2(\rho/\rho_{0})^{8/3}$ MeV [160].

The Dirac spinor in nuclear matter $u(\mathbf{p},s)$ can be decomposed in the complete basis of free Dirac spinor (both positive-energy $u_{0}(\mathbf{p},s)$ and negative-energy solutions $v_{0}(\mathbf{p},s)$) [156],

[TABLE]

with the expansion coefficients $a(p)$ and $b(p)$ normalized to $a^{2}(p)+b^{2}(p)=1$. In first-order perturbation theory, the correction to the positive-energy wave function is [176]

[TABLE]

where $\Lambda_{-}(p)=(p_{0}-H_{D})/(2p_{0})$ is the Casimir projection operator, $H_{D}=\bm{\alpha}\cdot\mathbf{p}+\beta M$ is the single-particle Dirac Hamiltonian. In the small velocity limit and ignoring the correction to the large component, one has the first-order correction to $u(\mathbf{p},s)$,

[TABLE]

The mixing of negative-energy states in nuclear matter can be illustrated by the so-called Z-diagram in Fig. 2. Representing the matrix element of each dashed line as $-(U_{s}/M)\bm{\sigma}\cdot\mathbf{p}$ and the energy denominator $(2M)^{-1}$, the correction to the energy can be calculated to have the same value as in Eq. (51). This correction is included naturally in the dressed relativistic spinor $u$ (38) when one evaluates the full energy (50). If one wants to calculate this in a non-relativistic basis, where the spinor $u_{0}$ is used, an additional term (with matrix element $-(U_{s}/M)\bm{\sigma}\cdot\mathbf{p}$) has to be added to the non-relativistic energy $E_{0}$. This matrix element is caused by the difference between the relativistic dressed spinor and the non-relativistic spinor $\delta u$: $u=u_{0}+\delta u=au_{0}+bv_{0}$, with the momentum-dependent expansion coefficients $a$ and $b$ of the positive- and negative-energy free spinors. Such a mixing of the free spinors with positive and negative energies in the dressed spinor $u$ is interpreted in the non-relativistic framework as an admixture of nucleon-antinucleon excitations. In other words, in leading order the relativistic effect caused by the strong scalar potential can be expressed in the non-relativistic framework as a three-body interaction caused by an intermediate nucleon-antinucleon excitation.

In Ref. [296], a three-body interaction was constructed from the two-body Paris interaction, and among different contributions, the effect of the nucleon-antinucleon excitation term was confirmed with the above discussed properties [54]. The same three-body interaction has been used together with the AV18 two-body interaction [77] in a non-relativistic BHF calculation and the results have been compared with RBHF using the Bonn B interaction [155]. After including the $2\sigma$-exchange three-body interaction with the intermediate nucleon-antinucleon excitation in the BHF framework, the results are very close to those of full RBHF theory. Of course, one should also notice that in the given three-body interaction [296], there are other important contributions such as the $\sigma\omega$ nucleon-antinucleon excitation, Roper resonances, etc. [155]. The other contributions cancel each other more or less and the net effect is slight attraction.

In the modern nuclear forces derived from chiral effective field theory, three-body interactions, as well as higher many-body interactions, appear on the same footing as the two-body interaction in the chiral expansions. It is therefore interesting to compare relativistic effects and the chiral three-body interaction. In Ref. [178], RBHF and non-relativistic BHF calculations with Bonn B have been compared with non-relativistic BHF calculations with N3LO without and with the three-body interaction [66]. The results are shown in Fig. 3. By comparing the results of BHF with Bonn B and with N3LO, it can be seen that these two interactions give very similar results in the same framework. By comparing the changes from BHF to RBHF and the changes from N3LO to N3LO plus three-body interaction, it is evident that the effect of a relativistic treatment is very similar to the effect of the three-body interaction in the non-relativistic framework.

3.3 RBHF formalism for finite nuclei

Using the Thompson choice, the relativistic Bethe-Goldstone (RBG) equation for finite nuclei has the form,

[TABLE]

The quantum numbers $a,b,c,\dots$ represent the single-particle states being the solutions of RBHF equations. The indices $k,l,m,\dots$ will be used for an arbitrary complete relativistic single-particle basis. The Pauli operator $Q$ has the same definition as in BHF theory, i.e., $Q(c,d)=1$ for $e_{c}>e_{F}$ and $e_{d}>e_{F}$, otherwise $Q(c,d)=0$.

Adopting the choice in between the gap choice and the continuous choice, the single-particle potential of RBHF theory is calculated similar as in Eq. (16),

[TABLE]

where some uncertainty exists in $e^{\prime}$ and its choice has been discussed in Refs. [37, 232]. This is exactly the same as Eq. (16) with only different labels for quantum states, except that now the single-particle energies $e_{a}$ and $e_{b}$ may correspond also to states with negative energies, to states in the Dirac sea in the relativistic framework.

With this single-particle potential, one can solve the RHF equation,

[TABLE]

which is also formally the same as the non-relativistic HF equation (11). The matrix elements of the kinetic energy are given in Eq. (34a), and the RHF equation can be solved by the standard diagonalization. The RBHF equations for finite nuclei (56,57,58) will be solved iteratively until convergence is reached. Different from the RBHF equations for nuclear matter where different components of single-particle potential must be decomposed before the solution can be obtained, the decomposition of single-particle potential is not necessary for finite nuclei. Only the matrix elements of single-particle potential (57) are needed, which can be calculated directly from the $G$-matrix.

This difference between the calculations in finite nuclei presented here and those in nuclear matter discussed in Section 3.2.2 is caused by the fact that the solution of the RBG equations (56) in finite nuclei involves not only the matrix elements between pairs of states $|ab\rangle$ with positive energies ($e_{a}>0,e_{b}>0$), but also matrix elements between all possible pairs ($e_{a}>0,e_{b}<0$), ($e_{a}<0,e_{b}<0$) etc., which are needed for the solution of the corresponding RHF equation (58). In the usual nuclear matter calculations discussed in Section 3.2.2, the RBG equations (44) are only solved for pairs with positive energies. The matrix elements of the $G$-matrix for all other pairs, which are needed for the solution of the RHF equations (35), are directly calculated from the potential $U$, i.e., one needs the decomposition of $U$ into the various relativistic channels, scalar, vector, etc., and this decomposition is not unique. Of course, a full solution of the RBG equations in nuclear matter for all pairs with arbitrary energies would solve this problem. However, this is complicated and usually not done. The fact that one solves in finite nuclei the full set of RBG equations increases the numerical efforts considerably.

The total energy in the RBHF framework is (14)

[TABLE]

To solve the RBHF equations, one must start with a complete initial basis $\{|k\rangle\}$, as the solution $\{|a\rangle\}$ is not known from the beginning. The RBHF single-particle states are expressed as linear combinations of the basis states

[TABLE]

with $D_{ka}$ the expansion coefficients.

As the requirement of the completeness of the basis, $|k\rangle$ runs not only over the positive-energy states but also over the negative-energy states in the Dirac sea [229]. One should not mix this with the no-sea approximation adopted in the RBHF or RHF calculations, which refers to neglecting the Dirac sea in calculating densities and currents. For the RBHF calculation, the input Bonn NN interaction was constructed without considering antinucleon degrees of freedom [165], and this will be kept in the many-body calculation.

For the case with spherical symmetry, the eigenfunctions of the Dirac equation can be written as

[TABLE]

where the radial, orbital angular momentum, total angular momentum, and magnetic quantum numbers are denoted by $n,l,j,$ and $m$, respectively, and the spinor spherical harmonics are

[TABLE]

In practice, the relativistic quantum number $\kappa$ is often used, i.e., $\kappa=\pm\left(j+1/2\right)$, for $j=l\mp 1/2$. The orbital angular moment for the lower component in Eq. (61), $\tilde{l}=2j-l$, is different from the one for the upper component. The wave functions $F(r)$ and $G(r)$ satisfy the radial equation

[TABLE]

where $\Sigma=V+S$ and $\Delta=V-S$ are the sum and the difference of vector and scalar potentials. The initial basis is chosen as the Dirac Woods-Saxon (DWS) basis which is obtained by solving the radial Dirac equation (63) with the Woods-Saxon potentials $\Sigma(r)$ and $\Delta(r)$ [229].

Being consistent with the no-sea approximation, the index $c$ in Eq. (57) is restricted to the $A$ occupied states in the Fermi sea. The RHF equation (58) in the basis $\{|k\rangle\}$ reads

[TABLE]

As both the RBG equation (56) and the single-particle potential (57) are defined in the RBHF single-particle basis $\{|a\rangle\}$, it requires a nested iteration procedure.

Start with a set of trial single-particle states, e.g., a discrete DWS basis [229] with the single-particle wave functions $|k\rangle$ and the corresponding energies $e_{k}$. 2. 2.

Calculate the matrix elements of the kinetic energy $T_{k^{\prime}k}$ and the antisymmetrized two-body matrix elements of the bare interaction $\bar{V}_{k^{\prime}l^{\prime}kl}$ in this basis. 3. 3.

Solve the the RBG equation (56) in this basis by matrix inversion. This yields a set of $G$-matrix elements $\bar{G}_{k^{\prime}l^{\prime}kl}$. 4. 4.

For the fixed $G$-matrix, solve the non-linear RHF equation (58) by a nested inner iteration:

1. (a)

Use these $G$-matrix elements to calculate the single-particle potentials $U_{k^{\prime}k}$ in Eq. (57). With the basis transformation coefficient $D_{kk^{\prime}}$, it becomes,

[TABLE]

In the first iteration, $D_{kk^{\prime}}$ is a unity matrix. 2. (b)

Solve the RHF equation (58) with these matrix elements $T_{k^{\prime}k}$ and $U_{k^{\prime}k}$:

[TABLE]

and a new set of single-particle states $\{|a\rangle\}$ with single-particle energies $e_{a}$ is obtained. 3. (c)

Calculate a new relativistic single-particle density by occupation of the A orbits in the local Fermi sea ( $0<e_{a}\leq e_{F}$ ) for the next step of the RHF-iteration and a new potential $U$ with the same $G$-matrix, and continue this inner RHF-iteration till convergence. 5. 5.

Compare the converged RHF solution $\{|a\rangle\}$ with the initial basis $\{|k\rangle\}$ in which the $G$-matrix is calculated. If $\{|k\rangle\}=\{|a\rangle\}$ then the whole RBHF iteration is converged. Otherwise the RHF basis $\{|a\rangle\}$ will be used in the next iteration in Step 2. Alternatively, basis transformation can be used instead of calculating the matrix elements as in Step 2,

[TABLE]

3.4 Local density approximation on RBHF and its promotion to CDFT

From the early 1990s, one of the links between the RBHF theory and the covariant density functional theory involves the local density approximation (LDA), which was first introduced in the non-relativistic framework by Brueckner et al. [297]. The $G$-matrix in a finite nucleus is assumed depend on the local density $\rho(\mathbf{r})$ and it is the same as the $G$-matrix in nuclear matter with the same density. Then the $G$-matrix is parameterized as a density-dependent (DD) effective interaction in the HF framework. In order to account for the discrepancy between the RHF results and experimental data, on top of the effective interaction derived from the $G$-matrix, higher-order correction are then parameterized [298, 38, 299, 300].

With the success of RBHF theory in describing nuclear matter, it is natural to describe finite nuclei in the same framework with the LDA. In general, an effective Lagrangian is constructed in a similar way as the bare Bonn interaction in Eqs. (17), but now with DD coupling strengths,

[TABLE]

By fitting to the results of RBHF for nuclear matter, the effective interaction, i.e., the DD coupling strengths $g(\rho)$, can be determined. In practice, one can define the effective interaction in the RMF or in the RHF approximation. Accordingly the obtained effective interactions will be used in DDRMF or DDRHF calculations. There are many choices to fit the RBHF results, e.g., fits to the single-particle potential or fits to the $G$-matrix elements.

Although self-consistent RBHF theory for finite nuclei has been improved recently [36, 37], self-consistent RBHF calculations for all nuclei in the nuclear chart are still far above the numerical capabilities of the most modern computers. Definitely it is expected that RBHF results for both nuclear matter and finite nuclei will provide more and more insight for the nuclear density functional.

The LDA in the RBHF theory was initiated by the Brooklyn group [301, 302]. The relativistic $G$-matrix in nuclear matter was calculated using the interactions HEA [303] and HM2 [304] by taking into account relativistic effects in first-order perturbation theory. The single-particle potential calculated by RBHF is then fitted to determine the effective interaction in the RHF framework with the mesons, $\sigma,\omega$, and $\delta$ [301]. The density-dependence of the effective interactions is parameterized by analytical functions, and nuclear matter saturation curve can be well reproduced [301, 302].

Marcos et al. [305] improved the input by self-consistent RBHF calculation for asymmetric nuclear matter [162]. The total energy and different components of the single-particle potential calculated by RBHF are fitted in the RHF framework by an effective interaction with four mesons, $\sigma,\omega,\pi$, and $\rho$, including form factors. Further work using self-consistent RBHF results for nuclear matter [159] as an input can be found in Ref. [306].

In the above works, the single-particle potential and optionally the total energy has been used in the fit to determine the effective interaction. The $G$-matrix elements can also be used to determine the effective interaction, and in principle they should contain more information [307].

Adopting a parametrization of the $G$-matrix calculated by RBHF using the Bonn C interaction in Ref. [160] as the form $G=V+\Delta G$, the first part $V$ was determined by fitting to the matrix elements of the original bare interaction, and the second part $\Delta G$ containing four mesons, $\sigma,\omega,\delta$, and $\rho$, together with the first part were then fitted to the $G$-matrix [307]. The correlation term, $\Delta G$, mainly reduces sizeably the repulsive $\omega$-exchange [307]. By discarding the form factors in the parametrization, the resulting effective interaction  [307] became simpler, but it was still difficult in DDRHF calculations for the study of finite nuclei. Therefore, a simple and accurate, in the sense of reproducing the RBHF results, effective interaction is important in describing properties of finite nuclei.

With the input of RBHF [162] and the non-relativistic variational method [308], Gmuca [309, 310, 311] found that effective interactions with a $\sigma$-$\omega$ model with nonlinear $\sigma$-couplings, and a $\sigma$-$\omega$ model with a DD $m_{\sigma}$ [312] give similar results for nuclear matter, but quite different results for finite nuclei. Extending the nonlinear $\sigma$ model to include a nonlinear $\omega$ term, the effective RMF interaction is found to reproduce better the EoS derived from RBHF [310]. In Ref. [311], the finite nuclei 16O and 40Ca were studied using the effective interaction determined by RBHF calculations using the Bonn A, B, and C interactions [160] and the Groningen interaction [313]. Nonlinear effective interactions including $\sigma,\omega,\rho$, and $\pi$ mesons were further developed by Savushkin et al. [314].

Self-consistent RBHF calculations for nuclear matter by Brockmann and Machleidt using the newly developed Bonn A, B, and C interactions presented more detailed and more systematic results [160]. Jong and Malfliet improved the treatment of the singularity in the scattering matrix [313]. These results became very popular and were used directly to construct effective interactions.

A DD effective interaction with $\sigma$-$\omega$ mesons was developed in Refs. [111, 160] by fitting to the single-particle potential of nuclear matter derived from RBHF calculations with the Bonn A interaction. An extension to the DDRHF framework can be found in Ref. [193, 194]. Furthermore, the isovector $\rho$ meson without tensor was taken into account with a fixed coupling constant [315]. The Groningen group [313] developed a DDRHF effective interaction parameterized in the invariant Lorentz space by fitting to the $G$-matrix derived from RBHF calculations [187, 316]. In Ref. [187] an explicit density dependence was introduced in the fitting procedure, therefore such a parametrization can be applied to finite nuclei easily [317].

Up to this point, DDRMF and DDRHF do not include the so-called rearrangement term. This term appears if the single-particle potential is defined by the variation of the total energy with respect to the single-particle state,

[TABLE]

For a usual non-density-dependent two-body interaction, this will give the results of Eq. (12). However, for a interaction with DD coupling strengths, the variation will give an additional term,

[TABLE]

where the last term is the so-called rearrangement potential. In the HF total energy (10), the rearrangement term does not appear explicitly, but it has significant influence on the single-particle energies and wave functions and thus will change the total energy indirectly.

In (R)BHF theory, as discussed in Section 3.1, the definition of single-particle potential is not unique. It is chosen usually for a better convergence of the hole-line expansion. With the choice of Eq. (16), which is similar to the HF potential (12), there is no rearrangement term.

The concept of a “rearrangement term” in (R)BHF theory is similar but more complicated. It was first used by Brueckner to refer to the difference between the binding of a $A$-nucleon system and a ($A-1$)-nucleon system plus the 1 nucleon removal energy, due to the change of the Pauli operator in the BG equation when removing a particle [318]. Considering the single-particle energy derived in perturbation theory, the leading correction from Eq. (16) corresponds to another definition of the single-particle energy as [319, 320]

[TABLE]

where the occupation number $n_{i}$ is the expectation value of the number operator for the state $i$. This is the energy required to remove a particle from the system, leaving a hole in the state $i$. Therefore it satisfies the separation energy theorem of Hugenholtz and van Hove [321]. Another type is sometimes referred as the “orbital rearrangement”, in contrast to the “Brueckner rearrangement”, caused by the change of the self-consistent potential (or the wave function) [322, 323, 324, 325]. See also a discussion in Ref. [326] for these two rearrangement terms in BHF theory and the relation to Brandow’s definition of single-particle energy [327, 263], and the renormalized BHF theory [328, 329].

Back to the LDA for RBHF in nuclear matter, the missing rearrangement term in the calculation of finite nuclei using the DD effective interactions was investigated in detail in Refs. [330, 196]. The DDRMF effective interaction used [331] was based on the single-particle potential of RBHF in nuclear matter [160], and two density dependence schemes were used. For the vector density dependence (VDD), the baryon density $\rho_{b}=\rho_{0}$ (the [math]-component of vector) density) is used and for both the $\sigma$-coupling and $\omega$-coupling. For the scalar density dependence (SDD), the scalar density $\rho_{s}$ is used for $\sigma$-coupling while vector density for $\omega$-coupling. Both schemes including the rearrangement term improved the description of binding energy and charge radius simultaneously. Later, the effect of the rearrangement term was studied [332] using the DDRMF effective interaction based on RBHF calculations for nuclear matter in a full basis [279, 189], and similar results were found.

Most of the studies discussed so far are based on RBHF calculation for symmetric nuclear matter, and the isovector mesons were not included or only included with fixed phenomenological coupling constants. With RBHF results for asymmetric nuclear matter [333, 334], the DDRMF effective interaction with $\sigma,\omega,\rho$, and $\pi$ mesons were determined and applied to finite nuclei in Ref. [182]. However, the resulting coupling strength $g_{\rho}(\rho)$ decreased too much with the density and even became negative [182]. Using a parametrization of the $G$-matrix of Ref. [187, 316], it was demostrated that the single-particle potential determined by analyzing the momentum dependence in Ref. [160] might predict a “wrong sign” for the isovector dependence of $U_{s}$ and $U_{0}$ [181]. The RBHF calculation for asymmetric nuclear matter of the Giessen group adopted the projection method [286], and therefore it did not have the sign problem in the LDA study [197]. Another slightly different approach was adopted by the Tübingen group which used the so-called substracted $T$-matrix method to study asymmetric nuclear matter with RBHF [184, 186], and the results were later used for an LDA study [335]. In Ref. [199], by choosing all the mesons as density-dependent and determining the coupling strengths by fitting to the RBHF results of asymmetric nuclear matter, the DDRHF model with $\sigma,\omega,\rho,\pi$, and $\delta$ mesons yielded the best reproduction of the RBHF results in nuclear matter and provided a good description for finite nuclei.

As shown in the above discussion, RBHF calculations have made substantial progress and deepened our understanding for relativistic ab initio calculations of the nuclear system over the past decades. However, it is still a challenge to realize a fully self-consistent RBHF for nuclear matter in the full Dirac space with high precision. An exact treatment and the estimation of the approximations for the angle-average of total momentum [290] and Pauli operator [230, 231, 288] in the RBHF framework are also quite important.

On the other hand, the ab initio RBHF calculations provides valuable information and stimulate the development of nuclear (covariant) density functional theory. In the following, some of the examples on the guidance of the CDFT through RBHF will be discussed. Possible future opportunities to improve CDFT through RBHF calculation in finite nuclei will be discussed in Section 4.

The importance of the nonlinear coupling for the $\omega$ meson has been pointed out in the parametrization of effective interactions in nuclear matter [310]. By comparing the single-particle potentials derived from the RBHF calculations [160], Sugahara and Toki showed the necessity to introduce a nonlinear coupling for the $\omega$ meson [106].

Suggested by the microscopic ab initio RBHF calculations for nuclear matter [162, 160, 282, 286] and the related LDA studies [305, 111, 193, 331, 194, 187, 316], the DD effective interaction, later often referred as TW99, was proposed in Ref. [112]. After that, based on these general considerations, phenomenological CDFTs with density dependent coupling constants have been developed. By a careful adjustment of the parameters new DD effective interactions were developed which reproduced with high precision the available data not only for ground state properties, but also for collective rotational and vibrational excitations. They also have been successfully applied for density functional applications beyond mean field. For reviews see Refs. [18, 119, 122, 19]

Among these DDRMF functionals, DD-ME$\delta$ is another example how the ab initio RBHF studies can guide the development of nuclear density functionals [260]. Before that, most of relativistic functionals do not include the scalar-isovector $\delta$ meson, because a general fit to the binding energies and charge radii of finite nuclei is not sensitive to this degree of freedom. One needs ab initio information. Guided by the proton-neutron effective mass splitting in neutron matter calculated by RBHF with the Bonn A interaction [336, 186], the strength of the $\delta$ coupling can be determined, as shown in Fig. 4. Although the inclusion of the $\delta$ meson does not improve the accuracy of the properties of finite nuclei such as masses and radii (for example, an earlier parameter set DD-ME2 which does not include the $\delta$ meson outperforms DD-ME$\delta$ in a global description of the properties of finite nuclei [127]), the proton-neutron effective mass splitting is correctly incorporated by DD-ME$\delta$, and its EoS at higher densities is in a better agreement with experimental data derived from heavy-ion reactions than that of DD-ME2 [260].

The first DDRHF effective interaction PKO1, which has comparable precision for nuclear ground-state properties with the DDRMF ones, was developed in Ref. [212]. By including the Fock terms and the contribution of $\pi$ meson, it improves the descriptions of the nucleon effective mass and its isospin and energy dependence, while it describes well for bulk properties like binding energies and radii.

However, although some of the studies of LDA from RBHF in nuclear matter included $\pi$ or the $\rho$ tensor [315, 195, 199], their values could not be constrained from the properties of nuclear matter and were fixed to phenomenological values. In such a case, the RBHF study for finite nuclei is indispensable and we will show in Section 4 how nuclear density functional can be further improved through the information of finite nuclei.

3.5 Self-consistent RBHF calculations for finite nuclei

In this Section, we will introduce the details of the self-consistent RBHF calculations for finite nuclei, including the convergence check, the treatment of center-of-mass motion, the self-consistent basis, and the choice of the single-particle potential.

The numerical settings for the RBHF calculation for finite nuclei has been given in Ref. [37] in detail. For the nucleus 16O, the realistic bare $NN$ interaction Bonn A which is adjusted to the $NN$ scattering data [165] is used as the only input. The DWS basis, the solution of the Dirac equation with a Woods-Saxon potential, is then obtained by solving the spherical Dirac equation (63) in a box with the box size $R=7$ fm and mesh size $dr=0.05$ fm [229]. The procedure to solve the full RBHF equations has been introduced in Section 3.3. As one works always in the RHF basis during the iteration, the Pauli operator in Eq. (56) and all its relativistic structure is fully taken into account. In particular, there is no angle averaging involved, as such an averaging is involved in most Brueckner calculations for nuclear matter [230, 231]. The RBG equation (56) 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 [271].

3.5.1 Convergence check

It is well known that the bare $NN$ interaction contains a strong repulsive core and a strong tensor part which can connect the nucleons below the Fermi surface to the states with high momentum in the continuum. In most of the non-relativistic ab initio calculations, such a bare interaction is usually renormalized to a soft one while keeping the low-momentum properties unchanged. The renormalized interaction is then applied to the nuclear-structure calculations, see, e.g., the $V_{{\rm low}~{}k}$ [337] or the similarity renormalization group method [338]. The RBHF study for finite nuclei developed in Refs. [36, 37] uses directly the bare Bonn interactions [165], and therefore it is crucial to investigate the corresponding convergence of the RBHF calculations with respect to the basis space.

The total energy $E$ and the charge radius $r_{c}$ of 16O calculated by RBHF theory are shown in Fig. 5 as a function of the energy cut-off $\varepsilon_{\rm cut}$. In usual RMF (similarly in RHF) calculations with a DWS basis expansion for nuclei like 16O, a cut-off $\varepsilon_{\rm cut}=300$ MeV would be more than enough to achieve convergence with very high precision [229]. However, as seen in Fig. 5, this is not the case for RBHF calculations. To achieve a precision $\Delta E<1$ MeV, one needs a cut-off of at least $\varepsilon_{\rm cut}=1100$ MeV. This is due to the short-range correlation induced by the strong repulsive core of the bare interaction, which involves contributions of the states with high momentum in the continuum.

Other cut-offs, such as the box size, angular momentum, etc., have been carefully checked in Ref. [37] and will not be repeated here. For 16O, the final basis space used in the calculation is as following: the single-particle angular momentum cut-off $l_{\rm cut}=20\hbar$, the energy cut-off in the Fermi sea $\varepsilon_{\rm cut}=1100$ MeV and that in the Dirac sea $\varepsilon_{\rm cut}=-1700$ (with the condition that at least 2 negative-energy states are included in each block), the box size $R=7$ fm, the total angular momentum cut-off $J_{\rm cut}=6\hbar$ for $jj$-coupled two-body matrix elements, and the pair cut-off for the single-particle angular momentum and energy is applied when $J>3~{}\hbar$ as $l_{1}+l_{2}<20~{}\hbar$, $\varepsilon_{1}+\varepsilon_{2}<1100$ MeV. With these conditions, the total size of the two-body matrix elements stored in the computer is about 256 GB, and the computation time is about several days depends on specific machines. With the development of computational resources and possible optimization of parallelizing the RBHF code, the study of heavy nuclei such as 208Pb by the relativistic ab initio calculations is promising in the foreseen future.

3.5.2 Center-of-mass motion

For light finite nuclei, the treatment of the spurious center-of-mass (c.m.) motion is very important. 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}$. It has been shown in Refs. [339, 145] that, for large values of the particle number $A$, the projected energy is obtained in a good approximation by removing the c.m. 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 c.m. energy in Eq. (73) is removed after the variation Ref. [107]. This is a projection after variation (PAV). A strict treatment would be to exclude this term also in the (R)HF equation, i.e., to carry out a projection before variation (PBV) [340, 341]. In none of these cases, i.e., PAV and PBV, the c.m. term is included in the solution of the BG equation (56). These two approaches have been discussed in the non-relativistic BHF calculation long time ago, which were found to give similar results [342]. However, in Ref. [342], the BHF equations were not solved fully self-consistently, since the BG equation was only solved at the beginning and the $G$-matrix was not recalculated during the iteration.

In Table 3, the total energy, charge radius, and c.m. correction (73) for 16O calculated by RBHF for PBV and PAV are shown. It can be seen that although the values of the c.m. correction energy are similar for these two approaches, the corresponding total energies are quite different. PBV gives $9$ MeV more binding than PAV. This difference is much larger than that reported in Ref. [342]. From Eq. (73), one can conclude that the wave functions given by PBV and PAV are similar to each other, since the c.m. correction energies are similar. The difference in the total energy then can only be due to the difference in $G$-matrix for the two cases.

In order to understand this more clearly, at each step of the iteration the total energy is shown in the left panel of Fig. 6. It can be seen that there is no significant difference between PBV and PAV during the first RBHF iteration step, where the $G$-matrix $G_{1}$ is calculated from the initial DWS basis. One may say that with the same interaction, PAV can be viewed as a good approximation of PBV. While in the next RBHF iterations, the $G$-matrices $G_{2},\,G_{3},$ and $G_{4}$ are calculated each time in the converged RHF basis of the previous step, and we observe that the energy of PBV becomes smaller than that of PAV. The reason for such a sudden change can be understood from the right panel of Fig. 6, where the single-particle energies of the $s$ and $p$ blocks are given after the first ($G_{1}$) and the last ($G_{4}$) RBHF iterations. It can be seen that in the RBHF calculation the PBV approach leads in general to lower single-particle energies than the PAV, especially for the high-lying states. In the conventional RHF calculations, PBV and PAV approaches give similar results because only the occupied states are involved and they are similar for PAV and PBV. Whereas for the RBHF calculations, the difference in the single-particle spectra, in particular in the high-lying states, will lead to different $G$-matrices in next iterations, $G_{2},G_{3}$, and $G_{4}$. Lower unoccupied states generally coincide with more attractive $G$-matrix elements of occupied states. As a result, the $G$-matrix in PBV is more attractive and the total binding energy becomes larger.

3.5.3 Self-consistent basis

In the above discussion for the c.m. correction, two approaches, PAV and PBV, make a big difference if the RBHF equations are solved in the self-consistent RHF basis. This indicates the importance of using the self-consistent basis in the (R)BHF calculations for finite nuclei. Another way to see the importance of the self-consistent RHF basis is to examine the dependence of the final results on the chosen initial basis.

At the beginning, one does not know the final RBHF solution and therefore a trial basis is needed. To take into account the relativistic structure of the Pauli operator in Eq. (56), a relativistic basis is preferred. In Refs. [36, 37], the relativistic DWS basis [229] has been chosen as the initial trial basis. In comparison with the harmonic oscillator basis [93], it has advantages 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.

To show the importance of the self-consistent RHF basis, two types of calculations have been performed, i.e., solving the RBHF equations with a fixed initial basis [36] or with the self-consistent RHF basis [37]. In both calculations, the DWS basis is obtained by solving the radial Dirac equation  (63) with the Woods-Saxon potentials.

The total energies of both calculations at each iteration step are plotted in Fig. 7. The filled symbols stand for the self-consistent RBHF calculations, and the open ones for the fixed DWS basis calculations, and different types of symbols represent different DWS potential depths $V_{0}$. Similar to Fig. 6, the first RBHF iteration is represented by $G_{1}$. In Fig. 6 or Fig. 7, different calculations may have different number of iteration steps, but to have better comparison it has been adjusted slightly in the figures to make the number of iterations the same without losing too much precision. For the first RBHF iteration $G_{1}$, the self-consistent RBHF calculations are identical to the ones with fixed DWS basis for various $V_{0}$ values, because in both cases the same basis is used. After the first iteration $G_{1}$, there appears a distinct difference between the self-consistent calculations and those with the fixed basis. Moreover, the self-consistent RBHF calculations get converged with respect to the potential depths $V_{0}$, while those with a fixed basis do not. As expected, self-consistency is very important to get unambiguous results.

3.5.4 Choice of the single-particle potential

One uncertainty in (R)BHF theory is the choice of the single-particle potential for the particle states in Eq. (57). Different choices have been discussed in BHF for finite nuclei and more details can be found in Ref. [271]. Those were non-relativistic investigations and therefore one only had matrix elements $\langle a|U|b\rangle$ for the single-particle states $|a\rangle$ and $|b\rangle$ in the Fermi sea and above the Fermi level. In the relativistic framework, it becomes more complicated as the single-particle states $|a\rangle$ and $|b\rangle$ can have negative energies (in the Dirac sea) and the definition of single-particle potential for these states is another open question.

In Eq. (57), the single-particle states $|a\rangle$ and $|b\rangle$ in the Dirac sea were treated as occupied. This choice seems to be reasonable since in the Bethe-Goldstone equation (56) the intermediate states $|c\rangle$ and $|d\rangle$ are only allowed to be states above the Fermi surface $e_{c},\,e_{d}>e_{F}$. From this point of view, the single-particle states in the Dirac sea are “occupied” hole states.

In Ref. [37], these states were treated as unoccupied, i.e. in a similar way as the states above the Fermi level, and the formula is slightly different from that of Eq. (57) as

[TABLE]

where the index $i$ runs over the occupied states in the Fermi sea (no-sea approximation).

As discussed in Ref. [232], there is no “right” or “wrong” choice for the single-particle potential in (R)BHF theory, since (R)BHF theory can be viewed as the two-hole-line expansion in the more general hole-line expansion (or the Brueckner-Bethe-Goldstone expansion) [132] and as the expansion goes to higher order the results become independent of the choice of $U$ [270]. On the other hand, there do exist “better” choices of $U$ as they will affect the convergence rate of the hole-line expansion. It has been shown that the definition for hole states ($0<e\leq e_{F}$) in Eq. (16) cancels a certain large amount of higher-order diagrams and this accelerates the convergence of the hole-line expansion and improves the BHF approximation [141, 138]. However, there is no similar proof for the particle states nor for the states in the Dirac sea.

The effect of different choices are shown in Table 4 [37, 232] including,

Formula I, Eq. (74), and for $e_{a},e_{b}>e_{F}~{}\text{or}~{}<0$ the gap choice $\langle a|U|b\rangle=0$ is chosen. 2. 2.

Formula I, Eq. (74), and $e^{\prime}$ fixed as the energy of the lowest occupied states, $e_{\nu 1s1/2}$. 3. 3.

Formula I, Eq. (74), and $e^{\prime}$ fixed as the energy of the highest occupied states, $e_{\pi 1p1/2}$. 4. 4.

Formula II, Eq. (57), and $e^{\prime}$ fixed as the energy of the highest occupied states, $e_{\pi 1p1/2}$.

For Formula I, it can be seen that the gap choice gives the least bound system, while fixing $e^{\prime}$ as an energy among the occupied states gives $1.4$ to $1.7$ MeV per nucleon more binding. This is consistent with the non-relativistic BHF calculations. In the framework of the Bethe-Brueckner-Goldstone expansion (or hole-line expansion, see for instance Ref. [132]), the gap choice and the continuous choice have been discussed in nuclear matter in Ref. [270]. 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. They lie between the results of these two choices at the BHF level [270]. From this point of view, the choice of $e^{\prime}$ as an energy among the occupied states is a reasonable choice.

For Formula II, fixing $e^{\prime}=e_{\pi 1p1/2}$ gives $0.4$ MeV per nucleon more binding than Formula I and is in better agreement with the data. On the other hand, the rms charge radius is by $0.03$ fm smaller, and the SO splitting is smaller by $0.1$ MeV. The two formulas will have another major effect on the spin symmetries in the Dirac sea, and this will be shown in Section 3.7.

3.6 Ground-state properties for finite nuclei

In this Section, the ground-state properties for 4He, 16O, and 40Ca are studied by RBHF. Here we will show RBHF results for 16O and 40Ca using formula II (57) for the choice of the single-particle potential, instead of using formula I (74) as it was done in Ref. [37]. The formula II is preferred as the treatment of the single-particle states in the Dirac sea is similar to that of occupied states, and the resulting spin symmetry properties are more consistent with the findings in the studies of CDFT (see the discussion of next Section). The results of 4He in Ref. [37] will not be updated as RBHF theory is more focusing on medium-heavy and heavy nuclei, and for 4He there are exact solutions by solving the Faddeev-Yakubovsky equation, as for instance in Ref. [346]. Without specification, the center-of-mass correction is treated by the PBV method.

In Fig. 8 we show results of the RBHF calculations with the interactions Bonn A, B, and C [165], the energy per nucleon and the charge radius of 16O, and compare them with the renormalized BHF [191] and the relativistic effective density approximation (EDA) [191] using the same interactions. Similar as for the results of nuclear matter, we have connected the results of Bonn A, B, and C to indicate schematically the “Coester lines”. The experimental data [343, 344] is denoted by the black star.

By comparing the results of BHF and EDA with RBHF, it can be seen that the relativistic effects improve the description considerably, both for the energy and the radius. The “relativistic Coester line” is much closer to the experimental data, consistent with the results in nuclear matter. By comparing the results of EDA and RBHF, it can be seen that self-consistentcy is very important and it further improves the descriptions. One interesting point to be noted is that the three points of RBHF (with the Bonn A, B, C interactions) are closer to each other than those for BHF or EDA. This is similar as the results for nuclear matter, as shown in Fig. 1, where the three saturation points given by RBHF with Bonn A, B, and C are closer to each other than those of BHF. It should be mentioned that both calculations of BHF and EDA have adopted renormalized occupation probabilities, which takes into account certain rearrangement effects [328]. This has not been done for RBHF, and therefore the comparison in Fig. 8 is not fully on the equal grounds. If the results of the investigations of renormalized BHF hold also for RBHF, renormalized RBHF should further improve the description for finite nuclei, especially for the radii.

In Fig. 9, the single-particle energies for 16O calculated by RBHF using the interactions Bonn A, B, and C [165] are shown in comparison with the experimental data [345]. The results of RBHF are in a reasonable agreement with the data, especially for the spin-orbit splittings. However, the $p$-levels are slightly low. It might be due to the lack of more complicated configurations such as particle-vibration coupling [252] in the RBHF framework, where only the ladder diagrams have been included.

In Table 5 we list 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 [37, 232] with the interactions Bonn A, B, and C [165], in comparison with the experimental data [343, 344, 347, 345]. The RBHF results are different form Ref. [37] as they are updated with formula II for the single-particle potential (57) instead of formula I (74). We also include the corresponding results from the density-dependent relativistic Hartree-Fock (DDRHF) with the phenomenological interactions PKO1 [212] and PKA1 [207], the non-relativistic BHF results [348] with $V_{{\rm low}-k}$ derived from Argonne $v_{18}$ [77], the coupled-cluster (CC) method [349], the in-medium similarity renormalization group (IM-SRG) method [350], the no-core shell model (NCSM) [351], and the self-consistent Green’s function (SCGF) method [352] with N3LO [66], nuclear lattice effective field theory (NLEFT) [353] with N2LO [354], and the quantum Monte-Carlo (QMC) method [355] with the local chiral force N2LO [356]. All the results listed in Table 5 do not include bare three-body interactions.

For the results of IM-SRG and NCSM, the uncertainties induced by the extrapolation procedure are given in parentheses. For the result of NLEFT, the combined statistical and extrapolation errors are given. For QMC, the first error is statistical, and the second is based on the EFT expansion uncertainty. For RBHF, the errors of extrapolation with respect to the box size were estimated as 0.3 MeV for the total energy and 0.02 fm for the charge radius [37]. Other major uncertainties arise for the RBHF method due to cut-offs with respect to the single-particle basis space, the choice of starting energy, and three hole-line contributions. These shall be investigated in more detail in the future.

In the above references of non-relativistic ab initio calculations with chiral interactions, some have reported also results with 3N interactions. However, there are still large uncertainties, and the results with 3N forces obtained by different many-body methods do not agree as well as the results with 2N interactions. For completeness, we also cite some of the values. For example, IM-SRG gives $E=-130.5(1)$ MeV, similar to $E=-130.8(1)$ MeV for SCGF; but NCSM gives $E=-143.7(2)$ to $-147.8(1)$ MeV when changing the flow parameter from $0.05$ to $0.08$; NLEFT gives $E=-138.8(5)$ MeV; QMC gives $E=-117(5)(16)$ MeV and $r_{c}=2.71(5)(13)$ when the cutoff parameter is chosen as 1.0 fm, but the results are affected by this parameter to some extent.

The phenomenological functional PKO1 includes the tensor correlations induced by the $\pi$-exchange and PKA1 includes in addition the tensor correlations by the $\rho$-exchange. As expected, the phenomenological results, which are both fitted to the experimental data [212, 207], are in much better agreements with data than all the ab initio results.

The total energies given by the RBHF calculation with Bonn A, B, and C interactions are underbound by $7.4,20.5,$ and $29.6$ MeV (or by $6\%,16\%,$ and $23\%$), respectively. The charge radius is smaller by $0.17,0.11,$ and $0.06$ fm (or by $6\%,4\%,$ and $2\%$) compared to the experimental value, respectively. In the self-consistent non-relativistic BHF calculations [348], 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 value of $E=-120.9$ MeV is obtained with CC [349] using the chiral $NN$ interaction N3LO, $E=-122.8$ MeV is obtained with IM-SRG [350], $E=-119.7(6)$ MeV is obtained with NCSM [351], and $E=-122$ MeV is obtained with SCGF [352]. Using the same interaction N3LO, the calculated results agree with each other quite well and are similar to the results obtained by RBHF with Bonn A. The one obtained with the NLEFT [353] using the interaction N2LO, $E=-121.4(5)$, is also similar to the above values. On the other hand, with the two-body force only, QMC [355] gives $E=-87$ MeV, much less binding than the others. This might be because the N2LO force used there is a local one and different from that of Ref. [354].

The spin-orbit splittings of the $1p$ proton shell obtained in RBHF theory with the Bonn A interaction $\Delta E_{\pi 1p}^{ls}=5.3$ MeV is slightly smaller than the data but still in reasonable agreement ($16\%$ deviation). Those with Bonn B and C, which have stronger tensor force, give a further decreased spin-orbit splitting.

Figure 10 shows the charge-density distributions of 16O calculated by RBHF theory with the interactions Bonn A, B, and C, in comparison with experimental data [357]. As already reflected in the radius (see Fig. 8), here the central density distributions are too large. While Bonn C (with the strongest tensor force) gives the largest radius, it also describes the density distribution best.

Similar problems have also been found in the non-relativistic ab initio calculations, where the radii are too small even when the 3N interaction is included [358]. Only by using revised chiral interactions, for instance NNLOsat [359], which also includes the radii of finite nuclei such as 16O in the fitting, the situation is improved. However, this interaction cannot describe well the NN scattering data above 35 MeV, and thus the concept of ab initio calculation using this interaction is still questionable. Further investigations are needed for a better understanding of this problem.

Besides the radii, other physical observables, such as the single-particle energies and spin-orbit splittings, are also closely related to the properties of the density distribution. In the non-relativistic BHF calculations, the single-particle energies of 16O are shown to be more bound than the data [348]. For example, the single-particle energies of proton in the $1s_{1/2},1p_{3/2},1p_{1/2}$ orbits given by Ref. [348] are $-72.8,-37.2,-24.2$ MeV, respectively. On the other hand, this implies that the nucleons concentrate too much in the center and the mean-field potential is deeply overbound there. As can be expected, a large slope of the density leads to too large spin-orbit splittings: the proton $1p$ spin-orbit splitting given by BHF is 13 MeV (Table 5), considerably deviating from the experiment with 6.3 MeV. In RBHF, the situation has been improved much: the central density increases only slightly, and the single-particle energies are, on average, only slightly overbound as compared to the data (Fig. 9). Obviously, in a relativistic theory, the slope of density distribution is not the only decisive factor, because the spin-orbit splittings given by RBHF with the Bonn interactions are smaller than the data even though the density slope is larger.

As there are many investigations for LDA based on RBHF using the Bonn A interaction, it is interesting to compare the self-consistent and the LDA results. In Fig. 11, the energy per nucleon and charge radius calculated by RBHF using Bonn A interaction are compared with LDA derived from RBHF in nuclear matter using the same interaction. The LDA results include: effective interactions defined in the nonlinear relativistic mean field (NLRMF) model [311], the nonlinear relativistic Hartree-Fock (NLRHF) model [314], the DDRMF model [111, 193, 194, 315, 330, 332, 182, 335, 199], and the DDRHF model [193, 194, 315, 195, 199]. The filled (open) symbols are used for calculations with (without) the rearrangement term.

As ambiguities exist in the mappings from RBHF calculations for nuclear matter to the effective interactions for finite nuclei in LDA, the results for finite nuclei are very different. In this case, the self-consistent RBHF calculation for finite nuclei [36, 37, 232] is indispensable to provide a benchmark.

In Fig. 11, the difference between LDA results with and without the rearrangement term can also be seen. The filled symbols including the rearrangement term are generally located on the right-hand side of the corresponding open ones without the rearrangement term. In other words, the rearrangement term improves the description of the radius substantially. It also improves the description of the binding energy, but to a smaller extent.

Both 4He and 40Ca are calculated in RBHF theory with the Bonn A interaction as well [37]. In Table 6, the results for 4He are shown. The c.m. correction is taken into account with both PBV and PAV as introduced in Section 3.5.2. The corresponding results are compared with the experimental data [343, 344, 360], as well as those from DDRHF with PKO1 [212] and PKA1 [207], the Faddeev-Yakubovsky (FY) equation [361] with CD-Bonn [362], FY [346] with N4LO [68], no-core shell model (NCSM) [363] with N3LO [66], nuclear lattice effective field theory (NLEFT) [353] with N2LO, and BHF [348] with $V_{{\rm low}-k}$ derived from Argonne $v_{18}$ [77]. In all the above calculations, and also for the later results for 40Ca, only the two-body bare $NN$ interactions are used.

As expected, the center-of-mass correction plays an important role for such light nuclei. The binding energy given by PBV is much overbinding, and the one given by PAV is slightly underbinding. The radii of both are much larger than the data, and this is also observed in other calculations, including DDRHF with PKO1 and PKA1, FY equation, and NCSM.

In Table 7, the total energy, charge radius, matter radius, and proton spin-orbit splitting for the $1d$ orbit are listed for 40Ca calculated by RBHF with the interaction Bonn A [165], in comparison with the experimental data [343, 344, 345]. The RBHF results are different form the values given in Ref. [37] as they are updated with formula II for the single-particle potential (57), instead of formula I (74). The corresponding results from DDRHF with PKO1 [212] and PKA1 [207], BHF [348], NCSM [364], and CC [365] with $V_{{\rm low}-k}$ derived from Argonne $v_{18}$ [77], CC [366] and IM-SRG [350] with N3LO [66] are also included.

Similar as for 16O, the total energy of 40Ca calculated by RBHF with the interaction Bonn A is underbound by 36 MeV (or by 11%) and the charge radius is smaller by 0.26 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 of the data. Most of the non-relativistic results, because of the missing three-body force, give a too-large binding energy and a too-small radius, except for the CC and IM-SRG methods with N3LO which reproduce the experimental binding energy well. In the calculation by IM-SRG, the result with inclusion of chiral 3N interaction has also been reported as $E=-376.1$ [350], giving more overbinding than that with 2N interactions only. Again, as has been mentioned in the discussion of 16O, there are still large uncertainties in the 3N interactions and more investigations are called for to get a better understanding.

In Fig. 12, we show the localized single-particle potential and the single-particle energy levels calculated by RBHF [37] with Bonn A, in comparison with the experimental data [345]. In RBHF, the single-particle potential is indeed nonlocal, as seen in Eq. (57). Therefore, the one shown in Fig. 12 is obtained for a specific state $|a\rangle$ by the radial Dirac equation (63),

[TABLE]

In this way, the single-particle potential can be illustrated intuitively, but of course, it depends on the chosen state $|a\rangle$. In Fig. 12, the proton (neutron) $1s_{1/2}$ wave function with positive energy is used to obtain the localized single-particle potential. Similar to the case of 16O, RBHF with Bonn A gives reasonable SO splittings. Comparing with the experimental data, the gap in the Fermi surface is relatively large, which might be due to the lack of high-order configurations like the particle-vibration coupling [252].

In Fig. 13, the charge-density distribution calculated by RBHF with Bonn A interaction [37] is shown, in comparison with the experimental data [357]. Similar to 16O, the central density given by RBHF for 40Ca is too large (about $1.4$ times larger than the experimental value at $r=0$ fm).

In Fig. 14, the energy and the radius of 40Ca obtained by self-consistent RBHF with Bonn A are compared with various LDA results. Similar to Fig. 14, the LDA results are organized as: effective interactions defined in the NLRMF model [311], the NLRHF model [314], the DDRMF model [111, 193, 194, 315, 330, 332, 182, 335, 199], and the DDRHF model [193, 194, 315, 195, 199]. Filled (open) symbols are used for calculations with (without) the rearrangement term. The self-consistent result (red point) sits in the middle of those LDA results without the rearrangement term. Among them, the LDA results of Refs. [182, 314] are very close to the self-consistent one.

3.7 Pseudospin and spin symmetries

The spin symmetry and its breaking in the nuclear system is very important for nuclear physics. The large spin-orbit (SO) splitting, which was introduced by Goeppert-Mayer [7] and Haxel et al. [8] in 1949, formed the ground for the nuclear shell model. Twenty years later, the so-called pseudospin symmetry was proposed to explain the near degeneracy between two single-particle states with the quantum numbers $(n,l,j=l+1/2)$ and $(n-1,l+2,j=l+3/2)$ [367, 368]. Such two states are regarded as pseudospin doublets with the pseudospin quantum numbers $(\tilde{n}=n-1,\tilde{l}=l+1,j=\tilde{l}\pm 1/2)$. It was found that the angular momentum of the pseudospin doublets $\tilde{l}$ is nothing but the orbital angular momentum of the lower component of the Dirac spinor (61), and the pseudospin symmetry is exact when the sum of vector and scalar potentials $V+S$ vanishes [369], or a more general condition $d(V+S)/dr=0$ which can be approximately fulfilled in exotic nuclei [370, 371]. Since then, pseudospin symmetry as a relativistic symmetry has been realized and much work has been done for investigating its origin and its properties using phenomenological single-particle Hamiltonians, relativistic mean field theory, or relativistic Hartree-Fock (RHF) theory [372, 373, 374, 212, 213, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384].

In the Dirac equation, there exist solutions not only with positive energy but also with negative energy, i.e., the states in the Dirac sea. It was pointed out in Ref. [385] that the spin symmetry in the Dirac sea has the same origin as the pseudospin symmetry in the Fermi sea. In other words, the spin doublets in the Dirac sea have the quantum numbers $(n,\tilde{l},j=\tilde{l}\pm 1/2)$, and the spin symmetry breaking term is proportional to $d(V+S)/dr$, similar to the pseudospin symmetry in the Fermi sea. See the related works in Refs. [386, 387, 388, 389, 390, 391] and reviews [392, 124].

In most of these investigations, the pseudospin symmetry in the Fermi sea or the spin symmetry in the Dirac sea has been studied by starting from a phenomenological single-particle Hamiltonian, or the relativistic density functional [101, 18, 119, 122, 19]. It is therefore interesting to see to what extent the spin symmetry in the Dirac sea is found in RBHF calculations starting from the bare nucleon-nucleon interaction. As the single-particle properties of states in the Dirac sea is concerned, the two formulas on the choice of single-particle potential in the Dirac sea, namely, formula I (74) and formula II (57) should be investigated [232].

The single-particle spectrum and the effective single-particle potential $\Delta(r)=V(r)-S(r)$ in the Dirac sea calculated by RBHF theory with different choices of the single-particle potential $U$ are shown in Fig. 15. Similar to Eq. (75), the effective single-particle potential $\Delta(r)$ is obtained by the radial Dirac equation (63) for a specific state $|a\rangle$,

[TABLE]

The localized potentials shown in Fig. 15 are calculated from the wave function $\nu 1s_{1/2}$ in the Dirac sea (or $\nu 1\tilde{p}_{1/2}$ if labelled with the angular momentum of the lower component in Eq. (61)). The single-particle levels are grouped by the angular momentum $\tilde{l}$ of the lower component in the Dirac spinor (61) with negative energy, thus, $\tilde{s},\tilde{p},\tilde{d},\dots$ means $\tilde{l}=0,1,2,\dots$. The potentials in both panels are not approaching 0 when $r\to\infty$ as usually found in the RMF study [385] because of the nonlocality of the RBHF single-particle potential $U$ in Eq. (57).

By comparing the two panels in Fig. 15, it can be seen that formula II gives a deeper single-particle potential in the Dirac sea, and the spectrum is higher by $100\sim 200$ MeV. Even though the two choices give only $0.4$ MeV per nucleon difference for the binding energy (see Table 4), the difference in the Dirac sea is much more significant. This is understandable as the difference between these two choices lies in the definition of states in the Dirac sea only, and the effect on the states in the Fermi sea is indirect. Beside the shift of single-particle levels upward as a whole with formula II in Fig. 15, the SO splittings given by II are also generally smaller.

To see the SO splittings more clearly, in Fig. 16 the SO splittings $\Delta E_{\rm ls}=e_{j_{<}}-e_{j_{>}}$ versus the average energy of the spin doublets $e_{\rm av}=(e_{j_{<}}+e_{j_{>}})/2$ are shown, with $j_{<}=\tilde{l}-1/2$ and $j_{>}=\tilde{l}+1/2$. The results are compared with those of the phenomenological relativistic density functionals PKO1 [206] and PKDD [107]. It can be seen that the SO splittings calculated by RBHF using formula II are much smaller than those of I, thus the spin symmetry is better conserved and is in better agreement with phenomenological relativistic density functional findings. This may give some hints that formula II of the single-particle potential (57), where the states in Dirac sea are treated as occupied states, is preferred.

In short, the prediction of the spin symmetry in the single-particle spectrum in the Dirac sea is supported by the relativistic ab initio calculations by starting from the bare nucleon-nucleon interaction. In future, it would also be interesting to continue this investigation considering particle-vibration coupling, as it has been shown to have an important effect on the pseudospin symmetry [393].

4 RBHF Calculations for Neutron Drops

4.1 Neutron drops

A neutron drop is an ideal system, composed of a finite number of neutrons and confined in an external field to keep the neutrons bound. Similar to nuclear matter, it is not a real system existing in nature, but by studying it one can obtain rich knowledge of nuclear structure and various properties of droplets [291]. Since in such systems there exists only the neutron-neutron interaction, the equations for neutron drops are much easier to solve and therefore they have been investigated by many ab initio methods [394, 395, 396, 397, 233], as well as by DFT [394, 398]. In this way, different ab initio methods can be benchmarked and useful information can be extracted for the nuclear density functional.

The neutron drops were first studied by the quantum Monte-Carlo method [394] for $N=7$ and $8$ using the two-nucleon interaction Argonne $v_{18}$ [77] and the three-nucleon interaction Urbana IX  [399]. It was found that the commonly used Skyrme functionals overestimate the central density of these drops and the spin-orbit splitting of the $N=7$ neutron drop [394]. In Ref. [400], the ground-state energy was studied for a $N=6$ neutron drop, and the neutron pairing energy was discussed by comparison with Ref. [394]. Later, more systematic studies were performed for larger $N$ values with different external fields and different interactions using the quantum Monte-Carlo method [401, 402, 396, 403]. The studies with the modern high-precision chiral 2N interaction N3LO [66] and the 3N interaction N2LO [65] have been benchmarked with different ab initio methods, including the no-core shell model [32] and the coupled-cluster theory [31], and it was found that the results were consistent with each other [397]. However, by comparing these ab initio calculations, one found a significant dependence on the selected interactions, especially on the 3N interactions [396, 397, 403].

Recently, the effect of the tensor force in neutron drops was studied by ab initio relativistic Brueckner-Hartree-Fock (RBHF) theory with the Bonn A interaction [233, 234]. An evidence of the tensor force was found in the evolution of the spin-orbit splittings with increasing neutron number. This finding can provide a guidance for determining the strengths of effective tensor force in the nuclear medium.

Various non-relativistic and relativistic density functionals have been used to study neutron drops, and a strong linear correlation between the rms radii of neutron drops and the neutron-skin thickness of 208Pb and 48Ca has been pointed out in Ref. [398]. Because of the uncertainty in the isovector part of the density functionals, there is also an uncertainty in the results of neutron drops for different functionals. In Ref. [404], ab initio calculations for neutron drops were used to further constrain the density functionals by the microscopic results of effective field theories.

4.2 General properties

Without further specification, all the results of neutron drops shown below are calculated in an external harmonic oscillator field

[TABLE]

Figure 17 shows the total energy divided by $\hbar\omega N^{4/3}$ for $N$-neutron drops ($N$ from $4$ to $50$ for Bonn A and from $4$ to $28$ for Bonn B and C) in a HO trap ($\hbar\omega=10$ MeV) calculated by RBHF theory [233, 234]. The factor $\hbar\omega N^{4/3}$ is based on the consideration that, in the Thomas-Fermi approximation [145], the total energy for a non-interacting $N$-fermion system in a HO trap is given by

[TABLE]

In other words, the energy below the line $E/\hbar\omega N^{4/3}\approx 1.082$ corresponds to the binding due to the nuclear interactions. This intrinsic binding energy grows linearly with the HO strength $\hbar\omega$. The energy of the drops grows with the neutron number $N^{4/3}$, in contrast to the nuclear case where the binding energy grows with mass number $A$. For the cases of open shells, the filling approximation is used. The results are compared with quantum Monte-Carlo (QMC) calculations [402, 396] based on the interactions AV8’ + UIX, AV8’, and AV8’ + IL7, with no-core shell model (NCSM) calculations [397, 396] based on chiral 2N + 3N forces, chiral 2N forces, and the JISP16 force, and finally with calculations using relativistic density functionals [398, 206].

It can be seen that the results of Bonn A, B, and C are very similar. This is because the main difference among these three Bonn interactions is the strength of the $T=0$ tensor force [165], which has no influence on the $T=1$ neutron-neutron states. It is also consistent with the results in pure neutron matter, where the EoS calculated by RBHF with these three interactions are very close to each other [405]. Therefore, in later discussions only the results of Bonn A will be shown.

By comparing with the QMC and NCSM calculations, the results of RBHF with the interaction Bonn A are similar to the results with the JISP16 interaction, and AV8’ + IL7 (for $N\leq 14$), and getting closer to AV8’ (for $N\geq 20$). This similarity is favourable as JISP16 is a phenomenological nonlocal $NN$ interaction which can reproduce the scattering data and works well for light nuclei [406, 407]. On the other hand, AV8’ + IL7 gives a better description of light nuclei up to A = 12 than AV8’ or AV8’ + UIX, but it gives too much over-binding for pure neutron matter at higher densities [408, 396].

In comparison with the relativistic density functional calculations in the right panel of Fig. 17, four types of functionals have been chosen as the following,

RH nonlinear meson-exchange models: NL3 [104], PK1 [107]; 2. 2.

RH density-dependent meson-exchange models: DD-ME2 [114], PKDD [107]; 3. 3.

RH nonlinear point-coupling models: PC-PK1 [117]; 4. 4.

RHF density-dependent meson-exchange models: PKO1 [206] (which includes the tensor force from the $\pi$ meson).

As there is no pairing in the RBHF calculations, pairing is also excluded in the relativistic density functional calculations. It can be seen that the binding energies (attraction between neutrons) given by RBHF are generally larger than those given by DD-ME2 and PKDD. For $N=8$, RBHF is close to PKDD, but getting closer to PK1 from $N=14$ to 26, and closer to PC-PK1, NL3, and PKO1 from $N=28$ to 36. From $N=20$ on, the results of RBHF and DD-ME2 are close to a horizontal line, while the others have a small tendency of increasing. The microscopic results obtained by RBHF can be a guidance for future density functionals. For example, the neutron-neutron interaction might be in general too repulsive in DD-ME2, whereas in NL3 it might be too attractive when the neutron number N is small and then becomes too repulsive as N becomes large. This is consistent with the findings in the equation-of-state of neutron matter, where NL3 gives overbinding for small densities but underbinding for larger densities [113].

Some clues on the shell structure are seen in Fig. 17. For example, the energies of neutron drops at $N=8,20,40$ decrease suddenly comparing with their neighbours. This can be understood easily as the external HO potential has the magic numbers 8, 20, and 40.

In order to see the shell structure more clearly, the two-neutron separation energies for the above calculations are presented in Fig. 18. The HO magic numbers 8, 20, and 40 are clearly shown in all calculations. Beside the above magic numbers, the results of RBHF indicate sub-shell closures at $N=16$ and $N=32$, similar as the results of AV8’ + IL7. However, the $N=32$ gap given by AV8’ + IL7 might be too strong. The sub-shell closure at $N=32$ is not significant for AV8’, and does not exist for AV8’ + UIX. For the $N=28$ sub-shell closure, the results of Bonn A and AV8’ + UIX show a small hint, while AV8’ and AV8’ + IL7 do not show it. On the other hand, all the relativistic density functionals only show the HO magic numbers 8, 20, and 40, but no clear sub-shell closures for $N=16,28$ or $32$.

The left panel of Fig. 19 shows the density distributions of $N$-neutron drops in a HO trap ($\hbar\omega=10$ MeV) calculated by RBHF theory using the interaction Bonn A. With the given HO strength, the neutron density gets saturated around $0.14\sim 0.17$ fm*-3*.

A comparison between the density given by RBHF with Bonn A and other ab initio and density functional calculations was given in Ref. [404], which is shown in the right panel of Fig. 19. These calculations include: QMC with AV8’ + UIX (purple circles) [396], NCSM with JISP16 (green squares) [396]; density functionals SLy5 (black dash-dot) [341], SkM* (orange dash) [409], UNEDF0 (gray dash-dot) [410], YGLO (blue line) [411], KIDS (green line) [412], and ELYO (red line) [413]. In this figure the density distribution of $N=8$ neutron drop is taken as an example.

It can be seen that the density distribution given by RBHF with Bonn A is very similar to that of JISP16. Indeed (see Fig. 17) the energy of these two calculations are also close to each other. On the other hand, AV8’ + UIX shows more repulsion than Bonn A from the energy point of view (Fig. 17). This is consistent with the density distribution. As the neutrons feel less attraction and move more outside, AV8’ + UIX has lower central density. Among the selected functionals, the results of SLy5 and YGLO are close to RBHF with Bonn A.

Figure 20 shows the rms radii of $N$-neutron drops in a HO trap ($\hbar\omega=10$ MeV) calculated by RBHF theory using the interaction Bonn A. In panel (a) the results are compared with the QMC calculations based on the interaction AV8’ + UIX [402], with the NCSM calculations [397, 396] based on the chiral 2N + 3N and the JISP16 forces. In panel (b) these results are compared with calculations based on relativistic density functionals. The black line $R_{N}=2.118N^{1/6}$ fm is obtained by solving for free Fermions in a $\hbar\omega=10$ MeV HO trap using the Thomas-Fermi approximation as

[TABLE]

For $M=938.926$ MeV and $\hbar\omega=10$ MeV, one finds a factor $2.118$ fm in front of $N^{1/6}$. The black line $R_{N}=1.862N^{1/6}$ fm is obtained by fitting to the results of Bonn A from $N$ = 6 to 50.

In general, all the calculated radii fulfill the relationship $N^{1/6}$ as a function of $N$. In all the selected calculations, RBHF with Bonn A gives the smallest radii. By comparing with other calculations in Fig. 17 (a) and Fig. 20 (a), it can be seen that, while it gives the smallest binding energies, AV8’ + UIX also gives the largest radii. The energies given by JISP16 are similar to those of Bonn A, but the corresponding radii are larger. On the other hand, the radii of relativistic density functionals are much larger than those of RBHF, even though some of their binding energies are larger than RBHF before $N=20$ as shown in Fig. 17 (b). It is known that the relativistic density functionals without density-dependence in the isovector channel (Nl3, PK1, and PD-PK1) show too large neutron radii in realistic nuclei [113], and this can be observed for the neutron drops too.

Since the neutron matter properties are in close relation to the symmetry energy $a_{\rm sym}$ and its slope parameter $L$, Table 8 lists the radius of the $N=50$ neutron drop calculated by RBHF theory using Bonn A and the corresponding symmetry energy and the slope parameter calculated at nuclear matter saturation [414, 184, 415]. They are compared with the results of relativistic density functionals. It can be seen that in general the radius of a neutron drop is large if $a_{\rm sym}$ or $L$ is large, although in detail small discrepancies exist. For example, DD-ME2 gives the smallest $a_{\rm sym}$ and $L$, and its radius is indeed the smallest among those of relativistic functionals, but it is still larger than that of Bonn A. The radius of PC-PK1 is the largest, and its $a_{\rm sym}$ or $L$ are large, but not the largest. They are slightly smaller than those of NL3 and PK1.

The correlation between the slope parameter $L$ and the neutron-skin is well known [416]. In Ref. [398], a strong linear correlation has been found between the neutron skin thickness $\Delta r_{np}$ and the rms radius $R_{N}$ of $N$-neutron drops in an external HO field. Figure 21 shows the linear correlation between the neutron-skin thickness of 48Ca and the radius of $N=20$ neutron drops in a $\hbar\omega=10$ MeV HO external field. The black circle and square symbols are calculated with different non-relativistic and relativistic density functionals, and the blue line is obtained by fitting to those results [398]. The inner (outer) colored regions depict the 95% confidence (prediction) intervals of the linear regression.

The red star in Fig. 21 is calculated by RBHF theory using Bonn A. The datum of the neutron-skin thickness of 48Ca is obtained by measuring the electric dipole polarizability in Ref. [417]. It can be seen that the neutron-skin thickness of 48Ca given by RBHF $\Delta r_{np}=0.14$ fm is located within the error bar of the experimental data, which is also consistent with the $0.12\leq\Delta r_{np}\leq 0.15$ fm given by the coupled-cluster calculations using the interaction NNLOsat [418].

Apart from the linear correlation between $\Delta r_{np}$ of 48Ca and the radius of $N=20$ neutron drops in Fig. 21, similar correlations can be found in other cases, for example for $\Delta r_{np}$ of 208Pb or the radius of other $N$. Using these linear correlations, the experimental data of neutron skins of 48Ca and 208Pb can be mapped to the data of radii of neutron drops with different numbers of $N$ [398], and the results are shown with green symbols in Fig. 22. In this way, the study on the neutron skin of heavy nuclei can be linked to the study of the radius of neutron drops, while the latter is much easier to be accessed by different ab initio calculations.

In Fig. 22, the radii for $N=20,14$, and $8$ neutron drops calculated by RBHF theory using the interaction Bonn A (red dashed lines) are shown. They are compared with the data (green symbols) determined from the linear correlations with the neutron-skin thicknesses of 208Pb and 48Ca [398], and other ab initio calculations (blue dotted lines) [402, 396, 403]. For $\Delta r_{np}$ of 208Pb, the data come from different measurements with antiprotonic atoms [419] (circle), pion photoproduction [420] (square), and electric dipole polarizability [421] (diamond), respectively. For $\Delta r_{np}$ of 48Ca, the datum comes from the measurement of the electric dipole polarizability [417] (triangle). For the local chiral forces N2LO from Refs. [356, 422], the results include a two-body force with a cutoff $R_{0}=1.0$ and $1.2$ fm, and a two-body plus three-body force ($2N+3N$) with a cutoff $R_{0}=1.2$ fm [403]. Theoretical uncertainties are denoted by the blue colored regions. There is no particular reason to choose $N=20,14$, and $8$ neutron drops, as long as the central density of the neutron drop does not differ too much from the saturation density ($\approx 0.16$ fm*-3*) [398].

It can be seen that the radii obtained in the RBHF calculations with the interaction Bonn A are in a good agreement with the data determined from the linear correlations with the neutron-skin thicknesses. In comparison with other ab initio calculations, AV8’ + UIX shows more repulsion and gives larger radii, as expected from the energies shown in Fig. 17. For the $2N$ local chiral forces N2LO, the softer interaction with a cut-off radius $R_{0}=1.2$ fm gives a smaller radius and the harder one with $R_{0}=1.0$ fm gives a larger radius. When including the $3N$ force for N2LO, the radius gets larger by 0.05 fm and is similar to that of AV8’ + UIX.

4.3 Spin-orbit splitting and tensor effects

One important purpose of studying the neutron drops with RBHF theory is to pave the way to extract information for nuclear density functionals from ab initio calculations. Among the various properties of the neutron drops, the spin-orbit splitting is of particular interest, as, in general, it is difficult to find significant features in experimental data which are only connected to tensor forces and therefore suitable for an adjustment of their parameters. As shown in Ref. [233], the evolution of SO splittings in neutron drops gives a clear hint of the tensor force in the nuclear density functional, which has been the topic of long debates in the past [130].

In Fig. 23, we present the SO splittings of $N$-neutron drops in a HO trap for the $1p,\,1d,\,1f$, and $2p$ doublets calculated by RBHF theory using the Bonn A interaction. It shows the evolution of the various SO splittings with neutron number. In the left panel, the SO splittings are compared with the results obtained by various phenomenological relativistic mean-field (RMF) density functionals, including the nonlinear meson-exchange models NL3 [104] and PK1 [107], the density-dependent meson-exchange models DD-ME2 [114] and PKDD [107], and the nonlinear point-coupling model PC-PK1 [117].

For the microscopic RBHF results a clear pattern can be seen: The SO splitting of a specific orbit with orbital angular momentum $l$ decreases as the next higher $j=j_{>}=l+1/2$ orbit is filled and reaches a minimum when this orbit is fully occupied. When the number of neutron increases further,the $j=j_{<}=l-1/2$ orbit begins to be occupied and the SO splitting grows up again.

This pattern is similar to the monopole effect of the tensor force acting between neutrons and protons in nuclei as explained by Otsuka et al. [126]. In their explanation, the tensor force produces an attractive interaction between a proton in a SO aligned orbit with $j=j_{>}=l+1/2$ and a neutron in a SO anti-aligned orbit with $j^{\prime}=j^{\prime}_{<}=l^{\prime}-1/2$, and the tensor force produces a repulsive interaction between the same proton and a neutron in a SO aligned orbit with $j^{\prime}=j^{\prime}_{>}=l^{\prime}+1/2$ (or both in a SO anti-aligned orbit with $j^{\prime}=j^{\prime}_{<}=l^{\prime}-1/2$). It was also mentioned in Ref. [126] that a similar mechanism, but with a smaller amplitude, exists for the tensor interaction between the same type of nucleons with $T=1$.

The behavior of the SO splittings given by RBHF in Fig. 23 can be understood in a qualitative way: Consider, for instance, the decrease of the $1d$ SO splitting from $N=20$ to $N=28$. Above $N=20$ additional neutrons are filled into the $1f_{7/2}$ orbit (SO aligned). They show a repulsive interaction with the $1d_{5/2}$ neutrons (SO anti-aligned) and an attractive interaction with the $1d_{3/2}$ neutrons (SO anti-aligned). The repulsive (attractive) interaction increases (decreases) the corresponding single-particle energy. Therefore, by filling in neutrons into the $1f_{7/2}$ shell the $1d_{5/2}$ orbit is shifted upward and the $1d_{3/2}$ is shifted downward, reducing the $1d$ SO splitting. Above $N=28$ the neutrons are filled into $2p_{1/2}$ and $1f_{5/2}$. They interact with the $1d$-neutrons in the opposite way and increase again the SO-splitting for the $1d$ configuration.

However, in the left panel of Fig. 23, it turns out, that this specific evolution of SO splitting is not significant for any of the phenomenological RMF functionals, which do not include a tensor term. On the other hand, as shown in the right panel of Fig. 23, the same calculation but with the RHF density functional PKO1 [206], which includes the tensor force induced by the pion coupling through the exchange term, shows a similar pattern as the RBHF calculations, but the amplitude is smaller. Therefore, in order to investigate the effects of the tensor force, the pion coupling has been multiplied with a factor $\lambda$ without readjusting the other parameters of this functional. As expected, the evolution of the SO splitting is influenced by the strengths of the tensor force significantly. For $\lambda=1$, one has the results of the original functional PKO1, which already shows the right pattern, but the size of the effect is slightly too small. This can be understood by the fact that it is difficult to fit the strengths of the tensor force just to bulk properties such as binding energies and radii [217, 423]. The general feature of these SO splittings found in the RBHF calculations can be well reproduced with PKO1 simply by multiplying a factor $\lambda=1.3$ in front of the pion coupling.

4.4 Connections to CDFT

Investigations of nuclear matter using RBHF theory already promotes many successful developments of CDFT. Investigations of finite nuclear systems will provide further information to deepen our understanding of nuclear structure and to develop nuclear density functionals of higher quality. In particular, one of the important guides is the systematic and specific pattern due to the effects of the tensor force found in the evolution of SO splittings in neutron drops.

The tensor force in non-relativistic nuclear density functionals is a topic of long debates and it is still not fully settled [130]. The study of the tensor force in CDFT is still at its infancy and has attracted much attraction in the past decades [206, 207, 209, 424, 423, 425, 426, 211, 427, 428, 215, 216, 429, 382, 430, 218, 217]. In Ref. [206], the first RHF density functional PKO1 has been developed, which has a comparable precision of describing finite nuclei as other successful RMF functionals. With the inclusion of exchange terms and the contribution of the pion, the effect of the tensor force can be studied, and it improves the agreement with experimental single-particle levels [209]. The tensor force from the exchange of the $\rho$-meson was further developed in the RHF functional PKA1 in Ref. [207].

Even though it is well known that the tensor force mainly originates from the $\pi$-exchange, it is not an easy task to identify different contributions quantitatively and study the relativistic origin of the tensor force on an equal footing as in the non-relativistic context. The origin of the tensor force in CDFT is studied using a non-relativistic reduction and the one associated with the Fock diagrams of Lorentz scalar ($\sigma,\delta$) and vector ($\omega,\rho$) couplings are also discussed [216]. A more detailed study has been done very recently to identify the tensor force up to order $1/M^{2}$ in each meson-nucleon coupling using RHF theory [217]. With this newly developed formalism, one is eventually able to make a quantitative comparisons with the corresponding tensor force in the non-relativistic density functionals.

However, in most of the previous investigations on the influence of the tensor force for CDFT (or for non-relativistic DFT), the main focus was to compare with experimental single-particle levels [423]. For example, Fig. 24 shows the energy difference $\Delta E=E_{\pi 1h_{11/2}}-E_{\pi 1g_{7/2}}$ in $Z=50$ isotopes as a function of neutron excess calculated by the RHF functionals PKO1, PKO2, and PKO3, and by the RMF functional PKDD [209]. As the tensor force is (basically) not included in PKO2 and PKDD but included in PKO1 and PKO3, it can be seen that the tensor force can improve the agreement with the data. Indeed, this is also consistent with the non-relativistic findings [431, 432].

There might be several reasons why the inclusion of tensor components in the RHF calculations shown in Fig. 24 is still far from reproducing the experimental results. (i) The particle-vibrational coupling, which goes beyond mean field and is not taken into account here, has a strong influence on the single-particle energies (see below). (ii) In the PKO functionals, the pion-nucleon coupling strength has an exponential density dependence. The pion-nucleon coupling strength is fixed to the free value for vanishing density and it drops down very fast in the nuclear interior. This leads to a reduced tensor strength, because it has been shown in Ref. [423] that a fit to experimental masses and radii does not like a tensor term. In Fig. 25, we show the same energy splittings as those in Fig. 24 in results of RHF calculations with three functionals containing a fixed pion-nucleon coupling constant $\lambda f^{2}_{\pi}$ in addition to the conventional DD-ME functional [423]. In these calculations, only the Fock term of pion-nucleon coupling is taken into account and the values of $\lambda$ are chosen as [math], $0.5$, and $1.0$, respectively, and then the remaining parameters are adjusted with the same fitting protocol as DD-ME2 [114]. For $\lambda=0$, we have the usual RH theory, which does not show substantial increasing of the energy splitting. While the splittings are much too large for $\lambda=1.0$, one finds reasonable values for $\lambda=0.5$, i.e., the inclusion of a tensor term allows a very reasonable fit to the single-particle levels in this region and to the remaining ground-state properties.

On the other hand, as mentioned above, the particle-vibrational coupling (PVC) plays an important role in spherical nuclei [252, 393]. Going beyond mean field, this effect can be included consistently without any new fitting parameter, and it produces reasonable values not only for the single-particle gap in magic configurations [252] but also for the width of giant resonances [433]. In Ref. [428], the energy splittings in the Sn region were calculated with the RMF functional NL3*, which does not include a tensor force, both at the mean-field level and beyond mean-field with inclusion of PVC. It was found that the description of these energy differences can also be improved substantially when the particle-vibration coupling is included, as shown in Fig. 26.

Recent experiments found in the $N=20$ isotones below 40Ca that there are considerable changes in the $1f$ and $1p$ spin-orbit splittings, in particular between 36S and 34Si. On the basis of various relativistic and non-relativistic density functionals, Karakatsanis et al. [430] investigated several reasons for these experimental observations, including the difference in the isospin dependence of relativistic and non-relativistic spin-orbit couplings, tensor forces, pairing correlations, and particle-vibrational coupling. They found that all these effects had influences on the $1f$ and $1p$ spin-orbit splittings, of course depending in size on the various parameters entering in these theories. Since these parameters have been phenomenologically adjusted in one way or another, they found it very difficult to isolate specific reasons, such as tensor forces, in the nuclear density functional by just looking into experimental single-particle levels.

In such a situation, parameter free, microscopic ab initio calculations are definitely needed to get a clearer understanding of such phenomena and to learn more about the various parameters in the density functionals. For instance, as discussed in Section 4.3, the effect of the tensor force can be studied clearly in RBHF calculations of the evolution of spin-orbit splittings in neutron drops. As a relativistic description, this theory includes the proper parameter-free spin-orbit coupling. Such calculations produce meta-data. On the present level, they neither include pairing correlations, nor beyond-mean-field effects such as particle-vibrational couplings. This allows to adjust the tensor force to these meta-data without the ambiguity of additional effects, since they are neither included in the present concept of DFT nor in the present RBHF calculations.

5 Summary and Perspectives

Relativistic Brueckner-Hartree-Fock (RBHF) theory was investigated already in the 1980s for nuclear matter, and it has shown much improvement as compared with the non-relativistic calculations. Besides the development of RBHF theory and the important progress in understanding nuclear structure on the basis of ab initio calculations using the bare two-body nucleon-nucleon interaction, there are two other benefits of RBHF: First, it shows the relativistic origin of the three-body force in the non-relativistic framework. It is well known that the non-relativistic ab initio methods cannot describe the saturation property for nuclear matter using only two-body NN interactions and a three-body force has to be introduced with new adjustable parameters to reproduce the empirical values. Of course, the microscopic origin of the full three-body force in nuclei is much more complicated, and there are also discrepancies between experimental data and the RBHF results, which should be studied in future; still, the relativistic origin of three-body forces is an important component. Second, the studies of RBHF theory for nuclear matter stimulated the development of ab initio covariant density functional theories in the past. Examples include the successful density-dependent relativistic mean-field and density-dependent relativistic Hartree-Fock models, or the study of the importance of scalar-isoscalar channel by DD-ME$\delta$, and so on.

Due to the complexity of the theory and the limitation of computer resources, a fully self-consistent RBHF study for finite nuclei has not been realized until recently [36, 37]. Before that, RBHF studies for finite nuclei stayed at the level of the local density approximation (LDA). These LDA studies have deepened our understanding of the nuclear system and of the structure of covariant density functional theories. However, as an indirect way to study finite nuclei, by a mapping from results of nuclear matter, it suffers from considerable ambiguities. In order to solve this problem and also from the point of view of completeness, a self-consistent RBHF study for finite nuclei was called for. Starting from the Dirac Woods-Saxon basis and using the self-consistent RHF basis during the iteration, the coupled system RBHF equations (the relativistic Bethe-Goldstone equation and RHF equation) have been solved directly in a finite nuclear space [36, 37]. Convergence is studied carefully. The results show much improvement over the non-relativistic investigations, which confirms the importance of relativistic effects in finite nuclear systems. Moreover, full RBHF studies in finite nuclei not only solve the non-uniqueness problem of LDA, but also overcome the difficulty in decomposing the Lorentz structure of the single-particle potential which causes until today lots of troubles and uncertainties in the RBHF studies of nuclear matter.

The study of finite nuclei has further significance for the development of nuclear density functionals. For example, there are shell effects and tensor-force effects, which do not exist in nuclear matter.

The spin-symmetry of the single-particle levels in the Dirac sea of 16O has been investigated by RBHF theory and compared with results of covariant density functional theories. The studies on the spin- and pseudospin-symmetries beyond phenomenological functionals give us more insights of these symmetries and enables a better connection between CDFTs and RBHF calculations.

The investigation of neutron drops in the framework of RBHF theory provides many interesting insights. As an ideal system, a neutron drop is easy to be calculated and serves as a benchmark of different methods and interactions. By looking into the evolutions of spin-orbit splitting as the number of neutrons increases, a clear signal of the tensor force is revealed. This can provide guidance for the development of nuclear density functionals. For example, a non-relativistic functional with tensor terms, SAMi-T, has already been developed along this direction. It shows promising results for the descriptions of both ground-state and excited-state properties [235]. Counterparts in the relativistic framework are also in progress [217, 434].

With all these quite impressive successes, the way towards an ab initio covariant density functional theory is still at the beginning, much work is still needed.

First, the development of RBHF theory for nuclear matter is far from finished. The Lorentz-structure decomposition of the self-energy has not been fully solved, and certain approximations, such as angle-averaging of the Pauli operator, are existing. The study of nuclear matter by RBHF theory has important applications to nuclear astrophysics and for the construction of better nuclear density functionals through the LDA. Therefore, a fully self-consistent calculation with high precision is needed.

Second, even though an application of RBHF theory for finite nuclear systems has now been realized, the computational demand is vast and limits the further extension of the theory. Renormalization procedures such as $V_{{\rm low-}k}$ or the similarity renormalization group (SRG), which are commonly used in the non-relativistic ab initio context, would be very useful to reduce this computational cost and would also allow applications for heavier systems.

Third, through RBHF applications in nuclear matter and in finite nuclei, the descriptions have been improved over the non-relativistic ones, but they are still not perfect. As the lowest order of the hole-line expansion, the (R)BHF theory is not supposed to be the full story yet, and the next order (three hole-line) is known to have non-negligible contributions in the non-relativistic case. Therefore, how much the three hole-line contributes in the relativistic Brueckner theory is an open and interesting question. There is a possibility that, even with higher orders of the hole-line expansion added upon the current RBHF theory, discrepancies between theoretical and experimental results may still exist. Then one would have to consider the possibility to use bare three-body forces in the relativistic framework too. If the discrepancy between theoretical results with two-body interaction and experimental data can be attributed to the contribution of three-body interaction, as the relativistic effect already improves much of the description, the contribution of three-body force in RBHF should be smaller than in BHF.

Fourth, with higher orders of the hole-line expansion, diagrams beyond the ladder type are needed in principle. Therefore, how to go beyond the ladder diagrams in a systematic and controllable way is another open and interesting question. The connection to the recent progress in the method of the so-called diagrammatic quantum Monte-Carlo [435, 436] may be one of the possible ways.

Last but not least, the only input for present RBHF theory, i.e., the relativistic nucleon-nucleon interactions, which can reproduce the two-body scattering data and the deuteron properties, are much less developed than their non-relativistic counterparts. Up to now, the most commonly used relativistic NN interactions are the Groningen [162] and the Bonn interactions [165], which are in some sense outdated comparing with the non-relativistic modern and high-precision interactions. Therefore, efforts to develop relativistic chiral NN interactions are also highly demanded and work in this direction is in progress [71, 72]. Of course, one needs in this context a new power counting system [71], because in a relativistic theory all the terms of the $1/M$-expansion are summed up automatically. In the leading order of the relativistic chiral interaction, a much better description of the ${}^{1}S_{0}$ and ${}^{3}P_{0}$ phase shifts was achieved comparing with the non-relativistic case, and descriptions in other partial waves are similar to non-relativistic case [71]. There are indications that this improves the convergence properties of relativistic chiral perturbation theory.

Acknowledgements

We would like to express our gratitude to all the collaborators and colleagues who contributed to the investigations presented here, in particular to L. S. Geng, J. N. Hu, B. Long, X. L. Ren, H. Tong, S. B. Wang, S. Q. Zhang, and P. W. Zhao. This work was partly supported by the National Key R&D Program of China (Contract No. 2017YFE0116700 and No. 2018YFA0404400), Natural Science Foundation of China under Grants No. 11335002, and No. 11621131001, the Overseas Distinguished Professor Project from Ministry of Education of China No. MS2010BJDX001, the JSPS Grant-in-Aid for Early-Career Scientists under Grant No. 18K13549, the JSPS-NSFC Bilateral Program for Joint Research Project on Nuclear mass and life for unravelling mysteries of the r-process, the DFG (Germany) cluster of excellence “Origin and Structure of the Universe” (www.universecluster.de), the Funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 654002. H.L. would like to thank the RIKEN iTHES project and iTHEMS program.