Why does the sign problem occur in evaluating the overlap of HFB wave functions?

For the overlap matrix element between Hartree-Fock-Bogoliubov states, there are two analytically different formulae: one with the square root of the determinant (the Onishi formula) and the other with the Pfaffian (Robledo's Pfaffian formula). The former formula is two-valued as a complex function, hence it leaves the sign of the norm overlap undetermined (i.e., the so-called sign problem of the Onishi formula). On the other hand, the latter formula does not suffer from the sign problem. The derivations for these two formulae are so different that the reasons are obscured why the resultant formulae possess different analytical properties. In this paper, we discuss the reason why the difference occurs by means of the consistent framework, which is based on the linked cluster theorem and the product-sum identity for the Pfaffian. Through this discussion, we elucidate the source of the sign problem in the Onishi formula. We also point out that different summation methods of series expansions may result in analytically different formulae.


I. INTRODUCTION
The Hartree-Fock-Bogoliubov (HFB) theory gives a simple but profound basis for the nuclear many-body problem where the competition between the nuclear pairing and deformation plays a primary role in the determination of the ground state as well as the excited states. Especially, a combination of the HFB method with the technique of angular momentum projection allows the direct comparison between the theoretical calculations and the experimental data. The projected HFB states can produce more elaborate and accurate calculations although the simplicity of the HFB wave functions is kept from a mean-field point of view. In this way, not only the simple HFB state but also a superposition of different HFB states (i.e., the projected HFB state) have been extensively used for nuclear-structure studies. Behind this success of the HFB theory, there was a hidden problem for the overlap matrix element between the HFB states.
Half a century ago, a formula for the overlap matrix element was derived by Onishi and Yoshida [1] and is called the Onishi formula [2]. To derive the Onishi formula, we begin with the Thouless representation [2], [3] of the HFB wave functions, |φ (k) (k = 0, 1) defined as where c † 's are the creation operators and |− is the bare Fermion vacuum with c i |− = 0 (i = 1, · · · , N ). The dimension of the Fermion single-particle space is N . Z is an N ×N complex skew-symmetric matrix. The Thouless representation is a specific one of the Bogoliubov quasiparticle states. In this representation, the overall phase is fixed for two HFB wave functions |φ (0) and |φ (1) , respectively, as in Refs. [4], [5]. The overlap matrix element between these two HFB wave functions is defined as which can be expressed as This formula is known as the Onishi formula [1]. Due to the square root function, the Onishi formula is twovalued and does not give a definite sign if Z's are complex matrices. This indefiniteness of the sign assignment is referred to as the sign problem of the Onishi formula, which becomes quite serious in the application of the full angular momentum projection. So far, there are several approaches known to remedy the problem [4], [5], [6], [7], [8].
Among them, Robledo [5] has recently derived an alternative and ambiguity-free formula for the overlap matrix element by the Pfaffian as φ (0) |φ (1) where s N = (−) N (N +1)/2 and I is the N ×N identity matrix. This formula is proved with rather advanced techniques, that is, the Fermion coherent state and Grassmann integral. His proof is mathematically very elegant and interesting [5]. Moreover, these techniques led us to a relation to the generalized Wick's theorem and its related topics [9], [10], [11], [12], [13]. The proof by Robledo is, however, rather abstract and it somewhat keeps us from an intuitive understanding of the reason for the disappearance of the sign ambiguity.
In the present paper, we derive both formulae in Eqs. (3) and (4) directly from Eq.(2) and elucidate an origin of the sign problem. First, we expand the exponential operators in Eqs. (1,2). After handling the vacuum expectation values of the product of the creation-annihilation operators, the overlap matrix element can be, in principle, expressed as a polynomial of the matrix elements of Z. The overlap is, therefore, single-valued.
Next, we will consider two summation methods. We show that an expansion of the HFB wave function in Eq.(1) can be expressed by the Pfaffians and that the overlap matrix element can be, thereby, revealed by a finite series of the product-sum of the Pfaffians. This finite series can be summed up into Robledo's Pfaffian formula in Eq. (4). By this derivation, Robledo's Pfaffian formula is turned out to be obviously single-valued and to be free of the sign problem. We also show that the other summation method with the linked cluster theorem [14] brings us to the Onishi formula. We present that this summation concerning the connected diagrams involves an infinite series, which is in sharp contrast to the original finite series and that it gives rise to the square root function in the Onishi formula. We also clarify that the skew-symmetric property of the Thouless matrix Z can remove the sign problem from the Onishi formula completely.
The present paper is organized as follows. In Sec. II, we show a basic structure of the overlap matrix element through the series expansion of the HFB wave functions and present an alternative derivation of Robledo's Pfaffian formula. In Sec. III, we show a relation between the Onishi formula and the linked cluster theorem, and we discuss the origin of the sign problem. In Sec. IV, we give a conclusion. In the appendices, we summarize useful identities concerning the Pfaffian and show the derivation of the connected term.

II. OVERLAP FORMULA WITH THE PFAFFIAN
A. Basic structure of the overlap matrix element First, we show a basic mathematical structure of the overlap matrix element by expanding the exponential operators in Eq. (2). Defining pair-annihilation and paircreation operators,Â andB aŝ the HFB wave functions are shown by The overlap matrix element is simply denoted by By expanding exponential operator, the HFB wave function can be shown as where this series expansion terminates in order N/2 because the number of single particle states, namely, the dimension of the matrices Z 0 and Z 1 is N . The overlap matrix element is rewritten by Note that −|Â l l!B k k! |− vanishes if l = k becauseÂ and B are pair-annihilation and pair-creation operators, respectively.
Next, we investigate an expanded form of the 1 k operator is generally expressed by the 2k creation operators with the coefficients in terms of matrix elements of Z 1 as, As theÂ k operator is also similarly shown, the overlap matrix element can be straightforwardly expressed as a function of matrix elements of Z 0 * and Z 1 as, The matrix element in the third line of the above equation gives intricate restriction concerning p's, q's, p ′ 's, q ′ 's by taking the contractions. The above formula shows a very complicated structure regarding the matrix elements of Z 0 * and Z 1 . It is, however, quite evident that the overlap matrix element is a polynomial of the matrix elements of Z 0 * and Z 1 and has, thereby, no sign ambiguity.
In the subsequent subsections, we will show that Eq. (12) can be rewritten by the Pfaffians and will directly derive Robledo's Pfaffian formula from Eq. (12). Furthermore, in the next section, by handling Eq.(12) with the linked cluster theorem, we will derive the Onishi formula.

B. Overlap formula with product-sum of the Pfaffians
In this subsection, we consider to rewrite Eq.(12) by investigating the detailed structure of Eq. (11).
For example, the 2nd order term is expressed as where 2! is removed due to the additional condition p 1 < p 2 . Now let us change the integer indices p 1 , q 1 , p 2 , q 2 to new indices n 1 , n 2 , n 3 , n 4 with n 1 < n 2 < n 3 < n 4 , and let us obtain the coefficients b {n} in the following form; where n's run 1 to N under the condition n 1 < n 2 < n 3 < n 4 . These two kinds of indices have different conditions. The relation between them is classified into the following three cases: . Therefore, the coefficients b {n} can be given in terms of the Pfaffian as where we use Eq.(A4). In general, the 1 k!B k operator is expressed by the 2k creation operators with the coefficients in terms of matrix elements of Z 1 as, where the summations are performed with the restriction, p 1 < q 1 , · · · , p k < q k , p 1 < · · · < p k . The k! is removed due to the additional condition p 1 < · · · < p k . As the same procedure, we change the integer indices p 1 , q 1 , · · · , p k , q k (p 1 < q 1 , · · · , p k < q k , p 1 < · · · < p k ) to 2k different integer indices n 1 , · · · , n 2k (n 1 < · · · < n 2k ). As this condition is the same as Eqs.(A2,A3), we can introduce the Pfaffian as where n's run 1 to N under the restriction, n 1 < · · · < n 2k and the m × m skew-symmetric matrix Z 1 m is defined as The HFB wave functions are expressed as where i takes 0 or 1. The expansion of the HFB wave function can be generally shown by the Pfaffians. The overlap matrix element, φ (0) |φ (1) , is given in terms of the Pfaffians as By introducing an index set I = {n 1 , n 2 , · · · , n 2t } (n 1 < · · · < n 2t ), we can define the sub-matrix Z 1 2t as Z 1 I , which is defined as where i, j run 1 to 2t. With this notation, the overlap matrix element is compactly expressed as where I N 2t is a set, consisting of subsets with 2t elements in {1, 2, 3, · · · , N }. Next, we show that the Pfaffians in Eq.(22) can be expressed with single Pfaffian, thanks to the productsum identity of the Pfaffian, Eq.(B1) in Appendix B. We define two skew-symmetric matrices P and Q as The l.h.s. of Eq. (B1) is rewritten as which is Robledo's Pfaffian, except s N . As the dimension of Z 0 and Z 1 is N , m = 2N in Eq.(B1). The r.h.s. of Eq.(B1) becomes where I 2N 2r is a set, consisting of subsets with 2r elements in {1, 2, · · · , 2N }. Below, we will show that Eq.(25) can be reduced to Eq.(22) by the Pfaffian identities and classification of the indices in Eq.(25).
As the matrices P and Q have a bipartite structure, we divide the I into two parts I 1 and I 0 , We classify the I by the symmetry with the bipartite structure. One is symmetric concerning I 1 and I 0 and is denoted as I s . The I s with numbers of elements 2r is a sum of I 1 = {l 1 , l 2 , · · · , l r } and I 0 = {l 1 + N, l 2 + N, · · · , l r + N }. The {1, 2, 5, 6} is an example of I s . The other is asymmetric and is called I a , whose examples are {1, 3, 7, 8}, {2, 4, 5, 8}. Only these symmetric sets I s with even integer r can contribute to Eq.(25), while other sets, namely, symmetric sets I s with odd integer r and asymmetric sets I a give no contribution.
In the former case, I s with 2r (r:odd), matrix dimensions of Z 1 where we use Eq.(A7). Therefore, the number of elements of I s , giving non-vanishing Pfaffian, is 4t (t:integer).
In the latter case, if the numbers of elements of I a 1 and I a 0 are different, Pf (Q I a )=0 because Q I a is not a square matrix. For example, we again take the case of {1, 2, 4, 5} for N = 4 and r = 2. The numbers of elements of I a 1 and I a 0 are 3 and 1, respectively. The matrix Q I a is 1 × 3 and is not a square one. If the numbers of elements of I a 1 and I a 0 are the same and are 2t, P f (Q I a )=0 due to asymmetry of index. For instance, we again take the case of {1, 3, 7, 8} for N = 4 and r = 2. Its complementary set is {2, 4, 5, 6}. The Q I a is given by    q 2,2 q 2,4 q 2,5 q 2,6 q 4,2 q 4,4 q 4,5 q 4,6 q 5,2 q 5,4 q 5,5 q 5,6 q 6,2 q 6,4 q 6,5 q 6,6 because of the definition of Q, Eq.(23). Therefore, its Pfaffian is zero. In general, Q I a is shown in the bipartite form as where the diagonal block matrices are zero and the offdiagonal block matrix is denoted by C. By the identity Eq.(A8), it reduces to (−1) m(m−1)/2 Det [C] where m = N − 2t. Let I a 1 and I a 0 be {i 1 , i 2 , · · · , i m } (i 1 < i 2 < · · · < i m ) and {j 1 + N, j 2 + N, · · · , j m + N } (j 1 < j 2 < · · · < j m ), respectively. As some indices are the same and others are different, we re-sort the indices as {i 1 , i 2 , · · · , i m } → {i ′ 1 , · · · , i ′ k , · · · , i ′ m } and In this representation, the matrix C can be shown as while the diagonal block matrices in Eq.(28) are still zero. Therefore, P f (Q I a ) = 0 is proved. Next we consider I s with 4t (t:integer), which consists of I 1 s = {l 1 , l 2 , · · · , l 2t } and I 0 s = {l 1 +N, l 2 +N, · · · , l 2t + N }. For instance, such a case is {1, 2, 5, 6} for N = 4 and t = 1(r = 2). The Pfaffian of P is rewritten as where Z 1 where I I s 1 is a (N − 2t) × (N − 2t) identity matrix and we use Eq.(A8). Thus, the product-sum identity of the Pfaffian is reduced into As |I s | = 2|I s 1 | + 2N t, (−1) |I s | = 1. Thereby, the sign of r.h.s. of Eq.(32) becomes (−1) 1 2 N (N +1) , which is just Robledo's s N [5]. Therefore, the obtained overlap matrix element agrees with Robledo's Pfaffian expression as, Thus, we algebraically derived Robledo's Pfaffian formula, summing up the expansion terms by applying the product-sum identity of the Pfaffian. As this expansion forms a finite series and is a polynomial, it is evident that its summation is also single-valued and the obtained formula has no sign problem.

III. OVERLAP FORMULA WITH THE DETERMINANT
The Onishi formula was first obtained by Onishi and Yoshida [1], and Onishi and his collaborators used the linked cluster expansion [14] for the double-variational method [15], [16]. Here, we derive the Onishi formula based on the linked cluster theorem [14], which is more standard in view of the quantum many-body theory. Especially we focus on the origin of square-root function.

A. Expansion of overlap matrix element with the contractions
Let us begin with the overlap matrix element in Eq. (2), which is rewritten as where N is a dimension of the model space. This is the same equation as Eq.(10) while here we denote a vacuum expectation value simply as Ô ≡ −|Ô|− wherê O is an arbitrary operator. As Eq.(34) forms a finite series, it seems to be difficult to derive square-root function. However, by evaluating Eq.(34) with contractions, we naturally obtain an infinite series. Before we discuss the general term, we explicitly show several terms for k = 0 ∼ 3. The 0-th order term is unity.
For k = 1, taking the contractions, we can obtain the 1st order term as, where Y ≡ Z 0 * Z 1 . The contribution only from the connected diagrams is denoted with the suffix "c" attached to the expectation value. For this notation, the 1st order term is shown as, where the 1st power of Y is denoted as Y 1 for later convenience. The second-order term is shown by The contractions are classified into two groups. One is a disconnected term and is a product of 1st-order connected terms. The other is a connected term. The former one is shown by where 2! is the number of the repeated diagrams. The other is a connected term, which corresponds to where the coefficient 1 −2·2 is explained in Appendix C. Therefore, the second order term is given by The 3rd-order term has two disconnected terms as ÂB 3 c and 2! c ÂB c and one connected term as 3! c . One of the disconnected terms is shown by The other disconnected terms are shown by The connected term is shown by where the coefficient 1 −2·3 is explained in Appendix C. The 3rd order term is, therefore, given by (44) In general, we consider k-th order term, which has several disconnected terms and one connected term.
The disconnected terms are shown as ÂB k c , 1 , · · ·. The p q power of q-th connected term, 1 The k-th term is thereby expressed as In Eq.(46), the connected term is shown as where the coefficient 1 −2k is explained in Appendix C.

B. Onishi formula via the linked cluster theorem
According to the linked cluster theorem, the logarithm of the overlap matrix element is expressed with its connected diagrams as The contribution of connected diagrams for the overlap matrix element is given with the summation of all connected terms Eq.(47) as where this expression becomes an infinite series because the connected terms, which is a sum of the connected terms and the disconnected terms as shown in Eq.(46), are always zero for k > N/2 due to Fermi statistics. As a result, the overlap matrix element is also expressed by an infinite series through Eq.(48), although Eq.(34) is a finite series. This fact means that the overlap matrix element can be expressed in two analytically different ways.
Next, we continue to investigate Eq.(49). The eigenvalues of the Y matrix are denoted as e i (i = 1, · · · , N ). Then Tr(Y ) = N i=1 e i . The summation of the connected terms is shown as Here we consider the complex logarithmic function ln(1 − z) by its power series ln(1 − z) = −z − z 2 2 − z 3 3 − · · ·, which has the convergence radius |z| < 1. For the domain beyond the convergence radius, the logarithmic function can be defined by the analytic continuation, except on the singularity at z = 1. Note that, as for z = 1, the logarithmic function diverges, this point cannot be defined and is a singularity. The existence of this singularity was reported as the nodal line (a collection of zeros of the norm overlap) through the numerical investigation [7], which was found to be the major obstacles in the implementation of the Hara-Hayashi-Ring method [6].
where pv. means principal value and square root function appears. In the case of the inverse function of the logarithm as a complex function, that is, e lnz = z, the infinite-multivalued nature of complex logarithm vanishes. In this case, due to the factor of 1 2 , we have to discriminate the Riemann surface as pv.lnξ This is the well-known form of the Onishi formula [1]. Clearly, the present expression contains the square root function and it suffers from the sign problem. As discussed above, the origin of the square root is due to the infinite series expansion concerning the connected diagrams, which results in the logarithm with the 1 2 factor. Now we go back to Eq.(50), where we use the eigenvalues e i 's of the Y matrix. As this Y matrix is a product of two skew-symmetric matrices, the eigenvalues of such a matrices are doubly-degenerated [17]. Therefore, we can rewrite Eq.(50) as, where the eigenvalues of the Y matrix are pair-wisely denoted as {e i , e i } (i = 1, · · · , N/2) andξ = . This double degeneracy of the eigenvalues was also rediscovered by K. Neergård and E. Wüst [4], who proved it indirectly. This doubly degenerate nature cancels the 1 2 factor of the logarithm in Eq.(50) that is the direct origin of the sign problem. Eq. (51) can be, thereby, changed to e eÂeB c = e lnξ =ξ.
(54) Therefore, the skew-symmetric property of Z 0 and Z 1 matrices, which comes from Fermi statistics, can fundamentally and completely remove the sign problem from the Onishi formula. Thus, paying attention to such a double degeneracy of the eigenvalues, we have directly derived the sign-problem-free version of the Onishi formula from the definition of the overlap matrix element, Eq.(2).

IV. SUMMARY
We investigated two kinds of analytically different formulae for the overlap matrix elements between HFB wave functions. One is the Onishi formula that was derived half a century ago [1]. This formula is, in general, twovalued as a complex function and has the sign problem. The other is Robledo's Pfaffian formula [5], which has been derived recently. This formula is single-valued and is free of the sign problem. It is theoretically interesting to investigate why there exist two analytically different formulae and why the sign problem occurs only in the Onishi formula.
To understand both formulae more deeply, we began with the Thouless representation [2], [3] of the HFB wave function in Eq.(1) and the overlap matrix element in Eq.(2). By a naive series expansion of the exponential operators in Eq.(2), the overlap matrix element can be, in principle, expressed with a polynomial with respect to Z's matrix elements due to the Fermi statistics, as shown in Eq. (12). Thereby, the overlap is essentially singlevalued although such a simple expansion does not give rise to any useful formula. Hence, we investigate various summation methods of the series expansion.
First, by expanding the exponential operator in the HFB wave function in Eq.(1), we found that the HFB wave function can be expressed with the Pfaffians in Eq. (19), which is a finite series due to Fermi statistics. The overlap matrix element can be, thereby, rewritten by the product-sum form of the Pfaffians in Eq.(22), which is also a finite series as the naive expansion. Thanks to the product-sum identity of the Pfaffian Eq.(B1), we can sum up these Pfaffians, and finally, we succeeded in algebraically deriving Robledo's Pfaffian formula [5]. This derivation shows a relation between the finite series expansion and Robledo's Pfaffian formula [5], as well as its single-valued property and the sign-problem-free nature.
Next, starting with the overlap matrix element in Eq.(2), we evaluated the Onishi formula with the linked cluster expansion [14]. We investigated the summation procedure in the series expansions of the overlap matrix element in detail, where an infinite series of the connected diagrams shows up. This infinite summation can alter the analytical property. In fact, the summation for the connected diagrams leads to the logarithm with the factor 1 2 . As the overlap matrix element is given by the exponential function of the summation of the connected diagrams, this factor 1 2 results in the square root function, which can be considered as the origin of the sign problem. We also pointed out that the sign problem is completely clarified by paying attention to a mathematical fact that the eigenvalues of a product of two skew-symmetric matrices are always doubly-degenerated [17]. More details and further considerations about this aspect are to be discussed elsewhere [18]. This double degeneracy removes the factor 1 2 , and the sign problem disappears. Finally, through this study, we showed a case that the linked cluster theorem gives a different analytic expression for the matrix element from its original analytic property. It may be interesting to find other cases with regard to the sign problem caused by the use of the linked cluster theorem.
For an n × n (n = odd) skew-symmetric matrix, Pf [A] = 0. For a 2 × 2 skew-symmetric matrix, Pf [A] = a 12 . For a 4 × 4 skew-symmetric matrix, Pf [A] = a 12 a 34 − a 13 a 24 + a 14 a 23 . (A4) For a skew-symmetric matrix A with dimension 2n × 2n, the following relations hold as and where Q is an arbitrary 2n × 2n matrix. The Pfaffian of skew-symmetric block diagonal matrix A with dimension 2n× 2n becomes a product of the Pfaffians of n × n sub-matrices A 1 and A 2 as For a special block skew-symmetric matrix with n × n sub-matrices C, the following identity holds as, where I = {i 1 , i 2 , · · · , i 2r } is a subset with 2r elements of {1, 2, · · · , m} and || means sum of the elements, that is, |I| = i 1 + i 2 + · · · + i 2r . I is the complementary set of I concerning {1, 2, · · · , m}. This proof is given by papers in pure mathematics [19] and [20].