Calculating the norm matrix to solve the A-body Schr\"odinger equation within a set of non-orthogonal many-body states

TL;DR

This paper introduces a new, general method for calculating the norm matrix in many-body quantum physics, improving the solution of the A-body Schrödinger equation with non-orthogonal states, exemplified on Bogoliubov states.

Contribution

It presents an alternative, broadly applicable method for computing the norm matrix in many-body problems involving non-orthogonal states, addressing longstanding phase ambiguity issues.

Findings

Method successfully applied to Bogoliubov states

Numerical tests on a toy model demonstrate effectiveness

Provides a more general approach than previous solutions

Abstract

There are efficient many-body methods, such as the (symmetry-restored) generator coordinate method in nuclear physics, that formulate the A-body Schr\"odinger equation within a set of nonorthogonal many-body states. Solving the corresponding secular equation requires the evaluation of the norm matrix and thus the capacity to compute its entries consistently and without any phase ambiguity. This is not always a trivial task, e.g. it remained a long-standing problem for methods based on general Bogoliubov product states. While a solution to this problem was found recently in Ref. [L. M. Robledo, Phys. Rev. C79, 021302 (2009)], the present work introduces an alternative method that can be generically applied to other classes of states of interest in many-body physics. The method is presently exemplified in the case of Bogoliubov states and numerically illustrated on the basis of a toy model.

Equations

$$\mathcal{M}\equiv\big{\{}|\Phi_{k}\rangle,k=1\ldots,N\big{\}}\,.$$

$$|{\Psi }_n⟩\equiv \sum _{k=1}^N{f}_{nk}|{\Phi }_k⟩,$$

$$\mathcal{H}{\mathfrak{f}}_n=E_n\mathcal{N}{\mathfrak{f}}_n,$$

$${\mathcal{N}}_{kl}$$

$${\mathcal{H}}_{kl}$$

$$|{\Phi }_l⟩\equiv e^{iS[k,l]}|{\Phi }_k⟩,$$

$$\mathcal{M}[k,l]\equiv\big{\{}|\Phi_{kl}(\theta)\rangle\equiv e^{i\theta S[k,l]}|\Phi_{k}\rangle\,,\,\theta\in[0,1]\big{\}}\,,$$

$$n_{kl}[⟨\Theta |,\theta ]\equiv \frac{⟨\Theta |{\Phi }_{kl}(\theta )⟩}{⟨\Theta |{\Phi }_k⟩}$$

$$\begin{array}{rl}\frac{d}{d\theta }{n}_{kl}[⟨\Theta |,\theta ]& =i\frac{⟨\Theta |S[k,l]|{\Phi }_{kl}(\theta )⟩}{⟨\Theta |{\Phi }_k⟩}.\end{array}$$

$$n_{kl}[⟨\Theta |,\theta ]=e^{i{\int }_0^{\theta }d\varphi s_{kl}[⟨\Theta |,\varphi ]},$$

$$s_{kl}[⟨\Theta |,\theta ]\equiv \frac{⟨\Theta |S[k,l]|{\Phi }_{kl}(\theta )⟩}{⟨\Theta |{\Phi }_{kl}(\theta )⟩},$$

$$\text{Arg}(⟨{\Phi }_1|{\Phi }_1⟩)=\text{Arg}(⟨{\Phi }_1|{\Phi }_2⟩)\dots =\text{Arg}(⟨{\Phi }_1|{\Phi }_N⟩)=0,$$

$$\frac{{\mathcal{N}}_{1l}}{{\mathcal{N}}_{11}}=e^{i{\int }_0^{1}d\varphi s_{1l}[⟨{\Phi }_1|,\varphi ]},$$

$$\Re e{\int }_0^{1}d\varphi s_{1l}[⟨{\Phi }_1|,\varphi ]=0,$$

$$\frac{{\mathcal{N}}_{ml}}{{\mathcal{N}}_{m1}}=e^{i{\int }_0^{1}d\varphi s_{1l}[⟨{\Phi }_m|,\varphi ]},$$

$${\left(\begin{array}{c}\beta \\ {\beta }^{†}\end{array}\right)}_{[k]}$$

$$S[k,l]$$

$$+\frac{1}{2}{\left({\beta }^{†}\beta \right)}_{[m]}{\left(\begin{array}{cc}{S}^{11}[k,l]& S^{20}[k,l]\\ -S^{02}[k,l]& -S^{11\ast }[k,l]\end{array}\right)}_{[m]}{\left(\begin{array}{c}\beta \\ {\beta }^{†}\end{array}\right)}_{[m]}$$

$${\left(\begin{array}{cc}{S}^{11}[k,l]& S^{20}[k,l]\\ -S^{02}[k,l]& -S^{11\ast }[k,l]\end{array}\right)}_{[k]}$$

$$\frac{{\mathcal{N}}_{1l}}{{\mathcal{N}}_{11}}=e^{-\Im m\frac{1}{2}{\int }_0^{1}d\varphi \mathrm{Tr}\left(S^{02}[1,l{]}_{[1]}{R}_{1l}^{--}[⟨{\Phi }_1|,\varphi ]\right)},$$

$${\left(R_{1l}^{--}[⟨{\Phi }_1|,\varphi ]\right)}_{k_1{k}_2}\equiv \frac{⟨{\Phi }_1|{\beta }_{k_1}^{[1]}{\beta }_{k_2}^{[1]}|{\Phi }_{1l}(\theta )⟩}{⟨{\Phi }_1|{\Phi }_{1l}(\theta )⟩},$$

$$\mathcal{N}\equiv \left(\begin{array}{ccc}⟨{\Phi }_1|{\Phi }_1⟩& ⟨{\Phi }_1|{\Phi }_2⟩& ⟨{\Phi }_1|{\Phi }_3⟩\\ ⟨{\Phi }_2|{\Phi }_1⟩& ⟨{\Phi }_2|{\Phi }_2⟩& ⟨{\Phi }_2|{\Phi }_3⟩\\ ⟨{\Phi }_3|{\Phi }_1⟩& ⟨{\Phi }_3|{\Phi }_2⟩& ⟨{\Phi }_3|{\Phi }_3⟩\end{array}\right).$$

$${\mathcal{W}}_{[k]}$$

$${\bar{U}}_{[k]}(p,\bar{p})$$

$${\bar{V}}_{[k]}(p,\bar{p})$$

$${\left(\begin{array}{c}\tilde{\beta }\\ {\tilde{\beta }}^{†}\end{array}\right)}_{[k]}$$

Full text

Introduction. The A-body Schrödinger equation is most often (approximately) solved by representing it on a (truncated) orthonormal basis of the A-body Hilbert space $\mathcal{H}_{A}$. Some approaches, however, represent the A-body Schrödinger equation on a finite set of non-orthogonal states of ${\cal H}_{A}$. At the price of losing the orthogonality, it permits to select fewer states on the basis of their expected physical relevance. The set may even exceed ${\cal H}_{A}$ by employing states mixing vectors belonging to Hilbert spaces associated with different particle numbers, i.e. states that are genuine vectors of Fock space ${\cal F}$. This is for instance the case of the generator coordinate method (GCM) in use in nuclear physics Ring and Schuck (1980); Bender et al. (2003); Duguet (2014); Egido (2016) and of symmetry restoration calculations employed in nuclear physics Ring and Schuck (1980); Bender et al. (2003); Duguet (2014); Egido (2016) and quantum chemistry Jimenez-Hoyos et al. (2012a, b); Samanta et al. (2012); Rivero et al. (2013). In both methods, the Hamilton operator is diagonalized in the finite-dimensional vector space spanned by a set of non-orthogonal Bogoliubov product states. The secular equation to be solved requires the evaluation of the norm matrix constructed from overlaps between all members of the set. Furthermore, the recently developped particle-number-restored Bogoliubov coupled-cluster and particle-number-restored many-body perturbation theories Duguet and Signoracci (2017) also build on a manifold of non-orthogonal Bogoliubov states. The norm kernels at play are more general as they explicitly incorporate many-body correlations and reduce to the mere overlap between two non-orthogonal Bogoliubov states whenever such correlations are omitted. Last but not least, the efficient computation of overlaps constitutes a key element of quantum monte carlo (QMC) approaches, especially when they rely on more elaborate walkers and/or trial states than Slater determinants Guerrero et al. (1998); Roggero et al. (2013); Puddu (2003); Juillet et al. (2017).

The closed-form evaluation of the overlap between many-body states constitutes a long-standing problem. For instance, while the overlap between two non-orthogonal A-body Slater determinants has long been known to be computable as a determinant Blaizot and Ripka (1986), the efficient and unambiguous111In general, the well-known Onishi formula Onishi and Yoshida (1966) can provide the norm of the overlap but not its complex phase. calculation of the overlap between two arbitrary Bogoliubov quasiparticle states has only become available recently as the Pfaffian of a skew-symmetric matrix Robledo (2009, 2011); Avez and Bender (2012).

In this context, we presently propose a closed-form expression for the overlap between arbitrary many-body states and implement it in the context of GCM calculations based on a set of Bogoliubov states. All pertinent technical details related to the latter case can be found in Ref. Bally and Duguet (2017). The application to other classes of states in use in many-body physics remains to be investigated in the future.

Objective. We consider the situation where the static A-body Schrödinger equation is represented on a set of $N$ (a priori non-orthogonal) many-body states

[TABLE]

The eigenvectors of the Hamiltonian $H$ are expanded in terms of the states of $\mathcal{M}$ according to

[TABLE]

where $f_{nk}$ are complex numbers. The corresponding secular equation takes the form of a generalized eigenvalue problem

[TABLE]

where $\mathfrak{f}_{n}$ is the weight column matrix, i.e. $(\mathfrak{f}_{n})_{k}\equiv f_{nk}$, and where the norm and Hamiltonian hermitian matrices are made out of $N(N+1)/2$ independent elements

[TABLE]

respectively. Solving Eq. (3) gives access to energies $E_{n}$ and weigths $\mathfrak{f}_{n}$ for states $|\Psi_{n}\rangle$. Typically, such a generalized eigenvalue problem is addressed by first diagonalizing the norm matrix ${\cal N}$ before solving a standard eigenvalue problem in the (possibly smaller) orthonormal basis built from the norm eigenvectors with non-zero eigenvalues. A key step of the many-body calculation is thus the consistent computation of the $N(N+1)/2$ independent entries of the norm matrix, including their complex phase.

Master formulae. Given an arbitrary pair of states $(|\Phi_{k}\rangle,|\Phi_{l}\rangle)$ belonging to $\mathcal{M}$, we assume that a unitary transformation linking both states is either known or can be extracted under the form

[TABLE]

where $S[k,l]$ is a hermitian operator. Based on this sole hypothesis, an auxiliary set connecting both states is introduced as

[TABLE]

such that $|\Phi_{kl}(0)\rangle=|\Phi_{k}\rangle$ and $|\Phi_{kl}(1)\rangle=|\Phi_{l}\rangle$. Considering an arbitrary many-body bra $\langle\Theta|$, the quantity

[TABLE]

satisfying $n_{kl}[\langle\Theta|,0]=1$ is defined along the auxiliary manifold. Differentiating it with respect to $\theta$ leads to

[TABLE]

Assuming that $n_{kl}[\langle\Theta|,\theta]\neq 0$ along $\mathcal{M}[k,l]$, one divides both sides of Eq. (8) by it and integrates the corresponding first-order differential equation between [math] and $\theta$ to obtain

[TABLE]

where the off-diagonal linked-connected kernel Duguet (2015); Duguet and Signoracci (2017) of the operator $S[k,l]$ defined along the manifold $\mathcal{M}[k,l]$

[TABLE]

is unambiguous in the sense that it is independent of the relative phase between $|\Phi_{k}\rangle$ and $|\Phi_{l}\rangle$. Equations (9)-(10) constitute the master formulae repeatedly used below to access the norm matrix (Eq. (4a)).

Phase convention. The normalized states belonging to $\mathcal{M}$ are all individually defined up to a phase that must not influence the computation of observables. This freedom must be explicitly controlled such that the $N(N+1)/2$ independent entries to the norm matrix are computed consistently. This relates to fixing the relative phases of the states in a synchronized fashion, which effectively impacts the definition of each operator $S[k,l]$. This can be done by specifying the phase each member of $\mathcal{M}$ entertains with a common known state of reference generically denoted as $|\bar{\Phi}\rangle$. Among many possibilities, a natural and practical choice consists of picking this pivot state within $\mathcal{M}$ itself and requiring that all states of $\mathcal{M}$ have the same phase relative to it. Accordingly, we choose $|\bar{\Phi}\rangle\equiv|\Phi_{1}\rangle$, although any other state of $\mathcal{M}$ would be equally appropriate, and impose that

[TABLE]

given that $\langle\Phi_{1}|\Phi_{1}\rangle$ is real.

Algorithm. The evaluation of the $N(N+1)/2$ independent entries to the norm matrix follows three successive steps

The $N$ diagonal elements ${\cal N}_{kk}$ are trivially obtained by normalizing all members of the set, i.e. by imposing that $\langle\Phi_{k}|\Phi_{k}\rangle=1$ for $k=1,\ldots,N$. 2. 2.

The $N-1$ remaining elements of the first row are computed on the basis of the $N-1$ operators $S[1,l]$. Applying Eq. (9) for $k=1$, $\langle\Theta|\equiv\langle\Phi_{1}|$ and $\theta=1$, one obtains

[TABLE]

where $s_{1l}[\langle\Phi_{1}|,\theta]$ runs over the manifold $\mathcal{M}[1,l]$ and is defined for the bra $\langle\Phi_{1}|$. The phase convention stated above constrains the pure number entering the definition of $S[1,l]$ according to

[TABLE]

such that ${\cal N}_{1l}$ is effectively real. 3. 3.

Applying Eq. (9) again for $k=1$ but now setting $\langle\Theta|\equiv\langle\Phi_{m}|$, $1<m<l\leq N$, the remaining $(N-1)(N-2)/2$ independent overlaps ${\cal N}_{ml}$ are obtained consistently for $\theta=1$ via

[TABLE]

where the off-diagonal kernel of $S[1,l]$ still runs over the manifold $\mathcal{M}[1,l]$ but now involves the bra $\langle\Phi_{m}|$ rather than $\langle\Phi_{1}|$. Since ${\cal N}_{m1}={\cal N}_{1m}$ is among the $N-1$ overlaps computed in step 2, Eq. (13) completes the norm matrix. While the diagonal, the first row and the first column of the norm matrix are real, the remaining entries are a priori complex.

Each entry ${\cal N}_{ml}$, including its complex phase, has been powerfully expressed in terms of the integral of an off-diagonal kernel of the operator $S[1,l]$ along the auxiliary manifold connecting $|\Phi_{l}\rangle$ to $|\Phi_{1}\rangle$. In essence, the rationale of the method is to commute the computation of a pure overlap into the computation of the linked-connected kernel of an operator that has no phase ambiguity. Of course, the usefulness of the method relies on our ability to compute such an operator kernel, which itself depends on the nature of the many-body states making up $\mathcal{M}$. More specifically, the character of $|\Phi_{1}\rangle$ and $|\Phi_{l}\rangle$ determines the nature of the operator $S[1,l]$ driving their unitary connection along with our capacity to extract it and compute efficiently its off-diagonal kernels $s_{1l}[\langle\Phi_{m}|,\theta]$.

Application. We apply the above scheme on the basis of a set made out of arbitrary Bogoliubov states Ring and Schuck (1980); Blaizot and Ripka (1986). This particular choice is characteristic of state-of-the-art (symmetry-restored) GCM calculations in nuclear physics.

Each member $|\Phi_{k}\rangle$ of $\mathcal{M}$ is defined as a vacuum, i.e. $\forall\,\mu,\beta^{[k]}_{\mu}|\Phi_{k}\rangle=0$, of the set of quasi-particle operators $\{\beta^{[k]}_{\mu};\beta^{[k]\dagger}_{\mu}\}$. These creation and annihilation operators relate to those defining a basis $\{c_{p};c^{\dagger}_{p}\}$ of the $n$-dimensional one-body hilbert space ${\cal H}_{1}$ via a linear Bogoliubov transformation

[TABLE]

where the unitarity of the $2n\times 2n$ matrix ${\cal W}_{[k]}$ ensures the fermionic character of the quasi-particle operators. The above procedure actually defines $|\Phi_{k}\rangle$ only up to a complex phase Hara and Iwasaki (1979); Blaizot and Ripka (1986).

Given $|\Phi_{k}\rangle$ and $|\Phi_{l}\rangle$, the operator $S[k,l]$ parameterizing their unitary connection is a generic hermitian one-body operator reading, in the quasi-particle basis of an arbitrary third state $|\Phi_{m}\rangle$, as

[TABLE]

where $S^{00}[k,l]_{[m]}$ is a real number, $S^{11}[k,l]_{[m]}$ is a hermitian matrix whereas $S^{20}[k,l]_{[m]}$ and $S^{02}[k,l]_{[m]}$ denote skew-symmetric matrices satisfying $S^{02}[k,l]_{[m]}=(S^{20}[k,l]_{[m]})^{\ast}$. The non-trivial part of the operator $S[k,l]$ is uniquely extracted (Bally and Duguet, 2017) from the sole knowledge of the Bogoliubov transformations ${\cal W}_{[k]}$ and ${\cal W}_{[l]}$, i.e. it reads in the quasi-particle basis of $|\Phi_{k}\rangle$ as

[TABLE]

Following the algorithm layed down above, only the $N-1$ operators $S[1,l]$ are effectively needed. Their off-diagonal kernels $s_{1l}[\langle\Phi_{m}|,\theta]$ along the auxiliary manifold $\mathcal{M}[1,l]$ can be unambiguously computed on the basis of the off-diagonal Wick theorem Balian and Brézin (1969). Fixing the constant $S^{00}[1,l]_{[1]}$ entering $S[1,l]$ via the application of Eq. (12), the entries on the first row of the norm matrix (Eq.(11)) are thus obtained under the form

[TABLE]

where the elementary off-diagonal contraction defined as

[TABLE]

is computable from ${\cal W}_{[1]}$ and ${\cal W}_{[l]}$ (Bally and Duguet, 2017). Eventually, the norm matrix can be completed via a similar specification of Eq. (13) that requires the introduction of the elementary contraction $R^{--}_{1l}[\langle\Phi_{m}|,\phi]$. The Onishi formula Onishi and Yoshida (1966) is recovered for each entry by taking the norm of the corresponding expression Bally and Duguet (2017).

The overlap between two arbitrary Bogoliubov product states can be alternatively computed as the Pfaffian of a skew-symmetric matrix Robledo (2009, 2011); Avez and Bender (2012). This result relies on the Thouless representation Thouless (1960) of the Bogoliubov states and, as such, relates to a phase convention, i.e. $\forall\,k,\text{Arg}(\langle 0|\Phi_{k}\rangle)=1$, that differs from the one presently used. This is interesting to demonstrate that, while impacting individual overlaps, the overall phase convention does not influence output observables. All is needed is an internally consistent computation of the entries to the norm and Hamiltonian matrices.

Toy model. The method is numerically implemented on the basis of a set \mathcal{M}\equiv\big{\{}|\Phi_{1}\rangle,|\Phi_{2}\rangle,|\Phi_{3}\rangle\big{\}} of three different Bogoliubov states. The associated norm matrix reads as

[TABLE]

The Bogoliubov transformation associated with $|\Phi_{k}\rangle$ presently reads as

[TABLE]

where $L$ denotes a random $n\times n$ complex unitary matrix transforming the ten-dimensional, i.e. $n=10$, basis of ${\cal H}_{1}$ made out of five doubly-degenerated single-particle levels. Following ${\cal L}_{[k]}$, $\bar{{\cal W}}_{[k]}$ is a BCS transformation characterized by the set of real $2\times 2$ blocks of the form

[TABLE]

where $\bar{p}$ denotes the conjugated partner of $p$ and where $u^{2}_{p}[k]+v^{2}_{p}[k]=1$. The BCS occupations $v^{2}_{p}[k]$, $p=1,\ldots 5$, are decreasingly chosen in the interval $]0,1[$ to mimic a realistic fully paired system. Even though, for the sake of simplicity, fully paired states are presently considered, the validity of the method has been checked for Bogoliubov states containing an even or odd number of fully occupied single-particle states.

Results. The upper panel of Fig. 1 displays in the complex plane individual overlaps making up the norm matrix. Squares represent the results obtained from the Pfaffian method whereas circles denote those computed from the presently proposed method. Thick lines provide the auxiliary pathes followed from one overlap to the other, starting from $\langle\Phi_{1}|\Phi_{1}\rangle=1$. Consistently with the scheme exposed above, all the overlaps involving the pivot state $|\Phi_{1}\rangle$ are real, which is not the case for the Pfaffian method. Interestingly, the complex conjugate values $\langle\Phi_{2}|\Phi_{3}\rangle$ and $\langle\Phi_{3}|\Phi_{2}\rangle$ are consistently obtained by going through $|\Phi_{2}\rangle$ or $|\Phi_{3}\rangle$ first. Three circles help visualize that, while individual overlaps differ in both methods because of the distinct phase conventions used, they only do so by a complex phase. Eventually, the lower panel demonstrates that the eigenvalues of the norm matrix obtained from both methods are identical, thus showing the consistency of both calculations and the independence of the result on the phase convention used.

The derivation of Eq. (9) relied on the hypothesis that $n_{kl}[\langle\Theta|,\theta]\neq 0$ along $\mathcal{M}[k,l]$. In the present application, it is clear from Fig. 1 that this hypothesis is indeed fulfilled for all the overlaps involved. There however exists situations, e.g. global-gauge symmetry restoration calculations, in which it is not the case Bally and Duguet (2017). To overcome this apparent difficulty, one can perform an extra trivial Bogoliubov transformation of the quasi-particle operators of, e.g. $|\Phi_{k}\rangle$, among themselves

[TABLE]

where $K$ is a $n\times n$ unitary matrix. Maintaining the constraint from the phase convention, such a trivial Bogoliubov transformation modifies non-trivially the operator $S[k,l]$ and the auxiliary manifold linking both states without changing their overlap Bally and Duguet (2017). Figure 2 illustrates this powerful flexibility of the method that can be used to bypass rare problems associated with potential zeros of the overlap along the path linking both states.

Conclusions. The present paper proposes a powerful method to compute the overlap between many-body states belonging to a set $\mathcal{M}$ used to represent the A-body Schrödinger equation. Solving the corresponding secular equation requires the evaluation of the norm matrix and thus the capacity to compute its entries consistently and without any phase ambiguity. This is not always a trivial task, e.g. it remained a long-standing problem for methods based on general Bogoliubov product states in used in nuclear physics. While a solution to this problem was found recently Robledo (2009, 2011); Avez and Bender (2012), the presently proposed method provides an alternative that can be generically applied to other classes of states of interest in many-body physics.

The overlap of any two states belonging to $\mathcal{M}$ is given as the exponential of the integral of the off-diagonal linked-connected kernel of an operator along an auxiliary continuous set joining both states. Such a linked-connected kernel is free from any phase ambiguity. The operator in question, which needs to be known or extracted during the procedure, is the hermitian generator of a unitary transformation connecting both states.

In the present paper, the algebra is specified for sets made out of Bogoliubov states and numerically illustrated on the basis of a toy model. In this context, the overlap between arbitrary Bogoliubov states is computed, without any phase ambiguity, via elementary linear algebra operations. The numerical application nicely demonstrates the versatility of the method.