Exploring non-linear correlators on AGP

Let the $k$-body, or rank-$k$, Hilbert-space Jastrow operator Neuscamman (2013b) be written as

where the sums run over all orbitals, $\alpha$ is a symmetric tensor and is invariant under the exchange of its indices, and ${1}/{2^{k}}$ is introduced for convenience. Importantly, since $\mathbf{N}_{p}^{2}=2\mathbf{N}_{p}$ in the seniority-zero space, we impose that $\alpha$ is zero if any two indices are the same and eliminate these terms from the sum in Eq. (13).

From the $k$-body Jastrow operator we can define the $k$-body Jastrow-AGP (J_kAGP) wavefunction as

As we prove in Appendix. A, we do not need to include lower body Jastrow operators in this J_kAGP wavefunction, as they are contained inside $\boldsymbol{J}_{k}$. Thus, for example, we may write

where $\boldsymbol{J}_{2}^{\prime}$ is another two-body Jastrow operator.

Slater determinants built from the same orbitals used to define the Hilbert-space Jastrow operators are eigenfunctions of $\boldsymbol{J}_{k}$, that is

where $|\Phi_{\mu}\rangle$ is a determinant, $N_{p,\mu}$ is the occupation number of level $p$ in determinant $\mu$, and the indices $i_{1}\ldots i_{k}$ run over the occupied orbitals in $|\Phi_{\mu}\rangle$. From this, it follows that

For a complex $\alpha$, the magnitude of the coefficient of a given determinant is modified by the Hermitian part of $\boldsymbol{J}_{k}$, i.e. the real part of $\alpha$. The anti-Hermitian part of $\boldsymbol{J}_{k}$ (the imaginary part of $\alpha$) only adjusts the phase of each coefficient. It is thus apparent that for an $N$-body initial wavefunction $|\Psi\rangle$ in which all determinants $|\Phi_{\mu}\rangle$ have non-zero coefficients, $\exp(\boldsymbol{J}_{N})|\Psi\rangle$ can be made an exact eigenstate of a chosen Hamiltonian provided that $\boldsymbol{J}_{N}$ is non-Hermitian.

AGP is, in general, such an initial wavefunction, so $\exp(\boldsymbol{J}_{N})$ acting on it can be made exact for the seniority zero space. For any $k\leq N$, it follows from Eq. (16) and Eq. (17) that

In contrast, we can see how the linear wavefunction $|{\Psi}\rangle=\boldsymbol{J}_{k}|\text{AGP}\rangle$, which we refer to as J_kCI AGP Dutta et al. (2020), modifies AGP as follows

Thus, the way in which the linear and non-linear $\boldsymbol{J}_{k}$ correlators factorize the DOCI coefficients can be easily deduced from Eq. (18) and Eq. (19). Together with Eq. (10), we can readily see that

Indeed the two factorizations are different, and Eq. (20a) can be understood as the linear approximation to Eq. (20b) up to a constant shift of all amplitudes. An implication can be made, for example, about the sign of each coefficient; if $\alpha$ is real, the DOCI coefficient in Eq. (20b) can be negative only if some $\eta$’s are negative, whereas the same restriction is not implied in (20a).

Having shown analytically how $\boldsymbol{J}_{k}$ acts on AGP, we concern ourselves with optimizing the energy over JAGP in the subsequent sections.

Consider $k=1$ for which the J₁AGP wavefunction is just an AGP,

Evidently, the action of $\boldsymbol{J}_{1}$ just changes the geminal coefficient $\eta_{p}$ to $\eta_{p}\,\mathrm{e}^{\alpha_{p}}$. For $k=1$, J₁AGP can optimize an AGP, but it cannot correlate it; we need $k>1$ to go beyond an AGP mean-field. Unfortunately, except for $k=1$, the norm and density matrices over $\mathrm{e}^{\boldsymbol{J}_{k}}|\text{AGP}\rangle$ are not trivial to compute, and aside from quantum Monte Carlo methods, Foulkes et al. (2001); Casula et al. (2004); Neuscamman (2013a, 2016); Wei and Neuscamman (2018); Cohen et al. (2019) we are not aware of an efficient non-stochastic algorithm to compute the overlaps. This for examples makes the variational optimization of JAGP energy

out of reach for the existing non-stochastic methods.

However, a case can be made for the similarity transformed Hamiltonian using $\boldsymbol{J}_{2}$, that is

The two-body Jastrow correlator has a long history in electronic structure methods and has shown to be capable of recovering a significant portion of the correlation energy. Casula and Sorella (2003); Neuscamman (2013a); Genovese et al. (2019) Previous work in our group has shown that the similarity transformation of fermionic Hamiltonians using a two-body Jastrow operator is summable to all orders and the amplitudes can be solved by left projection via components of $\boldsymbol{J}_{2}$ acting on the HF reference Wahlen-Strothman et al. (2015); Wahlen-Strothman and Scuseria (2016). Following the same line of reasoning, we extend that formalism to the pairing $su(2)$ algebra Eq. (4) and show how to obtain the working equations for AGP in the next section.

The main idea relies on the fact that the similarity transformation of the Hamiltonian reduces the rank of $\boldsymbol{J}_{2}$ by one, thus the remaining one-body Jastrow operator can be subsequently absorbed into the reference. The same idea can be applied to the unitary transformation of the Hamiltonian with an anti-Hermitian $i\boldsymbol{J}_{2}$ operator

wherein $\alpha$ is assumed to be real. Based on the argument in the previous section, the Jastrow correlators with purely imaginary amplitudes only contribute to the phase of each Slater determinant in AGP and do not alter their magnitudes. This can also be seen from Eq. (18) by $\alpha\rightarrow i\alpha$. The phase contributions could be important for Hamiltonians that contain complex terms including those that break time-reversal symmetry. For other Hamiltonians, the mere change in phase is insufficient to make the wavefunction exact.

In this section, we extend the results of Ref. Wahlen-Strothman et al. (2015); Wahlen-Strothman and Scuseria (2016) to the AGP wavefunction and our $su(2)$ algebra. For pedagogical reasons, we keep this section self-contained. Suppose one had a generic Jastrow operator $\boldsymbol{J}$. In a manner similar to coupled cluster theory, we intend to solve for the ground state energy by similarity transforming the Hamiltonian,

Similarity transformation does not change the spectrum of $\boldsymbol{H}$, but $\bar{\boldsymbol{H}}$ is no longer Hermitian and $\bar{E}$ is not variational. Therefore we obtain the residual equations by left projecting the Schrödinger equation to get

where the expectation values are taken with respect to AGP. In traditional coupled cluster theory using particle-hole excitations, one expresses the similarity transformation by the Baker-Campbell-Hausdorff (BCH) expansion whose series naturally truncates at the 4th commutator for a two-body Hamiltonian. In contrast, the transformation using the Jastrow operator does not truncate, but it is summable to all orders for any Hamiltonian.

Our ultimate goal is to use a two- or higher-body Jastrow operators to similarity transform the Hamiltonian and provide correlation. However, it will prove useful to begin by working through the case in which we use a one-body operator, for which the algebra is much simpler. Therefore, consider the $\boldsymbol{J}_{1}$ operator

Using Eq. (4), it is easy to show that

Thus, the similarity-transformed Hamiltonian simply adjusts the coefficients of the $\mathbf{P}^{\dagger}$ and $\mathbf{P}$ operators.

Alternatively, we can use Eq. (21) to absorb $\exp(\boldsymbol{J}_{1})$ into AGP. Although we have seen how this can be done, we will extract the same result in an different fashion. Notice by Eq. (5) that

where we used the fact that $\exp(\boldsymbol{J}_{1})|\text{vac}\rangle=|\text{vac}\rangle$. From the algebra, Eq. (4), it follows

from which it is clear that $\mathrm{e}^{\boldsymbol{J}_{1}}|\text{AGP}\rangle$ leads to $\eta_{p}\rightarrow\eta_{p}\mathrm{e}^{\alpha_{p}}$ for all geminal coefficients in the ket, and $\langle\text{AGP}|\mathrm{e}^{-\boldsymbol{J}_{1}}$ results in $\eta_{p}\rightarrow\eta_{p}\mathrm{e}^{-\alpha_{p}}$ in the bra.

We now proceed to derive the equations for similarity transformation with $\boldsymbol{J}_{2}$. Define $\boldsymbol{J}_{1}(p)=\sum_{q}\alpha_{pq}\mathbf{N}_{p}/4$ and recall $\alpha_{pp}=0$, then as in Eq. (28) we obtain

where we used $[{\boldsymbol{J}_{1}(p)},{\mathbf{P}^{\dagger}_{p}}]=0$ to symmetrically distribute the exponentials in the right hand side of Eq. (31) on either side of each operator. In fact, using these equations, we can transform any many-body operator of the form $\mathbf{P}^{\dagger}_{p}...\mathbf{N}_{q}...\mathbf{P}_{r}...$ in a symmetric way. In general, it can be shown that

To obtain a closed-form expression for the energy and the residual equations, we need to further absorb each $\boldsymbol{J}_{1}(p)$ into AGP. Following the same steps as in Eq. (30), we can see that

This together with Eq. (32) show that we can evaluate the expectation value of any similarity transformed operator of the form $\mathbf{P}^{\dagger}_{p}...\mathbf{N}_{q}...\mathbf{P}_{r}...$ with $\boldsymbol{J}_{2}$ by scaling the geminal coefficients appropriately.

The same basic steps follow for even higher-body Jastrow operators; the initial similarity-transformation of the Hamiltonian leads to modified matrix elements along with exponentials of Jastrow operators whose ranks are reduced by one. However, only $\boldsymbol{J}_{2}$ permits easy evaluation of the remaining expectation values.

For unitary JAGP, as in Eq. (24), one follows the same steps as above but replaces $\alpha\rightarrow i\alpha$.

We now show how to apply the equations above to similarly transform the pairing Eq. (2). It follows form Eq. (32) that

For every $\mathbf{P}^{\dagger}_{p}\mathbf{P}_{q}$, define $|\widetilde{\mathrm{AGP}}_{pq}\rangle=\mathrm{e}^{\boldsymbol{J}_{1}(q)-\boldsymbol{J}_{1}(p)}|\text{AGP}\rangle$, and note that $|\widetilde{\mathrm{AGP}}_{pq}\rangle$ is an AGP whose geminal coefficients have been rescaled by $\tilde{\eta}_{r}=\eta_{r}\mathrm{e}^{(\alpha_{qr}-\alpha_{pr})/2}$ for all orbitals $r$ from Eq. (33). As such, we can obtain the energy,

Following similar steps as above, it is straightforward to get closed-form expressions for the residual equations, Eq. (26), as well. We simply use

which can then be normal ordered in terms of $\mathbf{P}^{\dagger}_{p}\ldots\mathbf{N}_{q}\ldots\mathbf{P}_{r}\ldots$.

As outlined in Sec. II, it is a property of AGP-based RDMs that $Z_{p}^{(1,1)}=\langle\text{AGP}|{\mathbf{N}_{p}}|\text{AGP}\rangle$ is sufficient to construct all higher-body density matrices, where the cost of building $Z^{(1,1)}$ is $\mathcal{O}(N_{\mathrm{grid}}M)$ using grid integration in the number-projected BCS representation of AGP. Of course the same principle applies to RDMs using $|\widetilde{\mathrm{AGP}}_{pq}\rangle$, except that $\widetilde{Z}^{(1,1)}_{r;pq}=\langle\widetilde{\mathrm{AGP}}_{pq}|\mathbf{N}_{r}|\widetilde{\mathrm{AGP}}_{pq}\rangle$ needs to be evaluated for every $pq$ in the Hamiltonian, as demonstrated for Eq. (35). The cost of building $\widetilde{Z}^{(1,1)}_{r;pq}$ for all $pqr$ is $\mathcal{O}(N_{\mathrm{grid}}M^{3})$ and is the most expensive step in evaluating the ST-J₂AGP energy. The cost of evaluating the residuals using the reconstruction formulae is $\mathcal{O}(M^{4})$. Without the reconstruction formulae, the cost would be $\mathcal{O}(M^{6})$.

In Fig. 1, we show the energy error and the fraction of correlation energy captured by ST-J₂AGP. All calculations were carried out using the energy optimized AGP for $\boldsymbol{H}$. The correlation energy is computed as the difference with that of HF. As shown in Fig. 1, ST-J₂AGP over-correlates near the recoupling regime, $G\approx G_{c}$ and smaller, including $G<0$. In the strongly attractive correlated limit, $G\gg G_{c}$, ST-J₂AGP is highly accurate and becomes exact for sufficiently large $G/G_{c}$. It is noteworthy, however, that AGP is an eigenstate of the pairing Hamiltonian for $G/G_{c}\rightarrow\infty$. For comparison purposes, we have plotted the results for variational J₂AGP where we carried out the optimization using an in-house DOCI code. We have also included results using a linear wavefunction $\boldsymbol{J}_{2}|\mathrm{AGP}\rangle$ in which the energy is determined variationally; this J₂CI-AGP has $\mathcal{O}(M^{6})$ scaling, or $\mathcal{O}(M^{8})$ without reconstruction formulae.

For larger systems, the over-correlation of ST-J₂AGP near the recoupling regime worsens while the linear variational theory remains well-behaved. Thus, although ST-J₂AGP has a lower cost than other AGP-based methods developed recently in our group, neither var-J₂AGP nor ST-J₂AGP yield as accurate energies for the pairing Hamiltonian. While correlators build from particle-hole excitations generally benefit from being put in an exponential operator, it is far from clear that the same is true of correlators built from Hilbert space Jastrow operators. Thus in the next section we discuss an exponential ansatz that we build from different generators.

Recall the definition of the anti-Hermitian pair-hopper operator Khamoshi et al. (2020)

where $\tau_{pq}$ is the amplitude and is antisymmetric. $\boldsymbol{\mathcal{T}}$ is derived from

where $\mathbf{K}_{pq}|\text{AGP}\rangle=0$ is the two-body killer of AGP Henderson and Scuseria (2019). Its exact expression is

Notice that $\mathrm{e}^{\boldsymbol{\mathcal{T}}}|{\text{HF}}\rangle$ is equivalent to paired unitary CC with generalized indices (UGCC) Nooijen (2000); Stein et al. (2014); Lee et al. (2019).

Define the unitary pair-hopper ansatz as $\mathrm{e}^{\boldsymbol{\mathcal{T}}}|\text{AGP}\rangle$ which we denote as uP-AGP. We want to solve for the amplitudes by minimizing the energy, that is

However, since the BCH expansion of the unitary transformed Hamiltonian does not truncate naturally, we choose to truncate the expansion at the highest order that gives us an affordable scaling, i.e. $\mathcal{O}(M^{6})$. Again, using the reconstruction formulae, the expectation value of all two-body or higher density matrices can be accessed at $\mathcal{O}(1)$ cost per element. This means the asymptotic scaling of the commutator expansion is dominated by the many-body contractions. For example, for a two-body Hamiltonian, the first commutator, $\langle{[{\boldsymbol{H}},{\boldsymbol{\mathcal{T}}}]}\rangle$ contains three-body contractions; the second commutator, $\langle{[{[{\boldsymbol{H}},{\boldsymbol{\mathcal{T}}}]},{\boldsymbol{\mathcal{T}}}]}\rangle$ contains four-body contractions, and so on. By truncating the expansion at the 4th commuter, we obtain a scaling of $\mathcal{O}(M^{6})$ for the energy. In other words, we can state

such that addition of every term in the expansion multiplies the leading scaling by a factor of $M$. Obtaining the equations above analytically is cumbersome but not prohibitive. We use our home-built algebraic manipulator software, drudge Zhao (2018), to generate this and the analytic gradient of the energy. As a result, we can compute both $E(\tau)$ and $\nabla E(\tau)$ at $\mathcal{O}(M^{6})$ cost albeit with a relatively large prefactor.

We minimize the energy as in Eq. (40) using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm Byrd et al. (1995) where we provide the gradients using our analytical expressions. We could alternatively find the solutions to $\nabla E(\tau)=0$ using non-linear root-finding algorithms, but we find that the former is typically more robust when an accurate initial guess is provided.

Because truncating the commutator expansion yields a non-variational energy expression, we must take care with the minimization of the energy. This becomes particularly important in larger systems for which the numerical algorithms may go below the exact energy and eventually blow up. Fortunately, this failure mode is easy to identify as it is accompanied by extremely large amplitudes in $\tau$. To avoid this problem, we introduce an $L_{2}$ regularization Nocedal and Wright (2006)

where $\lambda>0$ is the regularization parameter and is sufficiently small to allow a robust convergence. In effect, the regularization penalizes the energy optimization in proportion to the magnitude of the amplitudes, thereby preventing it from blowing up. Although this compromises the accuracy of final energy, we gain a smoother convergence. To find $\lambda$, we pick a small value e.g. $10^{-3}$, converge to an optimal solution for $\tau$, and use this solution as the initial guess of another optimization in which a smaller $\lambda$ is used. The first initial guess for $\tau$ can be set to zero. We repeat this process until we can no longer converge or the change in energy is negligible. The convergence threshold of the gradient is set to $10^{-8}$, and we have not observed significant dependence of the final results on the initial guess for $\tau$ or $\lambda$ in our calculations.

In Fig. 2 we plot the results for the pairing Hamiltonian in a small system in which we can compare the accuracy of our truncated commutator expansion with that of exact uP-AGP which we obtain from a full configuration interaction (FCI) code. As we can see in the plot, the agreement between the exact and truncated uP-AGP is excellent on the attractive side while at the same time being slightly more accurate than the linear wavefunction $\boldsymbol{\mathcal{T}}|\text{AGP}\rangle$ denoted as P-CI AGP. On the repulsive side, the results for the truncated uP-AGP is comparable to P-CI and becomes more accurate only at sufficiently small $G/G_{c}$.

For benchmarking purposes, we also plot HF-based unitary coupled cluster (UCC) Bartlett et al. (1989); Taube and Bartlett (2006) as well as that with generalized indices, UGCC, in Fig. 2. The results are obtained from a FCI code. As expected UCC works well in the repulsive regime, but it breaks down as $G$ approaches $G_{c}$. On the other hand, UGCC is remarkably more accurate than UCC, but it too eventually breaks down for large $G/G_{c}$. In contrast, our AGP based methods are well behaved in all correlation regimes.

The results for a larger system is show in Fig. 3 where we compare the accuracy of the truncated uP-AGP with ST-J₂AGP and P-CI. Although all of these methods have the same number of parameters, they clearly differ considerably in their accuracy. The observation that uP-AGP could be more accurate than P-CI in principle is a testament to its non-linear nature. Previously, it was shown that CI methods based on the generators of the algebra, namely J₂-CI, P-CI, and K-CI, yield identical results Dutta et al. (2020). The other advantage of the non-linear model is that it is agnostic to the linear dependencies of the ansatz. This is true for both ST-J₂AGP and the unitary pair-hopper. Although removing the linear dependencies by diagonalizing the metric can yield faster convergence, our numerical results show very little difference in the final energy.

The downside of the truncated uP-AGP is that it has a relatively large prefactor which is a result of the brute-force expansion of the BCH formula. In the next section we apply canonical transformation theory as another means for approximating the unitary transformed Hamiltonian.

In canonical transformation (CT) theory, as formulated originally by Yanai and Chan Yanai and Chan (2006, 2007), the commutator expansion of a unitary transformed Hamiltonian,

where $\boldsymbol{A}$ is an anti-Hermitian operator, is recursively summed by systematically replacing all higher-body operators obtained from each commutator by an approximate, lower-body operator. For example, if $\boldsymbol{A}$ and $\boldsymbol{H}$ are each two-body operators, one approximates $[{\boldsymbol{H}},{\boldsymbol{A}}]\rightarrow[{\boldsymbol{H}},{\boldsymbol{A}}]_{(1,2)}$ wherein the three-body operator resulting from $[{\boldsymbol{H}},{\boldsymbol{A}}]$ is approximated in terms of one- and two-body operators. This process yields an effective Hamiltonian, $\bar{\boldsymbol{H}}\approx\mathrm{e}^{-\boldsymbol{A}}\boldsymbol{H}\mathrm{e}^{\boldsymbol{A}}$, that, in its simplest form, can be written as

where we have retained up to two-body terms in the effective Hamiltonian. Eq. (44) in particular can be summed recursively by

where $\bar{\boldsymbol{H}}^{(0)}=\boldsymbol{H}$, and the sum is continued until the change in the parameters of the effective Hamiltonian is sufficiently small. The energy is then defined as $E=\langle{\Psi}|\bar{\boldsymbol{H}}|{\Psi}\rangle$ where $|{\Psi}\rangle$ is some mean-field reference. The minimization of the energy can be carried out by solving CC-style amplitude equations $\langle{\Psi}|\boldsymbol{\hat{A}}^{\dagger}_{i}\bar{\boldsymbol{H}}|{\Psi}\rangle$ which can be rearranged to give

where $\hat{\boldsymbol{A}_{i}}$ are the operator parts of $\boldsymbol{A}$.

In essence, the goal of CT theory is to make the effective Hamiltonian more and more diagonal by successively approximating higher order excitations. It can be viewed as a kind of renormalization group approach and is rooted in an earlier work by White White (2002) as well as the flow renormalization approach by Wegner Wegner (1994) and Glazek and Wilson Glazek and Wilson (1994). In this section, however, we very closely follow the formalism presented in Ref. Yanai and Chan (2006, 2007); Neuscamman et al. (2009, 2010); Yanai et al. (2012) by Yanai, Chan, and Neuscamman.

The chief difference between our approach and those of earlier papers is in the operator decomposition or downfolding scheme. Instead of invoking the cumulant decomposition of RDMs Kutzelnigg and Mukherjee (1999) as a means to get lower-body operators, we use AGP and our reconstruction formulae. Just as in cumulant decomposition, the reconstruction formulae are statements about the RDMs, but we can extended them to operators as an approximation.

For example, consider the three-body irreducible density matrix $Z^{(3,3)}_{pqr}$. From the reconstruction formulae Khamoshi et al. (2019), it follows that

where $\Lambda_{qp}=2\eta_{q}^{2}/(\eta_{q}^{2}-\eta_{p}^{2})$. For CT, we propose to approximate $\mathbf{N}_{p}\mathbf{N}_{q}\mathbf{N}_{r}$ by removing the expectation values to get

Note that since the reconstruction formulae apply only to irreducible RDMs, we limit our downfolding to irreducible operators, i.e. $\mathbf{P}^{\dagger}_{p}\ldots\mathbf{N}_{q}\ldots\mathbf{P}_{r}\ldots$ where all indices are different. Indeed, all reducible operators (RDMs) can be written as a sum of irreducible ones.

To carry out the CT procedure, we take the unitary pair-hopper $\mathrm{e}^{\boldsymbol{\mathcal{T}}}$ ansatz to be the generator of our transformation, and we seek an effective Hamiltonian of the form

where $w_{pq}$ is symmetric and $v_{pq}$ is Hermitian. Note that Eq. (49) is the general form of a two-body Hamiltonian in the seniority-zero space Henderson et al. (2015). In order to get a relatively more accurate approximation in CT, we would like to delay the downfolding as much possible. To this end, following Ref. Neuscamman et al. (2010), we split $\bar{\boldsymbol{H}}=\bar{\boldsymbol{H}}_{1}+\bar{\boldsymbol{H}}_{2}$ such that $\bar{\boldsymbol{H}}_{1}$ and $\bar{\boldsymbol{H}}_{2}$ are the one- and two-body parts of the effective Hamiltonian respectively, and the downfolding is delayed until four-body operators appear. As such, one obtains

Similarly, to get accurate residual equations, we split $R_{pq}=R_{1,pq}+R_{2,pq}$ and evaluate the residuals associated with the one- and two-body parts recursively

In both Eq. (50) and Eq. (51) we retain all two body-terms that naturally appear in the commutator expansion and absorb their coefficients into $w_{pq}$ and $v_{pq}$ in Eq. (49). All higher-body operators are downfolded directly into one-body with the intent of making the Hamiltonian more diagonal. The cost of evaluating the effective Hamiltonian and the residuals is $\mathcal{O}(M^{4})$.

Having the effective Hamiltonian and the residual equations, one can solve the amplitudes using a non-linear root-finding algorithm. In so doing, we can provide an approximation to the Jacobian matrix as follows

The cost of building the Jacobian is also $\mathcal{O}(M^{4})$ by using the reconstruction formula.

We apply the CT as described above to the pairing Hamiltonian Eq. (2). As encountered in other implementations of CT, we find that the optimization over the amplitude equations is ill-conditioned Neuscamman et al. (2009, 2010); Yanai et al. (2012). These numerical difficulties are attributed to near-zero eigenvalues of the Jacobian and are said to be similar to the intruder states in the second order perturbation theory Neuscamman et al. (2010). Several techniques have been proposed to mitigate these difficulties Neuscamman et al. (2010); Yanai et al. (2012). In our implementation, we find that shifting of the amplitude equations Yanai et al. (2012)

where $\mu,\nu$ are flattened indices, to be the most effective. This is analogous to the regularization introduced in Eq. (42). The shifting parameter $\lambda>0$ can be chosen ad hoc Yanai et al. (2012). While we may follow the same steps as in the previous section to find optimal values for $\lambda$, for convenience, we picked some sufficiently small values with which the equations can converge reliably.

The results for 40 levels, 20 pairs is shown in Fig. 4. We picked $\lambda=1$ for $G>0$ and a larger $\lambda=2.5$ for $G<0$. CT with the unitary pair-hopper on AGP (CT-uP-AGP) is considerably less expensive than the truncated uP-AGP as well as the linear P-CI AGP. However, the approximation made to higher order excitations compromises the accuracy of CT-uP-AGP. Systematic improvements can be made to CT-uP-AGP by delaying the downfolding even further, but we leave these improvements for future work.