Atomic photoionization cross sections beyond the electric dipole approximation

Photoionization is the physical process in which atoms, molecules, or solids emit electrons upon irradiation. It constitutes the basis for photoelectron spectroscopy (PES), also known as photoemission, an experimental technique routinely used for the characterization of materials.^1,2 By recording the kinetic energy of emitted electrons, PES provides a description of the occupied electronic structure of the sample. Besides information on the binding energies of the electronic levels, PE spectra measured using different photon energies can be used to distinguish between states with different atomic compositions.^3,4

Commonly, the comparison between experimental PE spectra and calculated electronic structures or total densities of states (DOS) is used to interpret the quality of a particular electronic structure method.^5–13 A step beyond such a simple comparison is to include in the calculation the probability of ionization of a particular state. This can be performed by multiplying the projected density of states with the corresponding photoionization cross section calculated for atomic orbitals within the dipole approximation, as performed, for example, in Refs. 14–18. Another possibility is to compute the dipole transition matrix elements using the Fermi golden rule, as in Refs. 19–23.

Photoionization cross sections have been calculated within the dipole approximation for different ranges of photon energies and using different mathematical descriptions for the initial and final states. Most notably, Yeh and Lindau computed the total and partial ionization cross sections for all atoms in the periodic table for photon energies between 10 and 1500 eV.^24,25 They used an algorithm introduced by Cooper and Manson,^26,27 where both initial and final states were described on a radial grid and the final state was relaxed in the presence of the positive ion.²⁸ Other attempts have also been made using plane waves,^{21,22,29–31} orthogonalized plane waves,^22,29–31 Coulomb waves,^21,32 or B-splines^33,34 to describe the final states. In the solid state community, time-reversed (TR) low energy electron diffraction (LEED) states^20,35,36 for the final state are used in a Green’s function method approach to obtain angle-resolved photoemission spectra (ARPES) of 3d- and 4f-materials.^20,37–39

The electric dipole (ED) approximation is based on the assumption that the electromagnetic field can be considered uniform across the spatial extent of the initial state wave function. This allows for the truncation of the series expansion of the spatial part of the field at zeroth order, i.e., e^ik·r ≈ 1. Here, k represents the wavevector of the electromagnetic radiation with wavelength λ, and we have k = 2π/λ. Thus, the dipole approximation is not valid when the wavelength λ is comparable to or smaller than the spatial extent of the initial-state wave function (or orbital), i.e., when the term k · r becomes comparable to or larger than 1. Going beyond, the electric dipole approximation becomes especially important when the ionizing photon energy is large and the ionized electronic states are highly delocalized, such as the valence orbitals of conjugated organic molecules or the valence bands of solids.

Non-dipole (ND) effects have been theoretically described and experimentally observed especially in terms of the angular distributions of the emitted photoelectrons.^23,40–47 In general, non-dipole photoionization cross sections and photoelectron angular distributions are calculated by including higher order terms, usually up to the first or second order, in the series expansion of e^ik·r as performed, for example, in Refs. 45 and 47–53. For a review on the topic, we refer the reader to Ref. 54. Not many studies include the full field operator in the photoionization cross section, but we note the work of Demekhin⁴¹ for atoms and of Seabra and co-workers²³ for molecules. Demekhin included the full field operator to calculate non-dipole cross sections for H, He, Be, Ne, Ar, and Kr using atomic orbitals relaxed in a scalar relativistic framework, while Seabra et al.²³ used Dyson orbitals and a plane wave (PW) final state to compute the error introduced by the dipole approximation for HF, H₂O, NH₃, and CH₄.

In the present study, we analyze the correction introduced by using the full field operator in the expression for the atomic photoionization cross section, following the theoretical derivation based on response theory from Ref. 55. By using Gaussian type orbitals (GTOs) and plane waves (PWs) for the initial and final states, respectively, we compute ND cross sections (denoted σ_ND throughout this article) for the atoms in the first four rows of the periodic table, alongside their most common ions, and estimate the magnitude of the ND correction as a function of photon energy for each system. Furthermore, we analyze the error introduced by the PW description for the outgoing photoelectron and find a kinetic energy threshold above which this approximation becomes applicable. It is often assumed that a PW description for the final state is correct at large photoelectron kinetic energies, but, to our knowledge, there is no systematic study of how the PW approximation performs for different kinetic energies, atoms, and atomic shells.

The article is organized as follows: the derivation of the expression for σ_ND is presented in Sec. II, the computational details are described in Sec. III, while results for selected atoms are presented in Sec. IV. The complete set of data for all atoms and ions included in this study is provided at https://www.theochem.kth.se/BED_atomic_cross_sections.

In this section, we present the derivation of the equation used to calculate atomic and molecular photoionization cross sections beyond the electric dipole approximation. We start from the expression of the absorption cross section derived from complex response theory by List and co-workers⁵⁵ and apply it to the photoionization process.

The photoabsorption cross section obtained as the linear response of the system of electrons to a time-dependent electromagnetic field is⁵⁵ 

where ℏ, e, ε₀, c, and m_e are the reduced Planck constant, elementary charge, vacuum permittivity, speed of light in vacuum, and electron mass, respectively. ω is the angular frequency of the electromagnetic field, k is the associated wavevector, and r_j is the position vector corresponding to electron j. $|0$ and $|n$ represent the many-body ground state and excited state, respectively. The summation is performed over all excited states and ω_n0 designates the angular frequency associated with the energy difference between the excited state and the ground state.

Equation (1) has been derived from response theory by List and co-workers,⁵⁵ but for the sake of completeness, we include the intermediate steps required to arrive at this expression in the supplementary material.

Equation (1) is a general expression applicable to single photon absorption where the initial and final states are described by the many-body wave functions |0⟩ and |n⟩, respectively. Considering that the electronic structure of the atoms we are interested in may be approximated by a single Slater determinant, the matrix element of the one-electron operator $e−ik⋅rjϵ⋅∇j$ can be re-written in terms of spin-orbitals,⁵⁶ 

where χ_i is an occupied atomic orbital and χ_f is the wave function corresponding to the photoelectron. We note that the same expression is also valid in the case of molecular orbitals.

Considering that ND corrections become more important at high photon energies, where the wavelength of the ionizing electromagnetic radiation becomes comparable to the spatial extent of (valence) atomic orbitals, it is reasonable to assume that the photoelectron, with a high kinetic energy, can be described as a periodic plane wave normalized over an arbitrary volume Ω,

with density of states

Using PWs and the corresponding free electron density of states, in the limit of a large volume Ω, the sum over final states in Eq. (2) becomes an integral over a continuum,

where the kinetic energy of the photoelectron (E_f = ℏω_f) satisfies energy conservation ω_f = ω_i + ω and the photoelectron wavevector k_e is correspondingly chosen. The advantage of using a PW for the final state is the fact that PWs are eigenfunctions of the momentum operator, and the expression for the photoionization cross section may be, therefore, further simplified,

Finally, by expressing the initial state wave function χ_i as a linear combination of Gaussian type orbitals (GTOs), the photoionization cross section becomes

where the summation is over all primitives used to describe the initial state and g_l is a GTO with l = l_x + l_y + l_z,

with R_x, R_y, and R_z the coordinates of the atom where the GTO is centred.

Equation (7) can be straightforwardly generalized to molecular orbitals, by replacing χ_i with a linear combination of atomic orbitals.

Finally, the matrix element to be calculated is

where $r=xe^x+ye^y+ze^z$ and $K=Kxe^x+Kye^y+Kze^z$ determined from the input parameters k_e and k, K = k_e − k.

The advantage of working with PWs and GTOs is the fact that the integrals in Eq. (9) can be computed analytically.⁵⁷ We provide analytic formulas up to l = 5 in the supplementary material. Furthermore, the dipole approximation is obtained simply by setting the photon wavevector k = 0, and hence, the dipole and ND photoionization cross sections are treated with the same formalism, and no further variables are introduced in the comparison.

The electronic structure of the atoms from the first four rows of the periodic table (H–Kr) was calculated with the Gaussian 16⁵⁸ quantum chemistry software, at the B3LYP level of theory⁵⁹ and using the universal Gaussian basis set (UGBS) by de Castro and co-workers.^60–68 The electronic structure calculations were performed using Cartesian Gaussian functions because the expression for the matrix elements involved in the photoionization cross section, Eq. (7), has been derived for this GTO variant. For the N atom, we have additionally computed cross sections using five other basis sets, namely, STO-3G,^69,70 4-31G,^71–74 6-31G,^71–80 cc-pVTZ,⁸¹ and cc-pVQZ.⁸² The basis set results are shown in Fig. S1 of the supplementary material.

In addition to the neutral atoms, we have also computed the most common ions corresponding to each atom. The atoms and ions included are listed in Table I, alongside their electronic configurations. The B3LYP-optimized electronic structures were used to compute the dipole and beyond dipole photoionization cross sections for a range of photon energies between 20 eV and 12 keV, with an energy step of 10 eV. A photoionization cross section was calculated for each atomic orbital in the relaxed electronic configuration using Eq. (7) and the computational setup shown in Fig. 1. The wavevector of the photon, k, was set in the y-direction, while the polarization vector was chosen in the z-direction. The norm of k is calculated from the photon energy ℏω. The norm of the photoelectron wavevector k_e is calculated from the kinetic energy KE, in turn, computed from energy conservation,

where ε_AO represents the computed negative eigenvalue of a particular atomic orbital.

The ND photoionization cross section is finally computed by averaging over all possible photoelectron directions generated on a sphere using the Lebedev quadrature.⁸³ We tested the Lebedev quadrature at different orders. Due to the spherical symmetry of the atomic system, the lowest order of the quadrature (i.e., 6) sufficed in the case of the dipole photoionization cross section. For ND, the order 50 was instead required to fully converge the cross sections. We have therefore computed both ND and dipole cross sections using the order 50 for the Lebedev quadrature.

In the case of open-shell atoms, we additionally averaged over the possible initial state configurations of the partially filled shell. For example, in the case of the C atom, the 2p cross section is computed as follows:

where the index $pxα$ means that the spin-up 2p_x orbital is occupied.

The total atomic photoionization cross section was calculated by summing over the cross sections of the occupied atomic shells. Dipole photoionization cross sections were calculated by the same algorithm with the only difference that the wavevector k was set to zero.

In the following (Sec. IV A), we analyze the performance of the PW approximation for the final state by comparing our calculated photoionization cross sections to calculated results from the literature and experimentally measured data (Sec. IV B). We then compare the dipole and ND cross sections to each other and calculate the magnitude of the relative correction introduced by going beyond the electric dipole approximation as a function of photon energy.

Figure 2 shows the total dipole photoionization cross section calculated using a PW final state (σ_PW) in comparison with the cross section calculated by Yeh and Lindau^24,25,84 using a final state relaxed in the presence of the positive ion (σ_YL). The comparison is shown as a function of photon energy, from 20 eV up to 1500 eV, for Li [Fig. 2(a)], C [Fig. 2(b)], Mn [Fig. 2(c)], and O [Fig. 2(d)]. The figure also depicts the ratio σ_PW/σ_YL as a function of photon energy.

What becomes clear from the figure is that σ_PW converges toward σ_YL as the photon energy is increased, i.e., as the kinetic energy of the photoelectron increases. In the vicinity of the ionization threshold, the PW approximation performs poorly, but it improves as the photon energy is increased. The photon energy above which the PW approximation becomes applicable depends on the particular element, becoming larger as the atomic number increases. For example, going approximately 100 eV above the 1s threshold is required in the case of the Li atom, but more than 700 eV above the 2s threshold is instead required for the Mn atom. The poor performance in the vicinity of the ionization threshold is expected since the outgoing electron, in this case, has low kinetic energy and is not free electron like.

In addition, we have also analyzed selected partial dipole photoionization cross sections and how they compare to the results obtained by Yeh and Lindau.^24,25,84 These comparisons are shown in Fig. 3 for C [Fig. 3(a)] and O [Fig. 3(b)] in the photon energy window between 20 eV and 1500 eV. As it can be seen from Fig. 3, the PW approximation becomes better further away from the ionization threshold, but different atomic orbitals behave differently. The agreement to the Yeh and Lindau results^24,25,84 is better in the case of s orbitals than in the case of p orbitals.

It is interesting to note that the cross sections calculated for the 2s orbitals of both the C and O atoms present a pronounced dip (Cooper minimum), at ℏω ≈ 90 eV and ℏω ≈ 160 eV, respectively (marked by arrows in Fig. 3). The reason for this minimum is that the transition matrix element is zero at that particular photon energy. Because the 2s atomic orbital has one node, it is expressed as a linear combination of GTOs with both positive and negative coefficients. Since the integrals involving s-type GTOs and a PW are always real and positive (see the analytical formulas for the integrals in the supplementary material), there are both negative and positive contributions to the transition matrix element. For a particular value of the photoelectron wavevector, the positive and negative contributions cancel out, giving rise to the minimum in the photoionization cross section. All other orbitals which have nodes present similar minima. General equations for different types of atomic orbitals described as linear combinations of GTOs may be worked out. We include one example for a very simple 2s orbital in the supplementary material.

Minima arising from the cancellation of the positive and negative contributions to the integrand have been first described by Cooper,²⁶ using final states relaxed in the presence of the positive ion. Subsequently, Cooper minima have been observed experimentally, most notably in the case of Ar, due to the ionization of the 3s subshell.⁸⁵ 

We should note that a perfect cancellation between the positive and negative parts of the transition matrix element does not arise for Slater type orbitals (STOs) in combination with a PW final state, as exemplified in the supplementary material for a 2s STO.

To further elucidate if the PW approximation for the final state is applicable, we have compared the non-dipole cross section (σ_ND) to experimental data and calculations performed over a much larger photon energy range. Figure 4 shows this comparison for photon energies in the range of 20 eV to 12 keV for Li [Fig. 4(a)], N [Fig. 4(b)], Ar [Fig. 4(c)], and O [Fig. 4(d)]. The experimental data were compiled for nine atoms from selected reference sources by Berkowitz.⁸⁶ The measurements are typically performed in gas phase and at room temperature, as described in Refs. 89–91. The calculated data were obtained within the dipole approximation using the Hartree-Fock-Slater method (the same algorithm used by Yeh and Lindau) and are stored on the National Institute of Standards and Technologies (NIST) database.⁸⁸ 

The figure shows that far from the ionization threshold the cross sections calculated using a PW final state (beyond the electric dipole approximation) compare very well with the experimental data. Additionally, the values of σ_ND are very close to the calculated dipole cross sections σ_NIST. As it will be discussed in Sec. IV B, the reason for this similarity is that the relative correction introduced by ND to the total photoionization cross section does not amount to more than ∼5%.

Figure 5 shows the comparison between partial cross sections calculated beyond the electric dipole using a PW final state and the first order corrected photoionization cross sections computed by Cooper.⁴⁷ The two results are in good agreement in the case of the s orbitals, once far enough from the ionization threshold. In the case of the 4s orbital of Kr, the comparison is hindered by the presence of one of the Cooper minima in σ_ND at approximately 5 keV photon energy. The cross sections calculated using a PW final state and the full field operator are larger than the first order corrected cross sections for the p and d orbitals. This difference may partly be assigned to the larger error introduced by the PW approximation for these types of orbitals (see Fig. 3) and partly to the use of the full field operator in σ_ND vs. the first order correction in Ref. 47.^53,54

We have also computed the photoelectron angular distribution at a photon energy ℏω = 5206 eV, for the Ar 1s orbital, shown in Fig. 6. The angular distribution of photoelectrons is typically measured using linearly polarized light by placing an electron analyser at a fixed polar angle θ, as depicted in Fig. 6(b). The number of electrons emitted at different azimuthal angles is then recorded by rotating the electron analyser around the z-axis and performing measurements at different values of φ. We note that the polar angle θ is defined as the angle between the polarization direction of the photon (ϵ) and the photoelectron direction (k_e), while the azimuthal angle φ [shaded in gray in Fig. 6(b)] is the angle between the projection of k_e in the xy-plane and the photon direction (k).^47,92,93 If the photoionization process would take place according to the dipole approximation, the number of emitted photoelectrons would not depend on φ and the graph in Fig. 6(a) would be a straight horizontal line. The deviation of the curve from a straight line is therefore a measure of nondipolar asymmetry. The differential cross section for photoionization (dσ/dΩ) using linearly polarized light, taking into account quadrupole terms, may be expressed as^47,92

where Ω is the solid angle, β is the anisotropy parameter, δ and γ are the asymmetry parameters which quantify the nondipolar effects, and P₂(cos θ) is the second Legendre polynomial.

At the magic angles θ = 54.7° and 180° − θ, P₂(cos θ) is zero and the differential cross section then depends only on the azimuthal angle φ,⁹³ 

We have calculated σ_ND for the 1s atomic orbital of Ar at different φ angles by rotating the k_e vector around the z-axis, as depicted in Fig. 6(b). The photon polarization was kept in the z-direction and the wavevector of the photon was kept in the y-direction. We normalized the value of the calculated cross section to the value calculated for φ = 0°, as performed for the experimental data.⁹² The experimental and calculated curves match rather well, and our calculated angle-dependent σ_ND cross section is able to reproduce the magnitude of nondipolar effects observed experimentally.

To summarize this section, by comparing our atomic photoionization cross sections calculated using a PW final state to previously calculated results from the literature and experimental data, we can conclude that the PW approximation for the final state is applicable at high kinetic energies, i.e., far from the ionization threshold. However, the kinetic energy above which the approximation is applicable depends on the atom and type of orbitals since the comparison becomes worse for larger atoms and for orbitals with larger l quantum numbers. Using the PW final state, the cross sections computed for orbitals with one or more nodes present minima at photon energies (generally close to the ionization threshold) where the positive and negative contributions to the transition matrix element cancel out. Finally, our calculated σ_ND is able to reproduce the nondipolar angular distribution of photoelectrons at high photon energies.

Having analyzed the advantages and limitations of a PW final state, we now move to examine the magnitude of the correction introduced by going beyond the electric dipole approximation. Figure 7 shows the correction obtained for total and partial photoionization cross sections of selected atoms as a function of photon energy, between 20 eV and 12 keV. The correction is represented as

The partial cross sections are grouped together based on the orbital type (s, p, d) because the values computed from different shells (1s, 2s, etc.), but the same l quantum number, are very similar to each other.

As expected, the ND correction increases as the photon energy is increased. In the case of the s partial cross sections, the values reach only up to 5% for a photon energy ℏω = 12 keV. At this photon energy the wavelength of the photon (1.95 a.u.) begins to become comparable to the size of the orbitals (∼2 a.u., for 2s orbitals). The total cross sections (not shown) follow closely the behavior of the s-type orbitals. This is expected since the last shell to ionize will be the one which will dominate, as illustrated in Fig. 3, where the 1s cross section is 2-3 orders of magnitude larger than the cross sections of the 2s and 2p orbitals. For all the elements of the first four rows of the periodic table, at 12 keV, the last shell to ionize is either a 1s or a 2s shell, resulting in the high similarity between the correction for s partial photoionization cross section and the total one.

In the case of p and d atomic orbitals, the correction introduced by ND is slightly larger, between 5% and 10% for the three atoms shown in Fig. 7. One last aspect to note is the appearance of sharp peaks, generally at low photon energies, notably visible in the case of the Mn atom [marked with arrows in Fig. 7(c)]. These are related to the Cooper minima previously discussed, resulting from the cancellation between the positive and negative contributions to the transition matrix element which involves an orbital with at least one node.

We further compare our calculated ND corrections to the results of Demekhin.⁴¹ Figure 8 shows the comparison between the ratio σ_ED/σ_ND calculated using a PW final state and the same ratio obtained in Ref. 41, using atomic orbitals within a scalar relativistic framework. We obtain the same ratio as Demekhin⁴¹ for the s orbitals of H, He, and Be, and for the 2p orbitals of Ne. Our calculated ratio is larger than the one in Ref. 41 in the case of Kr 2p, Ar 1s > Ar 2s > Ar 2p, and Ne 1s > N 2s (where the atoms and subshells are ordered from the largest difference to the lowest). Considering this ordering, the difference between our results and those of Demekhin⁴¹ may be attributed partly to the fact that our electronic structure calculation does not include relativistic effects.

Finally, we collect the corrections introduced by going beyond the electric dipole approximation for all atoms at the largest photon energy included in the study, i.e., ℏω = 12 keV. The results are represented as a function of atomic number in Fig. 9. The corrections to the total cross section are depicted in gray [Fig. 9(b)], while the l-shells with the largest correction to the cross sections are depicted in color [Fig. 9(a)]. The atomic number at which a new shell starts to fill is marked by a dotted line. Note that the maximum correction may be obtained for a different orbital than the last to be filled. For example, the largest correction obtained for the Na atom corresponds to a photoionization of the 2p orbitals, rather than the 3s which is the last to be filled.

The first feature to notice in Fig. 9(a) is that there are three “edges”: the first at Z = 1, the second at Z = 5, and the third at Z = 21. Interestingly, these edges arise when making a transition from l to l + 1, i.e., going from s to p, and, respectively, from p to d, where the corrections to the p and d cross sections are larger.

The larger sensitivity of the p and d orbitals to the ND correction may be understood by visualizing the effect of the ND correction in the computational setup, as represented in Fig. 9(c). The dipole approximation is obtained by setting the photon wavevector to zero, meaning the transition matrix elements that have to be calculated involve the initial state (s, p, or d orbital) and a PW of wavevector k_e directed along the z-axis. Going beyond the dipole approximation (k ≠ 0) translates into having to calculate transition matrix elements between the initial state and a new PW now directed along the z′-axis, with a new K = k_e − k. This new PW has both a different magnitude and a different direction, compared to k_e. In the case of the s orbital, which is spherically symmetric, only the difference in magnitude affects the transition matrix element, while in the case of the p and d orbitals, both magnitude and direction affect the transition matrix element, making this type of orbitals more sensitive to the dipole approximation.

Another feature of the graph in Fig. 9(a) is that each edge is followed by a linear decrease in the magnitude of the ND correction. Considering that all these are valence orbitals of relatively similar binding energies (between 5 and 20 eV for the majority of elements) and that relativistic effects may be neglected, this decrease may be explained by the decrease in the spatial extent of the particular type of orbitals when going through the series. This is illustrated for the 2p elements in Fig. 9(d).

In summary, we have derived an equation to compute photoionization cross sections beyond the electric dipole approximation from response theory, using Gaussian type orbitals for the initial state and a plane wave with appropriate kinetic energy for the final state. We compared our calculated cross sections to the experimental data and to the calculated dipole photoionization cross sections from the literature. This comparison showed that the PW approximation for the final state is applicable far from the ionization threshold. However, how far is sufficient depends on the atom and the type of atomic orbital, the comparison becoming worse for the larger atoms and for orbitals with larger l quantum numbers.

Given that the dipole and beyond dipole photoionization cross sections can be computed analytically and on the same footing, we could estimate the relative correction introduced by going beyond the electric dipole approximation. As expected, this correction increases with photon energy, but it is, in general, rather small, ∼5%–10% at the largest photon energy we included (12 keV). The fact that we obtain such a small correction is not surprising, considering the small spatial extent of the atomic orbitals. Even though the ND corrections for the atoms are quite small, we expect larger corrections for more delocalized states as, for example, molecular orbitals of conjugated molecules or valence states of solids.

Finally, the description for the atomic photoionization cross sections may be improved by replacing the PW final state with a linear combination of PWs relaxed in the presence of the positive ion. This is, however, rather complicated because the relaxation should be performed for a small kinetic energy window centred around the kinetic energy obtained from energy conservation. Considering all states up to this kinetic energy is simply not feasible.

Supplementary material contains (a) the complete derivation of the linear response function corresponding to a time-dependent perturbation $V^$ and a generic operator Ô; (b) derivation of the equation for the photoabsorption cross section from complex response theory; (c) analytical formulas for the integrals in Eq. (9); (d) photoionization cross sections calculated using different basis sets for N; (e) the analytical form of a transition matrix element between a simplified 2s orbital (described as a linear combination of two GTOs) and a plane wave.

Financial support from the Knut and Alice Wallenberg Foundation (Grant No. KAW-2013.0020) and the Swedish Research Council (Grant No. 621-2014-4646) is acknowledged. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC and HPC2N. I.E.B. is very grateful to Faris Gel’mukhanov, Nanna Holmgaard List, Barbara Brena, and Igor Di Marco for all the insightful discussions. O.E. also acknowledges support from eSSENCE.