$$
\mid\Psi\rangle = \exp\left(-\hat{T}\right)\mid\Phi\rangle\ .
$$
We will discuss how to define the operators later in this work. This simple
ansatz carries enormous power. It leads to a non-perturbative many-body
theory that includes summation of ladder diagrams , ring
diagrams, and an infinite-order
generalization of many-body perturbation theory.
Before we proceed with the derivation of the Coupled cluster equations, let us repeat some of the arguments we presented during our FCI lectures. In our FCI discussions, we rewrote the solution of the Schroedinger equation as a set of coupled equationsin the unknown coefficients \( C \). Let us repeat some of these arguments. To obtain the eigenstates and eigenvalues in terms of non-linear equations is not a very practical approach. However, it serves the scope of linking FCI theory with approximative solutions to the many-body problem like Coupled cluster (CC) theory
$$
\langle \Phi_0 | \hat{H} -E| \Phi_0\rangle + \sum_{ai}\langle \Phi_0 | \hat{H} -E|\Phi_{i}^{a} \rangle C_{i}^{a}+
\sum_{abij}\langle \Phi_0 | \hat{H} -E|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}=0,
$$
or
$$
E-E_0 =\Delta E=\sum_{ai}\langle \Phi_0 | \hat{H}|\Phi_{i}^{a} \rangle C_{i}^{a}+
\sum_{abij}\langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab} \rangle C_{ij}^{ab},
$$
where the energy \( E_0 \) is the reference energy and \( \Delta E \) defines the so-called correlation energy.
The single-particle basis functions could be the results of a Hartree-Fock calculation or just the eigenstates of the non-interacting part of the Hamiltonian.
$$
E-E_0 =\Delta E^{HF}=\sum_{abij}\langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}.
$$
$$
\Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+
\sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}.
$$
This equation determines the correlation energy but not the coefficients \( C \).
We need more equations. Our next step is to set up
$$
\langle \Phi_i^a | \hat{H} -E| \Phi_0\rangle + \sum_{bj}\langle \Phi_i^a | \hat{H} -E|\Phi_{j}^{b} \rangle C_{j}^{b}+
\sum_{bcjk}\langle \Phi_i^a | \hat{H} -E|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+
\sum_{bcdjkl}\langle \Phi_i^a | \hat{H} -E|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}=0,
$$
as this equation will allow us to find an expression for the coefficents \( C_i^a \) since we can rewrite this equation as
$$
\langle i | \hat{f}| a\rangle +\langle \Phi_i^a | \hat{H}|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{bj\ne ai}\langle \Phi_i^a | \hat{H}|\Phi_{j}^{b} \rangle C_{j}^{b}+
\sum_{bcjk}\langle \Phi_i^a | \hat{H}|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+
\sum_{bcdjkl}\langle \Phi_i^a | \hat{H}|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}=EC_i^a.
$$
$$
C_{i}^{a}=\frac{\langle i | \hat{f}| a\rangle}{\epsilon_i-\epsilon_a}.
$$
The observant reader will however see that we need an equation for \( C_{jk}^{bc} \) and \( C_{jkl}^{bcd} \) as well. To find equations for these coefficients we need then to continue our multiplications from the left with the various \( \Phi_{H}^P \) terms.
$$
\langle \Phi_{ij}^{ab} | \hat{H} -E| \Phi_0\rangle + \sum_{kc}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{k}^{c} \rangle C_{k}^{c}+
$$
$$
\sum_{cdkl}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{kl}^{cd} \rangle C_{kl}^{cd}+\sum_{cdeklm}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{klm}^{cde} \rangle C_{klm}^{cde}+\sum_{cdefklmn}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{klmn}^{cdef} \rangle C_{klmn}^{cdef}=0,
$$
and we can isolate the coefficients \( C_{kl}^{cd} \) in a similar way as we did for the coefficients \( C_{i}^{a} \).
A standard choice for the first iteration is to set
$$
C_{ij}^{ab} =\frac{\langle ij \vert \hat{v} \vert ab \rangle}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}.
$$
$$
\Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+
\sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}.
$$
The coefficients \( C \) result from the solution of the eigenvalue problem. The energy of say the ground state is then
$$
E=E_{ref}+\Delta E,
$$
where the so-called reference energy is the energy we obtain from a Hartree-Fock calculation, that is
$$
E_{ref}=\langle \Phi_0 \vert \hat{H} \vert \Phi_0 \rangle.
$$
Popular methods are
All these methods start normally with a Hartree-Fock basis as the calculational basis.
The ansatz for the wavefunction (ground state) is given by
$$
\begin{equation*}
\vert \Psi\rangle = \vert \Psi_{CC}\rangle = e^{\hat{T}} \vert \Phi_0\rangle =
\left( \sum_{n=1}^{A} \frac{1}{n!} \hat{T}^n \right) \vert \Phi_0\rangle,
\end{equation*}
$$
where \( A \) represents the maximum number of particle-hole excitations and \( \hat{T} \) is the cluster operator defined as
$$
\begin{align*}
\hat{T} &= \hat{T}_1 + \hat{T}_2 + \ldots + \hat{T}_A \\
\hat{T}_n &= \left(\frac{1}{n!}\right)^2
\sum_{\substack{
i_1,i_2,\ldots i_n \\
a_1,a_2,\ldots a_n}}
t_{i_1i_2\ldots i_n}^{a_1a_2\ldots a_n} a_{a_1}^\dagger a_{a_2}^\dagger \ldots a_{a_n}^\dagger a_{i_n} \ldots a_{i_2} a_{i_1}.
\end{align*}
$$
$$
\begin{equation*}
E_{\mathrm{CC}} = \langle\Phi_0\vert \overline{H}\vert \Phi_0\rangle,
\end{equation*}
$$
where \( \overline{H} \) is a similarity transformed Hamiltonian
$$
\begin{align*}
\overline{H}&= e^{-\hat{T}} \hat{H}_N e^{\hat{T}} \\
\hat{H}_N &= \hat{H} - \langle\Phi_0\vert \hat{H} \vert \Phi_0\rangle.
\end{align*}
$$
$$
\begin{equation*}
0 = \langle\Phi_{i_1 \ldots i_n}^{a_1 \ldots a_n}\vert \overline{H}\vert \Phi_0\rangle.
\end{equation*}
$$
The similarity transformed Hamiltonian \( \overline{H} \) is expanded using the Baker-Campbell-Hausdorff expression,
$$
\begin{align*}
\overline{H}&= \hat{H}_N + \left[ \hat{H}_N, \hat{T} \right] +
\frac{1}{2} \left[\left[ \hat{H}_N, \hat{T} \right], \hat{T}\right] + \ldots \\
& \quad \frac{1}{n!} \left[ \ldots \left[ \hat{H}_N, \hat{T} \right], \ldots \hat{T} \right] +\dots
\end{align*}
$$
and simplified using the connected cluster theorem
$$
\begin{equation*}
\overline{H}= \hat{H}_N + \left( \hat{H}_N \hat{T}\right)_c + \frac{1}{2} \left( \hat{H}_N \hat{T}^2\right)_c
+ \dots + \frac{1}{n!} \left( \hat{H}_N \hat{T}^n\right)_c +\dots
\end{equation*}
$$
The coupled cluster wavefunction is now given by
$$
\begin{equation*}
\vert \Psi_{CC}\rangle = e^{\hat{T}_1 + \hat{T}_2} \vert \Phi_0\rangle
\end{equation*}
$$
where
$$
\begin{align*}
\hat{T}_1 &=
\sum_{ia}
t_{i}^{a} a_{a}^\dagger a_i \\
\hat{T}_2 &= \frac{1}{4}
\sum_{ijab}
t_{ij}^{ab} a_{a}^\dagger a_{b}^\dagger a_{j} a_{i}.
\end{align*}
$$
If we truncate our equations at the CCSD level, it corresponds to performing a transformation of the Hamiltonian matrix of the following type for a six particle problem (with a two-body Hamiltonian):
| \( 0p-0h \) | \( 1p-1h \) | \( 2p-2h \) | \( 3p-3h \) | \( 4p-4h \) | \( 5p-5h \) | \( 6p-6h \) | |
| \( 0p-0h \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | 0 | 0 | 0 | 0 |
| \( 1p-1h \) | 0 | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | 0 | 0 | 0 |
| \( 2p-2h \) | 0 | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | 0 | 0 |
| \( 3p-3h \) | 0 | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | 0 |
| \( 4p-4h \) | 0 | 0 | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) |
| \( 5p-5h \) | 0 | 0 | 0 | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) |
| \( 6p-6h \) | 0 | 0 | 0 | 0 | \( \tilde{x} \) | \( \tilde{x} \) | \( \tilde{x} \) |
$$
\Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+
\sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}.
$$
In Coupled cluster theory it becomes (irrespective of level of truncation of \( T \))
$$
\Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle t_{i}^{a}+
\sum_{abij}\langle ij | \hat{v}| ab \rangle t_{ij}^{ab}.
$$
There are several interesting features:
However, Coupled Cluster theory is
We will now approximate the cluster operator \( \hat{T} \) to include only \( 2p-2h \) correlations. This leads to the so-called CCD approximation, that is
$$
\hat{T}\approx \hat{T}_2=\frac{1}{4}\sum_{abij}t_{ij}^{ab}a^{\dagger}_aa^{\dagger}_ba_ja_i,
$$
meaning that we have
$$
\vert \Psi_0 \rangle \approx \vert \Psi_{CCD} \rangle = \exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle.
$$
$$
\hat{H}=\hat{H}_N+E_{\mathrm{ref}},
$$
with
$$
\hat{H}_N=\sum_{pq}\langle p \vert \hat{f} \vert q \rangle a^{\dagger}_pa_q + \frac{1}{4}\sum_{pqrs}\langle pq \vert \hat{v} \vert rs \rangle a^{\dagger}_pa^{\dagger}_qa_sa_r,
$$
we obtain that the energy can be written as
$$
\langle \Phi_0 \vert \exp{-\left(\hat{T}_2\right)}\hat{H}_N\exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle =
\langle \Phi_0 \vert \hat{H}_N(1+\hat{T}_2)\vert \Phi_0\rangle = E_{CCD}.
$$
$$
E_{CCD}=E_{\mathrm{ref}}+\frac{1}{4}\sum_{abij}\langle ij \vert \hat{v} \vert ab \rangle t_{ij}^{ab},
$$
where the latter is the correlation energy from this level of approximation of CC theory.
Similarly, the expression for the amplitudes reads
$$
\langle \Phi_{ij}^{ab} \vert \exp{-\left(\hat{T}_2\right)}\hat{H}_N\exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle = 0.
$$
$$
\begin{align}
0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab} & \nonumber \\
+\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}+\hat{P}(ij\vert ab)\sum_{kc} \langle kb \vert \hat{v} \vert cj \rangle t_{ik}^{ac} & \nonumber \\
+\frac{1}{4}\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ij}^{cd}t_{kl}^{ab}+\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{ac}t_{jl}^{bd}& \nonumber \\
-\frac{1}{2}\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{dc}t_{lj}^{ab}-\frac{1}{2}\hat{P}(ab)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{lk}^{ac}t_{ij}^{db},&
\tag{1}
\end{align}
$$
where we have defined
$$
\hat{P}\left(ab\right)= 1-\hat{P}_{ab},
$$
where \( \hat{P}_{ab} \) interchanges two particles occupying the quantum numbers \( a \) and \( b \).
$$
\hat{P}(ij\vert ab) = (1-\hat{P}_{ij})(1-\hat{P}_{ab}).
$$
Recall also that the unknown amplitudes \( t_{ij}^{ab} \)
represent anti-symmetrized matrix elements, meaning that they obey the same symmetry relations as the two-body interaction, that is
$$
t_{ij}^{ab}=-t_{ji}^{ab}=-t_{ij}^{ba}=t_{ji}^{ba}.
$$
The two-body matrix elements are also anti-symmetrized, meaning that
$$
\langle ab \vert \hat{v} \vert ij \rangle = -\langle ab \vert \hat{v} \vert ji \rangle= -\langle ba \vert \hat{v} \vert ij \rangle=\langle ba \vert \hat{v} \vert ji \rangle.
$$
The non-linear equations for the unknown amplitudes \( t_{ij}^{ab} \) are solved iteratively. We discuss the implementation of these equations below.
It is useful to make approximations to the equations for the amplitudes. The standard method for solving these equations is to set up an iterative scheme where method's like Newton's method or similar root searching methods are used to find the amplitudes. Itreative solvers need a guess for the amplitudes. A good starting point is to use the correlated wave operator from perturbation theory to first order in the interaction. This means that we define the zeroth approximation to the amplitudes as
$$
t^{(0)}=\frac{\langle ab \vert \hat{v} \vert ij \rangle}{\left(\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b\right)},
$$
leading to our first approximation for the correlation energy at the CCD level to be equal to second-order perturbation theory without \( 1p-1h \) excitations, namely
$$
\Delta E_{\mathrm{CCD}}^{(0)}=\frac{1}{4}\sum_{abij} \frac{\langle ij \vert \hat{v} \vert ab \rangle \langle ab \vert \hat{v} \vert ij \rangle}{\left(\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b\right)}.
$$
$$
\begin{equation}
0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab}+\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}.
\tag{2}
\end{equation}
$$
$$
\tau_{ij}^{ab}= \left(\omega-\epsilon_a-\epsilon_b\right)t_{ij}^{ab},
$$
and inserting
$$
1=\frac{\left(\omega-\epsilon_c-\epsilon_d\right)}{\left(\omega-\epsilon_c-\epsilon_d\right)},
$$
in the intermediate sums over \( cd \) in Eq. (2), we can rewrite the latter equation as
$$
\tau_{ij}^{ab}(\omega)= \langle ab \vert \hat{v} \vert ij \rangle + \frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle \frac{1}{\omega-\epsilon_c-\epsilon_d}\tau_{ij}^{cd}(\omega),
$$
where we have indicated an explicit energy dependence. This equation, transforming a two-particle configuration into a single index, can be transformed into a matrix inversion problem. Solving the equations for a fixed energy \( \omega \) allows us to compare directly with results from Green's function theory when only two-particle intermediate states are included.
$$
\begin{equation}
0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab}+\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}.
\tag{3}
\end{equation}
$$
This equation is solved the same way as we would do for Eq. (2). The final step is then to include all terms in Eq. (1).