A convenient implementation of the overlap between arbitrary Hartree-Fock-Bogoliubov vacua for projection

Overlap between Hartree-Fock-Bogoliubov(HFB) vacua is very important in the beyond mean-field calculations. However, in the HFB transformation, the $U,V$ matrices are sometimes singular due to the exact emptiness ($v_i=0$) or full occupation ($u_i=0$) of some single-particle orbits. This singularity may cause some problem in evaluating the overlap between HFB vacua through Pfaffian. We found that this problem can be well avoided by setting those zero occupation numbers to some tiny values (e.g., $u_i,v_i=10^{-8}$). This treatment does not change the HFB vacuum state because $u_i^2,v_i^2=10^{-16}$ are numerically zero relative to 1. Therefore, for arbitrary HFB transformation, we say that the $U,V$ matrices can always be nonsingular. From this standpoint, we present a new convenient Pfaffian formula for the overlap between arbitrary HFB vacua, which is especially suitable for symmetry restoration. Testing calculations have been performed for this new formula. It turns out that our method is reliable and accurate in evaluating the overlap between arbitrary HFB vacua.


Introduction
The Hartree-Fock Bogoliubov (HFB) approximation has been a great success in understanding interacting many-body quantum systems in all fields of physics. However, the beyond mean-field effects (e.g., the nuclear vibration and rotation) are missing in the HFB calculations. Methods that go beyond mean-field, such as the Generator Coordinate Method(GCM) and the projection method, are expected to take those missing effects into consideration and present better description of the many-body quantum system. In the beyond meanfield calculations, operator matrix elements and overlaps between multi-quasiparticle HFB states are basic blocks. These matrix elements and overlaps can be evaluated using the generalized Wick's theorem (GWT) [1,2], or equivalently using Pfaffian [3,4,5,6,7], or using the compact formula in Ref. [8]. However, in the efficient calculations (e.g., see [5]), all of the matrix elements and overlaps require the value of the overlap between HFB vacua.
Thus, the reliable and accurate evaluation of the overlap between HFB vacua is very important for the stabil-ity and the efficiency of the beyond mean-field calculations. Especially in cases near to the Egido pole [9], the overlap between HFB vacua is very tiny, and a small error could lead to a large uncertainty of the matrix elements. In the past, numerical calculations of the overlap were performed with the Onishi formula [10]. Unfortunately, the Onishi formula leaves the sign of the overlap undefined due to the square root of a determinant. Several efforts have been made to overcome this sign problem [11,12,13,14,15,16]. In 2009, Robledo proposed a different overlap formula with the Pfaffian rather than the determinant [17]. This formula completely solves the sign problem but requires the invert of the matrix U in the Bogoliubov transformation. To avoid the singularity of U, the formula for the limit when several orbitals are fully occupied is given in Ref. [18]. Simultaneously, the limit when some orbits are exact empty was also considered to reduce the computational cost. Meanwhile, various Pfaffian formulae for the overlap between HFB vacua have been proposed by several authors [3,6,7]. In Ref. [7], the overlap formula does not require the invert of U, but the empty orbits in the Fock space should be omitted.
In practical calculation, one should first identify the singularity of the matrices U and V in the Bogoli-ubov transformation. This can be easily tested with the Bloch-Messiah theorem (see details in Ref. [19]). The matrices U and V can be decomposed as U = DŪC and V = D * V C. Here, D and C are unitary matrices.Ū and V refer to the BCS-transformation and are constructed from the occupation numbers u i , v i with 0 ≤ u i , v i ≤ 1 and u 2 i + v 2 i = 1 (see Eq. (7.9) and Eq. (7.12) in Ref. [19]). The limits of fully occupied (u i = 0, v i = 1) and fully empty (u i = 1, v i = 0) levels have been carefully treated in Ref.s [6,7,18] to avoid the collapse of the overlap computation.
However, we note that in most realistic cases the v i 's can be extremely close to 0 or 1 but not exact 0 or 1. Strictly speaking, these levels with such extreme emptiness or occupation should be considered but may lead to exotic values ( extremely huge or extremely tiny) of the Pfaffian in the proposed formulae. What is worse, the Pfaffian values are easily out of the scope of the double precision data type and cause the computation collapsed.
Careful treatment must be made to avoid such data overflow. In this paper, we implement an accurate and reliable calculation for the overlap between arbitrary HFB vacua in a uniformed way. For the cases of (u i = 1, v i = 0) and (u i = 0, v i = 1), we treat them as the cases of (u i = 1, v i = 10 −8 ) and (u i = 10 −8 , v i = 1), respectively. Because v 2 i (u 2 i ) = 10 −16 is actually zero relative to u 2 i (v 2 i ) = 1 in practical calculations. Thus this treatment does not change the HFB vacuum at all. Therefore, without losing the generality, we assume that all levels in the Fock space are partly occupied, but some of their u i , v i values are allowed to be extremely close to 0 or 1. Ideally, U, V are nonsingular in our assumption, and we can derive a new formula for the overlap between the HFB vacua based on the work of Bertsch and Robledo [7]. This formula is especially convenient for the symmetry restoration. Numerical calculations have been carried out for heavy nuclear system to test the precision of the new formula by comparing with the Onishi formula.
In section 2, the formalism of the new overlap formula is given. Section 3 provides an example of numerical calculation. A summary is given in section 4.

The overlap between the HFB vacua
We denoteĉ † i andĉ i as the creation and annihilation operators defined in a M-dimentional Fock-space. The Hartree-Fock-Bogoliubov(HFB) transformation is Here, we assume U and V are nonsingular matrices, and their shapes are M × M. The HFB vacuum (unnormalized) can be written as where |− is the true vacuum. By definition, one haŝ The second HFB vacuum |φ ′ is defined in the same way, but the prime, ' ′ ', is attatched to the corresponding symbols to show difference. The overlap between |φ and|φ ′ is given by where, Following the technique of Bertsch and Robledo [7], one can obtain The shape of the matrix in Eq. (5) is 2M × 2M, and no empty levels are omitted. For the norm overlap φ|φ , it is real and positive. From Eq.(5) and the Bloch-Messiah theorem, one can get Denoting M/2 i=1 v i by N, the normalized quasiparticle vacuum, |ψ , can be written as Then, one finds that In the symmetry restoration, the general rotational operator, involving the spin and particle number projection, may be written aŝ whereR(Ω) is the rotation operator, and Ω refers to the three Euler angles α, β, γ. e −iNφ n and e −iẐφ p are 'gauge' rotational operators induced by the neutron and proton number projection.N andẐ are neutron and proton number operators, respectively. φ n and φ p are "gauge" angles for neutron and proton, respectively. Ξ refers to (Ω, φ n , φ p ). The matrix element ψ|R(Ξ)|ψ ′ needs to be calculated. Let's define the general rotational transformation for symmetry restoration, where where By comparing Eq.(11) with Eq. (1), one can obtain the rotated overlap by replacing U ′ and V ′ in Eq. (8) with D(Ξ)U ′ and D * (Ξ)V ′ , respectively. Thus where This formula is essentially the same as the one proposed by Bertsch and Robledo [7], but we will transform it into a new form. Supposing that there is a Ξ 0 satisfying N pf (Ξ 0 ) 0, we have where and P is Therefore, one can get where, the coefficient C is actually independent of Ξ 0 , and can be written as Here, ∆ is a phase determined by In Eq.(18), we have used the the Bloch-Messiah theorem and the following equation Eq.(18) looks more convenient to be implemented and may save some computing time in contrast to Eq. (13), where extra evaluation of V T D(Ξ)V ′ * is required for each mesh point in the integral of projection.
For comparison, let us present a brief introduction of the overlap of the Onishi formula [10]. The unitary and canonical transformation of the quasi-particles under ro-tationR(Ξ) can be written as, where The Onishi formula is then expressed as (see Ref. [20]), where, M n and M p are the numbers of neutron and proton orbits in the Fock space, respectively. The value of det[X(Ξ)] is a complex number, and the sign of the square root is left undefined. Extra efforts must be made to determine the sign before the application of the Onishi formula. For instance, in the Projected Shell Model [21] without particle number projection, the overlap between the BCS vacua is real and positive, thus there is no sign ambiguity and the Onishi formula works.

Numerical test of the overlap formulae
Although the sign problem is solved in Eq.(13) and Eq. (18). One can imagine that N, N ′ are extremely tiny numbers by definition. Thus pf[M(Ξ)] is also very tiny, but pf[W(Ξ)] should be huge. Numerical accuracy of Eq.s (13) and (18) needs to be carefully tested. It is believed that the Onishi formula is accurate except for its undetermined sign. So, it is helpful to compare the numerical values of the overlaps using Eq.(13), Eq. (18) and Eq.(24).
To demonstrate the accuracy and the reliability of the Eq.s (13) and (18), numerical calculations are performed for the typical example of the deformed heavy nucleus 226 Th. For projection, we should take |ψ = |ψ ′ , and then ∆ = 1. The U, V matrices are obtained from the Nils-son+BCS method. The single particle levels are generated from the Nilsson Hamiltonian with the standard parameters [22]. The single-particle model space contains 5 neutron major shells with N = 4 − 8 and 5 proton major shells with N = 3 − 7, i.e., the Fock space has 145 neutron levels (M n = 290) and 110 proton levels (M p = 220). The numbers of the active neutrons and protons are 96 and 70, respectively. The quadrupole deformation is taken to be ǫ 2 = 0.2. Here, we only consider the axial symmetry for simplicity.
In the no pairing case, the BCS vacuum becomes a pure slater determinant, which is a challenge for Eq. (18) because all v i 's above the fermi surface are zero. Consequently, N = 0 and W(Ξ) is meaningless due to the singularity of V. Here, we set v i = 10 −8 for those v i = 0 orbits to avoid the collapse of calculation. Therefore we have In all calculations, we found that R < 10 −12 with double precision. This confirms that a small change of v i from zero to 10 −8 almost does not affect numerical accuracy. However, it is crucial to keep Eq. (18) valid. Yet notice that N pf (Ξ) in Eq. (18) is obtained from a product of tiny and giant numbers. The same calculations has also been done with Eq.(13), and we also get R < 10 −12 . Thus we have presented an alternative way of using Eq. (13), where we set v i = 10 −8 for those empty orbits rather than omitting them [7]. Once the overlap is available, it is straightforward to perform the symmetry restoration. The deformed BCS vacuum of 226 Th has been projected onto good particle number and spin. Therefore, one can test how precise the numerical calculations with Eq. (18) satisfy whereP N ,P Z , andP I MK are neutron-number, protonnumber, and spin projection operators, respectively. For the above vacuum state without pairing (i.e. the ground state slater determinant), the particle numbers of both neutrons and protons are good. Indeed, our particle number projection (using 16 mesh points in the integral) shows that ψ n |P N |ψ n = 1 (N = 96), or 0 (N 96) with numerical errors less than 10 −13 . Calculations for the protons also have the same accuracy. This again shows the reliability of Eq. (18). Angular momentum projection is also performed on the same state in addition to the particle number projection. The amplitude of ψ|P NPZPI 00 |ψ with (N = 96, Z = 70) is plotted as a function of spin I in Fig.2. In the integral of the spin projection, 100 mesh points are taken, and the range of spin is 0 ≤ I ≤ 70, and we indeed reproduced Eq.(26) with numerical error around 10 −12 .
We also have tested Eq. (18) in the projection of the triaxially deformed vacuum with normal pairing, which seems more convenient to use Eq. (18). With the present method, similar accuracy has also been achieved. Spin I Figure 2: The amplitude of projection, ψ|P NPZPI 00 |ψ , as a function of spin I at N=96 and Z=70 using Eq. (18). |ψ is the axially deformed BCS vacuum but without pairing.

Summary
Following the strategy of Bertsch and Robledo [12], we have proposed a new formula of the overlap between HFB vacua by using the Pfaffian identity and assuming that the inverse of the V matrix exists. This formula is especially convenient and efficient in the symmetry restoration, and has the same high accuracy as the Onishi formula as well as the correct sign. The reliability of the present formula has been tested by carrying out the calculations of the overlap and the quantum number projection for the heavy nucleus 226 Th. In the testing calculations, one has to face with two numerical problems: (1) The extreme (huge or tiny) quantities are certainly encountered, and we have properly treated this situation to avoid data overflow (see the text). (2) For those empty orbits with v i = 0, which make Eq. (18) invalid, one can change v i to a small nonzero value (e.g., 10 −8 ) to avoid this singularity of V matrix. It turns out that such treatments work very well. Testing calculations have confirmed that the present formula is even applicable to the pure slater determinant without losing the numerical accuracy. Thus it is promising that Eq.(18) may be applicable between arbitrary HFB vacua.