Ab initio nuclear many-body perturbation calculations in the Hartree-Fock basis

Ab initio nuclear many-body perturbation calculations in the Hartree-Fock basis

B. S. Hu (ºú°Øɽ) State Key Laboratory of Nuclear Physics and Technology, School of Physics, Peking University, Beijing 100871, China    F. R. Xu (Ðí¸¦ÈÙ) State Key Laboratory of Nuclear Physics and Technology, School of Physics, Peking University, Beijing 100871, China    Z. H. Sun (ËïÖкÆ) State Key Laboratory of Nuclear Physics and Technology, School of Physics, Peking University, Beijing 100871, China    J. P. Vary Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA    T. Li (Àîͨ) State Key Laboratory of Nuclear Physics and Technology, School of Physics, Peking University, Beijing 100871, China
July 14, 2019

Starting from realistic nuclear forces, the chiral NLO and JISP16, we have applied many-body perturbation theory (MBPT) to the structure of closed-shell nuclei, He and O. The two-body NLO interaction is softened by a similarity renormalization group transformation while JISP16 is adopted without renormalization. The MBPT calculations are performed within the Hartree-Fock (HF) bases. The angular momentum coupled scheme is used, which can reduce the computational task. Corrections up to the third order in energy and up to the second order in radius are evaluated. Higher-order corrections in the HF basis are small relative to the leading-order perturbative result. Using the anti-symmetrized Goldstone diagram expansions of the wave function, we directly correct the one-body density for the calculation of the radius, rather than calculate corrections to the occupation propabilities of single-particle orbits as found in other treatments. We compare our results with other methods where available and find good agreement. This supports the conclusion that our methods produce reasonably converged results with these interactions. We also compare our results with experimental data.

21.60.De, 21.30.Fe, 21.10.Dr, 21.10.Ft
thanks: frxu@pku.edu.cn

I Introduction

A fundamental and challenging problem in nuclear structure theory is the calculation of finite nuclei starting from realistic nucleon-nucleon () interactions. The realistic nuclear forces, such as CD-Bonn Machleidt (2001), Nijmegen Stoks et al. (1994), Argonne V18 (AV18) Wiringa et al. (1995), INOY Doleschall (2004) and chiral potential Entem and Machleidt (2003); Machleidt and Entem (2011), contain strong short-range correlations which cause convergence problems in the calculations of nuclear structures. To deal with the strong short-range correlations and speed up the convergence, realistic forces are usually processed by certain renormalizations. A traditional approach is the G-matrix renormalization in the Brueckner-Bethe-Goldstone theory Brueckner (1955a); Goldstone (1957); Bethe et al. (1963) in which all particle ladder diagrams are summed. Recently, a new class of renormalization methods has been developed, including Bogner et al. (2002, 2003), Similarity Renormalization Group (SRG) Bogner et al. (2007), Okubo-Lee-Suzuki Ôkubo (1954); Suzuki and Lee (1980); Suzuki (1982a); Suzuki and Okamoto (1983); Suzuki (1982b); Suzuki and Okamoto (1994) and Unitary Correlation Operator Method (UCOM) Roth et al. (2005, 2010). The renormalizations soften realistic interactions and generate effective Hamiltonians, while all symmetries and observables are preserved in the low-energy domain. The renormalization process also generates effective multi-nucleon interactions (sometimes called ”induced” interactions) that are typically dropped for four or more nucleons interacting simultaneously. We will neglect three-nucleon and higher multi-nucleon interactions both ”bare” and ”induced”. There is another class of “bare” forces which are sufficiently soft that they can be used without renormalization, e.g., the JISP interaction which is obtained by the -matrix inverse scattering technique Shirokov et al. (2004, 2005, 2007). These interactions can often be used directly for nuclear structure calculations.

A renormalized interaction should retain its description of the experimental phase shifts up to a cutoff. At the same time, the renormalized interaction provides better convergence in nuclear structure calculations without involving parameter refitting or additional parameters. The calculations based on realistic forces are called ab initio methods when they retain predictive power and accurate treatment of the first principles of quantum mechanics. There have been several ab initio many-body methods, such as No-Core Shell Model (NCSM) Navrátil et al. (2000a, b, 2009); Caprio et al. (2013); Barrett et al. (2013), Green’s Function Monte Carlo (GFMC) Pieper et al. (2001, 2004); Pervin et al. (2007); Marcucci et al. (2008) and Coupled Cluster (CC) Hagen et al. (2008, 2009, 2010). However, due to the limit of computer capability, the NCSM and GFMC calculations are currently limited to light nuclei (e.g., O), while the CC calculations are limited to nuclei near double closed shells.

While renormalization methods typically address short-range correlations, the Hartree-Fock (HF) approach is used to treat long-range correlations. However, the conventional HF method that takes only one Slater determinant describes the motion of nucleons in the average field of other nucleons and neglects higher-order correlations. For a phenomenological potential, one can adjust parameters to improve the agreement of the HF results with data. For realistic interactions, one needs to go beyond the HF approach to include the intermediate-range correlations which are missing in the lowest order HF approach. The many-body perturbation theory (MBPT) is a powerful tool to include the missing correlations Shavitt and Bartlett (2009); Coraggio et al. (2003); Hasan et al. (2004); Roth et al. (2006). The perturbation method starts from a solvable mean-field problem and derives a correlated perturbed solution. The most well-known perturbation expansions are the Brillouin-Wigner (BW) Brillouin (1932); Wigner (1935) and Rayleigh-Schrödinger (RS) Rayleigh et al. (1894); Schrödinger (1926) methods. MBPT calculations are usually performed with an order-by-order expansion represented in the form of groups of diagrams Shavitt and Bartlett (2009). The diagrams of MBPT proliferate as one goes to higher orders but some techniques, such as those introduced by Bruekner Brueckner (1955b), lead to useful cancellations of entire classes of diagrams. This leads to the linked-diagram theorem which simplifies greatly perturbation calculations up to high orders. Goldstone first proved the theorem valid to all orders in the non-degenerate case Goldstone (1957). Later, the theorem was extended to the degenerate case Brandow (1967); Johnson and Baranger (1971); Kuo et al. (1971); Sandars (2007). The linked-diagram theorem in the degenerate case is often referred to as the folded-diagram method.

Some recent works Coraggio et al. (2003); Hasan et al. (2004); Roth et al. (2006) show that the MBPT corrections to HF can significantly improve calculations which were based on realistic forces. The authors used different renormalization schemes, , OLS and UCOM, and obtained the convergence of low-order MBPT calculations Coraggio et al. (2003); Hasan et al. (2004); Roth et al. (2006). In the present work, we perform similar MBPT calculations with the SRG-renormalized chiral NLO potential Entem and Machleidt (2003); Machleidt and Entem (2011) and the “bare” JISP16 interaction Shirokov et al. (2004, 2005, 2007). We also calculate the MBPT corrections to the nuclear radius with the anti-symmetrized Goldstone (ASG) diagrams of the one-body density (up to the second order). We note that, in Ref. Coraggio et al. (2003), the same ASG diagrams for the corrections to energy were used for the corrections to the radius. In Refs. Hasan et al. (2004); Roth et al. (2006), corrections to the radius were approximated through corrections to occupation probabilities. In order to reduce computational task, we calculate the diagrams in the angular momentum coupling representation. Our MBPT corrections to energy are up to the third order, while our MBPT corrections to the radius are up to the second order.

Ii Theoretical framework

ii.1 The effective Hamiltonian

The intrinsic Hamiltonian of the -nucleon system used in this work reads


where the notation is standard. The first term on the right is the intrinsic kinetic energy, and is the interaction including the Coulomb interaction between the protons. We do not include a three-body interacton. In the present work, two different interactions have been adopted for comparison. One is the chiral potential NLO developed by Entem and Machleidt Entem and Machleidt (2003). Another one is the “bare” interaction JISP16 Shirokov et al. (2004, 2005, 2007).

The NLO potential is renormalized by using the SRG technique to soften the short-range repulsion and short-range tensor components. The SRG method is based on a continuous unitary transformation that suppresses off-diagonal matrix elements and drives the Hamiltonian towards a band-diagonal form Bogner et al. (2007). The process leads to high- and low-momentum parts of the Hamiltonian being decoupled. This implies that the renormalized potential becomes softer and more perturbative than the original one. In principle, the SRG method generates three-body, four-body, etc., effective interactions. We neglect these induced terms for the purposes of examining the similarities and differences of results with NN interactions alone. After the renormalization, the Coulomb interaction between protons is added.

The “bare” JISP16 interaction is obtained by the phase-equivalent transformations of the -matrix inverse scattering potential. The parameters are determined by fitting to not only the scattering data but also the binding energies and spectra of nuclei with Shirokov et al. (2007). In the JISP16 potential, the off-shell freedom is exploited to improve the description of light nuclei by phase-equivalent transformations. Polyzou and Glockle Polyzou and Glöckle (1990) have shown that changing the off-shell properties of the two-body potential is equivalent to adding many-body interactions. Therefore, the phase-equivalent transformation can minimize the need of three-body interactions. The “bare” JISP16 interaction has been used extensively and successfully in configuration interaction calculations of light nuclei Maris and Vary (2013); Shirokov et al. (2014a) and in nuclear matter Shirokov et al. (2014b).

ii.2 Spherical Hartree-Fock formulation

With the effective Hamiltonian established, we first perform the HF calculation and then calculate the MBPT corrections to the HF result. For simplicity of computational effort, we limit our investigations here to the spherical, closed-shell, nuclei He and O. These systems are sufficient to gain initial insights into the convergence rates of the ground-state energy and radius with these realistic interactions.

The spherical symmetry preserves the quantum numbers of the orbital angular momentum (), the total angular momentum () and its projection () for the HF single-particle states. In the spherical harmonic oscillator (HO) basis , the HF single-particle state can be written as


where the labels are standard with and for the radial quantum number of the HO basis and isospin projection, respectively. The HF wave function for the -body nucleus is then represented by an anti-symmetrized Slater determinant constructed with the HF single-particle states. By varying the HF energy expectation value (with respect to the coefficients ), we obtain the HF single-particle eigen equations,


where represents the HF single-particle eigen energies, and designates the matrix elements of the HF single-particle Hamiltonian given by


where and are the matrix elements of the two-body effective Hamiltonian and one-body density, respectively. They can be written




where is the occupation number of the HF single-particle orbit, i.e., (occupied) or (unoccupied).

In practice, we diagonalize the following equation to solve the HF single-particle eigenvalue problem


This is a nonlinear equation with respect to variational coefficients . In the spherical closed shell, the HF single-particle eigenvalues are independent of the magnetic quantum number , which leads to a degeneracy. In this case, we can rewrite the eigenvalues by omitting , i.e., and . Then we can simplify Eq. (7) in the angular momentum coupled representation as followsRoth et al. (2006),


with and one-body density matrix


where is the number of the occupied magnetic subshell, i.e., (occupied) or 0 (unoccupied).

ii.3 Rayleigh-Schrödinger perturbation theory

We can separate the -nucleon Hamiltonian Eq. (1) into a zero-order part and a perturbation ,


The exact solutions of the -nucleon system are


For the zero-order part, we write


If we choose the HF single-particle Hamiltonian Eq. (4) as , the zero-order energy is simply the summation of the single-particle energies up to the Fermi level. In the present work, we only investigate the ground states of closed-shell nuclei. For simplicity, we denote the ground-state energy and wave function by and , respectively, omitting the subscript. For the ground state , we formulate the Rayleigh-Schrödinger perturbation theory (RSPT), as follows,


where is called the resolvent of . Here we use intermediate normalization


Arranging the above expressions according to the perturbation orders of , we have


The first-, second-, third-order corrections are


Similarly, the wave function can be written in the perturbation scheme






for the first- and second-order corrections to the wave function, respectively. We can use the diagrammatic approach to describe various terms in RSPT. The ASG diagrams are the most commonly-used method of the diagrammatic representation.

ii.4 Diagrammatic expansion for Rayleigh-Schrödinger perturbation theory in the Hartree-Fock basis

If we choose the HF Hamiltonian as an auxiliary zero-order one-body Hamiltonian , many of the ASG diagrams are cancelled Shavitt and Bartlett (2009). Only a small number of low-order ASG diagrams for RSPT remain. In this subsection, we give the remaining AGS diagrams for the energy and wave function written in the standard perturbation theory Kuo and Osnes (1990). We consider corrections up to third order for the energy and second order for the wave function. To evaluate other observables that can be expressed by one-body operators, we calculate the corrections up to second order for the one-body density. It has been shown that the corrections up to third order for the energy in the HF basis give well-converged results for soft interactions Tichai et al. (2016). Spherical HF (SHF) produces degenerate single-particle states, so we can evaluate the vacuum-to-vacuum linked diagrams in angular momentum coupled representation Kuo et al. (1981) which is computationally efficient.

Figure 1: The first-, second-, and third-order ASG diagrams of energy corrections in the RS expansion Coraggio et al. (2003).

Fig. 1 displays the ASG diagrams corresponding to the first-, second- and third-order corrections to the energy in RSPT. The vertices, i.e., the dashed lines, represent in Eq. (1). The diagrams (a) and (b) are for and , respectively, while the diagrams (c), (d) and (e) sum up for . The zero-order energy is the simple summation of the HF single-particle energies up to the Fermi level, i.e., , where represents the HF single-particle energy. The summation of the and gives the HF energy, i.e., , since the initial Hamiltonian is entirely expressed in relative coordinates Bozzolo et al. (1988); Hasan et al. (2004).

ii.4.1 Corrections to the one-body density

MBPT corrections to the wave function bring configuration mixing. The convergence can be discussed in order-by-order perturbation calculations. Any observable that is expressed by one-body operators can be calculated by using the One-Body Density Matrix (OBDM). By definition, the local one-body density operator in an -body Hilbert space is written as Cockrell et al. (2012)


where is the unit vector in the direction , and is the spherical harmonic function.

We can write the density operator in the second quantization representation in the HO basis as




The ’s are the radial components of the HO wave function. We use the Condon-Shortley convention for the Clebsch-Gordan coefficients. Since we are dealing with a spherically symmetric system (K=0), we can obtain a simple form,


By introducing the normally-ordered product relative to the SHF ground state , the local one-body density operator can be written as


where gives the HF density, while brings corrections to the density. is the density matrix elements , and indicates the normally-ordered product of the creation and annihilation operators. It is required that all annihilation and creation operators which take to zero when acting on it are to the right of all other operators which do not take to zero. The expectation value of the density is obtained with the corrected wave function through Eq. (29). In the present work, we consider the first- and second-order wave function corrections.

Figure 2: ASG diagrams for the first- and second-order corrections to the wave function Shavitt and Bartlett (2009). The panel (a) is for the first-order correction, while (b) (c) … (i) are for the second-order correction.

The ASG diagrams for the first- and second-order corrections to the wave function Shavitt and Bartlett (2009) are displayed in Fig. 2. The first-order wave function diagram, i.e., panel (a) in Fig. 2, produces the second-order correction to the density. While diagrams (b) and (c) of the second-order wave function correction produce second-order corrections to the density, other diagrams of the second-order wave function correction contribute to higher-order corrections to the density. The first- and second-order wave function corrections which correct the density up to the second order can be written as


The total wave function that corrects the density up to the second order is


Then, the corrected density is written as


where , and . They are displayed using the language of the diagram in Fig. 3. Dashed lines with cross contribute to the reduced matrix elements .

The detailed formulae of the density correction terms in the angular momentum coupled scheme are written as


where is Wigner 6-j symbol. The letters indicate occupied single-particle levels in (i.e., hole states), the letters for unoccupied levels (i.e., particle states). or is the energy of particle or hole state, respectively. States or includes the quantum numbers of the orbital angular momentum , total angular momentum , isospin projection quantum number , and additional quantum number , i.e., or . We define an anti-symmetrized two-particle state (unnormalized) coupled to a good angular momentum with a projection ,


(a)   (b)   (c)   (d)  

Figure 3: ASG diagrams for the second-order corrections to the density.

ii.4.2 Root-mean-square radii

The root-mean-square (rms) radius is an important global indicator for the change of the density distribution arising from correlations beyond HF. The squares of the rms radii for point-like proton, neutron and nucleon (matter) distributions are the averaged values of the operators Kamuntavičius (1997), respectively,


with the c.m. position . The charge radius obtained from the point-proton radius using the standard expression Ekström et al. (2015)


where , , fm. The point-proton or point-neutron rms radius operator is a two-body operator. The squares of the rms radii can be calculated either from the translational invariant local density or directly using the two-body operators [ i.e., Eqs. (40), (41) and (42) ]. Since we adopt MBPT with intermediate normalization [ i.e., Eqs. (17) ], the perturbed wave function is unnormalized. In the present work, we use the one-body local density to calculate the radius, as


The wave function is written in the laboratory HO coordinate, starting from an anti-symmetrized Slater determinant which contains the component of the center-of-mass (c.m.) motion. Consequently, the local one-body density calculated with the wave function includes contribution from the c.m. motion. The c.m. correction to the radius can be approximated as follows. Eq. (42) gives


If the cross term is neglected, we have


Similarly for the proton radius,


This gives an approximate c.m. correction to the point-proton rms radius,


where is the point-proton rms radius calculated by Eq (44). Then the rms radius of the point-proton distribution is obtained by


Iii Calculations and discussions

In this section, we apply the method outlined in Section II to two light closed-shell nuclei, He and O. The SRG-softened chiral NLO and the “bare” JISP16 interactions are adopted for the effective Hamiltonians.

iii.1 Calculations with chiral NLO interaction

The SHF is carried out within the HO basis. The HO basis is truncated by a cutoff according to the number , where indicates how many major HO shells are included in the truncation. After the SHF calculation, the MBPT corrections are calculated in the SHF basis. In the present calculations, the basis spaces employed take =7, 9, 11 and 13. We verify that such a truncation is sufficient for the converged calculations of the ground state energies for these magic nuclei He and O.

Figure 4: HF-MBPT calculations of ground-state energy through third order as a function of oscillator parameter with the chiral NLO potential Entem and Machleidt (2003); Machleidt and Entem (2011) renormalized by SRG at different softening parameters . The dashed line represents the experimental ground-state energy.
Figure 5: Point-proton rms radius of He as a function of oscillator parameter with different . The chiral NLO potential Entem and Machleidt (2003); Machleidt and Entem (2011) is softened by the SRG method.
SRG flow parameter (fm)
1.5 2.0 2.5 3.0
Expt. Audi et al. (2012) -28.296 -28.296 -28.296 -28.296
NCSM  Jurgenson et al. (2013) -28.20 -28.41 -27.43 -26.80
SHF -25.754 -21.864 -15.854 -10.278
PT2 -1.788 -5.088 -9.652 -13.783
PT3 -0.391 -0.899 -1.523 -1.953
SHF+PT2+PT3 -27.933 -27.850 -27.029 -26.013
Table 1: Ground-state energy (in MeV) of He, analyzed in order-by-order HF-MBPT calculations with NLO softened at different SRG-softening parameter values (). PT2 and PT3 represent the second- and third-order corrections to energy, respectively. We take and MeV.
SRG flow parameter (fm)
1.5 2.0 2.5 3.0
Expt. 1.477 1.477 1.477 1.477
SHF 1.677 1.652 1.714 1.816
PT2 0.007 0.001 -0.021 -0.065
-0.226 -0.222 -0.227 -0.235
SHF+PT2+ 1.458 1.431 1.466 1.516
Table 2: Point-proton rms radius (in fm) of He in the HF-MBPT calculations with NLO softened at different SRG-softening parameter values. PT2 designates the second-order correction to the radius. and MeV are taken. The experimental point-proton rms radius is obtained using Eq. (43) with the experimental charge radius taken from Angeli and Marinova (2013).

Fig. 4 shows the MBPT calculated ground-state energy of He. The calculations were done with the chiral NLO interaction which was renormalized by SRG. We see that good convergence of the calculated energy by virtue of independence from the oscillator parameter and is obtained at least for the truncations and 13. We note that the dependence on the parameter displays behavior similar to NCSM calculations Bogner et al. (2008); Jurgenson et al. (2013). The softening parameter fm seems to be insufficient to produce an interaction soft enough for good convergence in MBPT. Jurgenson et al., have investigated the SRG evolution with the softening parameter in He at MeV Jurgenson et al. (2009, 2013). They found that fm can reasonably reproduce the experimental He ground-state energy with the -only interaction (without requiring a three-body force).

Fig. 5 shows the radius calculations at different with fm. Tables 1 and  2 give the details of the HF-MBPT calculations with different values. We see that both second- and third-order corrections to energy decrease with decreasing . This is easily understood because MBPT mainly treats intermediate-range correlations and these correlations are weakened with decreasing . With sufficiently small