Bivariational Principle for an Antisymmetrized Product of Nonorthogonal Geminals Appropriate for Strong Electron Correlation

Most computational methods for the electronic structure theory of molecules and solids are based on the orbital (for molecules) or band (for solids) picture helgaker_book ; raghavachari ; headgordon ; helgaker . In this picture, the ground-state wavefunction is approximated by a Slater determinant of the single-electron eigenfunctions (orbitals) from a mean-field model (e.g., Hartree-Fock, Kohn-Sham). The electrons are assumed to move independently, with each electron feeling only the average repulsion from the other electrons in the system. In many systems this is a reasonable first approximation and a single Slater determinant wavefunction is a good starting point to develop the wavefunction. The orbital picture is particularly reliable for equilibrium thermodynamic properties of organic molecules but is much less accurate for chemical transition states, electronic excited states, and inorganic substances containing $d$-block and $f$-block elements.

Strongly-correlated systems require many Slater determinants for even a qualitatively correct physical description.garnet2 ; anisimov This is a clear indication that Slater determinants are not an efficient basis and the mean-field behaviour is not independent electrons. Conventional density-based and wavefunction-based methods are usually unreliable for strong correlation, and the methods that are appropriate are usually computationally expensive. Large molecules and complex materials with thousands of valence electrons can be routinely modelled with the orbital picture, but no such tools exist when the orbital picture fails. Our ultimate goal is to develop new models for strongly correlated systems with hundreds of valence electrons.

Because modelling strongly correlated substances is difficult, condensed-matter physicists usually model these systems by introducing model Hamiltonians that capture the key qualitative features of the system in question. We have recently shown how these model Hamiltonians can also be used to provide quantitative predictions for real physical systems johnson ; limacher ; kasia_hubbard ; pawel_geminal . The basic idea is to (1) find a model Hamiltonian that reproduces the key qualitative features of the system of interest, (2) use the wavefunction-form of that model Hamiltonian to model the substance.

In our previous workjohnson:2020 ; fecteau:2020 ; fecteau:2021 ; johnson:2021 ; moisset:2022 ; fecteau:2022 we have chosen wavefunctions that are eigenvectors of a model Hamiltonian,

$\displaystyle\hat{H}_{model}(\bm{\eta})\ket{\Psi_{model}(\bm{\eta})}$

$\displaystyle=E_{model}\ket{\Psi_{model}(\bm{\eta})}$

(1)

and minimized the expectation value of our target Hamiltonian with respect to the parameters, $\bm{\eta}$, in the model Hamiltonian,

$\displaystyle E_{GS}\approx\min_{\bm{\eta}}\braket{\Psi_{model}(\bm{\eta})}{\hat{H}}{\Psi_{model}(\bm{\eta})}.$

(2)

Because these model Hamiltonians are exactly solved by the algebraic Bethe Ansatz (ABA) korepin_book , we refer to the eigenfunctions of the model Hamiltonian as on-shell Bethe vectors. A variational method based on on-shell Bethe vectors has been implemented previously and shown to be quite accurate provided the correct on-shell vector is chosen. The relevant details are summarized briefly in section III.

The requirement that the wavefunction be an eigenfunction of the model Hamiltonian may be relaxedjohnson ; limacher ; moisset:2022 yielding off-shell Bethe vectors. Variational optimization of off-shell Bethe vectors is much more difficult (and perhaps often impossible), so our previous work used the projected Schrödinger equation to establish a system of nonlinear equations that can be solved for the parameters in the off-shell Bethe vectors,

$\displaystyle\braket{\Phi}{\hat{H}}{\Psi(\bm{\eta})}=E\braket{\Phi}{\Psi(\bm{\eta})}.$

(3)

The primary goal of this paper is to present a bivariational principle for off-shell Bethe vectors; this is presented in section IV.1. The basic idea is to rewrite the variational principle as,

$\displaystyle E_{GS}=\min_{\begin{subarray}{c}\braket{\tilde{\Psi}(\bm{\eta})}{\Psi(\bm{\eta})}=1\\ \ket{\Psi(\bm{\eta})}=\ket{\tilde{\Psi}(\bm{\eta})}\end{subarray}}\braket{\tilde{\Psi}(\bm{\eta})}{\hat{H}}{\Psi(\bm{\eta})}.$

(4)

That is to say, we take two states $\ket{\Psi(\bm{\eta})}$ and $\ket{\tilde{\Psi}(\bm{\eta})}$, and minimize (4) while attempting to ensure that the two states are the same. This will be explained in section V.2.

This approach is applicable to any ABA solvable model Hamiltonian, though we focus for the moment on the Richardson-Gaudin family of Hamiltonians gaudin2 ; rich1 ; rich2 ; rich3 ; dukelsky3 . The necessary background material is presented in section II. After discussing the on-shell and off-shell formulations, we discuss the extension to nonorthogonal orbitals in section IV.3. The bivariational principle for off-shell RG states and a different approach to the antisymmetrized geminal power wavefunction are presented in sections V and VI, respectively. Notes on the computational algorithm are presented in section VII, and we discuss possible generalizations in section VIII.

The systems we wish to solve are described by the Coulomb Hamiltonian,

$\displaystyle\hat{H}_{C}=\sum_{ij}h_{ij}\sum_{\sigma}a^{\dagger}_{i\sigma}a_{j\sigma}+\frac{1}{2}\sum_{ijkl}V_{ijkl}\sum_{\sigma\tau}a^{\dagger}_{i\sigma}a^{\dagger}_{j\tau}a_{l\tau}a_{k\sigma},$

(5)

where $a^{\dagger}_{i\sigma}$($a_{i\sigma}$) are the operators for creating (destroying) an electron in spatial orbital $i$ with spin projection $\sigma$. In Eq. (5), the integrals are expressed in physicists’ notation

	$\displaystyle h_{ij}$	$\displaystyle=\int d\mathbf{r}\phi^{*}_{i}(\mathbf{r})\left(-\frac{1}{2}\nabla^{2}-\sum_{I}\frac{Z_{I}}{|\mathbf{r}-\mathbf{R}_{I}|}\right)\phi_{j}(\mathbf{r})$		(6)
	$\displaystyle V_{ijkl}$	$\displaystyle=\int d\mathbf{r}_{1}d\mathbf{r}_{2}\frac{\phi^{*}_{i}(\mathbf{r}_{1})\phi^{*}_{j}(\mathbf{r}_{2})\phi_{k}(\mathbf{r}_{1})\phi_{l}(\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$		(7)

With a wavefunction ansatz $\ket{\Psi}$ the energy is computed,

$\displaystyle E[\Psi]=\frac{\braket{\Psi}{\hat{H}_{C}}{\Psi}}{\braket{\Psi}{\Psi}}=\sum_{ij}h_{ij}\sum_{\sigma}\gamma^{\sigma}_{ij}[\Psi]+\frac{1}{2}\sum_{ijkl}V_{ijkl}\sum_{\sigma\tau}\Gamma^{\sigma\tau}_{ijkl}[\Psi]$

(8)

in terms of the normalized 1- and 2-electron reduced density matrices (RDMs),

	$\displaystyle\gamma^{\sigma}_{ij}[\Psi]$	$\displaystyle=\braket{\Psi}{a^{\dagger}_{i\sigma}a_{j\sigma}}{\Psi}$		(9)
	$\displaystyle\Gamma^{\sigma\tau}_{ijkl}[\Psi]$	$\displaystyle=\braket{\Psi}{a^{\dagger}_{i\sigma}a^{\dagger}_{j\tau}a_{l\tau}a_{k\sigma}}{\Psi}.$		(10)

Because it is extremely difficult to determine, much less optimize, $E[\Psi]$ for arbitrary wavefunction forms, most practical approaches to the quantum many-body problem simplify the form of either the Hamiltonian or the wavefunction. We do not want to surrender our ability to accurately model real physical systems, so we will only approximate the wavefunction in this paper.

The particular model wavefunctions we consider are built from fully-paired states: they are geminal products. hurley ; mcweeny3 ; parr ; parks1 ; surjan4 ; surjan7 ; silver4 ; mehler ; silver5 ; kutz2 ; miller1 ; coleman1 ; coleman2 . It is convenient to work with objects that create/destroy pairs:

$\displaystyle S^{+}_{i}=a^{\dagger}_{i\uparrow}a^{\dagger}_{i\downarrow},\quad S^{-}_{i}=a_{i\downarrow}a_{i\uparrow},\quad S^{z}_{i}=\frac{1}{2}\left(a^{\dagger}_{i\uparrow}a_{i\uparrow}+a^{\dagger}_{i\downarrow}a_{i\downarrow}-1\right)$

(11)

with the structure,

$\displaystyle[S^{z}_{i},S^{\pm}_{j}]$

$\displaystyle=\pm\delta_{ij}S^{\pm}_{i},\quad[S^{+}_{i},S^{-}_{j}]=2\delta_{ij}S^{z}_{i}.$

(12)

The pairing scheme may be engineered to be more general bajdich2 ; limacher ; scuseria1 ) though the algebraic structure of the states does not change. The operators $S^{\pm}_{i}$ create/destroy a pair in state $i$, while $S^{z}_{i}$ counts the number of pairs in state $i$. It is convenient to work with the number operator

$\displaystyle\hat{n}_{i}=2S^{z}_{i}+1.$

(13)

The most general geminal mean-field wavefunction is the antisymmetric product of interacting geminals (APIG)silver2 ; silver1 for $N_{P}$ pairs in $N_{orb}$ spatial orbitals,

$\displaystyle\ket{\text{APIG}}=\prod^{N_{P}}_{\alpha=1}\sum^{N_{orb}}_{i=1}g^{i}_{\alpha}S^{+}_{i}\ket{\theta}.$

(14)

The state $\ket{\theta}$ represents the vacuum; it ordinarily represents the physical vacuum but it is only required to be a vacuum with respect to the removal of electron pairs. For APIG, the amplitudes $g^{i}_{\alpha}$ have no additional structure. Equation (14) can be expanded as a linear combination of $\binom{N_{orb}}{N_{P}}$ Slater determinants

$\displaystyle\ket{\text{APIG}}=\sum_{\begin{subarray}{c}m_{i}=\{0,1\}\\ \sum^{N_{orb}}_{i=1}m_{i}=N_{P}\end{subarray}}\left|\mathbf{C}(\mathbf{m})\right|^{+}\left(S^{+}_{1}\right)^{m_{1}}\left(S^{+}_{2}\right)^{m_{2}}\dots\left(S^{+}_{N_{orb}}\right)^{m_{N_{orb}}}\ket{\theta}.$

(15)

Each expansion coefficient is the permanent of an $N_{P}\times N_{P}$ matrix consisting of the geminal amplitudes of the occupied (i.e., $m_{i}=1$) orbital pairs,

	$\displaystyle\left|\mathbf{C}(\mathbf{m})\right|^{+}=\left|\begin{array}[]{cccc}g^{i_{1}}_{1}&g^{i_{2}}_{1}&\dots&g^{i_{N_{P}}}_{1}\\ g^{i_{1}}_{2}&g^{i_{2}}_{2}&\dots&g^{i_{N_{P}}}_{2}\\ \vdots&&\ddots&\vdots\\ g^{i_{1}}_{N_{P}}&g^{i_{2}}_{N_{P}}&\dots&g^{i_{N_{P}}}_{N_{P}}\end{array}\right|^{+}$		(20)
	$\displaystyle 1=m_{i_{1}}=m_{i_{2}}=\dots=m_{i_{N_{P}}}.$

Evaluating the permanent of a matrix is #P-hard valiant , so the APIG wavefunction is computationally intractable. For certain special forms of $g^{i}_{\alpha}$, however, the permanents may be efficiently evaluated johnson ; limacher ; borchardt ; carlitz .

It has long been recognized that for antisymmetric products of geminals (Eq. (14)) kutz1 and doubly-occupied configuration interaction wavefunctions (Eq. (15)) weinhold ; weinhold:1967b ; cook have especially simple formulas for the reduced density matrices. Specifically, the energy expression becomes

$\displaystyle E=2\sum_{i}h_{ii}\gamma_{i}+\sum_{i\neq j}(2V_{ijij}-V_{ijji})D_{ij}+\sum_{ij}V_{iijj}P_{ij}$

(21)

in terms of the normalized RDM elements

	$\displaystyle\gamma_{i}$	$\displaystyle=\frac{1}{2}\frac{\braket{\Psi}{\hat{n}_{i}}{\Psi}}{\braket{\Psi}{\Psi}}$		(22a)
	$\displaystyle D_{ij}$	$\displaystyle=\frac{1}{4}\frac{\braket{\Psi}{\hat{n}_{i}\hat{n}_{j}}{\Psi}}{\braket{\Psi}{\Psi}}$		(22b)
	$\displaystyle P_{ij}$	$\displaystyle=\frac{\braket{\Psi}{S^{+}_{i}S^{-}_{j}}{\Psi}}{\braket{\Psi}{\Psi}}.$		(22c)

The 1-RDM $\gamma_{i}$ is diagonal, while the 2-RDM has two non-zero pieces: the diagonal correlation function $D_{ij}$ and the pair correlation function $P_{ij}$. Notice that the diagonal terms refer to the same matrix element $P_{ii}=D_{ii}$ so to avoid double counting this term is assigned to $P_{ii}$. Elements of the 1-RDM represent the probability of occupying the $i$th orbital. The elements $D_{ij}$ represent the probability of occupying the $i$th and $j$th orbitals simultaneously, while elements of $P_{ij}$ represent the probability of transferring a pair from level $i$ to level $j$.

Since we are interested in systems that are well-described by pairing wavefunctions (14), we need Hamiltonians whose eigenstates are of this form. One such Hamiltonian is the reduced Bardeen-Cooper-Schrieffer (BCS) Hamiltonian bardeen:1957a ; BCS ; gaudin2 ; rich1 ; rich2 ; rich3 ,

$\displaystyle\hat{H}_{BCS}=\frac{1}{2}\sum^{N_{orb}}_{i=1}\varepsilon_{i}\hat{n}_{i}-\frac{g}{2}\sum_{ij}S^{+}_{i}S^{-}_{j}.$

(23)

The reduced BCS Hamiltonian is believed to provide qualitative insights into strongly correlated systems like superconducting nanograins dukelsky3 ; dukelsky4 ; kruchinin ; vonDelft . Richardsonrich1 ; rich2 ; rich3 showed that the eigenvectors of the reduced BCS Hamiltonian are products of geminals:

$\displaystyle\ket{\{u\}}=\prod^{N_{P}}_{a=1}S^{+}(u_{a})\ket{\theta}$

(24)

where the operators $S^{+}(u_{a})$ create a pair delocalized over the entire single particle space, i.e.

$\displaystyle S^{+}(u_{a})=\sum^{N_{orb}}_{i=1}\frac{S^{+}_{i}}{u_{a}-\varepsilon_{i}}.$

(25)

The complex numbers $\{u_{a}\}^{N_{P}}_{a=1}$ are called by various names, including rapidities, quasimomenta, spectral parameters, and pairing energies. The parameters in the reduced BCS Hamiltonian, $\{\varepsilon_{i}\}^{N_{orb}}_{i=1}$ and $g$, represent single-particle energies and the pairing strength respectively (cf. Eq. (23)). The two-electron states generated by Eq. (25) are called Cauchy geminals because $1/(u_{a}-\varepsilon_{i})$ are the elements of a Cauchy matrix. The form of the solution, Eq. (24)-(25), is also referred to as the rational solution for the XXX (isotropic) Richardson-Gaudin model.

The wavefunction form (24) is a particular example of the ABA: acting with the Hamiltonian (23) on the state (24) and collecting terms yields

$\displaystyle\hat{H}_{BCS}\ket{\{u\}}=\sum_{a}u_{a}\ket{\{u\}}-\frac{g}{2}\sum_{i}S^{+}_{i}\sum_{a}\Lambda(u_{a})\prod_{b\neq a}S^{+}(v_{b})\ket{\theta}$

(26)

one term proportional to (24) and $N_{P}$ linearly independent terms which vanish provided that

$\displaystyle\Lambda(u_{a})$

$\displaystyle=\frac{2}{g}+\sum^{N_{orb}}_{i=1}\frac{1}{u_{a}-\varepsilon_{i}}+\sum_{b\neq a}\frac{2}{u_{b}-u_{a}}=0.$

(27)

The set of equations (27) are Richardson’s equations and must be solved numerically. Many algorithms are available.rombouts ; faribault:2011 ; faribault_DG ; stijn_PRC ; guan:2012 ; pogosov:2012 ; claeys:2015

Computationally practical expressions for the RDM elements have been derived previously, so we will only present the main ideas. The engine driving all of these results is Slavnov’s theorem slavnov ; zhou ; claeys:2017b ; belliard:2019 , which gives a formula for the inner product of two RG states, Eqs. (24)-(25), where one of the sets of rapidities (here $\{u\}$) is a solution of Richardson’s equations (27), with the other set arbitrary. Specifically,

$\displaystyle\braket{\{u\}}{\{v\}}$

$\displaystyle=\frac{\prod_{a,b}\left(u_{a}-v_{b}\right)}{\prod_{a<b}\left(v_{a}-v_{b}\right)\left(u_{b}-u_{a}\right)}\det J\left(\{u\},\{v\}\right)$

(28)

where the matrix $\mathbf{J}$ has elements

$\displaystyle J_{ab}=\frac{1}{u_{a}-v_{b}}\left(\sum_{i}\frac{1}{\left(u_{a}-\varepsilon_{i}\right)\left(v_{b}-\varepsilon_{i}\right)}-2\sum_{c\neq a}\frac{1}{\left(u_{a}-u_{c}\right)\left(v_{b}-u_{c}\right)}\right).$

(29)

Choosing $\{v\}=\{u\}$ gives the norm of an RG state,

$\displaystyle\braket{\{u\}}{\{u\}}$

$\displaystyle=\det G$

(30)

where the elements of the matrix G are

$\displaystyle G_{ab}$

$\displaystyle=\begin{cases}\sum^{N}_{i=1}\frac{1}{(u_{a}-\varepsilon_{i})^{2}}-2\sum^{N_{P}}_{c\neq a}\frac{1}{(u_{a}-u_{c})^{2}}&a=b\\ \frac{2}{(u_{a}-u_{b})^{2}}&a\neq b\end{cases}.$

(31)

This expression for the norm was originally derived by Richardson by an alternative method rich3 . The Gaudin matrix, $G$, is the Jacobian of Richardson’s equations. I.e.,

$\displaystyle G_{ab}$

$\displaystyle=\frac{\partial\Lambda(u_{a})}{\partial u_{b}}$

(32)

where $\Lambda(u_{a})$ is one of Richardson’s equations, cf. Eq. (27). The expressions for the RDM elements are corollaries to Slavnov’s theorem obtained from the form factor approach. For the 1-RDM, we can use the commutation relations (12) to move $\hat{n}_{i}$ to the right until it acts on the vacuum, obtaining

$\displaystyle\frac{1}{2}\braket{\{u\}}{\hat{n}_{i}}{\{v\}}$

$\displaystyle=\sum_{a}\frac{\braket{\{u\}}{S^{+}_{i}}{\{v\}_{a}}}{v_{a}-\varepsilon_{i}}.$

(33)

Here $\ket{\{v\}_{a}}$ denotes the $(N_{P}-1)$-pair state with the rapidity $v_{a}$ removed, and the scalar product $\braket{\{u\}}{S^{+}_{i}}{\{v\}_{a}}$ is called a form factor. Similarly, the 2-RDM elements are obtained from

$\displaystyle\frac{1}{4}\braket{\{u\}}{\hat{n}_{i}\hat{n}_{j}}{\{v\}}$

$\displaystyle=\sum_{a\neq b}\frac{\braket{\{u\}}{S^{+}_{i}S^{+}_{j}}{\{v\}_{a,b}}}{(v_{a}-\varepsilon_{i})(v_{b}-\varepsilon_{j})}$

(34)

and

$\displaystyle\braket{\{u\}}{S^{+}_{i}S^{-}_{j}}{\{v\}}$

$\displaystyle=\sum_{a}\frac{\braket{\{u\}}{S^{+}_{i}}{\{v\}_{a}}}{v_{a}-\varepsilon_{j}}-\sum_{a\neq b}\frac{\braket{\{u\}}{S^{+}_{i}S^{+}_{j}}{\{v\}_{a,b}}}{(v_{a}-\varepsilon_{j})(v_{b}-\varepsilon_{j})}$

(35)

where $\{v\}_{a,b}$ is the set $\{v\}$ without $v_{a}$ and $v_{b}$. The expression (35) appears asymmetric in $i$ and $j$ since all the terms arise from moving $S^{-}_{j}$ to the right, but from many numerical testsfaribault1 ; faribault2 ; johnson:2020 ; moisset:2022 it is clear that it is in fact symmetric as it should be.

The geminal-creation operators, Eq. (25), have simple poles whenever the rapidity coincides with a single-particle energy $\{\varepsilon\}$. The residues of the poles are the pair creators

$\displaystyle S^{+}_{i}$

$\displaystyle=\lim_{v\rightarrow\varepsilon_{i}}(v-\varepsilon_{i})S^{+}(v).$

(36)

The form factors that appear in Eqs. (33)-(35) can therefore be evaluated as the corresponding residues of the scalar product, $\braket{\{u\}}{\{v\}}$

$\displaystyle\braket{\{u\}}{S^{+}_{i}}{\{v\}_{a}}=\lim_{v_{a}\rightarrow\varepsilon_{i}}(v_{a}-\varepsilon_{i})\braket{\{u\}}{\{v\}}$

(37)

$\displaystyle\braket{\{u\}}{S^{+}_{i}S^{+}_{j}}{\{v\}_{a,b}}=\lim_{v_{a}\rightarrow\varepsilon_{i}}\lim_{v_{b}\rightarrow\varepsilon_{j}}(v_{a}-\varepsilon_{i})(v_{b}-\varepsilon_{j})\braket{\{u\}}{\{v\}}.$

(38)

Expressions for $\gamma_{i}$, $D_{ij}$, and $P_{ij}$ are directly obtained by making the substitution $\{v\}\rightarrow\{u\}$ and normalizing. The form factors become ratios of determinants differing by one (or two) columns

	$\displaystyle\frac{\braket{\{u\}}{S^{+}_{i}}{\{u\}_{a}}}{\braket{\{u\}}{\{u\}}}$	$\displaystyle=\frac{\det G^{i}_{a}}{\det G}$		(39)
	$\displaystyle\frac{\braket{\{u\}}{S^{+}_{i}S^{+}_{j}}{\{u\}_{a,b}}}{\braket{\{u\}}{\{u\}}}$	$\displaystyle=\frac{\det G^{ij}_{ab}}{\det G}.$		(40)

The matrix $G^{i}_{a}$ is the matrix in (31) with the $a$th column replaced by the $i$th vector

$\displaystyle\textbf{r}_{i}=\begin{pmatrix}\frac{1}{(u_{1}-\varepsilon_{i})^{2}}\\ \frac{1}{(u_{2}-\varepsilon_{i})^{2}}\\ \vdots\\ \frac{1}{(u_{N_{P}}-\varepsilon_{i})^{2}}\\ \end{pmatrix}$

(41)

while $G^{ij}_{ab}$ is (31) with the $a$th column replaced by the $i$th version of (41) and the $b$th column replaced by the $j$th version of (41). It is not difficult to verify that the double column replacement can be simplified to a $2\times 2$ determinant of single column replacements,

$\displaystyle\frac{\det G^{ij}_{ab}}{\det G}=\frac{\det G^{i}_{a}}{\det G}\frac{\det G^{j}_{b}}{\det G}-\frac{\det G^{i}_{b}}{\det G}\frac{\det G^{j}_{a}}{\det G}.$

(42)

and rather than compute the $N_{P}\times N_{orb}$ determinants $G^{i}_{a}$, we can solve the $N_{orb}$ sets of linear equations

$\displaystyle G\textbf{x}=\textbf{r}_{i}.$

(43)

From Cramer’s rule, the solutions of these linear equations directly give the ratios of determinants

$\displaystyle\textbf{x}_{a}=\frac{\det G^{i}_{a}}{\det G}.$

(44)

The RDMs for on-shell RG states are thus constructed from solutions of linear equations. No determinants are computed numerically.