# Derivation of general polarization propagator methods

Most higher-order response methods use some form of the polarization propagator, which is what I intend to derive here. It can then be used to derive other methods, such as TDHF/RPA and SOPPA.

We have derived earlier a response function, or frequency dependent polarizability,

where ${A}$ is the applied perturbation, and ${B}$ is the observable, and both are assumed to be Hermitian. ${\omega_{0n}}$ is the excitation energy for the change between states ${0}$ and ${n}$. It should be clear that the response function has poles when ${\omega}$ — the applied field frequency – equals to the excitation energy ${\omega_{0n}}$. Finding these poles is precisely the goal of polarization propagator methods. In the polarization propagator approach, the above equation has ${\eta}$ set to 0, and the response function (the `propagator’), defined as:
$\langle \langle B;A \rangle \rangle_{\omega} \equiv \sum\limits_{n\neq0} \left\{ \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega - \hbar\omega_{0n}} + \frac{\langle 0 \vert A \vert n\rangle \langle n \vert B \vert 0 \rangle}{-\hbar\omega - \hbar\omega_{0n}} \right\} \ \ \ \ \ (2)$

Now we want to describe the propagator in terms of commutators between ${A}$ and ${B}$. Make the observation that ${\frac{ab}{c+d} = \frac{ab}{c} - \frac{d}{c}\left(\frac{ab}{c+d}\right)}$, and applying to the first term of the above yields:
$\sum\limits_{n\neq0} \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega - \hbar\omega_{0n}} = \sum\limits_{n\neq0} \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega} -\sum\limits_{n\neq0} \frac{-\hbar\omega_{0n}\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\hbar\omega\left(\hbar\omega - \hbar\omega_{0n}\right)} \ \ \ \ \ (3)$

Do the same for the second term and combine, recognizing that the ${n=0}$ term vanishes in the first part (thus we get a sum over all ${n}$), and making use of the fact that ${1 = \sum\limits_{n}\vert n\rangle\langle n\vert }$ and ${H\vert n \rangle = E_n\vert n\rangle}$ and ${\hbar\omega_{0n} = E_n - E_0}$:
$\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \langle 0 \vert \left[B,A\right] \vert 0 \rangle + \frac{1}{\hbar\omega} \sum\limits_{n\neq0} \left\{\frac{\langle 0 \vert B \vert n\rangle\langle n \vert \left[H,A\right] \vert 0\rangle}{\hbar\omega - \hbar\omega_{0n}} + \frac{\langle 0 \vert \left[H,A\right] \vert n\rangle\langle n \vert B \vert 0\rangle}{-\hbar\omega - \hbar\omega_{0n}}\right\} \ \ \ \ \ (4)$

Which is to say that
$\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \langle 0 \vert \left[B,A\right] \vert 0 \rangle + \frac{1}{\hbar\omega}\langle \langle B;\left[H,A\right] \rangle \rangle_{\omega} \ \ \ \ \ (5)$

Or, as we will use it:
$\hbar\omega\langle \langle B;A \rangle \rangle_{\omega} = \langle 0 \vert \left[B,A\right] \vert 0 \rangle - \langle \langle \left[H,B\right] A \rangle \rangle_{\omega} \ \ \ \ \ (6)$

As you may have started to see, we can define the propagator iteratively in terms of commutator expectation values of ever-increasing complexity. This is what is known as the so-called ‘‘moment expansion’’ of the propagator. Thus by iteration:
$\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \left\{ \langle 0 \vert \left[B,A\right] \vert 0 \rangle + \left(\frac{-1}{\hbar\omega}\right)\langle 0 \vert \left[\left[H,B\right],A\right] \vert 0 \rangle + \left(\frac{-1}{\hbar\omega}\right)^2\langle 0 \vert \left[\left[H,\left[H,B\right]\right],A\right] \vert 0 \rangle + \cdots \right\} \ \ \ \ \ (7)$

We introduce the ‘‘superoperator’’ (analogous to the Liouville operator in Statistical Mechanics), which acts on operators to give their commutator:

With this definition, we have the power series
$\langle \langle B;A \rangle \rangle_{\omega} = \frac{1}{\hbar\omega} \sum\limits_{n=0}^{\infty} \left(\frac{-1}{\hbar\omega}\right)^n \langle 0 \vert \left[\hat{H}^nB,A\right]\vert 0\rangle \ \ \ \ \ (9)$

At this point we make two useful observations. First, recognize that
$\langle 0 \vert \left[\hat{H}B,A\right]\vert 0\rangle = -\langle0\vert \left[B,\hat{H}A\right]\vert 0\rangle \ \ \ \ \ (10)$

and so ${\hat{H}}$ can be applied to ${A}$ instead of ${B}$ insofar as we introduce a factor of ${(-1)^n}$. Furthermore, note that the power series is equivalent to
$\frac{1}{1-x} = 1 + x + x^2 + x^3 + \cdots \ \ \ \ \ (11)$

Making use of these two observations (and using ${\hat{1}X = X}$ and ${\hat{H}^0 = \hat{1}}$, where ${\hat{1}}$ is the unit superoperator), we have
$\langle \langle B;A \rangle \rangle_{\omega} = \langle 0 \vert \left[B,\left(\hbar\omega\hat{1} - \hat{H}\right)^{-1}A\right]\vert 0\rangle \ \ \ \ \ (12)$

Which is merely a cosmetic change at this point, as the superoperator resolvent is defined by the series expansion. We need to find a matrix representation of the resolvent, which implies that we find a complete basis set of operators. To do this, we are going to develop an operator space, where ${\hat{H}}$ is defined by its effect on operators instead of vectors. Introducing the notation
$\left(X\vert Y\right) = \langle 0 \vert \left[X^{\dagger},Y\right] \vert 0 \rangle \ \ \ \ \ (13)$

and it follows that ${\left(Y\vert X\right) = \left(X\vert Y\right)^*}$. As defined, we now have
$\langle \langle B;A \rangle \rangle_{\omega} = \left(B^{\dagger}\vert \left(\hbar\omega\hat{1} - \hat{H}\right)^{-1}\vert A \right) \ \ \ \ \ (14)$

Which is formally exact, albeit useless until we develop approximations. However, the form of the above equation does look similar to ordinary vector spaces in Hartree-Fock, etc. methods. Truncation of a basis in linear vector space ${V}$ to ${n}$ elements produces a subspace ${V_n}$, and truncation of a general vector corresponds to finding its projection onto the subspace. It follows, then, that we need to find a projection operator ${\rho}$, associated with the truncated basis. If the basis (${\mathbf{e}}$, say) is orthonormal we write
$\rho = \sum\limits_i e_ie_i^* = \mathbf{ee}^{\dagger} \ \ \ \ \ (15)$

which in a complete basis gives:
$\rho = \sum\limits_i \vert e_i\rangle \langle e_i\vert = \mathbf{1} \ \ \ \ \ (16)$

If it is not an orthonormal basis, we must include the metric matrix ${\mathbf{S} = \mathbf{e}^{\dagger}\mathbf{e}}$ (or L"{o}wdin basis ${\mathbf{\overline{e}}\mathbf{S}^{-1/2}}$):
$\rho = \mathbf{eS}^{-1}\mathbf{e}^{\dagger} = \sum\limits_{ij} e_i (S^{-1})_{ij}e_j^* \ \ \ \ \ (17)$

When using a truncated basis in operator space, two kinds of projections are useful (Löwdin, 1977, 1982),
$A' = \rho A \rho, \qquad A'' = A^{1/2} \rho A^{1/2} \ \ \ \ \ (18)$

which are the outer projection and inner projection, respectively, onto space ${V_n}$ defined by ${\rho}$. Note that $AB = C$ does not imply $A'B' = C'$ does not imply $A''B'' = C''$. Plugging the metric into ${A''}$:
$A'' = A^{1/2}\mathbf{eS}^{-1}\mathbf{e}^{\dagger}A^{1/2} \ \ \ \ \ (19)$

and we define
$\mathbf{f} \equiv A^{1/2}\mathbf{e} = \left(A^{1/2}\mathbf{e}_1 \quad A^{1/2}\mathbf{e}_1 \quad \cdots \right) \ \ \ \ \ (20)$

We assume that ${A}$ is Hermitian and positive-definite, so that ${A^{1/2}}$ can be defined. Note that ${\mathbf{S} = \mathbf{e}^{\dagger}\mathbf{e} = \left(A^{-1/2}\mathbf{f}\right)^{\dagger}\left(A^{-1/2}\mathbf{f}\right) = \mathbf{f}^{\dagger}A^{-1}\mathbf{f} \implies A'' = \mathbf{f}\left(\mathbf{f}^{\dagger}A^{-1}\mathbf{f}\right)^{-1}\mathbf{f}^{\dagger}}$. Because ${A}$ is arbitrary, replace it with ${A^{-1}}$, and since ${\mathbf{f}^{\dagger}A\mathbf{f} = \mathbf{A}}$ with ${A_{ij} = \langle f_i \vert A\vert f_j\rangle}$:
$\left(\mathbf{A}^{-1}\right)'' = \mathbf{f}\left(\mathbf{f}^{\dagger}A\mathbf{f}\right)^{-1}\mathbf{f}^{\dagger} = \mathbf{fA}^{-1}\mathbf{f}^{\dagger} \ \ \ \ \ (21)$

As the basis ${V_n \rightarrow V}$, the inner projection ${\rightarrow \mathbf{A}^{-1}}$, else it is simply a finite basis approximation to the inverse. This is the operator inverse in terms of a matrix inverse. Since ${\mathbf{e}}$ was an arbitrary basis defining ${V_n}$, let ${\mathbf{f}}$ define n-dimensional subspace ${V_n'}$. Thus:
$\mathbf{A}^{-1} \approx \mathbf{eA}^{-1}\mathbf{e}^{\dagger} = \sum\limits_{ij} e_i(\mathbf{e}^{\dagger}\mathbf{Ae})^{-1}_{ij}\mathbf{e}^*_j \ \ \ \ \ (22)$

Thus the inner projection leads to an approximation for the projector. Let us define the (as of yet undefined) operator basis:
$\mathbf{n} = (\mathbf{n}_1 \quad \mathbf{n}_2 \quad \mathbf{n}_3 \quad \cdots) \ \ \ \ \ (23)$

Given that the binary product (compare with ${\langle x\vert y\rangle = x^*y}$ for vectors)
$X^{\dagger} \cdot Y = (X\vert Y) = \langle 0 \vert \left[X^{\dagger},Y\right] \vert 0 \rangle \ \ \ \ \ (24)$

then for our resolvent superoperator we have
$\hat{R}(\omega) = (\hbar\omega\hat{1} - \hat{H})^{-1} = \mathbf{n}R(\omega)\mathbf{n}^{\dagger} = \sum\limits_{r,s} \mathbf{n}_r[R(\omega)]_{rs}\mathbf{n}_{s}^{\dagger} \ \ \ \ \ (25)$

where ${\mathbf{n}_r}$ and ${\mathbf{n}_s^{\dagger}}$ are the analogues of ${\mathbf{e}}$ and ${\mathbf{e}^*}$ in operator space. Finally, if
$R(\omega) = M(\omega)^{-1}, \qquad M(\omega) = \mathbf{n}^{\dagger} (\hbar\omega\hat{1} - \hat{H})\mathbf{n} \ \ \ \ \ (26)$

then we have
$\langle \langle B;A \rangle \rangle_{\omega} = B\cdot R(\omega)\cdot A = B \cdot \mathbf{n} R(\omega) \mathbf{n}^{\dagger} \cdot A \ \ \ \ \ (27)$

which is the key to calculating approximations to response properties. The matrix M is determined once we have chosen an operator basis. This approximation depends on two things: 1) the basis ${\mathbf{n}}$ (and its truncations), and 2) the reference function, that is not the exact ground state. Any approximations to these two things are where we get out various response methods.