Many-body perturbation theory

Morten Hjorth-Jensen, National Superconducting Cyclotron Laboratory and Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA & Department of Physics, University of Oslo, Oslo, Norway

July 6-24 2015

Table of contents

Many-body perturbation theory
Interpreting the correlation energy and the wave operator

Many-body perturbation theory

We assume here that we are only interested in the ground state of the system and expand the exact wave function in term of a series of Slater determinants $$ \vert \Psi_0\rangle = \vert \Phi_0\rangle + \sum_{m=1}^{\infty}C_m\vert \Phi_m\rangle, $$ where we have assumed that the true ground state is dominated by the solution of the unperturbed problem, that is $$ \hat{H}_0\vert \Phi_0\rangle= W_0\vert \Phi_0\rangle. $$ The state \( \vert \Psi_0\rangle \) is not normalized, rather we have used an intermediate normalization \( \langle \Phi_0 \vert \Psi_0\rangle=1 \) since we have \( \langle \Phi_0\vert \Phi_0\rangle=1 \).

The Schroedinger equation is $$ \hat{H}\vert \Psi_0\rangle = E\vert \Psi_0\rangle, $$ and multiplying the latter from the left with \( \langle \Phi_0\vert \) gives $$ \langle \Phi_0\vert \hat{H}\vert \Psi_0\rangle = E\langle \Phi_0\vert \Psi_0\rangle=E, $$ and subtracting from this equation $$ \langle \Psi_0\vert \hat{H}_0\vert \Phi_0\rangle= W_0\langle \Psi_0\vert \Phi_0\rangle=W_0, $$ and using the fact that the both operators \( \hat{H} \) and \( \hat{H}_0 \) are hermitian results in $$ \Delta E=E-W_0=\langle \Phi_0\vert \hat{H}_I\vert \Psi_0\rangle, $$ which is an exact result. We call this quantity the correlation energy.

This equation forms the starting point for all perturbative derivations. However, as it stands it represents nothing but a mere formal rewriting of Schroedinger's equation and is not of much practical use. The exact wave function \( \vert \Psi_0\rangle \) is unknown. In order to obtain a perturbative expansion, we need to expand the exact wave function in terms of the interaction \( \hat{H}_I \).

Here we have assumed that our model space defined by the operator \( \hat{P} \) is one-dimensional, meaning that $$ \hat{P}= \vert \Phi_0\rangle \langle \Phi_0\vert , $$ and $$ \hat{Q}=\sum_{m=1}^{\infty}\vert \Phi_m\rangle \langle \Phi_m\vert . $$

We can thus rewrite the exact wave function as $$ \vert \Psi_0\rangle= (\hat{P}+\hat{Q})\vert \Psi_0\rangle=\vert \Phi_0\rangle+\hat{Q}\vert \Psi_0\rangle. $$ Going back to the Schr\"odinger equation, we can rewrite it as, adding and a subtracting a term \( \omega \vert \Psi_0\rangle \) as $$ \left(\omega-\hat{H}_0\right)\vert \Psi_0\rangle=\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle, $$ where \( \omega \) is an energy variable to be specified later.

We assume also that the resolvent of \( \left(\omega-\hat{H}_0\right) \) exits, that is it has an inverse which defined the unperturbed Green's function as $$ \left(\omega-\hat{H}_0\right)^{-1}=\frac{1}{\left(\omega-\hat{H}_0\right)}. $$

We can rewrite Schroedinger's equation as $$ \vert \Psi_0\rangle=\frac{1}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle, $$ and multiplying from the left with \( \hat{Q} \) results in $$ \hat{Q}\vert \Psi_0\rangle=\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle, $$ which is possible since we have defined the operator \( \hat{Q} \) in terms of the eigenfunctions of \( \hat{H} \).

These operators commute meaning that $$ \hat{Q}\frac{1}{\left(\omega-\hat{H}_0\right)}\hat{Q}=\hat{Q}\frac{1}{\left(\omega-\hat{H}_0\right)}=\frac{\hat{Q}}{\left(\omega-\hat{H}_0\right)}. $$ With these definitions we can in turn define the wave function as $$ \vert \Psi_0\rangle=\vert \Phi_0\rangle+\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle. $$ This equation is again nothing but a formal rewrite of Schr\"odinger's equation and does not represent a practical calculational scheme. It is a non-linear equation in two unknown quantities, the energy \( E \) and the exact wave function \( \vert \Psi_0\rangle \). We can however start with a guess for \( \vert \Psi_0\rangle \) on the right hand side of the last equation.

The most common choice is to start with the function which is expected to exhibit the largest overlap with the wave function we are searching after, namely \( \vert \Phi_0\rangle \). This can again be inserted in the solution for \( \vert \Psi_0\rangle \) in an iterative fashion and if we continue along these lines we end up with $$ \vert \Psi_0\rangle=\sum_{i=0}^{\infty}\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle, $$ for the wave function and $$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle, $$ which is now a perturbative expansion of the exact energy in terms of the interaction \( \hat{H}_I \) and the unperturbed wave function \( \vert \Psi_0\rangle \).

In our equations for \( \vert \Psi_0\rangle \) and \( \Delta E \) in terms of the unperturbed solutions \( \vert \Phi_i\rangle \) we have still an undetermined parameter \( \omega \) and a dependecy on the exact energy \( E \). Not much has been gained thus from a practical computational point of view.

In Brilluoin-Wigner perturbation theory it is customary to set \( \omega=E \). This results in the following perturbative expansion for the energy \( \Delta E \) $$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle= $$ $$ \langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I+ \hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I+\dots\right)\vert \Phi_0\rangle. $$ $$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle=$$ $$ \langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I+ \hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I+\dots\right)\vert \Phi_0\rangle. $$ This expression depends however on the exact energy \( E \) and is again not very convenient from a practical point of view. It can obviously be solved iteratively, by starting with a guess for \( E \) and then solve till some kind of self-consistency criterion has been reached.

Actually, the above expression is nothing but a rewrite again of the full Schr\"odinger equation.

Defining \( e=E-\hat{H}_0 \) and recalling that \( \hat{H}_0 \) commutes with \( \hat{Q} \) by construction and that \( \hat{Q} \) is an idempotent operator \( \hat{Q}^2=\hat{Q} \). Using this equation in the above expansion for \( \Delta E \) we can write the denominator $$ \hat{Q}\frac{1}{\hat{e}-\hat{Q}\hat{H}_I\hat{Q}}= $$ $$ \hat{Q}\left[\frac{1}{\hat{e}}+\frac{1}{\hat{e}}\hat{Q}\hat{H}_I\hat{Q} \frac{1}{\hat{e}}+\frac{1}{\hat{e}}\hat{Q}\hat{H}_I\hat{Q} \frac{1}{\hat{e}}\hat{Q}\hat{H}_I\hat{Q}\frac{1}{\hat{e}}+\dots\right]\hat{Q}. $$

Inserted in the expression for \( \Delta E \) leads to $$ \Delta E= \langle \Phi_0\vert \hat{H}_I+\hat{H}_I\hat{Q}\frac{1}{E-\hat{H}_0-\hat{Q}\hat{H}_I\hat{Q}}\hat{Q}\hat{H}_I\vert \Phi_0\rangle. $$ In RS perturbation theory we set \( \omega = W_0 \) and obtain the following expression for the energy difference $$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{W_0-\hat{H}_0}\left(\hat{H}_I-\Delta E\right)\right\}^i\vert \Phi_0\rangle= $$ $$ \langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)+ \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)+\dots\right)\vert \Phi_0\rangle. $$

Recalling that \( \hat{Q} \) commutes with \( \hat{H_0} \) and since \( \Delta E \) is a constant we obtain that $$ \hat{Q}\Delta E\vert \Phi_0\rangle = \hat{Q}\Delta E\vert \hat{Q}\Phi_0\rangle = 0. $$ Inserting this results in the expression for the energy results in $$ \Delta E=\langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I+ \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I+\dots\right)\vert \Phi_0\rangle. $$

We can now this expression in terms of a perturbative expression in terms of \( \hat{H}_I \) where we iterate the last expression in terms of \( \Delta E \) $$ \Delta E=\sum_{i=1}^{\infty}\Delta E^{(i)}. $$ We get the following expression for \( \Delta E^{(i)} \) $$ \Delta E^{(1)}=\langle \Phi_0\vert \hat{H}_I\vert \Phi_0\rangle, $$ which is just the contribution to first order in perturbation theory, $$ \Delta E^{(2)}=\langle\Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\vert \Phi_0\rangle, $$ which is the contribution to second order. $$ \Delta E^{(3)}=\langle \Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\Phi_0\rangle- \langle\Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\langle \Phi_0\vert \hat{H}_I\vert \Phi_0\rangle\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\vert \Phi_0\rangle, $$ being the third-order contribution.

Interpreting the correlation energy and the wave operator

In the shell-model lectures we showed that we could rewrite the exact state function for say the ground state, as a linear expansion in terms of all possible Slater determinants. That is, we define the ansatz for the ground state as $$ |\Phi_0\rangle = \left(\prod_{i\le F}\hat{a}_{i}^{\dagger}\right)|0\rangle, $$ where the index \( i \) defines different single-particle states up to the Fermi level. We have assumed that we have \( N \) fermions. A given one-particle-one-hole (\( 1p1h \)) state can be written as $$ |\Phi_i^a\rangle = \hat{a}_{a}^{\dagger}\hat{a}_i|\Phi_0\rangle, $$ while a \( 2p2h \) state can be written as $$ |\Phi_{ij}^{ab}\rangle = \hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_j\hat{a}_i|\Phi_0\rangle, $$ and a general \( ApAh \) state as $$ |\Phi_{ijk\dots}^{abc\dots}\rangle = \hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{c}^{\dagger}\dots\hat{a}_k\hat{a}_j\hat{a}_i|\Phi_0\rangle. $$

We use letters \( ijkl\dots \) for states below the Fermi level and \( abcd\dots \) for states above the Fermi level. A general single-particle state is given by letters \( pqrs\dots \).

We can then expand our exact state function for the ground state as $$ |\Psi_0\rangle=C_0|\Phi_0\rangle+\sum_{ai}C_i^a|\Phi_i^a\rangle+\sum_{abij}C_{ij}^{ab}|\Phi_{ij}^{ab}\rangle+\dots =(C_0+\hat{C})|\Phi_0\rangle, $$ where we have introduced the so-called correlation operator $$ \hat{C}=\sum_{ai}C_i^a\hat{a}_{a}^{\dagger}\hat{a}_i +\sum_{abij}C_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_j\hat{a}_i+\dots $$ Since the normalization of \( \Psi_0 \) is at our disposal and since \( C_0 \) is by hypothesis non-zero, we may arbitrarily set \( C_0=1 \) with corresponding proportional changes in all other coefficients. Using this so-called intermediate normalization we have $$ \langle \Psi_0 | \Phi_0 \rangle = \langle \Phi_0 | \Phi_0 \rangle = 1, $$ resulting in $$ |\Psi_0\rangle=(1+\hat{C})|\Phi_0\rangle. $$

In a shell-model calculation, the unknown coefficients in \( \hat{C} \) are the eigenvectors which result from the diagonalization of the Hamiltonian matrix.

How can we use perturbation theory to determine the same coefficients? Let us study the contributions to second order in the interaction, namely $$ \Delta E^{(2)}=\langle\Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\vert \Phi_0\rangle. $$

The intermediate states given by \( \hat{Q} \) can at most be of a \( 2p-2h \) nature if we have a two-body Hamiltonian. This means that second order in the perturbation theory can have \( 1p-1h \) and \( 2p-2h \) at most as intermediate states. When we diagonalize, these contributions are included to infinite order. This means that higher-orders in perturbation theory bring in more complicated correlations.

If we limit the attention to a Hartree-Fock basis, then we have that \( \langle\Phi_0\vert \hat{H}_I \vert 2p-2h\rangle \) is the only contribution and the contribution to the energy reduces to $$ \Delta E^{(2)}=\frac{1}{4}\sum_{abij}\langle ij\vert \hat{v}\vert ab\rangle \frac{\langle ab\vert \hat{v}\vert ij\rangle}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}. $$

If we compare this to the correlation energy obtained from full configuration interaction theory with a Hartree-Fock basis, we found that $$ E-E_0 =\Delta E= \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}, $$ where the energy \( E_0 \) is the reference energy and \( \Delta E \) defines the so-called correlation energy.

We see that if we set $$ C_{ij}^{ab} =\frac{1}{4}\frac{\langle ab \vert \hat{v} \vert ij \rangle}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}, $$ we have a perfect agreement between FCI and MBPT. However, FCI includes such \( 2p-2h \) correlations to infinite order. In order to make a meaningful comparison we would at least need to sum such correlations to infinite order in perturbation theory.

Summing up, we can see that