deletions | additions       diff --git a/Sternheimer.tex b/Sternheimer.tex  index 388e6f7..5407ec2 100644  --- a/Sternheimer.tex  +++ b/Sternheimer.tex 
...
\section{Sternheimer}  \label{sec:sternheimer}  David, Xavier In textbooks, perturbation theory is formulated in terms of sum over  states and response functions. These are useful theoretical  constructions that permit a good description and understanding of the  underlying physics. However, this is not a good description for  numerical applications, since it involves the calculation of a large  number of eigenvectors, infinite sums and the representation of  functions that depend on more than one spatial variable.  % (see for example  %\cite{Lup2008PRB}).  A very interesting approach, that essentially solves the problems  mentioned above, is the reformulation of standard perturbation theory  in terms of differential equations for the variation of the  wave-functions (to a given order), this is what is named in the  literature as Sternheimer equation~\cite{Ste1951PR}. Although a  perturbative technique, it avoids the use of empty states, and has a  quite good scaling with the number of atoms. This method has already  been used for the calculation of many response  properties~\cite{Bar2001RMP} like atomic vibrations (phonons),  electron-phonon coupling, magnetic response, etc. In the domain of  optical response, this method has been mainly used for static  response, although a few first-principles calculations for low  frequency (far from resonance) (hyper)polarisabilities have  appeared~\cite{Sen1987PRA,Kar1990CPL,Gis1997PRL,Iwa2001JCP}.  Recently, an efficient reformulation of the Sternheimer equation in a  super operator formalism was presented~\cite{Wal2006PRL}. When  combined with a Lanczos solver, it allows to calculate very  efficiently the first order polarisability for the whole frequency  spectrum. However, the generalization of this method to higher orders  is not straightforward.  Here, we present a modified version of the Sternheimer equation that is  able to cope with both static and dynamic response in and out of  resonance~\cite{And2007JCP}. The method is suited for linear and  non-linear response; higher order Sternheimer equations can be  obtained for higher order variations. For second-order response,  however, we can apply the \emph{2n\,+\,1 theorem} to get the  coefficients directly from first-order response variations. This  theorem states that the \(n\)th-order variations of the wave functions  (and the lower-order ones) are enough to obtain the $2n+1$ derivative  of the energy~\cite{Gon1989PRB,Bar2001RMP}.\footnote{We must keep in mind that the \(n\)th  derivative of the energy accounts for \((n-1)\)th order response, as  one order accounts for the outgoing field. For example, the static  polarisability is proportional to the second-order derivative of the  energy and corresponds to the linear (order 1) response  term.} The theorem also holds for the  dynamic case~\cite{Dal1996PRB}.  To derive the Sternheimer formalism we start by considering a  monochromatic perturbative field \(\lambda  \delta{\vv}(\vec{r})\cos\left(\omega{t}\right)\). If we assume that  the magnitude, \(\lambda\), is small, we can use perturbation theory to  expand the Kohn-Sham wave-functions in powers of \(\lambda\). To  first-order the KS orbitals read  %\begin{multline}  \begin{equation}  \label{eq:psi}  \bar{\varphi}_{\io}(\vec r, t) =  e^{-{\imi}\left(\epsilon_\io+\lambda\delta\epsilon_{\io}\right){t}}\Big[  \varphi_\io(\vec r) + \frac\lambda2e^{{\imi}\omega{t}}  \delta\varphi_{\io}(\vec r, \omega) +\frac\lambda2{e^{-{\imi}\omega{t}}}  \delta\varphi_{\io}(\vec r, -\omega)\Big] \,,  \end{equation}  %\end{multline}  where $\varphi_\io(\vec r)$ are the wave-functions of the static  Kohn-Sham Hamiltonian $\op{\HH}$ obtained by taking $\lambda=0$  \begin{equation}  \op{\HH}\varphi_\io(\vec r) =\epsilon_\io\varphi_\io(\vec r)  \,,  \end{equation}  and $\delta\varphi_{\io}(\vec r, \omega)$ are the first order variations of  the time-dependent Kohn-Sham wave-functions.  From Eq.~\eqref{eq:psi} and the definition of the time-dependent density,  Eq.~\eqref{eq:tdksden}, we can obtain the time-dependent density  %\begin{multline}  \begin{equation}  \label{eq:exprho}  \bar{n}\dert = n\der + \frac\lambda2  e^{{\imi}\omega{t}}\delta{n}(\vec r, \omega)  +\frac\lambda2{e^{-{\imi}\omega{t}}}\delta{n}(\vec r, -\omega)\,,  \end{equation}  %\end{multline}  \noindent with the definition of the first-order variation of the density  %\begin{multline}  \begin{equation}  \label{eq:varrho}  \delta{n}(\vec r, \omega) = \sum_\io^{\rm occ.} \Big\{  \left[\varphi_\io\der\right]^*\delta\varphi_{\io}(\vec r, \omega)\\  + \left[\delta\varphi_{\io}(\vec r, -\omega)\right]^*\varphi_\io\der  \Big\}\ .  \end{equation}  %\end{multline}  By replacing the expansion of the wave-functions \eqref{eq:psi} in the  time-dependent Kohn-Sham equation \eqref{eq:tdks}, and picking up the  first-order terms in $\lambda$, we arrive at a Sternheimer equation  for the variations of the wave-functions  %\begin{multline}  \begin{equation}  \label{eq:sternheimer}  \left\{\op{\HH} - \epsilon_\io\pm\omega +  \imi\eta\right\}\delta\varphi_{\io}(\vec r, \pm\omega) = \\  -\mathrm{P}_c\,\delta{\op{\HH}}(\pm\omega) \varphi_\io(\vec r)  \,,  \end{equation}  %\end{multline}  with the first order variation of the Kohn-Sham Hamiltonian  %\begin{multline}  \begin{equation}  \label{eq:h1}  \delta{\op{\HH}}(\omega)=  \delta{\op{\vv}}(\vec{r})  +\mint{r'} \frac{\delta{n}(\vec{r}',\omega)}{|\vec{r}-\vec{r}'|}  +\mint{r'} f_{\rm xc}(\vec r, \vec r', \omega)\,\delta{n}(\vec{r'}, \omega)  \ .  \end{equation}  \(\mathrm{P}_c\) is the projector onto the unoccupied subspace and  \(\eta\) a positive in\-fi\-ni\-te\-si\-mal, essential to obtain the  correct position of the poles of the causal response function, and,  consequently obtain the imaginary part of the  po\-la\-ri\-za\-bi\-li\-ty.\footnote{Furthermore, using a small but  finite \(\eta\) allows us to solve numerically the Sternheimer  equation close to re\-so\-nan\-ces, as it removes the divergences of  this equation.} The projector \(\mathrm{P}_c\) effectively removes  the components of \(\delta\varphi_{\io}(\vec r, \pm\omega)\) in the  subspace of the occupied ground-state wave-functions. In linear  response, these components do not contribute to the variation of the  density\footnote{This is straightforward to prove by expanding the  variation of the wave-functions in terms of the ground-state  wave-functions, using standard perturbation theory, and then  replacing the resulting expression in the variation of the density,  Eq.~(\ref{eq:varrho}).}, and therefore we can safely ignore the  projector for first-order response calculation. This is important for  large systems as the cost of the calculation of the projections scales  quadratically with the number of orbitals.  The first term of \(\delta{\op{\HH}}(\omega)\) comes from the external  perturbative field, while the next two represent the variation of the  Hartree and exchange-correlation potentials. The exchange-correlation  kernel is a functional of the ground-state density $n$, and is  given, in time-domain, by the functional derivative  \begin{equation}  f_{\rm xc}[n](\vec r, \vec r', t - t') =   \frac{\delta v_{\rm xc}[n](\vec r, t)}{\delta n(\vec r', t)}  \ .  \end{equation}  Equations \eqref{eq:varrho} and \eqref{eq:sternheimer} form a set of  self-consistent equations for linear response that only depend on the  occupied ground-state orbitals.