\section{Complex scaling and resonances}
% Commands to unify notation (vectors, constants, etc.) in this section.

\newcommand{\ee}{\mathrm e} % 2.718281828...
\newcommand{\ii}{\mathrm i} % sqrt -1
\newcommand{\xc}{\mathrm{xc}} % exchange--correlation
\newcommand{\Ha}{\mathrm H} % hartree
\newcommand{\dee}{\mathrm d} % dx, dr
\newcommand{\fdiff}[2]{\frac{\delta #1}{\delta #2}} % functional derivative

In this section we discuss the calculation of resonant electronic states by
means of the complex-scaling method, as implemented in Octopus.  By
``resonant states,'' we mean metastable electronic states
of finite systems, such as atoms or molecules, with a characteristic energy
and lifetime.

Mathematically, resonances can be defined as poles of the scattering
matrix or cross-section at complex energies~\cite{PhysRev.56.750,Hatano01022008}.
If a pole is close to the
real energy axis, it will produce a large, narrow peak in the
cross-section of scattered continuum states.
One way resonances can arise is from application of an
electric field strong enough to ionize the system through tunnelling.
Resonant states may temporarily capture incoming electrons or electrons
excited from bound states, making them important intermediate states
in many processes.


The defining characteristic of a resonant state, often called a Siegert
state~\cite{PhysRev.56.750}, is that it has an
outgoing component but not an incoming one.
They can be determined by solving the
time-independent Schrödinger
equation with the boundary condition that the wavefunction must asymptotically
have the form
\begin{align}
  \psi(r) \sim \frac{\mathrm{e}^{\ii k r}}{r}\quad\textrm{as}\quad r\rightarrow\infty\ ,
\end{align}
where the momentum $k$ is complex and has a negative imaginary part.
This causes the state to diverge exponentially in space as
$r\rightarrow\infty$.  The state can further be ascribed a complex
energy, likewise with a negative imaginary part, causing it to decay
over time at every point in space uniformly.

Resonant states are not eigenstates of any Hermitian operator and 
in particular do not reside within the Hilbert space.  This precludes their direct
calculation with the standard computational methods
from DFT.  However, it turns out that a suitably chosen analytic
continuation of a Siegert state is localized, and this form can be used
to derive information from the state.
This is the idea behind the
\emph{complex-scaling method} \cite{aguilarcombes,Balslev:1971ez} where
states and operators are represented
by means of the transformation
\begin{align}
  \hat R_\theta\, \psi(\vec r) = \ee^{\ii N \theta / 2} \psi(\vec r \ee^{\ii\theta})\ ,
\end{align}
where $N$ is the number of spatial dimensions to which the scaling operation
is applied, and $\theta$ is a fixed \emph{scaling angle} which determines
the path in the complex plane along which the analytic continuation is taken.
The transformation maps the Hamiltonian to a non-Hermitian operator
$\hat H_\theta = \hat R_\theta \hat H \hat R_{-\theta}$.

The Siegert states $\psi(\vec r)$ of the original Hamiltonian are square-integrable
eigenstates $\psi_\theta(\vec r)$
of $\hat H_\theta$, and their eigenvalues $\epsilon_0 - \ii\Gamma/2$
define the energy $\epsilon_0$ and width $\Gamma$ of the
resonance~\cite{simon1973resonances,Reinhardt_1982,Ho19831}.

A typical example of a spectrum of the transformed Hamiltonian $\hat H_\theta$
is shown in Fig.~\ref{fig:cs-spectrum}, and the corresponding potential
and lowest bound and resonant states in Fig.~\ref{fig:cs-potential-wfs}.
The bound-state energies are unchanged while the continuum rotates by
$-2 \theta$ around the origin.  Finally, resonances appear as isolated
eigenvalues in the fourth quadrant once $\theta$ is sufficiently large
to ``uncover'' them from the continuum.
Importantly, matrix elements (and in
particular energies) of states are independent of $\theta$ as long as
the states are localized and well represented numerically --- this
ensures that all physical
bound-state characteristics of the untransformed Hamiltonian are retained.

Our implementation supports calculations with complex scaling 
for independent particles or
in combination with DFT and
selected
xc functionals~\cite{Larsen:2013cw}.
The energy functional in KS-DFT consists of several
terms that are all expressible as integrals of the density or the
wavefunctions with the kinetic operator and various potentials.
The functional is complex-scaled as per the prescribed method
by rotating the real-space integration contour
of every term by $\theta$ in the complex plane.
The DFT energy functional becomes
%
\newcommand{\rprime}[0]{\vec{r}'}
\begin{multline}
  E_\theta = \ee^{-\ii2\theta}
  \sum_n \int\dee\vec r\, \varphi_{\theta n}(\vec r) \left(-\frac12 \nabla^2\right)
\varphi_{\theta n}(\vec r)\\
  + \ee^{-\ii\theta} \frac12
\iint \dee \vec r \, \dee \rprime \,
\frac{n_\theta(\vec r)n_\theta(\rprime)}{\left| \vec r - \rprime \right|}\\
\quad+ E_\xc^\theta[n_\theta]
+ \int\dee \vec r\, v_{\mathrm{ext}}(\vec r \ee^{\ii \theta}) n_\theta(\vec r)\ ,
\end{multline}
%
with the now-complex electron density
\begin{align}
  n_\theta(\vec r) = \sum_n f_n \varphi_{\theta n}^2(\vec r)\ ,
\end{align}
with occupation numbers $f_n$, and complex-scaled KS states $\varphi_{\theta n}(\vec r)$.
Note that no complex conjugation is performed on the left component in
matrix elements such as the density or kinetic energy.
In order to define the complex-scaled xc potential, it is necessary to perform an analytic continuation procedure~\cite{Larsen:2013cw}.

In standard DFT, the KS equations
are obtained by taking the functional derivative of the energy functional with respect
to the density.  Solving the equations corresponds to searching
for a stationary point, with the idea that this minimizes the energy.
In our case, since the energy functional is complex-valued~\cite{WM07}, we cannot minimize the energy functional, but we can still search for stationary points to find the resonances~\cite{Whitenack_2010,WW11}.
The complex-scaled version
of the KS equations thereby becomes similar to the usual ones:
\begin{align}
  \left[-\frac12 \ee^{-\ii2\theta}\nabla^2 + v_\theta(\vec r)
    \right] \varphi_{\theta n}(\vec r) = \varphi_{\theta n}(\vec r) \epsilon_{\theta n}\ .
\end{align}
The effective potential $v_\theta(\vec r)$ is the functional derivative
of the energy functional with respect to the density $n_\theta(\vec r)$, and, therefore,
consists of the terms
\begin{align}
  v_\theta(\vec r) \equiv \fdiff{E}{n_\theta(\vec r)} = 
  v_\Ha^\theta(\vec r) + v_\xc^\theta(\vec r) 
  + v_{\mathrm{ext}}(\vec r \ee^{\ii\theta})\ ,
\end{align}
where $v_{\mathrm{ext}}(\vec r \ee^{\ii\theta})$ may represent atomic
potentials as analytically
continued pseudopotentials, and where the Hartree potential
\begin{align}
  v_\Ha^\theta(\vec r) &=
  \ee^{-\ii\theta}\int\dee \rprime\, \frac{n_\theta(\vec r')}{\left| \vec r'-\vec r\right|}
\end{align}
is determined by solving the Poisson equation defined by the complex density.
Together with the xc potential, 
\begin{align}
  v_\xc^\theta(\vec r) &= \fdiff{E_\xc^\theta[n_\theta]}{n_\theta(\vec r)},
\end{align}
this defines a self-consistency cycle very similar to ordinary
KS DFT, although more care must be taken to occupy the correct states,
as they are no longer simply ordered by energy.

Fig.~\ref{fig:cs-ionization-He} shows calculated ionization rates for the He~1$s$~state in a
uniform Stark-type electric field as a function of field strength.
In the limit of weak electric fields, the simple approximation
by Ammosov, Delone and Krainov (ADK)~\cite{adk}, which depends
only on the ionization potential, approaches the accurate reference calculation
by Scrinzi and co-workers~\cite{PhysRevLett.83.706}.
This demonstrates that the ionization rate is determined largely by the
ionization potential for weak fields.  As the local density approximation is known to produce inaccurate ionization potentials
due to its wrong asymptotic form at large distances, it necessarily yields
inaccurate rates at low fields.
Meanwhile exact exchange, which is known to produce
accurate ionization energies, predicts ionization rates
much closer to the reference calculation.  The key property 
of the xc functional that allows accurate determination of decay rates
from complex-scaled DFT therefore appears to be that it must yield
accurate ionization potentials, which is linked to its ability
to reproduce the correct asymptotic form of the potential at large distances
from the system~\cite{PhysRevA.49.2421}.
