Technical aspects of the evaluation of the overlap of Hartree- Fock-
Bogoliubov wave functions

Show preview

Hide preview

ar X
iv :1
10 2.
21 51
v1 [
nu cl-
th] 1
0 F eb
20 11
Technical aspects of the evaluation of the overlap of Hartree- Fock- Bogoliubov wave
functions
L.M. Robledo∗
Departamento Física Teórica, Universidad Autónoma de Madrid, E-28049 Madrid, Spain
Several technical aspects concerning the evaluation of the overlap between two mean field wave functions of the Hartree Fock Bogoliubov type, are discussed. The limit when several orbitals become fully occupied is derived as well as the formula to reduce the dimensionality of the problem when exactly empty orbitals are present. The formalism is also extended to deal with the case where the bases of each of the wave functions are different. Several practical results concerning the evaluation of pfaffians as well as the canonical decomposition of norm overlaps are also discussed in the appendices.
I. INTRODUCTION
In a recent publication [1] a new formulation, based on the pfaffian of a skew-symmetric matrix, has been pro- posed to compute the overlap between mean field wave functions of the Hartree Fock Bogoliubov (HFB) type, in- cluding its phase (or sign, for real overlaps). The result for the overlap is obtained by recursing to the powerful concept of fermion coherent states [2–6] and it involves the evaluation of a quantity called the pfaffian of a skew- symmetric matrix -see, for instance [7] for a definition of the pfaffian in a physical context- that is similar in spirit (linear combinations of products of matrix elements) to the determinant of a general matrix. In the derivation of the formula for the overlap [1] it is assumed that the two mean field wave functions (of the HFB type -see [8] for definition and properties-) can be related to a common reference one (usually chosen as the true particle vacuum) by means of their Thouless parametrization. It is very likely to find cases where the Thouless parametrization is ill defined because, for that particular case, it involves the inverse of a (near) singular matrix. This situation corresponds to the presence of particles (or quasiparti- cles in the general case) that have an occupancy of one, rendering the wave function orthogonal to the common (or reference) wave function. To handle those singular cases it was suggested in [1] to just change the common reference mean field wave function in order to modify the occupancies with the hope that none of them will be close to one. However, it is desirable to have an al- ternative for those cases where the change of reference wave function is either not possible or too cumbersome to carry out. Therefore, I have considered the formal limit of occupancies going to one and I have obtained a formula which is well defined in that limit and pro- vides a sound answer to such singular limit. By using the same kind of ideas I have also handled explicitly the case where some particles (or quasiparticles) have zero occupancy and therefore they do not contribute to the overlap. The formula obtained in this case involves ma- trices of smaller dimensions than the original ones and
∗ luis.robledo@uam.es
therefore should be regarded as a thrifty alternative to the original formula for situations where the wave func- tions are expanded in huge bases. Similar manipulations to reduce the size of matrices and/or deal with full occu- pancies, have been considered in Refs [9–12]. I also ad- dress the case where each of the mean field wave functions are expressed in different single particle basis related to each other by a general (not necessarily unitary) trans- formation. The result obtained is useful, for instance, to compute the overlap of the operators for spatial transfor- mations (as translations or rotations) between arbitrary mean field HFB wave functions. The result is general enough as to allow for transformations that do not map the single particle basis into itself (non complete basis under the transformation). This was already considered in [9] in a general framework but not considering the present formulation including the phase of the overlap. In [12] the implications of considering two different bases are also address, but there it is implicitly assumed that both bases share the same block structure defining the conjugate states. This is a limiting assumption that do not hold in general for time reversal violating (crank- ing, for instance) wave functions. Finally, some useful results concerning the formal evaluation of pfaffians are discussed in the appendices. These results can be of inter- est in other branches of physics where the use of pfaffians is becoming increasingly popular [13–15]. To cover also more practical aspects, the reader is referred to [16] for a thorough description of useful algorithms to compute nu- merically and symbolically the pfaffian of arbitrary skew symmetric matrices.
The relevance of the results present here is a direct consequence of the increasing popularity of the so called “beyond mean field methods” in nuclear physics [10, 17– 21] that demand the evaluation of both the modulus and phase of the overlaps between arbitrary HFB wave func- tions. A reliable determination of the sign of the norm can also be useful in to order to pin down the location of the zeros of the HFB overlaps [22]. This determina- tion would eventually be useful to get rid of the so called “pole problem” that plagues present beyond mean field calculations [23–25].
In section II the formulas pertaining the two limits considered are derived and their implications discussed.

2 We also explicitly show how to implement the change to a different common reference HFB state that could be an easy alternative in some cases. In section IIIthe case where the two HFB wave functions are referred to different bases is discussed. Finally, in appendices III, II, andI some relevant results required in the derivations are discussed.
II. EVALUATION OF THE OVERLAP IN
VARIOUS LIMITS
To begin with, it is convenient to recall the nota- tion used in [1] as well as the formulas obtained there. The goal is to evaluate both the modulus and phase of the overlap 〈φ0|φ1〉 between two HFB wave functions |φ0〉 and |φ1〉 which are given in terms of the Thouless parametrization of a general HFB wave function
|φi〉 = exp ( 1
2
∑ kk′
M (i) kk′a
+ k a
+ k′
) |0〉 (1)
The skew-symmetric matrices M (i) = (V (i)U (i)−1)∗, of dimension N × N (N is assumed to be an even number N = 2q as required for fermions), are related to the co- efficients U (i) and V (i) of the Bogoliubov transformation defining the quasiparticle annihilation operators
α (i) k =
∑ l
U (i)∗ lk al + V
(i)∗ lk a
+ l
associated to |φi〉. The above wave functions are not normalized to one, but as to have 〈0|φi〉 = 1 instead. As shown in [1] the overlap can be written as
〈φ0|φ1〉 = (−1)N(N+1)/2pfM (2) with
M =
( M (1) −I I −M (0)∗
) (3)
is a 2N×2N matrix. To obtain the results of the present section the Bloch-Messiah decomposition of the Bogoli- ubov amplitudes [8], namely U (i) = D(i)U¯ (i)C(i) and V (i) = D(i)∗V¯ (i)C(i), is used. In the previous expres- sions, D(i) and C(i) are given unitary matrices and U¯ (i)
and V¯ (i) are real matrices with special diagonal forms. By using this decomposition we can write
M (i) = D(i)M (i) C D
(i) T (4)
where the skew-symmetric matrix M (i) C is in “skew-
symmetric diagonal” (or canonical) form
M (i) C =
( 0 M¯ (i)
−M¯ (i) 0 )
(5)
The diagonal matrix M¯ (i)has matrix elements
M¯ (i) jk =
v (i) j
u (i) j
δjk. (6)
The extreme values of the ratios v (i) j /u
(i) j are infinity for
fully occupied levels (v(i) = 1) or zero for empty levels (v(i) = 0). For further developments it is convenient to single out those values and write
M¯ (i) =
( N¯ (i) 0 0 O¯(i)
) (7)
where the diagonal matrix O¯(i)contains the K(i) diagonal elements belonging to the extreme values, infinity or zero, mentioned above. The dimension of this matrix is K(i)× K(i). Each of the two limiting cases require different considerations and hence we will from now on considered them separately.
A. Limit of fully occupied levels
In this case, there are K(i) fully occupied levels in each of the HFB wave functions |φi〉 and the corresponding diagonal elements of the matrices M
(i) (the ones cor-
responding to O¯(i) in Eq. (7)) tend to infinity. This is a serious challenge, as the overlap Eq (2) as well as the norm of the |φi〉 diverge. The divergence has to be regularized and singled out of the overlap in order to cancel it out with the diverging factors coming from the norms of the HFB wave functions. To this end we write M
(i) C = R
(i)M (i) CRR
(i)T where we have introduced the “canonical regularized” (CR) matrix
M (i) CR =
 
N¯ (i) 0 0 IK(i)
−N¯ (i) 0 0 −IK(i)
 
as well as
R(i) =
 
IN−K(i) 0 0 0
0 IK(i) 0 0
0 0 IN−K(i) 0
0 0 0 O¯(i)
  (8)
In all the cases IK(i) represents the unit matrix of dimen- sionK(i). With the above definitions the matrix M of Eq (3), which enters the expression of Eq (2) for the overlap, is factorized as
M =
( R˜(1) 0
0 R˜(0)∗
)( M
(1) CR S
−ST −M (0)∗CR
)( R˜(1)T 0
0 R˜(0)+
)
where
R˜(i) = D(i)R(i)

3 and
S = R˜(1)−1 ( R˜(0)+
)−1 . (9)
Using now the property pf(BTAB) = det(B)pf(A) we obtain
pf(M) = det(D(1))det(D(0)∗)det(R(1))det(R(0)∗)pf(M˜)
with
M˜ =
( M
(1) CR S
−ST −M (0)∗CR
)
In the case of fully occupied levels, the diagonal matricesO¯(i) introduced in Eq. (8) tend to infinity. As a consequence, the determinants det(R(i)) diverge as do
some of the matrix elements of R˜(i). The later is not a problem as just the inverse of this matrix is required in Eq. (9)
R˜(i)−1 = R(i)−1D(i)−1
and R(i)−1 is a well defined quantity
R(i)−1 =
  IN 0 00 IN−K(i) 0
0 0 0K(i)
 
As a consequence of this structure, the matrix S of Eq (9) is, in this limit, a matrix where the last K(1) rows and last K(o) columns are set strictly to zero. We can use this property together with the special structure of
the matrices M (i) CR to reduce the size of the matrices to
be considered in the evaluation of the pfaffian. This will be addressed in the next subsections in a slightly differ- ent context. The only truly diverging quantities, namely the determinants det(R(i)), cancel out when we compute
the normalized overlap 〈ϕ0|ϕ1〉/ √|〈ϕ0|ϕ0〉|2|〈ϕ1|ϕ1〉|2
as, from the previous formulas, |〈ϕi|ϕi〉|2 is proportional to |det(R(i))|2.
B. Limit of fully empty levels
Another situation often encounter in numerical ap- plications is when many orbitals have zero occupancies v2 = 0 and therefore their contribution to the overlap is zero. To disentangle those contributions and reduce, in this way, the computational cost of the evaluation of the norm it is convenient to consider the limit of fully empty levels. This limit has been considered by other authors [9–12] in the past but in the more traditional formula- tion of the overlap not including the phase. In the limit of fully empty levels I will show that the evaluation of the overlap can be reduced to considering matrices with dimension equal to the number of orbitals with non-zero
occupancy. This is an advantage as the number of non- zero occupancy levels is usually much smaller than the total dimensionality of the basis used. I will use the same notation as before, but in this case, the O¯(i) are diago- nal matrices whose diagonal matrix elements are made to tend to zero. The number of diagonal elements, K(i), represent the number of empty levels in each HFB wave function |φi〉. In this case it is convenient to reorder the single particle basis by using unitary similarity transfor- mations P23 (see appendix I) permuting blocks 2 and 3,
to write M (i) C = P
(i) 23 M
(i) CRP
(i)T
23 with
M (i) CR =
 
0 N¯ (i)
−N¯ (i) 0 0 O¯(i)
−O¯(i) 0
  =
( N (i) 0
0 O(i)
)
where the matrix N (i) has dimension 2(N−K(i))×2(N− K(i)) and O(i) has dimension 2K(i)× 2K(i). The matrix carrying out the permutation P
(i) 23 depends on the basis
as the size of each block N¯ (i) is not necessarily the same. The impact of this dependence will show up in the eval-
uation of the pfaffian where the determinant of P (i) 23 is
required. Now, the new unitary matrix D (i) R = D
(i)P (i) 23
is introduced to write
M (i) = D (i) R M
(i) CRD
(i) T R . (10)
Now a technical detail related to the way how to handle
the situation where the bipartite structure ofM (0) CR is dif-
ferent from the one ofM (1) CR has to be considered. The dif-
ferent bipartite structure happens when the dimensions K(0) and K(1) differ, or in other words, that the number of empty levels in each mean field wave function is differ- ent. The optimal situation in terms of simplicity is that when both K(i) are the same. This is the strategy we will use in the following by assuming the same bipartite
structure in both M (i) CR and D
(i) R with a common dimen-
sion KS which is the smaller of K (0) and K(1). With
this in mind, we endow the matrices M (i) CR with the new
bipartite structure
M (i) CR =
( N
(i) C 0
0 O(i)
) (11)
where N (i) C are 2(N −KS)× 2(N−KS) matrices and the
O(i) have the same dimension 2KS × 2KS. In the same way, we will use below the same bipartite structure for
the matrix D (i) R
D (i) R =
( D
(i) 11 D
(i) 12
D (i) 21 D
(i) 22
) (12)
where D (i) 11 is a square matrix of dimension 2(N −KS)×
2(N −KS), D(i)12 is of dimension 2(N −KS)× 2KS , D(i)21

4 is of dimension 2KS × 2(N −KS) and finally D(i)22 is of dimension 2KS× 2KS. If we now use Eqs. (11) and (12) to reconstruct the matrix M (i) of Eq. (10) we realize
that in the O(i) → 0 limit the sub-matrices D(i)12 and D(i)22 do not enter the final expression of M (i). We can use this freedom to chose those “arbitrary” matrices in such a way as to simplify some of the results to be obtained
below. A possible, and convenient, choice is D (i) 12 = 0 and
D (i) 22 = I2KS . The matrix obtained with this choice will
be denoted by D (i)
R with
D (i)
R =
( D
(i) 11 0
D (i) 21 I2KS
) (13)
and has the nice property of having a simple inverse
D (i)−1
R =
( D
(i)−1 11 0
−D(i)21D(i)−111 I2KS
)
involving only the inverse of the matricesD (i) 11 which have
a moderate dimensionality. This is the property guiding the choice made, as we are implicitly assuming that N − KS ≪ KS and it has to be kept in mind that the cost of most of the matrix operations grow as the cubic power of their dimension. With all these definitions the matrix M which enters
the formula for the overlap is written as,
M =
( D
(1)
R 0
0 D (0)∗
R
)( M
(1) CR U
−UT −M (0)∗CR
)( D
(1)T
R 0
0 D (0) +
R
)
where
U = D (1)−1
R
( D
(0)+
R
)−1 (14)
(please note that the unitary character of the D (i) R matri-
ces is lost with the introduction of the D (i)
R ones). Using now the properties of the pfaffian we obtain
pf(M) = det(D (1)
R )det(D (0)∗
R )pf(M˜) (15)
= det(D (1) 11 ) det(D
(0) ∗ 11 )pf(M˜) (16)
with
M˜ =
( M
(1) CR U
−UT −M (0)∗CR
) . (17)
Let us now analyze the structure of the block matrix U entering the definition of M˜. Using Eq (14) together with Eq (13) allows to write
U =
( U11 U12 U21 U22
)
with
U11 = ( D
(0)+ R 11 D
(1) R 11
)−1 (18)
U12 = −U11D(0)+R 21 (19) U21 = −D(1)R 21U11 (20) U22 = I+D
(1) R 21U11D
(0) + R 21 (21)
With this definition, the matrix M˜ of Eq (28) acquires, in the limit where O(i) → 0, the block structure
M˜ =
 
N (1) C 0 U11 U12 0 0 U21 U22
−UT11 −UT21 −N (0)∗C 0 −UT12 −UT22 0 0
  .
This block structure is still not beneficial for the simpli- fication of the corresponding pfaffian and we need to use the exchange matrices defined in appendix I to exchange blocks 2 and 3
M˜R = P23M˜P T 23 =
 
N (1) C U11 0 U12
−UT11 −N (0)∗C −UT21 0 0 U21 0 U22
−UT12 0 −UT22 0
  .
Using the formulas of appendix II for the pfaffian of a bipartite matrix we obtain
Pf(M˜R) = Pf
( 0 U22
−UT22 0
)
× Pf [(
N (1) C X12
−XT12 −N (0) ∗C
)]
where
X12 = U11 + U12U −1 22 U21 (22)
The first pfaffian in the right hand side of the above ex- pression is simply given by (−1)KS detU22 whereas the second pfaffian can be computed using again the expres- sion for the pfaffian of a bipartite matrix. Collecting all the terms together we finally obtain
Pf(M˜R) = (−1)KS det(U22)Pf ( N
(1) C
) × Pf
( −N (0) ∗C +XT12N (1)−1C X12
) .
Taking into account that det(P23) = 1 we finally obtain
pf(M) = (−1)KS det(D(1)11 ) det(D(0) ∗11 ) det(U22)Pf ( N
(1) C
) × Pf
( −N (0) ∗C +XT12N (1)−1C X12
) (23)
The advantage of this ugly result over the general ex-
pression is that the dimensionality of the N (1) C , N
(0) C , D
(i) 11
and U11 matrices is 2(N − KS) which is much smaller

5 than the one of the original problem (2N). The only big matrix in Eq. (22) is the inverse of U22 with dimension 2KS, however, this does not pose a challenge as its spe- cial structure Eq. (21) is very well adapted to the use of the Sherman-Morrison formulas for the determinant and inverse of this kind of special matrices [26].
C. Using a different reference vacuum
From the above discussion, it is clear that the structure and properties of the M (i) matrices is intimately related to the reference vacuum used to express the HFB wave functions |φi〉. This suggests to use another reference vac- uum, instead of the true vacuum implicitly assumed in the previous discussions, with the hope that the new ma-
trices M (i)
will acquire a structure where there would be no fully occupied quasiparticles and the number of empty ones will be very large and comparable to the size of the basis. Matrices with that properties will not require the use of the “fully occupied limit” formulas and will benefit from the reduction in computational burden of the “fully empty limit” results. To be more specific, let us consider the example where the |φi〉 correspond to two HFB wave functions with the same average of the number of par- ticles and different quadrupole deformation parameters
q (i) 2 . It looks rather intuitive that a more appropriate ref- erence vacuum than the true vacuum could be the HFB wave function |φ¯〉 with the same average of the number of particles and a deformation parameter q¯2close to both
q (i) 2 (for instance the mean value
1 2 (q
(0) 2 +q
(1) 2 )). It is to be
expected that, with respect to this reference vacuum |φ¯〉, the new U
(i) amplitudes will be very close to the identity
matrix whereas the new V (i)
amplitudes will be small. In other words, the HFB wave functions |φi〉 will be rep- resented by a linear combination of quasiparticle exci- tations of the reference state |φ¯〉 with small amplitudes that quickly decrease with the number of quasiparticle excitations (i.e. the amplitudes of the four quasiparticle excitations much smaller that the amplitudes of the two quasiparticle ones). The expected properties of the new
U (i)
and V (i)
amplitudes imply that most of the quasi- particle excitations will correspond to the “fully empty” limit discussed in the previous subsection (with the as- sociated advantages of having to deal with matrices of small size) and far from the problematic fully occupied limit of subsection IIA. Let us consider the new reference vacuum |φ¯〉 with the
associated creation α¯+and annihilation α¯ quasiparticle operators which are defined in terms of a single particle basis of creation and annihilation operators by means of linear combinations involving the U¯ and V¯ amplitudes(
α¯
α¯+
) =
( U¯+ V¯ +
V¯ T U¯T
)( a
a+
) = W¯+
( a
a+
) .
The same relation holds true for the quasiparticle oper-
ators α(i) and α(i) + with amplitudes W (i). Using the unitarity of the matricesW (i) and W¯ we can express the set of quasiparticle operators α(i) and α(i) + in terms of the α¯ and α¯+ ones as(
α(i)
α(i) +
) =W (i)+W¯
( α¯
α¯+
) = W¯ (i) +
( α¯
α¯+
)
Using Thouless theorem we can also express the |φ¯i〉 wave functions ( satisfying 〈φ¯|φ¯(i)〉 = 1 and therefore different from the previous |φi〉 by a normalization factor) in terms of the |φ¯〉 reference vacuum
|φ¯i〉 = exp ( 1
2
∑ kk′
M¯ (i) kk′ α¯
+ k α¯
+ k′
) |φ¯〉. (24)
with M¯ (i) = (V¯ (i)U¯ (i)−1)∗ (not to be confused with the diagonal matrix of Eq (25)). These two matrices can be easily computed once the W¯ coefficients have been given and it is even possible to give an analytical expression [27] in terms of M¯ and M (i)
M¯ (i) = Q(M¯T +M (i))(I+ M¯+M (i))−1(Q+)−1. (25)
with Q = (I+ M¯M¯+)−1/2. All the formulas given above (and below) are equally valid for the wave functions given in the form of Eq (24) with the amplitudes of Eq. (25).
III. DIFFERENT SINGLE PARTICLE BASES
It is very common that the HFB wave functions |φi〉 are defined in terms of different single particle basis, with creation and annihilation operators a+k (i) and ak(i) that will carry indexes (i) to indicate the HFB wave function they belong to. Those bases are usually not complete and therefore they span different subspaces of the full Hilbert space. As a consequence, the formulas obtained above cannot be used because they implicitly rely on a common basis for the two HFB wave functions [9]. The strategy to overcome this problem is to find a bigger subspace en- compassing both subspace and use an orthogonal basis defined in that bigger subspace for the two HFB wave functions |φi〉. In an early consideration of this prob- lem [9] I used the whole Hilbert space as the common subspace. Another possibility, explored in this paper, is to consider the subspace generated by the union of the two subspaces. In this case, special care has to be taken with the resulting basis, union of the of the two original bases, as it can be redundant (i.e. it can contain linearly dependent vectors). Let me start by considering the two basis {a+k (0), k =
0, . . . , N(0)} and {a+k (1), k = 0, . . . , N(1)} which are de- fined in terms of single particle creation operators (typi- cally those of the harmonic oscillator basis) and with di- mensions N(0) andN(1), respectively, which do not neces- sarily coincide. It is worth introducing the set of creation

6 operators
A+µ =
{ a+k (0) µ = k k = 1, . . . , N(0) a+l (1) µ = l +N(0) l = 1, . . . , N(1)
embracing the two sets of creation operators of the bases. They satisfy the commutation relations {Aµ, A+ν } = Nµν and {Aµ, Aν} = {A+µ , A+ν } = 0 where the overlap matrix N (dimension (N(0) + N(1)) × (N(0) + N(1)) is given in terms of the rectangular matrix Tk k′ = {ak(0), a+k′(1)} = 〈−|ak(0)a+k′(1)|−〉 by the expression
N = (
I(0) T
T+ I(1)
) . (26)
The overlap matrix is hermitian, semi-positive definite and therefore can be diagonalized by a unitary transfor- mation D, i.e. N = DnD+where the diagonal matrix n of the eigenvalues is of dimension (N(0)+N(1))× (N(0)+ N(1)). In order to deal with the zero (or smaller than a given threshold) eigenvalues case (which signals the ap- pearance of linearly dependent basis states) it is conve- nient to introduce the notation
n =
( n¯ 0
0 ǫ
) .
where n¯ is a diagonal matrix with the eigenvalues dif- ferent from zero and ǫ is the diagonal matrix of dimen- sion Nǫcontaining those eigenvalues with value zero (or smaller than a numerical threshold). It is convenient in the ensuing developments to consider that the matrix ǫ is different from zero and therefore can be inverted. At the end of the calculations ǫ will be made to tend to zero to obtain the final result. Taking this regulariza- tion scheme into account, we can define the square root of the overlap matrix N 1/2 = Dn1/2D+ and its inverse N−1/2 = Dn−1/2D+ that are required to define the op- erators
Bµ = ∑ ν
N−1/2µν Aν
as well as the inverse relation Aµ = ∑
ν N 1/2µν Bν . The creation and annihilation operators B+µ and Bν satisfy canonical commutation relations {Bµ, B+ν } =(N−1/2NN−1/2)
µν = δµν . They are introduced to ex-
press the HFB wave functions of Eq. (24) in the standard way as
|φi〉 = exp   ∑ µµ′
1
2 N˜
(i) µµ′B
+ µ B
+ µ′
  |0〉
with the matrices of dimension (N(0) + N(1)) × (N(0) + N(1))
N˜ (i) = N 1/2+M˜ (i)N 1/2 ∗ (27)
given in terms of the extended matrices (also of dimen- sion (N(0) +N(1))× (N(0) +N(1)))
M˜ (0) =
( M (0) 0
0 0
) , M˜ (1) =
( 0 0
0 M (1)
) . (28)
As the operators B+µ and Bν satisfy canonical commu- tation relations, it is now possible to apply the standard formalism already developed in Ref [1] to write
〈φ0|φ1〉 = SN(0)+N(1)pf(M˜)
where the matrix M˜ entering the argument of the pfaffian is given in terms of the matrices defined in Eq (28) as
M˜ =
( N˜ (1) −I I −N˜ (0) ∗
) .
Using the results of appendix III concerning the eigenval- ues and eigenvectors of the norm matrix N (to simplify the notation in the following we consider the size of the two basis to be equal N(0) = N(1) = N) and using the definition of Eq (27) we have
N˜ (i) = D
( n 1/2 + 0
0 n 1/2 −
) D+M˜ (i)D∗
( n 1/2 + 0
0 n 1/2 −
) DT
= DN˜ (i) D D
T
which defines the matrices N˜ (i) D as
N˜ (i) D =
( n 1/2 + 0
0 n 1/2 −
) D+M˜ (i)D∗
( n 1/2 + 0
0 n 1/2 −
) .
Please note that these matrices are well defined when some of the eigenvalues of the norm overlap matrix N go to zero (i.e. some of the elements of the diagonal matrix n. are zero). Using the explicit form of the matrix D given in appendix III in terms of the matrices E and F entering the Singular Value Decomposition (SVD) of the
matrix T and defining the auxiliary matricesE± = n 1/2 ± E
and F± = n 1/2 ± F we get
N˜ (0) D =
1
2
( E+ 0
0 E−
)( M (0) −M (0) −M (0) M (0)
)( ET+ 0
0 ET−
)
and
N˜ (1) D =
1
2
( F+ 0
0 F−
)( M (1) M (1)
M (1) M (1)
)( FT+ 0
0 FT−
)
Finally, using known properties of the pfaffian we can
express the overlap in terms of the N˜ (i) D as
〈φ0|φ1〉 = (−1)Npf ( N˜
(1) D −I I −N˜ (0) ∗D
) (29)

7 which is the final expression for the overlap. To finish this section it is worth considering the limit
where the two bases are connected by means of an unitary transformation (see appendix III). In this case, the SVD
of T is trivial and we have n 1/2 + =
√ 2I, n
1/2 − = 0, E = I
and F = T . Using these values we have
N˜ (1) D =
( TM (1)T T 0
0 0
)
and
N˜ (0) D =
( M (0) 0
0 0
) .
As a consequence of the zero eigenvalues of the overlap
matrix, the matrices N˜ (i) D acquire a bipartite structure
where only the upper diagonal block is different from zero. In this case, the matrix in the pfaffian in Eq. (29) becomes a block matrix
 TM (1)T T 0 −I 0
0 0 0 −I I 0 −M (0) ∗ 0 0 I 0 0
  .
The pfaffian of this matrix can be simplified by using the results of appendix I and, after exchanging blocks 2 and 3, bring the matrix to block diagonal form. Once in block diagonal form, the pfaffian can be reduced to the product of the pfaffian of each of the diagonal blocks by using a simplified version of the results of appendix II for the pfaffian of a bipartite matrix. The final results is
〈φ0|φ1〉 = (−1)N(N+1)/2pf ( TM (1)T T −I
I −M (0)∗ )
as expected (see Ref [1]).
IV. CONCLUSIONS
We have analyzed the formula for the evaluation of the norm overlap of two different HFB wave functions in the situation where some of the occupancies of the quasipar- ticle levels are one and the standard approach leads to indeterminacies what have to be singled out in order to obtain a meaningful answer. In the case where there are fully empty single particle levels the overlap formula is well behaved but this limit is also addressed as a (sig- nificant) reduction in the computational burden can be obtained if the number of empty levels is large enough. Finally, the case where each of the two HFB wave func- tions are expressed in different single particle basis is ad- dressed and the formalism to compute the overlap in this situation is developed. A common basis, union of the other two is considered, and special attention has to be
paid to the problem of redundancy of the enlarged sub- space. The formulas given in this paper are a practical complement of the general one given in [1] and should be useful for a practical implementation of the calculation of the overlaps of two different HFB wave functions in the most general case.
ACKNOWLEDGMENTS
This work was supported by MICINN (Spain) under research grants No. FPA2009-08958, and No. FIS2009- 07277, as well as by Consolider-Ingenio 2010 Programs CPAN CSD2007-00042 and MULTIDARK CSD2009- 00064.
I. REORDERING OF MATRICES AND ITS
IMPACT IN THE PFAFFIAN
In some situations we will have to make use of eventual peculiarities of the block structure of the skew-symmetric matrix in order to simplify the final expression of the pfaffian. To do so, it is convenient to know how to reorder rows and columns of a matrix as well as the impact of such reordering in the pfaffian. A useful set of matrices is the one of the matrices E(i, j) that exchanges columns i and j of a matrix when multiplied to the right hand side of that matrix, i.e.( . . . ai−1 ai ai+1 . . . aj−1 aj aj+1 . . . ...
... ...
... ...
... ...
... ...
) E(i, j) =
( . . . ai−1 aj ai+1 . . . aj−1 ai aj+1 . . . ...
... ...
... ...
... ...
... ...
)
The matricesE(ij) are characterized by the matrix ele- ments
E(i, j)kl = δkl − δkiδli − δkjδlj + δkjδli + δkiδlj and are unit matrices where the elements i and j of the diagonal are set to zero and the elements i, j and j, i are set to one. Another useful set of matrices is the one of the S(i, j)
matrices such that, when multiplied to the right hand side of a matrix, moves the column j of the matrix to the position of column i (i < j) and then shifts column i to position i+ 1, column i+ 1 to position i+ 2, and so on, up to column j − 1 that is shifted to column j, i.e.( . . . ai−1 ai ai+1 ai+2 . . . aj−1 aj aj+1 . . . ...
... ...
... ...
... ...
... ...
...
) S(i, j) =
( . . . ai−1 aj ai ai+1 ai+2 . . . aj−1 aj+1 . . . ...
... ...
... ...
... ...
... ...
...
)

8 In terms of matrix elements they are given by
S(i, j)kl = δkl − j∑
s=i
δksδls +
j−1∑ s=i
δksδls+1 + δkjδli.
These matrices are unit matrices where the 1 in position i, i is shifted to position i, i+1 , the 1 in position i+1, i+1 is shifted to position i + 1, i + 2 and so on up the 1 in position j, j that is shifted to position j, i. The determi- nants of the two kind of matrices are easy to determine and they are given by
det(E(i, j)) = −1 and
det(S(i, j)) = (−1)j−i. The successive application of the matrices E(i+ k, j+ k) for k = 0 up to k = N defines a matrix
PN (i, j) = N−1∏ k=0
E(i + k, j + k)
that exchanges a set of N columns at once( . . . ai . . . ai+N−1 . . . aj . . . aj+N−1 . . . ...
... ...
... ...
... ...
... ...
) PN (i, j) =
( . . . aj . . . aj+N−1 . . . ai . . . ai+N−1 . . . ...
... ...
... ...
... ...
... ...
)
Applying the matrix PTN to the left of the matrix the corresponding exchange of rows is produced. As a con- sequence
PTN (i, j)
 
. . . ...
... ...
· · · AN (i, i) · · · AN (i, j) · · · ...
... . . .
... ...
· · · AN (j, i) · · · AN (j, j) · · · ...
... ...
. . .
  PN (i, j) =
 
. . . ...
... ...
· · · AN (j, j) · · · AN (j, i) · · · ...
... . . .
... ...
· · · AN (i, j) · · · AN (i, i) · · · ...
... ...
. . .
 
(1)
where AN (i, j) are sub-matrices of dimension N × N whose first element is located in the row i and column j of the matrix where the AN are embedded. The result obtained, together with that of the pfaffian of a bipar- tite matrix will be useful to reduce some of the pfaffians encountered in the main body of the paper. Obviously detPN = (−1)N . Unfortunately, this trick can not be
applied when the number of columns to be “exchanged” is not the same. In such case, we have to consider the more general shift operation S(i, j). To see how it works let us consider a matrix A were there are three groups of columns denoted by L, C, and R such that the first group goes from column i to columns i+NL− 1 (i.e. NL columns), the group C goes from column i + NL up to columns i +NL +NC − 1 (i.e. NC columns) and finally the group R from column i + NL + NC up to column i+NL+NC +NR− 1 (i.e. NR columns). Schematically, the columns of the matrix A could be represented as
A = (. . . |L|C|R| . . .) Now we want to exchange the group L with the group R; to do so the group R of columns is moved to the position occupied by L using the product of matrices
P (1) LR =
NR−1∏ k=0
S(i+ k, i+ k +NL +NC)
giving
(. . . |L|C|R| . . .)P (1)LR = (. . . |R|L|C| . . .) Now, the group of columns C is moved to the position of the group L by means of the following product of “shift” matrices
PLC =
NC−1∏ k=0
S(i+ k +NR, i+ k +NR +NL)
Using this matrix we obtain
(. . . |L|C|R| . . .)P (1)LRPLC = (. . . |R|L|C| . . .)PLC = (. . . |R|C|L| . . .)
The matrix exchanging the set of columns L with the
set R will be denoted PLR = P (1) LRPLC . By applying P
T LR
to the left of the matrix the set of rows L and R are exchanged. As a consequence
PTLR(i, j)
 
. . . ...
... ...
· · · ANL×NL(i, i) · · · ANL×NR(i, j) · · · ...
... . . .
... ...
· · · ANR×NL(j, i) · · · ANR×NR(j, j) · · · ...
... ...
. . .
  PLR(i, j) =
 
. . . ...
... ...
· · · ANR×NR(j, j) · · · ANR×NL(j, i) · · · ...
... . . .
... ...
· · · ANL×NR(i, j) · · · ANL×NL(i, i) · · · ...
... ...
. . .
 
where ANL×NL(i, i) is a sub-matrix withNL rows andNL columns starting at row i and column i, ANL×NR(i, j) is

9 a sub-matrix with NL rows and NR columns starting at row i and column j = i + NL + NC , and so on. This result generalizes the one of Eq (1) for sub-matrices of different sizes. The impact of such exchange of rows and columns in the pfaffian is the product of the de- terminants of the S matrices involved, i.e. det(PLR) =
det(P (1) LR) det(PLC). With det(P
(1) LR) = (−1)(NL+NC)NR
and det(PLC) = (−1)NLNC the total phase is det(PLR) = (−1)NLNR+NCNR+NLNC that reduces to (−1)N when NL = NR = N , a result that is independent of NC .
II. THE PFAFFIAN OF A BIPARTITE MATRIX
To derive some of the results obtained in this paper, it is often required to compute the pfaffian of a bipartite skew-symmetric matrix with the general structure
S =
( M Q
−QT N
) (1)
where M and N are skew-symmetric matrices and Q is a general rectangular matrix. Using Aitken’s formula it is possible to diagonalize the bipartite matrix using a congruence transformation(
I 0
QTM−1 I
)( M Q
−QT N
)( I −M−1Q 0 I
) =
( M 0
0 N +QTM−1Q
)
or the equivalent expression( I −QN−1 0 I
)( M Q
−QT N
)( I 0
N−1QT I
) =
( M +QN−1QT 0
0 N
)
to be used if M−1 does not exist. This is a very con- venient block diagonalization formula as it involves con- gruence transformations that allow to use the property pf(PTRP ) = det(P )pf(R) of the pfaffian to obtain the next two identities
pf(S) = pf(M)pf(N +QTM−1Q) (2)
= pf(M +QN−1QT )pf(N) (3)
III. THE OVERLAP MATRIX
The overlap matrix N has the bipartite structure
N = (
I(0) T
T+ I(1)
)
where the rectangular matrix T , with matrix elements Tk k′ = {ak(0), a+k′ (1)} = 〈0|ak(0)a+k′(1)|0〉, is the matrix
of the overlaps between the elements of the two basis considered. As in the body of the paper, N(0) and N(1) denote the dimensions of the each of the basis and it is assumed for definiteness that N(0) ≥ N(1).The matrices I(0) and I(1) stand for the identity matrices of dimen- sions N(0) and N(1), respectively. For the developments considered in this paper, the analysis of the spectral de- composition of the overlap matrix is required in order to handle properly the occurrence of very small or zero eigenvalues of the overlap. The analysis is based on the Singular Value Decomposition (SVD) [28] of the matrix T
T = E+∆F (1)
where E and F are square unitary matrices of dimen- sions N(0) × N(0) and N(1) × N(1) respectively and ∆ is a rectangular matrix of dimension N(0) × N(1) with the “diagonal structure”
∆ =
( ∆¯
0
)
where ∆¯ is a real and positive square diagonal matrix with dimension N(1)×N(1). It is convenient to introduce a rectangular “identity matrix” I(01) of dimension N(0) × N(1) with a structure similar to the matrix ∆, namely
I(01) =
( I(1)
0
)
and such that ∆ = I(01)∆¯. This matrix also has the
property I+(01)I(01) = I(1).
Using the SVD of T defined in Eq. (1) we can finally write
N = D¯ (
I(0) ∆
∆+ I(1)
) D¯+
where
D¯ =
( E+ 0
0 F+
) .
The matrix in the middle can be easily brought to diag- onal form
D¯(0)+
( I(0) ∆
∆+ I(1)
) D¯(0) =
( 2I(0) + I(01)(∆¯− I(1))I+(01) 0
0 I(1) − ∆¯
)
by means of the D¯(0) transformation
D¯(0) = 1√ 2
( I(0) −I(01) I + (01) I(1)
)

10
Finally, introducing the matrix D = D¯D¯(0) we obtain the complete diagonalization of the overlap matrix
N = D (
2I(0) + I(01)(∆¯− I(1))I+(01) 0 0 I(1) − ∆¯
) D+.
(2) The semi-positive character of the matrix N implies that all its eigenvalues are positive or zero and therefore (lower block) the elements of the diagonal matrix ∆¯ cannot ex- ceed one. In the case of bases with equal dimensions N(0) = N(1) the above result becomes
N = D (
I+ ∆¯ 0
0 I− ∆¯
) D+ = D
( n+ 0
0 n−
) D+ (3)
where the diagonal matrices n± = I± ∆¯ have been intro- duced. As the matrix ∆¯ is also positive semidefinite it is clear that the zero norm eigenvalues of N are associated to values of ∆¯ equal to one. It is helpful to consider the special case when T is an
unitary matrix. It corresponds to the common situation
where the two bases are connected through a unitary transformation as, for instance, the ones associated to symmetry operations like rotations in real space, etc act- ing on closed basis. For unitary T , the SVD factors can be chosen as E+ = T , ∆ = I and F = I in the formulas above so that
D = 1√ 2
( T −T I I
) (4)
and
N = D ( 2I 0
0 0
) D+ (5)
with N -fold degenerate eigenvalues 2 and 0.
To finish this appendix, just mention that a similar (but less general) treatment of the norm overlap was con- sidered in Ref [29] but without resorting to the powerful concept of the SVD.
[1] L.M. Robledo, Phys. Rev. C79, (2009). [2] John R. Klauder, and Bo-Sture Skagerstam, Coherent
states: applications in physics and mathematical physics (World Scientific, Singapore, 1985).
[3] F.A. Berezin, The Method of Second Quantization (Aca- demic Press, New York, 1966).
[4] J.-P. Blaizot, and G. Ripka, Quantum Theory of Finite Systems (MIT Press, Cambridge, Massachusetts, and London, England, 1985).
[5] J. W. Negele, and H. Orland, Quantum Many Particle Systems (Addison Wesley, Redwood City, 1988).
[6] Y. Ohnuki, and T. Kashiwa, Prog. Theor. Phys. 60, 548 (1978).
[7] E.R. Caianiello, Combinatorics and renormalization in Quantum Field Theory (W.A. Benjamin, Massachusetts, 1973).
[8] P. Ring, and P. Schuck, The Nuclear Many Body problem (Springer-Verlag, Berlin, 1980).
[9] L.M. Robledo, Phys. Rev. C50, 2874 (1994). [10] J.M. Yao, J. Meng, P. Ring, and D. Pena Arteaga, Phys.
Rev. C79, 044312 (2009). [11] P. Bonche, J. Dobaczewski, H. Flocard, P.-H. Heenen,
and J. Meyer, Nucl. Phys. A510, 466 (1990). [12] Q. Haider, and D. Gogny, J. Phys. G18, 993 (1992). [13] M. Bajdich, L. Mitas, and L.K. Wagner, Phys. Rev B77,
115112 (2008). [14] C.K. Thomas, and A. A. Middleton, Phys. Rev E80,
046708 (2009). [15] J.-M. Stéphan, S. Furukawa, G. Misguich, and V.
Pasquier, Phys. Rev B80, 184421 (2009). [16] C. Gonzalez-Ballestero, L.M. Robledo, and G.F. Bertsch,
submitted to Comp. Phys. Comm, arXiv:1012.5022. [17] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod.
Phys. 75, 121 (2003). [18] R. R. Rodríguez-Guzmán, J. L. Egido, and L. M. Rob-
ledo, Nucl. Phys. A709, 201 (2002). [19] S. Perez-Martin, and L.M. Robledo, Phys. Rev. C76,
064314 (2007). [20] K. W. Schmid, Prog. Part. Nucl. Phys. 52, 565 (2004). [21] M. Bender, and P.-H. Heenen, Phys. Rev. C 78, 024309
(2008). [22] M. Oi, and N. Tajima, Phys. Lett. B606, 43 (2005). [23] F. Dönau, Phys. Rev. C58, 872 (1998). [24] M. Anguiano, J. L. Egido, and L. M. Robledo, Nucl.
Phys. A696, 467 (2001). [25] D. Lacroix, T. Duguet, and M. Bender, Phys. Rev. C79,
044318 (2009). [26] William H. Press, Saul A. Teukolsky, William T. Vet-
terling, and Brian P. Flannery, Numerical Recipes: The Art of Scientific Computing (3rd Edition), (Cambridge, 2007).
[27] J.P. Blaizot, and H. Orland, Phys. Rev. C24, 1740 (1981).
[28] G. H. Golub and C. F. Van Loan, Matrix Computations (Johns Hopkins University Press, Baltimore, 1996).
[29] J.F. Berger, and D. Gogny, Nucl. Phys. A333, 302 (1980).