InternationalX-ray absorption spectroscopy and related techniquesTables for Crystallography Volume I Edited by C. T. Chantler, F. Boscherini and B. Bunker © International Union of Crystallography 2022 |
International Tables for Crystallography (2022). Vol. I. Early view chapter
https://doi.org/10.1107/S1574870722001598 ## Finite-difference method for the calculation of X-ray spectroscopies
All of the ingredients of the finite-difference method (FDM) for the calculation of X-ray near-edge absorption and emission as well as X-ray resonant diffraction peaks are presented. Solution of the relativistic Schrödinger equation is performed on a dense grid of points inside a sphere surrounding a cluster of atoms containing the absorber (or resonant) atom. The most important point is that this simple scheme allows a free description of the potential shape, which most often offers a good reproduction of the experimental spectra without additional parameters. It is shown how the center of the atoms, where a classical expansion in spherical waves takes place, is embedded in the grid. The continuity with the region of the outer sphere is performed in the same way and ensures convenient normalization of the states. Thanks to the appropriate algorithms, the calculation of all states of photoelectrons can now be performed using standard personal computers. The core–unoccupied state transitions are then calculated, finally giving the absorption or emission cross sections or the resonant diffracted beam intensities. Some examples with light and heavy chemical elements, both magnetic and nonmagnetic, demonstrate the possibility of FDM, which is at the heart of the Keywords: finite-difference method. |

The finite-difference method (FDM) is a general technique (Grossmann *et al.*, 2007) for solving partial differential equations. Its first formulation for solving the Schrödinger equation was given in the 1930s (Kimball & Shortley, 1934). FDM requires significant computing power and its progress has followed improvements in computer performance. One of its first satisfactory applications in the field of solid-state physics involves the calculation of binding energy states (Nieminen & Puska, 1983). We first developed it for scattering simulations of electrons (Joly, 1992) and positrons (Joly, 1994). FDM was applied for the first time to simulate X-ray absorption near-edge structure (XANES) at the end of the 1990s (Joly, 1997; Joly *et al.*, 1999) in the *FDMNES* code presented in Bunău *et al.* (2021). Its success came from the fact that it was the first satisfactory attempt to realize such simulations with a free shape potential, often resulting in a significant improvement of the calculation in term of accuracy. Other methods now allow calculation with full potentials, among others pseudo-potential techniques or the full linear augmented plane-wave method. Nevertheless, because of its simplicity and recent improvements in its computational efficiency, FDM remains widely used in the X-ray absorption spectroscopy (XAS) community.

Depending on the absorption edge, on the chemical species of the absorbing atom and to a lesser extent on the material itself, a density-functional theory (DFT) approach, which is theoretically valid only for ground states, can be sufficient to conveniently simulate an absorption spectrum. It is most often the case at the *K* and *L*_{1} edges of any elements and at the *L*_{2,3} edges of heavy elements. In this case, let us call 〈ϕ_{g}| and 〈ϕ_{f}| the mono-electronic core and photoelectron (or final) states. They are normalized (〈ϕ_{g}|ϕ_{g}〉 = 〈ϕ_{f}|ϕ_{f}〉 = 1) and the absorption cross section is given by where α is the fine-structure constant. , and ℏω are the core-state, final-state and incoming photon energies, respectively. is the density of state at . Note that the factor () and the extra term in the energy conservation, , give the effect of the intra-atomic screening of the core hole. For simplicity, as often performed, we omit both terms in the following. is the electromagnetic operator. Neglecting the magnetic interaction term, and with expansion up to the quadrupole electric term, it is given by where **r** is the position, **ɛ** is the polarization and **k** is the wavevector of the beam.

Harmonic representation is especially helpful to define the electronic states. We use Bessel functions, , and we know that are solutions of the Schrödinger equation when the potential *V* is zero, and that they form a complete basis. is the modulus of the electron wavevector in atomic units. *r* and are the radial and angular positions. The indices *f* are now defined by the quantum numbers *L*_{f} = (*l*_{f}, *m*_{f}) and the energy . The factor (κ/π)^{1/2} is such that the scalar product of the eigenstates gives the vacuum density of state in Rydberg^{−1}.

The asymptotic completeness theorem (Deift & Simon, 1976) states that a set of eigenfunctions 〈φ^{f}| of the Hamiltonian *H* = −Δ + *V* such that lim_{V→0}φ^{f} = *J*^{f} is a complete basis set with normalization condition . Thus, such states are automatically normalized by the density of states and one can just writewhere is now the Kronecker symbol. In other words, for a specific photon energy the summation is only on states *f* at the same energy .

Nonresonant inelastic X-ray scattering (NRIXS), or X-ray Raman spectroscopy, and resonant elastic X-ray scattering (REXS) are also governed by electronic transitions between core and unoccupied states. Thus, their response functions contain the same matrix elements, . This gives a possible unified way to calculate them numerically. Knowing ϕ_{g} and φ^{f}, the evaluation of the matrix elements and their summation do not present difficulties. ϕ_{g} is localized and is calculated using standard atomic procedures. The most important part thus concerns calculation of the photoelectron states. This is the purpose of the finite-difference method.

The main idea in FDM is to approximate the values of given continuous functions at a set of discrete points in a specific space. In the context of XAS, these functions are the wavefunctions corresponding to the states probed by the photoelectron, and the space is a volume limited by a sphere around the absorbing atom. Its radius is such that the outgoing photoelectron wave does not back-scatter towards the absorbing atom. This radius typically ranges from 5 to 8 Å.

The first step is thus to build a 3D grid of points. This grid can simply be cubic. Each point has six first-neighbour points and six second neighbours. We will see that the mesh can also be hexagonal. In other DFT techniques, the unknowns are the coefficients of the expansion of the wavefunctions in a specific basis. Here, the unknowns are simply the values of the wavefunction at each grid point *i*: .

The Laplacian in the Schrödinger equation is obtained by approximating the wavefunction around point *i* by a polynomial of second or fourth order. Its expression depends on the inter-point distance *d*. The smaller it is, the better the simulation. It is given bywhere Δ_{ij} are the coefficients of the Laplacian, with *j* indexing the *n* neighbouring points. When using a cubic grid and a second-order development, one just needs the *n* = 6 first neighbours, Δ_{ij} = 1/*d*^{2} and Δ_{ii} = −6/*d*^{2}. Equivalent precision can be reached with a fourth-order polynomial but with a less dense grid. The latter is preferred because it requires less computing power. In this case one also needs the second neighbours in the three directions for both orientations. Thus now *n* = 12, Δ_{ij} = 4/3*d*^{2} and −1/12*d*^{2} when *j* is the first and second neighbour, respectively, and Δ_{ii} = = −5/2*d*^{2}.

Finally, the discrete Schrödinger equation in atomic units (Bohr and Rydberg) at all points *i* iswhere *V*_{i} is the potential at point *i*. It is also discretized on the same grid. *V* is calculated using standard approaches with an adapted energy-dependent, exchange–correlation potential. We note that apart from the discretization, the potential does not need any approximation on its shape. This conceptual simplicity is the great advantage of FDM over other methods. Here, we took the potential as local, but it can be noted that FDM is perfectly adapted to nonlocal potentials, such as for instance in the case of a correlated material, which in the FDM form is writtenwith *j* not being limited here to the first neighbours. This possibility is not yet used in our approach.

When the material contains heavy atoms, it is most often necessary to perform a relativistic calculation. Following Wood & Boring (1978), we explicitly solve only the large component orbital instead of the four determined by the Dirac–Slater calculations. This reformulation gives a couple of Schrödinger-like equations, including the spin–orbit effect, closely akin to but improving on the Pauli equation, withwhere the potential, *V* = *V*(**r**, σ), and the wavefunction, φ = φ(**r**, σ), are dependent on the spin σ. **S** is a vector whose components are the Pauli matrices, and *B* = *B*(**r**, σ) = . The last term, which depends on **S** and is proportional to the potential gradient, represents the spin–orbit effect. By its *S*_{x} and *S*_{y} matrices, it connects the spin-up and spin-down components of the wavefunction. We also see that the gradient operator is needed; thus, the relativistic terms also give an off-diagonal effect. At fourth order, in a cubic grid, the FDM gradient is given by ∇_{ij} = ±2/3*d* and ∓1/12*d* for the first- and second-neighbouring points, respectively. Finally, at each grid point we obtain a couple of linear equations making the connection between the up and down wavefunctions,

As for any theory, there is great interest in using the symmetry of the system. Here, we work in direct space and consider a cluster of atoms limited by some radius. We thus have to consider its point group. Using the symmetry operations, plane symmetry, inversion and rotation axis (plus time reversal for magnetic systems), we need to calculate the wavefunctions in only a limited section of the original volume (see Fig. 1).

Due to the symmetry, off-diagonal Laplacian terms for points on a symmetry plane or rotation axis may be applied several times to the same neighbouring points, yielding the *H*_{ij} off-diagonal term. It is easy to check in Fig. 1 that the same equation at a neighbouring point which is not on the symmetry plane gives an off-diagonal *H*_{ji} ≠ *H*_{ij}. To recover the Hermitian matrix, we set the wavefunctions times the fraction, *v*_{i}, of the volume of the point (*d*^{3} for cubic grids) as unknowns instead of the simple wavefunction. For example, for a point on a symmetry plane *v*_{i} = 1/2. To do this, we just multiply the equation by *v*_{i}. In the nonrelativistic case, this givesWe have also multiplied and divided by the fraction, *v*_{j}, of the point *j*. Finally, using the new notation *v*_{i}φ_{i} → φ_{i}, we simply obtainIn this way we obtain *H*_{ji} = *H*_{ij}.

We previously gave the expression of the Laplacian for a cubic grid. More general formulae can be obtained for noncubic or non-uniform grids (Joly, 1992). Up to now, for XAS we have only used cubic and hexagonal grids of points (see Fig. 1). The latter is better for trigonal or hexagonal systems because the value of a missing neighbouring point in the Laplacian can be found versus an equivalent point in the grid. Of course, this value must then be multiplied by the character corresponding to the representation and the symmetry operation used. For a hexagonal grid there are six first and six second neighbours in the basal plane perpendicular to the threefold axis. The corresponding Laplacian coefficients are Δ_{ij} = 8/9*d*^{2} and Δ_{ij} = −1/18*d*^{2}, respectively. With the two first and two second neighbours along the axis, one obtains *n* = 16.

In a magnetic system, when neglecting the spin–orbit effect, up and down states are calculated independently and the point-group symmetry fully applies. When the spin–orbit effect cannot be neglected, time reversal, , has to be considered. In addition to the usual character, we must use .

We use a real-space cluster approach considering all of the atoms inside a so-called outer sphere. The calculation is performed in all of this volume, with the absorbing atom often, but not necessarily, at the centre. As seen in the following, the volume itself is divided into two zones: around the atomic cores limited by small atomic spheres and between the atoms.

We saw that we can obtain the solutions in the outer-sphere region from the known *J*^{f} solutions in vacuum. There, states *f* must be of the formwhere ; is the outgoing Hankel function. are the scattering amplitudes of the whole cluster. Without any atoms they are zero by the completeness theorem. The expansion in spherical harmonics defining the electronic states and the scattering amplitude is typically limited, as usual, by *l*_{max} ≃ κ*R*, where *R* is the cluster radius. Let us call *n*_{s} = (*l*_{max} + 1)^{2} the corresponding number of harmonics.

We need to write the analytic form of the embedding. One uses the continuity of the wave and of its derivative. For this purpose, one first writes for the neighbouring points just outside the outer sphere (see Fig. 2)

Secondly, one does the same for all of the points inside the outer sphere which have a first neighbour outside, and takes a Fourier-like transform, for a specific *L*′, summing over all these points,

Doing this for the *n*_{s} harmonics *L*′, we recover a number of equations equal to the number of outer-sphere scattering amplitudes. Note that the form an *n*_{s} × *n*_{s} nearly diagonal matrix because of the orthogonality of the spherical harmonics.

For systems with the spin–orbit effect we must also consider the spin σ; thus, the *f* states are doubled *f* = (*l*_{f}, *m*_{f}, σ_{f}) and we obtain spin-flip scattering amplitudes and write equation (14) in both spin situations,

Around the atomic core the wavefunctions have strong spatial oscillation, and a good reproduction of this would need very small inter-point distances. We have used a non-uniform grid in the past, but it is more convenient to use an expansion in spherical harmonics which does not necessarily need a spherical approximation of the potential shape. This full potential is given by , and one obtains a set of solutionswhere are the solutions of the coupled radial Schrödinger equations and are the atomic amplitudes obtained by continuity with the interatomic area. It is important to note that these expansions of both the potential and the wavefunction in spherical harmonics are performed in smaller atomic spheres than those used in classical multiple-scattering theory (MST). Indeed, here the corresponding radius typically ranges from 0.3 to 1 Å and increases with atomic number. Consequently, the potential can most often be limited to the first and isotropic term and indeed, to date, we have not found cases in which the following terms in the potential expansion give a real improvement. In the following we thus keep this approximation; the wavefunction then becomes diagonal and needs only the classical single-index *L* expansionwhere *b*_{l} is the solution of the standard radial Schrödinger equation.

The continuity with the FDM area is obtained in the same way as that between the outer sphere and the interior of the sphere. We use the continuity inside the atomic sphere to obtain the neighbouring values, with *B*_{Lj} = *b*_{l}*Y*_{L}. We also perform the Fourier-like transform on the points at the boundary of the ion core,

One thus obtains *n*_{a} = (*l*_{max} + 1)^{2} new equations corresponding to the *n*_{a} unknown atomic amplitudes. *l*_{max} is calculated here versus the atomic radius.

When considering the spin–orbit effect, starting from equations (7) and (8) we form a couple of Schrödinger-like radial equations. Their solutions take the formwhere indexes a solution. When neglecting the spin–orbit effect, *s* = σ. To make the match between the core and the grid of points we apply equation (19) for both spins.

From the FDM Schrödinger equation (equations 5, 7 and 8) applied to the *n*_{p} grid points and the Fourier-like transform for embedding in the outer sphere (equations 14 and 15) and with the atoms (equations 19 and 20) we form a large system of linear equations. It is an *n* × *n* matrix with *n* = , where *a* indexes the atoms in the cluster. *n* is multiplied by 2 when considering the spin–orbit effect.

In the XANES energy range, *d* = 0.25 Å is typically found to be sufficient at fourth order to obtain a good agreement. Bourke *et al.* (2016) also apply the FDM approach in the extended EXAFS range. For this purpose they use decreasing *d* values when the energy increases, down to 0.08 Å at 500 eV beyond the edge.

It is convenient to order the FDM points from the centre of the cluster to its outer sphere. In this way, the matrix is roughly band-type. When no symmetry applies *n* can be higher than 10 000, and the resolution of such a system is computationally demanding. Nevertheless the FDM matrix is sparse, and very efficient algorithms exist to solve it (Guda *et al.*, 2015).

Its solutions are the wavefunctions (times the volume fraction) at all of the node points, but also the atomic and outer-sphere amplitudes. From the atomic amplitudes, one can easily calculate the XANES, REXS or NRIXS signal response by applying the corresponding formula.

Applying the optical theorem and using the same normalization on the radial solutions of the Schrödinger equation, one obtains the following connection to the multiple-scattering amplitude, ,

This equation also applies in the spin–orbit case, with *L* → (*l*, *m* + σ − *s*, *s*). The different uses of the multiple-scattering theory to obtain response functions can thus also be applied when calculated by FDM. For example, one easily obtains the atomic electron density of state, as well as the atomic spin and orbital moments by the appropriate summations. One can also calculate the crystal orbital overlap population () between atoms *A* and *A*′ related to the multiple-scattering amplitude between the corresponding atoms, , with being the overlap integral. From the sign of (or its sum over *L* and *L*′) one has direct access to the bonding and antibonding states between atoms *A* and *A*′.

From the work of Schwitalla & Ebert (1998), we know that time-dependent DFT (TDDFT) simulations at the *L*_{2,3} edges of 3*d* elements provide a significant improvement over DFT. They first used a DFT simulation, on top of which they simulated the atomic multi-channel transition using an intra-atomic kernel. Bunău & Joly (2012*a*) followed the same scheme, with the DFT ingredient being provided by the FDM.

In this scheme, the DFT susceptibility is given in the spherical harmonics basis by and thus depends on the FDM atomic amplitudes. Formally, η is a positive number tending to zero. To mimic the finite core-hole lifetime it can be taken as equal to the equivalent level width. For simplicity we do not consider the spin and omit radial normalizations. The kernel, , node of this theory, allowing the crossing of transition channels between the different core and photoelectron states, is written in the same basis. One obtains the TDDFT susceptibility by a simple matrix expression . The matrix is not diagonal in core states and is introduced into the XANES formula in the same way as the multiple-scattering amplitude but with a double sum on core states: The interesting aspect of this theory is that it is fast in comparison with other *ab initio* multi-electronic techniques such as the Bethe–Salpeter equation, dynamical mean field theory or the multi-configuration method. It improves the ratio between the *L*_{2} and *L*_{3} edges, but performed here only with a local kernel it does not give the excitonic features. It is also independent of the calculation method used to obtain the DFT susceptibility. We recently analyzed the improvement provided by FDM with respect to the MST when coupled with TDDFT in various compounds (Bunău & Joly, 2012*b*).

The *FDMNES* code (Bunău & Joly, 2009; Bunău *et al.*, 2021) can use FDM or MST on demand. The validity of the code was first tested by a comparison of the techniques. To improve the validity of the test, we imposed the muffin-tin (MT) approximation on the potential shape (a spherical average within the atoms and a constant potential between them), without overlap, when working in FDM mode. Indeed, as is most often the case, MST uses the MT approximation in *FDMNES*. In this way, we have first shown that both approaches are equivalent in this case (Joly, 2001). More importantly, it has been possible to show that most of the limitation of the MT approximation is due its constant potential in the interatomic region. Roughly speaking, when the material is metallic and dense, when the structure is disordered or when the core-hole width is important, the MT approximation most often gives good results.

Here, we present some typical examples of X-ray spectroscopies calculated using FDM. In two cases FDM shows a clear improvement thanks to its free-shape potential. In the NdMg case it is mandatory because the observed effect comes from the aspherical potential around the atoms.

The edge of low-*Z* elements very often needs a free-shape potential. In the case of organic molecules it is often absolutely necessary. Here, we give an example in a condensed matter system: TiO_{2} in its rutile structure. The complete study and the effect of co-doping with iron and nitrogen can be found in Kaspar *et al.* (2012). Here, we have just recalculated the linear dichroism at the O *K* edge; that is, the difference in signal when the polarization is along the *c* axis and perpendicular to it. Fig. 3 shows the FDM result, the data from the reference and what we obtain using the MST muffin-tin calculation (with a 15% MT radius overlap). The free-shape calculation is evidently better.

High-*Z* elements such as mercury can also be sensitive to the full-shape potential. Manceau *et al.* (2015) showed that in mercury complexes recorded at high resolution, that is with a lower broadening effect than in standard measurements, it was impossible to obtain a convenient accord using the muffin-tin approximation. In contrast, a satisfactory agreement was obtained when working with the FDM procedure. Moreover, it was then possible to make an assignment of the different features in the spectra in term of (*l*, *m*) projected densities of state. For example, in Fig. 4 we show the data and simulation with and without convolution.

XMCD simulations at the *K* edge are not so easy to perform because they need to take the spin–orbit effect into account in the calculation of the photoelectron states, and the signal is always small because the 4*p* probed states are only slightly magnetically polarized. We illustrate a typical result obtained at the Fe *K* edge for a simple metallic foil. FDM simulations are satisfactory from a 4.5 Å cluster radius. In Fig. 5 we show the agreement between an SCF calculation using a 7 Å cluster radius compared with the data from Mathon *et al.* (2004). We can check that this is a case in which the MST muffin-tin approximation is good and that the best agreement is close to the Fermi energy where the magnetically polarized 3*d* states are found. At a higher energy most of the effect is a derivative effect that is more difficult to reproduce quantitatively.

Orbital ordering is a typical case in which the resulting aspherical potential can be described well by FDM. This possibility is interesting because the purpose is to find a way to observe this phenomenon independently of the associated structural or magnetic loss of symmetry. For example, in the case of LaMnO_{3} we proved that the REXS signal observed at the (300) reflection at the Mn *K* edge is simply due to Jahn–Teller distortion and not to the Coulomb effect induced by the 3*d* Mn orbital ordering (Benfatto *et al.*, 1999).

More interestingly, in NdMg Bunău *et al.* (2010) have shown that a quadrupolar orbital ordering associated with a peculiar spin ordering can be seen in the nonmagnetic reflection at the Nd *L*_{2} and *L*_{3} edges. The simulated signal including both the spin–orbit effect and orbital ordering is ten times stronger that that with only the spin–orbit effect (see Fig. 6). The signal with no spin–orbit effect and no orbital ordering is zero.

### References

Benfatto, M., Joly, Y. & Natoli, C. R. (1999).*Phys. Rev. Lett.*

**83**, 636–639.Google Scholar

Bourke, J. D., Chantler, C. T. & Joly, Y. (2016).

*J. Synchrotron Rad.*

**23**, 551–559.Google Scholar

Bunău, O., Galéra, R. M., Joly, Y., Amara, M., Luca, S. E. & Detlefs, C. (2010).

*Phys. Rev. B*,

**81**, 144402.Google Scholar

Bunău, O. & Joly, Y. (2009).

*J. Phys. Condens. Matter*,

**21**, 345501.Google Scholar

Bunău, O. & Joly, Y. (2012

*a*).

*Phys. Rev. B*,

**85**, 155121.Google Scholar

Bunău, O. & Joly, Y. (2012

*b*).

*J. Phys. Condens. Matter*,

**24**, 215502.Google Scholar

Bunău, O., Ramos, A. Y. & Joly, Y. (2021).

*Int. Tables Crystallogr. I*, https://doi.org/10.1107/S1574870720003304 .Google Scholar

Deift, P. & Simon, B. (1976).

*J. Funct. Anal.*

**23**, 218–238.Google Scholar

Grossmann, C., Roos, H.-G. & Stynes, M. (2007).

*Numerical Treatment of Partial Differential Equations.*Berlin, Heidelberg: Springer.Google Scholar

Guda, S. A., Guda, A. A., Soldatov, M. A., Lomachenko, K. A., Bugaev, A. L., Lamberti, C., Gawelda, W., Bressler, C., Smolentsev, G., Soldatov, A. V. & Joly, Y. (2015).

*J. Chem. Theory Comput.*

**11**, 4512–4521.Google Scholar

Joly, Y. (1992).

*Phys. Rev. Lett.*

**68**, 950–953.Google Scholar

Joly, Y. (1994).

*Phys. Rev. Lett.*

**72**, 392–395.Google Scholar

Joly, Y. (1997).

*J. Phys. IV Fr.*

**7**, C2-111–C2-115.Google Scholar

Joly, Y. (2001).

*Phys. Rev. B*,

**63**, 125120.Google Scholar

Joly, Y., Cabaret, D., Renevier, H. & Natoli, C. R. (1999).

*Phys. Rev. Lett.*

**82**, 2398–2401.Google Scholar

Kaspar, T. C., Ney, A., Mangham, A. N., Heald, S. M., Joly, Y., Ney, V., Wilhelm, F., Rogalev, A., Yakou, F. & Chambers, S. A. (2012).

*Phys. Rev. B*,

**86**, 035322.Google Scholar

Kimball, G. E. & Shortley, G. H. (1934).

*Phys. Rev.*

**45**, 815–820.Google Scholar

Manceau, A., Lemouchi, C., Rovezzi, M., Lanson, M., Glatzel, P., Nagy, K. L., Gautier-Luneau, I., Joly, Y. & Enescu, M. (2015).

*Inorg. Chem.*

**54**, 11776–11791.Google Scholar

Mathon, O., Baudelet, F., Itié, J.-P., Pasternak, S., Polian, A. & Pascarelli, S. (2004).

*J. Synchrotron Rad.*

**11**, 423–427.Google Scholar

Nieminen, R. M. & Puska, M. J. (1983).

*Phys. Rev. Lett.*

**50**, 281–284.Google Scholar

Schwitalla, J. & Ebert, H. (1998).

*Phys. Rev. Lett.*

**80**, 4586–4589.Google Scholar

Wood, J. H. & Boring, A. M. (1978).

*Phys. Rev. B*,

**18**, 2701–2711.Google Scholar