Geminal wave functions
APG, APSG and APIG
Using the general geminal creation operator (cfr. \(\eqref{eq:general_geminal_creator}\)), we can start proposing geminal wave function models for systems with an even number of electrons. One of these is the antisymmetrized product of geminals, or APG, a term first coined by Kutzelnigg (Kutzelnigg 1964). This wave function model is written as \[\begin{equation} \label{eq:APG} \require{physics} \ket{\text{APG}} % = \prod_{\gamma}^{N_P} % \hat{\Gamma}_\gamma^+ % \ket{\text{vac}} % = \prod_{\gamma}^{N_P} \qty( % \sum_{PQ}^M % G^{PQ}_\gamma % \hat{a}_P^\dagger \hat{a}_Q^\dagger % ) \ket{\text{vac}} % \thinspace , \end{equation}\] where \(N_P\) is the number of electron pairs. We can see that this geminal wave function model consists of \(N_P\) possibly different geminals.
Naturally, we could then come up with the following wave function model for an \(N\)-electron state, built up from these electron pair functions: \[\begin{equation} \label{eq:APSG} \Psi(\vb{r}_1, \cdots, \vb{r}_N) % = \mathscr{A} \qty[ % \omega_1(\vb{r}_1, \vb{r}_2) % \omega_2(\vb{r}_3, \vb{r}_4) % \cdots \omega_{N/2}(\vb{r}_{N-1}, \vb{r}_N) % ] % \thinspace , \end{equation}\] in which we have introduced \(N/2\) pair functions \(\set{\omega_R(\vb{r})}\) instead of the usual \(N(N-1)/2\). The operator \(\mathscr{A}\) is the antisymmetrizer that takes care of the resulting antisymmetry of the total wave function model. If for \(R \neq S\), we require the strong orthogonality condition \(\eqref{eq:strong_orthogonality}\) to be fulfilled, the resulting wave function model is called the antisymmetric product of strongly orthogonal geminals (APSG) P. R. Surján (2016). Imposing the strong orthogonality constraint, the geminal theory is fimplied. In practice, however, the strong orthogonality condition is too severe, and the APSG model does not account for interorbital correlations. (Miller et al. 1977)
We will refrain from calling these two-electron functions geminals, and reserve the term geminals for the algebraic description using two elementary creation operators.
Using the same set of spatial orbitals for every electron pair, the APIG (antisymmetric product of interacting geminals) (Silver 1969) wave function is written as \[\begin{align} \ket{\text{APIG}} &= \prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma \ket{\text{vac}} \label{eq:APIG_generalized_pair} \\ &= \prod_\gamma^{N_P} \qty( \sum_p^K G^p_\gamma \hat{P}^+_p ) \ket{\text{vac}} \label{eq:APIG} \thinspace , \end{align}\] in which the geminal coefficients \(G^p_\gamma\) are collected in the APIG coefficient matrix \(\vb{G}\): \[\begin{equation} \vb{G} = \qty( \begin{array}{ccccc|ccc} G^1_1 & \cdots & G^i_1 & \cdots & G^{N_P}_1 & G^{N_P+1}_1 & \cdots & G^K_1 \\ \vdots & \ddots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ G^1_\gamma & \cdots & G^i_\gamma & \cdots & G^{N_P}_\gamma & G^{N_P+1}_\gamma & \cdots & G^K_\gamma \\ \vdots & \ddots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ G^1_{N_P} & \cdots & G^i_{N_P} & \cdots & G^{N_P}_{N_P} & G^{N_P+1}_{N_P} & \cdots & G^K_{N_P} \end{array} ) \thinspace . \end{equation}\] This is an \(N_P \times K\)-matrix, where the left block is of dimensions \(N_P \times N_P\) and the right block is of dimensions \(N_P \times (K-N_P)\). Each row thus signifies a geminal \(\gamma\) (and which weight was given to each spatial orbital \(p\)), and every column thus collects the weights of the spatial orbital \(p\) in every geminal.
Writing out equation \(\eqref{eq:APIG}\) explicitly, we get \[\begin{equation} \ket{\text{APIG}} = \sum_{pq \cdots r}^K \thinspace G^p_1 G^q_2 \cdots G^r_{N_P} \thinspace \hat{P}^+_p \hat{P}^+_q \cdots \hat{P}^+_r \ket{\text{vac}} \thinspace , \end{equation}\] in which the summation sign encompasses \(N_P\) sum labels, and there are \(N_P\) factors in the coefficients and electron pair creators. By the Pauli principle, all the operators \(\hat{P}^+_p \hat{P}^+_q \cdots \hat{P}^+_r\) have to be different, which leads to \[\begin{equation} \ket{\text{APIG}} = \sum_{pq \cdots r}^{\in \binom{K}{N_P} N_P!} \thinspace G^p_1 G^q_2 \cdots G^r_{N_P} \thinspace \hat{P}^+_p \hat{P}^+_q \cdots \hat{P}^+_r \ket{\text{vac}} \thinspace , \end{equation}\] since we choose \(N_P\) distinguishable operators out of \(K\) possible ones without repetition. The summation sign means that the summation runs over all possible non-vanishing operator strings, where the string encompasses operators labelled with indices \(pq \cdots r\). Furthermore, pair creation operators commute, which leads to the constatation that many of operator strings are equal. Since we will permute all indices \(pq \cdots r\) (and the operators ), the `coefficient’ that belongs to a certain unique operator string becomes a . The number of terms that remain is \[\begin{equation} \binom{K}{N_P} \thinspace , \end{equation}\] since there are \(N_P!\) unique permutations, and we can write \[\begin{equation} \ket{\text{APIG}} = \sum_{\vb{m}}^{\binom{K}{N_P}} |\vb{G}(\vb{m})|_+ \qty( \prod_p^K \qty( \hat{P}^+_p )^{m_p} ) \ket{\text{vac}} \thinspace , \end{equation}\] where we have introduced a pair-occupancy vector \(\vb{m} = \qty(m_1, m_2, \ldots, m_K)\), which collects the pair occupancies of every spatial orbital \(1 \cdots K\) with the limitation that \[\begin{equation} \sum_p^K m_p = N_P \thinspace , \end{equation}\] and \(|\vb{G}(\vb{m})|_+\) is the permanent of the \(N_P \times N_P\)-matrix containing only the columns (i.e. the spatial orbitals) of \(\vb{G}\) for which the spatial occupation number is equal to one (\(m_p = 1\)). In the end, this means that APIG is an approximation to DOCI, in the sense that the DOCI coefficients must be able to be written as permanents Limacher et al. (2013). This means that any APIG flavor has in its expansion every seniority zero ONV, but the simplifications in the flavors lead to more and more severe restrictions on the coefficients. The coefficients don’t each have full flexibility to vary (or we would have the linear variation method inside the DOCI space), since they depend on one another in some non-linear fashion.
In other words, we the overlap between a doubly-occupied ONV \(\ket{\vb{m}}\) and the APIG wave function, i.e. the corresponding expansion coefficient in the seniority-zero wave function, is: \[\begin{equation} \braket{\vb{m}}{\text{APIG}} = |\vb{G}(\vb{m})|_+ \thinspace . \end{equation}\] Written out explicitly for a few cases, we have: \[\begin{equation} \label{eq:APIG_overlap_paired} \braket{\Phi_0}{\text{APIG}} = \qty| \begin{array}{ccc} G^1_1 & \cdots & G^{N_P}_1 \\ \vdots & \ddots & \vdots \\ G^1_{N_P} & \cdots & G^{N_P}_{N_P} \end{array} |_+ \thinspace , \end{equation}\] \[\begin{equation} \label{eq:APIG_overlap_single_pair_excitation} \braket{\Phi_i^a}{\text{APIG}} = \qty| \begin{array}{ccccc} G^1_1 & \cdots & G_1^a & \cdots & G^{N_P}_1 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ G^1_{N_P} & \cdots & G^a_{N_P} & \cdots & G^{N_P}_{N_P} \end{array} |_+ \thinspace , \end{equation}\] and \[\begin{equation} \label{eq:APIG_overlap_double_pair_excitation} \braket{\Phi_{ij}^{ab}}{\text{APIG}} = \qty| \begin{array}{ccccccc} G^1_1 & \cdots & G_1^a & \cdots & G_1^b & \cdots & G^{N_P}_1 \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ G^1_{N_P} & \cdots & G_{N_P}^a & \cdots & G_{N_P}^b & \cdots & G^{N_P}_{N_P} \\ \end{array} |_+ \thinspace , \end{equation}\] respectively. In equation \(\eqref{eq:APIG_overlap_single_pair_excitation}\), the \(i\)-th column has been replaced with the \(a\)-th column of the coefficient matrix, and in equation \(\eqref{eq:APIG_overlap_double_pair_excitation}\) (with (\(i<j\))), the \(i\)-th column has been replaced with the \(a\)-th column of the coefficient matrix, and similarly for \(j\) and \(b\).
Some more commutators
The commutators \(\eqref{eq:geminal_commutators_first}\) through \(\eqref{eq:comm_N_geminal}\) can repeated, but now for a string of \(N_P\) geminal creators: \[\begin{align} & \comm{\hat{P}^\circ_p}{\prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma} = \qty( \sum_\gamma^{N_P} G^p_\gamma \qty(\prod_{\delta \neq \gamma}^{N_P} \hat{\Gamma}^+_\delta)) \hat{P}^+_p \\ & \comm{\hat{P}^+_p}{\prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma} = 0 \\ & \comm{\hat{P}^-_p}{\prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma} = -2 \qty( \sum_\gamma^{N_P} \sum_{\delta = \gamma + 1}^{N_P} G^p_\gamma G^p_\delta \hat{\Gamma}^+_1 \cdots \hat{\Gamma}^+_{\gamma - 1} \qty(\prod_{\subalign{\lambda &= \gamma + 1 \\ &\neq \delta}}^{N_P} \hat{\Gamma}^+_\lambda)) \hat{P}^+_p \notag \\ & \phantom{\comm{\hat{P}^-_p}{\prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma} = } - 2 \qty(\sum_\gamma^{N_P} G^p_\gamma \qty(\prod_{\delta \neq \gamma}^{N_P} \hat{\Gamma}^+_\gamma)) \hat{P}^\circ_p \\ & \comm{\hat{N}_p}{\prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma} = 2 \qty( \sum_\gamma^{N_P} G^p_\gamma \qty(\prod_{\delta \neq \gamma}^{N_P} \hat{\Gamma}^+_\delta)) \hat{P}^+_p \\ & \comm{\hat{N}}{\prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma} = 2 N_P \prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma \label{eq:comm_N_geminal_string} \thinspace . \end{align}\] From \(\eqref{eq:comm_N_geminal_string}\) follows that any APIG wave function, i.e. of the form \(\eqref{eq:APIG}\) contains \(2 N_P\) (2 for every pair) electrons: \[\begin{equation} \hat{N} \ket{\text{APIG}} = 2 N_P \ket{\text{APIG}} \thinspace . \end{equation}\] Since the commutator of the seniority operator with the pair operators vanishes, we can furthermore show that any APIG wave function has seniority zero: \[\begin{equation} \label{eq:seniority_zero_APIG} \hat{\Omega} \ket{\text{APIG}} = 0 \thinspace , \end{equation}\] as is expected.
AGP
A conceptually simpler wave function model than APG is then given by AGP, or the antisymmetric geminal power: \[\begin{equation} \label{eq:AGP} \ket{\text{AGP}} % = \qty( % \hat{\Gamma}^+ % )^{N_P} \ket{\text{vac}} % = \qty( % \sum_{PQ}^M % G^{PQ} % \hat{a}_P^\dagger \hat{a}_Q^\dagger % )^{N_P} \ket{\text{vac}} % \thinspace , \end{equation}\] in which every geminal is constrained to be the same. Introducing the pair wave function \(\omega(\vb{r}_1, \vb{r}_2)\) that is related to the geminal creator in \(\eqref{eq:AGP}\), we can write the AGP model in first quantization as: \[\begin{equation} \Psi_{\text{AGP}}(\vb{r}_1, \cdots, \vb{r}_N) % = \mathscr{A} \qty[ % \omega(\vb{r}_1, \vb{r}_2) % \omega(\vb{r}_3, \vb{r}_4) % \cdots \omega(\vb{r}_{N-1}, \vb{r}_N) % ] % \thinspace , \end{equation}\] where we see that the same pair function \(\omega(\vb{r}_1, \vb{r}_2)\) is used \(N_P\) times, and we let \(\mathscr{A}\) be the antisymmetrizer that takes care of the overall antisymmetry of the total wave function. (P. R. Surján 2016)
If we let every generalized pair \(\gamma\) be equal in the APIG model, i.e. requiring \[\begin{equation} \hat{\Gamma}^+ = \sum_p^K G_p \hat{P}^+_p \thinspace , \end{equation}\] the APIG wave function reduces to \[\begin{equation} \ket{\text{AGP}} = \qty( \hat{\Gamma}^+ )^{N_P} \ket{\text{vac}} \thinspace , \end{equation}\] being an AGP wave function: an antisymmetric geminal power. (Paul Andrew Johnson 2014) The coefficient matrix then becomes \[\begin{equation} \vb{G} = \qty( \begin{array}{ccc|ccc} G_1 & \cdots & G_{N_P} & G_{N_P+1} & \cdots & G_K \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ G_1 & \cdots & G_{N_P} & G_{N_P+1} & \cdots & G_K \\ \end{array} ) \thinspace , \end{equation}\] which is an \(N_P \times K\)-matrix with equal rows, with \(K\) unique coefficients. The simplifications to equations \(\eqref{eq:APIG_overlap_paired}\), \(\eqref{eq:APIG_overlap_single_pair_excitation}\), and \(\eqref{eq:APIG_overlap_double_pair_excitation}\) are \[\begin{equation} \braket{\Phi_0}{\text{AGP}} = N_P! \thinspace \prod_i^{N_P} G_i \thinspace , \end{equation}\] \[\begin{equation} \braket{\Phi_i^a}{\text{AGP}} = N_P! \thinspace G_a \prod_{j \neq i}^{N_P} G_j \thinspace , \end{equation}\] and \[\begin{equation} \braket{\Phi_{ij}^{ab}}{\text{AGP}} = N_P! \thinspace G_a G_b \prod_{k \neq i, j}^{N_P} G_k \thinspace . \end{equation}\]
Simplifications
The evaluation of a permanent is #P-hard (Valiant 1979), which means that for large dimension the computation becomes intractable, as can be seen from Figure \(\ref{fig:permanent}\). This naturally leads to a search of approximations to APIG, in which the coefficients are written in special forms, in order to overcome the particular computational difficulties.
The AP1roG wave function
The AP1roG (antisymmetric product of 1-reference-orbital geminals) wave function (Limacher et al. 2013) will be introduced as \[\begin{align} \ket{\text{AP1roG}} &= \prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma \ket{\text{vac}} \\ &= \prod_\gamma^{N_P} \qty( \hat{P}^+_\gamma + \sum_{a = N_P + 1}^K G_\gamma^a \hat{P}^+_a ) \ket{\text{vac}} \label{eq:AP1roG} \thinspace , \end{align}\] where \(a\) labels virtual orbitals, so that the coefficient matrix becomes \[\begin{equation} \vb{G} = \qty( \begin{array}{cccc|ccc} 1 & 0 & \cdots & 0 & G^{N_P+1}_1 & \cdots & G_1^K \\ 0 & 1 & \cdots & 0 & G^{N_P+1}_2 & \cdots & G_2^K \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 & G_{N_P}^{N_P+1} & \cdots & G^{K}_{N_P} \end{array} ) \thinspace , \end{equation}\] in which the left \(N_P \times N_P\)-block is an identity matrix, and the right \(N_P \times (K - N_P)\)-block contains the AP1roG coefficients. The indices of the coefficients then are \[\begin{align} &a=N_P+1, \ldots, K \label{eq:coefficients_ap1rog_matrix_1} \\ &i=1, \ldots, N_P \label{eq:coefficients_ap1rog_matrix_2} \thinspace . \end{align}\]
That AP1roG is a special case of APIG can be seen by recognizing that the coefficients can be written as \[\begin{equation} G^p_{\gamma, \text{AP1roG}} = \delta_{p \gamma} + \qty( \sum_k^{i-1} \delta_{k, N_P} ) G^p_{\gamma, \text{APIG}} \thinspace . \end{equation}\] Furthermore, the AP1roG wave function can be written as a pair coupled-cluster doubles wave function Stein, Henderson, and Scuseria (2014): \[\begin{equation} \ket{\text{AP1roG}} = \exp( \sum_\gamma^{N_P} \sum_{a=N_P+1}^K G^a_\gamma \hat{P}^+_a \hat{P}^-_\gamma ) \ket{\Phi_0} \thinspace , \end{equation}\] of which a detailed derivation can be found below.
For AP1roG, we know the derivatives: \[\begin{equation} \ket{ % \pdv{ \Psi(\vb{G}) }{ G_i^a } % } \equiv % \pdv{G_i^a} \ket{\Psi(\vb{G})} % = \hat{P}^+_a \hat{P}^-_i \ket{ \Psi(\vb{G}) } % \thinspace , \end{equation}\] which can be confirmed to be equal to \[\begin{equation} \ket{ % \pdv{ \Psi(\vb{G}) }{ G_i^a } % } % = \hat{P}^+_b \qty[ % \prod_{\gamma \neq i}^{N_P} \qty( % \hat{P}^+_\gamma % + \sum_{a = N_P + 1}^K G_\gamma^a \hat{P}^+_a % ) % ] \ket{\text{vac}} \end{equation}\] as in (De Baerdemacker 2017). Since these derivatives are so simple, we can quickly calculate overlaps with pair-excited reference determinants: \[\begin{align} & \braket{\Phi_0}{ % \pdv{ \Psi(\vb{G}) }{ G_i^a } % } = 0 \\ % & \braket{\Phi_j^b}{ % \pdv{ \Psi(\vb{G}) }{ G_i^a } % } = \delta_{ij} \delta_{ab} \\ % & \braket{\Phi_{jk}^{bc}}{ % \pdv{ \Psi(\vb{G}) }{ G_i^a } % } = (1 - \delta_{jk}) % (1 - \delta_{bc}) \qty( % \delta_{ij} \delta_{ab} G_k^c % + \delta_{ij} \delta_{ac} G_k^b % + \delta_{ik} \delta_{ab} G_j^c % + \delta_{ik} \delta_{ac} G_j^b % ) % \thinspace . \end{align}\]
The APr2G wave function
Let’s take the Richardson-Gaudin Hamiltonian (Paul Andrew Johnson 2014) \[\begin{equation} \hat{\mathcal{H}}_{\text{RG}} = \sum_p^K \varepsilon_p \hat{P}^\circ_p - g \sum_{pq}^K \hat{P}^+_p \hat{P}^-_q \thinspace, \end{equation}\] which is also called the reduced BCS Hamiltonian, and introduce the following geminal operators: \[\begin{align} & \hat{\Gamma}^+_\gamma = \sum_p^K \frac{\hat{P}^+_p}{\lambda_\gamma - \varepsilon_p} \\ & \hat{\Gamma}^-_\gamma = \sum_p^K \frac{\hat{P}^-_p}{\lambda_\gamma - \varepsilon_p} \\ & \hat{\Gamma}^\circ_\gamma = \frac{1}{g} - \sum_p^K \frac{\hat{P}^\circ_p}{\lambda_\gamma - \varepsilon_p} \thinspace , \end{align}\] which just boils down to a simplification of the APIG model: \[\begin{equation} G^p_\gamma = \frac{1}{\lambda_\gamma - \varepsilon_p} \thinspace . \end{equation}\] It is important to note that \[\begin{equation} \hat{\Gamma}^-_\gamma \neq \qty(\hat{\Gamma}^+_\gamma)^\dagger \thinspace . \end{equation}\]
Let us furthermore introduce the Bethe ansatz as an ansatz for the eigenvectors of the RG-Hamiltonian: \[\begin{equation} \ket{\set{\lambda}} = \prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma \ket{\text{vac}} \thinspace , \end{equation}\] and calculate \[\begin{equation} \hat{\mathcal{H}}_{\text{RG}} \ket{\text{vac}} = - \frac{1}{2} \sum_p^K \varepsilon_p \thinspace , \end{equation}\] and the commutator \[\begin{equation} \comm{\hat{\mathcal{H}}_{\text{RG}}}{\hat{\Gamma}^+_\gamma} = \lambda_\gamma \hat{\Gamma}^+_\gamma + \sum_p^K \hat{P}^+_p \qty( 1 - 2g \hat{\Gamma}^\circ_\gamma ) \thinspace . \end{equation}\] From this commutator follows the nested commutator \[\begin{equation} \comm{\comm{\hat{\mathcal{H}}_{\text{RG}}}{\hat{\Gamma}^+_\gamma}}{\hat{\Gamma}^+_\delta} = \frac{2g}{\lambda_\delta - \lambda_\gamma} \sum_p^K \hat{P}^+_p \qty( \hat{\Gamma}^+_\gamma - \hat{\Gamma}^+_\delta ) \thinspace , \end{equation}\] from which we can see that any higher-order nested commutator in this fashion vanishes, which substantially simplifies further calculations. By moving \(\hat{\mathcal{H}}_{\text{RG}}\) to the right, we will, by brute force (Paul Andrew Johnson 2014), calculate that \[\begin{equation} \hat{\mathcal{H}}_{\text{RG}} \ket{\set{\lambda}} = E \ket{\set{\lambda}} \thinspace , \end{equation}\] with \[\begin{equation} E = - \frac{1}{2} \sum_p^K \varepsilon_p + \sum_\gamma^{N_P} \lambda_\gamma \thinspace , \end{equation}\] if and only if Richardson’s equations are fulfilled: \[\begin{equation} \label{eq:richardson} \forall \lambda_\gamma: \qquad - \frac{1}{g} - \sum_p^K \frac{1}{\lambda_\gamma - \varepsilon_p} + 2 \sum_{\delta \neq \gamma}^{N_P} \frac{1}{\lambda_\gamma - \lambda_\delta} = 0 \thinspace . \end{equation}\] A Bethe wave function in which the parameters satisfy Richardson’s equations are called on-shell.
Inspired by the procedure described in the previous section, let us take the following geminal: \[\begin{equation} \hat{\Gamma}^+_\gamma = \sum_p^K \frac{G^p_\gamma \hat{P}^+_p}{\lambda_\gamma - \varepsilon_p} \thinspace , \end{equation}\] which is a realization of the APIG-model, in which the coefficients are \[\begin{equation} G^p_{\gamma, \text{APIG}} = \frac{G^p_\gamma}{\lambda_\gamma - \varepsilon_p} \thinspace , \end{equation}\] such that \[\begin{align} \ket{\text{APr2G}} &= \prod_\gamma^{N_P} \hat{\Gamma}^+_\gamma \ket{\text{vac}} \\ &= \prod_\gamma^{N_P} \qty( \sum_p^K \frac{G^p_\gamma \hat{P}^+_p}{\lambda_\gamma - \varepsilon_p} ) \ket{\text{vac}} \thinspace , \end{align}\] in which APr2G stands for antisymmetric product of rank-two geminals. For the overlap with reference determinants, we find \[\begin{equation} \braket{\Phi_0}{\text{APr2G}} = \qty| \begin{array}{ccc} \frac{G^1_1}{\lambda_1 - \varepsilon_1} & \cdots & \frac{G^{N_P}_1}{\lambda_1 - \varepsilon_{N_P}} \\ \vdots & \ddots & \vdots \\ \frac{G^1_{N_P}}{\lambda_{N_P} - \varepsilon_1} & \cdots & \frac{G^{N_P}_{N_P}}{\lambda_{N_P} - \varepsilon_{N_P}} \end{array} |_+ \thinspace , \end{equation}\] \[\begin{equation} \braket{\Phi_i^a}{\text{APr2G}} = \qty| \begin{array}{ccccc} \frac{G^1_1}{\lambda_1 - \varepsilon_1} & \cdots & \frac{G^a_1}{\lambda_1 - \varepsilon_a} & \cdots & \frac{G^{N_P}_1}{\lambda_1 - \varepsilon_{N_P}} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \frac{G^1_{N_P}}{\lambda_{N_P} - \varepsilon_1} & \cdots & \frac{G^a_{N_P}}{\lambda_{N_P} - \varepsilon_a} & \cdots & \frac{G^{N_P}_{N_P}}{\lambda_{N_P} - \varepsilon_{N_P}} \end{array} |_+ \thinspace , \end{equation}\] and \[\begin{equation} \braket{\Phi_{ij}^{ab}}{\text{APr2G}} = \qty| \begin{array}{ccccccc} \frac{G^1_1}{\lambda_1 - \varepsilon_1} & \cdots & \frac{G^a_1}{\lambda_1 - \varepsilon_a} & \cdots & \frac{G^b_1}{\lambda_1 - \varepsilon_b} & \cdots & \frac{G^{N_P}_1}{\lambda_1 - \varepsilon_{N_P}} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ \frac{G^1_{N_P}}{\lambda_{N_P} - \varepsilon_1} & \cdots & \frac{G^a_{N_P}}{\lambda_{N_P} - \varepsilon_a} & \cdots & \frac{G^b_{N_P}}{\lambda_{N_P} - \varepsilon_b} & \cdots & \frac{G^{N_P}_{N_P}}{\lambda_{N_P} - \varepsilon_{N_P}} \end{array} |_+ \thinspace , \end{equation}\] respectively. Now, this doesn’t seem like much of an improvement, but we can see that the matrices of which the permanents have to be computed are Cauchy matrices, for which Borchardt’s theorem (Borchardt 1857) holds, which simplifies the calculation of the permanents.
If we determine \(\set{\lambda_\gamma}\) and \(\set{\varepsilon_p}\) from Richardson’s equation \(\eqref{eq:richardson}\), we call the resulting geminals on-shell. If they do not satisfiy Richardson’s equation, then we call them off-shell. In a sense, solving Richardson’s equations is much like a Hartree-Fock procedure. What is meant by this is that in Hartree-Fock, we find a set of optimal spin orbitals for a single Slater determinant, which in a sense fixes \(\hat{P}^+_p\), whereas solving Richardson’s equations determines the parameters, \(\set{\lambda_\gamma}\) and \(\set{\varepsilon_p}\), such that in a way the analoguous operator \(\frac{\hat{P}^+_p}{\lambda_\gamma - \varepsilon_p}\) is fixed.