# 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,

$\Pi(BA\vert \omega) = \lim_{\eta \rightarrow 0} \left(\frac{1}{\hbar}\right) \sum\limits_{n\neq0} \left\{ \frac{\langle 0 \vert B \vert n\rangle\langle n \vert A \vert 0\rangle}{\omega + i\eta - \omega_{0n}} - \frac{\langle 0 \vert A \vert n\rangle \langle n \vert B \vert 0 \rangle}{\omega + i\eta + \omega_{0n}} \right\} \ \ \ \ \ (1)$

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:

$\hat{H}B = \left[H,B\right], \qquad \hat{H}^2B = \left[H,\left[H,B\right]\right], \qquad \hat{H}^3B = \left[H,\left[H,\left[H,B\right]\right]\right], \qquad \cdots \ \ \ \ \ (8)$

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.