A new and practical formulation for overlaps of Bogoliubov vacua

In this letter we present a new expression for the overlaps of wavefunctions in Hartree-Fock-Bogoliubov based theories. Starting from the Pfaffian formula by Bertsch et al (Phys. Rev. Lett. 108,042505 (2012)), an exact and computationally stable formula for overlaps is derived. We illustrate the convenience of this new formulation with a numerical application in the context of the particle-number projection method. This new formula allows for substantially increased precision and versatility in chemical, atomic, and nuclear physics applications, particularly for methods dealing with superfluidity, symmetry restoration and uses of non-orthogonal many-body basis states.


Introduction.
A very successful approach in the context of many-body theory is to incorporate correlation effects through symmetry breaking followed by restoration of symmetries [2]. The intuitive picture gained from the symmetry-broken solutions has provided much insight into the symmetry and beauty of nature. For instance, many nuclei can be described as intrinsically deformed where both axial and reflection symmetry are broken [3]. Moreover, combining shape degrees of freedom with angular momentum orientation, additional symmetries arise that can also be broken with spectacular manifestations in nuclear spectra [4]. For each broken symmetry new collective modes emerge, which complete a picture of nuclear structure and nuclear spectra with an increasing precision. As starting points, mean-field theories play a large role because of their ability to find approximate solutions to the many-body Schrödinger equation at a relatively low computational cost [5]. A typical example is the Hartree-Fock-Bogoliubov (HFB) method where the wavefunction solution is written as a product of independent quasiparticles. Nuclear Physics applications of HFB coupled with symmetry-restoration methods have allowed to accurately describe a vast swath of nuclear properties such as binding energies, mean-square radii, deformation and spectra. Under certain circumstances, such as for instance in the case of shape coexistence, it becomes necessary to include additional correlations between the quasiparticles. Correlations beyond the mean-field can be included, for example, using the framework of the Generating Coordinate Method (GCM) [5][6][7] where the wavefunction is written as a linear combination of different non-orthogonal mean-field configurations. In practice, the solution of GCM equations and the application of symmetry-restoration methods require the precise computation of overlap functions. The modulus of the overlap between two HFB vacua can be computed with the Onishi formula [8] thus leaving ambiguity on its sign. In the past, several studies have been dedicated to overcome this issue [9][10][11][12][13][14][15][16]. In [17], Robledo solved this issue by deriving an expression for the overlap given in terms of a Pfaffian. In a later publication [1], Bertsch and Robledo extended the Pfaffian formulation to overlaps between odd-A systems and overlaps for operators involved in symmetry restoration methods. Let us consider two HFB vacua |Φ and |Φ with even number parity (for an even-A nucleus) and the associated overlap O ≡ Φ|Φ . The formula by Bertsch and Robledo [1] gives O as: where N is the (even) number of single-particle basis states (n, l, j, m) and U, V (U , V ) are matrices of the Bogoliubov transformation associated with |Φ (|Φ ) [5].
Using the Bloch-Messiah decomposition [18], one can write U = DŪ C and V = D * V C where D and C are both unitary matrices and D defines the so-called canonical basis associated with |Φ .Ū andV can be chosen as diagonal and skew-symmetric, respectively. The matrix V is written in terms of N/2 blocks of dimension 2 × 2 with elements (v i , −v i ) where v 2 i is the occupation probability of the canonical basis state i (the matrix elements u i ofŪ are such that u 2 i + v 2 i = 1). We adopt the usual phase convention u i > 0 and v i > 0 [5]. For v i = 1, the level i is fully occupied whereas it is empty for v i = 0. v i and C in (1) are obtained by the Bloch Messiah decomposition of (U , V ). In practice, due to the tiny values of v i , v i for the least occupied levels, the computation of the overlap (1) can become unstable. Indeed, in this context, both the denominator i,i v i v i and the Pfaffian in the numerator of Eq. (1) have small values that can become out of the scope of the double precision data type. A potential solution to this issue could be to discard levels for which v i , v i ∼ 0 in the computation of the overlap. But in practice, it can become difficult to check the reliability of such an approach. Indeed, one would need to check that the omitted levels have a negligible contribution in the value of the overlap and this would be done by increasing the number of levels considered. Eventually one might again run into instability of the numerical computation. It is especially important to obtain precise values of overlaps in the context of GCM where the basis states are not orthogonal and the overlaps can be arbitrary small. Other solutions to this issue has been arXiv:2010.08459v1 [nucl-th] 16 Oct 2020 proposed in [19] where small values of v i are replaced by an ad hoc tiny numerical parameter . However, it would be advantageous and more convenient for systematic calculations to be able to bypass the introduction of such a parameter.
In this letter, we present a new and practical formulation of the overlap (see Eq. (10) in the case of an even-A system and Eq. (12) for an odd-A system) which allows for a precise and stable computation and is also amenable to controllable truncations. This letter is organized as follows. We first show the main steps involved in the derivation of the formula for the overlap between two HFB vacua for an even-A system and then show the expression for the overlap in the case of an odd-A system. We then illustrate the numerical stability offered by this formulation by computing, as a function of the number of canonical basis states included, the matrix element of the particle-number projection operator.
We start the derivations from Eq. 1. First, we want to point out that this equation can not be used directly when including empty levels since in that case, the expression involves products of zero and infinity. The equation is however generally valid if one assumes a tiny occupation for the empty levels. In our final expression Eqs. (10,12) these factors cancel out analytically and one may safely evaluate the expression in the limit → 0.
We denote M the matrix argument of the Pfaffian in Eq. (1). Using the Bloch-Messiah decomposition for U, V and U , V we obtain: From Eq. (2) and the relation we can write: In order to avoid the numerical instability that arises in the computation of the overlap directly from Eq. (1) we factorize the norms of the HFB vacua out of pf(M ). This is achieved in the following steps by first writing the matrix argument of the Pfaffian on the right-hand side of Eq. (4) as: Introducing the diagonal matrix and a similar matrix Λ with the occupation number v i , we can rewrite (5) as : From the expression above, we can now write Eq. (4) as: where Eq. 9 has been obtained from Eq. 8 using once again the relation ( (1), we now arrive at an expression of the overlap as: where σ is the N × N tridiagonal skew-symmetric matrix with elements 1 and -1 [20]. The expression (10) is exact and allows for a stable numerical computation of overlaps independently of how tiny the occupation number v i might be (this includes the limit of vanishing occupations e.g. in the case of no pairing, where occupation numbers are exactly 0 for levels with energy greater than the Fermi energy). Moreover, the formula enables to systematically reduce the computing cost by including only levels for which the occupation number v i (v i ) is greater than a given value η and accordingly truncating the matrix in Eq. (10). The quality of the truncation can then be checked a posteriori by decreasing η to arbitrary small values. Let us assume that for a given value of η, n (n ) canonical basis states fulfill the criteria v i ≥ η (v i ≥ η).
In that case, the values of the overlap in this truncated space is given by [20]: The extension for odd-A systems is straightforward and the derivation goes along the same line (see [20]). Let us consider the overlap between two odd-A states written as one quasi-particle creation operators acting on even-number-parity HFB vacua that is, β † a |Φ and β † a |Φ . The corresponding overlap Φ|β a β † a |Φ is given as: with matrix related to the even part of the system, related to the connection between the even and the odd part and related to the odd particles. The index a (a ) in the equations above is used to denote the column vectors associated with the two blocked quasi-particles. As in the case of the overlap between even-number-parity vacua, the expression for the overlap in the odd case is numerically stable and allows for controllable truncations. In this context, the truncated expression for the overlap can be written as [21] (see [20]): Example. As a proof of principle, we now show results for the computation of a matrix element involved in the particle-number projection method for an even-A nucleus. We denote this matrix element O PN : where e −iN φ is the gauge rotational operator for the particle number projection,N the particle number operator and φ the gauge angle. We focus on projecting the number of neutrons and generate the occupation numbers with a HFB calculation of 192 Pb using a basis of 16 major oscillator shells, the SLY4 functional [22] and a separable Gaussian pairing interaction [23]. We rewrite O PN as the overlap between the vacuum |Φ and a rotated vacuum |Φ with the corresponding Bogoliubov-transformation matrices U = U e iφ and V = V e −iφ . We set φ = π/10 and for convenience, we randomly generate the unitary matrices D and C. We calculate O PN with Eq.  cally unstable. In such a case the computation limit is reached at n = 323 for which v i = 0.097593. With the new formula all states can be included without encountering any numerical instability. The ability to include all states is expected to be particularly important for the description of the radial tail of the wavefunction. A correct description of the wavefunction at large distances is in turn critical for many applications involving scattering and reaction processes [24] such as alpha decay [25] or the computation of nucleon-nucleus optical potential [26,27].
Conclusion. In this Letter we have presented a new, exact and numerically stable formulation of overlap functions that appear in beyond mean-field theories such as the GCM method. We have illustrated the benefits offered by this formulation by computing a matrix element for the particle-number projection method. We want to emphasize that the derivations presented here is applicable to overlap functions in the context of both even and odd number parity. The new formula offers substantial improvements over current methods since exact values can now be obtained even in large spaces without being limited by the ability of representing very small or large numbers. It adds no extra computational effort or complexity and may be truncated in a systematic way, which allows for a smooth convergence towards the exact value. Finally, the new formula opens the door towards precision calculations of nuclear reactions and nuclear structure where one takes advantage of the power of symmetry breaking.