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{aligned} \psi(r) \sim \frac{\mathrm{e}^{{\mathrm i}k r}}{r}\quad\textrm{as}\quad r\rightarrow\infty\ ,\end{aligned}\] 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 complex-scaling method \cite{aguilarcombes,Balslev:1971ez} where states and operators are represented by means of the transformation \[\begin{aligned} \hat R_\theta\, \psi(\vec r) = {\mathrm e}^{{\mathrm i}N \theta / 2} \psi(\vec r {\mathrm e}^{{\mathrm i}\theta})\ ,\end{aligned}\] where \(N\) is the number of spatial dimensions to which the scaling operation is applied, and \(\theta\) is a fixed 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 - {\mathrm i}\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
\[\begin{gathered} E_\theta = {\mathrm e}^{-{\mathrm i}2\theta} \sum_n \int{\mathrm d}\vec r\, \varphi_{\theta n}(\vec r) \left(-\frac12 \nabla^2\right) \varphi_{\theta n}(\vec r)\\ + {\mathrm e}^{-{\mathrm i}\theta} \frac12 \iint {\mathrm d}\vec r \, {\mathrm d}{\vec{r}'}\, \frac{n_\theta(\vec r)n_\theta({\vec{r}'})}{\left| \vec r - {\vec{r}'}\right|}\\ \quad+ E_{\mathrm{xc}}^\theta[n_\theta] + \int{\mathrm d}\vec r\, v_{\mathrm{ext}}(\vec r {\mathrm e}^{{\mathrm i}\theta}) n_\theta(\vec r)\ ,\end{gathered}\]
with the now-complex electron density \[\begin{aligned} n_\theta(\vec r) = \sum_n f_n \varphi_{\theta n}^2(\vec r)\ ,\end{aligned}\] 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{aligned} \left[-\frac12 {\mathrm e}^{-{\mathrm i}2\theta}\nabla^2 + v_\theta(\vec r) \right] \varphi_{\theta n}(\vec r) = \varphi_{\theta n}(\vec r) \epsilon_{\theta n}\ .\end{aligned}\] 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{aligned} v_\theta(\vec r) \equiv {\frac{\delta E}{\delta n_\theta(\vec r)}} = v_{\mathrm H}^\theta(\vec r) + v_{\mathrm{xc}}^\theta(\vec r) + v_{\mathrm{ext}}(\vec r {\mathrm e}^{{\mathrm i}\theta})\ ,\end{aligned}\] where \(v_{\mathrm{ext}}(\vec r {\mathrm e}^{{\mathrm i}\theta})\) may represent atomic potentials as analytically continued pseudopotentials, and where the Hartree potential \[\begin{aligned} v_{\mathrm H}^\theta(\vec r) &= {\mathrm e}^{-{\mathrm i}\theta}\int{\mathrm d}{\vec{r}'}\, \frac{n_\theta(\vec r')}{\left| \vec r'-\vec r\right|}\end{aligned}\] is determined by solving the Poisson equation defined by the complex density. Together with the xc potential, \[\begin{aligned} v_{\mathrm{xc}}^\theta(\vec r) &= {\frac{\delta E_{\mathrm{xc}}^\theta[n_\theta]}{\delta n_\theta(\vec r)}},\end{aligned}\] 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}.