Application of sampling theory in modelling of continuum processes: photoionization cross-sections of atoms

Photoionization processes play an important role in chemical physics as a fundamental probe of the structural and dynamical properties of many-body systems. Photoelectron spectroscopy is often applied to identify the electronic structure of materials in different fields as well as in different contexts. For example, photoelectron spectroscopy is used in astrophysics, aeronomy, radiation chemistry, environmental and atmospheric chemistry, metrology, surface science and catalysis and new material development in order to obtain detailed information about the energy levels and electronic wavefunctions of molecules and ions, as well as their interactions with radiation.

The development of theoretical models of electrodynamical processes in complex systems is driven primarily by the development in experimental techniques. Synchrotron radiation sources and the recent emergence of high-brightness x-ray free electron laser sources has stimulated renewed research activity in the area of molecular photoionization [1–3], which impacts the achievable resolution of methods of structure determination in molecules and clusters. A theoretical method suitable for the calculation of electrodynamical processes in molecules or clusters ultimately depends on the availability of the appropriate transition matrix elements involving continuous states, obtained either by detailed calculation or by some approximation scheme.

One of the unresolved issues in scattering theory in general, and photoionization theory in particular, involves the representation of the photoelectron wavefunction in a continuum state. The available theoretical methods for calculating the transition matrix can be broadly characterized as either employing a projection of the continuum wavefunction onto a square integrable (L²) basis set or the construction of a true continuum wavefunction with appropriate boundary conditions [4, 5]. For molecules, the lack of spherical symmetry and the need to generate delocalized bound-state functions as initial states has led to a highly developed technology for generating electronic structures using L² basis sets, typically involving Gaussian functions. This technology is almost wholly independent of any method that attempts to construct a true scattering type solution, so performing all determinations of electrodynamic processes for molecules within an L² basis seems to be preferable if it can be made to be feasible and accurate.

The simplest model within the assumptions of strong orthogonality and the use of a single-centre expansion of the photoelectron wavefunction to impose appropriate spherical boundary conditions far from the interacting system involves either the use of a plane or Coulomb wave representation. Some recent developments implement [6–8] this simple model along with correlated Dyson orbital, representing the initial and final state of the system. This method provides a reasonably accurate molecular photoinization cross sections of polyatomic systems. The exchange interaction of the photoelectron with the residual core in this model is accounted for by the use of a partial effective charge of the core. While this procedure is motivated by physical intuition the selection of the effective charge parameters is largely empirical.

A mathematical technique that exploits L² wavefunctions is based on Stieltjes moment theory, which is not directly related to the explicit representation of the continuum state. In this method, the final state of the photoelectron wavefunction is sampled by the discretized representation of the complete spectrum. The discrete transition matrix elements obtained from these L² states are used to extract the oscillator strength density in the continuum. In this way, the Stieltjes method avoids the need for the solution of the scattering equations by extracting cross sections from spectral moments of an oscillator strength distribution. The discretized samples of the continuum states are the excited states calculated using any L² method. Any electronic structure method that generates a discretized representation of the complete spectrum can be used, in principle, to calculate the photoionization cross section, of any other continuum process, using the Stieljes imaging method. This approach has been led by Langhoff and collaborators [9–12]. Further developments and applications of this method based on various ab initio electronic structure theory have been performed recently [13–16]. A number of other approaches based on L² discretizaton methods have also been reported, including the complex coordinate methods [17] and the reactance or K-matrix method [18].

The accurate representation of the true continuum states in the field of residual molecular region can be achieved by allowing non-orthogonality of the photoelectron wavefunction and occupied orbitals [19], and exploiting variational techniques or by considering extra bound electron configurations in the wavefunctions, such as is implemented in R-matrix methods [20]. The variational principle provides a powerful way of solving a wide range of differential and integral equations involved in calculating collision and photoionization cross sections [21]. The Kohn variational principle is probably the most widely used variational principle in describing collision phenomena [22], though there are at least two major difficulties associated with implementations of this method. The first of these is the occurrence of anomalous singularities in the reactants or K-matrices that make these techniques difficult to apply in large scale calculations. The second problem concerns computational difficulties associated with the evaluation of multi-centre integrals involving both bound and continuum functions [23–26]. These issues have been addressed since the late 1980s and have been further developed into the complex Kohn method for electron-molecule scattering in recent years.

McKoy and co-workers [27–32] have developed methods that are based on, or in some way related to, the Schwinger variational principle or extensions of it that rely on the Lippmann–Schwinger scattering equation. The major disadvantage of this method is the occurrence of the Green's function in this variational principle which is known in closed form only in some special cases, such as the free-particle case. This method does have some advantages, however, including automatic incorporation of the correct boundary conditions.

For simple atoms and diatomic molecules continuum states can be reconstructed using L² method based on B-splines [33]. The latter are known to give accurate representation of continuum wavefunctions within the 'box' after normalization factor is adjusted [34]. Although the results obtained via B-splines demonstrate high accuracy in reconstructing continuum states, the B-splines themselves are hardly applicable for description of molecular systems beyond few atoms in size.

In the present work, however, we have applied sampling theory methods originally developed for the signal processing and information theory to construct matrix elements involving continuum state amplitudes from discrete representations of the complete set of eigenfunctions of the electronic spectrum. We also investigate the perspectives of using a Gaussian basis set for the description of the matrix elements that involve continuum final state. The calculations were carried out for hydrogen and sodium atoms. These systems are well-studied and are widely used as a benchmark for photoionization calculations. Another model, discussed in detail, involves a hydrogenic initial state and a final state that is free particle. This model allows a simple analytic solution and may be used for illustrative purposes. The proposed approach significantly simplifies the computational labour. The motivation behind it is an application to molecular systems where sampling of a discrete representation of the complete spectrum is effectively mandatory and in majority of cases involves a Gaussian basis set.

The photoionization cross-section in the dipole approximation is given by [35]

$$\begin{eqnarray}&&\sigma (\omega )=2\pi \alpha \omega \displaystyle \sum _{f}| {\langle i| \vec{r}| f\rangle | }^{2},\end{eqnarray} \tag{ 1 }$$

where α is the fine structure constant, ω is the photon frequency, $| i\rangle$ is the initial state of the system, and $| f\rangle$ is the final state. Conservation of energy, ${\hslash }\omega ={E}_{{f}}-{E}_{{i}}$, is assumed to be fulfilled while the effect of recoil is ignored. In the case of atoms and molecules, the initial state in equation (1) is a bound state that can be calculated using number of well established methods. For simple systems, such as light atoms and molecules, K-matrix [36, 37] and R-matrix [38] methods can be employed to calculate the final state wavefunction, as well as convergent close-coupling [39, 40]. In general, however, these methods are limited by the rapid growth of computational complexity with the increase in the number of atoms and technical difficulties that arise through the loss of spherical symmetry. Hybrid basis set approaches have been proposed to tackle this problem [41], which aim to provide accurate representations of the continuum wavefunction within a finite box using B-spline functions and appropriate boundary conditions. This flexibility is achieved, however, at the expense of a very large basis set dimension and the need to calculate new types of interaction integral, which complicates its implementation within existing quantum chemistry packages.

Rather than attempting to construct an accurate representation of a molecular continuum state wavefunction, it seems preferable to deduce the values of matrix elements involved in equation (1) by some other means that is more consistent with the L² methods that are used to generate the target bound-states. As will be shown below, this can be achieved with the help of conventional L² Gaussian basis set that is widely used in the molecular structure calculations that form the practical basis of much of quantum chemistry. Such calculations are of particular importance for investigation of photoionization in medium to large size molecules which takes place in a number of applications such as the radiation damage processes that influence the quality of single molecule imaging algorithms [42].

Briefly, these L² methods solve the generalized matrix eigenvalue equations

$$\begin{eqnarray}&&{{\bf{Hc}}}_{\epsilon }={E}_{\epsilon }{{\bf{Sc}}}_{\epsilon },\end{eqnarray} \tag{ 2 }$$

where ${H}_{\mu \nu }=\langle \mu | \hat{H}| \nu \rangle$ is a matrix element of Hamiltonian, $\hat{H}$ of the system, $| \mu \rangle$ and $| \nu \rangle$ are the basis states, ${c}_{\mu \epsilon }$ is an element of the eigenvector, ${{\bf{c}}}_{\epsilon }$, ${E}_{\epsilon }$ is an eigenvalue and ${S}_{\mu \nu }=\langle \mu | \nu \rangle$ is an overlap matrix elements involving $| \mu \rangle$ and $\nu \rangle$. The eigenstates of the Hamiltonian are $| \epsilon \rangle ={\sum }_{\mu }{c}_{\mu ,\epsilon }| \mu \rangle$. Within the Hartree–Fock approximation, the system of equations defined by (2) is referred to as the Hartree–Fock–Roothaan equations or, simply, the Roothaan equations [43].

Use of L² discretized 'pseudo-states' as a basis to extract information about continuum processes is not a new idea. For example Heller, Reinhardt, and Yamani [44] applied this method using a matrix representation of the free-particle and hydrogenic Hamiltonians in Laguerre–Sturmian basis. In these cases, basis set matrix diagonalization can be carried out analytically, allowing a comprehensive study to be performed. Here, we generalize their approach in a manner that does not rely on a particular choice of basis set but which does require, as is always the case in molecular electronic structure calculations, a numerical diagonalization of the matrix representation of the Hamiltonian.

Suppose one needs to estimate a matrix element $\langle i| \hat{A}| k\rangle$ of some operator $\hat{A}$ for the system of interest, for example an atom or a molecule. The final state of the photoelectron in a single electron approximation, $| k\rangle$, has positive energy, so that

$$\begin{eqnarray}&&\hat{H}| E\rangle =E| E\rangle ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\langle E^{\prime} | E\rangle =\delta (E^{\prime} -E),\end{eqnarray} \tag{ 4 }$$

where $\delta (x)$ is the Dirac delta function, $E={k}^{2}/2\gt 0$ is the asymptotic kinetic energy of the continuum state. In addition, we assume that (2) has been solved in some L² finite basis set, yielding the solutions

$$\begin{eqnarray}&&\hat{H}| m\rangle ={E}_{m}| m\rangle ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\langle n| m\rangle ={\delta }_{{nm}},\end{eqnarray} \tag{ 6 }$$

where ${\delta }_{{nm}}$ is Kronecker delta function. The continuum states $| m\rangle$ are identified as the 'positive-energy' solutions of the above equation, for which ${E}_{m}={k}_{m}^{2}/2$. Functions $| m\rangle$ are not the actual eigenstates $| E\rangle$, but rather some wavepacket representation of them, so we shall refer to these discretized positive energy states as 'pseudo-states'. The following derivation establishes the relation between the matrix elements of these pseudo-states and the matrix elements that involve actual eigenstates of corresponding positive energy $| {E}_{m}\rangle$. The identity (6) can then be rewritten as

$$\begin{eqnarray}&&{\int }_{0}^{\infty }\langle n| E\rangle \langle E| m\rangle {\rm{d}}{E}\approx {\delta }_{{nm}}.\end{eqnarray} \tag{ 7 }$$

In the case of the free-particle Hamiltonian the states $| k\rangle$ form a complete set and equation (7) is a strict equality. In other cases, however, the partitioning of the spectrum into 'positive' and 'negative' energy parts renders (7) an approximation. It proves to be, however, an excellent approximation if $| m\rangle$ and $| n\rangle$ are discrete states that are classified, on the basis of their respective energies, as being of positive-energy type. The integral in (7) can be re-written by changing the integration variable as

$$\begin{eqnarray}&&{\int }_{0}^{\infty }\langle n| E\rangle \langle E| m\rangle {\rm{d}}{E}={\int }_{0}^{\infty }\langle n| E\rangle \langle E| m\rangle \displaystyle \frac{{\rm{d}}{E}(\mu )}{{\rm{d}}\mu }{\rm{d}}\mu ,\end{eqnarray} \tag{ 8 }$$

where $E(\mu )$ is a smooth interpolation of points $({E}_{m},m)$ such that for integer values $\mu =m$ from 1 to N, the interpolation $E(\mu )$ takes values $E(m)={E}_{m}$ that correspond to positive-energy solutions of (5). The above equation can approximated by the following quadrature if one drops out part of the integral from 0 to $\infty$:

$$\begin{eqnarray}&&{\int }_{0}^{\infty }\langle n| E\rangle \langle E| m\rangle {\rm{d}}{E}\approx \displaystyle \sum _{l=1}^{N}\langle n| {E}_{l}\rangle \langle {E}_{l}| m\rangle {\omega }_{E}^{(l)},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\omega }_{E}^{(l)}={\left.\displaystyle \frac{{\rm{d}}{E}}{{\rm{d}}\mu }\right|}_{\mu =l},\end{eqnarray} \tag{ 10 }$$

where μ $\epsilon$ $[1,N]$. Note that in the right hand side of (9) ${\rm{d}}\mu \to {\rm{\Delta }}\mu =1$ has being included. Substituting equation (9) into equation (7) we arrive at the system of equations

$$\begin{eqnarray}&&\displaystyle \sum _{l=1}^{N}\langle n| {E}_{l}\rangle \langle {E}_{l}| m\rangle {\omega }_{E}^{(l)}={\delta }_{{nm}},\end{eqnarray} \tag{ 11 }$$

where $\langle n| {E}_{l}\rangle$ are unknowns. The solutions of these equations are of the form

$$\begin{eqnarray}&&\langle n| {E}_{l}\rangle =\displaystyle \frac{{\delta }_{{nl}}}{\sqrt{{\omega }_{E}^{(l)}}}={\delta }_{{nl}}{\left.\sqrt{\displaystyle \frac{{\rm{d}}\mu }{{\rm{d}}{E}}}\right|}_{E={E}_{l}}.\end{eqnarray} \tag{ 12 }$$

Equation (12) allows one to evaluate the matrix element $\langle i| \hat{A}| {E}_{n}\rangle$ if one knows the value of $\langle i| \hat{A}| n\rangle$:

$$\begin{eqnarray}&&\langle i| \hat{A}| {E}_{m}\rangle =\displaystyle \sum _{n=1}^{N}\langle i| \hat{A}| n\rangle \langle n| {E}_{m}\rangle =\langle i| \hat{A}| m\rangle {\left.\sqrt{\displaystyle \frac{{\rm{d}}\mu }{{\rm{d}}{E}}}\right|}_{E={E}_{m}},\end{eqnarray} \tag{ 13 }$$

where $\mu =\mu (E)$ is a smooth interpolation on the points $\{{E}_{n},n\}$. The above result generalizes a method described by [44] for a discrete set of wavepackets $\{| n\rangle \}$ obtained by discretizing a Hamiltonian in a finite basis set. The quality of the result obtained using (13) depends on the completeness of the basis set $\{| n\rangle \}$, which was established in the particular case of a complete Laguerre–Sturmian basis in [44, 45]. Similar relation is used in B-spline based methods [34] to adjust the normalization of the positive energy pseudo-states.

For the application of photoionization processes in (13) we take $\hat{A}=\hat{E}1$, where $\hat{E1}$ is the electric dipole transition operator. It is also convenient to use wavenumber normalization of the continuum state wavefunctions. For bound-continuum electric dipole matrix elements and the transition from a bound-state to a continuum pseudo-state of corresponding energy, we have

$$\begin{eqnarray}&&\langle i| \hat{E1}| {k}_{m}\rangle =\langle i| \hat{E1}| m\rangle {\left.\sqrt{2\pi \displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}}\right|}_{k={k}_{m}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\langle i| \hat{E1}| k\rangle | }_{k={k}_{m}}\equiv \langle i| \hat{E1}| {k}_{m}\rangle ,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{s}}{{k}{\rm{d}}{k}}=\displaystyle \frac{{\rm{d}}\mu }{{\rm{d}}{E}},\end{eqnarray} \tag{ 16 }$$

where k_m is the wavenumber corresponding to the eigenvalue E_m of pseudo-state $| m\rangle$, $s=s(k)$ is a smooth interpolation on the points $\{{k}_{m},m\}$ and $| i\rangle$ is the initial bound state of the atom or molecule. Note that ${\rm{d}}{s}/{\rm{d}}{k}{| }_{k={k}_{n}}{\rm{\Delta }}k$ can be interpreted as the number of pseudo-states in the interval ${\rm{\Delta }}k$ in the vicinity of $k={k}_{n}$ (see for example figure 1) and that, therefore, ${\rm{d}}{s}/{\rm{d}}{k}{| }_{k={k}_{n}}$ represents the density of pseudo-states in the vicinity of k_n. If the square-integrable basis set tends towards completeness, the density of pseudo-states tends to infinity for any wavenumber k.

Zoom In Zoom Out Reset image size

In order to calculate $\langle i| \hat{E1}| k\rangle$ for all values of k and not just those values equal to the wavenumbers k_m of pseudo-states, one can use a simple spline or other kind of polynomial interpolation. A more powerful way of calculating $\langle i| \hat{E1}| k\rangle$ within the approximations that have been imployed uses the Kramer sampling theorem [46]. Within this formalism, the matrix element for electric dipole transitions can be written as

$$\begin{eqnarray}&&\langle i| \hat{E1}| k\rangle =\mathop{\mathrm{lim}}\limits_{N\to \infty }\displaystyle \sum _{m=1}^{N}\langle i| \hat{E1}| m\rangle \langle m| k\rangle ,\end{eqnarray} \tag{ 17 }$$

where the completeness condition $\sum | m\rangle \langle m| =1$ has been used. Using (14) the above equation can be written as

$$\begin{eqnarray}&&\langle i| \hat{E1}| k\rangle =\mathop{\mathrm{lim}}\limits_{N\to \infty }\displaystyle \sum _{m=1}^{N}\langle i| \hat{E1}| {k}_{m}\rangle S(k,{k}_{m}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&S(k,{k}_{m})=\langle m| k\rangle {\left./\sqrt{2\pi \displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}}\right|}_{k={k}_{m}},\end{eqnarray} \tag{ 19 }$$

where $s=s(k)$ is defined as in (14). We shall refer to the functions $S(k,{k}_{m})$ as 'Kramer sampling functions' and to the functions $\langle m| k\rangle$, which are the projections of pseudo-states onto the continuum states, as 'sampling functions'. The Kramer sampling theorem formulated in the form of (18) allows one to calculate the matrix elements $\langle i| \hat{E1}| k\rangle$ for any value of the wavenumber k if one knows the values of matrix elements for some discrete set of points ${k}_{m}$ called sampling points. In a particular case when the sampling points are equally spaced $s(k)=k$ and the continuum states are free-particle spherical outgoing s-type waves confined in the spherical box of radius ${r}_{{\rm{\max }}}=2\pi$, the Kramer sampling formula reduces to well-known Whittaker–Shannon–Kotel'nikov sampling theorem and $S(k,{k}_{m})\sim \mathrm{sinc}(\pi k-\pi {k}_{n})$.

In the general case, the functions $S(k,{k}_{m})$ are unknown. In the limit of a complete L² basis functions, the functions defined in (19) satisfy the conditions

$$\begin{eqnarray}&&S({k}_{n},{k}_{m})={\delta }_{{nm}}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\int {S}^{* }(k,{k}_{m})S(k,{k}_{n}){\left.\displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}\right|}_{k={k}_{n}}{\rm{d}}{k}={\delta }_{{nm}}.\end{eqnarray} \tag{ 21 }$$

Some examples of these functions, their properties, as well as more rigorous and detailed formulation of Kramer sampling theorem can be found in [47, 48]. It is shown in [47] that the replacement of the exact functions $S(k,{k}_{m})$ given by (19) by functions which correspond to a different set of continuum states, leads to the same result in the limit of the infinite number of sampling points. In the present case we restrict our focus to atoms, so it is more convenient to use sampling functions derived from Bessel–Hankel transform [48] given by

$$\begin{eqnarray}&&\tilde{S}(y,{y}_{m})=\displaystyle \frac{2{({{yy}}_{m})}^{1/2}{J}_{\nu }(y)}{{J}_{\nu }^{\prime }({y}_{m})({y}^{2}-{y}_{m}^{2})},\end{eqnarray} \tag{ 22 }$$

where y_n are solutions of ${J}_{\nu }(y)=0$. We generate a set of function $\tilde{S}(y,{y}_{m})$ using Bessel–Hankel transform since this transform has a spherical Bessel function as a core. This corresponds to the radial part of solution of the free particle Hamiltonian in spherical coordinates. As it will be shown in this and subsequent sections, this is a good choice not only for free particle case, but for much wider class of Hamiltonians with physically realistic potentials. To construct suitable sampling function from (22) one can use the relation

$$\begin{eqnarray}&&S(k,{k}_{m})=\tilde{S}(y(k),{y}_{m})\sqrt{{N}^{2}\displaystyle \frac{{\rm{d}}{y}}{{\rm{d}}{k}}}{\left./\sqrt{\displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}}\right|}_{k={k}_{m}},\end{eqnarray} \tag{ 23 }$$

where $y=y(k)$ is a smooth interpolation of points $\{{k}_{n},{y}_{n}\}$, ${N}^{2}\approx 1/\pi$ is the normalization constant of $S(y,{y}_{m})$. This function can be shown to satisfy (20), but in general, it does not satisfy (21). In the above equation function the limit of dy/dk is $\pi {\rm{d}}{s}/{\rm{d}}{k}$ for $k\gg {k}_{1}$. For small wavenumbers, dy/dk depends strongly on the interpolation procedure employed and varies dramatically between different orders of Bessel function in (22). In calculations it is more convenient to construct sampling functions using

$$\begin{eqnarray}&&S(k,{k}_{m})=\tilde{S}(y(k),{y}_{m})\sqrt{\displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}}{\left./\sqrt{\displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}}\right|}_{k={k}_{m}},\end{eqnarray} \tag{ 24 }$$

where, as above, $y=y(k)$ and $s=s(k)$ provide smooth interpolations of the points $\{{k}_{n},{y}_{n}\}$ and $\{{k}_{n},n\}$, respectively. For the function defined by (24), condition (20) is satisfied exactly, while (21) holds only approximately, though to good accuracy, according to our numerical calculations. Summarizing these results, equations (18), (19), (22) and (24) may be used to construct an approximation to the matrix element of the electric dipole operator according to

$$\begin{eqnarray}&&\langle i| \hat{E1}| k\rangle =\displaystyle \sum _{n=1}^{N}\langle i| \hat{E1}| n\rangle \langle n| k\rangle ,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\langle n| k\rangle =\displaystyle \frac{2{(y(k){y}_{n})}^{1/2}{J}_{\nu }(y(k))}{{J}_{\nu }^{\prime }({y}_{n})({y}^{2}(k)-{y}_{n}^{2})}\sqrt{2\pi \displaystyle \frac{{\rm{d}}{s}}{{\rm{d}}{k}}}.\end{eqnarray} \tag{ 26 }$$

We denote the function $\langle n| k\rangle$ in (26) a 'Bessel–Hankel sampling function', to reflect its relation to the Kramer sampling function $S(k,{k}_{m})$ and the Bessel–Hankel transformation from which it is derived. The sign of $\langle n| k\rangle$ is undetermined, so for later convenience we impose the convention that $\langle r| n\rangle \gt 0$ for $r\to 0$.

It is convenient to start our analysis of photoionization using a simple model, in which initial state is a ground-state of hydrogen atom and final state is a free particle. This model has being discussed in [49] and treats potential energy as a perturbation. Its advantage is that most of the results can be obtained analytically, allowing direct comparison with our numerical approach. Due to spherical symmetry of the problem it is convenient to adopt spherical coordinates and integrate over the angular part to obtain

$$\begin{eqnarray}&&\langle 1s| r\cos (\theta )| {{kp}}_{{\rm{fp}}}\rangle =-\sqrt{\displaystyle \frac{1}{3}}\langle 1s\parallel r\parallel {{kp}}_{{\rm{fp}}}\rangle ,\end{eqnarray} \tag{ 27 }$$

where the reduced matrix element $\langle 1s\parallel r\parallel {{kp}}_{{\rm{fp}}}\rangle$ involves radial integration only. In the free-particle basis, the reduced electric-dipole transition matrix element is

$$\begin{eqnarray}&&\langle 1s\parallel r\parallel {{kp}}_{{\rm{fp}}}\rangle =32\displaystyle \frac{{k}^{2}}{{(1+{k}^{2})}^{3}},\end{eqnarray} \tag{ 28 }$$

where $| 1s\rangle$ represents the ground state of hydrogen atom and $| {{kp}}_{{\rm{fp}}}\rangle$ is p-type free particle final state with wavenumber $k=\sqrt{2{E}_{{f}}}$. At the same time, one can apply equation (25) to calculate the matrix element (28) with functions $| m\rangle$ obtained via diagonalization of the free particle Hamiltonian ${\hat{H}}_{{\rm{fp}}}=-{{\rm{\nabla }}}^{2}/2$ using p-type Gaussian basis functions

$$\begin{eqnarray}&&\langle 1s\parallel r\parallel {{kp}}_{{\rm{fp}}}\rangle =\displaystyle \sum _{m}\langle 1s\parallel r\parallel m\rangle \langle m| {{kp}}_{{\rm{fp}}}\rangle ,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\langle r| m\rangle =\displaystyle \sum _{a=1}^{N}{C}_{{ma}}\displaystyle \frac{{2}^{\tfrac{11}{4}}{\alpha }_{a}^{\tfrac{5}{4}}}{\sqrt{3}{\pi }^{\tfrac{1}{4}}}r\exp [-{\alpha }_{a}{r}^{2}],\end{eqnarray} \tag{ 30 }$$

where ${k}_{m}=\sqrt{2{E}_{m}}$, and the expansion coefficients, C_ma, are obtained by solving the Roothaan equations, (2), in the p-type Gaussian basis set. A total of N = 30 basis set exponents with ${\alpha }_{a}=\alpha {\beta }^{a-1}$, $\alpha =0.01$, and $\beta =1.5$ were used; Convergence and linear dependence properties of this basis set are discussed in [50]. Note that in the case of a free particle all energies E_m are positive in the final state. Since we know the exact solution for the final state wavefunction it is possible to calculate the sampling functions $\langle m| {{kp}}_{{\rm{fp}}}\rangle$, using (30)

$$\begin{eqnarray}\langle m| {{kp}}_{{\rm{fp}}}\rangle & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{r}\langle m| r\rangle r\sqrt{2\pi {kr}}{J}_{3/2}({kr}),\\ & = & \displaystyle \sum _{a=1}^{N}{C}_{{ma}}{(8\pi )}^{\tfrac{1}{4}}\displaystyle \frac{{k}^{2}}{\sqrt{3}{\alpha }_{a}^{\tfrac{5}{4}}}{{\rm{e}}}^{-\tfrac{{k}^{2}}{4{\alpha }_{a}}}.\end{eqnarray} \tag{ 31 }$$

Some of the sampling functions (31) are shown in figure 2. For the sake of brevity we use a simplified notation, so that $\langle m| {{kp}}_{{\rm{fp}}}\rangle \equiv \langle m| k\rangle$. Substitution of sampling functions (31) in (29) leads to the representation of the matrix elements using discrete positive energy solutions. The results are presented in figure 3. Note that the exact free particle continuum states were used to calculate sampling functions. As one can see the agreement between matrix elements calculated using equations (28) and (29) is quite good up to $k\approx 10$ a.u., which corresponds to energies of about ${E}_{{f}}=600\,\mathrm{eV}$. For higher values of k, the interpolation procedure is unsatisfactory and inaccurate. This reveals a fundamental limitation about the particular basis set, rather than the use of L² amplitudes to extract information about continuum processes. In order to expand a rapidly oscillating continuum function using pseudo-states (30) one needs a very dense set of ${\alpha }_{a}$ corresponding to $\beta \to 1$, but that leads to a well known problem of the linear dependence of the basis [50]. While this problem may be circumvented entirely using the Laguerre–Sturmian basis, the Gaussian basis is intrinsically limited in this regard.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The method proposed in the previous section is now applied to the calculation of the matrix element defined by (28). In equation (29) we take the sampling functions $\langle m| {{kp}}_{{\rm{fp}}}\rangle \equiv \langle m| k\rangle$ given by equation (26). Figure 1 represents the smooth interpolation, s(k), of the set of points $\{{k}_{n},n\}$ so that $s({k}_{n})=n$. Function y(k) interpolates the set of points $\{{k}_{n},{y}_{n}\}$, where y_n is nth zero of the Bessel function of the order ν, so that $y({k}_{n})={y}_{n}$. Figure 4 compares the Bessel–Hankel sampling functions with m = 1 and parameter $\nu =0.5,1.5,2.5$ with the exact $\langle 1| k\rangle$ given by (31). One can see that those functions almost coincide exactly for $\nu =1.5$, while for $\nu =0.5,2.5$ there are small deviations. It is clearly seen that the zeros of the exact function (30) only approximately coincide with eigenvalues k_n, while those of Bessel–Hankel sampling functions coincide exactly. Another important difference is that the number of zeros for the exact sampling function (28) equals to the number of eigenvalues N, while Bessel–Hankel sampling functions have an infinite number of zeroes. We expect, therefore, both functions to coincide in the limit of a complete basis set with a dense distribution of eigenvalues and for the fitting to be accurate for a finite-dimensional representation except in the high-energy limit.

Zoom In Zoom Out Reset image size

Figure 5 presents the result of interpolation of the reduced electric dipole matrix element $\langle 1s\parallel r\parallel {{kp}}_{{\rm{fp}}}\rangle$ using (29) with Bessel–Hankel sampling functions. The fitting deviates from the exact solution as the wavenumber k increases, but replicates the behaviour of the exact fitting (29) to a very good approximation, as shown in figure 3. The best results are obtained for $\nu =1.5$, which has the same asymptotic behaviour as the exact p-type free particle wavefunction. However for s-type ($\nu =1/2$) and d-type ($\nu =5/2$) fitting functions give a fitting that is essentially identical to that obtained with $\nu =1.5$. In fact the difference between s,d, and p-type fittings is noticeable only for small kinetic energies, where the sampling functions are proportional to ${k}^{(\nu +1/2)}$, and can not be distinguished at the present scale.

Zoom In Zoom Out Reset image size

In order to demonstrate the efficiency of the proposed method for calculation of photoionization matrix elements, we consider next the hydrogen atom. The exact expression for the reduced electron dipole matrix element is given by [35]

$$\begin{eqnarray}&&\langle 1s\parallel r\parallel {{kp}}_{H}{\rangle }^{2}={2}^{8}\displaystyle \frac{2\pi k}{{(1+{k}^{2})}^{5}}\displaystyle \frac{{{\rm{e}}}^{-\tfrac{4\arctan (k)}{k}}}{1-{{\rm{e}}}^{-2\pi /k}}.\end{eqnarray} \tag{ 32 }$$

To obtain the hydrogen continuum pseudo-states we used p-type even tempered Gaussian basis set with N = 35 basis set functions with ${\alpha }_{a}=\alpha {\beta }^{a-1}$, $\alpha =0.001$, and $\beta =1.5$. The six lowest-energy p-type states obtained using this basis set have negative energies and correspond to bound states. Therefore unlike the example described in the previous section, where all the states had positive energies, the negative-energy states are excluded from the summation in (29). This leads to the expression for wavepacket representation of hydrogen photoionization matrix elements given by

$$\begin{eqnarray}&&\langle 1s\parallel r\parallel {{kp}}_{H}\rangle =\displaystyle \sum _{(m=1| {E}_{m}\gt 0)}^{N}\langle 1s\parallel r\parallel m\rangle \langle m| k\rangle ,\end{eqnarray} \tag{ 33 }$$

where the sampling functions $\langle m| k\rangle$ are given by equation (26) with $\nu =1/2$, and m = 1 denotes the lowest positive energy pseudo-state. The result for hydrogen is shown in figure 6. As in the previous example of free particle final state, one can notice rapidly accumulating error of the interpolation for values of $k\approx 10$ a.u. due to the reasons outlined in previous section. Also some new deviations of interpolation from an exact solution (32) can be observed for interpolation on the range of eigenvalues ${k}_{1}\lt k\lt {k}_{4}$. Deviation in that region was barely noticeable in figure 5 since the sampling function, especially with $\nu =3/2$ replicated the behaviour of exact functions $\langle m| k\rangle$ to a very good approximation. Kramer theorem [46] indicates that if one would have a large number of sampling points $\langle 1s\parallel r\parallel m\rangle$ in the limit of small kinetic energy, the deviation in that region would vanish; but for finite basis sets one does not enjoy the same close match between the sampling scheme and the asymptotic behaviour of the matrix elements. For a small number of Gaussian basis functions we always have only a few pseudo-states with small eigenvalues for which ${k}_{m}\gg 1$, so the interpolation experiences some undersampling. One could expect the same thing to happen at the other end of the energy range, but Gaussian solutions tend to fail in describing rapid oscillation of continuum states $| k\rangle$ before that happens. We also performed interpolation with sampling functions (22) with $\nu =3/2,5/2$, but it did not increase the quality of the interpolation. This occurs because Bessel–Hankel sampling functions with $\nu =3/2$ do not have the same asymptotic behaviour near the origin as the actual sampling functions $\langle m| k\rangle$ for hydrogen (see figure 4(b)).

Zoom In Zoom Out Reset image size

For more complex systems with many electrons additional terms due to electron–electron interaction arise in the Hamiltonian of the system. For such systems no analytic solutions exists. Discrete states can be calculated to a very good precision by diagonalizing the Hamiltonian in a sufficiently large basis set. While the few lowest negative energy solutions are usually associated with bound state orbitals, the positive energy states are used mainly in the finite summations over virtual states that arise in perturbation theory. Using methods developed in this work we can extract some information, such as photoionization matrix elements, from these positive energy solutions, which we shall show in the case of the sodium atom. We use the simplest Hartree–Fock method for calculations of its ground state electronic structure and employ the frozen core approximation for calculations of the valence electron states (see for example [53]). In this approximation the states of the valence electron are calculated in the Hartree–Fock potential of the closed-shell electrons of sodium ion. Continuum pseudo-states were calculated with the same primitive Gaussian basis set as was used for hydrogen atom. Approximately 3% accuracy in the energy of 3s orbital can be achieved using this method, although more sophisticated and precise approaches give considerably better precision [51, 54].

The result of interpolation on the photoionization matrix element in the length form of the dipole operator for sodium is shown in figure 7. One can see that the fitted curve follows the finite difference solution of Hartree–Fock equations [51] very accurately up to the lowest energy eigenvalue. Figure 8 represents photoionization cross-section obtained by substituting interpolated matrix element into expression (1). Although our solution accurately matches length gauge Hartree–Fock results of [51], it deviates considerably from experimental results from [52] given by the blue points. The reasons for such a discrepancy are not entirely clear and [51] discusses some possibilities in more detail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have demonstrated how the Kramer sampling theorem can be used to recover photoionization matrix elements from conventional L² discretization of electronic spectra, such as one obtains from the positive energy solutions of the Roothaan equations. Although a Gaussian basis set does not give an appropriate continuum wavefunctions, it still can be used to recover the matrix elements involving continuum final state. This procedure has being demonstrated and discussed for hydrogen and sodium but can be easily applied to more complex atoms and molecules. The main advantage of proposed method over the spline interpolation is that it is globally defined withing the eigenvalue range of up to few hundred eV. The proposed method has a close connection with Shannon–Nyquist sampling theorem widely used in digital signal sampling and processing. Therefore some important refining techniques that are widely used in signal processing can be adopted to refine the evaluation of continuum processes in complex quantum systems.

Here, we have solely considered atomic radiative transitions, but the Gaussian basis set methods of quantum chemistry have been developed for molecular applications. Decomposition of molecular states into their atomic constituents is always possible, so that the atomic sampling methods developed here may be applied to each atomic component and convolved to produce molecular information by superimposing quantum mechanical amplitudes. The implementation of these ideas is the subject of ongoing research that will be reported in future publications.

The authors acknowledge the support of the Australian Research Council through the Centre of Excellence for Coherent x-ray Science and Advanced Molecular Imaging. AK is grateful to Daniel Lewis for a useful discussion during the work on this manuscript.