Symmetry restoration in mean-field approaches

TL;DR

This paper reviews symmetry restoration techniques in mean-field approaches for quantum many-body systems, discussing group-theory methods, approximations, and applications in nuclear and mesoscopic systems, highlighting unresolved issues.

Contribution

It provides a comprehensive overview of symmetry restoration methods, including detailed expressions and practical applications, addressing computational challenges and unresolved problems.

Findings

Symmetry restoration improves comparison with experimental data.

Approximate projection methods are effective in small systems.

Applications span nuclear models and energy density functional theory.

Abstract

The mean-field approximation based on effective interactions or density functionals plays a pivotal role in the description of finite quantum many-body systems that are too large to be treated by ab initio methods. Some examples are strongly interacting medium and heavy mass atomic nuclei and mesoscopic condensed matter systems. In this approach, the linear Schrodinger equation for the exact many-body wave function is mapped onto a non-linear density-dependent one-body potential problem. This approximation, not only provides computationally very simple solutions even for systems with many particles, but due to the non-linearity, it also allows for obtaining solutions that break essential symmetries of the system, often connected with phase transitions. In this way, additional correlations are subsumed in the system. However, the mean-field approach suffers from the drawback that the corresponding wave functions do not have sharp quantum numbers and, therefore, many results cannot be compared directly with experimental data. In this article, we discuss general group-theory techniques to restore the broken symmetries, and provide detailed expressions on the restoration of translational, rotational, spin, isospin, parity and gauge symmetries, where the latter corresponds to the restoration of the particle number. In order to avoid the numerical complexity of exact projection techniques, various approximation methods available in the literature are examined. Applications of the projection methods are presented for simple nuclear models, realistic calculations in relatively small configuration spaces, nuclear energy density functional theory, as well as in other mesoscopic systems. Further, unresolved problems in the application of the symmetry restoration methods to the energy density functional theories are highlighted in the present work.

Tables

Table 1. Table 1: Various symmetries and the corresponding parts of the projection operators.

symmetry	shift	generator	group operator
	g	$\hat{S}$	$R\left(g\right)$
rotation (1D)	$\alpha$	ang. mom. ${\hat{J}}_z$	$e^{-i\alpha {\hat{J}}_z}$
rotation (3D)	$\alpha ,\beta ,\gamma$	ang. mom. $\widehat{𝑱}$	$e^{-i\alpha {\hat{J}}_z}{e}^{-i\beta {\hat{J}}_y}{e}^{-i\gamma {\hat{J}}_z}$
translation (3D)	$𝒂$	momentum $\widehat{𝑷}$	$e^{-i𝒂\cdot \hat{P}}$
rotation in gauge space	$\varphi$	part. number $\hat{N}$	$e^{-i\varphi \hat{N}}$
isorotation (1D)	$\phi$	isospin ${\hat{T}}_3$	$e^{-i\phi {\hat{T}}_3}$
isorotation (3D)	$\phi ,\theta ,\psi$	isospin $\widehat{𝑻}$	$e^{-i\phi {\hat{T}}_3}{e}^{-i\theta {\hat{T}}_2}{e}^{-i\psi {\hat{T}}_3}$
parity	$\varphi =0,\pi$	${\hat{N}}_{-}$	$e^{-i\varphi {\hat{N}}_{-}}$

Table 2. Table 2: Table  1 continued.

symmetry	weight function	eigenvalue	proj. operator
	$D^{I\ast }\left(g\right)$	I	in equation
rotation (1D)	$\frac{1}{2\pi }{e}^{i\alpha M}$	$M$	(62)
rotation (3D)	$\frac{2I+1}{8{\pi }^2}{D}_{MK}^{I\ast }(\alpha ,\beta ,\gamma )$	$I,M,K$	(64)
translation (3D)	$\frac{1}{{\left(2\pi \right)}^3}{e}^{i𝒂\cdot P}$	$𝑷$
rotation in gauge space	$\frac{1}{2\pi }{e}^{i\varphi N}$	$N$	(56)
isorotation (1D)	$\frac{1}{2\pi }{e}^{i\phi T_3}$	$T_3=\frac{1}{2}\left(N-Z\right)$
isorotation (3D)	$\frac{2T+1}{8{\pi }^2}{D}_{T_3,T_3^{\prime }}^{T\ast }(\phi ,\theta ,\psi )$	$T,T_3,T_3^{\prime }$
parity	$\frac{1}{2}$,$\frac{p}{2}$	$p$	(57)

Equations

$$|{\Psi }^{\pm }⟩$$

$$|{\Psi }^{\text{mbox}left}⟩$$

$$|{\Psi }^{\text{mbox}right}⟩$$

$$\begin{array}[]{rcl}\hat{H}|\Psi^{+}\rangle&=&E^{+}|\Psi^{+}\rangle\\ \hat{H}|\Psi^{-}\rangle&=&E^{-}|\Psi^{+}\rangle\\ \end{array}$$

$$\begin{array}[]{rcl}\langle\Psi^{\mbox{\scriptsize{left }}}|\hat{H}|\Psi^{\mbox{\scriptsize{left }}}\rangle&=&E^{\mbox{\scriptsize{loc}}}\\ \langle\Psi^{\mbox{\scriptsize{right}}}|\hat{H}|\Psi^{\mbox{\scriptsize{right}}}\rangle&=&E^{\mbox{\scriptsize{loc}}}\\ \langle\Psi^{\mbox{\scriptsize{left }}}|\hat{H}|\Psi^{\mbox{\scriptsize{right}}}\rangle&=&-{\textstyle\frac{1}{2}}\delta{E}\\ \langle\Psi^{\mbox{\scriptsize{right}}}|\hat{H}|\Psi^{\mbox{\scriptsize{left }}}\rangle&=&-{\textstyle\frac{1}{2}}\delta{E}\\ \end{array}$$

$$\delta E=E^{-}-E^{+},E^{\text{mbox}loc}=\frac{1}{2}\left(E^{+}+E^{-}\right),$$

$$\Big{(}\langle\Psi^{+}|,\langle\Psi^{-}|\Big{)}\hat{H}\left(\begin{array}[]{l}|\Psi^{+}\rangle\\ |\Psi^{-}\rangle\end{array}\right)=\left(\begin{array}[]{ll}E^{+}&0\\ 0&E^{-}\\ \end{array}\right)$$

$$\Big{(}\langle\Psi^{\mbox{\scriptsize{left }}}|,\langle\Psi^{\mbox{\scriptsize{right}}}|\Big{)}\hat{H}\left(\begin{array}[]{l}|\Psi^{\mbox{\scriptsize{left }}}\rangle\\ |\Psi^{\mbox{\scriptsize{right}}}\rangle\end{array}\right)=\left(\begin{array}[]{rr}E^{\mbox{\scriptsize{loc}}}&-{\textstyle\frac{1}{2}}\delta{E}\\ -{\textstyle\frac{1}{2}}\delta{E}&E^{\mbox{\scriptsize{loc}}}\\ \end{array}\right),$$

$$|{\Phi }^{\pm }⟩=\frac{1}{\sqrt{2\pm 2ϵ}}\left(|{\Phi }^{\text{mbox}left}⟩\pm |{\Phi }^{\text{mbox}right}⟩\right)\text{mbox}forϵ=⟨{\Phi }^{\text{mbox}left}|{\Phi }^{\text{mbox}right}⟩.$$

$$\hat{\Pi }|{\Phi }^{\text{mbox}left}⟩=|{\Phi }^{\text{mbox}right}⟩\text{mbox}and\hat{\Pi }|{\Phi }^{\text{mbox}right}⟩=|{\Phi }^{\text{mbox}left}⟩,$$

$$ϵ=⟨{\Phi }^{\text{mbox}left}|{\Phi }^{\text{mbox}right}⟩=⟨{\Phi }^{\text{mbox}left}|\hat{\Pi }|{\Phi }^{\text{mbox}left}⟩=⟨{\Phi }^{\text{mbox}right}|\hat{\Pi }|{\Phi }^{\text{mbox}right}⟩=⟨{\Phi }^{\text{mbox}right}|{\Phi }^{\text{mbox}left}⟩,$$

$$\Big{(}\langle\Phi^{+}|,\langle\Phi^{-}|\Big{)}\hat{H}\left(\begin{array}[]{l}|\Phi^{+}\rangle\\ |\Phi^{-}\rangle\end{array}\right)=\left(\begin{array}[]{ll}{\cal E}^{+}&0\\ 0&{\cal E}^{-}\\ \end{array}\right)$$

$$\Big{(}\langle\Phi^{\mbox{\scriptsize{left }}}|,\langle\Phi^{\mbox{\scriptsize{right}}}|\Big{)}\hat{H}\left(\begin{array}[]{l}|\Phi^{\mbox{\scriptsize{left }}}\rangle\\ |\Phi^{\mbox{\scriptsize{right}}}\rangle\end{array}\right)=\left(\begin{array}[]{rr}{\cal E}^{\mbox{\scriptsize{loc}}}&\Delta\\ \Delta&{\cal E}^{\mbox{\scriptsize{loc}}}\\ \end{array}\right),$$

$${\mathcal{E}}^{\pm }$$

$$\delta \mathcal{E}$$

$${\mathcal{E}}^{\text{mbox}loc}$$

$$\Delta$$

$$\Phi^{S_{z}=0}_{\rm UHF}(1,2)=\frac{1}{\sqrt{2}}\left|\begin{array}[]{cc}\varphi_{l}({\bm{r}}_{1})\alpha_{1}&\;\varphi_{l}({\bm{r}}_{2})\alpha_{2}\\ \varphi_{r}({\bm{r}}_{1})\beta_{1}&\;\varphi_{r}({\bm{r}}_{2})\beta_{2}\end{array}\right|,$$

$${\hat{P}}_{\mathrm{spin}}^s=(1-{\varpi }_{12})/2,$$

$$2\sqrt{2}{\hat{P}}_{\mathrm{spin}}^s{\Phi }_{\mathrm{UHF}}^{S_z=0}(1,2)=|{\phi }_l(1){\overline{\phi }}_r(2)⟩-|{\overline{\phi }}_l(1){\phi }_r(2)⟩.$$

$$2{\hat{P}}_{\mathrm{spin}}^s{\Phi }_{\mathrm{UHF}}^{S_z=0}(1,2)$$

$$\chi (s=0,S_z=0)=({\alpha }_1{\beta }_2-{\alpha }_2{\beta }_1)/\sqrt{2}.$$

$${\Psi }_{\text{mbox}PRJ}^{\mathrm{s}}(1,2)=N_{+}\sqrt{2}{\hat{P}}_{\mathrm{spin}}^s{\Phi }_{\mathrm{UHF}}^{S_z=0}(1,2),$$

$$N_{+}=1/\sqrt{1+S_{lr}^2},$$

$$S_{lr}=\int {\mathrm{d}}^2\mathbf{r}{\phi }_l(\mathbf{r}){\phi }_r(\mathbf{r}).$$

$$E_{\text{mbox}PRJ}^{\mathrm{s}}=N_{+}^2\left[h_{ll}+h_{rr}+2S_{lr}{h}_{lr}+J_{lr}+K_{lr}\right],$$

$$E_{\text{mbox}HF}^{\mathrm{s}}=h_{ll}+h_{rr}+J_{lr}.$$

$$E_{\mathrm{HF}}^{\mathrm{t}}=h_{0n,0n}+h_{1n,1n}+J_{0n,1n}-K_{0n,1n}.$$

$$\hat{H}=-G{\hat{S}}_{+}{\hat{S}}_{-},$$

$${\hat{S}}_{+}=({\hat{S}}_{-}{)}^{†}=\sqrt{\frac{\Omega }{2}}[a_j^{†}{a}_j^{†}{]}_{J=0},$$

$${\hat{S}}_0=\frac{1}{2}(\hat{N}-\Omega ),$$

$${\hat{\mathbf{S}}}^2={\hat{S}}_{+}{\hat{S}}_{-}+{\hat{S}}_0^2-{\hat{S}}_0,$$

Full text

††: J. Phys. G: Nucl. Phys.

1. 1 Introduction

2. 2 Symmetry breaking in simple illustrative models

3. 2.1 Doubly symmetric potential well

4. 2.2 Dissociation of the natural molecular hydrogen and other similar two-dimensional artificial dimers

5. 2.3 The seniority model

6. 3 Symmetry Restoration - General Formalism

7. 3.1 Projection Operator - Mathematical Basis

8. 3.1.1 Group-Theory Approach

9. 3.1.2 Generator Coordinate Method

10. 3.1.3 Generic properties of projection

11. 3.2 Projection Methods for Various Symmetries

12. 3.2.1 Particle-number

13. 3.2.2 Linear Momentum

14. 3.2.3 Parity

15. 3.2.4 One-dimensional rotation

16. 3.2.5 Three-dimensional rotation

17. 3.3 Symmetry restoration of the HFB wave function

18. 3.3.1 Rotated Norm and Energy Kernels

19. 3.3.2 Variation after projection

20. 4 Approximate projection methods

21. 4.1 The Lipkin method

22. 4.2 The Lipkin-Nogami method

23. 4.3 The Kamlah method

24. 4.4 Applicability and Applications

25. 4.4.1 First-order terms

26. 4.4.2 Second-order terms

27. 4.4.3 Higher-order terms

28. 5 Projection methods in simple nuclear models and ab initio calculations

29. 5.1 The pairing-plus-quadrupole model

30. 5.2 Projection in small configuration spaces

31. 5.3 Projection methods in ab initio calculations

32. 6 Projection methods and nuclear density functional theory

33. 6.1 Difficulties encountered in restoring symmetries with nuclear EDF

34. 6.1.1 Basis not closed under the symmetry action

35. 6.1.2 Self-energy and Pauli principle

36. 6.1.3 Density-dependent prescription

37. 6.1.4 Non-integer powers of the density

38. 6.1.5 Future implementations

39. 6.2 Applications of symmetry restoration with nuclear EDFs

40. 6.2.1 Non relativistic EDFs

41. 6.2.2 Relativistic EDFs

42. 6.3 Approximate projection for nuclear EDF

43. 7 Projection methods in other mesoscopic systems beyond atomic nuclei

44. 7.1 The general hierarchical methodology.

45. 7.2 Quantum dots

46. 7.2.1 The microscopic many-body Hamiltonian.

47. 7.2.2 Mean-field equations for electrons: UHF wave functions using the Pople-Nesbet equations or Slater determinants with displaced Gaussians.

48. 7.2.3 Combining spin and total-angular-momentum restorations.

49. 7.2.4 More on spin restoration.

50. 7.2.5 Molecular symmetries of the UHF wave functions and magic angular momenta.

51. 7.2.6 A tour of the literature and Wigner molecules.

52. 7.3 Trapped ultracold neutral atomic gases and ions

53. 7.4 Spin-projected UHF, Hubbard systems, and natural molecules

54. 7.5 Other electronic systems

55. 7.6 Other emerging directions

56. 7.6.1 Relation to entanglement and quantum information science.

57. 7.6.2 Time evolution in finite systems out of equilibrium.

58. 8 Projected statistics

59. 8.1 Symmetry restoration at finite temperature

60. 8.2 Thermo-field dynamics

61. 9 Summary, conclusions, and perspectives

62. 10 Acknowledgments

63. A Overlaps and matrix elements between HFB states

1 Introduction

Mean-field approaches play a central role in the description of quantum many-body problems in areas like quantum chemistry, atomic, molecular, condensed matter, and nuclear physics. The simplicity of the associated wave functions, both in the fermion and the boson cases, is the reason behind the popularity of the mean-field approaches. These product-type wave functions allow, on the one hand, an easy implementation of the symmetrization principle of quantum mechanics required for identical particles, and, on the other hand, permit the application of techniques used in the field theory, like the Wick’s theorem, which enormously simplify the evaluation of the matrix elements of the many-body operators.

The optimal mean field, generating the single-particle orbitals, is usually determined through the application of the variational principle. In the fermion case, the variational principle, performed in the space of Slater determinants, leads to the familiar Hartree-Fock (HF) method. It is common to find situations where short-range attractive interactions induce correlations leading to the superfluidity or superconductivity phenomena that are well described by the BCS theory. The quasiparticles introduced in the BCS theory can be combined with the concepts of the HF theory to give the Hartree-Fock-Bogoliubov (HFB) mean-field theory, which is widely employed in nuclear physics.

Another facet of the mean-field approach is revealed within the Density Functional Theory (DFT). By using one-body densities as efficient and relevant degrees of freedom, DFT aims to map the exact wave functions of many-body systems onto the product states, which leads to the dynamical equations becoming formally identical to those given by the HF method. Although the foundations and approximations leading to HF and DFT methods are different, the similar structure of dynamical equations allows us to use for both, HF and DFT, the common name of mean-field approach.

One of the most salient characteristics of mean-field approaches is the fact that solutions often spontaneously break symmetries of the Hamiltonian. This is the case, for instance, in the BCS and HFB theories, where the associated mean-field wave functions, which are the vacua of the corresponding quasiparticle operators, do not represent states with good particle number. The spontaneous symmetry-breaking mechanism provides a way to incorporate nontrivial dynamic correlations on top of the simple Slater determinants, while preserving the simplicity of the mean-field description.

It is the nonlinearity of the mean-field equations that favors the spontaneous symmetry-breaking mechanism, and may constitute the simplest description of the symmetry-breaking effects occurring at a more fundamental level. In this way, the mean-field symmetry breaking leads to interesting perspectives to understand the physics of a given problem. For example, it leads to an easy and efficient description of different collective effects, such as the appearance of rotational bands being the result of the rotational-symmetry breaking, so common in nuclear or molecular physics. Symmetry breaking can also be a useful concept in the presence of stationary degenerate symmetry-conserving states, whereupon the famous Jahn-Teller effect becomes effective.

Nevertheless, advantages of the spontaneous mean-field symmetry breaking come at a price: the resulting wave functions are not invariant or covariant with respect to broken-symmetry groups and, therefore, they cannot be labeled with the symmetry quantum numbers such as angular momentum, parity, etc. This represents a serious drawback if quantities like electromagnetic transition probabilities, with their selection rules, are to be computed. Another drawback of the symmetry-breaking mechanism is connected with sharp transitions observed between the symmetry-conserving and symmetry-breaking solutions, which may occur as a function of some parameters of the Hamiltonian. Such sharp transitions are typical for infinite systems, but cannot characterize finite many-body or mesoscopic systems.

A way to overcome the disadvantages of the symmetry-breaking mean field is to restore the broken symmetries, which is the subject matter of this review article. A general idea of such a restoration is to take linear combinations, properly weighted, of wave functions obtained by applying the elements of the symmetry group to the symmetry-breaking (often called “deformed") mean-field wave function. As a consequence of the symmetry of the Hamiltonian, these “rotated” deformed wave functions are degenerate, and linear combinations, thereof, are expected to have lower energies. In addition, if the weights of the linear combination are chosen according to the rules of group theory, the obtained wave functions become invariant or covariant with respect to the underlying symmetry group and can be labeled with proper quantum numbers. This procedure is denoted in the literature as “symmetry restoration" or “projection". As it leads to linear combinations of product states, it can be understood as introducing correlations beyond the mean-field approach.

The theory behind the projection method is rooted in group theory as the weights of the linear combination are given by the irreducible representations of the symmetry group. In most of the cases, the symmetry group is a continuous Lie group (rotation, translation, particle-number gauge, etc.) while in other cases it is a discrete group (parity).

Once the structure of the projector operator is fixed, two alternatives are available to determine the deformed intrinsic state. The simplest one is to restore the symmetry of the deformed state obtained after solving the HF, HFB, or Kohn-Sham equations. This procedure is called “projection after variation”, where variation refers to the minimization of the mean-field energy. Another approach, fully self-consistent and variational, called “variation after projection", determines the deformed state through minimization of the projected energy, separately for each quantum numbers. In this way, different deformed states are obtained for a given quantum number. It turns out that extra flexibility brought in by the variation after projection method is able to smear out the sharp transitions mentioned before.

To compute basic quantities involved in the symmetry restoration, one takes advantage of the generalized Wick’s theorem, which allows for calculating matrix elements of operators between mean-field states. This theorem can be applied to projection, because the rotated mean-field states are mean-field states again – they simply correspond to rotated (quasi)particles. This property is a direct consequence of the Thouless theorem, and of the fact that the Lie algebras of the relevant symmetry groups can be represented in terms of one-body operators.

The program to perform the projected calculations can be directly implemented when the problem is defined in terms of a Hamiltonian operator. However, this is not always the case, and in many applications the Hamiltonian is replaced by a density functional in the spirit of the Hohenberg-Kohn or Kohn-Sham approach. In this situation, one is forced to introduce some sort of prescriptions on how to compute the energy kernels. However, these recipes are plagued with conceptual problems that were not resolved yet in the most general case.

Further, it is also common to use different (effective) interactions for each of the three contributions to the energy coming from a two-body operator, namely the direct, exchange, and pairing contributions (the most typical case is probably the use of the Slater approximation for the Coulomb exchange contribution and the neglect of the Coulomb antipairing field). In this case, a naive use of the generalized Wick’s theorem can lead to spurious contributions and specific ways to deal with this problem have to be devised. These difficulties represent serious impediments for the practical implementation of the symmetry-restoration methods in nuclear physics.

Since the effective interactions that are most often used to construct density functionals usually contain density-dependent terms, the corresponding mean-field approaches are referred to as based on the energy density functionals (EDFs), see Duguet (2014a); Schunck (2019); Furnstahl (2020) for further reading. The main benefit of such approaches is their applicability over the entire Segré chart. Very successful effective interactions were developed over the years with great success in describing bulk nuclear properties. The success of the EDF methods motivated the introduction of the symmetry restoration (beyond mean field step) aimed to gain access to symmetry-conserving observables and to increase the accuracy of the bulk properties.

The authors of this review are fully aware of the fact that nuclear physics specialists may find it difficult to follow the sections on electronic wave-function-based methods and vice versa. The two domains were evolving separately indeed. However, as far as the symmetry restoration is concerned, apart from apparent differences in the notation and language, the similarities are abundant. Obviously we were confronted with the question of unification of notation. We came to the conclusion that a single review would never change what has been deeply encrusted in the historical development of both domains. Therefore, we think that promoting a unified notation and language would not be a good idea, because then the text would become alien to everybody. We thus decided to always keep the original terminology and as clearly as possible explain what various names and symbols mean.

The present review article is organized in the following manner. Basic ideas regarding the symmetry breaking and restoration are presented in section 2. Section 3 outlines the general formalism of the symmetry restoration and section 4 discusses the applicability of the approximate projection methods. Symmetry restoration methods in simple nuclear models are discussed in section 5, nuclear DFT in section 6, and other mesoscopic systems in section 7. A brief description of the projection in statistical approaches is given in section 8 and finally, we provide a summary and concluding remarks in section 9.

2 Symmetry breaking in simple illustrative models

In this section, we introduce the subject matter of symmetry breaking and restoration by presenting three very simple examples. Firstly, in section 2.1, we discuss a solution of a one-particle problem in one dimension, where neither mean-field nor many-body complications appear. Secondly, in section 2.2, we present a two-particle problem in two dimensions, where one can illustrate the role of the mean-field approximation. And third, in section 2.3, we discuss the case of a many-body setting. These three simple models are exactly solvable, which allow us to analyze the problem of symmetry breaking and restoration in the quantum mechanical context and to clearly delineate the role of approximations that unavoidably have to be made in realistic situations.

2.1 Doubly symmetric potential well

Consider the doubly symmetric potential well Sakurai (1994), that is, a one-dimensional infinite potential well of width $2a$ with a step-like potential barrier of width $2b$ and height $V$ placed in the middle. To link this example to nuclear-physics scales of mass, distance, and energy, let us use the parameters of $\hbar^{2}/2m=20$ MeV fm2, $a=10$ fm, $b=1$ fm, and $V=40$ MeV. In this model, exact wave functions can be very easily determined; those of the two lowest eigenstates are plotted in figure 1(a).

The model is symmetric with respect to the middle of the well, and thus the eigenstates are either symmetric or antisymmetric, $\hat{\Pi}|\Psi^{\pm}\rangle=\pm|\Psi^{\pm}\rangle$, where $\hat{\Pi}$ is the inversion operator $x\rightarrow{-x}$. The parameters of the model are chosen in such a way that the two lowest states reside predominantly within the left and right well, and not in the barrier region.

The two lowest eigenstates of opposite parity can be expressed as linear combinations of two localized configurations, that is,

[TABLE]

for

[TABLE]

see figure 1(b). That is, the localized configurations are wave packets built of the two lowest eigenstates of the system. In these configurations, the particle resides either in the left or in the right well. The four states are pairwise orthogonal, $\langle\Psi^{+}|\Psi^{-}\rangle=0$ and $\langle\Psi^{\mbox{\scriptsize{left }}}|\Psi^{\mbox{\scriptsize{right}}}\rangle=0$, and both, the pair of exact states, $|\Psi^{+}\rangle$ and $|\Psi^{-}\rangle$, and that of localized wave packets, $|\Psi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Psi^{\mbox{\scriptsize{right}}}\rangle$, span the same subspace of the two lowest eigenstates.

If we denote the exact Hamiltonian of the doubly symmetric potential well by $\hat{H}$, we obviously have

[TABLE]

and

[TABLE]

for

[TABLE]

where in our model the exact eigenstates are split in energy by $\delta{E}=75.6$ keV and the average energies of both wave packets $E^{\mbox{\scriptsize{loc}}}$ are, of course, the same and located exactly in the middle between the two eigenenergies. In the matrix notation, equations (4) and (5) can be represented as

[TABLE]

and

[TABLE]

respectively.

It is now very important to realize that by breaking the symmetry of the problem, we can build a very reasonable model of the localized wave packets, see figure 1(c). Indeed, by keeping only the left or right potential well, we obtain the left and right broken-symmetry states, $|\Phi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Phi^{\mbox{\scriptsize{right}}}\rangle$. The broken-symmetry states are the exact eigenstates in the modified potential wells, but at the same time they approximate exact localized states in the original doubly symmetric potential well. In the scale of figure 1, they cannot really be distinguished from the exact wave packets $|\Psi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Psi^{\mbox{\scriptsize{right}}}\rangle$. Note that $|\Phi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Phi^{\mbox{\scriptsize{right}}}\rangle$ are stationary in the modified potential wells, whereas we use them to model non-stationary wave packets $|\Psi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Psi^{\mbox{\scriptsize{right}}}\rangle$ of the original doubly symmetric potential well.

At this point, we arrive at the very heart of the subject matter of this article. Namely, the symmetry-broken solutions, which pertain to a different problem than the original one, can serve us as approximate solutions of the original problem. This is achieved by restoring their symmetry, that is, by considering the normalized symmetric and antisymmetric combinations of the normalized states $|\Phi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Phi^{\mbox{\scriptsize{right}}}\rangle$, see figure 1(d),

[TABLE]

Since the inversion transforms the two broken-symmetry states one into another,

[TABLE]

states (9) have correct symmetry properties of $\hat{\Pi}|\Phi^{\pm}\rangle=\pm|\Phi^{\pm}\rangle$ with their normalization factors that depend on the overlap $\epsilon$ between the approximate broken-symmetry states. Note that since the inversion is a hermitian operator, the overlap $\epsilon$,

[TABLE]

is real.

For one-dimensional representations of the symmetry group, as is the case of the parity symmetry discussed in this section (see also discussion in section 3.2), the symmetry restoration of the broken-symmetry states, equation (9), automatically brings the Hamiltonian matrix to its diagonal form. Then, using the symmetry condition $[\hat{H},\hat{\Pi}]=0$, the analogues of equations (7) and (8) read

[TABLE]

and

[TABLE]

respectively, where the approximate eigenenergies ${\cal E}^{\pm}$ and their splitting $\delta{\cal E}$,

[TABLE]

depend on the matrix elements of the Hamiltonian,

[TABLE]

where, by the argument analogous to that used when deriving equation (11), $\Delta$ is real.

Figure 2 summarizes the logic of the construction presented above and depicts energies of all discussed states. In the left panel, we show how the pair of exact eigenstates, $|\Psi^{+}\rangle$ and $|\Psi^{-}\rangle$, is transformed into the pair of exact localized wave packets, $|\Psi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Psi^{\mbox{\scriptsize{right}}}\rangle$. In the right panel of figure 2, we show how the pair of approximate localized states, $|\Phi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Phi^{\mbox{\scriptsize{right}}}\rangle$, is by the symmetry restoration transformed into the pair of approximate eigenstates $|\Phi^{+}\rangle$ and $|\Phi^{-}\rangle$.

It is gratifying to see that the model energies of approximate localized states ${\cal E}^{\mbox{\scriptsize{loc}}}$ are only 1.6 keV higher than the average energies of the localized exact wave packets $E^{\mbox{\scriptsize{loc}}}$ (6). It is even more gratifying to see that energies of the symmetry-restored states, are only 1.8 and 1.7 keV above their exact counterparts. The energy splitting between the two symmetry-restored states, which equals 75.5 keV, almost exactly reproduces the exact result.

We should mention at this point that the exact localized wave packets, $|\Psi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Psi^{\mbox{\scriptsize{right}}}\rangle$, which represent non-stationary solutions of the Schrödinger equation in the doubly symmetric potential well, evolve in time in such a way that after the time of $T=\pi\hbar/\delta{E}$, the left wave packet will appear on the right-hand side and vice versa. For the selected parameters of the model, this left-right quantum oscillation time is very short, $T\simeq 3\times 10^{-21}$ s, and thus, a localized particle created in the left or right well will not really keep its identity. However, if the barrier width is increased from 2 to 45 fm, this oscillation time becomes $T\simeq 40$ days, and the particle created in one of the wells would remain there as a classical system would do. In the case of the wide barrier, the approximate localized states, $|\Phi^{\mbox{\scriptsize{left}}}\rangle$ and $|\Phi^{\mbox{\scriptsize{right}}}\rangle$, become extremely good representations of the localized exact wave packets, and the symmetry restoration becomes an extremely efficient method to obtain perfect approximations of the symmetry-conserving exact eigenstates.

The simple model discussed in this section shows that the symmetry breaking and restoration is a useful concept of describing the physical reality of quantum mechanics, and that it is not inherent to complicated many-body systems. Nevertheless, it is for these complicated and difficult systems that it finds its most prominent and successful applications.111Analyses performed in schematic models can be found, e.g., in Robledo (1992); Yannouleas and Landman (2002a, b). In particular, in nuclear and molecular physics, there is overwhelming evidence that symmetry-restored mean-field states provide for a global understanding of multiple phenomena and experimental observations.

At this point, to relate the symmetry restoration to the rigorous DFT Hohenberg and Kohn (1964); Kohn and Sham (1965); Bartlett (2010); Becke (2014); Jones (2015), a few comments are in order. The basis for existence theorems of exact DFT is the variational principle, whereby one reaches the exact ground-state of the system and its density. Within our simple example above, it would mean that the DFT is bound to yield the exact, symmetry conserving, positive-parity ground-state wave function $\Psi^{+}(x)=\langle{x}|\Psi^{+}\rangle$ and its density $\rho(x)=|\Psi^{+}(x)|^{2}$, and not the localized wave functions $\Psi^{\mbox{\scriptsize{left}}}(x)=\langle{x}|\Psi^{\mbox{\scriptsize{left}}}\rangle$ or $\Psi^{\mbox{\scriptsize{right}}}(x)=\langle{x}|\Psi^{\mbox{\scriptsize{right}}}\rangle$ and their respective densities.

However, it is obvious that densities (be they average or maximum) of exact and localized wave functions differ by about a factor of two, compare figures 1(a) and (b). This simple observation creates an important issue for systems, like nuclei, for which the equilibrium local density (the so-called saturation density of about 0.16 fm*-3*) is an important physical parameter determined by the nature of the underlying interaction. Indeed, for such systems we build (or derive) functionals that describe infinite saturated systems, which leads to the local density approximation, or finite self-bound systems, within a single potential well (typical for a drop of matter). Such functionals then have minimum energies at saturation density and thus can properly work only for localized wave functions and not for the exact symmetry-conserving ones.

To bring the discussion above away from the simple example, which we introduced only to illustrate basic concepts, and towards a realistic case, consider a positive-parity ground state of an alpha-particle emitting nucleus. Before the decay, the density of nucleons is almost constant within the nucleus and equal to the saturation density. After the decay, the exact parity-conserving wave function would correspond to a symmetric combination of a recoil nucleus moving right, with the alpha particle moving left, and that of the recoil moving left and the alpha moving right. It is obvious that such a state cannot be modeled by the same density functional as that used to model the nucleus before the decay, because densities are now twice smaller than the saturation density. However, it is also obvious that a symmetry-broken state, e.g., the one with the recoil moving left and alpha moving right, is entirely within the remit of that functional, because both subsystems do have similar local densities, not very different from the saturation density.

The case of the alpha-emitting nucleus illustrates crucial points of nuclear DFT, whereupon the symmetry breaking plays a fundamental role. It also tells us that the symmetry restoration is an equally fundamental piece of the description. Indeed, after modeling the DFT state of the recoil moving left and alpha moving right, we must symmetrize the obtained solution, because the alpha-particle detectors will, of course, never see any left-right asymmetry of the decay process. In this sense, the DFT description of many-body systems gives us immediate access to physical localized states describing specific configurations, by which we mean specific arrangements of constituents of composite objects. However, it is now clear that these configurations should never be confused with exact eigenstates, as they simply represent specific wave packets thereof, whereas a reasonable modeling of the exact eigenstates is then accessible via the symmetry restoration.

There remains, nevertheless, one troubling element of the link between the DFT and symmetry restoration. Indeed, to restore the symmetry, we need to have access not only to the average energies of the localized broken-symmetry states $|\Phi^{\mbox{\scriptsize{left}}}\rangle$ or $|\Phi^{\mbox{\scriptsize{right}}}\rangle$, which are within the remit of DFT, but also to the overlaps, $\epsilon=\langle\Phi^{\mbox{\scriptsize{left}}}|\Phi^{\mbox{\scriptsize{right}}}\rangle$, and matrix elements, $\Delta=\langle\Phi^{\mbox{\scriptsize{left}}}|\hat{H}|\Phi^{\mbox{\scriptsize{right}}}\rangle$, thereof, neither of which is. Within the nuclear-DFT applications, there is overwhelming evidence that $\epsilon$ and $\Delta$ can be evaluated using the corresponding Kohn-Sham states and generalized Wick’s theorem. This gives us a rich and reasonable description of numerous experimental data. However, such an approach constitutes a hybrid mix of the DFT and wave-function approaches and, to our knowledge, it has as yet no justification in any solid formalism. It appears that the many-body-physics community has executed a spectacular triple Axle jump into a pool without really verifying whether the water is there or not. Nevertheless, the obtained excellent results indicate that we may rather worry about finding a justification than about questioning the method itself.

The reader is begged to excuse us for the partly simplistic and partly philosophical narrative of this introductory section. We thought that exposing these basic notions could constitute a useful background of the following sections, where we move right on to the forefront description pertaining to the subject matter of this review. However, the advanced discussion that is coming up should not obscure the vision of the forest behind trees.

2.2 Dissociation of the natural molecular hydrogen and other similar

two-dimensional artificial dimers

The second illustrative example is taken from the areas of condensed-matter physics and chemistry, where the long-range interparticle Coulombic repulsion plays an essential role in inducing symmetry breaking at the level of mean-field treatments of interacting electrons. In this context, the simplest nontrivial examples involve at a minimum a pair of electrons. Here we will describe the case of an artificial two-dimensional (2D) quantum dot molecule (QDM), specifically a manmade system denoted as H2-QDM. The H2-QDM system consists of two electrons trapped inside two parabolic quantum dots (each with a harmonic-potential confinement specified by a frequency $\hbar\omega_{0}$) separated by an interdot distance $d$ and an interdot barrier $V_{b}$. The overall confinement potential is effectively that of a 2D two-center oscillator. Increasing the interdot separation $d$, or the interdot barrier $V_{b}$, generates a process bearing analogies to the dissociation of the natural H2 Hydrogen molecule. This process offers an immediate illustration of the symmetry dilemma facing the mean-field approaches. This expression [coined by Löwdin Lykos and Pratt (1963)] succinctly conveys the fact that imposing the Hamiltonian symmetries on the mean-field treatment [exemplified by the restricted Hartree Fock (RHF) method for electrons; see section 7], provides wave functions with the proper symmetries, but often the corresponding total energy is higher than that obtained when the symmetry requirements are relaxed [using the unrestricted Hartree Fock (UHF) method for electrons; see section 7]. In the double-dot example here, the UHF lowers the total energy by breaking in the ensuing wave functions the symmetries of total spin and of the parity along the separation axis.

The stretched natural H2 in both the RHF and the UHF mean-field levels and the ensuing correct-symmetry/higher-energy versus lower-energy/wrong-symmetry dilemma are described in detail in Chapter 3.8.7. of Szabo and Ostlund (1989). In wave-function-based approaches describing electronic systems, the symmetry dilemma can be overcome by using symmetry restoration.222In the Kohn-Sham DFT treatment of electronic systems, the symmetry dilemma remains an open question, with a main obstacle being the large self-interaction error; see, e.g., Perdew et al. (1995). An explicit illustration to this effect for the case of the artificial H2-QDM was provided in Yannouleas and Landman (2001, 2002b, 2002a). Before proceeding further, we mention two points that will be crucial in grasping the overall picture behind the mathematical formalism and the numerical details below. Namely, (I) The UHF equations do not provide a symmetry-broken solution in all circumstances. Given a specific many-body state, whether a ground state or an excited one, symmetry breaking appears at well-defined regions of the parameters characterizing the many-body electronic problem. In the absence of symmetry breaking, the UHF solutions coincide with the RHF ones, and thus subsequent symmetry restoration has no effect. Anticipating the specific numerical results below, we mention that, in the case of the singlet state of the H2-QDM, smaller values of Coulombic repulsion, interdot separation, and interdot barrier tend to suppress symmetry breaking. (II) Two electrons correspond to a closed electronic shell in both two and three dimensions. However, unlike the nuclear experience where shell closures prevent symmetry breaking (associated with nuclear shape deformations and nuclear pairing), the two-particle shell closure in electronic systems is not immune to symmetry breaking (associated with electron localization).

RHF and UHF numerical results for the singlet state of the H2-QDM are displayed in figure 3 for a case with trapping frequency of $\hbar\omega_{0}=5$ meV, interdot distance of $d=30$ nm and barrier of $V_{b}=4.95$ meV. In the RHF result (left column), a single bonding molecular orbital ($\sigma_{g}$) is occupied by both the spin-up and spin-down electrons. In the UHF case (right column), however, the total spin and parity symmetries are broken, and the spin-up electron occupies an atomic-like ($1s$) orbital located on the left quantum dot, while the spin-down electron occupies another different atomic-like ($1s^{\prime}$) orbital centered on the right quantum dot. For the total energies, the RHF result is $E_{\rm RHF}=13.68$ meV, while the UHF energy yields $E_{\rm UHF}=12.83$ meV, amounting to a lowering in energy of 0.85 meV. Furthermore, the energy of the triplet state is $E_{\rm UHF}=13.01$ meV, and thus the singlet state conforms to the requirement that for two electrons at zero magnetic field the singlet is always the ground state; in sharp contrast, the RHF molecular-orbital solution fails to fulfill this exact requirement.

Next we show how to construct an improved wave function by using the spin projection method to restore the broken symmetry of the UHF step. For clarity, we explicitly write down the UHF determinant, cf. equation (174),

[TABLE]

where $\varphi_{l}({\bm{r}})$ and $\varphi_{r}({\bm{r}})$ are the $1s$ (left) and $1s^{\prime}$ (right) localized orbitals of the UHF solution displayed in the right column of figure 3, and $\alpha$ and $\beta$ denote the up and down spins, respectively. We also found useful to introduce a shorthand notation for the determinant above as $\sqrt{2}\Phi^{S_{z}=0}_{\rm UHF}(1,2)\equiv|\varphi_{l}(1)\overline{\varphi}_{r}(2)\rangle$, where a bar over a space orbital denotes a spin-down electron; absence of a bar denotes a spin-up electron.

$\Phi^{S_{z}=0}_{\rm UHF}(1,2)$ is an eigenstate of the projection $S_{z}$ of the total spin ${\bm{S}}={\bm{s}}_{1}+{\bm{s}}_{2}$, but not of ${\bm{S}}^{2}$. A many-body wave function which is an eigenstate of ${\bm{S}}^{2}$ with eigenvalue $s(s+1)$ (here $s=0$) can be generated by applying the singlet-state spin projection operator, i.e.,

[TABLE]

where $\varpi_{12}$ is an operator that interchanges opposite spins of the two electrons; see the general formula in equation (180) below.

Upon the spin-symmetry restoration, the singlet state of two electrons (with $s=0$) is given by:

[TABLE]

We note that the beyond-mean-field projected many-body wave function (20) is a linear superposition of two Slater determinants, in contrast to the single-determinant wave functions of the RHF and UHF methods.

A further expansion of the determinants in equation (20) produces the equivalent expression

[TABLE]

where $\chi(s=0,S_{z}=0)$ is the spin eigenfunction for the singlet state, and is given by

[TABLE]

The expression in the right-hand side of equation (21) is similar to the Heitler-London Heitler and London (1927) or valence-bond wave function. However, at variance with the Heitler-London approach, which employs the left and right orbitals of the fully separated atoms, this expression involves the UHF orbitals that are optimized self-consistently at any separation $d$ and potential barrier height $V_{b}$.

Considering the normalization of the spatial part, one arrives at the following improved wave function of the singlet state:

[TABLE]

where the normalization constant is given by

[TABLE]

$S_{lr}$ being the overlap integral of the spatial orbitals $\varphi_{l}({\bm{r}})$ and $\varphi_{r}({\bm{r}})$,

[TABLE]

We stress again that the improved wave function $\Psi^{\rm s}_{\mbox{\scriptsize{PRJ}}}(1,2)$ exhibits all the symmetries of the original many-body Hamiltonian.

The total energy of the symmetry-restored singlet state $\Psi^{\rm s}_{\mbox{\scriptsize{PRJ}}}(1,2)$ [equation (23)] reads

[TABLE]

where $h$ is the single-particle contribution to the total Hamiltonian defined in section 7.2.1, and $J$ and $K$ are the direct and exchange matrix elements of the $e-e$ repulsion $e^{2}/\kappa r_{12}$, where $\kappa$ is the material dielectric constant. In equation (26), the subscripts $l$ and $r$ indicate that the matrix elements are associated with the left and right space orbitals, $\varphi_{l}({\bm{r}})$ and $\varphi_{r}({\bm{r}})$, respectively. For comparison, we provide also here the corresponding formula for the HF total energy, which is valid for both the RHF [with $\varphi_{l}({\bm{r}})=\varphi_{r}({\bm{r}})=\varphi({\bm{r}})$] and UHF [with either $\varphi_{l}({\bm{r}})\neq\varphi_{r}({\bm{r}})$ or $\varphi_{l}({\bm{r}})=\varphi_{r}({\bm{r}})=\varphi({\bm{r}})$] cases, i.e.,

[TABLE]

For the triplet state, the UHF solution (with $S_{z}=1$) does not break any symmetries for all set of values of the parameters $d$, $V_{b}$, and strength of the $e$-$e$ repulsion. For the triplet HF determinant, one has $\Phi^{\rm t}_{\rm RHF}(1,2)=\Phi^{\rm t}_{\rm UHF}(1,2)=|\varphi_{\rm 0n}({\bm{r}})\varphi_{\rm 1n}({\bm{r}})\rangle$, with the absence of bars over the $\varphi$ orbitals reflecting that both electrons have $\alpha$ (up) spins. The indices "0n" and "1n" denote zero-node and one-node space orbitals, respectively, and, for the set of parameters used here, they have the general form of molecular orbitals, similar to the $\Psi^{+}$ and $\Psi^{-}$ in figure 1(a); for the nodeless $\varphi_{\rm 0n}$, see also the MO orbital (modulus square) in the left column of figure 3. In this case, the spin-projection operator, $\hat{P}_{\rm spin}^{{\rm t},S_{z}=1}\equiv 1$, has no effect [no opposite spins to be interchanged, see the general formula in equation (180) below], and the projected wave function coincides with the starting HF determinant. As a result the energies for the triplet state in all three approximations are equal, i.e., $E^{\rm t}_{\mbox{\scriptsize{PRJ}}}=E^{\rm t}_{\mbox{\scriptsize{UHF}}}=E^{\rm t}_{\rm RHF}=E^{\rm t}_{\rm HF}$, with

[TABLE]

Like in the case of the natural H2, the energy difference, $\Delta\varepsilon=E^{\rm s}-E^{\rm t}$, between the singlet and the triplet states of the artificial H2-QDM, must vanish from below as the separation increases. Equivalently one can keep the interdot distance $d$ constant and vary the height of the interdot barrier $V_{b}$.

Figure 4 illustrates the evolution of $\Delta\varepsilon$ as a function of $V_{b}$ in zero magnetic field, and for all three successive approximation steps, i.e., the RHF (molecular-orbital theory, top line), the UHF (broken symmetry, middle line), and the projection approach (symmetry restoration, bottom line).

An inspection of figure 4 shows that both the UHF and projection-technique wave functions describe the energetics of the separation limit ($\Delta\varepsilon\rightarrow 0$ for $V_{b}\rightarrow\infty$) rather well, while the RHF approach fails. Furthermore, in this limit, the projected (symmetry restored) singlet wave function [see equations (20) and (21)] has the advantage of reproducing the fully spin-entangled two-fermion wave function introduced by David Bohm Bohm and Aharonov (1957); Bohm (1951) as a simpler example for demonstrating the Einstein-Podolsky-Rosen quantum-mechanical entanglement paradox Einstein et al. (1935). This exemplifies the potential for utilizing the symmetry restoration methodologies in the context of the emerging field of quantum information science Cloët (2019).

2.3 The seniority model

So far we discussed symmetry breaking and restoration for one or two particles. However, an essential point of the symmetry breaking discussed in this article is the approximate treatment of correlations in a many-body system by introducing the mean-field approximation, that is, by describing the many-body system in terms of a product state $|\Phi\rangle$ of uncorrelated particles (or quasi-particles) moving in a single-particle potential with a broken symmetry. To elaborate this aspect in more detail, we now briefly discuss the seniority model introduced by Kerman (1961) as an example. We consider $N$ fermions in a degenerate single $j$-shell (with $\Omega=j+\frac{1}{2}$) interacting through a monopole pairing force with the corresponding many-body Hamiltonian $\hat{H}$,

[TABLE]

where $G$ is the strength of the interaction and the operator,

[TABLE]

where $a_{j}^{{\dagger}}$ creates a single particle in the j-th shell and $\hat{S}_{+}$ creates a Cooper-pair of particles coupled to angular momentum $J=0$. Together with the operator $\hat{S}_{0}$ connected with the particle-number operator $\hat{N}$ by the relation

[TABLE]

operators $\hat{S}_{+},\hat{S}_{0},$ and $\hat{S}_{-}$ form the algebra of the group $SU(2)$ of the quasi-spin. It has the Casimir operator $\mathbf{\hat{S}}^{2}$,

[TABLE]

with the eigenvalue $S(S+1)$. The Hamiltonian (29) can be expressed as

[TABLE]

and thus it is diagonal in the basis characterized by the quantum numbers $S,S_{0}$ or $S,N$. The particle vacuum $|-\rangle$ with $N=0$ is given by $S=\frac{\Omega}{2}$ and $S_{0}=-\frac{\Omega}{2}$. Starting from this vacuum $|-\rangle$, and applying the raising operator $\hat{S}_{+}$, one finds the exact ground states of the system with even particle numbers $N$ Hara (1967),

[TABLE]

This is a condensate of $n=\frac{N}{2}$ Cooper pairs. In terms of the original fermions this is a highly correlated state.

We now use the mean-field approximation, i.e., we start with the BCS-state

[TABLE]

where $\alpha_{jm}^{{\dagger}}=u_{m}a_{jm}^{{\dagger}}-v_{m}a_{j\bar{m}}$. $\bar{m}$ denotes the time-reversed state. $|\Phi\rangle$ has the form of a generalized product state and it is the vacuum for the quasiparticle operators $\alpha_{jm}^{{\dagger}}$. We, therefore, can apply the Wick theorem (see equation (7.47) of Ref. Ring and Schuck (1980)) for the two-body operator (29). Neglecting, as usual in the BCS approximation, the terms $a^{\dagger}a$ we obtain the BCS Hamiltonian

[TABLE]

with the gap parameter $\Delta=G\langle\Phi|S_{+}|\Phi\rangle$. The BCS state $|\Phi\rangle$ is an eigenstate of $\hat{H}_{BCS}$ and the BCS coefficients $u_{m}$ and $v_{m}$ are determined by the diagonalization of the mean-field Hamiltonian $\hat{H}_{BCS}$.

It is evident that $\hat{H}_{BCS}$ breaks the particle-number symmetry; therefore one looks for eigenfunctions of $\hat{H}_{BCS}^{\prime}=\hat{H}_{BCS}-\lambda\hat{N}$, where the Lagrange parameter $\lambda$ is determined by the subsidiary condition $\langle\Phi|\hat{N}|\Phi\rangle=N$. Since the BCS Hamiltonian depends on the gap parameter $\Delta$, one ends up with a non-linear problem, which has to be solved iteratively. In the seniority model discussed here, the single-particle part (${\sim}a^{{\dagger}}a$) of $\hat{H}_{BCS}$ vanishes and thus does not depend on the quantum number $m$, and hence the BCS amplitudes $u_{m}$ and $v_{m}$ do not depend on $m$ either.

The BCS Hamiltonian $\hat{H}_{BCS}$ breaks the particle-number symmetry, and so it also breaks the gauge symmetry, see section 3.2.1 below. Through this symmetry breaking, it is possible to represent the wave function as a product state in terms of quasiparticles, cf. equation (35). In the exact solution, these quasiparticles are not independent and there are additional correlations, which are not taken into account in the product state. However, we can bring these correlations back by the restoration of the symmetry. This is achieved by the particle-number projection, that is, by neglecting in the wave function $|\Phi\rangle$ all the contributions with particle numbers different from $N$. For this purpose we express the BCS state $|\Phi\rangle$ (35) in terms of the operator $\hat{S}_{+}$ as

[TABLE]

where we have used the fact that the squares of fermion creation operators are equal to zero and thus only the first two terms of the exponential remain, as in equation (35). The particle-number projection then leads to

[TABLE]

We find that the restoration of the symmetry leads to the exact solution (34) of the system. Of course, this is a very specific model and the fact, that symmetry restoration brings us back to the exact solution, depends definitely on the fact, that the operator $\hat{S}_{+}$ does not depend on the particle number nor on other properties of the model such as the strength parameter $G$. Nonetheless, also in more general cases, where the symmetry restoration does not lead to the exact solution, we will find that by restoring symmetries, one can improve the mean-field approximation considerably.

3 Symmetry Restoration - General Formalism

As discussed in the previous section, mean-field solutions may break symmetries that the original many-body Hamiltonian obeys. For strong symmetry breaking, approximate methods can be used to evaluate the observable quantities to a good accuracy, and these are discussed in section 4. For weak symmetry breaking, wave functions defined in the intrinsic frame of reference should have their broken symmetries restored. In the sixties and seventies of the last century, considerable efforts were made to decouple the total Hamiltonian in terms of intrinsic and collective degrees of freedom Lipkin et al. (1955); Lipkin (1956); Villars (1957); Lipkin (1958); Villars and Cooper (1970). Owing to the Galilean invariance, this method is well defined only for the simplest case of the linear-momentum, where the center of mass coordinate separates exactly from the intrinsic degrees of freedom, described, for instance, by Jacobi coordinates. This is used for few-body systems, i.e., in very light nuclei Navratil et al. (2000). For heavier nuclei this technique becomes quite cumbersome because the wave functions are far from being product states in these coordinates. Only recently a method based on Quantum Monte Carlo techniques has been introduced, which allows the exact treatment of Galilean invariance even for density functionals with finite range Massella et al. (2020). Because of its numerical complexity, however, applications are also limited to light nuclei and it cannot be generalized to other symmetry violations, where Galilean invariance does not apply.

Moreover, the popular terminology referring to the intrinsic and laboratory reference frames can be either confusing, or useless, or both. Indeed, in the language of the symmetry restoration, no reference frame is ever changed, namely, both the broken-symmetry and restored-symmetry wave functions reside in the same unique Hilbert space with one unique reference frame conveniently predefined and used. Then, the restored-symmetry wave function is obtained from the broken-symmetry one by acting on it with a specific symmetry-restoration operator, which does not change any reference frame either, but rather rotates the broken-symmetry wave function in the predefined reference frame.

In this article, we use the notion of an active rotation, whereby not the reference frame but the states are rotated. However, even if we used a passive rotation scheme, whereby the wave function stays the same and the reference frame is rotated, this would not have been equivalent to any change of the reference frame from intrinsic to laboratory. Although we may occasionally slip into the traditional terminology of the intrinsic and laboratory frames, the reader should always follow the correct description by translating the term "wave function in the intrinsic (laboratory) frame" as "broken-symmetry (restored-symmetry) wave function".

For a continuous symmetry group, an infinite number of solutions, which are degenerate in energy, are obtained by applying the elements of the Lie group on the intrinsic state. For instance, all Nilsson intrinsic states that differ only by rotation in space have the same energy. Peierls and Yoccoz (1957) employed a linear superposition of these degenerate states with the weight functions determined through a second variational procedure. In this way, the variational procedure is performed in two steps. In the first one, the intrinsic state is determined, and in the second one, the energy is minimized within the subspace of states, invariant under the symmetry group, which can be projected from that intrinsic state. This double variational approach is referred to as projection after variation.

In a single-step variational procedure, the variation is performed by simultaneously considering symmetry-restored trial states that can be obtained from any intrinsic state. In the case of rotational symmetry, this corresponds to first projecting the intrinsic wave function onto a state with well-defined angular-momentum and then performing the variational procedure over the intrinsic states. This variation after projection method, proposed originally by Zeh Zeh (1965), for a simple degenerate model leads to practically exact results Sheikh and Ring (2000). The difference between two-step and single-step variational approaches becomes quite obvious for the case of particle permutational symmetry, with the method of Peierls and Yoccoz giving rise to an approximate anti-symmetrized Hartree solution, while the method of Zeh leading to the correct HF solution.

Projection methods developed to restore the symmetries can be divided into those pertaining to abelian and non-abelian symmetry groups Löwdin (1967); MacDonald (1970); Ring and Schuck (1980). Restoration of the linear-momentum symmetry, gauge symmetry associated with the particle number, and parity symmetry pertain to the abelian groups. For this class, projection operators have expected mathematical properties of idempotency and hermiticity. For the class of non-abelian symmetry groups, which includes three-dimensional rotations, and corresponds to restoration of angular-momentum or isospin, the projection operators do not have properties of idempotency and hermiticity. However, as we discuss below, such operators project out the relevant quantum numbers from the intrinsic wave function, and that is what is important for physical applications.

In the following, we first construct the projection operators using group-theory and generator-coordinate methods. As is evident, all the projection operators can be expressed as basis transformations in some representation. Projection methods for various symmetries, such as linear momentum, three-dimensional angular momentum, particle number, and parity are discussed in section 3.2.

The projection formalism using the generalized HFB basis is discussed in section 3.3. In this section, the expressions for the norm and the Hamiltonian kernels between the HFB transformed and the initial basis are obtained using the generalized Wick’s theorem. Methods for performing variation after projection are then discussed in this section. In particular, it is demonstrated that variation after projection of an arbitrary symmetry-projected energy leads to HFB-like equations, as is the case for the unprojected energy, with the only difference that the pairing and HF fields get modified and depend on the projected quantum numbers as well.

3.1 Projection Operator - Mathematical Basis

The mathematical structure of the projection operators is constructed in this subsection using the group-theory approach and the generator coordinate method.

3.1.1 Group-Theory Approach

Given a symmetry group, a projection operator, $\hat{P}^{\lambda}$ can be constructed such that for an arbitrary wave function $|\Phi\rangle$, $\hat{P}^{\lambda}|\Phi\rangle$ is the component of $|\Phi\rangle$ belonging to the irreducible representation characterized by the quantum number(s) $\lambda$. In the following, we construct the mathematical structure of such an operator using the group-theory approach and closely follow the textbooks of Hamermesh Hamermesh (1962) and Gilmore Gilmore (2008a).

Let us suppose that the states $|\lambda\mu\rangle$, enumerated by quantum number(s) $\mu$, span an irreducible representation of the group defined by the quantum number(s) $\lambda$. Within the irreducible representation, the group transformations $\{{\hat{R}}(g)\}$ act in the following way

[TABLE]

where $g$ denotes the group element. The matrix functions $D^{\lambda}_{\mu\nu}(g)$ are the continuous single-valued representations of the group, which for the group of three-dimensional rotations are called Wigner $D$-functions Edmonds (1957); Varshalovich et al. (1988). They obey the orthogonality theorem of representation theory

[TABLE]

where $n_{I}$ denotes the dimension of the irreducible representation.

For finite groups, $V$ is the order of the group and the sum in equation (40) runs over all group elements, whereas for continuous compact groups, $V$ is the volume of the parameter space of the group, and the sum over the group elements should be replaced by an integral over the group parameters. For example, for the rotational group in three dimensions, where the group elements depend on three Euler angles $\Omega=(\alpha,\beta,\gamma)$, the volume is $V=\int d\Omega=8\pi^{2}$ and the dimension of the representation characterized by the angular momentum $I$ is $n_{I}=2I+1$. Moreover, for continuous non-compact groups (like the translational group), representations are labeled by continuous parameters and the Kronecker deltas must be replaced by Dirac deltas.

Multiplying equation (39) by $D^{\lambda\ast}_{\mu\nu}(g)$ and summing over the group elements, we have

[TABLE]

This allows us to define operators $\hat{P}^{\lambda}_{\mu\nu}$,

[TABLE]

such that

[TABLE]

and

[TABLE]

It is evident from the above equation that the diagonal operators, $\hat{P}^{\lambda}_{\mu\mu}$, project out the $\mu$-th columns of the $\lambda$-th irreducible representation, and obey

[TABLE]

It means that they are hermitian and idempotent, i.e., projection operators in the strict mathematical sense.

For non-abelian groups, the wave functions $|\Phi^{\lambda}_{\mu\mu}\rangle\equiv\hat{P}^{\lambda}_{\mu\mu}|\Phi\rangle$ are no longer tensors with respect to this group, i.e., under the operation $\hat{R}(g)$ they do not behave as wave functions $|\lambda\mu\rangle$ shown in equation (39). However, owing to the completeness of the matrix functions $D^{\lambda}_{\mu\nu}(g)$ in the space of regular functions of $g$, operators $\hat{P}^{\lambda}_{\mu\mu}$ obey the following important relation,

[TABLE]

that is, they provide for the so-called resolution of unity. This means, that every wave function $|\Phi\rangle$ can be uniquely split into a sum of components $|\Phi^{\lambda}_{\mu\mu}\rangle$, $|\Phi\rangle=\sum_{\lambda\mu}|\Phi^{\lambda}_{\mu\mu}\rangle$, even if these components are not necessarily orthogonal or normalized, or even if some of them can vanish.

3.1.2 Generator Coordinate Method

The expression for the projection operator in equation (42) can also be obtained as a special case of the generator coordinate method Hill and Wheeler (1953); Griffin and Wheeler (1957); Ring and Schuck (1980). In this method, Peierls and Yoccoz (1957) started from the symmetry breaking wave function $|\Phi\rangle$ and diagonalized the many-body Hamiltonian in the collective subspace spanned by a set of generating functions

[TABLE]

Using the ansatz

[TABLE]

the variation of the energy expectation value with respect to the weight function $f(g)$ leads to a generalized eigenvalue problem with the generating functions as non-orthogonal basis states. To restore the symmetry, the weight function is chosen in such a way that the resulting many-body wave function transforms as in equation (39). It can be shown that the collective subspace is invariant under the symmetry transformations, i.e.,

[TABLE]

where $(-g+g^{\prime})$ is a short hand notation for the group element $\hat{R}^{-1}(g)\hat{R}(g^{\prime})$. This implies that the projector $\hat{P}$ onto this subspace commutes with the symmetry operator $\hat{R}(g)$ and it is possible to find simultaneous eigenstates of $\hat{P}\hat{H}\hat{P}$ and the symmetry operator. Thus, a function $f(g)$ exists that not only minimizes the energy but also has the proper symmetry. This function $f(g)$ can be found by expanding it in terms of the representations of the group, characterized by eigenvalues of the Casimir operators.

3.1.3 Generic properties of projection

Over the years projection techniques were introduced for various symmetry groups. In the following, we use a generic form of the projection operator (42),

[TABLE]

Assuming that the Hamiltonian $\hat{H}$ commutes with the symmetry operator $\hat{R}(g)$, the projected energy is given by

[TABLE]

or equivalently as Sheikh and Ring (2000)

[TABLE]

where the rotated wave function $|g\rangle$ is defined as a suitably normalized generating function (47),

[TABLE]

and two auxiliary functions $x(g)$ and $y(g)$ are defined as

[TABLE]

Similar expressions can be found for other operators, such as the multipole operators and for electromagnetic transition matrix elements Dobaczewski et al. (2009).

For the one-dimensional case, $I$ represents directly the quantum number on which one projects. In the general case, $I$ represents several quantum numbers, as for instance $K_{x},K_{y},K_{z}$ for projection onto the linear momentum or $I,M,K$ in the case of three-dimensional rotations.

3.2 Projection Methods for Various Symmetries

In Tables 1–2, we show examples of symmetry restoration corresponding to several symmetry groups, with five specific cases discussed in sections 3.2.1–3.2.5 in more detail. In addition, the isospin symmetry restoration being formally identical to the angular-momentum restoration, sections 3.2.4 and 3.2.5, is not explicitly covered in this section, apart from listing the relevant expressions in Tables 1 and 2.

3.2.1 Particle-number

Rotations in gauge space are defined by

[TABLE]

where $\hat{N}$ corresponds to the number operator for neutrons or protons. Wave functions with good particle number are only multiplied with a phase. This is no longer true in the case of HFB- or BCS-wave functions, being linear combinations of states with different particle numbers. In this case one has to restore the symmetry by projection onto a good particle number. This is an abelian symmetry group and from Table 2 we find the corresponding operator

[TABLE]

Particle-number projection, in particular in the BCS-case is relatively simple. It was applied using different methods Bayman (1960); Dietrich et al. (1964); Fomenko (1970); Janssen and Schuck (1981); Egido and Ring (1982a), in particular, also in the framework of approximate projections, see section 4.

3.2.2 Linear Momentum

In the nuclear theory, localized single-particle states are employed that are not eigenstates of the momentum operator. Conversely, the plane waves that are eigenstates of the momentum operator cannot describe a localized system of particles.

Apart from very light nuclei, in most of the practical applications the nuclear wave functions are based in one way or another, on wave functions in localized potentials. They violate the translational invariance of the underlying Hamiltonian. As shown in Tables 1 and 2, this is a three-dimensional group with the three operators $\hat{{\bm{P}}}=\sum_{i}{\bm{p}}_{i}$ of the total momentum as generators.

In principle, one should carry out a projection onto good momentum Peierls and Yoccoz (1957). This is technically rather difficult and therefore, there are very few examples in the literature, where such an exact projection was carried out Schmid (2001, 2002a, 2002b, 2003); Rodríguez-Guzmán and Schmid (2004a). In most of the applications approximate methods are used. As discussed in section 4, it can be shown that the corrections introduced in this way decrease with numbers of particles, $A$ as $1/A$ and therefore in most of the applications conservation of the linear momentum is only taken into account approximately (see section 4).

3.2.3 Parity

The parity projection operator is connected to a discrete symmetry. It is similar to the projection of spin singlet and triplet states (19), and it can be written as Egido and Robledo (1991),

[TABLE]

where $p=\pm 1$ and $\hat{\Pi}=\exp{(-i\pi\hat{N}_{-})}$ is the standard parity (inversion) operator. Here $\hat{N}_{-}=\sum_{k}^{{}^{\prime}}a^{\dagger}_{k}a_{k}$ is a restricted summation over all states $k$ with negative parity. We find $\hat{\Pi}^{2}=1$ and $\hat{P}^{p}$ is a true projection operator.

3.2.4 One-dimensional rotation

The one-dimensional rotation by the angle $\alpha$ around the $z$-axis is given by

[TABLE]

where $\hat{J}_{z}$ denotes the $z$ component of the angular-momentum operator.333In this review, we use dimensionless linear-momentum, angular-momentum, and spin operators, which represent generators of the corresponding symmetry groups. Whenever the corresponding quantum-mechanical operators and/or their eigenvalues are considered, we tacitly assume that the Planck-constant prefactors $\hbar$ are also included. The irreducible representations of this group are given by

[TABLE]

where $M$ denotes the eigenvalues of $\hat{J}_{z}$. Expressing function $f(\alpha)$ of equation (48) in terms of these irreducible representations is equivalent to a Fourier transformation,

[TABLE]

Using the $M$-th component $D^{M\ast}(\alpha)$ as the weight function in equation (48), we obtain

[TABLE]

with the projector

[TABLE]

This is in full agreement with the general expression (42) for the projection operator derived in section 3.1.1 from group-theory considerations. We note here, that the one-dimensional symmetry restoration projects deformed states on good projections $M$ of angular momentum $\hat{J}_{z}$, but each projected states is still a mixture of components with good total angular momenta $I\geq{M}$.

3.2.5 Three-dimensional rotation

As an example of a non-abelian group, we consider rotations in three dimensions, cf. recent comprehensive review in Ref. Bally and Bender (2021). They are characterized by the Euler angles $\Omega=(\alpha,\beta,\gamma)$. The corresponding group element is given by the operator

[TABLE]

and we have to introduce the generalized projection operators (42):

[TABLE]

where

[TABLE]

are Wigner $D$-functions Edmonds (1957); Varshalovich et al. (1988). Following Lamme and Boeker (1968), we define a complete and orthogonal set of many-body wave functions $|IMi\rangle$ which are eigenstates of the angular-momentum operators $\mathbf{\hat{J}}^{2}$ and $\mathbf{\hat{J}}_{z}$, and $i$ combines all the remaining quantum numbers in the many-body Hilbert space. Using equation (39), we have

[TABLE]

and owing to the completeness relations for the states $|IMi\rangle$, we can express the generalized projectors (64) as

[TABLE]

This again shows that only the diagonal term $\hat{P}^{I}_{MM}$ is a true projector onto the sub-space of the Hilbert space with the quantum numbers $I$ and $M$. However, individual states $\hat{P}^{I}_{MM}|\Phi\rangle$ cannot be identified with the basis states $|IMi\rangle$, i.e., they do not obey equation (66).

For a better understanding of the additional444It is better to call $K$ additional quantum number than to call it projection of angular momentum on the intrinsic axis, which is the term frequently used. Indeed, states $\langle IKi|$ are “bra” representations of “ket” states $|IMi\rangle$, and both correspond to projections of angular momentum on the same predefined quantization axis. quantum number $K$ in the projector $\hat{P}^{I}_{MK}$, we start from a deformed intrinsic (symmetry-breaking) wave function $|\Phi\rangle$, and, in analogy to equation(69) we expand the weight function $f(\Omega)$ in equation (48) in terms of the complete set of Wigner functions

[TABLE]

If we restrict us in this sum to the terms with fixed $I$ and $M$, we obtain, by using equations (48) and (64), the projected many-body state

[TABLE]

From (67) it is evident that $|\Psi^{I}_{M}\rangle$ has good quantum numbers $I$ and $M$ for arbitrary expansion coefficients $g^{I}_{MK}$ (for further details see Zeh (1967)).

In contrast to abelian groups, here the weight function $f^{I}_{M}(\Omega)$ for the generator-coordinate-method ansatz (48) is not completely defined by the symmetry group. The coefficients $g^{I}_{MK}$ have to be determined by the dynamics of the system, i.e., by diagonalizing the many-body Hamiltonian or by minimizing the projected energy.

Only in special cases, this additional diagonalization is not necessary, e.g., if the intrinsic (symmetry-breaking) wave function $|\Phi\rangle$ is symmetric with respect to rotations around the intrinsic $z$-axis, i.e., if $\hat{J}_{z}|\Phi\rangle=K_{0}|\Phi\rangle$, then $\hat{P}^{I}_{MK}|\Phi\rangle=0$ for $K\neq K_{0}$ and there is only one coefficient $g^{I}_{MK_{0}}$, which is determined by the normalization. A simple case is the intrinsic state $|\Phi\rangle_{K=0}$ of the ground state of an axially deformed even-even nucleus. Here we find the projected states

[TABLE]

which do obey equation (66),

Further, by integrating in equation (64) over the Euler angles $\alpha$ and $\gamma$, and by using equation (62), the generalized projector $\hat{P}^{I}_{MK}$ can be decomposed into three steps:

[TABLE]

We can now describe this result in two different ways:

In the traditional language of passive transformations between the intrinsic and the laboratory reference frames, we start with the projector, $P^{K}=\delta(\hat{J}_{z}-K)$, onto the quantum number $K$, corresponding to the component of the angular momentum $\hat{\bm{J}}$ along the intrinsic $z$-axis, then we have a rotation by angle $\beta$ around the $y$-axis from the $z$ axis in the intrinsic system to the $z$-axis in the laboratory frame, and finally $P^{M}=\delta(\hat{J}_{z}-M)$ projects on an eigenstate with quantum number $M$ in the laboratory frame. 2. 2.

In the language of symmetry-broken and symmetry-restored states, the first operation projects out the symmetry-broken state on the good quantum number $K$ corresponding to the $z$ quantization axis of a predefined reference frame. The second step projects on the good total-angular-momentum quantum number $I$, but it does that by a rotation along the $y$ axis, and thus mixes again the previously restored projections of the angular momentum. Then, the third-step projector is required to restore the projection $M$ on the $z$ quantization axis of a predefined reference frame, and to give the fully symmetry-restored wave function.

The final wave function $|\Psi^{I}_{M}\rangle$ in equation (68) is a quantum-mechanical superposition of all these different orientations and, as usual in the generator coordinate method, the weight functions are related to the corresponding probability amplitudes Ring and Schuck (1980).

It is important to emphasize, that the concept of generator coordinates, which corresponds here to the projection onto the subspaces determined by the symmetry group, deals only with the coordinates of the $A$ particles in the corresponding wave functions in the intrinsic or in the laboratory frame, i.e., before or after the symmetry restoration. The collective coordinates, in this case the Euler angles, enter only in a parametric way. In none of these considerations one has to introduce "redundant" coordinates and no spurious states are involved. One stays, from the beginning to end, completely in the quantum-mechanical framework and no "requantization" is necessary. However, the results are not characterized by orbits in the collective subspace, but by probabilities corresponding to different orientations.

It needs to be added that there is a prize to pay as these calculations become relatively complicated. Therefore, although these concepts were around since more than half a century, many of these calculations, in particular, those with realistic applications are possible only nowadays using modern computing resources Bally et al. (2014). Some of the applications still have to wait for more advanced implementations Romero et al. (2019).

3.3 Symmetry restoration of the HFB wave function

In most of the projection studies, one starts with a broken-symmetry mean-field wave function $|\Phi\rangle$ of the HFB type Ring and Schuck (1980); Schunck (2019). These states can be expressed as vacua of quasiparticles $\alpha_{k}|\Phi\rangle=0$ which are connected by a unitary transformation to the operators $a^{\dagger}_{n},a_{n}$ of an arbitrary reference basis

[TABLE]

There is a one-to-one correspondence between the HFB wavefunction and the corresponding one-body densities,

[TABLE]

and Wick’s theorem allows us to evaluate the matrix elements of arbitrary many-body operators $\langle\Phi|\hat{O}|\Phi\rangle$ in terms of these densities, see Appendix A. In the conventional HFB theory, wave function $|\Phi\rangle$, i.e., the quasiparticle amplitudes $U_{nk}$ and $V_{nk}$, are determined by a variation of the energy $E=\langle\Phi|\hat{H}|\Phi\rangle$ with respect to the densities. This leads to two single-particle fields

[TABLE]

and the corresponding HFB-equations

[TABLE]

where ${\cal H}^{\prime}={\cal H}-\lambda{\cal N}$, $E$ is the diagonal matrix of quasiparticle energies and the chemical potential $\lambda$ is determined by adjusting the average particle number.

3.3.1 Rotated Norm and Energy Kernels

Standard expressions for the norm overlap are given in Appendix A, whereas those that show the explicit dependence on the densities $\rho$ and $\kappa$ are given in Sheikh and Ring (2000):

[TABLE]

Here $R_{g}$ is the matrix representing the group element $g$ in the original basis

[TABLE]

For $|g\rangle$ given by equation (53), the generalized Wick’s theorem allows us to express matrix elements of the form $\langle\Phi|\hat{O}|g\rangle$ in terms of the corresponding transition densities Onishi and Yoshida (1966); Balian and Brezin (1969); Hara and Iwasaki (1979); Ring and Schuck (1980):

[TABLE]

It is important to note that here all the matrix elements are expressed in terms of the intrinsic densities $\rho$ and $\kappa$ and the matrix representation $R_{g}$ of the group element $g$.

As an example, let us consider a Hamiltonian with a two-body interaction of the form

[TABLE]

for which we obtain the transition matrix element:

[TABLE]

with the rotated fields

[TABLE]

In principle, the evaluation of projected matrix elements is relatively straightforward. One only has to replace the normal density matrices $\rho$, $\kappa$, and $\kappa^{*}$ by the transition densities $\rho(g)$, $\kappa(g)$, and $\overline{\kappa}^{*}(g)$ and integrate over the parameter space of the group. In practice, however, depending on the dimension of the single-particle space and the number of mesh-points in parameter space this can require a large computational effort, in particular, for triaxial nuclei. At each point $g$ in parameter space one has to invert a large (often complex) matrix $A_{g}$ (76) with the dimension of the single-particle space [see for instance Yao et al. (2014)].

As the Madrid group Anguiano et al. (2001) showed, it may happen in regions of level-crossings that $\langle\Phi|\hat{R}(g)|\Phi\rangle$ vanishes at certain values of $g$. This leads to poles in certain parts of the Hamiltonian matrix element (81), e.g., in ${\rm Tr}\left(\Gamma_{g}\rho_{g}\right)$, see discussion in section 6.1.2.

3.3.2 Variation after projection

With the techniques discussed in the previous section, it is relatively simple to carry out a projection after variation. However, such a procedure is not variationally optimal. Therefore, the method of variation after projection was proposed Zeh (1965); Yoccoz (1966); Rouhaninejad and Yoccoz (1966), where the mean-field wave function $|\Phi\rangle$ is determined by minimizing the projected energy, i.e., by solving the equation:

[TABLE]

For the exact solution of this problem, the following two methods were proposed:

Gradient Method

A particularly powerful method to minimize the projected energy with respect to the product state is the gradient method, which was introduced in Mang et al. (1976), and which was applied for variation-after-projection calculations in Egido and Ring (1982a, b). In this method, in the neighborhood of an arbitrary point $|\Phi_{0}\rangle$, the manifold of the HFB wave functions $|\Phi\rangle$ is parameterized by the Thouless theorem,

[TABLE]

where operators $\alpha_{k}^{{\dagger}}$ are the quasiparticle operators with respect to the quasiparticle vacuum $|\Phi_{0}\rangle$, i.e., $\alpha_{k}|\Phi_{0}\rangle=0$. The gradient of the projected energy with respect to parameters $Z_{kk^{\prime}}$ is given by

[TABLE]

These matrix elements can be evaluated using the generalized Wick’s theorem, in a similar way as it is done for the average energy. Following the direction of steepest descent on this manifold by a step size of $\eta$, we obtain in the next step of the iteration, the following wave function

[TABLE]

The resulting HFB coefficients of $|\Phi_{1}\rangle$ have to be orthogonalized [for details see Egido and Ring (1982a)]. Calculating the new projected energy $E_{1}^{I}$ and changing the step size $\eta$ in the next step accordingly, the minimum of the projected energy surface can be found without diagonalizing any matrix. This method is particularly useful, if one wants to minimize the energy surface with additional constraints. In this case the method of Lagrange multipliers is used, where the total gradient is projected onto the gradient along the hyper-surface determined by the constraining operator. The speed of convergence of this method can be considerably improved by using the conjugate gradient method Egido et al. (1995).

It is evident that the gradient method can only be applied for cases, where the final solution corresponds to a minimum in the energy surface. In all the applications of the Covariant Density Functional Theory, see section 6.2.2, because of the "no-sea approximation", the solution of the corresponding mean-field equations do not correspond to a minimum, but rather to a saddle point on the energy surface. In this case, the projected HFB equations are more useful.

Projected HFB Equations

In a similar way as the normal HFB-equations in equation (74) were derived by a variation of the unprojected energy $E=\langle\Phi|\hat{H}|\Phi\rangle$ with respect to the intrinsic densities $\rho$ and $\kappa$, projected HFB-equations were derived in Sheikh and Ring (2000) by a variation of the projected energy $E^{I}$ in equation (52). These equations have the same form as the unprojected HFB equations (74). However, in this case the HF and pairing fields depend on the quantum number $I$:

[TABLE]

In order to write down these quantities explicitly, we need the analytic form of the projected $E^{I}$ in equations (52) and (81) in terms of the intrinsic densities $\rho$ and $\kappa$ as given in equations (76)–(79). The detailed expressions for $\hat{h}^{I}$ and $\hat{\Delta}^{I}$ were derived in Sheikh et al. (2002) and are not repeated here.

It needs to be mentioned that several other methods were developed to perform the particle-number projection Ring and Schuck (1980) in the BCS case. In particular, the method of residuum integrals was introduced to perform the exact particle-number projection Dietrich et al. (1964) before the variation. In this approach, the particle-number projected method is cast into a set of non-linear equations, which are similar in structure to those of BCS equations.

4 Approximate projection methods

The general formalism of symmetry restoration, presented in the previous section, leads to expressions that involve multi-dimensional integrals of norm overlaps and Hamiltonian transition matrix elements over the symmetry group parameters, see equation (51). In the standard numerical approach these integrals are replaced with finite sums,555An alternative method based on solving linear equations was recently proposed in Johnson and O’Mara (2017); Johnson and Jiao (2019). which entails evaluations of the norm overlaps and Hamiltonian transition matrix elements for each mesh point in the multi-dimensional space of group parameters. For each mesh point, the norm overlap is evaluated using the Onishi or Pfaffian method and the Hamiltonian transition matrix elements are expressed in terms of the one-body rotated densities by employing the generalized Wick’s theorem. In general, the numerical cost of calculating one point of the projected integrand is somewhat larger than the cost of performing one iteration of the self-consistent method, required to determine the broken-symmetry state being projected.

Estimates of the number of mesh points vary a lot depending on the mass or deformation of the nucleus, and on the precision desired for the final results. On the one hand, for a strong symmetry breaking case, a larger number of eigenstates of the symmetry generators are contained within the broken-symmetry state, and thus more mesh points are required to resolve them. On the other hand, conserved symmetries allow us to limit the integration domain and thus to decrease the number of mesh points.

Typically, for a medium heavy nucleus with moderate deformation, about 10 integration points are needed to project on good proton and neutron numbers, and about 50 integration points, to project on good angular momentum of an axial nucleus. Restoration of symmetry then requires a numerical expense somewhat larger than that required to perform 5000 self-consistent iterations. It is evident that even with such a modest mesh size, the restoration of symmetries lead to numerical cost largely exceeding the typical 100 iterations required to converge the broken-symmetry state itself. For a triaxial state, where a three-dimensional integration over the Euler angles is needed, the number of integration points would increase to 12,500,000 and thus would become unmanageable. Therefore, up to now, calculations of this type were restricted to lighter systems only, see Bally et al. (2014).

In the early days of projection theory, well before the above tour-de-force achievements were envisioned and when the adequate computing power was not yet available, several approximate methods for symmetry restoration were proposed and implemented. Here we discuss in detail the most popular one, based on the so-called Lipkin method Lipkin (1960) or Kamlah expansion Kamlah (1968), along with the variant of the former one proposed by Nogami Nogami (1964). In fact, the Lipkin and Kamlah ideas were basically identical, although Kamlah did not apparently know about, and he did not cite the much earlier work of Lipkin.666However, Kamlah did cite Goodfellow and Nogami (1966) who had cited Lipkin (1960).

The main objective of both approaches is to obtain approximate expressions for symmetry-projected energies and to employ them in the implementation of the variation-after-projection approach. However, the two approaches differ on the physical quantities to be described: Kamlah is primarily concerned with determining the projected energy, which is then varied; whereas Lipkin aims to model the entire spectrum of collective states related to a given broken symmetry. Another difference between them is that Kamlah does and Lipkin does not consider the effects of collective motion brought about by the so-called pushing or cranking terms.

The baseline of the Lipkin and Kamlah approaches is the observation that the average values and matrix elements of operators calculated between symmetry-projected states always involve transition matrix elements of operators between symmetry-transformed states, equation (81). It is thus obvious that a meaningful approximation of the latter may lead to a useful approximation of the former.

Finally, we should stress the fact that although the Lipkin and Kamlah approaches are primarily concerned with identifying the variation-after-projection symmetry-breaking state, they do not actually determine the projected state. For that, an explicit projection of the variation-after-projection symmetry-breaking state is always necessary. Only then one can calculate correct transition probabilities respecting all symmetry properties of the transition operators.

4.1 The Lipkin method

The main idea of the Lipkin approach Lipkin (1960) can be formulated as a proposal to flatten the spectrum of projected energies $E^{I}$ (51) Dobaczewski (2009); Wang et al. (2014); Gao, Dobaczewski and Toivanen (2015),

[TABLE]

where

[TABLE]

are average values of the so-called Lipkin operator ${\hat{K}}$, evaluated between the projected states, and the constant on the right-hand side of equation (89) does not depend on labels $I$ of the projected states. Equivalently, as it is evident from equation (51), one can flatten the reduced energy kernels, i.e.,

[TABLE]

where the constant on the right-hand side, which is the same constant as in equation (89), does not depend on the group element $g$. Indeed, since the analogue of equation (51) relates $E^{I}-K^{I}$ to ${\hat{H}}-{\hat{K}}$, we can plug there the numerator of equation (91) expressed though its denominator, which immediately gives equation (89). Otherwise, by plugging equation (89) into the same analogue of equation (51), we see that the function on the group defined by$\langle\Phi|({\hat{H}}-{\hat{K}})\hat{R}(g)|\Phi\rangle-\langle\Phi|\hat{R}(g)|\Phi\rangle\times const.$ is orthogonal to all functions $D^{I}(g)$, and thus must be equal to zero, proving equation (91).

How to find the Lipkin operator ${\hat{K}}$ that does the job as desired? The strategy is obvious, namely, since the quantum numbers $I$ are related to the symmetry generators, one can build ${\hat{K}}$ as a function of the symmetry generators. For example:

• •

For the particle-number symmetry, projected energies depend on the number of particles, so ${\hat{K}}$ may depend on the particle-number operators, ${\hat{N}}$ and ${\hat{Z}}$, for neutrons and protons:

[TABLE]

where $N_{0}$ and $Z_{0}$ are numbers of protons and neutrons of the state we want to describe.

• •

For the translational symmetry, projected energies depend on total momenta, so ${\hat{K}}$ may depend on the components of the total momentum, ${\hat{P}}_{x}$, ${\hat{P}}_{y}$, and ${\hat{P}}_{z}$:

[TABLE]

• •

For the rotational symmetry, projected energies depend on total angular momentum, so ${\hat{K}}$ may depend on the components of the total angular momentum, ${\hat{J}}_{x}$, ${\hat{J}}_{y}$, and ${\hat{J}}_{z}$:

[TABLE]

• •

For the rotational symmetry and axial nucleus oriented along the $z$ axis, projected energies depend on total angular momenta, so ${\hat{K}}$ may depend on the total angular momentum $\bm{{\hat{J}}}^{2}$:

[TABLE]

where $I$ is the total angular momentum of the state we want to describe.

In equations (92)–(95) and below, constants $k$ represent adjustable parameters (the Lipkin parameters), which should best approximate $E^{I}$ by $K^{I}$, equation (89).

The main idea behind building the Lipkin operators is to have the best possible description of spectra $E^{I}$ of projected energies in terms of averages of group generators $K^{I}$. We stress that we do not deal here with real spectra of the system, but with energies $E^{I}$ of symmetry-conserving components $\hat{P}^{I}|\Phi\rangle$ derived from the symmetry-breaking state $|\Phi\rangle$. In fact, among all the projected states, we are interested in only one of them, that with energy $E^{I_{0}}$ and wave function $\hat{P}^{I_{0}}|\Phi\rangle$. Of course, if the flattening, equations (89) or (91), is perfect – this does not matter; however, if it is not perfect, we better do the best possible job for the one state $I_{0}$ that we want to describe. Then, the Lipkin operator constructed so that

[TABLE]

gives us obviously

[TABLE]

In fact, for the examples of the Lipkin operators presented in equations (92) and (95), condition (96) is fulfilled, whereas those in equations (93) or (94) apply to states at rest ($\bm{P}=0$) or non-rotating ($\bm{J}=0$), respectively.

At this point, by evaluating equation (91) at $g=0$, we obtain the “magic” Lipkin formula:

[TABLE]

namely, the projected energy $E^{I_{0}}$ can be obtained as an average value of ${\hat{H}}-{\hat{K}}$ calculated for the symmetry-breaking state $|\Phi\rangle$ without performing any projection at all. Of course, we can benefit from the magic formula only if we can find appropriate Lipkin operators that correctly flatten the spectrum, and the precision of it is dictated by the precision of the flattening.

Therefore, the main thrust of the method now lies in finding the numerical coefficients in equations (92)–(95) that define the Lipkin operators in terms of the symmetry generators. Before going into details of specific applications, let us introduce a generic form of the Lipkin operator as a linear combination of different terms:

[TABLE]

Following the original idea of Peierls and Yoccoz (1957), we now evaluate equation (91) at $M+1$ group elements, $g_{i}$, for $i=0,\ldots,{M}$, $g_{0}=0$, which leads to a set of linear equations that determine the Lipkin parameters $k_{m}$,

[TABLE]

where we extended the list of symmetry generators by defining,

[TABLE]

Then, coefficients in (100) are defined by the following reduced kernels:

[TABLE]

and the Lipkin parameters can be obtained by inverting matrix $A$:

[TABLE]

In doing so, we can always adjust the choice of group elements $g_{i}$ so as to obtain a non-singular matrix $A$. It is noted here that one can simply ignore the value of $k_{0}$ given by equation (104). Indeed, since it is by definition equal to $E^{I_{0}}$, equation (101), one can always evaluate it from the magic Lipkin formula (98).

Based on the Peierls-Yoccoz prescription to determine the Lipkin parameters, one has to calculate a few norm overlaps and energy transition matrix elements – the same ones that are required for the execution of the full projection, equation (51). Their number is, however, significantly smaller than that required for a full projection, and thus the Lipkin method is computationally much less intensive.

The biggest advantage of the Lipkin method manifests itself when we attempt to obtain the variation-after-projection solution. Indeed, an exact implementation leads then to difficult programming and lengthy calculations, cf. section 3.3. On the other hand, variation of the projected energy obtained from the magic Lipkin formula (98) is as easy as a direct variation of the energy of the symmetry-breaking state. Clearly, the Lipkin method gives only the projected energy, whereas, if average values of other observables are to be calculated, the full projection has to be anyhow performed. Then, the Lipkin method allows for obtaining variation-after-projection results at the expense of a single projection-after-variation calculation, which still constitutes a substantial gain in computing time and efficiency.

One should stress one important aspect of the Lipkin method, namely, when varying state $|\Phi\rangle$ that appears in the magic formula (98) to obtain the variation-after-projection result, one should treat the Lipkin parameters $k_{m}$ as constants that do not undergo variation. Indeed, even if their values parametrically depend on $|\Phi\rangle$ through equations (102)–(104), their role is to provide the best flattening of the final spectrum, so during the variation they should be kept fixed and equal to the final variation-after-projection values. For iterative solutions of self-consistent equations, the algorithm to achieve such a goal can easily be implemented. Indeed, it is enough to treat Lipkin parameters at every iteration as constants, and at the end of every iteration recalculate them using equations (102)–(104). Once the convergence is reached, such an algorithm yields the desired result. In the following sections, for specific cases we discuss solutions of equations (102)–(104).

4.2 The Lipkin-Nogami method

In a series of papers, Nogami and collaborators Nogami (1964); Nogami and Zucker (1964); Nogami (1965); Goodfellow and Nogami (1966); Pradhan et al. (1973) developed a variant of the Lipkin method that replaces the calculation of norm overlaps and Hamiltonian transition matrix elements by a calculation of several average values. This replacement can be derived by first rewriting equations (91), (99), and (101) as

[TABLE]

Since for any $g$, and in particular for $g$ near 0, the group operators $\hat{R}(g)$ are equal to the exponents of linear combinations of symmetry generators ${\hat{K}}_{i}$, it follows from equation (105) that for any $i$ we have,

[TABLE]

Now the Lipkin parameters can be evaluated in an analogous way to how it was done in the Lipkin approach (104):

[TABLE]

but for

[TABLE]

Had the Lipkin method been exact, the Lipkin-Nogami expressions would also be same, and would lead to exact results. Otherwise, the Lipkin and Lipkin-Nogami methods may give different results, and it is a priori difficult to say which one is superior. Nevertheless, if the calculation of transition matrix elements, and not only of the average values, is available, the Peierls-Yoccoz method is certainly easier to implement. Indeed, in case the Hamiltonian is a 2-body operator and the Lipkin operator is an $n$-body operator, the Peierls-Yoccoz method requires calculating transition matrix elements of these 2-body and $n$-body operators only, whereas the Lipkin-Nogami calls for calculating averages of $n$+2-body and $2n$-body operators. In spite of that, at second order ($n=2$) and for the particle-number projection, the Lipkin-Nogami method was applied quite widely, see section 4.4.2.

4.3 The Kamlah method

The principal idea of the Kamlah expansion Kamlah (1968) is that the Hamiltonian transition matrix elements can by efficiently expanded into a series of derivative operators ${\cal K}_{m}$ acting on the norm overlaps:

[TABLE]

This expansion is supposed to work best in the limit of strong symmetry breaking, for example, at large deformations. In this limit, the norm overlaps and Hamiltonian transition matrix elements are both strongly peaked near $g=0$, and therefore, the expansion of the former in a series of derivatives of the latter may have a chance to converge rapidly.

Since for every continuous group, polynomials of symmetry generators ${\hat{K}}_{m}$ can always be represented by derivatives ${\cal K}_{m}$ with respect to the group parameters,

[TABLE]

Kamlah expansion (110) is strictly equivalent to the Lipkin flattening condition (91) applied for the Lipkin operator of equations (99) and (101). Since, in addition, Kamlah proposes to determine coefficients $k_{m}$ by evaluating derivatives at $g=0$, his method gives equations for $k_{m}$ that are strictly equivalent to the Lipkin-Nogami method (106).

There are, nevertheless, two important differences. First, the Kamlah proposal involves variation of the projected energy (98) “as it is”, i.e., a variation over symmetry-breaking states $|\Phi\rangle$ should also involve variation of $k_{m}$, see discussion in section 4.4.2. Second, the Kamlah expansion may contain terms that are not invariants of the symmetry group, and therefore, they are not really within the realm of the Lipkin method. This latter property mostly relates to the so-called pushing and cranking terms discussed in detail in section 4.4.1.

4.4 Applicability and Applications

The Lipkin, Kamlah, and Lipkin-Nogami methods, discussed in the previous sections, all rely on polynomial expansions of collective spectra or reduced kernels. This principal assumption creates two main limitations of these approaches. First, obviously, the expansions have to be carried out up to a sufficiently high order, see section 4.4.3. And second, and most importantly, these methods cannot really be applied to spectra that have a non-analytical dependence on the quantum numbers. Unfortunately, the latter situation occurs in two physically meaningful cases, namely, when particle-numbers are restored in (semi)magic nuclei Dobaczewski and Nazarewicz (1993); Wang et al. (2014) and when the angular-momentum is restored in weakly-deformed systems Gao, Dobaczewski and Toivanen (2015).

4.4.1 First-order terms

*The Fermi energy:

*The simplest application of the Lipkin, Kamlah, or Lipkin-Nogami methods concerns the approximated restoration of the particle-number symmetry up to the first order in the particle number (92), i.e.,

[TABLE]

where we can treat numbers of protons and neutrons separately. Then, the Lipkin-Nogami equations (107) and (109) give

[TABLE]

where brackets $\langle\rangle$ denote average values calculated for the particle-number-symmetry-breaking state $|\Phi\rangle$. The variation-after-projection equation, which is derived from (98), now reads

[TABLE]

for

[TABLE]

According to Lipkin’s methodology, variation over $|\Phi\rangle$ has to be carried out at constant Lipkin coefficient $k_{1}$, and according to the Kamlah’s methodology, expression for $k_{1}$ (114) should be inserted into (115) and then varied.

In this sense, at first order, the Lipkin-Nogami and Kamlah prescriptions lead to the same result. Moreover, the Lipkin coefficient $k_{1}$ can now be reinterpreted as a Lagrange multiplier $\lambda_{1}$, that is, as a Fermi energy, which has to be adjusted so as to obtain the correct average particle number. Then, the Lipkin-Nogami expressions (113) and (114) simplify tremendously, and give

[TABLE]

The expression for $k_{0}$ is thus compatible with (116) and that for $k_{1}$ stems from (115), provided $\delta_{\Phi}|\Phi\rangle=({\hat{N}}-N_{0})|\Phi\rangle$ is an allowed variation.

In a similar way, we can evaluate the Lipkin expressions (104), which gives

[TABLE]

where

[TABLE]

are reduced kernels of the shifted particle-number operator, evaluated at gauge angle $\phi_{i}$, and $h(\phi_{i})=h_{i}$ are the analogous reduced kernels of the Hamiltonian (103). Again it is beneficial to carry out variation (115) with the average particle number kept correct, $\langle{\hat{N}}\rangle=N_{0}$, which gives $n_{1}(0)=0$. In this case, expression (119) reduces again to $k_{0}=\langle\hat{H}\rangle$ and expression (120) stems from (115), provided the finite-difference derivatives are allowed as variations $\delta_{\Phi}|\Phi\rangle$.

We conclude, that the Lipkin, Lipkin-Nogami, and Kamlah symmetry restoration at first-order are completely equivalent to using Lagrange multipliers for adjusting average values of symmetry generators.

*The pushing model and Thouless-Valatin mass:

*For the restoration of translational symmetry, at first order the Lipkin operator reads,

[TABLE]

where (93) is generalized to the case of a nucleus moving with the average total momentum $\bm{P}_{0}$. We note here that the lowest-order invariant of the translational group is equal to ${\hat{\bm{P}}}^{2}$, and therefore, the flattening of the spectrum requires using the second-order Lipkin operator, see section 4.4.2. Therefore, the first-order model (122) rather pertains to the Kamlah approach.

Since components of the momentum operator ${\hat{\bm{P}}}$ commute, we can treat them independently. Then, following the derivations presented for the particle number, we conclude that the restoration of translational symmetry is, at first order, equivalent to performing minimization of the total energy, constrained to the given momentum $\bm{P}_{0}$,

[TABLE]

where the vector Lipkin coefficient $\bm{k}_{1}$ acquires interpretation of the Lagrange multiplier $\bm{v}$, i.e., of the velocity of the system.

For the translational symmetry, variation of the total energy $E^{\bm{P}_{0}}$ constrained to the momentum $\bm{P}_{0}$ is particularly simple. Indeed, suppose we found the state $|\Phi_{\bm{P}_{0}=\bm{0}}\rangle$, which is at rest, $\langle\Phi_{\bm{P}_{0}=\bm{0}}|\hat{\bm{P}}|\Phi_{\bm{P}_{0}=\bm{0}}\rangle=\bm{0}$, and fulfills the variation-after-projection equation $\delta{}E^{\bm{P}_{0}=\bm{0}}=0$ for $\bm{v}=\bm{0}$. Then, the Galilean invariance,

[TABLE]

where $\hat{\bm{R}}$ is the center-of-mass coordinate and $Am$ is the total mass of the system, allows us to boost state $|\Phi_{\bm{P}_{0}=\bm{0}}\rangle$ to momentum $\bm{P}_{0}$,

[TABLE]

so that

[TABLE]

We see that the restoration of translational symmetry at first order, that is, the pushing model, correctly reproduces all classical-motion relations. In particular, from the analog of the Lipkin-Nogami expression (118) we obtain the velocity vector as,

[TABLE]

which gives the mass

[TABLE]

that, in the translational case, correctly reproduces the true mass of the system, ${\cal{M}}=Am$ Thouless and Valatin (1962).

*The cranking model and Thouless-Valatin moment of inertia:

*Restoration of the rotational symmetry at first order leads to the so-called cranking model, which was introduced originally in a semiclassical time-dependent picture of a system rotating with constant angular velocity $\bm{\omega}$ around a fixed axis parallel to the angular momentum $\bm{J}$ Inglis (1954, 1956). This model was very successfully used in nuclear physics to describe a multitude of physical phenomena related to collective rotation Bohr and Mottelson (1998); Ring and Schuck (1980). In this case, the Lipkin operator reads,

[TABLE]

Since components of the angular-momentum operator ${\hat{\bm{J}}}$ do not commute, we cannot treat them independently. equation (131) has thus to be understood as corresponding to a nucleus having a fixed projection $J_{0}=|\bm{J}_{0}|$ of the angular momentum on a quantization axis oriented along the Lipkin coefficient vector, $\bm{k}_{1}=\bm{\omega}$. In systems with approximate axial symmetry along the $z$-axis, for instance, for the ground-state bands in well deformed even-even nuclei, the rotational axis is perpendicular to the symmetry axis and usually chosen along the $x$-axis. In this case, we have $J_{0}=\langle J_{x}\rangle=\sqrt{I(I+1)}$ with integer values of $I$ (for odd systems, see Ring et al. (1974)). Different directions of $\bm{k}_{1}$ then mean a freedom of choosing an arbitrary direction of the quantization axis. This defines the so-called tilted-axis cranking model Kerman and Onishi (1981); Frauendorf (1993, 2001); Shi et al. (2013), where the vector of the average angular momentum $\langle{\hat{\bm{J}}}\rangle$ is arbitrarily oriented with respect to the principal axes of the mass distribution of the rotational-symmetry-breaking state.

Following the derivations presented for the momentum operator, we conclude that the restoration of rotational symmetry is, at first order, equivalent to performing minimization of the total energy, constrained to the given projection of the angular momentum $J_{0}$,

[TABLE]

where the vector Lipkin coefficient $\bm{k}_{1}$ acquires interpretation of the Lagrange multiplier $\bm{\omega}$, that is, of the angular velocity of the system.

The principal difference between translational and rotational symmetry is the fact that for rotations there is no analogue of the Galilean invariance (124), and one cannot simply boost a non-rotating state to higher rotational frequencies without changing its structure. Indeed, with increasing rotational frequency, the quantum analogues of the classical Coriolis and centrifugal forces set in, and modify the state. Therefore, a constrained minimization of the total energy has now to be explicitly performed.

Recall that the average momentum is exactly proportional to the translational velocity, equation (127), with a constant proportionality factor (mass). Although the average angular momentum $\langle\hat{\bm{J}}\rangle$ has to be parallel to the angular frequency $\bm{\omega}$ Kerman and Onishi (1981), the proportionality constant (moment of inertia) can vary along the rotational band. Therefore, we define two important local, that is, frequency-dependent characteristics of the band, which are called the first ${\cal{J}}^{(1)}$ (kinematic) and the second ${\cal{J}}^{(2)}$ (dynamic) moments of inertia,

[TABLE]

respectively.

In parallel with equation (130), the Lipkin expression allows us to determine the kinematic moment of inertia,

[TABLE]

The dynamic moment of inertia is identical to this value only at $\omega=0$. For all other values of $\omega$, it corresponds to the Thouless-Valatin moment of inertia ${\cal{J}}^{(2)}={\cal{J}_{\mbox{\scriptsize{TV}}}}$ that can be derived in linear response theory Thouless and Valatin (1962).

*Isocranking:

*To perform the approximate restoration of the isospin symmetry, Frauendorf and Sheikh Frauendorf and Sheikh (1999) introduced the concept of cranking in isospin with a focus on the isovector proton-neutron pair field without Coulomb interaction. In this isospin-symmetry conserving framework, the direction of the isocranking axis in the isospin space is not relevant. For a specific choice of this direction, the isocranking about the isospin $z$-axis corresponds to the standard use of the neutron and proton Fermi energies that are not equal to one another. Satuła and Wyss Satuła and Wyss (2001b, a); Głowacz et al. (2004) used the same terminology of isocranking in their isocranking model with Coulomb interaction included, where dynamical consequences of the isocranking about the axis tilted with respect to the isospin $z$-axis became relevant, and where they could investigate both isovector and isoscalar pair fields.

In the isospin context, the Lipkin operator reads as

[TABLE]

where arrows denote vectors in the isospace (isovectors) and symbol "$\circ$" denotes their scalar product. Since components of the isospin operator ${\hat{\vec{T}}}$ do not commute, we cannot treat them independently. Thus equation (135) has to be understood as corresponding to a nucleus having a fixed projection $T_{0}=|\vec{T}_{0}|$ of the isospin on the isoquantization axis oriented along the isovector Lipkin coefficient $\vec{k}_{1}$. Moreover, in even (odd) systems, projections $T_{0}$ can only equal to integer (half-integer) numbers. Different directions of $\vec{k}_{1}$ then mean a freedom of choosing an arbitrary direction of the isoquantization axis. The Lipkin coefficient $\vec{k}_{1}$ is interpreted as the isovector Fermi energy $\vec{\lambda}$ Sato et al. (2013); Sheikh et al. (2014), which fixes the average values of components of the isospin $\langle{\hat{\vec{T}}}\rangle$.

We note that the standard definition of the isospin implies that its $z$ component is equal to half of the neutron excess, $T_{z}={\textstyle\frac{1}{2}}(N-Z)$. Therefore, the $z$ component of the isovector Fermi energy $\lambda_{z}$ along with the standard isoscalar Fermi energy $\lambda$ simply fix the neutron $N$ and proton $Z$ numbers. When the isocranking axis is tilted away from the $z$ direction, one must use the formalism where proton and neutron components of single-particle states are mixed Sato et al. (2013); Sheikh et al. (2014). Such a situation occurs when the isospin-symmetry-breaking terms are added to the nuclear Hamiltonian Bączyk et al. (2018). Figure 5 shows energies of states in 48Cr isocranked to $\langle{}T_{x}\rangle=0$, 2, 4, 6, and 8, while keeping $\langle{}T_{z}\rangle=0$ Sato et al. (2013). We see that one obtains a perfectly rigid isorotational band (a sequence of states in 48Cr with increasing isospin $T$). The obtained Thouless-Valatin moment of inertia then corresponds to the symmetry energy coefficient $a_{I}$ in the symmetry energy, $E_{I}(N,Z)={\textstyle\frac{1}{2}}a_{I}(N-Z)/A$, i.e., to $a_{I}=2A\times 1.39(2)/4=33.4(5)$ MeV.

4.4.2 Second-order terms

Although at first order the Lipkin or Kamlah approach provides correct understanding of the collective effects, including the proper determination of the collective mass, it does not really fulfill Lipkin’s requirement of flattening the spectrum of projected states. Indeed, already from the example of the translational motion, we see that the main component of this dependence may rather be quadratic then linear. In this section, we thus examine the Lipkin operators expanded up to second-order terms in the symmetry generators. This allows us to model the spectra in terms of the quadratic Casimir operators of the corresponding symmetry groups.

*Particle-number restoration:

*The main focus of the Lipkin-Nogami method Nogami (1964); Nogami and Zucker (1964); Nogami (1965); Goodfellow and Nogami (1966); Pradhan et al. (1973) was up to now on the approximate restoration of the particle-number symmetry, whereby the Lipkin operator is postulated as,

[TABLE]

Assuming again that $k_{1}$ is always adjusted so as to obtain the correct particle number, $\langle{\hat{N}}-N_{0}\rangle=0$, the Lipkin (104) Wang et al. (2014) and Lipkin-Nogami (107) Valor et al. (1996) methods give, respectively,

[TABLE]

and

[TABLE]

Here $\Delta\hat{N}\equiv\hat{N}-N_{0}$ is the shifted particle-number operator, $h(\phi_{2})$ and $n_{i}(\phi_{2})$ are the Hamiltonian (103) and shifted particle-number (121) kernels, and we used the traditional notation of $\lambda_{2}\equiv{}k_{2}$ for the second-order Lipkin-Nogami coefficient. For an HFB vacuum $|0\rangle$, an alternative and equivalent expression for $\lambda_{2}$ was derived in Sandhu and Rustgi (1978) as

[TABLE]

where $|4\rangle\equiv\alpha^{+}_{\mu}\alpha^{+}_{\nu}\alpha^{+}_{\mu^{\prime}}\alpha^{+}_{\nu^{\prime}}|0\rangle$ stands for all four-quasiparticle states. After evaluating all required matrix elements, one obtains Flocard and Onishi (1997); Stoitsov et al. (2003)

[TABLE]

where the potentials

[TABLE]

can be calculated in full analogy to $\Gamma$ and $\Delta$ by replacing the $\rho$ and $\kappa$ in terms of which they are defined by $\rho(1-\rho)$ and $\rho\kappa$, respectively. In the case of the seniority pairing interaction with strength $G$, equation (140) simplifies to Pradhan et al. (1973):

[TABLE]

where $k>0$ denotes the summation over one state of each canonical (or Kramers-degenerate) pair of single-particle states.

Evaluation of the Lipkin coefficient $k_{2}$, as in equation (137), is fairly simple, see section 4.4.3, but it was implemented only in Wang et al. (2014). A rigorous evaluation of the Lipkin-Nogami coefficient $\lambda_{2}$, equation (138) or (140), is for realistic Hamiltonians rather cumbersome, so it was rarely implemented in full, see, e.g., Valor et al. (1996). A practical workaround, which was used quite often, see, e.g., Stoitsov et al. (2003, 2007); Kortelainen et al. (2010) was to use the seniority-pairing expression equation (143) with the effective strength $G\equiv{}G_{\mbox{\scriptsize{eff}}}=-\frac{\bar{\Delta}^{2}}{E_{\mbox{\scriptsize{pair}}}}$, determined from the pairing energy $E_{\mbox{\scriptsize{pair}}}$ and the average pairing gap $\bar{\Delta}$.

Figure 6 shows comparison of the Lipkin-Nogami and Kamlah approaches within a simple solvable two-level model Zheng et al. (1992); Dobaczewski and Nazarewicz (1993). In the strong-pairing regime, both approaches give excellent reproduction of the exact results, however, for the half-filled shell [figure 6(b)], at or below the critical pairing strength both fail. This is so because the kink in the dependence of the exact energies on the particle number, which is a characteristic feature of a shell gap, cannot be reproduced by the quadratic Dobaczewski and Nazarewicz (1993) or higher-order Wang et al. (2014) form of the Lipkin operator. Away from the shell gap [figure 6(a)], the Lipkin-Nogami approach works well for all pairing strengths, but the Kamlah approach fails below the critical pairing strength. This latter feature can be attributed to the fact that by exploiting errors of the approximation, the exact minimization over the Lipkin coefficient $\lambda_{2}$ brings the approximate projected energy below the exact result.

*Peierls-Yoccoz mass:

*For the translational symmetry, at second order, the general form of the Lipkin operator reads

[TABLE]

As already discussed in section 4.4.1, the first-order terms define the pushing model, wherein the nucleus moves in space with average momentum $\bm{P}_{0}$ and velocity $\bm{v}\equiv\bm{k}_{1}$. This motion leads to an increase of energy that is quadratic in the momentum, and the corresponding proportionality coefficient is called the Thouless-Valatin mass (130), which is correctly equal to the translational mass $mA$.

The role of the Lipkin coefficients $k_{2n}$ is different – they are meant to flatten the spectrum of energies projected on momentum eigenstates in the direction of $n=x$, $y$, or $z$. Their values thus characterize the momentum distributions within the translational-symmetry-breaking state, and have nothing to do with the physical motion of the system. This is particularly evident when we consider the state at rest, $|\Phi_{\bm{P}_{0}=\bm{0}}\rangle$, which can be obtained by simply conserving the time-reversal symmetry. This state does not move, so the Lipkin coefficients $k_{2n}$ cannot describe inertia, which is the reaction of the system under boost.

Coefficients $k_{2n}$ calculated for the $n=x$, $y$, or $z$ directions can differ from one another Gao, Dobaczewski and Toivanen (2015). Indeed, along the longer or shorter principal axis of the mass distribution, the momentum distribution is narrower or wider, respectively, and the corresponding projected energy components can thus differ from one another.

Nevertheless, historically, quantities ${\cal M}_{\mbox{\scriptsize{PYn}}}={\textstyle\frac{1}{2}}k_{2n}^{-1}$, corresponding to translational Lipkin coefficients $k_{2n}$, are called Peierls-Yoccoz Peierls and Yoccoz (1957) or Yoccoz Ring and Schuck (1980) masses. Assuming that those corresponding to the $n=x$, $y$, or $z$ directions are independent from one another, they can be calculated in the Lipkin or Lipkin-Nogami approach, respectively, as:

[TABLE]

or

[TABLE]

cf. equations (137) or (138). Here $\Delta\hat{P}_{n}={\hat{P}}_{n}-P_{n0}$ are shifted momentum operators for $n=x$, $y$, or $z$ and $p_{in}(\phi_{2n})$ are their reduced kernels calculated at distances $\phi_{2n}$.

For conserved time reversal, average values of all odd powers of momentum are equal to zero, and thus equation (4.4.2) reduces to

[TABLE]

In addition, if we consider a spherical nucleus, where the Lipkin coefficients corresponding to three directions $n=x$, $y$, or $z$ are equal, and the Lipkin operator takes the form ${\hat{K}}=k_{2}\hat{\bm{P}}^{2}$, we then have

[TABLE]

We see that equality of $k_{2n}$ and $k_{2}$ requires independence of the three directions, $\langle\hat{P}_{n}^{2}\hat{P}_{m}^{2}\rangle=\langle\hat{P}_{n}^{2}\rangle\langle\hat{P}_{m}^{2}\rangle$ for $n\neq{}m$, which was assumed when deriving equations (4.4.2) and (147). Further, within the Gaussian Overlap Approximation Ring and Schuck (1980), we have $\langle\hat{\bm{P}}^{4}\rangle=3\langle\hat{\bm{P}}^{2}\rangle^{2}$, and equation (148) simplifies to

[TABLE]

Figure 7 shows comparison of the Peierls-Yoccoz and Thouless-Valatin (exact) masses calculated in doubly magic nuclei. We see that the former are never equal to the latter, because they represent different quantities. Indeed, the Peierls-Yoccoz masses characterize the curvatures of energies projected from the symmetry-breaking states at rest, whereas the Thouless-Valatin masses characterize the increase of the energy when the symmetry-breaking states are boosted to non-zero momenta. In addition, the figure shows the Peierls-Yoccoz masses evaluated for energies minimized before and after including the Lipkin operator. One clearly sees that the self-consistent inclusion of the Lipkin correction does modify the curvatures of projected energies. In practice, however, differences between the Peierls-Yoccoz and Thouless-Valatin translational masses do not exceed 10% and vary smoothly with nuclear masses. More discussion can be found in Bender et al. (2000); Dobaczewski (2009).

Finally, let us mention here a popular way of implementing the so-called center-of-mass correction in self-consistent calculations Bender et al. (2000); Dobaczewski (2009), which amounts to combining the magic formula (98) and the second-order Lipkin operator (144) with the Thouless-Valatin mass into an approximate expression for the projected energy,

[TABLE]

The Peierls-Yoccoz mass, which should appear in this formula is not evaluated but approximated here by the true inertial mass of the system. Further approximations of the average value $\langle\Phi|\hat{\bm{P}}^{2}|\Phi\rangle$ are also frequently used Bender et al. (2000); Dobaczewski (2009).

*Peierls-Yoccoz moment of inertia:

*For the rotational symmetry, at second order, the Lipkin operator reads

[TABLE]

As already discussed in section 4.4.1, the first-order terms define the cranking model, wherein the nucleus rotates in space with average angular momentum $\bm{J}_{0}$ and frequency $\bm{\omega}\equiv\bm{k}_{1}$. This motion leads to an increase of energy, which is (in deformed nuclei) approximately quadratic in the angular momentum, and the corresponding proportionality coefficient is called the Thouless-Valatin moment of inertia (134). On the other hand, the Lipkin coefficients $k_{2n}$ are meant to flatten the spectrum of energies projected on angular-momentum eigenstates. However, since the three components of the angular momentum do not commute, the three first-order and three second-order terms in the Lipkin operator (151) are not independent from one another. Moreover, in axial nuclei, the rotational symmetry is not fully broken, namely, the total angular momenta are mixed, but the projections of the angular momentum on the symmetry axis continue to be a good quantum numbers. In this case, the second-order Lipkin operator in the form of equation (95) is more appropriate.

Expressions for Peierls-Yoccoz moments of inertia can be obtained in full analogy to those for Peierls-Yoccoz masses, equations (4.4.2)–(149). In particular, for a one-dimensional rotation of an axial nucleus about the $x$ axis perpendicular to the symmetry axis, the analogue of equation (147) reads

[TABLE]

In analogy with approximate expressions used for translations (150), also the so-called rotational-motion corrections Frauendorf et al. (1969); Bender et al. (2003) are sometimes implemented along with various approximations for the Peierls-Yoccoz moments of inertia.

*Isospin:

*Almost immediately after Nogami’s work Nogami (1964), an analogous method was suggested in the restoration of the isospin symmetry Kissener and Münchow (1966, 1967); Ghosh et al. (1975), whereupon the Lipkin operator,

[TABLE]

combines the isocranking term (135) with the second-order correction. Such approach, however, was later not too often employed, because the exact isospin restoration was early implemented Caurier et al. (1980a, b); Caurier and Poves (1982) and is now efficiently used Satuła et al. (2009, 2010, 2016).

4.4.3 Higher-order terms

In the Lipkin approach or Kamlah expansion, higher-order terms were studied for the particle-number restoration. Rodríguez et al. (2005) introduced the so-called reduced variation-after-projection method, which aimed to improve second-order Lipkin-Nogami or Kamlah results. The Lipkin method up to sixth order was implemented in Wang et al. (2014). This work showed that away from semi-magic nuclei, the Lipkin method converges already at 4th order, with the 6th order corrections bringing almost no change. Near the semi-magic nuclei, however, the non-analytical dependence of energy in function of the particle number did not allow for obtaining well converged results.

5 Projection methods in simple nuclear models and ab initio calculations

The basic objective of the projection methods is to include the many-body correlations beyond the mean-field level, and in order to treat them accurately, the best is to perform them in the $\bm{r}$-space Baye and Heenen (1984). However, projection in the $\bm{r}$-space leads to solution of non-local potential problem, even with zero-range effective interaction. The non-local potential problem is quite prohibitive to solve even with the modern computing facilities. This has also to do with the fact that even the bare three-dimensional HFB problem is very difficult to solve in $\bm{r}$-space. Most of the modern and efficient 3D methods in $\bm{r}$-space are restricted the HF- or HF+BCS calculations Ryssens et al. (2019), whereas full 3D-HFB calculations are carried out in oscillator space  Egido et al. (1995). Two-dimensional HFB equations were solved in $\bm{r}$-space using the spline functions Kegley Jr. et al. (1996); Pei et al. (2008) and there were also attempts to solve the three-dimensional HFB problem using the multi-wavelet method Pei et al. (2014). As evident from these studies, the solution of the full HFB equations in the three-dimensional $\bm{r}$-space requires, at present, a too high computational effort to be applied in realistic applications. For the above reason, most of the application of projection methods were carried out by expanding the single-particle wave functions in a finite basis. This leads to relatively simple matrix operations. The basis states, normally chosen, are solutions of harmonic oscillator or Woods-Saxon potentials. A very efficient method relies on solving the HFB equations on the basis of eigenstates of the HF Hamiltonian Gall et al. (1994); Bender et al. (2003) that can be obtained in the $\bm{r}$-space.

For axial and triaxial mean-field calculations, the oscillator basis has the advantage that these wave functions are separable in the coordinates. This requires considerably less memory for the storage of the basis states and, in addition, such a basis can be deformed. Through a careful choice of the basis deformation parameters, this allows us to reproduce, in a relatively small basis of roughly 20 major oscillator shells, the exact mean-field results of the system with sufficient accuracy, i.e., with roughly 500 keV in the total nuclear binding energies even for heavy nuclei Dobaczewski et al. (2002). For the nucleus 208Pb with an experimental binding energy of 1636.446 MeV this corresponds to an accuracy of 0.3 permille. We also have to keep in mind, that the binding energy, being a very sensitive difference between two large quantities, the repulsive kinetic energy and the attractive potential energy, requires often much more accuracy than other quantities that are mostly determined by the properties of the valence shells, as for instance deformation energies etc.

However, for the application of projection methods, the use of a deformed basis leads to severe problems (see discussion in section 6.1.1), because the basis violates the corresponding symmetry from the beginning. By this reason, most of the applications of angular-momentum projection were carried out in a spherical basis.

In most of the realistic nuclear models, for instance, in the DFT approach, the projection methods were applied after variation and include correlations only partially beyond the mean-field level. Projection calculations before variation become exceedingly difficult, in particular, for the case of three-dimensional projection. Such studies were performed with simpler model Hamiltonians using a few major oscillator shells around the Fermi surface McCullen et al. (1964); Lawson (1961); Baranger and Kumar (1968).

The justification for using such a simpler approach, which considers only a restricted set of active orbits around the Fermi surface (what is often referred to in the literature as “the valence space"), is that we are mostly interested in properties associated to low energy excited states. Also, the pairing part of the effective interaction is dominated by single-particle states around the Fermi surface. Even the correlations leading to deformation are dominated by contributions from a valence space of a few major oscillator shells. It is known from several studies that states that are far above or below the Fermi surface do not contribute to the correlations beyond the mean-field level. This is clearly evident in the case of particle-number projection, where the states far from the Fermi surface have occupations close to either one or zero and, therefore, do not contribute to pairing correlations.

Such valence spaces are often too large for what is called “full configuration-interaction" calculations777The term refers to those calculations that consider as elements of the Fock state all mean field states built by taking into account all possible distribution of particles among the active orbits. The typical example being the interacting shell model Caurier et al. (2005)., but still small enough even for sophisticated projection techniques. A famous case is the Baranger-Kumar valence space Baranger and Kumar (1968) containing the shells $N=3,4$ for protons and $N=4,5$ for neutrons and used for the description of the well deformed Rare Earth region of the nuclear chart. Its configuration space contains 72 neutron and 50 proton levels. This leads for mean-field calculations to matrices with dimension 72 and 50, which are easy to handle. On the other side in a typical configuration-interaction calculation we have in the middle of the shells roughly $5\times 10^{34}$ configurations.

Of course, the model interaction has to be adjusted carefully to the underlying configuration space. If the interaction in such a valence space is properly chosen, much of the physics, in particular, the interplay between collective degrees of freedom and single-particle degrees of freedom can be well described in such a space, only global properties such as total binding energies, radii, or multipole moments have to be treated in the full space.

5.1 The pairing-plus-quadrupole model

The applications are further simplified by the use of separable interactions. From the configuration-interaction calculations Dufour and Zuker (1996); Caurier et al. (2005), one knows that an effective interaction in such a restricted space can be represented as a sum of terms separable in the particle-hole and particle-particle channel, i.e.,

[TABLE]

where,

[TABLE]

that is, $\hat{Q}_{\lambda\mu}$ and $\hat{P}_{\lambda\mu}$ are multipole operators of multipolarity $\lambda$ in the particle-hole and particle-particle channel, respectively, with $\chi_{\lambda}$ and $G_{\lambda}$ being the corresponding coupling constants. The Hamiltonian of the model is

[TABLE]

where the doubly magic spherical core enters the calculations in terms of the single-particle energies $\varepsilon_{n}$ in the operator

[TABLE]

The parameter of such models are the single-particle energies $\varepsilon_{n}$ and the coupling constants $\chi_{\lambda}$ and $G_{\lambda}$.

In the configuration-interaction calculations with interactions derived from ab initio calculations Dufour and Zuker (1996), it is found while deriving separable representations of these effective forces that the quadrupole part with $\lambda=2$ plays the essential role in the particle-hole channel and the monopole part with $\lambda=0$ dominates the particle-particle channel. These two parts do not depend much on the underlying bare nucleon-nucleon force. Therefore, it is easy to understand that the pairing plus quadrupole model, introduced in the sixties Baranger and Kumar (1968), was pretty successful to describe the bulk of the important long-range correlations in nuclei in a very efficient way. Of course, it can be improved without great difficulties through the inclusion of additional separable terms, such as octupole ($\hat{Q}_{3}$) and hexadecapole ($\hat{Q}_{4}$) operators or the quadrupole pairing term $\hat{P}_{2}$. Some of the advantages of such models are:

1. (i)

one works in a spherical basis which preserves, apart from translational invariance, all the symmetries; 2. (ii)

the basis and the corresponding single-particle matrices are relatively small and, therefore, projection techniques can be applied easily; 3. (iii)

if the parameters are carefully adjusted, one obtains excellent results which can be compared with the experimental data and the basic physics can be understood in this way.

Of course, such model have restrictions as quantities influenced by physics outside of the valence shells cannot be described. These include, most importantly, the total binding energy of the nucleus and also the phenomena driven by so-called intruder states coming from the excluded higher shells or coming from the core at large deformations. The fission process and superdeformed Sheline et al. (1972) configurations belong to this category. There are also limitations for phenomena driven by special forces not considered in the model Hamiltonian, for instance, such parts of the tensor force Otsuka et al. (2020) which are not included in the phenomenologically adjusted single-particle energies.

In figure 8, we show self-consistent cranking results for the pairing gap for protons and neutrons in the nucleus 168Hf as a function of the average angular momentum, $I$ calculated with and without particle-number projection before the variation. Without projection, we observe in the neutron pairing collapse at spin $I=24\hbar$. The proton gap is also somewhat quenched after $I=10\hbar$. On the other hand, for the case of number projection the pairing correlations are reduced smoothly.

Angular-momentum projection in triaxial systems is much more complicated. So far there are very few results for variation after projection. In most of the cases the variation of the angular-momentum projected energy is restricted to a few external parameters as, for instance, deformation parameters. In this case, one needs for each deformation only to evaluate the projected energy in equation (51) after a constrained mean-field calculation. This is considerably simpler than the evaluation of the projected gradient in each step of the iteration.

In figure 9, we show as an example [taken from Hayashi et al. (1984)], the energy surface of the nucleus 188Os as a function of the Bohr quadrupole deformation parameters $\beta$ and $\gamma$. The unprojected energy does not depend on the triaxiality $\gamma$, and a mean-field calculation without constraint would be relatively unstable. However, projection on the angular momentum $I=0$ leads to a clear triaxially deformed minimum with $\gamma\approx 30^{\circ}$.

Several advanced nuclear models were developed based on the pairing plus quadrupole-quadrupole interaction and the projection theory. Here, we discuss the approach of triaxial projected shell model (TPSM) Sheikh and Hara (1999), which was extensively employed in recent years to investigate the high-spin band structures in triaxial nuclei. In this method, intrinsic basis states are obtained by solving the triaxial Nilsson potential with the expected deformation values for the system under investigation. Explicit three-dimensional angular-momentum projection method is then employed to project out the states with good angular-momentum. Apart from the projected vacuum state, multi-quasiparticle states are also projected to the laboratory frame of reference in this approach. In the most recent version of the TPSM approach, for even-even systems, quasiparticle basis states employed are:

[TABLE]

where $|\Phi\rangle$ is the vacuum state and $\alpha^{+}_{n}~{}(\alpha^{+}_{p})$ are neutron (proton) quasiparticle operators.

The projected basis of equation (158) is then used to diagonalize the pairing-plus-quadrupole Hamiltonian, equations (154)–(157), consisting only of monopole ($\lambda=0$) and quadrupole ($\lambda=2$) terms in the particle-hole and particle-particle chanel, respectively. The $QQ$-force strength $\chi_{2}$ is adjusted such that the physical quadrupole deformation $\beta$ is obtained as a result of the self-consistent mean-field HFB calculation Hara and Sun (1995).

As an illustrative example, we present the TPSM results recently obtained in Jehangir et al. (2017) for the high-spin band structures observed in 156Dy. For this system, rich band structures have been populated in Majola et al. (2015). The TPSM calculated energies for various band structures are compared with the experimental data in figure 10. It is quite evident from the figure that TPSM results reproduce the data quite accurately. At first glance it looks somewhat surprising that angular-momentum projection from a fixed mean-field employed in the TPSM approach can reproduce the data so nicely up to quite high-spin. However, it is to be reminded that in this approach, angular-momentum projection is performed not only from the zero-quasiparticle, but also from multi-quasiparticle states. Further, configuration mixing of these quasiparticle states is performed through diagonalization of the shell model Hamiltonian Hara and Sun (1995), which we consider takes into account minor perturbations in the mean-field potential with increasing spin.

We would like to add that considering only mixing from two-quasiparticle states without projection corresponds to the Tamm-Dancoff approximation (TDA) that is known to show a number of disadvantages as compared to the random phase approximation (RPA) Ring and Schuck (1980), which can be justified as the small amplitude limit of time-dependent-mean field theory. In order to take into account ground state correlations beyond the projected two-quasiparticle TDA, several versions of a symmetry conserving RPA were introduced in Federschmidt and Ring (1985); Sangfelt et al. (1987); Kyotoku et al. (1990); Dukelsky et al. (2019). Applications of these theories, however, were restricted, so far, to very light systems or to simple group-theoretical models.

5.2 Projection in small configuration spaces

Over the years many applications of various projection methods were carried out in relatively small configuration spaces and for light nuclei. This was done in the first place due to the fact that many of these techniques require a considerable numerical effort and that it was or is impossible to apply them in large spaces and for heavy nuclei. A second reason was also the fact that exact configuration-interaction calculations can be carried out in such small spaces, which allow a comparison of the projection methods with the exact solution in the corresponding space.

A variation after projection on particle number and angular momentum was carried out already in the eighties by the Tübingen group Schmid et al. (1984a). This method was called VAMPIR (variation after mean field projection in realistic model spaces). First it was restricted to axially symmetric intrinsic HFB wavefunctions in the valence space containing the 0s a 0d orbits (the sd-shell valence space) and a phenomenological interactions adjusted to this space were used. The results were compared with configuration-mixing calculations. At low spins the results are reasonable. However, in this approach the intrinsic wavefunctions obey time-reversal symmetry. This means that they are not flexible enough as to contain all the physics at high spin and the resulting moments of inertia are too small. The missing physics can be considered by mixing to the ground state band with projected two-quasiparticle configurations. This model Schmid et al. (1984c, b) was called MONSTER (model handling many number- and spin-projected two-quasiparticle excitations with realistic interactions and model spaces). In Hammaren et al. (1985) a similar method was used for odd-mass nuclei in the mass $A=130$ region with great success. These configuration-mixing of projected many-quasiparticle configurations are based on a fixed HFB function for the ground state. Changes in the deformation etc. cannot be included here. Therefore a new method was developed by Schmid et al. (1986), the "excited VAMPIR". Here the intrinsic wavefunction for the lowest energy configuration for each angular momentum (the yrast levels) is determined in a first step. In the next step, for each angular momentum, a set of excited states is determined by new intrinsic wavefunctions, Schmidt-orthogonalized with respect to all the earlier solutions at this spin. Thus one successively constructs an optimal configuration space for the $A$-nucleon problem. Over the years, all these methods were applied in relatively small configuration spaces with realistic effective forces. They compare well with the corresponding exact configuration-mixing calculations. For reviews, see Schmid and Grümmer (1987); Schmid (2004). The problem is, as in the case of the configuration-mixing calculations, the determination of the effective interaction. It depends on the specific configuration space. In many cases one starts with an effective interaction obtained by the Brückner method Day (1967); Ring and Schuck (1980) and additional parameters are adjusted to experimental data.

In Gao, Horoi and Chen (2015) calculations with variation after projection on spin, isospin, and mass number were carried out in the even-even nuclei in the sd-shell using the well known and successful USDB interaction, which was adjusted to configuration-interaction calculations in this space Brown and Richter (2006). The binding energies turn out to be very close to the exact configuration-interaction results. The differences are very small in cases, where the number of parameters in the projected wave functions is close or larger than the dimension of the configuration-interaction space. For the opposite case one finds energy differences of up to 500 keV. Angular-momentum projection is very important for these results. Calculations with angular-momentum projection lead in most of the cases to triaxial deformations of the intrinsic states. This is not the case for the calculation without $\mathit{J}$-projection. This is in agreement with systematic mean-field investigations with the Gogny-force without projection in large model spaces Delaroche et al. (2010).

Intrinsic wave functions with different deformations can be included in the generator coordinate method. In Gao and Horoi (2009) this method is extended by adding at each deformation not only the lowest Hartree-Fock configuration, but in addition a large number $n$p-$n$h configurations and projecting all these states on good angular momentum. This is a projected configuration-interaction calculation. The results are compared with exact solutions for nuclei in the sd- and pf-shell and it is found that not only the energies, but also the quadrupole moments and the B(E2) transition probabilities are in very good agreement with the exact configuration-interaction calculations in the spherical basis.

It is an interesting question, whether it is possible to define an intrinsic deformation for the exact wavefunctions in the laboratory frame obtained by configuration-mixing calculations. In cases, where one finds good agreement between the projected states and the exact solution not only for the energy, but also for other operators, the projected mean field state is close to the exact solution and in such a case it is possible to define the intrinsic deformation of the exact wavefunction as the deformation of the intrinsic state in the projected theory. It turns out that in some cases the above mentioned calculations that the intrinsic deformation is not well determined, because the projected energy surfaces are very flat and the resulting projected wave functions are nearly identical, even if the intrinsic deformations are different.

The Monte-Carlo Shell Model of the Tokyo group Otsuka et al. (2001) (see also Mizusaki and Imada (2004)) uses a fixed configuration space and it provides an exact solution identical to the conventional calculations mixing spherical configurations in an oscillator basis. Using the Hubbard-Stratonovich transformation Hubbard (1959); Stratonovich (1958) the Monte-Carlo Shell Model uses a linear combination of angular-momentum and parity projected intrinsic wavefunctions

[TABLE]

where $\sigma$ runs over a large number of intrinsic Slater-determinants $|\Phi(\sigma)\rangle$ which are determined, together with the coefficients $g_{K}^{\sigma}$, by a stochastic Monte-Carlo sampling based on the Hamiltonian $\hat{H}$. Finally the coefficients $c_{\sigma}$ are calculated by the diagonalization of $\hat{H}$. Each intrinsic state has a certain deformation and the sum of the coefficients $|c_{\sigma}|^{2}$ with a specific deformation determine the weight of the intrinsic deformation in the exact state $|\Psi_{M}^{I\pi}\rangle$. In principle, one would not need the projection operators, because the importance sampling would automatically lead to eigenstates with good quantum numbers, i.e., would automatically carry out the integration over the Euler angles and the summation over $K$. However, as it turns out, the calculations are much more stable and faster by using the projection operators. Technically, the 3-dimensional projection is possible because it is applied here in relatively small configuration spaces, including only a few oscillator shells as is usual in configuration-mixing calculations. For recent reviews on modern extensions of this method to relatively large model spaces and to ab initio applications, the reader is referred to Shimizu et al. (2012, 2017).

5.3 Projection methods in ab initio calculations

During the last few decades, ab initio calculations based on effective Hamiltonians, derived without further approximations from the bare nucleon-nucleon interaction adjusted to scattering data Wiringa et al. (1995); Entem and Machleidt (2003), gained more and more interest in nuclear structure physics. At short distances these realistic potentials have a very strong repulsive core, leading, in interaction mixing calculations, to the admixture of configurations which very high momenta. Therefore such interaction cannot be used for calculations in small model spaces. There are several methods in the literature Bogner et al. (2001, 2010) to deduce, without phenomenological parameters, from these bare forces, effective soft interactions, which can be used in small model spaces and which lead, in principle, to the same results as calculations in full space.

At the beginning these effective forces have been used for very small spaces and therefore for very light systems. Here the effective Hamiltonian were fully diagonalized by configuration mixing Barrett et al. (2013) and the results were in good agreement with exact numerical solutions of the these few-body problems Carlson et al. (2015).

For heavier systems this is no longer possible, because for calculations without core one needs a relatively large oscillator basis. Here approximate methods, based on mean field calculations with additional projection techniques, are used. The same applies to heavy nuclei where ab initio calculations are possible in valence shells in the neighborhood of magic configurations. However, already for medium heavy nuclei these valence shells become relatively large, such that full configuration-mixing calculations require considerably numerical expense both in computer time and memory. Here also mean field calculations together with projection methods based on the corresponding effective interactions are useful.

No-core shell model calculations using chiral ab initio forces have been carried out with symmetry adapted basis sets for the nuclei 6Li, 8Be, and 6He in Dytrych et al. (2013) and also for heavier nuclei in Dytrych et al. (2016); Launey et al. (2020, 2021). Such methods are also used in Dytrych et al. (2015) for ab initio calculations of electron-scattering on the nucleus 6Li. Within this symmetry-adapted framework one can achieve a considerable reduction of the dimensions as compared to large configuration-interaction calculations while retaining the accuracy of the results.

Broken symmetries and additional projection form also the basis of the Monte-Carlo Shell-Model discussed in equation (159). It has been applied for no-core shell model calculations in light nuclei and for valence shells in heavier nuclei. In recent years it has been used for effective forces derived from ab initio calculations Shimizu et al. (2012). Extensions have been successfully applied in Shimizu et al. (2017) for heavier nuclei in the neighborhood of the magic shells with N = 28, 40, and 50. In most of these calculations this model uses a basis with good particle number. Pairing correlations are included by the additional correlations in the wave function (159). However, for heavier nuclei where pairing correlations extend over a larger region around the Fermi surface, the corresponding valence shells are not large enough. Therefore, this method has been extended for calculations based on the Bogoliubov quasiparticles and the corresponding symmetry violation is restored by number projection in Shimizu et al. (2021). Recently there have been also ab initio investigations of cluster structures in nuclei based on using projection techniques Neff and Feldmeier (2008).

Coupled Cluster methods provide a very successful tool in many-body calculations in electronic systems. They are also applied in nuclear physics for ab initio calculations Hagen et al. (2014). Modern version uses basis sets with symmetry violations and additional projection, such as broken rotational symmetry Duguet (2014b) or broken gauge symmetry Signoracci et al. (2015); Duguet and Signoracci (2016); Qiu et al. (2019)

The generator-coordinate-method approaches together with symmetry projection in deformed mean-field wavefunctions have also been used for ab initio calculations with microscopic interactions derived by the Similarity Renormalization Group method for even calcium and nickel isotopes Hergert et al. (2014). Recently this method has been applied to the calculation of the nuclear matrix elements for neutrinoless double beta ($0\nu\beta\beta$) decay Yao et al. (2018, 2020). So far, this was possible only for the decay from 48Ti to 48Ca, but this method opens up new ways towards fully microscopic calculations in heavier isotopes, where this decay is under discussion. This is an example of the importance of reliable ab initio calculations in heavy nuclei not only for nuclear physics, but also for general physics, in particular the physics beyond the standard model. Earlier, the ($0\nu\beta\beta$) decay matrix elements have been calculated only with phenomenological forces of phenomenological density functionals.

6 Projection methods and nuclear density functional theory

The complexity of the bare nucleon-nucleon force and strong in-medium effects characteristic of nuclear physics lead us to consider phenomenological effective forces or relativistic Lagrangians to define the underlying intrinsic mean-field, see the reviews in Bender et al. (2003); Vretenar et al. (2005); Robledo et al. (2019); Schunck (2019). The associated HF or HFB equations are traditionally solved by expanding on a basis, as for instance the harmonic oscillator or Woods-Saxon basis, or by using mesh techniques in the coordinate representation (which is also equivalent to a basis expansion, see Baye and Heenen (1984)). The presence of a phenomenological density-dependent term containing the mean field density implies that such effective forces do not lead to mean-field average values that can be obtained from a Hamiltonian. Instead, they can be considered as a special case of an energy density functional (EDF), where part of the functional is obtained from a two-body interaction (and therefore is quadratic in the densities) and the rest is purely phenomenological. Density dependent term of this kind is also found in the Slater approximation to Coulomb exchange. These peculiarities of the traditional nuclear EDF approach lead to difficulties in the implementation of the methodology of symmetry restoration:

• •

The harmonic oscillator or Woods-Saxon bases can break spatial symmetries like rotational and/or translational. This is also the case for the mesh representation of the wave functions. This symmetry breaking has to be taken into account in the formalism explicitly.

• •

The density-dependent term of the EDF is only well defined for mean values. For energy overlaps a prescription, satisfying consistency constraints, is required. It turns out that popular alternatives lead to consider complex quantities that are often raised to non-integer powers, requiring additional considerations in the selection of the branch-cuts in the complex plane.

• •

Different interactions for the particle-hole and particle-particle channels are often used raising self-energy and violation of the Pauli principle issues.

The existence of these difficulties, and the fact that some of them are still not satisfactorily addressed, slowed down the application of symmetry restoration techniques in nuclear structure with EDFs. However, there are numerous indications that symmetry restoration is required to improve our qualitative and quantitative theoretical understanding of many nuclear properties. The results obtained so far go in this direction and some of them will be discussed below.

6.1 Difficulties encountered in restoring symmetries with nuclear EDF

6.1.1 Basis not closed under the symmetry action

Symmetry restoration requires to consider the action of a symmetry operation (rotation, translation, etc) on an intrinsic symmetry violating wave function. Often, the single-particle basis used for the intrinsic state is not closed under the action of the symmetry operator and the “rotated" basis does not span the same subspace of the Hilbert space as the original basis. In the rotational case, this happens, for instance, when the oscillator lengths along different spatial direction are not the same, or when a Cartesian mesh is used in coordinate representation Baye and Heenen (1984). In this situation, the standard formalism to compute overlaps between rotated HFB states is not valid and has to be generalized – see below. A common alternative in the applications of angular-momentum projection is to use spherical harmonic oscillator bases with a sufficiently large number of complete major shells. This strategy increases the computational cost, and becomes impractical in some situations like fission, where the large variety of shapes involved in the dynamics would require a huge rotationally invariant basis. In addition, it cannot be applied, e.g., to the case of spatial translations as required in the restoration of translational and/or Galilean invariance Schmid and Reinhard (1991); Rodríguez-Guzmán and Schmid (2004a, b). To understand the reason, let us consider a simple model in one dimension where the basis contains just one single Gaussian state $\varphi_{0}(x)=e^{-x^{2}}$. After translation $\varphi_{0}(x)\rightarrow\varphi_{0}(x-x_{0})=e^{x_{0}^{2}}e^{-x^{2}}\sum_{k=0,\infty}\frac{(-2xx_{0})^{k}}{k!}$ an infinite number of Gaussian wave functions $e^{-x^{2}}x^{k}$ is required to reconstruct the translated wave function. The situation does not improve if another basis, like Woods-Saxon or a discrete mesh in space Baye and Heenen (1984), is used.

As mentioned above, the generalized Wick’s theorem, cannot be used in this case as it is implicitly assumed in its derivation (see, for instance, Balian and Brezin (1969)) that both HFB wave functions in the overlap are expanded in the same basis. To circumvent this problem, the original finite basis (denoted by 1) is formally expanded, by using the orthogonal complement (denoted by 2), as to span the whole Hilbert space Bonche et al. (1990); Robledo (1994); Valor et al. (2000). The same procedure is applied to the “rotated" basis ($1^{\prime}$). The expanded $U$ and $V$ Bogoliubov amplitudes have a block diagonal form with the block corresponding to basis 2 having a simple form with an uniform occupancy of 0. With the introduction of the expanded basis, the traditional generalized Wick’s theorem can be applied, However, all quantities referring to basis 2 have to be reexpressed Robledo (1994); Valor et al. (2000) in terms of quantities defined in basis 1, leading to an extended Wick’s theorem that differs from the traditional one by the presence of the non-unitary overlap matrix between the two basis. Further details can be found in Robledo (1994); Valor et al. (2000); Marević et al. (2019); Marević and Schunck (2020).

6.1.2 Self-energy and Pauli principle

Reduced kernels of a two-body Hamiltonian between two HFB wave functions $|\Phi_{0}\rangle$ and $|\Phi_{1}\rangle$ can be expressed with the help of the generalized Wick’s theorem in terms of the transition density matrix $\rho^{01}$ and transition pairing tensors $\overline{\kappa}^{01*}$ and $\kappa^{01}$ as888Compare definitions in equations (77)–(79), which were valid only for the normalization specified in equation (53).

[TABLE]

where the three right-hand-side terms of equation (161) are referred to as direct, exchange, and pairing terms. As the overlap $\langle\Phi_{0}|a^{+}_{k_{1}}a^{+}_{k_{2}}a_{k_{4}}a_{k_{3}}|\Phi_{1}\rangle$ is a finite quantity, the right hand side of equation (161) must diverge when the overlap $\langle\Phi_{0}|\Phi_{1}\rangle$ vanishes. The same argument applies to $\langle\Phi_{0}|a^{+}_{k_{1}}a_{k_{3}}|\Phi_{1}\rangle$, etc and therefore $\rho^{01}_{k_{3}k_{1}}$, $\overline{\kappa}^{01*}_{k_{1}k_{2}}$ and $\kappa^{01}_{k_{3}k_{4}}$ are also divergent quantities when the overlap vanishes. The order of the pole present in all the quantities, including the Hamiltonian kernel, must be the same and equal to the order of the zero of the overlap. Therefore, in the right hand side of equation (161) there must be a cancellation of poles, so as to reduce the order of the pole in the products of transition densities and pairing tensors. The cancellation indeed takes place and is a direct consequence of the Pauli exclusion principle Stringari and Brink (1978). The same kind of cancellation also happens for the three-body and many-body operators. This represents a serious problem if some of the contributions (typically exchange and/or pairing) for some interaction terms are neglected. In Tajima et al. (1992) a somehow arbitrary “regularization" scheme was introduced to define “regularized" contributions for the direct, exchange and pairing terms. In Dönau (1998), the problem was also discussed, raising serious doubts about the validity of calculations with the pairing plus quadrupole Hamiltonian (see section 5), where the exchange and pairing parts of the quadrupole-quadrupole force, as well as the direct and exchange parts of the pairing force, are all neglected.

In Anguiano et al. (2001) the same problem was discussed in the context of the particle-number projection with the Gogny force, see also Dobaczewski et al. (2007); Duguet et al. (2009); Bender et al. (2009); Hupin et al. (2011); Duguet et al. (2015) for follow-up studies with the Skyrme forces. In traditional studies with the Gogny force, the Coulomb exchange is treated in the Slater approximation, and Coulomb and spin-orbit pairing are neglected. As a consequence, unphysical and inconsistent values of the projected energies were obtained. The use of variation after projection aggravated the problem, as the projected energy could get attractive contributions from the unphysical contributions. Due to the simple form of the overlaps involved in the particle-number projection, it was possible to trace back the problem to the existence of vanishing overlaps at the gauge angle $\varphi=\pi/2$, due to configurations with occupancy $v_{k}^{2}=1/2$. The inclusion of all the missing contributions solved the problem. The computational cost can increase by up to two orders of magnitude if the exchange and pairing Coulomb contributions are required in a calculation with contact central potentials.

To illustrate the problem, let us compute the general matrix element of equation (161) with $|\Phi_{1}\rangle=\exp(-i\varphi\hat{N})|\Phi_{0}\rangle$ The transition density matrix and pairing tensor take a very simple form in the canonical basis of the Bogoliubov transformation (see Mang (1975); Ring and Schuck (1980) for details). If $c^{+}_{k}$ represent the creation operators in the canonical basis of the HFB wave function $|\Phi_{0}\rangle$, we obtain

[TABLE]

When $\varphi=\pi/2$ ($e^{-i2\varphi}=-1$) and the occupancy of one of the states $k$ is $u_{k}^{2}=v_{k}^{2}=1/2$, some of the denominators in this expression go to zero. In this case, the overlap $\langle\Phi_{0}|\exp(-i\varphi\hat{N})|\Phi_{0}\rangle=\prod_{k>0}(u_{k}^{2}+v_{k}^{2}e^{-2i\varphi})$ also goes to zero. As long as there is a single pole, or two poles but with $k_{1}\neq k_{2}$ (or $k_{3}\neq k_{1}$) the overlap in the numerator of the l.h.s. of equation (162) remains finite. The only problematic case is $\langle\Phi_{0}|c^{+}_{k}c^{+}_{\bar{k}}c_{\bar{k}}c_{k}\exp(-i\varphi\hat{N})|\Phi_{0}\ \rangle$ where the right-hand side of equation (162) apparently has a pole of order two that cannot be canceled out by the norm overlap. Although each of the two terms have an order two pole, the sum turns out to be a order one pole. When multiplying by the norm overlap, the sum gives a finite contribution, namely $v_{k}^{2}\cdot e^{-2i\varphi}\cdot\prod_{m>0,m\neq k}(u_{m}^{2}+v_{m}^{2}e^{-2i\varphi})$. As long as the direct, the exchange, and the pairing terms of the Wick’s factorization of equation (161) are taken into account, no divergences appear in the particle-number projection formalism. This cancellation is also connected with the so-called “self interaction" problem Perdew and Zunger (1981); Lacroix et al. (2009): the sum of the three terms is required to give a zero two body energy for a Slater determinant containing just one particle.

Although the solution to the self-energy problem is simple, and it is routinely used with Gogny forces (see section 6.2), it is not possible to implement it for those nuclear EDFs where the particle-hole and particle-particle interactions are independent of each other. In this case, by construction, the particle-particle interaction do not contribute to the particle-hole channel and vice versa and therefore the pole cancellation cannot take place. The problem was extensively discussed in the literature Dobaczewski et al. (2007); Duguet et al. (2009); Bender et al. (2009); Hupin et al. (2011); Duguet et al. (2015) and attempts to renormalize the divergences were proposed in Tajima et al. (1992); Lacroix et al. (2009); Satuła and Dobaczewski (2014). Unfortunately, they only work for functionals depending on integer powers of the density and a clear renormalization criteria is missing in the most general cases. This difficulty is the major reason for the recent interest in replacing EDFs by mean values of true Hamiltonians Bennaceur et al. (2014); Sadoudi et al. (2013); Bennaceur et al. (2017).

Note that the arguments given above are rather general, and thus they also apply to restoration of other symmetries like angular momentum Zduńczuk, Satuła, Dobaczewski and Kosmulski (2007). They do not depend on the kind of force, or the kind of two body operator considered in the mean value. Only for strictly density-independent forces (see section 6.1.3) like Skyrme SV Beiner et al. (1975) and no pairing, all symmetries, including the isospin symmetry, can be consistently restored Satuła et al. (2010). Obviously, the projected calculation of mean values of one-body operators are never affected by the inconsistencies mentioned above.

6.1.3 Density-dependent prescription

Most of the EDF used nowadays include in one way or another a phenomenological density dependence. In the non-relativistic case, it is introduced to mock-up the saturation property of the nuclear interaction in a simple way. Although the saturation could also be obtained without any density dependence Beiner et al. (1975), this inevitably leads to an unrealistic low value for the effective mass Davesne et al. (2018). In the relativistic case, the saturation property is a relativistic effect due to the difference between the scalar density, the source of the attractive part of the force and the baryon density, the source of the repulsive part of the force. Here, the saturation can also be obtained without any density dependence Walecka (1974). Nonetheless, to obtain a realistic description of the nuclear incompressibility and surface properties, an additional density dependence is used in the coupling constants.

In both the cases, density-dependent forces are implemented in a state-dependent way, so that the density $\rho(\bm{r})=\langle\Phi|\hat{\rho}(\bm{r})|\Phi\rangle$ is used to evaluate the HFB energy associated with the HFB wave function $|\Phi\rangle$. For most of the non-relativistic functionals, the density-dependent contribution to the energy is strongly repulsive, it is usually proportional to $\delta(\bm{r}_{1}-\bm{r}_{2})$ (contact term) and depends on the intrinsic center of mass density $\rho((\bm{r}_{1}+\bm{r}_{2})/2)$ raised to some power $\alpha$ which is, in most of the cases, a non integer number (usually 1/3, although 1/6 is also a common choice). The typical form of the density-dependent part of functional reads

[TABLE]

where $t_{3}$ and $x_{3}$ are parameters and $\hat{P}_{\sigma}$ is the spin-exchange operator, see Bender et al. (2003) for details.

When the variational principle is used to derive the HF or HFB equations, the density-dependent part has also to be varied, leading to the so-called rearrangement term in the mean field, which in the non-relativistic case is given by

[TABLE]

Given the typical values of $t_{3}$, this rearrangement term represents an important component of the HF or HFB mean field and it cannot be overlooked. It enters in the expressions of many quantities like chemical potentials, quasiparticle excitation energies, etc..

In the evaluation of symmetry-restored energies based on HFB intrinsic states, instead of the mean values of the HFB theory we have to consider kernels of the Hamiltonian between different HFB states. In the case of density-dependent forces some kind of prescription has to be given for the evaluation of the Hamiltonian kernels. This requirement also arises in the implementation of the generator coordinate method.

There are a few fundamental requirements that must be used to guide the choice of the prescription:

• •

The projected energy has to be a real quantity.

• •

The projected energy has to remain invariant with respect to symmetry operations applied to the intrinsic state.

• •

In the strong deformation limit of the intrinsic state (see section 4.3), the projected theory should reduce to the traditional mean-field approach Robledo (2007).

Among the various proposals made over the years Bonche et al. (1990); Valor et al. (2000); Rodríguez-Guzmán et al. (2002b); Duguet and Bonche (2003); Schmid (2004), we only comment on the two that are in vogue:

Mixed-density prescription:

The “mixed" or “transition" density prescription, proposed already in early nineties Bonche et al. (1990) within the framework of the generator coordinate method, and extended to the symmetry restoration case in Valor et al. (2000); Rodríguez-Guzmán et al. (2002b), employs the transition density,

[TABLE]

where $\hat{\rho}$ is the one-body density operator. Then, $\rho_{q,q^{\prime}}$ is used in the evaluation of the Hamiltonian kernel

[TABLE]

where $\hat{H}\{\rho\}$ denotes any density-dependent Hamiltonian containing terms similar to $V_{DD}$ of equation (163). This prescription is inspired by contractions used in the generalized Wick’s theorem when it is applied to three-body forces. It satisfies all the consistency requirements mentioned above Rodríguez-Guzmán et al. (2002b); Egido and Robledo (2004). However, since the transition density is in general a complex quantity, this prescription requires additional rules when it is implemented for the density raised to a non-integer power.

Projected-density prescription:

The density to be used in the Hamiltonian kernel (166) is given by

[TABLE]

where $\hat{P}$ represents generic projection operator (50). This prescription produces a density-dependent term which is invariant under the broken symmetries. For instance, in the case of the rotational symmetry breaking, this prescription leads to density $\rho^{P}_{q,q^{\prime}}$ that is spherically symmetric, and because of this, it was advocated by Schmid (2004). This prescription satisfies the first two requirements mentioned above, but the third one is more problematic.

The projected density prescription is highly inconsistent with the underlying intrinsic mean field Robledo (2010) when spatial symmetries like the parity (see discussion in section 2.1) or rotational invariance are under consideration. An example of the catastrophic consequences can be seen in figure 11, where the parity-projected energy for the nucleus 224Ra is plotted as a function of the octupole moment $Q_{30}$ for the two prescriptions discussed here. For the projected-density prescription, with increasing octupole moment, the projected energy continuously decreases as a consequence of the mismatch between the projected density and the transition density.

We note here, that in the case of the particle-number projection, the use of the projected-density prescription produces reasonable results. This is related to the fact that the particle-symmetry restoration affects only the occupation numbers of single-particle canonical states and does not lead to a spatial mismatch between the left and right states for which the Hamiltonian kernel (166) is evaluated. The Madrid group routinely uses a hybrid prescription, that is, the mixed-density prescription is used for all symmetries apart from the particle-number symmetry restoration when the projected-density prescription is used. In addition, as it was shown in Valor et al. (2000), the Lipkin-Nogami approximation (see section 4.2) leads to results that are very similar to those obtained using the mixed-density prescription. We also note that the random phase approximation applied to the generator coordinate method Jancovici and Schiff (1964) can also be used as an argument in favor of the transition-density prescription Robledo (2007).

6.1.4 Non-integer powers of the density

The transition-density prescription requires to consider the transition density raised to some non-integer power. However, the transition density is a complex quantity, in general, and we are confronted with the problem of how to consider the evaluation of a non-integer power of a complex number. The choice of branch cut in this case can have an enormous impact on the value of the matrix element as a consequence of the large positive value of the strength of the density-dependent term.

In addition, in the symmetry restoration case, where integrals over the parameters of the symmetry group have to be carried out, the presence of branch cuts associated with the density-dependent part break the analyticity of the integrand and leads to spurious dependencies on the integration path in the complex plane of the symmetry group parameters. In a simple case of the particle-number restoration, this issue was discussed in great detail by Dobaczewski et al. (2007); Duguet et al. (2009). Similar problems are expected for the restoration of other symmetries. Treatment of the density dependence with non-integer powers of the density is still an unsolved problem and represents a serious limitation in the applicability of the symmetry restoration and configuration-mixing methodology within nuclear EDF.

6.1.5 Future implementations

In recent years, the idea to give up phenomenological density-dependent terms in favor of real (multibody) operators is becoming increasingly popular in non-relativistic density functional theory Bennaceur et al. (2014); Sadoudi et al. (2013); Bennaceur et al. (2017). The idea is to find a three-body interaction (not excluding higher order terms) that is able to mimic the saturation property induced by the density-dependent interaction and at the same time is not spoiling pairing properties. Plain use of the generalized Wick’s theorem would be the only thing required to compute the kernel overlap. The difficulty encountered with this approach is the large number of possible terms with their associated free parameters. In order to reduce the number of terms, several assumptions were made, but so far, the three-body interactions proposed are not very accurate in describing nuclear properties.

On the other hand, the influence of the poles in the projected energy surfaces seem to be relatively narrow and in many applications, in particular in a projection after variation, they can be regularized by using a wider mesh in the integration over the various angles. In this way, the results depend on the choice of the discretization mesh, but at the end the numerical errors connected with such procedures do not play an essential role, in particular for heavy nuclei. In any case, one finds in the literature many very successful applications of generator coordinate method and projection after variation (see below), where the potential impact of such poles in the physical observables do not spoil the physical interpretation.

6.2 Applications of symmetry restoration with nuclear EDFs

The importance of symmetry restoration in nuclear structure was noticed very early in nuclear physics Bayman (1960); Dietrich et al. (1964). Applications of this method to simple models were already discussed in sections 3.2.1 and 5. Other approaches used shell-model-like effective interactions defined in a restricted configuration space involving a limited number of orbits. Among them we can mention the work developed by the Tübingen group using very sophisticated many body techniques involving mixing of configurations through projecting out mean-field intrinsic states obtained in a variation after projection framework Schmid (2004). We will not dwell on these approaches as we are more concerned with realistic nuclear energy density functionals.

6.2.1 Non relativistic EDFs

The first application of symmetry restoration with an EDF possibly dates back to the work of Caurier and Grammaticos Caurier and Grammaticos (1977) where rotational bands of light nuclei were computed using several flavors of the Skyrme interaction and angular-momentum projection. A few years later, the parity-projected excited intrinsic configurations of 20Ne were studied Marcos et al. (1983) using the BKN interaction Bonche et al. (1976). In that work, the center of mass correction was also computed using symmetry restoration techniques. Reflection symmetry restoration was also used along with the Gogny D1S interaction in Egido and Robledo (1991) to describe the physics of parity doublets and octupole deformation. Parity restoration was also applied to cranking wave functions in order to study the emergence of octupole deformation at high spins Garrote et al. (1997, 1998). Parity projection of non axial intrinsic states was performed in several mercury and lead isotopes in a generator coordinate context Skalski et al. (1993). Systematic calculations of the excitation energies, E3 transition strengths and ground state octupole correlation energies were performed in a combined generator coordinate method and parity projection framework Robledo and Bertsch (2011); Robledo (2015a) with several flavors of the Gogny force, depicting the importance of the dynamical octupole correlations in nuclear structure.

To describe collective negative parity states, intrinsic wave functions breaking reflection symmetry are important ingredients and therefore octupole deformation becomes important in this case. At the mean field level only a few nuclei in the actinide and rare earth regions are octupole deformed. The number of nuclei accessible increases substantially if the intrinsic state is determined in the variation after parity projection scheme, irrespective of whether it is implemented exactly or in a restricted variational space Robledo and Bertsch (2011); Robledo (2015b). Octupole correlations are better described using the generator coordinate method with the octupole moment as one of the collective coordinates Bernard et al. (2016); Bucher et al. (2017); Bernard et al. (2017).

Particle-number restoration is another important application of symmetry restoration in nuclear structure as pairing correlations are known to be rather weak in atomic nuclei and therefore the use of an intrinsic mean field wave function (BCS or HFB wave function) is not easy to justify. Fluctuations of the order parameter associated with pairing correlations (for instance, the fluctuation of the particle number $\langle\Delta N^{2}\rangle$) as well as the corresponding gauge angle associated with particle-number restoration are important ingredients for a proper description of pairing correlations in the weak-pairing regime of nuclear physics. Most of the applications use intrinsic states obtained from a HFB calculation, supplemented with the Lipkin-Nogami procedure (see section 4.2) generalized for density dependent forces Valor et al. (1997). There are many examples of full variation after particle-number projection calculations with the Gogny force, mostly in the framework of a particle-number projected generator coordinate method with restoration of additional symmetries. There are also early examples Anguiano et al. (2001) aimed to understand the effect of particle-number projection on the moment of inertia of rotational bands or in the ground state correlation energy Anguiano et al. (2002). Particle-number projection in the variation after projection scheme increases pairing correlations and, as a consequence, decreases the moment of inertia, increasing thereby the excitation energy of rotational $2^{+}$ states.

The structure of some Sr Heenen et al. (1993) or Pb isotopes Heenen et al. (2001) was analyzed with Skyrme EDF and the particle-number projected generator coordinate method. In Stoitsov et al. (2003), a mass table from proton to neutron dripline was generated with the SLy4 EDF and using volume pairing and implementing Lipkin-Nogami method followed by a full particle-number projection. The procedure was implemented in a computer code that is publicly available Stoitsov et al. (2005). Similar calculations were performed by the Brussels group in their quest for an accurate mass model Samyn et al. (2004).

A full variation after particle-number projection calculation with a Skyrme functional plus a zero range pairing force was carried out in Stoitsov et al. (2007) using the formulation of the particle-number projection method of Sheikh and Ring (2000), solely involving functions of the standard density and the abnormal pairing tensor, see Sect 3. The results were compared to the ones obtained with the Lipkin-Nogami method followed by a subsequent projection on particle number and a substantial improvement is observed, specially for magic or near magic nuclei.

Spontaneous symmetry breaking of rotational invariance is a defining characteristic of the nuclear interaction. It leads to the fruitful concept of nuclear deformation Casimir (1935) that allows us to explain a variety of phenomenology like, for instance, the ubiquitous existence of rotational bands in nuclear spectra Bohr and Mottelson (1953, 1955). Most of nuclei in the Segré chart are thought to exhibit rotational symmetry breaking in some of their quantum states, either the ground or excited states. Although it is possible to extract a lot of information out of the intrinsic deformed states by using the strong deformation limit [see Mang (1975); Ring and Schuck (1980) and section 4.3], the existence of weakly deformed states and/or the coexistence of different types of deformations in a limited range of energies requires the explicit restoration of the rotational symmetry. This is also the case for calculations of electromagnetic transition strengths in weakly-deformed nuclei Robledo and Bertsch (2012).

One of the difficulties of the restoration of angular momentum is the three-dimensional integration over the Euler angles. Often, the assumption of axial symmetry is made to reduce the number of integrals to just one, facilitating the application of the method, see section 3.2.5. Axial angular-momentum projection along with particle-number projection was first carried out in Valor et al. (2000) with several parameterizations of the Skyrme EDF. This paper, together with Heenen et al. (1993), contains a detailed description of the evaluation of operator matrix elements and different prescriptions for the density-dependent part of the interaction. Further applications include the study of the impact of the rotational energy correction in fission barriers Bender, Heenen and Bonche (2004) or the analysis of collective low lying structures in Kr Bender, Bonche and Heenen (2006) or Pb isotopes Bender, Bonche, Duguet and Heenen (2004). Other applications include the study of bubble structures in 34Si Yao et al. (2012) (see figure 12) that shows the relevance of angular-momentum projection for the description of matter density. The intrinsic density shows a bubble structure at the center that is much less pronounced for the density obtained from the projected ground-state wave function.

Applications with the Gogny D1S energy density functional include the work of Rodríguez-Guzmán et al. (2000a) where the erosion of the $N=20$ magic number in the magnesium isotopic chain was addressed. The physics of super-deformation in sulfur isotopes was also analyzed Rodríguez-Guzmán et al. (2000c) in this framework. The calculations were further extended to consider the generator coordinate method with the quadrupole degree of freedom as generating coordinate to explain more quantitative features of the deformation of Mg isotopes Rodríguez-Guzmán et al. (2000b, 2002b), the erosion of the $N=28$ shell closure Rodríguez-Guzmán et al. (2002a) and triple shape coexistence of neutron deficient lead isotopes Rodríguez-Guzmán et al. (2004).

Calculation of mass tables of even-even nuclei including angular-momentum projection restricted to axial symmetry plus particle-number projection and including generator coordinate method for the axial quadrupole degree of freedom with Skyrme Bender et al. (2005); Bender, Bertsch and Heenen (2006); Bender et al. (2008) or Gogny forces Rodríguez et al. (2015) were presented in the literature. Other large scale studies include the excitation energy of the collective $2^{+}$ state and its B(E2) transition strength to the ground state Sabbey et al. (2007); Rodríguez et al. (2015). Applications of this methodology to the study of neutrino-less double beta decay are essential to extract relevant nuclear matrix elements in medium mass and heavy nuclei Rodríguez and Martínez-Pinedo (2010), see results presented in figure 13.

Angular-momentum projection of triaxial intrinsic states was first carried out in an early work by Caurier and Grammaticos Caurier and Grammaticos (1977) with the Skyrme SIII interaction. A restricted variation after projection was used with the radii and quadrupole moments as variational quantities. A very small configuration space was used and therefore the applications were restricted to very light nuclei. Quite a few years later, Bonche et al. (1991) used conveniently chosen combinations of intrinsic states oriented along different axes to carry out an approximate projection of triaxial intrinsic states on $I=0^{+}$. A full triaxial angular-momentum projection was carried out with Skyrme interactions in Baye and Heenen (1984), but restricted to HF wave functions. The HF wave functions were represented on a mesh of nodes in coordinate representation and special care was taken to define the rotation operator in that case. In Bender and Heenen (2008) triaxial angular-momentum plus particle-number projection of HFB states was performed with Skyrme SLy4 for light nuclei. The main conclusion was that the inclusion of triaxiality improves the description of rotational bands as compared to axial results. The calculations were subsequently extended to heavier systems in Yao et al. (2013). Using the finite range Gogny forces, the technology to project triaxial states was developed in Rodríguez and Egido (2010) and applied to the study of 44Si in Rodríguez and Egido (2011a) and to the waiting point nucleus 80Zr in Rodríguez and Egido (2011b).

Most of the angular-momentum projected calculations so far were of the projection after variation type with intrinsic states preserving time-reversal invariance. As a consequence, the rotational bands obtained were stretched with respect to experiment by a typical factor of 1.4 Rodríguez-Guzmán et al. (2002a); Li, Nikšić, Ring, Vretenar, Yao and Meng (2012), a consequence of implicitly using the Peierls-Yoccoz moment of inertia instead of the Thouless-Valatin one, see section 4.4.1. In order to overcome this difficulty, the use of time-reversal breaking intrinsic states is required.

Cranked HF states were employed as early as in Baye and Heenen (1984) with a simple Skyrme like interaction (BKN+Coulomb). They were also used as intrinsic states in a full three-dimensional angular-momentum projection with the Skyrme SLy4 force to analyze (i) band termination in nuclei around $A=44$ Zduńczuk, Satuła, Dobaczewski and Kosmulski (2007) and (ii) angular-momentum projection in 156Gd and 155Er Zduńczuk, Dobaczewski and Satuła (2007). Later, calculations with HFB cranking intrinsic states were carried out with the Gogny force in Borrajo et al. (2015); Egido et al. (2016) [see also Rodríguez (2016) for an analysis of the different moments of inertia obtained in the different approaches]. In Egido et al. (2016), collective and single-particle degrees of freedom were studied in the nucleus 44S. Intrinsic wave functions $|\Phi(\beta,\gamma,\omega)\rangle$ of the HFB type with quadrupole deformation parameters $\beta$ and $\gamma$, and obtained at different angular frequencies $\omega$ were used in a configuration mixing calculation including projection on good particle number and angular momentum (see figure 14). The consideration of cranking states improves notably the description of moments of inertia, whereas inclusion of states with deformation $\gamma$ in the full interval of $-60^{\circ}\leq\gamma\leq 120^{\circ}$ allows for considering two-quasiparticle excitation like those that are present in the experimental spectrum of 44S.

The use of the cranking frequency $\omega$ as a generator coordinate was considered in Shimada et al. (2015, 2016) without particle-number projection. In these calculations, the intrinsic cranking state was rotated along a single axis. In contrast, in Tagami and Shimizu (2016) intrinsic states were generated by using infinitesimal cranking frequencies $\omega_{i}$ along three different cranking axes. This technique was well suited to describe $\gamma$ bands and was applied to study these structures in 164Er. Applications to wobbling motion in odd-$A$ nuclei and to chiral doublet bands were presented in Shimada et al. (2018a) and Shimada et al. (2018b), respectively.

Odd mass nuclei also require time-reversal breaking of intrinsic states for their analysis. In this case, the existence of zeros in the overlaps between rotated wave functions seem to be more likely than when time-reversal is preserved Oi and Tajima (2005); Zduńczuk, Dobaczewski and Satuła (2007). It is therefore to be expected that the problems concerning self-energies and ill-behaved density-dependent forces should be more relevant for odd mass nuclei. This consideration led Bally et al. (2014) to consider for their first projected calculation of an odd mass nucleus a Skyrme interaction which is fully derived from a Hamiltonian Sadoudi et al. (2013). The results for 25Mg look reasonable, see figure 15. Surprisingly, the implementation with the density-dependent Gogny force Borrajo et al. (2015); Borrajo and Egido (2018) does not lead to any apparent inconsistency in the Mg isotopes considered. Some results in heavier systems and for high spin states were discussed in Shimada et al. (2018a) and magnetic moments in 45Sc were determined in de Groote et al. (2020).

Isospin symmetry is explicitly broken in the atomic nucleus due to the Coulomb interaction and the tiny differences observed in the different isospin channels of the nucleon-nucleon interaction originating from the different mass and charge of the $u$ and $d$ quarks. In addition, when working at the mean field level, isospin symmetry is also spontaneously broken by the mean field wave functions Engelbrecht and Lemmer (1970); Brink and Svenne (1970). Independently of the origin of the broken symmetry, working in a basis preserving isospin quantum numbers is advantageous to understand the impact of the different sources of isospin symmetry breaking in the nuclear wave function. The first application of isospin projection in a variation-after-projection framework was carried out in Caurier and Poves (1982) using a density-independent interaction with a Brink-Boecker central potential and applied to the study of Coulomb displacement energies. Later, the formalism was applied to Skyrme functionals in Satuła et al. (2009, 2010, 2012) for the case of HF wave functions not mixing protons and neutrons at the single-particle level.

The proper treatment of isospin related effects often involves the use of EDF mixing proton and neutron densities as well as cranking techniques in the isospin space (isocranking) Sato et al. (2013); Sheikh et al. (2014).

The traditional applications of symmetry restoration for EDF are combined with the use of the generator coordinate method with collective continuous coordinates (deformation parameters, pairing gaps, etc) so as to describe low-energy collective states. An alternative is to use multi-quasiparticle excitations as discrete generating coordinates (like in the Projected Shell Model discussed in section 5.1) to gain flexibility in generating the correlated wave function. Both the projected shell model and the Shell Model employ restricted configuration spaces that require the introduction of a core as well as effective charges. Recently, a No Core Configuration Interaction (NCCI) method was implemented along with the density-independent SV Skyrme interaction Satuła et al. (2016). The method uses the full configuration space, removing the need for a core and/or effective charges. The present implementation of NCCI used p-h excitations of Slater determinants projected to good angular momentum and isospin and was employed to study excitation spectra of several $N\approx Z$ nuclei as well as $\beta$-decay and exotic-processes matrix elements Konieczka et al. (2016, 2018).

6.2.2 Relativistic EDFs

The relativistic mean-field approach, also called Covariant Density Functional Theory, represents an alternative approach to describe the structure of the nucleus. Its main ingredient is the Dirac equation, which is used to determine the nucleon orbits. The potentials entering Dirac’s equation are deduced in different ways depending on the version of the relativistic model used. In most of these models, the potentials experienced by nucleons depend upon several meson fields Reinhard (1989); Ring (1996) which are determined through classical inhomogeneous Klein-Gordon equations, where the sources are given in terms of the nucleon densities and currents. The simplest version of this model Walecka (1974) cannot reproduce the right incompressibility of nuclear matter. Therefore a density dependence was introduced by non-linear meson couplings Boguta and Bodmer (1977) or by an explicit density dependence of the meson-nucleon coupling constants Lalazissis et al. (2005).

In deformed nuclei the classical meson fields $\phi_{i}({\bm{r}})$ are deformed. In order to restore symmetries the meson fields have to be quantized using bosonic creation and annihilation operators, $b_{i}^{\dagger}({\bm{r}})$ and $b_{i}({\bm{r}})$. The total wave function $|\tilde{\Phi}\rangle$ of the system becomes the product of a Slater determinant $|\Phi\rangle$ in Fermion space and a coherent state in boson space, that is,

[TABLE]

A variation of the corresponding energy with respect to the single-particle wave functions of the fermions and with respect to the meson field leads to the classical relativistic mean-field equations. For the angular-momentum projection discussed in section 3.1.1, one has to evaluate integrals of norm and Hamiltonian overlap kernels, equations (52) and (54), not only in the fermion space, but also in the boson space. Expressions pertaining to the boson space can be found in Balian and Brezin (1969). However, because of the numerical complexity of the Hamiltonian matrix elements with finite-range interactions of Yukawa-type, such calculations were not carried out so far.

A simple way to bypass these problems is to use the relativistic point coupling models. Here the meson propagators with the large meson masses are expanded in momentum space up to second order in $q/m_{i}$, where $q$ is the momentum transfer and $m_{i}$ are meson masses, and one ends up with a Lagrangian without mesons, but containing zero-range fermion interactions and zero-range derivative terms, in full analogy to the non-relativistic Skyrme functionals Bürvenich et al. (2002). The large repulsive density-dependent contact term of the non-relativistic case is not needed here. However, for a good description of nuclear matter properties a density dependence is introduced in the Lagrangian by either three- and four-body contact terms Bürvenich et al. (2002); Nikšić et al. (2008); Zhao et al. (2010) or by density-dependent coupling constants of the two-body contact terms Nikšić et al. (2008). This leads to coupling constants depending on integer powers of density only, and therefore the difficulties mentioned in section 6.1.4 related to a non-analytical dependence on complex densities are not present. The self-energy and self-pairing problem of section 6.1.2 is however present because, in practice, the Fock terms, which are also of zero range, are usually neglected.

A much more serious problem is, however, related to the pairing channel. Because of the extremely large relativistic scalar and vector fields, one cannot use the same force in the HF and in the pairing channel Kucharek and Ring (1991). As a result, in the relativistic point coupling models it is a common practice to use pairing interaction different from that in the particle-hole channel. Therefore, the self-energy problems associated with the violation of Pauli principle are present in these calculations and cannot be easily avoided. The conclusion is that, as in most of the non-relativistic cases, the results obtained by symmetry restoration can contain some spurious contamination and their stability with respect to the parameters of the calculations should always be carefully checked.

Point coupling models with many-body contact terms and without explicit density dependence were used in the mean-field + BCS approximation for beyond mean-field calculations with symmetry restoration (and configuration mixing) in Nikšić et al. (2006a), with angular-momentum projection of axially symmetric intrinsic states, and also in Nikšić et al. (2006b) with simultaneous angular-momentum and particle-number projection. Angular-momentum projection was carried out after variation. The particle number was treated in the intrinsic state with the Lipkin-Nogami method (see section 4.2) and after that an exact particle-number projection was carried out. In this way it was possible to provide a microscopic description of the $X(5)$ quantum-phase transition, which was introduced in a group-theory model by Iachello (2001). It describes a transition from spherical to axially symmetric deformed intrinsic shapes and it is realized, e.g., in the chain of Nd-isotopes between the spherical nucleus 142Nd and the axially deformed nucleus 152Nd. At 150Nd, a first order phase transition occurs and this nucleus has a very characteristic spectrum (see figure 16), which can be described in the $X(5)$ model with only two phenomenological parameters.

In Nikšić et al. (2007), collective states in 152Nd were determined within a fully microscopic relativistic mean-field calculation with subsequent angular-momentum and particle-number projection. The entire calculated spectrum was subsequently scaled to match the experimental $2^{+}$ energy, see section 6.2.1. The resulting spectrum is even in a better agreement with experiment than the group-theory spectrum, which indicates that the nucleus 150Nd is not exactly at the transition point of the $X(5)$ model. The transition rates were calculated in the full configuration space without effective charges and show excellent agreement with the experimental data.

The model based on Covariant Density Functional Theory was further extended to describe (i) triaxial intrinsic states Yao et al. (2009, 2010, 2011), (ii) reflection asymmetric states in Yao, Zhou and Li (2015); Zhou et al. (2016); Yao and Hagino (2016), (iii) admixtures of projected two-quasiparticle configurations Zhao et al. (2016) important for band-crossing phenomena in rotating nuclei, and (iv) nuclear matrix elements for $0\nu\beta\beta$ decay Song et al. (2014); Yao, Song, Hagino, Ring and Meng (2015). In Yao and Engel (2016), the $0\nu\beta\beta$ decay of 150Nd to 150Sm was studied by including octupole correlations in the description of the ground and lowest-lying $0^{+}$ collective excited states. In figure 17, the result obtained for the nuclear matrix element corresponding to the $0^{+}_{1}\rightarrow 0^{+}_{1}$ $0\nu\beta\beta$ transition is compared with predictions including octupole correlations and with those of other similar calculations Rodríguez and Martínez-Pinedo (2010). The inclusion of octupole correlations in the ground states of both mother and daughter nuclei reduces the value of $M^{0\nu}$ by 7%, which, however, is not enough to reduce discrepancies with non-relativistic calculations using a similar framework, or with other calculations using the quasiparticle random phase approximation or the interacting boson model.

6.3 Approximate projection for nuclear EDF

As we discussed in section 6.2, practical applications of projection and configuration-mixing methods within nuclear-EDF approaches are quite successful. However, large configuration spaces make applications of the variation after three-dimensional angular-momentum projection in heavy nuclei quite difficult. Many of the applications are therefore restricted to light nuclei, where one has in principle also other methods such as configuration-interaction calculations or coupled-cluster methods. Methods based on the mean-field approximation are assumed to work better in heavy systems, where other methods cannot be applied. Therefore, for heavy nuclei, approximate methods were developed for nuclear EDFs. They are based on the fact, that the overlap $\langle\Phi(q)|\Phi(q^{\prime})\rangle$ and Hamiltonian $\langle\Phi(q)|\hat{H}\{\rho_{q,q^{\prime}}\}|\Phi(q^{\prime})\rangle$ kernels between two different HFB wave functions, cf. equation (166), are sharply peaked at $q=q^{\prime}$.

For heavy systems, the Gaussian Overlap Approximation Ring and Schuck (1980) is well justified. It was shown Haff and Wilets (1973); Giraud and Grammaticos (1974) that under this approximation one can derive a collective Hamiltonian in collective variables $q$. It contains a potential energy $V(q)=\langle\Phi(q)|\hat{H}\{\rho_{q,q}\}|\Phi(q)\rangle$, a kinetic term with microscopically derived inertia parameters and zero-point corrections [for details see Libert et al. (1999)]. In the case of three-dimensional angular-momentum projection, one ends up with the rigid-rotor Bohr Hamiltonian Une et al. (1976), where angular momentum is automatically preserved. For the generator-coordinate-method ansatz that includes quadrupole-deformation parameters $\beta$ and $\gamma$, one finds in this approximation the five-dimensional rotation-vibration Bohr Hamiltonian Próchniak and Rohoziński (2009). A similar collective Hamiltonian can also be derived within the adiabatic time-dependent HF theory Baranger and Vénéroni (1978).

The advantage of above approximation is that one only has to solve the constrained mean-field equations on the energy surface characterized by the parameters $q$ and to determine the expectation values of certain operators, e.g., $\langle\Phi(q)|\hat{H}\{\rho_{q,q}\}|\Phi(q)\rangle$ or $\langle\Phi(q)|\hat{H}\{\rho_{q,q}\}\hat{J^{2}}|\Phi(q)\rangle$. One avoids the complicated matrix elements and the problem of singularities connected with those.

As an example, in figure 18, we show the results of benchmark calculations by Yao et al. (2014), where full three-dimensional angular-momentum and particle-number projected generator-coordinate-method calculations are compared with experiment and with the results of the corresponding five-dimensional Bohr Hamiltonian (see also another example in Delaroche et al. (2010)). The agreement between the two calculations for this complicated spectrum in the transitional nucleus 74Kr is excellent. Having in mind, that the generator-coordinate-method calculations for this spectrum required 200 CPU hours with one processor, it is easy to understand that nowadays one can find many applications based on this approximation. Unfortunately, as discussed in Rodríguez et al. (2015), the method also has its downside related to negative values of zero-point-energy corrections.

7 Projection methods in other mesoscopic systems beyond atomic nuclei

The last few decades witnessed extraordinary advances in experimental techniques leading to the fabrication of mesoscopic and nanoscopic many-body systems with unparalleled control and diversity over the finite number of constituent particles, temperature, interparticle interactions, dimensionality, particle density, statistics (fermions versus bosons), and spin Ellenberger et al. (2006); Hanson et al. (2007); Serwane et al. (2011); Zürn et al. (2012); Murmann et al. (2015); Tai et al. (2017); Wineland et al. (1987); Li et al. (2017); Noguchi et al. (2014). Such manmade systems can be viewed as artificial atoms and molecules, and they offer unprecedented opportunities for generating and observing novel and exotic many-body states and phenomena, as well as for testing fundamental aspects of quantum physics that are beyond the reach of the natural chemical and condensed-matter systems. These nanosystems include two-dimensional semiconductor Hanson et al. (2007); Yannouleas and Landman (2007) and graphene Romanovsky et al. (2009) quantum dots confining electrons and ultracold traps confining neutral atoms Serwane et al. (2011); Zürn et al. (2012); Murmann et al. (2015); Ramanathan et al. (2011) or ions Wineland et al. (1987); Li et al. (2017); Noguchi et al. (2014) in a variety of trap shapes. Among the rich physics studied in these systems, one can mention Wigner molecules (which extend Wigner crystals to the quantum regime), the connection to the fractional quantum Hall effect for high magnetic fields, Aharonov-Bohm phenomena and quantum space-time crystals in ring-shaped devices, wave function entanglement111A pure quantum state describing two or more particles is entangled if it is unfactorizable. A mixed state is entangled if it cannot be written as a mixture of factorizable pure states Wootters (1998); see further Eckert et al. (2002); Aspect (2004). , Schrödinger-cat-state superpositions in strings of ultracold ions in linear traps, and the elucidation of the nature of correlations in assemblies of strongly repelling electrons (long-range Coulomb interaction) or strongly interacting neutral atoms with both an attractive or repulsive contact interaction. Areas of potential applications include quantum information999For increasing interest in the intersection between nuclear physics and quantum information science, see also Cloët (2019) and computing, improved electronic and photonic devices, atomic clocks, metrology, etc.

This section provides an outline of the attainments of symmetry-restoration methods in the area of mesoscopic systems beyond atomic nuclei. For an extensive background exposition, we note an earlier review in Yannouleas and Landman (2007).

7.1 The general hierarchical methodology.

Symmetry restoration in electronic and atomic systems is a particular step in a multilevel hierarchical scheme, which produces a lower total energy at each step; see figure 19 that describes the successive levels in this hierarchy. One starts with the restricted Hartree Fock (RHF) whose wave function [a single Slater determinant for fermions or permanent for bosons Romanovsky et al. (2004, 2006))] preserves all the symmetries of the many-body Hamiltonian; in particular, it imposes the spatial symmetries of the many-Hamiltonian on each individual HF orbital. The next level is an unrestricted Hartree Fock (UHF) whose single determinant allows for the breaking of some (or all) of the Hamiltonian symmetries in an appropriate range of the two-body interaction. In this case, the UHF total energy is lower compared to the RHF one, but the UHF space orbitals do not reflect the space symmetries of the many-body Hamiltonian, a behavior that is often referred to as "Löwdin’s symmetry dilemma" Lykos and Pratt (1963).

In a subsequent step, the broken symmetry in the UHF solutions is restored and the symmetry dilemma is resolved. This level produces a multi-determinantal (or multi-permanent) wave function by applying the projection-operator technique on the UHF single determinant. This level, which is depicted as a single item in figure 19, consists of two substeps, namely the step of variation before projection and the step of variation after projection, see section 3. The variation-after-projection step produces lower energies in general, while the wave functions retain the same multi-determinantal structure as in the variation-before-projection step. The energy difference between these two projection variants decreases as the symmetry breaking becomes stronger. In the context of this section, an example of the variation-after-projection step is offered by Romanovsky et al. (2004, 2006) where the localized-particle orbitals [displaced Gaussians of equation (175) with variational parameters] were used to build an approximate UHF Slater determinant for fermions (or permanent for bosons).

The final level corresponds to a configuration-interaction treatment which in principle provides the exact many-body energies and wave functions. The RHF and UHF are mean-field approximations; the restoration of symmetry and the configuration interaction are often referred to as beyond-mean-field approaches.

We note here that the symmetry-restoration approaches discussed earlier within the context of nuclear physics imply a similar hierarchy as the one shown in figure 19.

7.2 Quantum dots

Advances in nanolithography and growth techniques enabled the fabrication of small semiconductor devices with dimensions in the nanoscale range; they are known in the literature as quantum dots and they play a central role in the modern field of nanotechnology. Here we focus on two-dimensional electrostatically controlled quantum dots Kouwenhoven et al. (1997). Quantum dots are often referred to as "artificial atoms" Kastner (1993); Kouwenhoven and Marcus (1998) due to their having a discrete single-particle spectrum arising from their finite size. Such a terminology invokes a 2D analogue of the physics of 3D electronic shells (whether closed or open) which is associated with the Mendeleev periodic table of natural elements Kouwenhoven and Marcus (1998).

However, it was rather early realized Yannouleas and Landman (1999, 2000c, 2001, 2002b, 2002a), through UHF calculations that in 2D quantum dots, the process of symmetry breaking is highly operative [see also Müller and Koonin (1996)]. This is unlike the case of natural atoms, where due to the overwhelming Coulomb attraction from the central nucleus, the extent of spherical-symmetry breaking is minimal Fertig and Kohn (2000). As a result, the physics of 2D quantum dots overlaps in several ways Yannouleas and Landman (2007) with the nuclear many-body problem, transposed however in the milli-eV (meV) energy range, instead of the mega-eV (MeV) range of atomic nuclei.

7.2.1 The microscopic many-body Hamiltonian.

Before proceeding with the description of the many-body Hamiltonian, a brief discussion concerning the spin-orbit coupling is informative. As usual in atomic and molecular physics, the $L$-$S$ coupling scheme, where $L$ stands for orbital angular momentum and $S$ for spin, is also used for condensed-matter nanosystems and trapped ultracold atoms. This allows that the restorations of the orbital angular momentum and spin can be carried out independently of each other, whereas the restoration of the combined (orbital and spin) total angular momentum ${\bm{J}}={\bm{L}}+{\bm{S}}$ ($J$-$J$ coupling) is pertinent in nuclei due to the strong spin-orbit interaction. In contrast, the spin-orbit coupling in atomic, molecular, and electronic systems considered here is often weak compared to the corresponding coupling in nuclei.

The spin-orbit in condensed-matter systems (like quantum dots) is treated in two varieties: (i) the Rashba type Bychkov and Rashba (1984) and (ii) the Dresselhaus type Dresselhaus (1955). The Rashba or Dresselhaus couplings can be included following the steps of restoration of the spin and angular-momentum broken symmetries. For an example of incorporating the Rashba and Dresselhaus spin-orbit couplings in the context of two-dimensional (2D) quantum dots, see the configuration-interaction calculations in Szafran et al. (2009).

The many-body Hamiltonian describing $N$ fermions or bosons interacting via a two-body potential $U({\bm{r}}_{i}-{\bm{r}}_{j})$ is given by

[TABLE]

where $\hat{H}_{\mbox{\scriptsize{sp}}}({\bm{r}}_{i},{\bm{p}}_{i})$ denotes the single-particle Hamiltonian, which depends on the position ${\bm{r}}_{i}$ and momentum ${\bm{p}}_{i}$ of the $i$th particle. For electrons and ultracold ions, the Coulomb repulsion is pertinent as the two-body interaction in equation (169), namely

[TABLE]

where $e$ is the elementary charge (we assume single-ionized ions) and $\kappa$ is the dielectric constant of the material in the case of semiconductor quantum dots; for trapped ultracold ions, $\kappa=1$. Besides the familiar electrons, examples of trapped ultracold ions are: Be+, Ca+, and Yb+.

We assume that, in practice, the system is two dimensional and thus, for ${\bm{r}}=(x,y,z)$, it is confined to the $z=0$ plane. The single-particle Hamiltonian in a perpendicular external field ${\bm{B}}=(0,0,B)$ is given by

[TABLE]

where $m$ and ${\bm{p}}$ denote mass and momentum of the particle, respectively, the external-confinement potential acting in the $z=0$ plane is denoted by $V(x,y)$, and in the symmetric gauge, the vector potential ${\bm{A}}({\bm{r}})$ is given by

[TABLE]

In the case of charged particles, ${\bm{B}}$ coincides with the natural magnetic field and $\eta=e/c$. For electrons confined within a quantum dot, the mass $m$ should be replaced by the effective mass $m^{*}$.

Unlike in nuclear physics, in finite two-dimensional systems like quantum dots, the magnetic field plays an important role because of their relatively large spatial size. This allows the full range of orbital magnetic effects to be explored for magnetic fields that are readily attained in the laboratory (less than 40 T). In contrast, for natural atoms and molecules, magnetic fields of extremely large strength (i.e., larger than $10^{5}$ T) are needed to produce novel phenomena related to orbital magnetism (beyond the perturbative regime). Such strong fields are known to occur only in astrophysical environments (e.g., on the surface of neutron stars) Ruder et al. (1994). A main orbital effect is the progressive spatial shrinking of the single-particle orbitals as the magnetic field increases; this behavior can be directly visualized from the analytic width $\lambda$ [equation (176)] of the displaced Gaussian wave function given in equation (175). Another orbital effect is the acquisition of a Peierls phase factor [see again equation (175)]. These orbital effects are prerequisites behind the appearance of celebrated magnetic-field-dependent phenomena, like the Aharonov-Bohm effect and the formation of quantized Landau levels supporting integer and fractional quantum Hall effects.

To model a single circular or elliptic quantum dot, or ultracold confining trap, or a molecule-like double well, the external-confinement potential $V(x,y)$ can assume various parametrizations. In the case of an elliptic confinement in a harmonic-oscillator potential, one has

[TABLE]

where $\omega_{x}$ ($\omega_{y}$) is the oscillator frequency in the $x$ ($y$) direction. When $\omega_{x}=\omega_{y}=\omega_{0}$, the elliptic confinement reduces to the circular (parabolic) one. The appropriate parametrization of $V(x,y)$ in the case of a double potential well is more complicated. Often a parametrization based on a 2D version of a two-center oscillator with a smooth necking is used. Details of the double-potential-well parametrization are described in Yannouleas and Landman (2002b); Li et al. (2009).

7.2.2 Mean-field equations for electrons: UHF wave functions using the

Pople-Nesbet equations or Slater determinants with displaced Gaussians.

The UHF many-body wave function for $N$ fermions is a single Slater determinant,

[TABLE]

where $\chi_{k}({\bm{x}})$ stands for one of the $k=1,\ldots,{N}$ spin orbitals, with the symbol ${\bm{x}}_{j}$ denoting both the space and spin coordinates of the $j$th particle.

We stress here that the spin orbitals are characterized by conserved values of the projection of spin on the $z$ axis. Therefore, they are expressed through spatial orbitals, $\varphi^{\alpha}_{n}({\bm{r}})$ and $\varphi^{\beta}_{\bar{n}}({\bm{r}})$, and spinors, $\alpha$ (up) and $\beta$ (down), as $\chi_{k=n}({\bm{x}})=\varphi^{\alpha}_{n}({\bm{r}})\alpha,n=1,\ldots,N^{\alpha}$ for spin-up fermions and $\chi_{k=\bar{n}+N^{\alpha}}({\bm{x}})=\varphi^{\beta}_{\bar{n}}({\bm{r}})\beta,\bar{n}=1,\ldots,N^{\beta}$ for spin-down fermions. Thus the UHF Slater determinants (174) are eigenstates of the projection $\hat{S}_{z}$ of the total spin on the $z$ axis with eigenvalue $S_{z}=(N^{\alpha}-N^{\beta})/2$, where $N^{\alpha}$ ($N^{\beta}$) denotes the number of spin up (down) fermions. However, except for the fully spin-polarized case of only spin-up ($N^{\alpha}=N$) or only spin-down ($N^{\beta}=N$) orbitals being occupied, these determinants are not eigenstates of the square of the total spin, $\hat{\bm{S}}^{2}$.

To specify the spin orbitals entering in the UHF Slater determinant (174), one usually solves the self-consistent Pople-Nesbet equations, which are described in Chapter 3.8 of Szabo and Ostlund (1989); see also Yannouleas and Landman (2007). To derive them, one minimizes the total energy $\langle\Phi_{\rm UHF}(1,\ldots,N)|\hat{{\cal H}}|\Phi_{\rm UHF}(1,\ldots,N)\rangle$ by varying the two sets of spatial orbitals $\{\varphi^{\alpha}_{n}({\bm{r}})\}$ and $\{\varphi^{\beta}_{\bar{n}}({\bm{r}})\}$ under the constraint that both sets consist of orthonormal functions. Because these two sets of spatial orbitals are allowed to be different, the Pople-Nesbet equations are also referred to as the approach of “different orbitals for different spins.” We note that each UHF spatial orbital (the output of the Pople-Nesbet equations) is allowed to break the rotational symmetry. On the other hand, for $N^{\alpha}=N^{\beta}=N/2$ and for rotational-symmetry-conserving spatial orbitals that are pairwise identical, $\varphi^{\alpha}_{n}({\bm{r}})=\varphi^{\beta}_{\bar{n}}({\bm{r}})\equiv\varphi_{k}({\bm{r}})$, the Slater determinant of equation (174) corresponds to the RHF approximation, see the hierarchy of approximations displayed in figure 19 and Chapter 3.4 in Szabo and Ostlund (1989).

An illustrative example of broken-symmetry UHF solutions is given in figure 20 for the case of $N=9$ electrons in a parabolic quantum dot at $B=0$. In the case of repulsive interactions [but also for high magnetic fields Yannouleas and Landman (2007)], the symmetry breaking results in particle localization and a lowering of the continuous rotational symmetry to a point-group one. The localized humps in the UHF densities in figure 20 result from the tendency of the particles to avoid each other due to their strong mutual repulsion. For high magnetic fields, a similar localization effect is related to the shrinking of the space orbitals, as mentioned above.

Two UHF isomers of localized electrons, in the notation introduced in section 2.2 denoted as (2,7) and (1,8), are displayed in figure 20. Such nested polygonal-ring isomers are denoted in general as $(n_{1},n_{2},\ldots,n_{r})$, with $n_{r}$ being the number of localized electrons in the $r$th ring. They may compete with each other in a similar way to the prolate and oblate nuclear shape deformations.

The localization of individual particles, revealed by using the self-consistent Pople-Nesbet equations Yannouleas and Landman (1999), suggests a convenient and physically transparent approximation for the broken-symmetry UHF mean-field solution. Namely, one can use a Slater determinant $\Phi_{\rm UHF}^{\mbox{\scriptsize{app}}}(1,\ldots,N)$ made out of non-orthogonal spatial orbitals having the form of displaced Gaussian functions localized at positions ${\bm{R}}_{j}$ Li et al. (2006), i.e.,

[TABLE]

where $\lambda$ can be used as a variational parameter. However, for strong magnetic fields, one can fix it as

[TABLE]

where $\omega_{c}=\eta B/m$ is the cyclotron frequency, familiar for electron systems under the influence of a magnetic field. The augmentation of the effective trap frequency from $\omega_{0}$ to $\widetilde{\omega}$ expresses the associated diamagnetic behavior, which is operative due to the large size of the quantum dot (compared to natural atoms). The phase in equation (175) is due to the gauge invariance of magnetic translations Peierls (1933); London (1937); Ditchfield (1974); Pausch and Klopper (2020) and is given by $\psi({\bm{r}},{\bm{R}}_{j};B)=(xY_{j}-yX_{j})/(2l_{B}^{2})$, with $l_{B}=\sqrt{\hbar/\eta B}=\sqrt{\hbar/(m\omega_{c})}$ being the magnetic length. This approximation proved to be accurate and, in addition, it bypasses the numerical effort involved in solving the self-consistent Pople-Nesbet equations.

At zero magnetic field and for both the cases of a contact potential and a Coulomb interaction, the resulting energy gain from symmetry breaking becomes larger for stronger repulsion. Controlling this energy gain (the strength of correlations) is the ratio $R_{\delta}$ (for a contact potential) and $R_{W}$ (for a Coulomb interaction) between the strength of the repulsive potential and the zero-point kinetic energy. Specifically, for a 2D trap, one has

[TABLE]

with $l_{0}=\sqrt{\hbar/(m\omega_{0})}$ being the characteristic harmonic-oscillator length. [The subscript $W$ stands for “Wigner” Wigner (1934).]

7.2.3 Combining spin and total-angular-momentum restorations.

When the fermions are not fully polarized, the symmetry-broken Pople-Nesbet UHF determinantal solutions do have the total spin projection $S_{z}$ as good quantum numbers. However, the total spin $\hat{\bm{S}}=\sum_{i=1}^{N}\hat{\bm{s}}_{i}$ is not preserved. A simple example is the UHF Slater determinant which describes the $S_{z}=0$ ground state of two electrons in a parabolic quantum dot for $R_{W}=2.40$ (and $B=0$). figure 21(a) displays the azimuthally symmetric RHF electron density, which contrasts with the symmetry-broken UHF one displayed in figure 21(b). The associated UHF determinant is given by equation (18) in section 2.2.

The next step in restoring the total spin is described in detail in section 2.2, resulting in a $s=0$, singlet state given by equation (23); it is a superposition of two UHF Slater determinants. We note that the spatial reflection symmetry (parity) is automatically restored along with the spin symmetry Fukutome (1981).

For the exact singlet state of the circular confinement, one can generate approximate projected wave functions with good total angular momentum, exhibiting the required azimuthally-uniform electron densities, by applying the product operator,

[TABLE]

The singlet spin-projection operator $\hat{P}^{s}_{\rm spin}$ (19) produces the anisotropic singlet wave function of equation (23), and then the total-angular-momentum projection operator,

[TABLE]

acts upon this two-determinant wave function to restore the total 2D angular momentum $L$.101010 The description of the general projection-operator formalism regarding symmetry restoration is presented in section 3. See also Tables 1 and 2, where the nuclear-style notation $\hat{J}_{z}$ is used for the 2D angular momentum $\hat{L}$, commonly used in the condensed-matter and atomic-physics literature. Here we repeat the definition in equation (179) for clarity and for the convenience of the reader. Note that ”1D” in Tables 1 and 2 refers to the dimension of the integral, and not to the geometry of the physical system. For a full exposition of the adaptation of the general formalism to the 2D electronic and atomic systems, see Yannouleas and Landman (2007). In equation (179), $\hat{L}=\sum_{i=1}^{N}\hat{l}_{i}$, $i=1,2,\ldots,N$, $\hbar\hat{L}$ is the two-dimensional total angular-momentum operator, and $\gamma$ is the azimuthal angle. This double projection describes all the lowest-energy states [yrast band Yannouleas and Landman (2000a)] with good total angular momentum $L=0,2,4,\ldots$. (The yrast-band states with odd values, $L=1,3,5,\ldots$, are generated via a projection of the fully-polarized UHF state.)

The evolution of the ground-state $(L=0)$ electron densities according to the successive approximations, RHF, UHF, spin projection (SP), and combined spin-and-angular-momentum projection (SAP) is illustrated in figure 21. The exact wave functions for two electrons in a parabolic confinement are available Yannouleas and Landman (2000a), and the corresponding ground-state electron density is plotted in figure 21(e). The successive lowering of the ground-state total energies is also displayed.

7.2.4 More on spin restoration.

The literature of spin restoration in systems other than nuclei has a more complicated history compared to that of 2D angular momentum. Löwdin introduced a spin projection operator through the expression Löwdin (1955b)

[TABLE]

where the index $s^{\prime}$ runs over the quantum numbers associated with the eigenvalues $s^{\prime}(s^{\prime}+1)$ of $\hat{\bm{S}}^{2}$ (in units of $\hbar^{2}$), with $\hat{\bm{S}}$ being the total spin operator. When applying $\hat{\bm{S}}^{2}$ on a Hartree-Fock determinant, one uses:

[TABLE]

where $\varpi_{ij}$ interchanges the spins of electrons $i$ and $j$ provided that they are different; $N_{\alpha}$ and $N_{\beta}$ denote the number of spin-up and spin-down electrons, respectively, while $N=N_{\alpha}+N_{\beta}$.

For a large number of electrons, a computationally practical implementation of Löwdin’s spin projection formalism was recently discussed in Pons Viver (2018).

The operator $\hat{P}_{\rm spin}(S)$ was used Yannouleas and Landman (2002b); De Giovannini et al. (2007) to describe the energy spectra and wave functions for electrons in 2D quantum dots. We note that, for $N\geq 3$, there are multiple spin eigenfunctions for a given value $S$ of the total spin, and this multiplicity is important for obtaining a complete set of excited states (e.g., the yrast band), in addition to the ground state. The spin multiplicities are tabulated in the so-called branching diagram Li et al. (2009); Pauncz (2000); Yannouleas et al. (2016); Salmon (1974). For example, for $N=3$ fermions, there are two spin eigenfunctions with $S=1/2$. For a spin projection $S_{z}=1/2$ and using the notation ${\cal S}(S,S_{z};i)$ (where the index $i$ is employed for the degeneracies), a pair of basis spin eigenfunctions that spans the associated two-dimensional spin space is given by

[TABLE]

[TABLE]

For $N>3$ electrons, the complete set of spin eigenfunctions can be specified by solving a Heisenberg Hamiltonian $H_{H}=\sum_{i,j}J_{ij}\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{j}$ Li et al. (2009); Yannouleas et al. (2016), whose purely-spin solutions are given as a superposition of spin primitives $|\sigma_{1}\sigma_{2}\ldots\sigma_{N}\rangle$, where $\sigma$ stands either for $\alpha$ (spin-up) or $\beta$ (spin-down). Then using the fact that each electron is associated with a localized space orbital, one can generate a corresponding Slater determinant of spin orbitals from each spin primitive, thus generalizing Shi et al. (2007); Dai et al. (2007) the two-electron Heitler-London-type expression of equation (23).

Naturally, the total-spin restoration111111In the context of a generalized UHF with a broader set of unrestrictions, restoration of both the total spin $\hat{\bm{S}}$ and its projection $S_{z}$ (in the case that the mean-field wave functions break both of these symmetries) has also been discussed in Fukutome (1981). can also be performed Yannouleas and Landman (2007); Fukutome (1981); Hashimoto (1982); Igawa (1995) by using the projection-operator formulas that restore the three-dimensional (3D) total angular momenta, see section 3.2.5.

7.2.5 Molecular symmetries of the UHF wave functions and magic angular momenta.

The 2D projected wave functions have good total angular momentum $L$, and the corresponding single-particle density is circular and azimuthally uniform. Thus any association of the projected wave functions to a point-group symmetry (which corresponds to a single-particle density that is not azimuthally uniform) is counterintuitive. Despite this expectation, the projected trial wave functions do embody and reflect hidden (or emergent) molecular point-group symmetries similar to the case of natural molecules. Specifically, the $C_{N}$ point-group symmetry of the “classical” crystalline configuration, which is accounted by the symmetry-broken mean-field determinants $\Phi_{\rm UHF}(1,\ldots,N)$, equation (174), is reflected in the fact that the symmetry-restored wave functions $\Psi^{\rm PRJ}_{L}$ are identically zero except for a subset of magic angular momenta $L_{m}$. Of course this vanishing of wave functions is not present for the exact configuration-interaction ones with non-magic angular momenta. In the context of the symmetry-restoration approach, correspondence with the full set of configuration-interaction wave functions is established by considering the vibrations of the molecular configurations. This additional step was described in Yannouleas and Landman (2010, 2011); see also Yannouleas and Landman (2000a). An analogous situation in nuclear theory is provided by the generator-coordinate method, see section 3.1.2, which is a generalization of the strict symmetry-restoration approach via projection techniques.

For the simpler case of $N$ repelling fermions with parallel spins on a ring [i.e., with $S_{z}=N/2$ and a $(0,N)$ configuration], the magic total angular momenta can be determined by considering the point-group symmetry operator $\hat{R}(2\pi/N)\equiv\exp(-i2\pi{\hat{L}}/N)$ that rotates on the ring simultaneously the localized particles by an angle $2\pi/N$. In connection to the state $\Psi^{\rm PRJ}_{L}$, the operator $\hat{R}(2\pi/N)$ can be invoked in two different ways, namely, either by applying it on the “intrinsic” (symmetry-broken) wave function $\Phi_{\rm UHF}(1,\ldots,N)$ or on the phase factor $\exp(i\gamma L)$ appearing in the “laboratory” (symmetry-restored) wave function $\Psi^{\rm PRJ}_{L}$ of equation (179), see Chapter 4-2c in Bohr and Mottelson (1998). One gets

[TABLE]

from the first alternative and

[TABLE]

from the second alternative. The $(-1)^{N-1}$ factor in equation (184) results from the fact that the $2\pi/N$ rotation is equivalent to exchanging $N-1$ rows in the determinant $\Phi_{\rm UHF}(1,\ldots,N)$. Now, if $\Psi^{\rm PRJ}_{L}\neq 0$, the only way that equations (184) and (185) can be simultaneously valid is if the condition $\exp(2\pi Li/N)=(-1)^{N-1}$ is fulfilled. This leads to the following sequence of magic angular momenta,

[TABLE]

for $N$ odd, and

[TABLE]

for $N$ even.

The physics associated with magic-angular-momentum yrast states was extensively explored in the literature of 2D quantum dots Yannouleas and Landman (2007, 2006a, 2003b); Ruan et al. (1995); Maksym (1996); Seki et al. (1996); Maksym et al. (2000). An important property is the enhanced energy stabilization (compared to the rest of the spectrum as described by configuration-interaction calculations) that they acquire in their neighborhood in the regime of strong interactions (i.e., for large $R_{W}$, $R_{\delta}$, see equation (177), or for large magnetic fields).

In the case of strong magnetic fields, when the many-body Hilbert space can be restricted to the lowest Landau level, it is customary to calculate the configuration-interaction energy spectra keeping only the interaction term in the many-body Hamiltonian. In such partial spectra, the magic angular momenta are associated with the so-called “cusp” states Yannouleas and Landman (2003b, 2004a); Jain (2007) [see figure 22(a)], which are precursors of the fractional quantum-Hall-effect bulk states Yannouleas and Landman (2003b, 2004a); Jain (2007). When the single-particle part of the Hamiltonian is also included (i.e., kinetic energy plus external confinement), one obtains global energy spectra as a function of the magnetic field $B$; see figure 22(b). In these global spectra, all the ground states correspond to cusp states, illustrating the role played by the magic angular momenta in enhancing energy stabilization.

For magnetic-field-free systems, this energy stabilization leads to a separation of energy scales between the rotational the vibrational motions (formation of a near-rigid rotor), which is a familiar prerequisite in the formulation of nuclear effective field theories Papenbrock and Weidenmüller (2015). An example of such a separation of energy scales is portrayed in figure 1 of Yannouleas and Landman (2000a) for the case of two electrons in a parabolic 2D quantum dot.

In the above derivation, we considered fully polarized fermions only, that is cases when $S=S_{z}=N/2$, where $S$ is the total spin and $S_{z}$ is its projection. Consideration with this methodology of the other spin values $S_{z}<N/2$ is straightforward; it requires, however, restoration of both the total spin ${\bm{S}}^{2}$ and the total angular momentum. An explicit example for $N=3$ fermions is discussed in Yannouleas and Landman (2003a).

7.2.6 A tour of the literature and Wigner molecules.

Naturally, there are several key differences between the physics of 2D quantum dots and that of atomic nuclei, which arise from the fact that the inter-particle interaction in quantum dots is repulsive, instead of attractive as in nuclei, and that quantum dots consist of one kind of fermions (electrons), instead of two kind of particles (protons and neutrons). As a result, symmetry breaking in quantum dots is associated with individual-electron localization in space in the intrinsic frame (leading to formation of mean-field crystalline configurations), rather than the familiar shape deformations of the nuclear central mean-field confining potential.

Such mean-field crystalline configurations in quantum dots [see, e.g., figure 20] are referred to as "Wigner molecules" Yannouleas and Landman (1999). After restoration of the angular-momentum symmetries, they are often referred to Yannouleas and Landman (2004b) as "rotating Wigner molecules"121212The rotating-Wigner-molecule wave functions are stationary, exhibiting a time-independent single-particle density. The term ”rotating” here refers to these wave functions having good quantal total angular momenta, unlike a symmetry-broken crystalline UHF wave function. Time-evolution of wave packets formed through the superposition of several rotating Wigner molecules Yannouleas and Landman (2017) reintroduces (beyond mean field) symmetry breaking and leads to the concept of a quantum space-time crystal Wilczek (2012); Li, Gong, Yin, Quan, Yin, Zhang, Duan and Zhang (2012); Yannouleas and Landman (2017) and to phenomena of quantum-mechanical revival Seideman (1999); Yannouleas and Landman (2017); Eriksson et al. (2018); see also Liu et al. (2018, 2019) for analogous superpositions that break and restore symmetries, experimentally achieved in natural molecules and atoms., exhibiting two special cases of "rotating electron molecules" Yannouleas and Landman (2003b) or "rotating boson molecules" Romanovsky et al. (2006). Localized corpuscular patterns or cluster structures arise also in symmetry-broken mean-field single-particle densities of lighter nuclei Ichikawa et al. (2012); Girod and Schuck (2013); Ebran et al. (2013, 2017); they are, however, associated with $\alpha$-particle multi-nucleon clustering.

Starting with the early 2000’s, the two-step methodology that combines symmetry breaking with subsequent symmetry restoration was employed extensively to investigate the physics of quantum dots. In particular, using the Löwdin projection for restoring the total spin, Yannouleas and Landman (2001) investigated the coupling and dissociation of two electrons in a double-well confinement (artificial H2 molecule), see section 2.2. In addition to the spin restoration, the formation of a two-electron rotating Wigner molecule was described in Yannouleas and Landman (2002b) by restoring simultaneously the total angular momentum in the case of two electrons confined in a parabolic (circularly symmetric) single-well quantum dot.

For the case of zero or low magnetic fields and using the two-step method, subsequent literature studied a larger number of electrons in parabolic quantum dots (in the range of $3\leq N\leq 10$). An explicit demonstration that the projected ground-state wave function has a lower energy compared to the UHF one was given in Mikhailov and Ziegler (2002) for $N=2-8$ fully spin-polarized electrons. For $N=3$ electrons, a detailed analysis of the lower point-group symmetries of the UHF broken-symmetry molecular solutions and their influence upon the angular-momentum-restored wave functions was also carried out Yannouleas and Landman (2003a). The richness of the physics embodied in the projected wave functions was illustrated in Yannouleas and Landman (2004b), where it was shown that the rotating Wigner molecule can attain two opposite limits depending on the parameters of the system. Namely the limit of a rigid 2D rotor is reached for strong Coulomb repulsion (e.g., $R_{W}=200>>1$) in the absence of an applied magnetic field; the rotational spectrum (yrast band) in this case exhibits energy levels $\propto L^{2}$. An opposite limit of a hyper floppy rotor is reached for smaller $R_{W}\sim 10$, but very high magnetic field (the lowest-Landau-level regime); in this case the rotational energies (yrast band) have an ${\cal A}L+{\cal B}/\sqrt{L}$ dependence on the total angular momentum $L$. The limit of a 2D rigid rotor for $R_{W}\rightarrow\infty$ and low magnetic field was also demonstrated for the case of $N=9$ and $N=8$ ultracold ions confined in a 2D ring-shaped trap Yannouleas and Landman (2017). The limit of the 2D rigid-rotor rotational spectrum extracted in the papers above is reminiscent of the projected-energy Kamlah expansion131313The Kamlah expansion needs to be used in conjunction with the $C_{N}$ Wigner-molecule lower symmetry; otherwise Müller and Koonin (1996) the multifaceted effects originating from the magic angular momenta are missed. Moreover in the lowest-Landau-level regime, use Müller and Koonin (1996) of the Kamlah expansion cannot reproduce the $1/\sqrt{L}$ energy component characteristic of the hyper-soft rotor Yannouleas and Landman (2004b); Li et al. (2006). in integer powers of $L$ for strong symmetry breaking in rotating nuclei [see section 4 and Chapter 11.4.4 of Ring and Schuck (1980)]; in the present cases, however, only the dominant term $\propto L^{2}$ survives for $R_{W}\rightarrow\infty$.

Using broken-symmetry UHF solutions and following the Löwdin prescription for the total-spin projection [see equation (180)] in connection with the construction of spin eigenfunctions presented in Smith Jr. (1964), the combined restoration of both total-spin and angular-momentum approach was applied in a systematic investigation De Giovannini et al. (2007, 2008) at zero and low magnetic field $B$ of the properties of 2D parabolic quantum dots with up to $N=12$ electrons. In particular for $B=0$, it was confirmed that Hund’s rules apply for weaker interaction with $R_{W}\leq 2$; for stronger interaction ($R_{W}>4$), Hund’s rules are violated signaling the dominance of a strong Wigner molecule Yannouleas and Landman (1999).

For completeness, we mention that collective modes associated with the spurious RPA states were used to restore the broken rotational symmetry of UHF solutions in parabolic quantum dots. The case of $N=2$ electrons was systematically studied Serra et al. (2003); Birman et al. (2013). This RPA-based approach, however, becomes computationally prohibitive for larger $N$, due to the increasing number of RPA modes that are required.

The symmetry-restoration methodology was also successfully used to describe aspects of the many-body physics of few electrons in the lowest Landau level. This level forms at very large magnetic fields $B\rightarrow\infty$, and it consists exclusively of all single-particle levels $\propto r^{\ell}e^{i\ell\phi}e^{-r^{2}/2\lambda_{c}^{2}}$ with zero radial nodes and arbitrary single-particle angular momentum $\ell$. These levels are degenerate with energy of $\hbar\omega_{c}/2$, where $\omega_{c}$ is the cyclotron frequency $\omega_{c}=eB/m^{*}c$, and $\lambda_{c}=l_{B}\sqrt{2}=\sqrt{2\hbar/(m^{*}\omega_{c})}$; see, e.g., the Appendix in Yannouleas and Landman (2007).

By constructing a Slater determinant out of the displaced Gaussian orbitals in equation (175) (with $\lambda=\lambda_{c}$), and projecting out the good total angular momentum $L$, one can derive analytic expressions for the rotating electron molecule Yannouleas and Landman (2002c) for any number $N$ of fully spin-polarized electrons (i.e., with $S=S_{z}=N/2$) and any $L$. Note that for large magnetic fields, the electrons in the ground state are fully spin-polarized. Analytic expressions were derived for both the cases of rotating electron molecules with $(0,N)$ Yannouleas and Landman (2002c) and $(1,N-1)$ Yannouleas and Landman (2003b) ring configurations. A numerical investigation of lowest-Landau-level rotating electron molecules exhibiting a configuration where the electrons are arranged in a configuration consisting of $r$ concentric regular polygons ($n_{1},n_{2},\ldots,n_{r}$, $N=\sum_{i}^{r}n_{i}$) was also presented Yannouleas and Landman (2004a). Corresponding analytic expressions for rotating bosonic molecules for $N$ spinless bosons in the lowest Landau level in a double-ring configuration, $(n_{1},n_{2})$ with $n_{1}+n_{2}=N$, were subsequently derived Yannouleas and Landman (2010).

Going beyond the rotating-electron or rotating-boson molecular states (which describe pure vibrationless rotations), a class of trial wave functions portraying combined rotations and vibrations of Wigner molecules associated with concentric polygonal rings was further introduced Yannouleas and Landman (2010, 2011). These trial functions, referred to as rovibrational molecular functions, are valid for both bosons and fermions and provide a correlated basis that spans the translationally invariant part of the lowest-Landau-level spectra for both the yrast and excited lowest-Landau-level states, and for both low and high angular momenta Yannouleas and Landman (2010). As a result, the restoration of broken symmetry approach can describe the totality of the lowest-Landau-level states and not only the cusp states which are associated with ground states exhibiting magic angular momenta that are precursors of the fractional quantum Hall effect states characterized by fractions $\nu$.

A major subject in the lowest-Landau-level physics was the emergence of actual broken-symmetry Wigner-solid crystal states. Such Wigner-solid crystals were expected to appear for smaller fractions $\nu\leq 1/5$. It was thus surprising that a Wigner-crystal regime was experimentally observed Zhu et al. (2010) in the neighborhood of $\nu=1/3$ in the case of very clean samples. An interpretation of these observations was achieved using linear superpositions (wave packets) of angular-momentum-restored wave functions (specifically the analytic ones of the rotating electron molecules). These superpositions involve summation over several cusp states with different magic angular momenta; they naturally break the rotational symmetry to exhibit explicitly the crystalline structure, without necessarily reverting back to the UHF level. The triggering agent for the pinning of the rotating Wigner molecule and the enforcing of symmetry breaking is the presence of residual impurities and disorder in the sample.

For non-fully spin polarized electrons, the symmetry restoration in the lowest Landau level must involve both the total spin $S$ and the angular momentum $L$. Such combined $S$ and $L$ projection leading to spin-dependent rotating electron molecules with $S<N/2$ was performed for $N=4-5$ localized electrons in the lowest Landau level Shi et al. (2007); Dai et al. (2007). The combined spin and space projection was also demonstrated for $N$ lowest-Landau-level electrons confined in a ring geometry Yang et al. (2008).

Of interest is the property that the edge states at zero-magnetic field in a circular graphene dot with a zig-zag termination form a collection mimicking the lowest-Landau-level manifold; these edge states appear due to the existence of two valleys in the single-particle spectrum of the zero-mass Weyl-Dirac graphene electron. The formation of rotating Wigner molecules in this novel lowest-Landau-level manifold was investigated using both configuration-interaction and projection techniques Wunsch et al. (2008); Romanovsky et al. (2009).

7.3 Trapped ultracold neutral atomic gases and ions

Ultracold trapped neutral atoms interact via a Dirac-delta contact potential, namely, in equation (169) one takes

[TABLE]

In equation (188) above, the strength parameter $g$ can take both negative (attractive interaction) and positive (repulsive interaction) values. Experimentally, this parameter can be varied continuously from the attractive to the repulsive regime; see, e.g., Zürn et al. (2012); Brandt et al. (2015). For ultracold neutral atoms, the magnetic field ${\bm{B}}$ in equation (169) can be mimicked with artificial synthetic fields Goldman et al. (2014) or the rotational frequency $\Omega$ of a rotating harmonic trap Romanovsky et al. (2006) [with an appropriate modification of the parameter $\eta$ in equation (169)]. 6Li (fermionic) and 87Rb (bosonic) are examples of trapped neutral atoms.

In the case of ultracold atoms and molecules, dipole-dipole two-body interactions were also experimentally realized. However, no symmetry-restoration investigations with dipolar interactions were reported as yet.

The restoration of angular momentum was employed Romanovsky et al. (2004) to investigate systems with a finite number $N$ of spinless neutral and charged bosons (ions) confined in a 2D harmonic trap. The broken-symmetry UHF-type orbitals were approximated as in equation (175), treating the positions and the widths of the displaced Gaussians as variational parameters (which corresponds to a variation-after-projection scheme). Wigner molecules were described for both neutral and charged bosons in the regime of strong interparticle repulsion. For the case of neutral bosons, the Wigner-molecule regime corresponds to a process of 2D fermionization, when the strong repulsion keeps the particles away from each other overtaking the propensity of bosons to bunch together due to statistics. This fermionization behavior is well known for strongly repelling strictly 1D bosons Girardeau (1960); in two dimensions, it has also been recently further verified via exact numerical calculations for two interacting bosons Mujal et al. (2018).

The system of $N$ spinless bosons in rotating harmonic and toroidal traps was also studied Romanovsky et al. (2006) using angular-momentum restoration techniques. Figure 23 illustrates the patterns in the single-particle and two-body conditional probabilities as the successive steps sketched in figure 19 are applied. The conditional probability is defined [see section 1.5 in Yannouleas and Landman (2007)] as the probability for finding a particle at position ${\bm{r}}$ given that another particle is located (fixed) at a reference point ${\bm{r}}_{f}$. To be noted is the fact, see Romanovsky et al. (2006), that the ground state has zero angular momentum only for small values of the rotational frequency $\Omega$ (or equivalently for small magnetic field $B$); for larger values of $\Omega$, the ground-state angular momentum increases in steps of $N$, which is the hallmark of the emergence of magic angular momenta.

An interesting application Romanovsky et al. (2008) of the methodology of projection techniques is the restoration of the rotational symmetry starting from broken-symmetry Gross-Pitaevskii solutions that describe vortices in rotating traps. In this case, the symmetry-restored many-body wave functions can be characterized as rotating vortex clusters. The presence of vortices is not visible in the single-particle densities of the rotating vortex clusters, which are homogeneous; to reveal the hidden vortices, one needs to employ the conditional probability distributions. The Gross-Pitaevskii vortex states are shown to be wave packets composed of such rotating-vortex-cluster states.

Finally, we mention that the restoration of angular-momentum (under the name of "continuous configuration-interaction") was applied Alon et al. (2004) in the case of an attractive Bose gas on a ring. For strong attraction, the full-configuration-interaction method converges very slowly with the increasing dimension of the employed discrete single-particle basis, and thus the restoration of broken symmetry is advantageous, yielding lower total energies compared to the nonconverged configuration-interaction ones.

7.4 Spin-projected UHF, Hubbard systems, and natural molecules

In addition to the description of novel strongly-correlated many-body phases and phenomena Yannouleas and Landman (2007, 2006a) for electrons in quantum dots and trapped ultracold neutral atoms and ions, the two-step method of symmetry breaking and symmetry restoration beyond nuclear physics is also being developed in the direction of a powerful computational approach that can rival in effectiveness the Kohn-Sham density-functional computational technique. This computational direction is mainly associated with the spin-projected UHF (restoration only of spin) in the context of condensed-matter Hubbard systems and natural molecules traditionally belonging to the field of chemistry. The specific approach used to restore the total spin is similar to that used for the 3D total angular-momentum projection in nuclear physics Ring and Schuck (1980). This has definite practical advantages Tsuchimochi and Van Voorhis (2015) for large-scale computations compared to the prescription suggested by Löwdin Löwdin (1955b); however, see Pons Viver (2018). In this respect, we mention that, unlike the symmetry-restoration techniques, the single-determinantal Kohn-Sham density-functional formalism fails to provide a complete description of magnetic states Kaplan (2018), as well as to describe properly the regime of strong static correlation [when there are several competing degenerate states, as in dissociation processes and the formation of Wigner molecules; see, e.g., Cohen et al. (2008); Yannouleas and Landman (2007, 2006a)].

The Hubbard model Hubbard (1963) is widely used to describe strongly interacting electrons in periodic potentials (associated with natural ionic lattices in condensed-matter systems) and most recently ultracold atoms trapped in artificial optical lattices. The Hubbard model Hamiltonian for $N$ fermions is given by

[TABLE]

where $<i,j>$ denotes summation over nearest-neighbor sites and $\sigma$ sums over the up ($\uparrow$) and down ($\downarrow$) spins, with $n_{i\sigma}=a_{i,\sigma}^{+}a_{i,\sigma}$. The parameters $t$ and $U$ are the hopping parameter and the on-site repulsion, respectively.

In particular, following an early publication Igawa (1995), the method of restoration of spin symmetries was applied with increasing numerical sophistication in the case of one-dimensional Hubbard chains Tomita (2004, 2009); Schmid et al. (2005); Rodríguez-Guzmán et al. (2013). Partially two-dimensional Hubbard geometries (e.g., 2$\times$2 and 2$\times$4 plaquettes) were also studied Rodríguez-Guzmán et al. (2012) using the spin-projected HF approach. Interestingly, it was shown that the exact ground-state in a four-site Hubbard plaquette can be recovered by the spin-projected wave function irrespective of the interaction strength Leprévost et al. (2014). A combination of spin projection with Gutzwiller-type double-occupancy screening was also applied to 2D Hubbard lattices Wahlen-Strothman and Scuseria (2016).

The restoration of the total-spin symmetry was employed further to describe the ground-state correlations and dissociation profiles of natural molecules Scuseria et al. (2011); Jiménez-Hoyos et al. (2012, 2013) familiar from quantum chemistry. Examples of investigated systems were O2, N2, H2O, [Cu2O2]2+ core, etc. A projected coupled cluster theory is being also developed for both natural molecules and the Hubbard model. It was shown Qiu et al. (2017) that spin restoration via projection techniques significantly improves unrestricted coupled-cluster-method results while reimposing the required good quantum numbers.

7.5 Other electronic systems

In the early 1980’s, it was discovered that the doubly-excited states of the Helium atom exhibit rovibrational spectra that reflect the formation of a highly nonrigid linear symmetric XYX "molecular" structure Kellman and Herrick (1980); Berry (1989), where the X’s stand for the two excited electrons and Y for the He nucleus. In addition to other methods Berry (1989), these rovibrational spectra and corresponding underlying collective wave functions were studied using the approach of 3D angular-momentum restoration Iwai and Nakamura (1989). Such molecular structures in highly excited atoms indicate that physical processes associated with symmetry breaking cannot be dismissed even in the case of natural atoms; they provide a bridge Yannouleas and Landman (2007, 2006a) to the regime of Wigner molecules in two-dimensional quantum dots.

Another notable application of projection techniques beyond nuclear physics is the use of a number-projection method, see section 8.1 and Rossignoli (1995), to calculate the canonical-ensemble, temperature-dependent free energy of metal clusters for describing temperature effects in electronic-shell energy contributions Frauendorf and Pashkevich (1996), and in particular for describing temperature attenuation in ionization potentials, electron affinities, and fission fragmentation of poly-cationic and poly-anionic clusters Yannouleas and Landman (1997, 2000b); Yannouleas et al. (2002).

Particle-number projection techniques were also used to investigate the properties of superconducting metal grains Fernández and Egido (2003, 2005).

Finally worth mentioning is the use of projection-operator techniques to describe the dynamic Jahn-Teller effect in natural molecules in the case of tunneling between equivalently distorted energy-minimum configurations of the adiabatic potential energy surface Hallam et al. (1992); Dunn et al. (2012). Naturally, due to the very large masses of the ionic cores, the explicit symmetry-broken wave-packet state localized within a single minimum can be observed in cases when tunneling is suppressed Bersuker (2016). This is analogous to the observation of pinned classical Wigner crystals of trapped ultracold ions Yannouleas and Landman (2017); Thompson (2015).

7.6 Other emerging directions

7.6.1 Relation to entanglement and quantum information science.

The emergence of modern quantum information theory is being built around exotic and counterintuitive theoretical concepts, such as entanglement Aspect (2004); Wootters (1998) and quantumness Piani et al. (2014); Modi et al. (2012), which reflect the complexity of the structure of the quantum wave functions (e.g., non-factorizability in the case of two or more particles), or of quantum measurement. The symmetry-breaking mean-field solutions are at a disadvantage in this area because they do not conserve symmetries of the many-body Hamiltonian. In this context, it was shown Zeng et al. (2014) that the broken-symmetry BCS wave function represents a class of wave functions where the required quantumness was lost. It is noticeable that the lost quantumness [in the form of proper description for the concurrence Wootters (1998) and quantum discord Ollivier and Zurek (2001)] is restored simultaneously with the restoration of the particle number symmetry in the projected BCS wave function.

Earlier, the ability of the total-spin, symmetry-restored wave function to describe properly the entanglement [in the form of concurrence and von Neumann entropy Wootters (1998)] for two electrons in a double quantum dot under the influence of an increasing magnetic field was also investigated in detail Yannouleas and Landman (2006b, a, 2007).

7.6.2 Time evolution in finite systems out of equilibrium.

Apart from the small-amplitude harmonic vibrations, broken-symmetry wave functions (single determinants or permanents) fail to describe the proper time-evolution behavior when propagated in time with the corresponding mean-field Hamiltonian Lichtner and Griffin (1976); Yannouleas and Landman (2017); see also Chapter 12.2.4 in Ring and Schuck (1980). This drawback of the mean-field treatment of finite systems was earlier discussed in the framework of heavy-ion collisions in nuclear physics Lichtner and Griffin (1976). It is easily overcome by expressing the broken-symmetry wave function as a wave packet (superposition) of symmetry-restored wave functions and evolving independently in time (by multiplying by a time-dependent phase) each component of the wave packet. Using this approach, other symmetry-broken wave packets (different from the UHF solutions) can be envisaged that exhibit single-particle densities with controlled periodicities in both space and time, as was recently discussed in the framework of implementing a quantum space-time crystal of ultracold atoms or ions in a ring-shaped trap Yannouleas and Landman (2017). If the initial wave packet reproduces the UHF or Gross-Pitaevskii broken-symmetry solution, revival and recurrence in-time behavior is generated Eriksson et al. (2018).

8 Projected statistics

Symmetry restoration can also be introduced in the context of quantal statistical ensembles. The formal treatment of statistical ensembles requires the replacement of mean values of operators by traces over the whole Fock space "weighted" by a density matrix operator $\hat{D}$ responsible for the probability distribution Kadanoff et al. (1994); Huang (1987); Shieve and Horwitz (2009); Attard (2015). The form of this operator depends on the problem at hand but it is typically defined as the exponential of the Hamiltonian plus some additional terms. When the problem is restricted to the mean-field level, the density matrix operator is the exponential of the one-body mean-field Hamiltonian including pairing fields. In the trace, all possible multi-quasiparticle excitations of the mean-field ground state are considered.

In those mean-field applications in nuclear physics that require the inclusion of pairing correlations, the density matrix operator to be used is the one of the grand canonical ensemble. This is required to accommodate the possibility of exchanging particles with the external "reservoir". When restricted to the mean-field approximation at finite temperature $T$, $\hat{D}$ is proportional to the exponential of the one body HFB Hamiltonian $H_{\mathrm{HFB}}$

[TABLE]

with $\beta=1/(k_{B}T)$ and $Z=\Tr[\hat{D}_{\mathrm{HFB}}]$ the partition function. As $\hat{H}_{\mathrm{HFB}}-\lambda\hat{N}$ is a quadratic form of creation and annihilation operators

[TABLE]

its exponential $\hat{D}_{\mathrm{HFB}}$ is the operator of a canonical transformation acting on the quasiparticle operators satisfying

[TABLE]

with $\mathcal{H}^{\prime}=\mathcal{H}-\lambda\mathcal{N}$. The matrices $\mathcal{H}$ and $\mathcal{N}$ have been introduced in equation 74. This identity allows the calculation of any statistical trace by using Gaudin’s theorem for Hartree-Fock states Gaudin (1960), which is the extension to statistical ensembles of Wick’s theorem. Its generalization to HFB states is straightforward Ring et al. (1984); Rossignoli and Ring (1994). Like Wick’s theorem, it allows us to write the trace of any operator times the HFB statistical density matrix as a contraction of the operator’s matrix elements with the density and pairing tensors for the statistical ensemble. Introducing the set of operators $a_{\mu}=(a_{1},\ldots,a_{N},a^{+}_{1},\ldots,a^{+}_{N})$ one gets for the contraction

[TABLE]

The zero temperature limit of this result can be used to derive a sort of generalized Wick’s theorem for multiquasiparticle overlaps Perez-Martin and Robledo (2007). If we denote by

[TABLE]

the matrix diagonalizing $\mathcal{H}^{\prime}$ (see equation 74 ) then

[TABLE]

where $f_{\mu}=1/(1+\exp(\beta E_{\mu}))$ are the Fermi statistical occupation factors given in terms of the quasiparticle energies $E_{\mu}$, eigenvalues of the HFB-equations, see equation 74. Inserting this in equation 192 we finally arrive to the standard definitions for the finite temperature density $\rho=V^{*}(1-f)V^{T}+UfU^{+}$ and $\kappa=UfV^{+}V^{*}(1-f)U^{T}$.

8.1 Symmetry restoration at finite temperature

In the context of finite temperature or statistical ensembles, symmetry restoration means that the density matrix operator $\hat{D}$ has to be able to select from all the states considered in the statistical trace only those with the desired set of quantum numbers. The easiest way to achieve this property is by sandwiching the statistical operator with the projector onto the required quantum numbers one wants to select. For instance, for particle-number restoration we have to replace the density matrix operator $\hat{D}$ by

[TABLE]

For non-abelian symmetry groups, like the one of angular-momentum projection the expression of the density matrix operator gets a bit more involved and the reader is referred to Rossignoli and Ring (1994) for the technical details. In general, the expression for $\hat{D}_{N}$ is rather involved, but it simplifies enormously if $\hat{D}$ is restricted to be the exponential of an one-body operator, as it is the case in the mean-field approximation to the exact $\hat{D}$. In the following and just to illustrate the method, we will restrict the discussion to the abelian case of particle-number projection. The projector $P^{N}$ is a linear combination of exponentials of one body operators. Therefore, $\hat{P}^{N+}\hat{D}_{\mathrm{HFB}}\hat{P}^{N}$ becomes a (very involved) linear combination of products of exponentials of one body operators. As in equation (191) those products are generators of canonical transformations

[TABLE]

where $\mathcal{T}_{i}$ are the matrices defining the one body operators in the exponents of $\hat{T}_{i}$,

[TABLE]

As a consequence of this property, Gaudin’s theorem can still be used just replacing the $\exp[-\beta\mathcal{H}^{\prime}]$ in equation (192) by the appropriate product of exponentials (see Rossignoli and Ring (1994) for details). The only difficulty in carrying out this program is in the evaluation of the entropy, required to evaluate the free energy $F=H-TS$. The standard definition of the entropy involves the logarithm of the density matrix operator. In the standard mean-field approximation, this logarithm can be evaluated analytically and the final expression for the entropy in terms of quasi-particle energies is straightforward. Unfortunately, the projected density matrix can not be expressed in general as the exponential of an operator and therefore the evaluation of the entropy becomes a very complicated task Esebbag and Egido (1993). In spite of these difficulties, the use of projected statistics proved to be advantageous over other techniques when applied in the spirit of projection after variation, that does not require the evaluation of the entropy Fanto et al. (2017). The intrinsic difficulty associated with the sign ambiguity in the evaluation of the partition function was addressed in Fanto et al. (2017) in a time reversal preserving scenario and further generalized using the Pfaffian method Robledo (2009) to the more general case involving time reversal breaking intrinsic states Fanto (2017).

The shell model Monte Carlo method Lang et al. (1993); Koonin et al. (1997a, b) is often used to evaluate the partition function of nuclei with relatively large configuration spaces. The shell model Monte Carlo requires the use of particle-number projection to carry out calculations in the more convenient canonical ensemble Alhassid et al. (1999) and therefore this constitutes another field of application of the techniques discussed here. Recently, the use of particle-number projection to carry out calculations in the canonical ensemble was also explored in Magnus et al. (2017).

8.2 Thermo-field dynamics

To finish this section, let us briefly mention an alternative to the traditional approach described above and known under the name of thermo-field dynamics. It consists in computing the traces of statistical operators by using mean values of pure states defined in an extended Fock space including twice the original degrees of freedom. This approach was introduced in the context of quantum field theory by Takahashi and Umezawa Takahashi and Umezawa (1996). In this approach, the statistical average of an operator $\hat{F}$ with probabilities $p_{n}$

[TABLE]

is replaced by the mean value of $\hat{F}$ with the wave function

[TABLE]

which is a linear combination of wave functions $|n\tilde{n}\rangle=|n\rangle\otimes|\tilde{n}\rangle$ defined in an extended Fock space which is the tensor product of the original Fock space with itself. The ket $|\tilde{n}\rangle$ represents a new set of states with identical characteristics as $|n\rangle$. In the extended Fock space all the operators are defined as $\hat{F}\otimes\mathbb{I}$ where $\mathbb{I}$ is the identity in the space spanned by $|\tilde{n}\rangle$. A new set of creation and annihilation single-particle operators $\tilde{a}^{+}_{k}$ and $\tilde{a}_{k}$ satisfying fermion canonical commutation relations and anti-commuting with all the elements of the original ${a}^{+}_{k}$ and ${a}_{k}$ set is required too. An advantage of the formalism is that $|\Phi\rangle$ can be written as a HFB state, vacuum of a set of quasiparticles defined in terms of the $\tilde{a}^{+}_{k}$ ,$\tilde{a}_{k}$, ${a}^{+}_{k}$ and ${a}_{k}$ by means of an appropriate BCS like transformation. Therefore, we can use verbatim all the zero temperature formalism developed before (including the generalized Wick’s theorem) to restore symmetries but taking into account properly the doubling of the single-particle Fock space. It is not clear, however, if this method represents any advantage over the traditional one due to the doubling of matrix sizes. The procedure is analogous to the construction of statistical ensembles by taking the trace over a subsystem of the whole Hilbert space of a pure state density matrix operator. Applications to nuclear physics in the context of symmetry restoration are given in Tanabe and Nakada (2005) but only formal expressions are developed in the mentioned reference.

9 Summary, conclusions, and perspectives

In recent years, research in mesoscopic many-body systems witnessed a discernable progress with the development of the state-of-the-art models and methods. In particular, ab initio methods and mean-field theory based on effective interaction are now widely used to elucidate rich and fascinating properties of these quantum many-body systems. On the one hand, applications of the ab initio methods are still restricted to lighter systems only, whereas the mean-field approaches can be applied to investigate mesoscopic systems of any size. In particular, in nuclear physics, the density functional theory was employed to investigate ground-state properties of all nuclear species predicted to exist in the Segré chart.

Spontaneous symmetry breaking mechanism, inherent to the mean-field based approaches, played an important role in our understanding of many-body systems. For instance, in rotating nuclei, the breaking of rotational symmetry led to the fundamental concept of deformation in nuclei. Nevertheless, the quantal fluctuations of the observables, absent in mean-field approaches, become quite important for mesoscopic systems.

To build quantal fluctuations on top of the mean-field solutions, several approaches were developed. A powerful method to include these fluctuations is through the restoration of the broken symmetries. In the prelude section of this review, the spontaneous symmetry breaking mechanism was illustrated through three simple examples. The main purpose of the present review was to provide an overview of the recent developments and, more importantly, to bring to focus the bottlenecks in the application of symmetry-restoration methods.

The general formalism of symmetry restoration, whose origin can be traced in group theory and generator coordinate method, was laid out in section 3. There we distinguished between the symmetry restoration for abelian groups (particle number, linear momentum, parity), which have mathematical properties of projection operators, versus the symmetry restoration for non-abelian groups, as is the case of the rotational symmetry, relevant for the spatial coordinates, spin, or isotopic spin. In the later case, the symmetry-restoration operators do not obey properties of a projection operator, but they nevertheless project out the relevant symmetry quantum numbers. Further, we discussed the fact that the symmetry restoration can be performed either before minimization or after minimization of the energy functional. In the former approach, commonly referred to as variation after projection, symmetry breaking states (often called “intrinsic states") are determined by application of the variational principle on the projected energy (i.e., energy computed with the projected wave functions). In this procedure, different intrinsic states are obtained for different quantum numbers of the restored symmetry. In the latter approach, referred to as projection after variation, the intrinsic state is determined without consideration of the subsequent projection.

In general, implementation of the symmetry restoration is numerically quite challenging, especially in realistic applications where several quantum numbers need to be restored simultaneously, and calculations need to be performed in the spirit of variation-after-projection method. Due to these numerical challenges, development of approximate projection methods were actively pursued by exploiting the sharp character of the overlap kernel when the intrinsic state strongly breaks the underlying symmetry. This resulted in the development of popular methods of Lipkin-Nogami and Kamlah, and were discussed in section 4. We showed that these methods lead to successful approximations, like the mean-field cranking model – a useful concept to understand the physics of rotational bands. Symmetry-restoration methods were also successfully applied to simple nuclear models, where the Hamiltonian is separable or the configuration spaces are limited to a few oscillator shells. These applications were discussed in section 5.

Although the mechanism of symmetry restoration can be consistently formulated for systems described in terms of a Hamiltonian operator, this is not the case for the energy density functionals, which are commonly employed in nuclear physics to provide a description of low energy observables all over the Segré chart. The complexity of nuclear interaction, and in-medium effects that characterize many-body systems, required the introduction of phenomenological density-dependent interaction terms, for which symmetry-restoration methods cannot be uniquely defined. This is further aggravated when separate interactions are being considered in particle-particle and particle-hole channels. Recently, there were several attempts to overcome these difficulties, but a satisfactory solution, covering both sources of problems, is still not available. A possible solution, which is being vigorously pursued, is to base the functionals on the Hamiltonian picture with explicit three-body terms. However, this approach is not yet sufficiently developed to give definite answers at this stage. Nevertheless, many symmetry restored calculations performed with present-day energy density functionals seem to provide reasonable and a consistent picture of low-energy nuclear phenomena, as was elucidated in section 6. The results, however, should be taken with a pinch of salt as they might be contaminated with spurious effects.

Mesoscopic condensed matter systems and the physics of atoms and molecules, as well as assemblies of trapped ultracold ions and neutral atoms, are mostly free from the above mentioned difficulties as the interaction is often just the Coulomb or a contact interaction between the constituents of the system. Applications of symmetry restoration to those areas share many technical details with the ones in nuclear physics, but there are also clear differences like the fact that rotational symmetry restoration can be carried out separately for the spatial coordinates and the spin in condensed matter physics. Many applications were presented where symmetry breaking and restoration represent an easy way to understand the complexity of the problem, see sections 6 and 7. Finally, in section 8.1 we demonstrated that the concept of symmetry restoration can be extended to the realm of quantum statistical mechanics where pure states are replaced by a set of quantum states with a prescribed probability distribution.

Based on the results presented in this work, it can be concluded that the method of symmetry restoration applied to mean-field wave functions provides a simple and fruitful mechanism to incorporate important dynamic correlations, while still using a simple framework of product wave functions. Furthermore, in this approach one stays within a fully quantum mechanical description from the beginning to end, and the classical picture of collective motion does not have to be invoked.

The generator coordinate method can be employed along with the symmetry restoration to provide a powerful framework to describe quantal fluctuations of relevant degrees of freedom around the mean-field values. In future, we expect development of more advanced configuration interaction approaches with symmetry projected states as the basic building blocks. For instance, projected multi-quasiparticle basis configurations can be constructed around the optimal mean-field, in the spirit of traditional shell-model approach, to describe the physics of excited configurations and also to incorporate many-body correlations in the ground state. Finally, inspired by the successes of symmetry breaking and restoration in mean-field approaches, in nuclear physics and quantum chemistry analogous methods were also recently implemented in the context of the so-called ab initio calculations, see, e.g., Duguet (2014b); Duguet and Signoracci (2016); Qiu et al. (2019); Yao et al. (2020) and section 5.3.

10 Acknowledgments

JD would like to thank Michał Bączyk for a collaboration on the doubly-symmetric-potential-well model. The work of JD was partly supported by the STFC Grants No. ST/M006433/1 and No. ST/P003885/1 and by the Polish National Science Centre under Contract No. 2018/31/B/ST2/02220. The work of LMR was partly supported by Spanish MINECO Grant No. PGC2018-094583-B-I00. PR acknowledges partial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. CY wishes to thank his co-authors at Georgia Tech, and especially Uzi Landman, head of the Center for Computational Materials Science. The work of CY was supported over the years by the Office of Basic Energy Sciences of the US D.O.E. (Grant No. FG05-86ER45234) and by the Air Force Office of Scientific Research (USA) (Grant No. FA9550-15-1-0519).

Appendix A Overlaps and matrix elements between HFB states

The restoration of the symmetries broken by intrinsic HFB states requires to consider matrix elements of various operators between different HFB states. This is a rather general statement because the action of any operator belonging to the symmetry group on a given HFB state is again another HFB state. The origin of this property lies on the fact that the algebra of one-body operators can be used to provide a representation of any matrix Lie algebra Gilmore (2008b). Therefore, the symmetry operations (which are members of the group spanned by the corresponding Lie algebra) are equal to exponentials of one body operators and the Thouless theorem Thouless and Valatin (1962); Mang (1975) applies.

The evaluation of the required matrix elements is best carried out with the help of the generalized Wick’s theorem. The theorem, that can be derived in many different ways Löwdin (1955a); Onishi and Yoshida (1966); Balian and Brezin (1969); Hara and Sun (1995); Bertsch and Robledo (2012); Ring and Schuck (1980), states that the matrix elements of an arbitrary operator $\hat{O}$ between arbitrary, non orthogonal, HFB states $|\Phi_{0}\rangle$ and $|\Phi_{1}\rangle$, $\langle\Phi_{0}|\hat{O}|\Phi_{1}\rangle/\langle\Phi_{0}|\Phi_{1}\rangle$, can always be written in terms of the sum of all possible two-quasiparticle contractions,

[TABLE]

where the $\beta_{\mu}$ and $\beta^{+}_{\mu}$ are annihilation and creation operators associated with $|\Phi_{0}\rangle$. The only non-trivial contraction is given in terms of the skew-symmetry matrix $C_{\mu\nu}$ which is the product of the inverse of $A=U^{+}_{0}U_{1}+V^{+}_{0}V_{1}$ times $B=U^{+}_{0}V_{1}+V^{+}_{0}U_{1}$, i.e., $C=A^{-1}B$. The $U_{0}$, $V_{0}$ and $U_{1}$, $V_{1}$ are the Bogoliubov transformation amplitudes of the corresponding HFB states. As an example, one of the terms entering the Hamiltonian overlaps is given by

[TABLE]

In the general case, where the matrix element of a product of $n$ creation and annihilation quasiparticle operators is required, the number of terms in the sum grows very quickly and is given by $(n-1)!!$. This is the so called combinatorial explosion problem (see Hu et al. (2014) for an example) that hampers applications where multi-quasiparticle excitations have to be considered. This difficulty can be avoided using the formulation of Bertsch and Robledo (2012) in terms of Pfaffians (see below).

The overlap $\langle\Phi_{0}|\Phi_{1}\rangle$ given by the Onishi formula Onishi and Yoshida (1966); Balian and Brezin (1969); Hara and Sun (1995),

[TABLE]

suffers from a sign indeterminacy that requires further consideration. The sign of the overlap affects the quantities to be integrated in the symmetry restoration or configuration-mixing methods. A wrong assignment of the sign even in a small integration interval can substantially change the value of the integral. The sign problem was addressed in the past using different strategies like continuity arguments or determining pairwise degenerate eigenvalues of a general matrix Neergård and Wüst (1983). However, the use of techniques of fermion coherent states allows us to avoid the sign problem by expressing the overlap in terms of the Pfaffian of a skew-symmetric matrix Robledo (2009, 2011); Bertsch and Robledo (2012); Avez and Bender (2012),

[TABLE]

for the skew-symmetric matrices $M^{(i)}=(V_{i}U_{i}^{-1})^{*}$ given in terms of the Bogoliubov amplitudes. Formula (201) assumes that both HFB states are normalized as $\langle 0|\Phi_{i}\rangle=1$ and that their Bogoliubov amplitudes are expressed in a common single-particle basis of dimension $N$, which defines the phase factor $s_{N}=(-1)^{N(N+1)/2}$; see Avez and Bender (2012) for a generalization relaxing this assumption. We also note that recently Mizusaki et al. (2018); Bally and Duguet (2018) have provided additional perspectives into the sign problem. The Pfaffian of a skew-symmetric matrix is a quantity similar to the determinant and shares with it many properties. The numerical evaluation of the Pfaffian can be carried out using the traditional algorithms of linear algebra with a cost similar to the one of the determinant González-Ballestero et al. (2011); Wimmer (2012).