Variational approximations to exact solutions in shell-model valence spaces: calcium isotopes in the pf-shell

TL;DR

This paper evaluates the accuracy of mean-field and beyond-mean-field methods in shell-model calculations for calcium isotopes, demonstrating near-exact results when advanced techniques and collective coordinates are used.

Contribution

It introduces a comprehensive comparison of self-consistent mean-field and advanced projection methods against exact shell-model solutions in the pf-shell.

Findings

High accuracy of beyond-mean-field methods with particle-number and angular-momentum restoration

Explicit inclusion of quadrupole and neutron-neutron pairing improves results

Excellent agreement with exact diagonalization for calcium isotopes

Abstract

We study the performance of self-consistent mean-field and beyond-mean-field approximations in shell-model valence spaces. In particular, Hartree-Fock-Bogolyubov, particle-number variation after projection and projected generator coordinate methods are applied to obtain ground-state and excitation energies for even-even and odd-even Calcium isotopes in the pf-shell. The standard (and non-trivial) KB3G nuclear effective interaction has been used. The comparison with the exact solutions -- provided by the full diagonalization of the Hamiltonian -- shows an outstanding agreement when particle-number and angular-momentum restorations are performed and both quadrupole and neutron-neutron pairing degrees of freedom are explicitly explored as collective coordinates.

Figures

Tables

Table 1. Table 1: Different PGCM approximation schemes used to compute the structure of calciums isotopes.

Label	Energy minimization	Collective coordinates
PGCM1	HFB	$({\beta }_2,\gamma )$
PGCM2	PNVAP	$({\beta }_2,\gamma )$
PGCM3	PNVAP	$({\beta }_2,\gamma ,{\delta }_{nn})$

Table 2. Table 2: Reduced transition probabilities, B ​ ( E ​ 2 ) 𝐵 𝐸 2 B(E2) , and electric spectroscopic quadrupole moments, Q 𝑄 Q , calculated with different methods using proton (1.5) and neutron (0.5) effective charges. The B ​ ( E ​ 2 ) 𝐵 𝐸 2 B(E2) are given in units of e 2 superscript 𝑒 2 e^{2} fm 4 whereas the Q 𝑄 Q are expressed in units of e 𝑒 e fm 2 .

	PGCM1	PGCM2	PGCM3	ISM
$B(E2:2_1^{+}\to 0_1^{+})$	12.7	12.7	12.7	11.5
$B(E2:2_1^{+}\to 0_2^{+})$	0.8	0.9	1.0	1.0
$B(E2:2_2^{+}\to 0_1^{+})$	0.0	0.0	0.0	0.0
$B(E2:2_2^{+}\to 0_2^{+})$	30.4	24.5	23.0	21.6
$B(E2:4_1^{+}\to 2_1^{+})$	5.6	5.7	2.5	2.0
$Q\left(2_1^{+}\right)$	+5.0	+4.8	+4.4	+4.1
$Q\left(2_2^{+}\right)$	-11.4	-9.8	-9.3	-8.6
$Q\left(4_1^{+}\right)$	+10.8	+9.1	+8.3	+7.5

Table 3. Table 3: Occupation numbers computed for the two lowest 0 + superscript 0 0^{+} , 2 + superscript 2 2^{+} and 4 + superscript 4 4^{+} states in 48 Ca using the ISM method and the different PGCM schemes.

$J^{\pi \sigma }$	orbit	PGCM1	PGCM2	PGCM3	ISM
$0_1^{+}$	$f_{7/2}$	7.84	7.81	7.81	7.80
	$p_{3/2}$	0.10	0.08	0.08	0.07
	$f_{5/2}$	0.04	0.09	0.09	0.11
	$p_{1/2}$	0.02	0.02	0.02	0.02
$0_2^{+}$	$f_{7/2}$	5.89	5.84	5.84	5.82
	$p_{3/2}$	1.44	1.64	1.75	1.80
	$f_{5/2}$	0.15	0.18	0.15	0.17
	$p_{1/2}$	0.52	0.34	0.26	0.21
$2_1^{+}$	$f_{7/2}$	6.80	6.77	6.77	6.77
	$p_{3/2}$	1.07	1.07	1.05	1.08
	$f_{5/2}$	0.09	0.12	0.14	0.11
	$p_{1/2}$	0.04	0.04	0.04	0.04
$2_2^{+}$	$f_{7/2}$	5.86	5.85	5.87	5.85
	$p_{3/2}$	1.44	1.63	1.71	1.78
	$f_{5/2}$	0.15	0.14	0.12	0.13
	$p_{1/2}$	0.55	0.38	0.30	0.24
$4_1^{+}$	$f_{7/2}$	6.66	6.59	6.83	6.83
	$p_{3/2}$	1.15	1.12	1.06	1.03
	$f_{5/2}$	0.13	0.23	0.08	0.11
	$p_{1/2}$	0.06	0.06	0.03	0.03
$4_2^{+}$	$f_{7/2}$	6.63	6.60	6.79	6.82
	$p_{3/2}$	0.43	0.37	0.17	0.09
	$f_{5/2}$	0.11	0.15	0.10	0.11
	$p_{1/2}$	0.83	0.88	0.94	0.98

Table 4. Table 4: Dimensions of the Hamiltonian matrices for ISM (value on the left) and PGCM 2 (value on the right) in even (top) and odd (bottom) calcium isotopes in the p ​ f 𝑝 𝑓 pf -shell.

$J$	42 or 58	44 or 56	46 or 54	48 or 52	50
0	4 $|$ 4	28 $|$ 7	137 $|$ 7	347 $|$ 9	468 $|$ 9
2	8 $|$ 7	94 $|$ 16	512 $|$ 14	1390 $|$ 20	1935 $|$ 20
4	6 $|$ 6	99 $|$ 21	615 $|$ 20	1755 $|$ 28	2468 $|$ 13
6	2 $|$ 2	59 $|$ 17	462 $|$ 17	1426 $|$ 26	2051 $|$ 8
$J$	43 or 57	45 or 55	47 or 53	49 or 51
1/2	12 $|$ 9	107 $|$ 10	415 $|$ 12	790 $|$ 12
3/2	25 $|$ 18	198 $|$ 21	764 $|$ 26	1484 $|$ 15
5/2	28 $|$ 21	253 $|$ 25	1005 $|$ 37	1965 $|$ 30
7/2	27 $|$ 21	271 $|$ 23	1121 $|$ 24	2215 $|$ 32
9/2	23 $|$ 22	252 $|$ 41	1091 $|$ 45	2214 $|$ 38
11/2	16 $|$ 16	211 $|$ 34	974 $|$ 45	2017 $|$ 35
13/2	8 $|$ 8	153 $|$ 33	783 $|$ 39	1669 $|$ 28
15/2	5 $|$ 5	105 $|$ 25	577 $|$ 33	1284 $|$ 21

Equations

$$\hat{H}=\sum _{ab}{t}_{ab}{c}_a^{†}{c}_b+\frac{1}{(2!{)}^2}\sum _{abcd}{\bar{v}}_{abcd}{c}_a^{†}{c}_b^{†}{c}_d{c}_c,$$

$$|{\Psi }_{\sigma }^{JMNZ\pi }⟩=\sum _{qK}{f}_{\sigma ;qK}^{JMNZ\pi }{P}_{MK}^J{P}^N{P}^Z{P}^{\pi }|\Phi (q)⟩,$$

$$\sum _{q^{\prime }{K}^{\prime }}\left({\mathcal{H}}_{qK,q^{\prime }{K}^{\prime }}^{\Gamma }-E_{\sigma }^{\Gamma }{\mathcal{N}}_{qK,q^{\prime }{K}^{\prime }}^{\Gamma }\right)f_{\sigma ;q^{\prime }{K}^{\prime }}^{\Gamma }=0,$$

$${\mathcal{H}}_{qK,q^{\prime }{K}^{\prime }}^{\Gamma }$$

$${\mathcal{N}}_{qK;q^{\prime }{K}^{\prime }}^{\Gamma }$$

$$\sum _{q^{\prime }{K}^{\prime }}{\mathcal{N}}_{qK,q^{\prime }{K}^{\prime }}^{\Gamma }{u}_{\lambda ;q^{\prime }{K}^{\prime }}^{\Gamma }=n_{\lambda }^{\Gamma }{u}_{\lambda ;qK}^{\Gamma },$$

$$|{\Lambda }_{\lambda }^{\Gamma }⟩=\sum _{qK}\frac{u_{\lambda ;qK}^{\Gamma }}{\sqrt{n_{\lambda }^{\Gamma }}}{P}_{MK}^J{P}^N{P}^Z{P}^{\pi }|\Phi (q)⟩,$$

$$|{\Psi }_{\sigma }^{\Gamma }⟩=\sum _{\lambda }{G}_{\sigma ;\lambda }^{\Gamma }|{\Lambda }_{\lambda }^{\Gamma }⟩,$$

$$\sum _{{\lambda }^{\prime }}⟨{\Lambda }_{\lambda }^{\Gamma }|\hat{H}|{\Lambda }_{{\lambda }^{\prime }}^{\Gamma }⟩G_{\sigma ;{\lambda }^{\prime }}^{\Gamma }=E_{\sigma }^{\Gamma }{G}_{\sigma ;\lambda }^{\Gamma }.$$

$$\left(\begin{array}{c}\beta (q)\\ {\beta }^{†}(q)\end{array}\right)={\mathcal{W}}^{†}(q)\left(\begin{array}{c}c\\ c^{†}\end{array}\right),$$

$$\mathcal{W}(q)\equiv \left(\begin{array}{cc}U(q)& V^{\ast }(q)\\ V(q)& U^{\ast }(q)\end{array}\right).$$

$${\hat{Q}}_{2\mu }$$

$$q_{20}$$

$$q_{21}$$

$$q_{22}$$

$$q_{20}$$

$$q_{22}$$

$$\begin{array}{rl}{\left[{\hat{P}}^{†}\right]}_{M_J{M}_T}^{JT}& =\sum _{\breve{a}}{\left[{\hat{P}}_{\breve{a}}^{†}\right]}_{M_J{M}_T}^{JT}\\ & =\frac{1}{\sqrt{2}}\sum _{\breve{a}}\sqrt{2j_a+1}{\left[c_{\breve{a}}^{†}{c}_{\breve{a}}^{†}\right]}_{M_J{M}_T}^{JT}\end{array}$$

$$\begin{array}{rl}{\left[c_{\breve{a}}^{†}{c}_{\breve{b}}^{†}\right]}_{M_J{M}_T}^{JT}=& \frac{1}{\sqrt{1+{\delta }_{j_a{j}_b}}}\sum _{\begin{array}{c}{m}_{j_a}{m}_{j_b}\\ m_{t_a}{m}_{t_b}\end{array}}⟨j_a{m}_{j_a}{j}_b{m}_{j_b}|J{M}_J⟩\\ & \times ⟨\frac{1}{2}{m}_{t_a}\frac{1}{2}{m}_{t_b}|T{M}_T⟩c_a^{†}{c}_b^{†}.\end{array}$$

$${\delta }_{nn}=\frac{1}{2}⟨\Phi (q)|{\left[\hat{P}\right]}_{01}^{01}+{\left[{\hat{P}}^{†}\right]}_{01}^{01}|\Phi (q)⟩.$$

$$\left|F^{\Gamma}_{\sigma}(q)\right|^{2}=\bigg{|}\sum_{K\lambda}G^{\Gamma}_{\sigma;\lambda}u^{\Gamma}_{\lambda;qK}\bigg{|}^{2}.$$

$$n_{\check{a}}^{\Gamma \sigma }=⟨{\Psi }_{\sigma }^{\Gamma }|{\hat{n}}_{\check{a}}|{\Psi }_{\sigma }^{\Gamma }⟩,$$

$${\hat{n}}_{\check{a}}=\sum _b{c}_b^{†}{c}_b{\delta }_{\breve{b}\breve{a}}{\delta }_{m_{t_b}{m}_t},$$

Full text

I Introduction

The study of the structure of the atomic nuclei from a microscopic point of view is a very complicated task mainly due to two interrelated problems. First, the interaction among the composite nucleons is rather intricate and for a long time could not be reliably derived from first principles. The recent advent of Effective Field Theories in nuclear physics Epelbaum et al. (2009); Hammer et al. (2019) has considerably improved the situation. The other complication arises from the number of particles in a nucleus that is usually too large to solve exactly the problem but not large enough to describe the system using statistical mechanics. Therefore, one of the keys to the theoretical description of the atomic nucleus is the design and implementation of efficient quantum many-body techniques.

Over the last decade, ab initio methods, which aim at solving as exactly as possible, and in a systemically improvable manner, the many-body Schrödinger equation, have made great progress Somà et al. (2013); Hagen et al. (2014); Hergert et al. (2016); Tichai et al. (2018). In certain specific cases, they can even reach nuclei with $A\sim 100$ nucleons Morris et al. (2018). Their accuracy in the description of experimental data, however, is not yet on par with what is obtained by phenomenological models such as the Interacting Shell Model (ISM) Caurier et al. (2005); Otsuka et al. (2001) or the Energy Density Functional (EDF) method Bender et al. (2003); Nikšić et al. (2011); Egido (2016); Robledo et al. (2018).

Within the ISM framework, the problem with $A$ nucleons is reduced to a system with $A^{\prime}=A-A_{core}$ particles that are interacting in a restricted valence space. The resulting Hamiltonian matrix is then directly diagonalized to obtain the exact solutions of the reduced problem. The valence space is normally made of a few one-body spherical harmonic oscillator orbits and $A_{core}$ is normally a given double-magic nucleus that defines an inert core. To this day, the ISM is the most successful method to compute low-lying nuclear spectra of nuclei not very far away from shell closures.

On the other hand, EDF methods consider all the particles in the system but the many-body problem is solved through a variational approximation. The most advanced EDF methods, the so-called Multi Reference EDF (MREDF) approaches, combine the Hartree-Fock-Bogolyubov (HFB) methodology with symmetry restoration and configuration mixing techniques Duguet (2014). Calculations using Skyrme, Gogny111In Gogny implementations, MREDF approaches are also known as symmetry conserving configuration mixing (SCCM) methods Egido (2016); Robledo et al. (2018) and covariant functionals have been widely applied to describe various nuclear observables in the whole nuclear chart Bender et al. (2003); Nikšić et al. (2011); Egido (2016); Robledo et al. (2018).

The broadness of the applicability of the variational approaches is one of the main positive aspects of these methods. However, we cannot judge in general whether the results provided by these methods are close or far away from the exact solutions. This indetermination leads to an uncertainty in the theoretical results that cannot be easily assessed. In this work, we benchmark several variational approximation schemes against the exact solutions using a shell-model valence space and Hamiltonian as a non-trivial playground. In particular, we study the even and odd calcium isotopes in the $pf$-shell using the standard KB3G nuclear interaction. The different approaches cover from the simplest HFB method to the most sophisticated projected generator coordinate method (PGCM) that includes particle-number and angular-momentum restorations as well as the configuration mixing along the quadrupole and pairing degrees of freedom.

Similar implementations to the present approach have been recently used to evaluate neutrinoless double-beta decay nuclear matrix elements Jiao et al. (2017); Yao et al. (2018); Jiao et al. (2018). Some of them consider a limited number of collective degrees of freedom as generating coordinates and none of them implements the particle-number variation after projection (PNVAP) method. Additionally, only even-even nuclei were computed so far with valence-space Hamiltonian-based GCM. A study of HFB and PNVAP potential energy surfaces along the axial quadrupole degree of freedom in $sd$- and $pf$-shell even-even nuclei has been also reported in Ref. Maqbool et al. (2011) but neither angular momentum restoration nor configuration mixing were performed in that work. Moreover, Hartree-Fock (HF) total energy surfaces along quadrupole degrees of freedom have been used to interpret ISM Hadyńska-Klȩk et al. (2016); Bounthong (2016) and Monte Carlo Shell Model (MCSM) Tsunoda et al. (2014); Togashi et al. (2018) calculations. prominently Also, simpler Hamiltonians like Pairing-plus-quadrupole (PPQ) interactions have been also used in the past within the projected shell model (PSM) framework Hara and Sun (1995); Sun (2016) and with GCM approximations including two quasi-particle excitations explicitly Chen and Egido (2016, 2017). Our goal is to go beyond these simplistic Hamiltonians although, for the moment, some single-particle degrees of freedom included in the above methods will be omitted. Finally, we note that the idea to build sophisticated variational approximations is at the core of the famous variation after mean-field projection in realistic model spaces (VAMPIR) Schmid and Grummer (1987); Schmid (2004) approach that includes the variation after projection (VAP) onto particle number and angular momentum simultaneously. The angular-momentum variation after projection although highly desirable is numerically very intensive in the general case, i.e. considering many states without any symmetry restrictions and projected onto other quantum numbers, and therefore is not considered here.

The paper is organized as follows. First, in Sec. II, we define the valence space, the Hamiltonian, and the different variational approaches used in this work. Then, in Sec. III, we show how these methods are implemented using the nucleus 48Ca as an example. From there, we extend the study to the rest of even and odd calcium isotopes. Finally, a summary and outlook is given in Sec. IV.

II Theoretical framework

II.1 Basic principles

Let us first consider a model space spanned by a set of harmonic oscillator single-particle states $\left\{\phi_{a}\right\}$, associated with the set of creation and annihilation operators $\left\{c_{a};c^{\dagger}_{a}\right\}$, which are characterized by their principal quantum number $n_{a}$, orbital angular momentum $l_{a}$, spin $s_{a}=1/2$, total angular momtum $j_{a}$ and its third component $m_{j_{a}}$, and isospin $t_{a}=1/2$ and its third component $m_{t_{a}}$.222We use as a convention $m_{t_{a}}=-1/2$ for proton single-particle states and $m_{t_{a}}=+1/2$ for neutron single-particle states. For the sake of clarity, we use the shorthand notation $a\equiv(n_{a},l_{a},s_{a},j_{a},m_{j_{a}},t_{a},m_{t_{a}})$.

Given this model space, the Hamiltonian333The inclusion of three- or higher-body forces would not pose any conceptual problem. of the system reads

[TABLE]

where $t_{ab}$ are one-body matrix elements and $\bar{v}_{abcd}$ are (antisymmetrized) two-body matrix elements.

In the present work, the model space under consideration is restricted to the valence space made of the $pf$-shell spherical orbits ($0f_{7/2}$, $1p_{3/2}$, $1p_{1/2}$ and $0f_{5/2}$). As for the Hamiltonian, we use the matrix elements of the KB3G interaction Poves et al. (2001), which in particular provides us with the single-particle energies $t_{ab}=\epsilon_{a}\delta_{ab}$.

In such a small model space, it is possible to solve exactly the many-body Schrödinger equation for any given numbers of valence protons, $Z$, and neutrons, $N$, by performing a direct diagonalization of the Hamiltonian matrix in the many-body basis made of all $A^{\prime}$-particles Slater determinant, with $A^{\prime}=Z+N$. This is the standard ISM framework and, in this work, the ISM solutions are obtained using the code ANTOINE developed by the Strasbourg-Madrid collaboration Caurier et al. (2005).

On the other hand, the number of configurations and, correspondingly, the dimensions of the Hamiltonian matrices increase combinatorially with the number of active particles and the size of the valence space. This fact forbids the calculation of exact solutions in the general case and, therefore, approximate methods are in order. One of the most widely used is the PGCM Ring and Schuck (1980); Bender et al. (2003); Robledo et al. (2018) within which the nuclear many-body wave functions are taken to be linear combinations of the form

[TABLE]

where $\left\{|\Phi(q)\rangle\right\}$ are the so-called intrinsic states, assumed here to be Bogolyubov quasi-particle states, $P^{N(Z)}$ is the projector onto good number of neutrons (protons), $P^{J}_{MK}$ is the angular momentum projection operator, $P^{\pi}$ is the parity projector. The precise form of the projection operators can be found somewhere else Bender et al. (2003); Ring and Schuck (1980); Bally (2014). For the sake of clarity, we regroup the set of quantum numbers under the label $\Gamma\equiv(JMNZ\pi)$. The coefficients $f^{\Gamma}_{\sigma;qK}$ are obtained by solving the Hill-Wheeler-Griffin (HWG) equations Ring and Schuck (1980)

[TABLE]

where we have defined the Hamiltonian and norm overlaps matrices

[TABLE]

The associated generalized eigenvalue problem can be transformed into a normal eigenvalue problem by removing the linear dependencies of the original set of projected states $\left\{P^{J}_{MK}P^{N}P^{Z}P^{\pi}|\Phi(q)\rangle\right\}$. Hence, the nuclear states given in Eq. (2) are expressed in the so-called natural basis that is obtained by finding first the eigenvalues and eigenvectors of the norm overlap matrix

[TABLE]

then, the orthonormal states are defined as

[TABLE]

where we set a threshold for the smallest norm overlap matrix eigenvalue, $n^{\Gamma}_{\lambda}>\epsilon$. The nuclear states, Eq. (2), can now be written as

[TABLE]

and the HWG equations

[TABLE]

Finally, the energies $E^{\Gamma}_{\sigma}$ computed with the PGCM method are approximations to the exact eigenvalues of $\hat{H}$. In addition, the coefficients $G^{\Gamma}_{\sigma;\lambda}$ can be used to evaluate the collective wave functions, occupation numbers, transition probabilities, electromagnetic moments, etc. Bender et al. (2003); Robledo et al. (2018).

We note in passing that our Hamiltonian-based method is free from the problems of self-interaction and self-pairing appearing in MR EDF calculations Anguiano et al. (2001); Lacroix et al. (2009); Bender et al. (2009); Duguet et al. (2009), which in particular implies that all the direct, exchange and pairing terms are included. Also, to compute the norm overlap between two quasi-particle vacua, we use the Pfaffian method Robledo (2009); Avez and Bender (2012) with the numerical routines published in Ref. Wimmer (2012). Other methods have been recently proposed Bally and Duguet (2018); Mizusaki et al. (2018).

II.2 Generation of the intrinsic states

The quality of the PGCM approach relies heavily on the richness of the underlying set of selected intrinsic wave functions $\left\{|\Phi(q)\rangle\right\}$. In the present case, all of them are Bogolyubov quasi-particle wave functions, i.e. they are vacua for a set of quasi-particle operators $\left\{\beta_{a}(q);\beta^{\dagger}_{a}(q)\right\}$ defined through unitary linear Bogolyubov transformations

[TABLE]

where

[TABLE]

The matrices $U(q)$ and $V(q)$ are variational parameters that are obtained by solving either HFB or particle-number variation after projection (PNVAP) equations (see below). There are several aspects that characterize a set of intrinsic states. First, the above transformations can be restricted and the quasi-particles would fulfill certain self-consistent symmetries. In the present work, we allow for general HFB transformations without any restriction except for the fact that the matrices $U(q)$ and $V(q)$ are kept real. In particular, angular-momentum, particle-number and parity symmetry breaking, as well as proton-neutron mixing are allowed. The latter is key to include proton-neutron pairing in the wave functions and it is normally ignored in most of the EDF calculations. However, it will be irrelevant in the study of the calcium isotopic chain within the $pf$-shell since only neutrons are involved in the valence space. Moreover, all the single-particle orbits in the $pf$-shell have negative parity and, therefore, only positive (negative) parity states can be obtained for even (odd) nuclei. As a consequence, parity is not broken by the HFB transformations and the parity projection is not needed.

Second, we distinguish between HFB and VAP approaches Ring and Schuck (1980) depending on the definition of the energy minimized in the variational process. In this study, we perform either HFB or PNVAP calculations. The PNVAP method is variationally better for ground state energies, provides larger pairing correlations and is free from spurious phase transitions in the pairing channel that occur in HFB approaches Anguiano et al. (2001); Rodríguez et al. (2005); Rodríguez and Egido (2007). Nevertheless, we will compare the results obtained by PGCM schemes using both HFB and PNVAP underlying states. In both cases, the gradient method Robledo and Bertsch (2011); Egido et al. (1995) complemented by a momentum term Ryssens et al. (2019) is used to solve the non-linear equations at hand.

The intrinsic wave functions of nuclei with an odd number of particles are obtained by the blocking method to ensure their correct (odd) number-parity symmetry Ring and Schuck (1980). The solution of the variational equations could depend on the initial blocked state (the seed wave function needed to start the calculation) Bally et al. (2014); Bally (2014). This is particularly important in HFB calculations where local minima that correspond to ground or excited states could be obtained by the gradient method. In the present work, we try several blocked initial configurations. Then, we take the configuration that gives the lowest energy after convergence as the actual solution. On the other hand, PNVAP calculations are much less dependent on the choice of the seed state and the gradient method leads normally to the absolute minimum. The choice to select only the lowest self-consistent quasi-particle state for given set of constraint may warrant, a priori, a good variational exploration of excitated states. But as we will see later on, we obtain an excellent agreement with the exact results also for those states.

The third important aspect about the definition of the set of intrinsic states is the selection of the collective coordinates, $q$. Hence, the variational equations are solved with constraints, i.e., the energy that is minimized is modified with Lagrange multipliers. These terms ensure that the conditions $\langle\Phi(q)|\hat{Q}|\Phi(q)\rangle=q$ are fulfilled ($\hat{Q}$ are the operators associated to the collective coordinates). The most important collective degrees of freedom present in realistic nuclear interactions are the lowest multipole deformations, particularly the quadrupole deformations, and pairing Dufour and Zuker (1996). The former can be explored by imposing constraints in the values of the mass quadrupole operators

[TABLE]

where $C=\frac{3r_{0}^{2}A^{5/3}}{8\pi}$, $r_{0}=1.2$ fm and $A$ is the total mass number (including core and valence space particles). Additionally, electromagnetic transitions and moments are calculated with the effective charge $e_{p}$ for protons and $e_{n}$ for neutrons. In the $pf$-shell, we choose the standard values $1.5$ and $0.5$ for protons and neutrons, respectively Caurier et al. (2005); Otsuka et al. (2001).

Pairing degrees of freedom can be explored by constraining the expectation value of a pair creation operator with respect to the intrinsic states. In this work, we use an operator that couples pairs within a given orbit $\breve{a}\equiv(n_{a},l_{a},j_{a},s_{a},t_{a})$ to a good total angular momentum $J$ and total isospin $T$ Talmi (1993); Hinohara and Engel (2014)

[TABLE]

where the creation operators are $JT$-coupled according to

[TABLE]

Both isoscalar ($T=0$, $J=1$) $pn$-pairing and isovector ($T=1$, $J=0$) $pp$-, $nn$-, and $pn$-pairing can be explored with these operators. In this work, we only study the $nn$-pairing channel ($T=1$, $M_{T}=1$, $J=0$, $M_{J}=0$) because only neutrons are present in the calcium chain in the $pf$-shell, i.e., the intrinsic wave functions can be constrained to

[TABLE]

In Table 1, we summarize and label the different PGCM schemes that are examined in the present work depending on the type of energy minimization scheme used (HFB or PNVAP) and the collective coordinates explored. In all cases, particle-number (proton and neutron) and angular-momentum (three Euler angles) projections were performed. The number of integration points taken to discretize the integrals over the gauge and Euler angles was large enough as to ensure a full convergence in the nominal expectation values of the particle-number and angular-momentum operators computed with the GCM wave functions (Eq. (2)). The PGCM calculations were performed using the newly developed software TAURUS Bally et al. (2019).

III Results

III.1 48Ca

In this section, we illustrate our methodology taking the nucleus 48Ca as an example. It is noteworthy that, as it will be demonstrated below, it represents one of the most difficult cases for our model.

As mentioned in the previous section, the first step in our method is the construction of a set of quasi-particle states through a series of constrained HFB/PNVAP calculations. It is important to point out that in a restricted valence space the range of admissible values for the constraints is much more limited than in a no-core implementation. Indeed, working with a handful of particles and single-particle states, it is not possible to build a many-body wave function that satisfies any arbitrary values of the constraints. For example, the largest $\beta_{2}$ value reachable in the model space is quite small compared to the values used in traditional EDF calculations. In the present work, the bounds of the constraints are determined heuristically.

The total energy as a function of the quadrupole degrees of freedom, $(\beta_{2},\gamma)$, is represented in Fig. 1. On the top panels, the HFB (Fig. 1(a)) and PNVAP (Fig. 1(b)) total energy surfaces (TES) are shown. As expected in this doubly-magic nucleus, in both cases the absolute minimum is located at the spherical configuration and the energy rises quickly with $\beta_{2}$ and is almost independent of $\gamma$; also we observe that the PNVAP surface is slightly softer. At the spherical point, the pairing collapses in the HFB calculation, which is thus reduced to a plain HF minimization. By contrast, including pairing correlations through the PNVAP permits to gain about 750 keV extra binding energy. However, the main differences between the two approaches are found whenever the particle-number and angular momentum projection ($J=0$) are performed on their respective intrinsic states. In the PNVAP+PNAMP case (Fig. 1(d)), we obtain an interval in $\beta_{2}\in[0.0,0.1]$ where the TES is almost flat independently of the $\gamma$ degree of freedom. The energy gain is small with respect to the spherical minimum obtained previously. On the other hand, the HFB+PNAMP results also show some $\gamma$-softness but the absolute minimum is found at larger prolate deformations. Moreover, the energy gain is larger (740 keV) due to the particle-number projection that is not present in the HFB approach. Nevertheless, such an energy is still 300 keV above the energy obtained by a PNVAP+PNAMP method. The exact ground state provided by the ISM diagonalization $(E(0^{+}_{1})=-7.569$ MeV) is still below these minima (see caption of Fig. 1).

We now study the role of the $nn$-pairing correlations. To better visualize the influence of this degree of freedom, we display in Fig. 2 the results performing PNVAP calculations with constraints on the axial quadrupole deformation $(\beta_{2},\gamma=0^{\circ}\text{ or }180^{\circ})$ and the isovector pairing content $(\delta_{nn})$ defined in Eq. (14). In Figure 2(a), we show the intrinsic $nn$-pairing energy Ring and Schuck (1980), $E_{pair}=-\tfrac{1}{2}\mathrm{Tr}\left(\Delta\kappa^{*}\right)$, calculated using the PNVAP intrinsic wave functions without any projections.444i.e. using the pairing density $\kappa_{ab}=\langle\Phi|c_{b}c_{a}|\Phi\rangle$ and the pairing field $\Delta_{ab}=\frac{1}{2}\sum_{cd}\bar{v}_{abcd}\kappa_{cd}$. We observe that the pairing energy grows monotonically with increasing values of $\delta_{nn}$ and is also almost independent of $\beta_{2}$ (contour lines are horizontal). Therefore, the effect of varying the parameter $\delta_{nn}$ is just a way of changing the pairing content of the intrinsic wave functions Vaquero et al. (2011). Concerning the total energies (Figs. 2(b)-(c)), we observe a similar behavior to the $(\beta_{2},\gamma)$ case when we compare the PNVAP and the angular momentum projected $(J=0)$ surfaces. In the former, a distinct minimum around the spherical point $(\beta_{2}=0)$ and $\delta_{nn}\in[1.5,3.5]$ is found. This surface is softer in the pairing direction than in the quadrupole direction. The angular momentum projection, however, produces a much flatter TES (Fig. 2(c)) although slightly lower energies are obtained for larger values of $\delta_{nn}$. The energy difference between the minimum before and after projection is rather small ($\sim 26$ keV). This behavior, together with the degeneracy also found in the angular momentum projected surface in the $(\beta_{2},\gamma)$ plane, suggests that all these degrees of freedom should be taken into account simultaneously in a GCM calculation.

We compare now the different PGCM calculations (see Table 1) performed for the nucleus 48Ca with the exact results provided by the ISM framework. In Figure 3, the lowest eigenvalues for angular momenta $J=$ 0, 2, 3, 4 and 6 are shown. We first notice that all of the PGCM approaches are variational approximations to the exact solution plotted on the right, i.e., the energies are all above the energies obtained from the diagonalization of the Hamiltonian. We clearly see an improvement of the results whenever PNVAP rather than HFB intrinsic wave functions are used (from PGCM1 to PGCM2). Additionally, the explicit exploration of the pairing correlations (PGCM3) provides even closer results to the exact solutions. In fact, the difference between the exact and PGCM3 ground-state energies is $119$ keV showing the ability of these methods to approach the exact values. The reproduction of excited states with a larger value of angular momentum, however, is not as good as for the ground state. In particular, we found that the $3^{+}$ state is particularly sensitive to the variational approach. This effect can be understood as our minimization schemes probe primarily the ground-state energy. Such a deficiency can be corrected within the PGCM framework by enlarging the initial set of quasi-particle wave functions with: a) intrinsically rotating states through cranking calculations Borrajo et al. (2015); Egido et al. (2016); and/or, b) multi-quasi-particle excitations obtained through the blocking mechanism Chen and Egido (2016, 2017). Both improvements will be explored in the future but are beyond the scope of the present work.

Electric quadrupole transitions and moments are also computed and the results of the most relevant ones are written in Table 2. We observe a slight overestimation of the exact results in most of the transitions and moments. Furthermore, the PGCM3 method shows, as expected, the best agreement, and lower values of the spectroscopic quadrupole moments are obtained whenever the pairing correlations are better described (PNVAP vs. HFB minimization).

We now analyze two complementary aspects of the nuclear wave functions obtained by the PGCM framework, namely, the collective wave functions (c.w.f.) and the occupation numbers of the spherical single-particle orbits. The c.w.f. are useful quantities to reveal the role of the different collective degrees of freedom explored with constrained calculations and they are routinely computed in MREDF methods based on Skyrme, Gogny, and covariant functionals. These functions give the most relevant contributions of the collective coordinates in each individual nuclear state Ring and Schuck (1980); Bender et al. (2003); Robledo et al. (2018) and are defined as

[TABLE]

In Figure 4, we plot the c.w.f. for $0^{+}_{1,2}$, $2^{+}_{1,2}$ and $4^{+}_{1,2}$ states in the $(\beta_{2},\gamma)$ plane computed with the PGCM2 method. We first note that the ground state is flat in almost the whole range of deformations although the maximum contribution is found around the spherical point, as we could expect for a doubly-magic nucleus. This behavior is consistent with the degeneracy of the projected TES (Fig. 1(d)). The $2^{+}_{1}$ c.w.f. is also rather smooth and we find a considerable mixing of oblate and $\gamma=20^{\circ}$ deformed states at the edges of the available $\beta_{2}$ deformations. The rest of the c.w.f. represented in Fig. 4 are also localized at (or nearby) the border of the $\beta_{2}$ coordinate. Moreover, they show a rather sharp behavior which means that only very few intrinsic states actually contribute to the building of the nuclear states. We recall that the limits in $\beta_{2}$-values are determined by the maximum value that gives a convergent PNVAP solution. Finally, prolate deformed c.w.f. with a similar structure are obtained for the $0^{+}_{2}$, $2^{+}_{2}$ and $4^{+}_{3}$ (not shown) states. However, the $4^{+}_{1}$, $4^{+}_{2}$ excited states are mainly located at specific $\gamma=60^{\circ},40^{\circ}$ values, respectively, and, to a lesser extent, at other $\beta_{2}>0.1$ deformations.

It is important to note that prolate or oblate configurations of the $2^{+}$ and $4^{+}$ states obtained here are consistent with their values for the spectroscopic quadrupole moments given in Table 2. For example, the $2^{+}_{2}$ state is predominantly prolate deformed and the $Q(2^{+}_{2})$-value is negative, and so on.

Similar conclusions can be extracted from the c.w.f. given by the PGCM3 method. In Fig. 5 we show such quantities in the $(\beta_{2},\delta_{nn})$ plane with only axial deformed intrinsic states to simplify the discussion. The ground state is almost flat around $\beta_{2}\in\left[-0.1,0.1\right]$ with a mild dependence on $\delta_{nn}$. Again, the $0^{+}_{2}$ and $2^{+}_{2}$ c.w.f. are functions located at the prolate edge with two sharp peaks, one at the largest and the other at the smallest pairing values. This topology is also seen in the $2^{+}_{1}$ c.w.f. that is found at the oblate edge of the plot, consistently with the $2^{+}_{1}$ state obtained by the PGCM2 approach discussed above. In fact, the only significant difference between the c.w.f. calculated with the two approaches is found in the $4^{+}_{1}$ state. Here, both oblate states and a sharp region around prolate ($\beta_{2}=0.14$) and small-pairing values are present in this c.w.f. while in the $(\beta_{2},\gamma)$ plane the oblate configuration is the predominant one.

Interestingly, these c.w.f. (Figs. 4- 5) largely differ from those usually obtained by EDF calculations where, in general, the distributions found are smoother and better behaved at the boundaries Bender et al. (2003); Nikšić et al. (2011); Egido (2016); Robledo et al. (2018). One possible explanation is that in EDF calculations we have access to many more single-particle orbits that can evolve with the collective degrees of freedom following a Nilsson-like behavior. Such a picture allows for extending the exploration of these collective coordinates up to much larger values with the c.w.f. decaying progressively to zero.

More information can be learnt about the content of the correlated wave functions by computing their occupation numbers, i.e. the number of particles in a given $m_{t}$-orbit $\check{a}\equiv(\breve{a},m_{t})$,

[TABLE]

where

[TABLE]

and comparing the values to the decomposition of the exact eigenstates as obtained by the ISM calculations. A similar analysis was made in Ref. Rodríguez et al. (2016) in the context of nuclear EDF methods. The results are shown in Table 3 for the nuclear states discussed above. The differences among the various PGCM schemes are small but the best agreement with the exact results is consistently given by the PGCM3 approach. In fact, this approach is very close to the ISM results and reproduces very well the particle-hole structure of the ground and excited states. In particular, the ground state shows almost a pure $f_{7/2}$ closed-shell configuration, the $0^{+}_{2}$ and $2^{+}_{2}$ states are 2p-2h excitations in the $f_{7/2}$-$p_{3/2}$ orbits, the $2^{+}_{1}$ and $4^{+}_{1}$ states display 1p-1h configurations also in these orbits, and, finally, the $4^{+}_{2}$ state exhibits a 1p-1h configuration in the $f_{7/2}$-$f_{5/2}$ orbits. Therefore, the excitations found in this nucleus have, within our particular choice of basis, a well-defined single-particle nature. In addition, it is interesting to note that those kind of particle-hole configurations can be grasped by exploring a priori collective degrees of freedom such as quadrupole deformations and pairing correlations. Actually, this observation together with the shapes of the collective wave functions displayed above seem to indicate that, in this small valence space, these variables are just useful tools to generate highly correlated intrinsic states and the meaning of ”collectivity” is blurred.

III.2 Calcium isotopic chain

We now turn our attention towards the results obtained along the whole calcium isotopic chain in the $pf$-shell. In Fig. 6(a) we represent the ground-state energies for even-even and odd-even calcium nuclei computed exactly and with the different variational schemes. We see that the results are almost indistinguishable at the scale that reflects the absolute ground-state energies. To better visualize the accuracy of the different variational methods, we plot the difference between the approximate and exact energies, $\Delta E_{g.s.}=E^{approx}_{g.s.}-E^{ISM}_{g.s.}$, in Fig. 6(b). Here, we observe a better agreement with the exact result at the edges of the valence space than in the middle. As we will discuss below, the dimensions of the Hamiltonian matrices are the largest in the mid-shell and there will be configurations that cannot be captured only exploring the collective degrees of freedom studied here. The worst approximation is given by the HFB method, as expected, where a slight odd-even staggering is found. The smallest (largest) difference is $0.24$ MeV ($1.76$ MeV) obtained in the nucleus 42Ca (52Ca). The staggering is magnified in the results provided by the PGCM1 method that uses precisely the HFB intrinsic wave functions to define the subspace within which the GCM states are constructed. In this case, the even nuclei are better described than their odd neighbors except for the isotopes at the edges (42-45,56-58Ca) where the PGCM1 approximation gives ground state energies very close to the exact solution (for 42,43,57,58Ca they are on top of the ISM results). Using the PGCM1 method, the largest difference is reduced to $1.01$ MeV found in the isotope 51Ca. Smaller and smoother energy differences are obtained by PNVAP approaches (PNVAP itself and PGCM2,3). In particular, the superior treatment of pairing correlations in PNVAP approaches allows them to describe even- and odd-mass isotopes on an equal footing, i.e. there is no odd-even staggering in those cases. By constrast, HFB is a poor approximation whenever pairing correlations are weak, e.g. odd number parity states due to the blocking effect. In the best case, namely, the PGCM3 scheme, we obtain a largest difference of $0.24$ MeV (at 51Ca) proving the ability of these variational methods to approach closely the exact ground states.

We can also analyze the capability of these methods to reproduce the low-lying excitation spectra. In that order, the excitation energies of the lowest states computed with ISM and PGCM methods are displayed as a function of the angular momentum in Fig. 7 and Fig. 8 for even and odd calcium isotopes, respectively. The overall agreement with the exact results is again very good and the quantitative behavior of the excitation energies is well reproduced. In fact, the ISM solutions are obtained exactly at the edges of the valence space even with PGCM1,2 approaches. The quantitative agreement deteriorates in the mid-shell but it is improved by the inclusion of the pairing fluctuations within the PGCM3 method. We also notice in that case a better description of the exact spectra in odd-mass nuclei. The effect of pairing on the final spectra can be studied by comparing the PGCM approximations based on the HFB solutions with those based on the PNVAP ones. For even-mass nuclei, the PGCM1 approach produces systematically larger excitation energies in even nuclei that are subsequently reduced in the PGCM2,3 methods. For odd-mass nuclei, the PGCM2 approach gives systematically larger excitation energies that are corrected once the pairing fluctuations are included. Consequently, the PGCM3 calculations provide almost the exact spectra in the odd calcium isotopes.

For even-mass isotopes, we display in Fig. 9 the systematics of the lowest excitations with an even value of $J$. Again, we observe a very good reproduction of the exact results by the PGCM methods. In particular, these approximations are able to reproduce the main trends, e.g., the peaks at $N=28$ and $32$ in the $2^{+}_{1,2}$ states associated to shell closures, the sudden decrease of the $0^{+}_{2}$ energy at $N=34$, or the bell-shaped $4^{+}_{1}$ excitation energies. However, the PGCM1,2 approximations have problems in describing the $0^{+}_{2}$ and $4^{+}_{1}$ excitation energies in the vicinity of 48Ca. This discrepancy is partially corrected by the addition of the pairing fluctuations (PGCM3) although there is some room left to improvements, which could be corrected through the inclusion of additional degrees of freedom such as the rotational frequency or the explicit multi-quasi-particle excitations. Nevertheless, the PGCM3 approximation is almost on top of the exact results in most of the excitation energies shown in Fig. 9.

The study of the occupation numbers performed above for 48Ca can be extended to all even calcium isotopes. For the sake of simplicity, we only represent in Fig. 10 the difference between the ISM and the PGCM2 methods for the two first $0^{+}$ and $2^{+}$ states. We observe that the difference for the ground and the $2^{+}_{1}$ states are almost negligible in the whole isotopic chain. Slightly larger differences are found for $0^{+}_{2}$ and $2^{+}_{2}$ in the mid-shell nuclei (46-52Ca), consistently with the discrepancies shown above for ground and excited state energies. However, these differences are smaller than 0.2 particles in most of the cases. Therefore, these results prove the ability of the PGCM methods to describe the underlying single-particle structure of the nuclear wave functions by exploring coordinates that depict, a priori, collective degrees of freedom.

These results suggest that, for all the values of the angular momentum $J$ examined here, the lowest eigenstates of the full many-body Hilbert $J$-subspaces can be well approximated by the variational mixing of few selected correlated many-body states. As described in Sec. II, our approximated wave functions are built as optimal superpositions of the symmetry-projected quasi-particle states included in our trial set. This is equivalent to the diagonalization of the Hamiltonian in the subspace spanned by the set of projected states $\left\{P^{J}_{MK}P^{N}P^{Z}P^{\pi}|\Phi(q)\rangle,\forall\,qK\right\}$. In practice, however, due to the linear redundancy among the projected states, the dimension of the subspace within which the Hamiltonian is effectively diagonalized is much smaller. In Table 4, we compare the dimensions of the exact diagonalization using Slater determinants coupled to $J$ with the dimensions of the GCM diagonalization using the natural basis states for even and odd calcium isotopes (the dimensions are symmetric with respect to particles and holes). For simplicity, we compare again only the PGCM2 and ISM methods. The number of intrinsic states, $|\Phi(\beta_{2},\gamma)\rangle$, is 50 in this case. Since $K$-mixing is also performed, the actual number of states mixed within the GCM framework for each angular momentum is $50\times(2J+1)$. We observe in Table 4 a large reduction of this number whenever the corresponding natural bases are extracted, i.e., a substantial redundancy given by the linear dependence of the projected states is obtained. It is important to point out that at the edges of the valence space (42-58Ca, 43-57Ca) the dimensions of the natural bases are the same or almost the same as the bases built with Slater determinants coupled to $J$.

Therefore, in those cases the ISM and PGCM methods are necessarily equivalent and that explains the perfect agreement between exact diagonalizations and the variational approaches. However, for systems in the mid-shell, the number of many-body states in the natural bases are two orders of magnitude smaller than the shell-model bases and still the description is very good. Some pros and cons of this procedure can be mentioned. The main advantage of the PGCM methods are their scalability with the number of configurations/shells in the valence space. In this sense, these methods can be used in a broader context. We also see that these techniques are straightforwardly extendable by including additional degrees of freedom (pairing fluctuations, cranking terms, explicit quasi-particle excitations), even though the ground states are rather well approximated by including only the most relevant variables (quadrupole and pairing). However, the addition of more constraints is also costly from a computational point of view and the method can be eventually as complicated as the full diagonalization. We also observe a poorer description of the exact results in those systems where the number of Slater determinants is very much larger than the dimension of the corresponding natural basis (e.g., states in 46-54Ca). Concerning the latter point, it is obvious that only few (out of the total) number of nuclear states for a given angular momentum can be obtained by the current implementations of the PGCM method. Nevertheless, the present study shows that these methods are able to reproduce very well the exact results by mixing a very reduced set of strongly correlated basis states.

IV Summary and outlook

In this article, we have studied the performance of several PGCM approximation schemes in reproducing the exact eigenstates of the shell model Hamiltonian KB3G in the $pf$-shell valence space. The PGCM is a versatile variational method that includes beyond-mean-field correlations through the configuration mixing of a set symmetry-projected intrinsic states. In this work, the intrinsic states considered are Bogolyubov quasi-particle states that explore the quadrupole and the pairing degrees of freedom, and they have been generated by solving a series of constrained HFB or PNVAP equations.

We used the even and odd calcium isotopes in the $pf$-shell valence space to benchmark the different PGCM implementations. In general, we obtained a very good qualitative and quantitative description of ground- and excited-state exact energies. Furthermore, the use of PNVAP intrinsic states and the simultaneous exploration of quadrupole and $nn$-pairing degrees of freedom improved significantly the agreement with the shell model results, especially in the odd systems.

We have also analyzed the PGCM correlated wave functions in terms of the collective coordinates and their occupation numbers. For the nucleus 48Ca, we have found a very flat ground-state collective wave function and sharply peaked excited-state collective wave functions. These distributions seem to indicate that the meaning of the collective degrees of freedom in these isotopes and valence spaces is not as clear as in well-deformed nuclei and no-core configuration spaces. Recently, a similar conclusion has been drawn after analyzing the $(\beta_{2},\gamma)$ values and their variances extracted from the Kumar invariants Poves et al. (2019) in spherical isotopes. The non-collective nature of the excitations in this isotopic chain could nonetheless be reproduced by the present PGCM implementations as the analysis of the occupation numbers pointed out.

Altogether, our PGCM calculations have demonstrated that a very good approximation to the lowest eigenstates of the shell-model Hamiltonian can be obtained through the diagonalization of the aforementioned Hamiltonian within the subspace spanned by a small set of selected correlated wave functions instead of very large $J$-coupled Slater determinant basis. This is in line with the conclusions drawn in Refs. Ripoche et al. (2017, 2018) in the context of a truncated configuration-interaction method. There were, however, some discrepancies for some states belonging to many-body Hilbert spaces with large dimensions (i.e. for mid-shell isotopes). In those cases, the exploration of additional degrees of freedom, such as the rotational frequency or the explicit multi-quasi-particle excitations, in the generation of the intrinsic states would be needed. This is a work in progress.

It would certainly be interesting to extend the present study by testing the quality of our variational approximations in a larger valence space that includes more than a single harmonic oscillator shell. Finally, the role of other pairing channels, in particular, the isoscalar and isovector proton-neutron pairing degrees of freedom will be studied in the near future.

Acknowledgements

We thank L. M. Robledo for proofreading the manuscript. We also thank him, A. Poves, and F. Nowacki for useful discussions. The work of TRR was supported by the Spanish MINECO under Programa Ramón y Cajal 11420, FIS-2014-53434-P and by the Spanish MICINN under PGC2018-094583-B-I00, and, in part, by the ExtreMe Matter Institute EMMI at the GSI-Darmstadt, Germany. TRR gratefully thanks the support from GSI-Darmstadt and TU-Darmstadt computing facilities.