Projected Coupled Cluster Theory

TL;DR

This paper introduces projected coupled cluster theory, combining symmetry projection with coupled cluster methods to better handle strong correlations, demonstrated on molecules and the Hubbard model.

Contribution

It develops a disentanglement framework for applying symmetry projection to coupled cluster wavefunctions, improving accuracy in strongly correlated systems.

Findings

Spin projection enhances unrestricted coupled cluster accuracy.

Projected coupled cluster restores quantum numbers effectively.

Energy convergence is slower than projected Hartree-Fock in the thermodynamic limit.

Abstract

Coupled cluster theory is the method of choice for weakly correlated systems. But in the strongly correlated regime, it faces a symmetry dilemma, where it either completely fails to describe the system, or has to artificially break certain symmetries. On the other hand, projected Hartree-Fock theory captures the essential physics of many kinds of strong correlations via symmetry breaking and restoration. In this work, we combine and try to retain the merits of these two methods by applying symmetry projection to broken symmetry coupled cluster wavefunctions. The non-orthogonal nature of states resulting from the application of symmetry projection operators furnishes particle-hole excitations to all orders, thus creating an obstacle for the exact evaluation of overlaps. Here we provide a solution via a disentanglement framework theory that can be approximated rigorously and systematically. Results of projected coupled cluster theory are presented for molecules and the Hubbard model, showing that spin projection significantly improves unrestricted coupled cluster theory while restoring good quantum numbers. The energy of projected coupled cluster theory reduces to the unprojected one in the thermodynamic limit, albeit at a much slower rate than projected Hartree-Fock.

Equations

$$|\Psi ⟩={\mathrm{e}}^U|\varphi ⟩$$

$$H{\mathrm{e}}^U|\varphi ⟩=E{\mathrm{e}}^U|\varphi ⟩$$

$$E$$

$$0$$

$${\mathrm{e}}^{-U}H{\mathrm{e}}^U|\varphi ⟩=\bar{H}|\varphi =E|\varphi ⟩$$

$$\bar{\mathcal{O}}={\mathrm{e}}^{-U}\mathcal{O}{\mathrm{e}}^U.$$

$$E$$

$$0$$

$$⟨L|=⟨\varphi |\left(1+Z\right)$$

$$⟨\varphi |\left(1+Z\right)\bar{H}=E⟨\varphi |\left(1+Z\right).$$

$$\mathcal{E}$$

$$0$$

$$⟨\varphi |{\mathrm{e}}^{-U}$$

$$⟨\mu |{\mathrm{e}}^{-U}$$

$$|\Psi ⟩=P|\psi ⟩$$

$$P=\frac{1}{V_{\Omega }}\int \mathrm{d}\Omega R(\Omega )w(\Omega )$$

$$E=\frac{⟨\chi |P^{†}HP|\psi ⟩}{⟨\chi |P^{†}P|\psi ⟩}=\frac{⟨\chi |PH|\psi ⟩}{⟨\chi |P|\psi ⟩},$$

$$\mathcal{N}(\Omega )$$

$$\mathcal{H}(\Omega )$$

$$E$$

$$P_S={\int }_0^{2\pi }\frac{\mathrm{d}\alpha }{2\pi }{\int }_0^{\pi }\frac{\sin (\beta )\mathrm{d}\beta }{2}{\int }_0^{2\pi }\frac{\mathrm{d}\gamma }{2\pi }R(\alpha ,\beta ,\gamma )$$

$$R(\alpha ,\beta ,\gamma )={\mathrm{e}}^{-\mathrm{i}\alpha S_z}{\mathrm{e}}^{-\mathrm{i}\beta S_y}{\mathrm{e}}^{-\mathrm{i}\gamma S_z}.$$

$$P_S=P_{s_z=0}{\int }_0^{\pi }\frac{\sin (\beta )\mathrm{d}\beta }{2}{\mathrm{e}}^{-\mathrm{i}\beta S_y}{P}_{s_z=0}$$

$$P_{s_z=0}|\psi ⟩$$

$$⟨\chi |P_{s_z=0}$$

$$\mathcal{N}(\beta )$$

$$\mathcal{H}(\beta )$$

$$E$$

$$|\Psi ⟩=P{\mathrm{e}}^U|\varphi ⟩.$$

$$E$$

$$0$$

$$PH{\mathrm{e}}^U|\varphi ⟩=EP{\mathrm{e}}^U|\varphi ⟩$$

$${\mathrm{e}}^{-U}PH{\mathrm{e}}^U|\varphi ⟩=E{\mathrm{e}}^{-U}P{\mathrm{e}}^U|\varphi ⟩$$

$$\bar{P}\bar{H}|\varphi ⟩=E\bar{P}|\varphi ⟩.$$

$$E$$

Full text

I Introduction

Single-reference coupled cluster theory (CC)Coester (1958); Čížek (1966); Paldus and Li (1999); Bartlett and Musiał (2007); Shavitt and Bartlett (2009) is often considered the gold standard of quantum chemistry. While relatively expensive, coupled cluster with single and double excitations (CCSD) is affordable on systems of moderate to large size, and adding a perturbative correction for triple excitations leads to results which, for weakly correlated system, are generally within chemical accuracy at polynomial cost – $\mathcal{O}(N^{6})$ for CCSD and $\mathcal{O}(N^{7})$ with perturbative triple excitation corrections, where $N$ is some measure of system size.

Unfortunately, single-reference coupled cluster’s reputation is tarnished somewhat by its inability to describe strongly correlated problems with more than a few strongly correlated electrons. Active space coupled cluster methods are some help here, but are not a panacea, and in general systems with many strongly correlated electrons cannot readily be described by symmetry-adapted coupled cluster techniques.

To remedy this problem, we are often forced to use a broken-symmetry mean-field reference such as unrestricted Hartree-Fock (UHF), so that both the reference determinant and indeed the coupled cluster wave function lack some or perhaps all of the symmetries of the exact wave function. In the thermodynamic limit, this symmetry breaking is in fact real and physical, but for finite systems it is artificial and should be avoided. The practical result is that while symmetry-broken coupled cluster may be energetically accurate, the loss of symmetry can yield poor results for properties other than total energies.

At the mean-field level, one can obtain the energetic benefits of symmetry breaking in a symmetry-adapted picture by using the projected Hartree-Fock (PHF) method.Löwdin (1955); Ring and Schuck (1980); Blaizot and Ripka (1985); Schmid (2004); Jiménez-Hoyos et al. (2012) The idea of PHF is simple enough. One allows symmetry breaking in the mean-field, and then projects the broken-symmetry determinant back onto the correct symmetry. This can be accomplished by writing the PHF wave function as a relatively short linear combination of degenerate and non-orthogonal broken symmetry determinants obtainable from one another by a symmetry rotation operator. PHF is a black-box technique without the necessity of picking active orbitals. Furthermore, it possesses an underlying reference determinant – in fact, a whole manifold of them – which in principle allows for combination with single-reference coupled cluster theory.

Ideally, we would like to combine PHF and coupled cluster, and we have made several attempts to do so.Qiu, Henderson, and Scuseria (2016, 2017); Degroote et al. (2016); Wahlen-Strothman et al. (2016); Gomez, Henderson, and Scuseria (2017); Hermes, Dukelsky, and Scuseria (2017) In our previous work, we have written the projected Hartree-Fock wave function in terms of particle-hole excitations acting on a symmetry-adapted reference, in which form we can readily combine PHF with fairly traditional symmetry-adapted coupled cluster theory. In this work, we take a complementary approach and work in the symmetry-broken picture. The idea is simply to carry out the symmetry projection of a broken-symmetry coupled cluster wave function. We will focus on spin symmetry in this manuscript, because it is the symmetry that spontaneously breaks in molecular systems, but the basic framework is more general. The main challenge for our projected coupled cluster theory is the non-orthogonality to which we have previously alluded. Orthogonal determinants, differing by particle-hole excitations, lose their orthogonality after the symmetry rotation, but many of the simplifications of traditional coupled cluster theory rely on this orthogonality. Practical approximations require some way of either evaluating or truncating the overlaps between nonorthogonal states. In this manuscript, we do so by what we call disentangled cluster operators. This formalism permits us to work with orthogonal particle-hole excitations and thereby truncate projected coupled cluster in a systematic way due to the decay of the disentangled cluster operators with excitation rank.

The projected coupled cluster theory (PCC) introduced in this manuscript has some nice features. If the cluster operators are set to zero, it reduces to PHF. If the projection operator is disregarded, it reduces to broken symmetry coupled cluster theory. In between, we have a theory with particle-hole excitations to all orders; this is the result of disentangling the action of symmetry projection and broken-symmetry coupled cluster, but is a model that can be rigorously approximated by truncation, very much in the spirit of traditional coupled cluster theory. Just as the size-extensive component of the PHF energy is that of the symmetry-broken mean-field, the size-extensive component of the PCC energy is that of the symmetry-broken coupled cluster. However, PCC reduces to broken symmetry coupled cluster at a much slower rate than PHF reduces to broken symmetry Hartree-Fock. Moreover, it should not be forgotten that the PCC wave function has correct symmetry. While the size extensive energetic component of PHF is the same as broken-symmetry mean field (though there is also a size-intensive component independent of system size in the thermodynamic limitJiménez-Hoyos et al. (2012)), the wave functions are different.

Previously, Duguet proposed a symmetry broken and restored coupled cluster theory.Duguet (2015) We have borrowed some concepts and nomenclature from his paper, but these two approaches are very different, as we shall examine below in more detail.

II Background

Before we can discuss the basic formulation of our symmetry-projected coupled cluster theory, we briefly review both traditional coupled cluster theory and symmetry projection to establish some basic concepts. Readers familiar with these ideas may freely skip the following two subsections, though we do use them to establish our notation.

II.1 Traditional Coupled Cluster Theory

In coupled cluster theory, the ground state wave function $|\Psi\rangle$ is approximated by the exponential ansatz

[TABLE]

where $|\phi\rangle$ is a single determinant and where $U$ is an excitation operator. Inserting this wave function ansatz into the Schrödinger equation yields

[TABLE]

and leads to energy and amplitude equations

[TABLE]

where $|\mu\rangle$ is an excited determinant. So-called unlinked terms appear in the foregoing equations but cancel out.Hirata (2011)

Coupled cluster theory is not conventionally described as we have done. More commonly, one multiplies both sides of Eqn. 2 by $\exp(-U)$ to get a similarity-transformed Schrödinger equation:

[TABLE]

where generically the operator $\bar{\mathcal{O}}$ means

[TABLE]

Now the energy and amplitudes defining $U$ are obtained from

[TABLE]

The energy and amplitude equations of Eqn. 6 together imply that $|\phi\rangle$ is a right-hand eigenstate of $\bar{H}$. Since $\bar{H}$ is non-Hermitian, its left-hand eigenstate $\langle L|$ is not the adjoint of its right-hand eigenstate, and we can parametrize it in a configuration interaction-like way, as

[TABLE]

where $Z$ is an excitation operator acting to the left, and thus a de-excitation operator acting to the right. As the left-hand eigenstate, we can solve for $Z$ by demanding that

[TABLE]

With a little effort, one can show that Eqns. 6 and 8 together can be summarized as

[TABLE]

where $Z_{\mu}$ and $U_{\mu}$ are individual amplitudes in the operators $Z$ and $U$. If the amplitudes defining $U$ satisfy Eqn. 9b, then for any amplitudes $Z$, the two energies $E$ and $\mathcal{E}$ are identical; otherwise, they may differ significantly. We would argue that $\mathcal{E}$ is the more correct definition of the coupled cluster energy in general.

We have already noted that Eqns. 6 and 9 yield equivalent formulations of coupled cluster theory. Both are equivalent to the non-similarity-transformed approach of Eqns. 3, as can be shown by making use of the facts that

[TABLE]

where the coefficients $c_{\mu\nu}$ depend on the amplitudes defining $U$ and where the sum on $\nu$ includes determinants of equal or lower excitation level to $\langle\mu|$, including the ground state determinant $\langle\phi|$.

Thus far, we have made no approximations; when $U$ contains all possible excitations, the theory as we have outlined it is exact. In practical calculations, $U$ and $Z$ must be truncated, and the resultant theories are named according to the excitation levels retained. For example, the simplest theory is coupled cluster doubles (CCD) and retains only the double-excitation parts of $U$ and $Z$, where CCSD keeps both the single- and double-excitation parts of $U$ and $Z$. Typically we truncate $Z$ at the same level to which we truncate $U$.

II.2 Symmetry Projection

The basic idea of symmetry projection is straightforward. From a broken symmetry wave function $|\psi\rangle$ the symmetry adapted component is simply

[TABLE]

where the projector $P$ is Hermitian and idempotent, and commutes with the Hamiltonian.

Some important symmetries, such as number and spin, are continuous. For these symmetries, perhaps the simplest way to write the projector is as an integral over symmetry-generated transformations:

[TABLE]

where $R(\Omega)$ is a one-body rotation operator, $w(\Omega)$ is a weight which depends on the eigenvalue being projected onto, and $V_{\Omega}=\int\mathrm{d}\Omega$ is the volume of the space being integrated over. If the symmetry is discrete, the integration is replaced by a summation. In practice, the integral is replaced by a weighted sum over a grid even for a continuous symmetry such as spin.

Given a projected wave function, one has several options for extracting the energy. In PHF we use a simple expectation value, but here we generalize to a biorthogonal expectation value as is used in coupled cluster theory (c.f. Eqn. 9a) because we will need this form for the projected coupled cluster energy. In this kind of biorthogonal approach, one might write

[TABLE]

where $\langle\chi|$ is some other broken symmetry wave function and where we have used the properties of the projection operator. In view of the integral form of the projection operator, the energy can be expressed in terms of reduced norm and Hamiltonian kernels $\mathcal{N}(\Omega)$ and $\mathcal{H}(\Omega)$:

[TABLE]

Note that $\langle\chi|R(\Omega)\,H|\psi\rangle$ is proportional to the norm kernel.

Our primary interest in this work is in spin projection. While the general case of spin projection is slightly more complicated than what we have outlined above,Jiménez-Hoyos et al. (2012) for the special case of spin projection onto a singlet state we can write a projection-like operator

[TABLE]

where the rotation operator is

[TABLE]

This projector simplifies to

[TABLE]

where $P_{s_{z}=0}$ is the projector onto eigenstates of $S_{z}$ with eigenvalue zero.

When working in the unrestricted framework in which wave functions are eigenfunctions of $S_{z}$ but not of $S^{2}$, the projector simplifies somewhat. In particular, we can use the facts that

[TABLE]

for spin unrestricted wave functions $|\psi\rangle$ and $\langle\chi|$ to write the norm kernel, Hamiltonian kernel, and projected energy of Eqn. 14 as simply

[TABLE]

III Projected Coupled Cluster Theory

Now that we have reviewed the background material, we are in a position to introduce our spin-projected unrestricted coupled cluster theory. In this section we make no approximations other than the truncation of the cluster operator. As we shall see, the exact theory is computationally too cumbersome for practical use, and we will introduce computationally tractable approximate versions in Sec. IV. These approximations are motivated by the smallness of the cluster amplitudes in broken symmetry coupled cluster and afford truncation in the same spirit as traditional single-reference coupled cluster theory.

III.1 Theory

The basic idea of our projected coupled cluster (PCC) theory is simple: the symmetry adapted wave function $|\Psi\rangle$ is obtained by projecting a broken-symmetry coupled cluster wave function as

[TABLE]

Just as for traditional coupled cluster, the PCC energy and amplitude equations are derived by inserting the wave function ansatz into the Schrödinger equation to get

[TABLE]

where $|\mu\rangle$ is an excited determinant. We have used the fact that $H$ and $P$ commute.

Alternatively, one could set up a similarity-transformed approach which requires similarity transformations of both the Hamiltonian and the projection operator. Starting from the Schrödinger equation

[TABLE]

we multiply on the left by $\exp(-U)$ to obtain

[TABLE]

or equivalently,

[TABLE]

Recall that $\bar{H}$ is the similarity-transformed Hamiltonian (c.f. Eqn. 5) while $\bar{P}$ is the similarity-transformed projection operator. The energy and the amplitudes defining $U$ are then

[TABLE]

Furthermore, by parameterizing the left hand state as in traditional coupled cluster theory, one can introduce the PCC analogue of the traditional coupled cluster energy functional, in terms of which the energy and amplitude equations are

[TABLE]

Just as in traditional coupled cluster, one can show that all three formulations are equivalent, though it is important to note that, also as in traditional coupled cluster, we have $\mathcal{E}=E$ only when the cluster operator $U$ satisfies the amplitude equations.

Traditional coupled cluster theory finds it convenient to work with the similarity-transformed approach because the similarity transformation of the Hamiltonian naturally terminates at $\mathcal{O}(U^{4})$, a consequence of the Hamiltonian being a two-body operator. The projection operator, however, is an $N$-body operator, and its similiarity transformation therefore does not terminate in a convenient way. Accordingly, we find it simpler to work in the formalism of Eqn. 21. Of course one could truncate the commutator expansion of $\bar{P}$ manually, though so far we have considered only the approximation $\bar{P}\approx P$. This approximation leads to results of significantly lower quality, so we will not discuss them.

While the similarity transformation of $P$ is cumbersome and essentially precludes the use of the similarity-transformed approach as it is traditionally expressed, we will find the need to evaluate the biorthogonal energy $\mathcal{E}$ from time to time. To facilitate this, we define

[TABLE]

where $\tilde{Z}$ is a modified excitation operator when acting to the left and contains a component which creates no excitations. This is essentially a consequence of Eqn. 10. To give a concrete example, when $U$ and $Z$ contain only single and double excitation operators, we may write

[TABLE]

where the subscripts “$\mathrm{cd}$” and “$\mathrm{cs}$” respectively mean the connected deexcitation part and the connected scalar part. In terms of $\tilde{Z}$ we may write the biorthogonal energy as

[TABLE]

where $\tilde{Z}_{\nu}$ are the amplitudes defining the operator $\tilde{Z}$ and where the summation over $\nu$ includes the reference determinant $\langle\phi|$ as well as excited determinants.

Our basic formulation is agnostic as to the precise way in which the symmetry projection is carried out. In practice, we find it most convenient to use the integral representation outlined in Sec. II.2. To do so, we need only define the requisite integral kernels, which are

[TABLE]

From these kernels, we can obtain all quantities needed to evaluate the projected energy and amplitude equations:

[TABLE]

We should make an important caveat at this point. In traditional coupled cluster theory the exponential can be truncated in a natural way. Thus, for example, the CCSD equations require only up to the quadruple excitation part of $\exp(U)$. This is so because we consider states $\langle\mu|$ which are no more than doubly-excited, and the Hamiltonian (as a two-body operator) cannot connect a doubly-excited bra state with more than a quadruply-excited ket state. Unfortunately, in projected coupled cluster the exponential cannot be so conveniently truncated, due to the presence of the rotation operator $R(\Omega)$. If we imagine defining rotated states

[TABLE]

then even when $\langle\mu|=\langle\mu(\Omega=0)|$ is doubly-excited with respect to the ground state, $\langle\mu(\Omega\neq 0)|$ may have components which are very highly excited indeed. Section IV discusses how we circumvent this difficulty.

III.2 Results

Thus far we have been as general as possible. Our PCC theory can be combined with any symmetry projection, and if the cluster operator is not truncated it provides the exact ground state wave function. Here, we provide results for the special case of spin-projected unrestricted coupled cluster with single and double excitations (SUCCSD), and consider only projection onto a spin singlet. We have adapted a full configuration interaction code to do SUCCSD calculations, which limits our exact SUCCSD to small systems. A production version of the theory requires additional approximations which we introduce in Sec. IV and which allow for a computationally efficient implementation; our exact results here thus provide a benchmark to compare against. To ensure accuracy, we have used exact grids to integrate the various kernels where we have as many grid points as basis functions. In practice the number of grid points needed for spin-projected UHF (SUHF) scales roughly as the square root of the number of basis functions and we expect fewer grid points are needed for SUCCSD.

In this section we also wish to address the importance of obtaining the amplitudes $U_{\mu}$ and potentially $Z_{\mu}$ from our PCC equations rather than from traditional coupled cluster. In analogy with the literature on projected Hartree-Fock theory, we will refer to a variation after projection (VAP) approach in which the amplitudes solve the SUCCSD equations, where a projection after variation (PAV) approach means that we first solve the traditional unrestricted CCSD (UCCSD) equations and simply evaluate the projected energy. Without further qualification, by the PAV energy we mean the energy expression of Eqn. 21a, but we may also refer to the linear response PAV (LR-PAV) in which we use the biorthogonal expectation value of Eqn. 26a. Recall that these two energy expressions yield different results when using UCCSD amplitudes but the same result if one solves for the $U$ amplitudes in the presence of the projection operator. In other words, the VAP formulation of SUCCSD yields the same energy from the two different energy formulae.

Bearing all this in mind, let us first consider results in the Hubbard model Hamiltonian.Hubbard (1963) The Hubbard model describes electrons on a lattice, and the Hamiltonian is simply

[TABLE]

where $\mu$ and $\nu$ are lattice sites and the notation $\langle\mu\nu\rangle$ means that we include only sites connected in the lattice. As the relative interaction strength $U/t$ increases, the system becomes more and more strongly correlated. Our calculations use a one-dimensional lattice in which adjacent sites are connected. We impose periodic boundary conditions, which is of more general interest and for which exact results are readily available through the Lieb-Wu algorithm.Lieb and Wu (1968)

Figure 1 shows results for the half-filled lattice with six and ten sites. We obtain the reference determinant $|\phi\rangle$ from a variation after projection spin-projected UHF (SUHF) calculation, so that $|\phi\rangle$ is symmetry broken for all $U/t$. One can see that the spin projection improves significantly upon UCCSD, regardless of whether the amplitudes solve the SUCCSD equations (VAP) or the UCCSD equations (PAV). Unsurprisingly, our overall best results are obtained from the VAP approach. In the strongly correlated limit, VAP and LR-PAV are comparable, and LR-PAV is in general superior to PAV without including the $Z$ amplitudes. We do not know why the LR-PAV results have a relatively pronounced peak for small $U/t$.

Let us now turn to a few molecular examples. The first example is four hydrogen atoms placed on a circle of radius 3.284 Bohr,Van Voorhis and Head-Gordon (2000) as depicted in Fig. 2. As the bond angle $\theta$ approaches $90^{\circ}$, the system becomes more strongly correlated as the ground and first excited states become nearly degenerate (and exactly degenerate precisely at the high symmetry point). Where restricted CCSD (RCCSD) has a cusp and incorrectly predicts the existence of a minimum at $\theta=90^{\circ}$, UCCSD remedies these failures at the cost of correct spin symmetry. Moreover, UCCSD undercorrelates somewhat. In contrast, SUCCSD is virtually atop the exact result, though we have been able to compute only the PAV result because the cc-pVDZ basis set which we use is too large for our full configuration interaction code. Incorporating linear response adds virtually nothing in this case.

Finally we consider the dissociation of N2 in the cc-pVDZ bas set, as shown in Fig. 3. Restricted CCSD overcorrelates badly at the dissociation limit and has an unphysical bump in the potential energy curve. Adding perturbative triples to obtain RCCSD(T) makes the situation even worse (not shown). If one uses an unrestricted reference the results are much improved. Results are better yet with SUCCSD, though again we have not been able to carry out the VAP calculations. For small $R$, the UHF reference does not break symmetry, so symmetry projection has no effect. This causes the break in the LR-PAV curve. The LR-PAV curve is smoothed by using the SUHF reference instead (not shown), though overall the reference-dependence is weak.

We should make one small caveat about our N2 results. Because the system was beyond the scope of our full configuration interaction code, we were unable to evaluate the Hamiltonian and norm kernels exactly. Accordingly, we truncated the exponential of the cluster operator in both and retained terms only through $\mathcal{O}(U_{2}^{3})$. While exact results would require us to keep terms through $\mathcal{O}(U_{2}^{7})$, we have found that the energy is in practice converged at the truncation we have shown, presumably because there are six strongly correlated electrons. We have kept the $U_{1}$ terms to infinite order, as they can be absorbed into transforming the Hamiltonian, discussed in detail in Appendix B.

Nonetheless, the inability to exactly compute the Hamiltonian and norm kernels is a serious limitation, and simply truncating the exponential is not always practical. Accordingly, we now introduce our disentangled cluster formalism, which allows us to circumvent this limitation.

IV The Disentangled Cluster Formalism

As we have seen, the chief difficulty in evaluating the kernels needed to do PCC calculations is that in principle the exponential must be expanded to all powers of $U$, which is not generally feasible. We now turn our attention to overcoming this obstacle. We will need to make some preliminary simplifications first.

We begin by replacing the rotation operator, which in general is a unitary Thouless transformation,Thouless (1960) with a particle-hole-style Thouless transformation applied to the bra:

[TABLE]

where $V_{1}(\Omega)$ is a single de-excitation operator. Having done so, the norm and Hamiltonian kernels become

[TABLE]

where we have introduced the similarity-transformed Hamiltonian

[TABLE]

Since $V_{1}(\Omega)$ is a one-body operator, $\bar{H}_{V_{1}}(\Omega)$ contains no more than two-body operators (assuming that $H$ itself had only one- and two-body operators). Henceforth we will omit the explicit $\Omega$-dependence to reduce notational clutter.

In addition to the ground-state reduced kernels, we also need excited kernels $\mathcal{N}_{\mu}$ and $\mathcal{H}_{\mu}$. To obtain these, we first write

[TABLE]

where $\mathcal{Q}_{\mu}$ is an excitation operator when acting to the left. Having done so, we can then use

[TABLE]

where

[TABLE]

The excited norm and Hamiltonian kernels are thus

[TABLE]

We can evaluate the action of $\tilde{\mathcal{Q}}_{\mu}$ straightforwardly, because $\exp(V_{1})\,R^{-1}$ is just an orbital transformation operator, given in the molecular orbital basis as

[TABLE]

where $\mathbf{V}_{\mathrm{ov}}$ contains the amplitudes $v^{i}_{a}$ defining $V_{1}$ and where $\mathbf{R}$ is the matrix representation of the rotation operator. For example, the rotation operator $R$ for unrestricted spin projection is $\exp(-\mathrm{i}\,\beta\,S_{y})$ so the matrix $\mathbf{R}^{-1}$ is $\exp(\mathrm{i}\,\beta\,\mathbf{S}_{y})$ where $\mathbf{S}_{y}$ is the matrix representation of the $S_{y}$ operator.

As concrete examples, the singly- and doubly-excited kernels are

[TABLE]

We must make several comments to clarify the foregoing. First, we have followed the usual notation that indices $i$ and $j$ denote orbitals occupied in $|\phi\rangle$, while $a$ and $b$ denote orbitals unoccupied in $|\phi\rangle$ and $p$, $q$, $r$, and $s$ are general. Determinants $\langle\phi_{i}^{a}|$ and $\langle\phi_{ij}^{ab}|$ are singly- and doubly-excited, respectively. Our convention for matrices is that row (column) indices of the matrix $\mathbf{T}$ correspond to upper (lower) indices of $T_{p}^{q}$. Finally, we have introduced auxiliary quantities such as $\tilde{\mathcal{N}}_{p}^{q}$. These are

[TABLE]

These objects are sparse, and their non-zero elements can be worked out by considering the action of the fermionic creation and annihilation operators on the reference bra determinant $\langle\phi|$. We include the non-zero components in Appendix A.

IV.1 The Disentangled Cluster Operator

The remaining difficulty in evaluating the various kernels we need is handling the state $\mathrm{e}^{V_{1}}\,\mathrm{e}^{U}|\phi\rangle$. If we cannot do so, all of the foregoing is in vain. To accomplish this task, we now introduce our disentangled cluster approximation. We will specialize to the case of coupled cluster doubles where $U=U_{2}$, with $U_{2}$ a double-excitation operator. Singles need not be explicitly included because the single-excitation operator $U_{1}$ can be absorbed into transforming the Hamiltonian and symmetry operators, as we discuss in Appendix B. Higher excitations could be included in principle, and the basic framework generalizes to include them, but at increasing algebraic and computational complexity.

The basic idea is to rewrite $\mathrm{e}^{V_{1}}\,\mathrm{e}^{U_{2}}|\phi\rangle$ as an exponential of cluster operators, denoted by $W_{k}$, acting on $|\phi\rangle$. That is, when acting to the right, $\mathrm{e}^{U_{2}}$ is an excitation operator and $\mathrm{e}^{V_{1}}$ is a deexcitation operator, and their product, acting on the reference determinant, can be expressed purely in terms of excitations together with a normalization constant. One has simply

[TABLE]

where

[TABLE]

Here, the subscript “$\mathrm{ce}$” means only the connected excitation part is retained. The sum over $n$ in defining $W_{k}$ must include all possible excitation levels. Equation 45 is a rigorous result, proven in Appendix C, and is central to the remainder of this paper. We refer to it as the disentanglement equation, and to the $W_{k}$’s as the disentengled cluster operators. Equations 45a and 45b define $W_{k}$ for $k>=1$; $W_{0}$ is a number, given by

[TABLE]

and can be evaluated by our formula for $W_{k}$ with $k=0$ but replacing the connected excitation part with the connected scalar part.

Introducing the disentangled clusters makes it possible to work in an orthogonal framework. The necessary kernels can all be evaluated with low excitation rank $W_{k}$:

[TABLE]

Here, $C_{i}$ is the $i$-fold excitation part of $\mathrm{e}^{W}$:

[TABLE]

We need only up to $W_{2}$ to evaluate the energy without linear response. Evaluating the energy with linear response or solving the amplitude equations requires up to $W_{4}$.

IV.2 Evaluating the Disentangled Cluster Operators

At this point, we have everything we need for projected CCD or CCSD, provided only that we can obtain the disentangled cluster amplitudes defining $W$. We will discuss several strategies for doing so. We emphasize again that because $W$ depends on $V_{1}$ and $V_{1}$ depends on the integration variable $\Omega$, $W$ is $\Omega$ dependent, though we have suppressed this dependence for brevity.

IV.2.1 Truncated Disentangled Cluster Operators

The simplest approach is to presume that $U_{2}$ is small and truncate the summation defining $W$. If we truncate the sum over $U_{2}$ in Eqn. 45 to obtain amplitudes $W$ from which the kernels are extracted, however, we run into a problem: when the rotated state $\langle\phi|R$ has small overlap with $|\phi\rangle$, then $V_{1}$ is large and the series defining $W$ converges slowly, with terms having large amplitudes and alternating signs. For example, in the half-filled Hubbard Hamiltonian for large $U/t$ the mean-field ground state has a Néel structure in which each site is occupied by one electron and the electrons on neighboring sites point in opposite directions. Spin projection at $\beta=\pi$ amounts to flipping each spin, so $|\phi\rangle$ and $\langle\phi|R$ are strictly orthogonal; thus $V_{1}(\beta=\pi)$ cannot even be defined but the limit of $V_{1}$ as $\beta$ approaches $\pi$ is infinite.

We can do better by truncating the exponential itself. For example, we could expand

[TABLE]

where $W_{0}^{(n)}$ is the component of $W_{0}$ which has $n$ powers of $U_{2}$. Note that except for the norm, all integrals in Eqn. 31 are integrals of the product of the norm kernel and another kernel; in these cases, the product is truncated as a whole. For example, we would write

[TABLE]

where similarly $\mathcal{H}^{(n)}$ and $\mathcal{N}^{(n)}$ are portions of the kernel which contain $n$ powers of $U_{2}$.

The process of expressing a certain order of kernel in terms of $U_{2}$ and $V_{1}$ is tedious and error-prone. The derivation of equations is facilitated by an in-home algebra manipulator,Zhao and Scuseria and this part of code is generated by the accompanying automatic code generator.Zhao and Scuseria For this work, the energy functional without response is implemented up to fifth order in $U_{2}$, while the one with response is implemented up to third order. For the energy functional without response, the energy can be evaluated through $U_{2}^{3}$ at $\mathcal{O}(N^{6})$ cost, through $U_{2}^{4}$ at $\mathcal{O}(N^{7})$ cost, and through $U_{2}^{5}$ at $\mathcal{O}(N^{8})$ cost. The energy functional with response can be evaluated up through $U_{2}^{2}$ at $\mathcal{O}(N^{6})$ cost and through $U_{2}^{3}$ at $\mathcal{O}(N^{8})$ cost. All of these costs are per grid point in the numerical integration grid.

IV.2.2 Algebraic Equations for Disentangled Cluster Operators

A second approach to obtaining the disentangled cluster operators generalizes work of Jeziorski,Jeziorski and Moszynski (1993) that we rediscovered independently. Let us define the virtual occupation operator as $n_{V}=\sum_{a}c^{\dagger}_{a}c_{a}$. It can be shown that for a $k$th level excitation $U_{k}$, or a $k$th level deexcitation $V_{k}$,

[TABLE]

As a result,

[TABLE]

Applying these results to both sides of Eqn. 45 separately leads us to

[TABLE]

On the other hand,

[TABLE]

Together, these equations show that

[TABLE]

which implies

[TABLE]

Amplitudes of $W_{k}$ can be extracted through left-multiplication with excited determinants. For example

[TABLE]

or, more abstractly,

[TABLE]

Notice that $J$ contains one- and two-body parts only. Thus, the expressions for $W$ look like the traditional coupled cluster equations, with $H$ replaced by $J$, and accordingly, $W_{k}$ is coupled with up to $W_{k+2}$. Note also that because $W_{0}$ is simply a number and these algebraic equations are derived from commutators, we cannot use them to evaluate or even approximate $W_{0}$.

IV.2.3 Differential Equation for Disentangled Cluster Operators

For the special case of unrestricted spin projection, we can solve for the disentangled cluster operators $W$ by solving a set of differential equations with the spin rotation angle $\beta$ as the independent variable. By integrating this set of ordinary differential equations (ODEs), we can in principle obtain exact untruncated amplitudes, though in practice other approximations must be made.

To proceed, we need the derivative of $V_{1}$ with respect to $\beta$. We define

[TABLE]

where the subscript $d$ means only the deexcitation part is retained and where this defines $X$. We derive this equation in Appendix D.

Now we note that

[TABLE]

Differentiating with respect to $\beta$ leads to

[TABLE]

whence

[TABLE]

Given $W_{1}$, we can integrate this equation to obtain $W_{0}$.

Similarly, we have

[TABLE]

Differentiating both sides leads to

[TABLE]

Here, summation over repeated indices is implied. Notice that the first term on the right-hand-side is

[TABLE]

Then the disconnected parts on both sides cancel and we have simply

[TABLE]

We can obtain a similar result for double excitations. Starting from

[TABLE]

we obtain

[TABLE]

Let us note a few important features of these differential equations for the disentangled cluster operators. First, at $\beta=0$ we have $V_{1}=0$, so $W_{0}=0$, $W_{1}=0$, and $W_{2}=U_{2}$. Second, the derivative of $W_{k}$ contains $W_{k+1}$ but no higher excitations. We can thus approximate $W_{k+1}$ for some excitation level $k$ and from there solve for the lower excitation levels. For this reason we prefer the differential equations to the algebraic equations given in the previous section, as the latter require both $W_{k+1}$ and $W_{k+2}$ to obtain $W_{k}$; morever, the algebraic equations do not allow us to compute $W_{0}$.

Here, we consider several approximations. If we make the approximation $W_{3}=0$ and integrate to obtain $W_{0}$ to $W_{2}$, we call the resulting method SUCCSD(SD) where SUCCSD stands for spin-projected unrestricted CCSD and “(SD)” emphasizes that we retain through $W_{2}$ in the differential equation. Similarly, SUCCSD(SDT) means we set $W_{4}=0$ and solve for $W_{0}$ to $W_{3}$. The cost per grid point of SUCCSD(SD) is $\mathcal{O}(N^{6})$, and that of SUCCSD(SDT) is $\mathcal{O}(N^{7})$. A third approach we consider is to solve for an approximate $W_{3}$ from the algebraic equations, by setting $W_{3}$ through $W_{5}$ to zero in the right-hand-side of Eqn. 58, thereby obtaining an approximation for $W_{3}$ in terms of $W_{1}$ and $W_{2}$. We call the resulting method SUCCSD(SDt), and it scales as $\mathcal{O}(N^{6})$. Because including linear response in the energy expression or solving the amplitude equations would require through $W_{4}$, in this work we use the ODE technique only to evaluate the energy in a PAV approach without response. This allows us to focus on assessing the error introduced by our approximations for $W_{k+1}$.

IV.3 Results

The key objects of practical implementations of our spin projected coupled cluster theory are the disentangled cluster operators $W_{k}$. All our approximations are based on assuming that $W_{k}$ for sufficiently high excitation level are negligible or at least readily approximated by the $W_{k}$ for lower excitation levels. It is thus critical to check that this is indeed the case. Figure 4 shows the behavior of $W_{k}$ for the ten-site half-filled periodic Hubbard lattice at several values of interaction strength $U/t$. As we have discussed, only $W_{2}$ is non-zero at $\beta=0$. As $\beta$ is increased, other amplitudes appear. In general, $W_{k}$ indeed decays as excitation level increases for all $\beta$ and for all interaction strengths. This justifies our various truncation schemes. Note, however, that $W_{0}$ actually increases as $\beta$ increases. We will discuss this increase later.

A different way of justifying our truncation schemes is by comparison with the untruncated results. This is done in Fig. 5. Unsurprisingly, SUCCSD(SDT) is closer to the exact SUCCSD than is SUCCSD(SD). More interestingly, SUCCSD(SDt) is closer to the exact result than is SUCCSD(SDT). While this may seem counterintuitive, it ultimately arises from the observation that the algebraic equations lead to a more accurate approximation for $W_{3}$ than does solving the differential equation assuming $W_{4}=0$. We should note that both SUCCSD(SD) and SUCCSD(SDt) scale the same as UCCSD, and both are significantly more accurate.

To bear out our explanation for the superiority of SUCCSD(SDt) over SUCCSD(SDT), Fig. 6 shows that SUCCSD(SDt) produces more accurate disentangled amplitudes (as a consequence of which it also produces more accurate kernels, as shown in the Supplementary Material). Figure 6 also shows the behavior of the order-by-order expansion of the disentangled amplitudes. At small $\beta$ where $V_{1}$ is small, this expansion works exceptionally well, but as $\beta$ increases its quality deteriorates rapidly.

One might be concerned that approximating the disentangled cluster operators spoils the symmetry projection. After all, approximating $W$ implies approximately handling $V_{1}$, hence approximately handly the rotation operator $R$ and hence approximating the symmetry projection. Figure 7 shows that, at least with our SUCCSD(SDt) approximation, this is not a concern. Here, we plot the error in the expectation value of $\langle S^{2}\rangle$ per electron in the periodic half-filled eighteen-site Hubbard lattice. We use the reference determinant of SUHF as our mean-field reference so the mean-field breaks symmetry everywhere.

To be clear about what we plot, the correct definition of the coupled cluster expectation value of an operator $\mathcal{O}$ is

[TABLE]

This is what we have plotted as “LR-UCCSD” where $\mathcal{O}=S^{2}$. Neglecting $Z$ gives us the curve we have plotted as “UCCSD.” For projected coupled cluster, the expectation value would similarly be

[TABLE]

Evaluating this projected coupled cluster expectation value in the presence of $Z$ would require excited kernels which we have not yet constructed, so we have plotted it without $Z$, as “SUCCSD(SDt).” If the projection is done exactly, the expectation value of $S^{2}$ would be exact even in the absence of $Z$.

While the mean-field is badly symmetry broken for large Hubbard parameter $U/t$, unrestricted CCSD greatly reduces the degree of symmetry breaking. Interestingly, including linear response has minimal effect in this case. Our SUCCSD(SDt) gives almost exact symmetry projection, with error per electron in the expectation value of $S^{2}$ on the order of $10^{-3}$ (and very slightly negative in places, which is possible because we do not use a Hermitian expectation value).

IV.4 The Thermodynamic Limit

It is known that spin projected UHF is not size-extensive, or, more precisely, that its size-extensive energetic component is the same as that of the broken symmetry mean-field.Jiménez-Hoyos et al. (2012) Here, we demonstrate numerically that SUCCSD likewise has no size extensive correction beyond UCCSD, but that the SUCCSD energy returns to the UCCSD value at a much slower rate than the SUHF energy returns to UHF. This bodes well for large finite systems, and we expect SUCCSD to provide useful improvements to UCCSD even when SUHF provides only marginal improvements upon UHF. We should note that our results here are generated for the SUCCSD(SDt) approximation to the exact SUCCSD, and in the PAV sense, but we expect the exact VAP SUCCSD to display broadly similar features although we have not proven as much.

From the differential equation one can see that the derivatives of the $W_{k}$ contain only connected terms. At $\beta=0$ where $V_{1}=0$, $W=U_{2}$ and is extensive. This suggests that $W$ at every $\beta$ should be extensive and $W$ should converge with respect to system size at a rate similar to UCCSD. Since the reduced Hamiltonian kernel contains only connected terms involving $W$, it, too, should be size extensive. This is indeed the case, as can be seen from Fig. 8. The Hamiltonian kernel per electron is already converged with respect to system size for the 22-site Hubbard ring for all interaction strengths, which is where UCCSD also converges. The SUHF Hamiltonian kernel also converges rapidly with respect to system size.

The problem, then, must be the norm kernel, which as we see converges very slowly. In fact, in the thermodynamic limit the norm kernel for SUHF becomes a delta function at $\beta=0$; thus the energy only samples $\beta=0$, and we therefore obtain the UHF result.

In SUCC theory, the norm kernel is $\langle\phi|R(\beta)|\phi\rangle\,\mathrm{e}^{W_{0}}$, where $\langle\phi|R(\beta)|\phi\rangle\equiv\mathrm{e}^{S_{0}}$ is the SUHF norm kernel. Figure 9 shows how $W_{0}$ and $S_{0}$ behave as a function of system size. It is clear that

[TABLE]

where $w_{0}$ and $s_{0}$ are constants, so that the coupled cluster norm kernel is

[TABLE]

However, $s_{0}$ and $w_{0}$ have opposite signs; as the number of electrons increases, the SUHF norm kernel goes to zero while $W_{0}$ goes to infinity. This makes it difficult to tell how the SUCC norm kernel behaves, although clearly it decays less rapidly than does the SUHF norm kernel. To help resolve this behavior, Fig. 10 shows that $s_{0}$ is negative for all $\beta$, while $w_{0}$ is positive and thus counteracts the decay of the mean-field kernel, but their sum is always negative except at $\beta=0$. This means that the norm kernel in the thermodynamic limit vanishes except at $\beta=0$, just as it does for SUHF. Thus, SUCCSD approaches UCCSD in the thermodynamic limit, but at a much slower rate than SUHF approaches UHF (because $s_{0}+w_{0}$ is less negative than $s_{0}$ alone). Figure 11 shows that SUCCSD adds considerable correlation atop of UCCSD even at large system sizes (see results for the 46-site lattice).

We close with a brief discussion of how the kernels should behave as the theory becomes even more complete. In the full coupled cluster limit, the wave function, though written in the language of symmetry broken Hartree-Fock, has good quantum numbers (for finite systems). Accordingly, the norm kernel is one for all $\Omega$ and the Hamiltonian kernel is the exact ground state energy. Our results here show that SUCCSD has flatter norm and Hamiltonian kernels than does SUHF, and as we continue to increase the level of correlation in the cluster operator, we would expect to see the kernels become even flatter.

V Discussion

V.1 Comparison With Duguet’s Symmetry Broken and Restored Coupled Cluster Theory

Recently, Duguet proposed a symmetry broken and restored coupled cluster theory.Duguet (2015) Here, we briefly discuss the differences between our two approaches.

The biggest difference is that we begin with an explicit wave function ansatz rather than with the kernels. As a result, the linear response and optimization of the cluster operator $U$ force us to introduce excited Hamiltonian and norm kernels which do not appear in Duguet’s theory. Second, Duguet’s approach includes solving a coupled cluster-like equation at each $\Omega$ where we solve for the cluster operator $U$ only once and use these amplitudes to generate the $\Omega$-dependent kernels. Finally, we introduce the disentangled cluster operators $W$, which conveniently encapsulate all information from the various kernels we require. We may solve a set of coupled differential equations for the disentangled cluster operators. While Duguet also proposes a differential equation, it is used only for the norm kernel.

V.2 Conclusions

This work combines coupled cluster and symmetry projection into one tool which seems to have considerable promise. The broken symmetry coupled cluster is already fairly accurate for many cases, and the symmetry restored coupled cluster is even better while retaining good symmetries. Even the simple projection after variation scheme is a useful improvement on broken symmetry coupled cluster, and results can be improved further by including linear response or by reoptimizing the broken symmetry cluster operator in the presence of the symmetry projector. For the small half-filled Hubbard rings for which we have been able to carry out the exact VAP SUCCSD calculations, errors are well under $0.001t$ per electron.

Symmetry projected coupled cluster reduces energetically to broken symmetry coupled cluster in the thermodynamic limit, just as projected Hartree-Fock reduces to the broken symmetry mean field. However, the symmetry projected coupled cluster does so more slowly, and useful improvements upon broken symmetry coupled cluster are possible even for fairly sizable systems.

An important quantity introduced in this work is the disentangled cluster operators $W$. These allow us to work in the same orthogonal excitation framework as traditional coupled cluster theory. The needed kernels can be evaluated from disentangled cluster operators of low excitation rank. The disentangled cluster operators decay as the excitation level increases and, for unrestricted spin projection, satisfy a set of simple differential equations which couple $W_{k}$ of different excitation levels. We have shown several truncation schemes for the approximate construction of the disentangled cluster operators.

Much work, of course, remains. We have by no means explored all possible ways of approximating the disentangled cluster operators, nor have we used them in the construction of excited kernels, which one must do if one wishes to make linear response and VAP calculations practical. Neither have we implemented the combination of spin projected coupled cluster based on a generalized Hartree-Fock reference which breaks both $S^{2}$ and $S_{z}$ symmetries, which we expect it to be significantly more accurate yet.Henderson and Scuseria Thus far we have considered only the ground state energy, and naturally properties, gradients, and excited states can all in principle be accessed by suitable modifications of traditional coupled cluster methods. Nonetheless, we are greatly encouraged by this early foray into the symmetry projection of broken symmetry coupled cluster theory.

Supplementary Material

See supplementary material for data corresponding to Figs. 6, 8, 9, and 10 at other values of Hubbard parameters $U/t$. We also show errors in the disentangled cluster operator $W_{2}$ and the norm and Hamiltonian kernels corresponding to different approximations to SUCCSD.

Acknowledgements.

This work was supported by the National Science Foundation under award CHE-1462434. G.E.S. is a Welch Foundation Chair (C-0036). We would like to acknowledge computational support provided by the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award DE-SC0012575.

Appendix A Non-Zero Components of Auxilliary Kernels

In Eqn. 43 we introduced the auxiliary kernels $\tilde{\mathcal{N}}$ and $\tilde{\mathcal{H}}$. Here we list their nonzero elements for convenience. They are

[TABLE]

Appendix B Including $U_{1}$ in the Disentangled Cluster Formalism

In our disentangled cluster formalism we incorporate $U_{1}$ by using it to transform the Hamiltonian and symmetry projectors rather than explicitly including it in our equations. Here, we show how we do this. We will use the singly-excited Hamiltonian kernel $\mathcal{H}_{i}^{a}$ to demonstrate the general idea.

The excited Hamiltonian kernel we need is

[TABLE]

where $\bar{R}_{U_{1}}$ and $\bar{H}_{U_{1}}$ are

[TABLE]

The orbital transformation matrix $\bm{\tau}$ is

[TABLE]

where $\mathbf{U}_{\mathrm{vo}}$ is the matrix of $U_{1}$ amplitudes. The practical result is that the integrals defining the Hamiltonian and rotation operator should be replaced by the transformed ones and, this having been done, the transformation matrix $\mathbf{T}$ of Eqn. 41 should be replaced by $\mathbf{T}\,\bm{\tau}$.

Appendix C Proof of the Disentanglement Equation

Here we prove our formula for the disentangled cluster operators given in Eqn. 45a.

We begin by writing the wave function $|\psi\rangle=\mathrm{e}^{V_{1}}\,\mathrm{e}^{U_{2}}|\phi\rangle$ in two different ways:

[TABLE]

The excitation operators $C_{k}$ are given by the linked excitation components of

[TABLE]

where by “excitation components” we mean the maximally contracted part (so that all indices on the $V_{1}$ are contracted with indices on $U_{2}$). The cluster operator $W_{k}$ is the fully connected part of $C_{k}$, which proves Eqn. 45a for $k>0$. The case $k=0$ needed for $W_{0}$ is a special case of the linked diagram theorem.

Appendix D Proof of the Derivative of $V_{1}$

In Eqn. 59 we write the derivative of the operator $V_{1}$ with respect to $\beta$. Here we prove that result.

Recall that the matrix $\mathbf{V}_{\mathrm{ov}}$ contains the amplitudes defining $V_{1}$. From Thouless’ theorem,

[TABLE]

where $\mathbf{R}_{\mathrm{oo}}$ and $\mathbf{R}_{\mathrm{ov}}$ are respectively the occupied-occupied and occupied-virtual blocks of $\mathbf{R}$ (and similarly, subscripts “vv” and “vo” would denote the virtual-virtual and virtual-occupied blocks). This means that

[TABLE]

The rotation matrix $\mathbf{R}$ is

[TABLE]

where recall that $\mathbf{S}_{y}$ is the matrix representation of the operator $S_{y}$. Its derivative with respect to $\beta$ is therefore

[TABLE]

where for convenience we have defined

[TABLE]

Inserting the derivative of $\mathbf{R}$ into Eqn. 80 gives us

[TABLE]

which is the occupied-virtual block of the operator $\mathrm{e}^{V_{1}}\,A\,\mathrm{e}^{-V_{1}}$. Accordingly,

[TABLE]

where the subscript “d” means we keep the deexcitation part only and where we have used that $A=-\mathrm{i}\,S_{y}$.