Atomic photoionization cross-sections beyond the electric dipole approximation

A methodology is developed to compute photoionization cross-sections beyond the electric dipole (BED) approximation from response theory, using Gaussian type orbitals and plane waves for the initial and final states, respectively. The methodology is applied to compute photoionization cross-sections of atoms and ions from the first four rows of the Periodic Table. Analysing the error due to the plane wave description of the photoelectron, we find kinetic energy and concomitant photon energy thresholds above which the plane wave approximation becomes applicable. The correction introduced by going beyond the electric dipole approximation increases with photon energy and depends on the spatial extension of the initial state. In general, the corrections are below 10% for most elements, at a photon energy reaching up to 12 keV.

A methodology is developed to compute photoionization cross-sections beyond the electric dipole (BED) approximation from response theory, using Gaussian type orbitals and plane waves for the initial and final states, respectively. The methodology is applied to compute photoionization cross-sections of atoms and ions from the first four rows of the Periodic Table. Analysing the error due to the plane wave description of the photoelectron, we find kinetic energy and concomitant photon energy thresholds above which the plane wave approximation becomes applicable. The correction introduced by going beyond the electric dipole approximation increases with photon energy and depends on the spatial extension of the initial state. In general, the corrections are below 10% for most elements, at a photon energy reaching up to 12 keV. a) iubr@kth.se

I. INTRODUCTION
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 with calculated electronic structures or total densities of states (DOS) is used to interpret the quality of a particular electronic structure method. [5][6][7][8][9][10][11][12][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 done by multiplying the projected density of states with the corresponding photoionization crosssection calculated for atomic orbitals within the dipole approximation, as performed, for example, in Refs. [14][15][16][17][18]. Another possibility is to compute the dipole transition matrix elements using the Fermi golden rule, as in Refs. [19][20][21][22][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 1,500 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. 28 Other attempts have also been made using plane waves, 21,22,[29][30][31] orthogonalized plane waves, 22,[29][30][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 4fmaterials. 20,[37][38][39] The electric dipole 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. Nondipole effects have been theoretically described and experimentally observed especially in terms of the angular distributions of the emitted photoelectrons. For overviews from various perspectives, we refer the reader to Refs. [23,[40][41][42][43][44][45][46] and articles therein.
In the present study, we analyse the correction introduced by using the full field operator in the expression for the atomic photoionization cross-section. For this purpose, we follow the theoretical derivation based on response theory from Ref. [47] to obtain an equation for the cross-section beyond electric dipole approximation, denoted by σ BED in the present work.
By using Gaussian type orbitals (GTOs) and plane waves (PWs) for the initial and final states, respectively, we compute BED cross-sections for the atoms in the first four rows of the Periodic Table, alongside their most common ions, and estimate the magnitude of the BED correction as a function of photon energy for each system. Furthermore, we analyse the error introduced by the PW description for the outgoing photoelectron and find a kinetic energy threshold above which this approximation becomes applicable.
The article is organized as follows: the full derivation of the expression for σ BED is presented in Section II, the computational details are described in Section III, while results for selected atoms are presented in Section IV. The complete set of data for all atoms and ions included in this study is provided at https://static.sys.kth.se/cbh/teochem/.

II. THEORY
In this section we present a detailed derivation of the equations used to calculate atomic and molecular photoionization cross-sections beyond the electric dipole approximation.
First, we define the electronic Hamiltonian which describes the interaction between the atomic or molecular system with a classical electromagnetic wave. Second, we derive a general equation for the absorption cross-section and, finally, we obtain an expression for the photoionization cross-section by approximating the wave function corresponding to the photoelectron by a plane wave (PW) and by describing the atomic orbitals as a linear combination of Gaussian type orbitals (GTOs).

A. The Interaction Hamiltonian
The interaction of a system of electrons within an atom or molecule with a time-dependent electromagnetic field, assuming the Born-Oppenheimer approximation, can be described by the non-relativistic Hamiltonian: [48][49][50] where N is the number of electrons, m e is the electron mass, e is the elementary charge,p j = −ih∇ j is the momentum operator,ŝ j is the spin operator, A(r j , t) is the time-dependent vector potential, φ(r j , t) is the time-dependent scalar field and B(r j , t) is the time-dependent magnetic field.
In the Coulomb gauge, for radiation without electrostatic sources, the vector potential is chosen divergence free and the constant scalar potential is imposed to vanish (i.e. ∇ · A=0 and φ(r, t) = 0). 50 Additionally, the Zeeman term is small compared to the others and may be neglected. By also neglecting the term proportional to A 2 (r j , t) as it does not contribute to one-photon processes, 51 the equation becomes: A(r j , t) ·p j .
In the Coulomb gauge, for a classical monochromatic electromagnetic wave with angular frequency ω, the vector potential, the electric field, and respectively the magnetic field are: 47 where is the unit vector of the polarization, k is the wavevector in the direction of propagation, A ω 1 (r j ) ≡ A ω 1 j implicitly contains the dependence on r j , and A 0 = E 0 /2iω in terms of the real amplitude of the electric field, E 0 . We note that the intensity of the incoming wave may be calculated from the time-averaged magnitude of the Poynting vector over one period, 52 i.e. I(ω) = ε 0 c E 2 0 /2, with ε 0 the vacuum permittivity. Finally, the interaction potential of Eq. (2) may be written as: where we have denoted:ν Using Eq. (6), the Hamiltonian takes the form:

B. Absorption Cross-Sections from Response Theory
The problem presented in the previous section may now be addressed using response theory, as discussed in detail in Ref. [53]. Briefly, response theory consists in describing the time evolution of the particle system |0(t) using a unitary transformation of the timeindependent ground state |0 , where the time-dependent unitary operatorP (t) is expanded in series in terms of the perturbationV (t). The evolution of the average value of a generic operatorÔ can be calculated using |0(t) . Specific response functions may be obtained by truncating the series at the desired order. A detailed derivation is provided in the Supplementary Material. Given that the process we aim to describe is photoionization due to single photon absorption, we focus on the linear response function only, which for a generic operatorÔ is: Here, |n is an excited state of the time-independent HamiltonianĤ 0 with corresponding eigenvalue E n , ω n0 = (E n − E 0 )/h, and γ n is the half-width broadening associated with the inverse lifetime of the excited state.
With the purpose of deriving an equation for the absorption cross-section, we express the rate of absorption within a volume V in terms of the work done by the electric field E(r, t) by generating a current density j(r, t) in the sample (in our case, limited to bound electrons): 52 The time average of the absorption rate over one period may be expanded as: 47,54 where α, β, γ are the absorption coefficients and σ(ω) = α(ω) is the linear absorption crosssection.
To obtain an expression for σ(ω) we use the non-relativistic current density operator: 55 where δ(r j − r) is the Dirac delta function and we have neglected the spin contribution and diamagnetic current.
By replacing the generic operatorÔ in Eq. (9) withĵ(r), we obtain the first order perturbation-induced current density: The first order transition rate may now be calculated using the first order expression of j from Eq. (13) and the expression for E(r, t) from Eq. (4): where, in the first step, the order of integration between sample and electronic degrees of freedom is interchanged and, in the second step, the Dirac delta function has been integrated out.
Considering that in the Coulomb gauge the vector potential is divergence free, Eq. (14) can be written where we identified the linear response function: 47 The absorption rate can now be written in terms of the linear response function: Finally, the time-averaging over one period of time, T , is carried out to arrive at dW dt where only the terms with ω 1 = −ω 2 of the sum integrate to non-zero values (equal to T ).
Comparing Eqs. (11) and (18), the linear absorption cross-section can be expressed as: With G(−ω; ω) from Eq. (16), we obtain: The first term corresponds to photon absorption, while the second term corresponds to stimulated emission. We keep only the absorption term and calculate the absorption crosssection in the limit γ n → 0: where we have used: Note that the dimension of δ(ω n0 − ω) is time.

C. Approximations for the Initial and Final States
Eq. (23) 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·r j · ∇ j can be re-written in terms of spin-orbitals: 56 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 BED 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. (24) becomes an integral over a continuum: where the kinetic energy of the photoelectron (E f =hω 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.
Eq. (29) 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 = xê x + yê y + zê z , and K = K xêx + K yêy + K zê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. (31) can be computed analytically. 57 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 BED photoionization cross-sections are treated with the same formalism and no further variables are introduced in the comparison.

III. COMPUTATIONAL DETAILS
The electronic structure of the atoms from the first four rows of the Periodic Table (H-Kr) was calculated with the Gaussian 16 58 quantum chemistry software, at the B3LYP level of theory 59 and using the universal Gaussian basis set (UGBS) by de Castro and co-workers. [60][61][62][63][64][65][66][67][68] The electronic structure calculations were performed using Cartesian Gaussian functions be- 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 the electronic configuration relaxed with B3LYP. The relaxed electronic structures were used to compute the dipole and beyond dipole photoionization cross-sections for a range of photon energies between 20 eV to 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. (29) and the computational set-up shown in Figure 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 energyhω. The norm of the photoelectron wavevector k e is where ε AO represents the computed eigenvalue of a particular atomic orbital.
The BED photoionization cross-section is finally computed by averaging over all possible photoelectron directions generated on a sphere using the Lebedev quadrature. 83 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 BED, the order 50 was instead required to fully converge the cross-sections. We have therefore computed both BED and dipole cross-sections using the 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 p α x means that the spin-up 2p x orbital is occupied. The total atomic photoionization cross-section was calculated by summing over the crosssections 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.

IV. RESULTS AND DISCUSSION
In the following, (Section A) we analyse 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. (Section B) We then compare the dipole and BED 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 to 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. 2a), C (Fig. 2b), Mn (Fig.   2c), and O (Fig. 2d). The figure also depicts the ratio σ PW /σ YL as a function of photon energy. ratio between σ P W and σ Y L is also depicted (dotted line). The labels 1s, 2s, and 2p mark the ionization thresholds for these orbitals. The PW cross-sections have been multiplied by a factor of 1/3 to average over the directions of the incoming photon, as performed in Refs. [24,25]. Note that the vertical axis is logarithmic.

A. Performance of the PW Approximation
What becomes clear from the figure is that σ PW converges towards σ 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 approximatively 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. to average over the directions of the incoming photon as performed in Refs. [24,25]. Note that the vertical axis is logarithmic. The inset shows the product between the radial parts of a 2s atomic orbital with a PW as a function of the distance to the nucleus (r).
In addition, we have also analysed 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 Figure 3 for C (Fig. 3a), and O (Fig. 3b) 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, athω ≈ 90 eV andhω ≈ 160 eV, respectively (marked by arrows in Figure 3). The reason for this dip 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. To further elucidate if the PW approximation for the final state is applicable, we have compared the beyond electric dipole cross-section (σ BED ) 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. 4a), N (Fig. 4b), Ar (Fig. 4c), and O (Fig. 4d). The experimental data was compiled for nine atoms from selected reference sources by Berkowitz. 85 The measurements are typically performed in gas phase and at room temperature, as described in Refs. [86][87][88]. The calculated data was obtained within the dipole approximation using the Hartree-Fock-Slater method (the same algorithm used by Yeh and Lindau) and is stored on the National Institute of Standards and Technologies (NIST) database. 89 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 σ BED are very close to the calculated dipole cross-sections σ NIST . As it will be discussed in the following section, the reason for this similarity is that the relative correction introduced by BED to the total photoionization cross-section does not amount to more than ∼5%. (stored on the NIFS database, 90 σ exp , blue squares), and with total dipole cross-sections calculated using the Hartree-Fock-Slater method (σ NIST , red circles) from the NIST database. 89 The comparison is shown for (a) Li, (b) N, (c) Ar, and (d) O. σ BED has been multiplied by a factor of 1/3 to average over the directions of the incoming photon, as performed in Refs. [24,25]. Note that the vertical axis is logarithmic.
We have also computed the photoelectron angular distribution at a photon energyhω = 5206 eV, for the Ar 1s orbital, shown in Figure 5. 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 Figure 5b. 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 ϕ. If the photoionization process would take  Fig. 5a 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: 91,92 where Ω is the solid angle, β is the anisotropy parameter, δ and γ are the asymmetry parameters which quantify the nondipolar effects, and P 2 (cos θ) is the second Legendre polynomial.
At the magic angles θ = 54.7 • and 180 • −θ, P 2 (cos θ) is zero and the differential cross-section then depends only on the azimuthal angle ϕ: 93 We have calculated σ BED for the 1s atomic orbital of Ar at different ϕ angles by rotating the k e vector around the z-axis, as depicted in Figure 5b. 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. 91 The experimental and calculated curves match rather well, and our calculated angle-dependent σ BED 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 number. Using the PW final state, the cross-sections computed for orbitals with one or more nodes have discontinuities at photon energies (generally close to the ionization threshold) where the positive and negative contributions to the transition matrix element cancel out.
These are inherent to the GTOs used for the initial state. Finally, our calculated σ BED is able to reproduce the nondipolar angular distribution of photoelectrons at high photon energies.

B. Magnitude of the BED Correction
Having analysed 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 6 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 BED 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 energyhω = 12 keV. This is not unexpected given that only 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 a 2s orbitals). The total cross-sections (not shown) follow closely the behaviour of the s-type orbitals. This is also expected, since the last shell to ionize will be the one which will dominate, as illustrated in Figure 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 BED are slightly larger, between 5-10% for the three atoms shown in Fig. 6. One last aspect to note is the appearance of sharp dips, generally at low photon energies, notably visible in the case of the Mn atom (marked with arrows in Fig. 6c). These are related to the discontinuities 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. a For the Ni and Co atoms, the B3LYP electronic configuration relaxed to a situation where the 4s orbital remained unoccupied irrespective of the initial guess. We have, therefore, used the calculated data for the ions instead.
Finally, we collected the corrections introduced by going beyond the electric dipole approximation for all atoms at the largest photon energy included in the study, i.e.hω = 12 keV. The results are represented as a function of atomic number in Figure 7. The corrections to the total cross-section are depicted in gray (Fig 7b), while the l-shells with the largest correction to the cross-sections are depicted in color (Fig. 7a). The dots are red if the maximum correction was obtained for an s orbital, blue for a p orbital, and brown for a d orbital. 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 Figure 7a 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 BED correction may be understood by visualizing the effect of the BED correction in the computational set-up, as represented in Figure 7c. 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 Figure 7a is that each edge is followed by a linear decrease of the magnitude of the BED correction. 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 Figure 7d.

V. CONCLUSIONS
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 experimental data and to 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 number.
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 BED 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
The