Tables for
Volume I
X-ray absorption spectroscopy and related techniques
Edited by C. T. Chantler, F. Boscherini and B. Bunker

International Tables for Crystallography (2022). Vol. I. Early view chapter

Finite-difference method for the calculation of X-ray spectroscopies

Y. Joly,a* A. Y. Ramosa and O. Bunăua

aUniversité Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble, France
Correspondence e-mail:

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 FDMNES software.

Keywords: finite-difference method.

1. Introduction

The finite-difference method (FDM) is a general technique (Grossmann et al., 2007[link]) for solving partial differential equations. Its first formulation for solving the Schrödinger equation was given in the 1930s (Kimball & Shortley, 1934[link]). 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[link]). We first developed it for scattering simulations of electrons (Joly, 1992[link]) and positrons (Joly, 1994[link]). FDM was applied for the first time to simulate X-ray absorption near-edge structure (XANES) at the end of the 1990s (Joly, 1997[link]; Joly et al., 1999[link]) in the FDMNES code presented in Bunău et al. (2021[link]). 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.

2. XANES and related spectroscopies

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 L1 edges of any elements and at the L2,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 [\sigma = 4\pi^{2}\alpha\hbar\omega S_{0}^{2}\textstyle\sum\limits_{fg}\left|\left\langle\phi_{f} \right|\hat{o}\left|\phi_{g}\right\rangle\right|^{2}\delta\left({\cal E}_{f}-{\cal E}_{g}+\Delta{\cal E}_{a}-\hbar\omega\right), \eqno (1)]where α is the fine-structure constant. [{\cal E}_{g}], [{\cal E}_{f}] and ℏω are the core-state, final-state and incoming photon energies, respectively. [\delta({\cal E}_{f}-{\cal E}_{g}+\Delta{\cal E}_{a}-\hbar\omega)] is the density of state at [{\cal E}_{f}]. Note that the factor [S_{0}^{2}] ([S_{0}^{2}\,\lt\, 1]) and the extra term in the energy conservation, [\Delta{\cal E}_{a}], give the effect of the intra-atomic screening of the core hole. For simplicity, as often performed, we omit both terms in the following. [\hat{o}] is the electromagnetic operator. Neglecting the magnetic inter­action term, and with expansion up to the quadrupole electric term, it is given by [\hat{o} = \boldvarepsilon\cdot{\bf r}\left(1+{{i} \over {2}}{\bf k}\cdot{\bf r} \right), \eqno (2)]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, [j_{l_{f}}], and we know that [J^{f} = (\kappa/\pi)^{1/2}j_{l_{f}}(\kappa r)Y_{L_{f}}(\hat{r})] are solutions of the Schrödinger equation when the potential V is zero, and that they form a complete basis. [\kappa = ({\cal E}_{f})^{1/2}] is the modulus of the electron wavevector in atomic units. r and [\hat{r}] are the radial and angular positions. The indices f are now defined by the quantum numbers Lf = (lf, mf) and the energy [{\cal E}_{f}]. 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[link]) states that a set of eigenfunctions 〈φf| of the Hamiltonian H = −Δ + V such that limV→0φf = Jf is a complete basis set with normalization condition [|\varphi^{f}\rangle\langle\varphi^{f}| = |\phi_{f} \rangle\langle\phi_{f}|\delta({\cal E}_{f})]. Thus, such states are automatically normalized by the density of states and one can just write[\sigma = 4\pi^{2}\alpha\hbar\omega\textstyle\sum \limits_{fg}\left|\left\langle\varphi^{f}\right| \hat{o}\left|\phi_{g}\right\rangle\right|^{2}\delta_{k}\left({\cal E}_{f}- {\cal E}_{g}-\hbar\omega\right), \eqno (3)]where [\delta_{k}({\cal E}_{f}-{\cal E}_{g}-\hbar\omega)] is now the Kronecker symbol. In other words, for a specific photon energy the summation is only on states f at the same energy [{\cal E}_{f}].

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, [\langle\varphi^{f}|\hat{o}|\phi_{g}\rangle]. 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.

3. FDM for absorption spectroscopies

3.1. Discrete formulation

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: [\varphi^{f}_{i} = \varphi^{f}({\bf r}_{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 by[\Delta\varphi^{f}_{i} = \textstyle \sum \limits_{j}\Delta_{ij}\varphi^{f}_{j}+\Delta_{ii}\varphi^{f}_{i}, \eqno (4)]where Δ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/d2 and Δii = −6/d2. 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/3d2 and −1/12d2 when j is the first and second neighbour, respectively, and Δii = [-\sum_{j}\Delta_{ij}] = −5/2d2.

Finally, the discrete Schrödinger equation in atomic units (Bohr and Rydberg) at all points i is[(-\Delta_{ii}+V_{i}-{\cal E}_{f})\varphi^{f}_{i}-\textstyle\sum\limits_{j}\Delta_{ij}\varphi^{f}_{j} = 0, \eqno (5)]where Vi 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 written[V\varphi^{f}_{i} = \textstyle\sum\limits_{j}V_{ij}\varphi^{f}_{j}, \eqno (6)]with j not being limited here to the first neighbours. This possibility is not yet used in our approach.

3.2. Relativistic case and spin–orbit effect

When the material contains heavy atoms, it is most often necessary to perform a relativistic calculation. Following Wood & Boring (1978[link]), 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, [(-\Delta+V-{\cal E}_{f}+\Omega)\varphi^{f} = 0, \eqno (7)]with[\Omega = -{{\alpha^{2}} \over {4} }\{(V-{\cal E}_{f})^{2}-B[\nabla V\cdot\nabla+i{\bf S}\cdot(\nabla V\times\nabla)]\} = 0, \eqno (8)]where 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, σ) = [1/[1-(\alpha^{2}/4)(V - {\cal E})^{2}]]. The last term, which depends on S and is proportional to the potential gradient, represents the spin–orbit effect. By its Sx and Sy 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/3d and ∓1/12d 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,[\eqalignno{(-\Delta_{ii}& +V_{i\sigma}-{\cal E}_{f}+\Omega_{i\sigma i \sigma})\varphi^{f}_{i\sigma}+\Omega_{i\sigma i,-\sigma}\varphi^{f}_{i,-\sigma} \cr &+\textstyle \sum \limits_{j\sigma^{\prime}}(-\Delta_{ij}+\Omega_{i\sigma j \sigma^{\prime}})\varphi^{f}_{j\sigma^{\prime}} = 0. &(9)}]

3.3. FDM and symmetry

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[link]).

[Figure 1]

Figure 1

FDM grids of points. Left: a cubic grid corresponding to point group 4mm. The volume corresponding to each point is shown by the thin lines. For this symmetry group only the points shown are necessary. To calculate the Laplacian when a neighbouring point is outside the grid, one uses the plane symmetry to obtain an equivalent point. Right: for trigonal and hexagonal symmetry, the grid of points in the plane perpendicular to the threefold axis is hexagonal. When neighbouring points are outside the area, one uses the rotation symmetry to obtain a point within the area.

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 Hij off-diagonal term. It is easy to check in Fig. 1[link] that the same equation at a neighbouring point which is not on the symmetry plane gives an off-diagonal HjiHij. To recover the Hermitian matrix, we set the wavefunctions times the fraction, vi, of the volume of the point (d3 for cubic grids) as unknowns instead of the simple wavefunction. For example, for a point on a symmetry plane vi = 1/2. To do this, we just multiply the equation by vi. In the nonrelativistic case, this gives[(-\Delta_{ii}+V_{i}-{\cal E}_{f})v_{i}\varphi_{i}-{\textstyle \sum\limits_{j}\Delta _{ij}}{{v_{i}} \over {v_{j}}}v_{j}\varphi^{f}_{j} = 0.\eqno (10)]We have also multiplied and divided by the fraction, vj, of the point j. Finally, using the new notation viφiφi, we simply obtain[(-\Delta_{ii}+V_{i}-{\cal E}_{f})\varphi_{i}-{\textstyle\sum\limits_{j}}\Delta_{ij} {{v_{i}} \over {v_{j}}}\varphi^{f}_{j} = 0. \eqno (11)]In this way we obtain Hji = Hij.

We previously gave the expression of the Laplacian for a cubic grid. More general formulae can be obtained for non­cubic or non-uniform grids (Joly, 1992[link]). Up to now, for XAS we have only used cubic and hexagonal grids of points (see Fig. 1[link]). The latter is better for trigonal or hexagonal systems because the [\varphi^{f}_{j}] 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/9d2 and Δij = −1/18d2, 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, [\hat{T}], has to be considered. In addition to the usual character, we must use [\hat{T}\varphi_{i,\sigma} = \varphi^{*}_{i,-\sigma}].

3.4. Cluster approach

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 Jf solutions in vacuum. There, states f must be of the form[\varphi^{f}({\bf r},{\cal E}_{f}) = J^{f}({\bf r}, {\cal E}_{f})-i\textstyle\sum\limits_{L}s_{L}^{f}({\cal E}_{f})H_{L}({\bf r},{\cal E}_{f}), \eqno (12)]where [H_{L} = ({k/\pi})^{1/2}h^{+}_{l}(\kappa r)Y_{L}({\hat r})]; [h^{+}_{l}] is the outgoing Hankel function. [s_{L}^{f}] 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 lmaxκR, where R is the cluster radius. Let us call ns = (lmax + 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[link]) [\varphi^{f}_{j} = J^{f}_{j}-i\textstyle\sum\limits_{L}s_{L}^{f}H_{Lj}. \eqno (13)]

[Figure 2]

Figure 2

Cluster with point group mm2 and its FDM grid of points. The atomic cores, where expansion in spherical waves is performed, are also shown. The points at the boundary of the outer sphere with a grey background are used to make a Fourier-like transform to ensure the continuity of the wavefunction and of its derivative. The empty circles show the points in the outer sphere which are neighbours of the previous ones calculated using the expansion in spherical harmonics.

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, [\textstyle \sum \limits_{i}\varphi^{f}_{i}Y^{*}_{L^{\prime}i} = \textstyle \sum \limits_{i}J^{f}_{i}Y^{*}_{L^{\prime}i }-i\textstyle \sum \limits_{L}s_{L}^{f}\textstyle \sum \limits_{i}H_{Li}Y^{*}_{L^{\prime}i}, \eqno (14)]

Doing this for the ns harmonics L′, we recover a number of equations equal to the number of outer-sphere scattering amplitudes. Note that the [\textstyle \sum_{i}H_{Li}Y^{*}_{L^{\prime}i}] form an ns × ns 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 = (lf, mf, σf) and we obtain spin-flip scattering amplitudes [s_{L\sigma}^{L_{f}\sigma_{f}}] and write equation (14)[link] in both spin situations, [\textstyle \sum \limits_{i}\varphi^{f}_{i\sigma}Y^{*}_{L^{\prime}i} = \textstyle \sum \limits_{i}J^{f}_{i}Y^{*}_{L^{\prime}i}\delta_{\sigma\sigma_{f}}-i\textstyle \sum \limits_{L}s_{L\sigma}^{f}\sum_{i}H_{Li}Y^{*} _{L^{\prime}i}. \eqno (15)]

3.5. Atomic cores

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 [V = \textstyle \sum_{L}V_{L}Y_{L}], and one obtains a set of solutions[\varphi^{f}({\bf r},{\cal E}_{f}) = \textstyle \sum \limits_{L}a_{L}^{f}({\cal E}_{f})\textstyle \sum \limits_{L^{\prime}}b_{LL^{\prime}}(r,{\cal E}_{f} )Y_{L^{\prime}}({\hat r}), \eqno (16)]where [b_{LL^{\prime}}] are the solutions of the coupled radial Schrödinger equations and [a_{L}^{f}] 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 expansion[\varphi^{f}({\bf r},{\cal E}_{f}) = \textstyle \sum \limits_{L}a_{L}^{f}({\cal E}_{f})b_{l}(r,{\cal E}_{f})Y_{L}({\hat r}), \eqno (17)]where bl 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, [\varphi^{f}_{j} = \textstyle \sum \limits_{L}a_{L}^{f}B_{Lj}, \eqno (18)]with BLj = blYL. We also perform the Fourier-like transform on the points at the boundary of the ion core, [\textstyle \sum \limits_{i}\varphi^{f}_{i}Y^{*}_{L^{\prime}i} = \textstyle \sum \limits_{L}a_{L}^{f}\textstyle \sum \limits_{i}B_{Li}Y^{* }_{L^{\prime}i}.\eqno (19)]

One thus obtains na = (lmax + 1)2 new equations corresponding to the na unknown atomic amplitudes. lmax is calculated here versus the atomic radius.

When considering the spin–orbit effect, starting from equations (7)[link] and (8)[link] we form a couple of Schrödinger-like radial equations. Their solutions take the form[\varphi^{f}_{\sigma}({\bf r},{\cal E}_{f}) = \textstyle \sum \limits_{Ls}a_{l,m +\sigma-s,s}^{f}({\cal E}_{f})b_{L\sigma s}(r,{\cal E}_{f})Y_{L}(\hat{r}), \eqno (20)]where [s = \pm{{1} \over {2}}] 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)[link] for both spins.

3.6. The FDM matrix

From the FDM Schrödinger equation (equations 5[link], 7[link] and 8[link]) applied to the np grid points and the Fourier-like transform for embedding in the outer sphere (equations 14[link] and 15[link]) and with the atoms (equations 19[link] and 20[link]) we form a large system of linear equations. It is an n × n matrix with n = [n_{p}+\textstyle\sum_{a}n_{a}+n_{s}], 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[link]) 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[link]).

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.

3.7. Connection to MST and time-dependent DFT

3.7.1. MST

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, [\tau_{LL^{\prime}}], [{{1} \over {2}}(\tau_{LL^{\prime}}-\tau_{L^{\prime}L}^{*}) = {\hat \tau}_{ LL^{\prime}} = -i\textstyle\sum\limits_{f}a_{L}^{f}a_{L^{\prime}}^{f*}. \eqno (21)]

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, [\rho_{L}({\bf r}) = -\Im({\hat \tau}_{LL})|b_{L} ({\bf r})|^{2} = \textstyle\sum\limits_{f}|a_{L}^{f}|^{2}|b_{L}({\bf r})|^{2}, \eqno (22)]as well as the atomic spin and orbital moments by the appropriate summations. One can also calculate the crystal orbital overlap population ([O_{LAL^{\prime}A^{\prime}}]) between atoms A and A′ related to the multiple-scattering amplitude between the corresponding atoms, [{\hat \tau}_{LAL^{\prime}A^{\prime}} = -i\textstyle\sum_{f}a_{LA}^{f}a_{L^{\prime}A^{\prime}}^{f*}], [O_{LAL^{\prime}A^{\prime}} = 2S_{LAL^{\prime}A^{\prime}}\textstyle\sum\limits_{f}a_{LA}^{f}a_{L^{\prime}A^{\prime}}^{f*}\int B_{LA}({\bf r})B_{L^{\prime}A^{\prime}}^{*}({\bf r})\,{\rm d}({\bf r}), \eqno(23)]with [S_{LAL^{\prime}A^{\prime}}] being the overlap integral. From the sign of [O_{LAL^{\prime}A^{\prime}}] (or its sum over L and L′) one has direct access to the bonding and antibonding states between atoms A and A′.

3.7.2. Time-dependent DFT

From the work of Schwitalla & Ebert (1998[link]), we know that time-dependent DFT (TDDFT) simulations at the L2,3 edges of 3d 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 (2012a[link]) 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 [\hat{\chi}_{0LL^{\prime}gg^{\prime}}\simeq\delta_{gg^{\prime}}{{i} \over {\pi}}\, {\textstyle \int \limits}{\rm d}{\cal E}_{f}{{\hat{\tau}_{LL^{\prime}}} \over {\hbar\omega-({\cal E}_{f}-{\cal E}_{g})+i\eta}} \eqno (24)]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, [{\hat K}], 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 [\hat{\chi} = (1-\hat{\chi}_{0}\hat{K})^{-1}\hat{\chi}_{0}]. The [\hat{\chi}_{LL^{\prime}gg^{\prime}}] 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: [\sigma = -4\pi^{2}\alpha\hbar\omega\textstyle \sum\limits_{gg^{\prime}}\sum\limits_{LL^{\prime}}\langle g|\hat{o}^{*}|b_{l}Y_{L}\rangle\Im(\hat{\chi} _{LL^{\prime}gg^{\prime}})\langle b_{l^{\prime}}Y_{L^{\prime}}|\hat{o}|g^{\prime}\rangle. \eqno (25)]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 L2 and L3 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, 2012b[link]).

4. Applications and testing

The FDMNES code (Bunău & Joly, 2009[link]; Bunău et al., 2021[link]) 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[link]). 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.

4.1. Low-Z materials

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: TiO2 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[link]). 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[link] 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.

[Figure 3]

Figure 3

Linear dichroism at the O K edge in rutile TiO2. Top: experiment from Kaspar et al. (2012[link]). Bottom: simulations with FDM (thick line) and muffin-tin MST (thin line). The dotted line is the polarization parallel to the c axis; the full line is that perpendicular.

4.2. High-Z materials

High-Z elements such as mercury can also be sensitive to the full-shape potential. Manceau et al. (2015[link]) 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[link] we show the data and simulation with and without convolution.

[Figure 4]

Figure 4

Hg L3 edge in Hg(Cys)4, a four-coordinated mercury thiolate. The dotted line is the high-resolution experiment; the FDM simulations are shown convoluted (full line) and not convoluted (dashed line). The data and study are from Manceau et al. (2015[link])

4.3. X-ray magnetic circular dichroism (XMCD)

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 4p 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[link] we show the agreement between an SCF calculation using a 7 Å cluster radius compared with the data from Mathon et al. (2004[link]). 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 3d states are found. At a higher energy most of the effect is a derivative effect that is more difficult to reproduce quantitatively.

[Figure 5]

Figure 5

XMCD at the Fe K edge in a simple metallic iron foil. The data (dotted line) are from Mathon et al. (2004[link]); the full line shows the FDM simulation and the dashed line shows the MST simulation.

4.4. Resonant diffraction in NdMg

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 LaMnO3 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 3d Mn orbital ordering (Benfatto et al., 1999[link]).

More interestingly, in NdMg Bunău et al. (2010[link]) have shown that a quadrupolar orbital ordering associated with a peculiar spin ordering can be seen in the nonmagnetic [(-{{1} \over {2}},0,{{5} \over {2}})] reflection at the Nd L2 and L3 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[link]). The signal with no spin–orbit effect and no orbital ordering is zero.

[Figure 6]

Figure 6

[(-{{1} \over {2}},0,{{5} \over {2}})] reflection at the Nd L3 edge in NdMg. Dotted line, experiment; full line, simulation using FDM with spin–orbit effect and quadrupolar orbital ordering; dashed line, the same simulation but without orbital ordering.


First citationBenfatto, M., Joly, Y. & Natoli, C. R. (1999). Phys. Rev. Lett. 83, 636–639.Google Scholar
First citationBourke, J. D., Chantler, C. T. & Joly, Y. (2016). J. Synchrotron Rad. 23, 551–559.Google Scholar
First citationBunău, O., Galéra, R. M., Joly, Y., Amara, M., Luca, S. E. & Detlefs, C. (2010). Phys. Rev. B, 81, 144402.Google Scholar
First citationBunău, O. & Joly, Y. (2009). J. Phys. Condens. Matter, 21, 345501.Google Scholar
First citationBunău, O. & Joly, Y. (2012a). Phys. Rev. B, 85, 155121.Google Scholar
First citationBunău, O. & Joly, Y. (2012b). J. Phys. Condens. Matter, 24, 215502.Google Scholar
First citationBunău, O., Ramos, A. Y. & Joly, Y. (2021). Int. Tables Crystallogr. I, .Google Scholar
First citationDeift, P. & Simon, B. (1976). J. Funct. Anal. 23, 218–238.Google Scholar
First citationGrossmann, C., Roos, H.-G. & Stynes, M. (2007). Numerical Treatment of Partial Differential Equations. Berlin, Heidelberg: Springer.Google Scholar
First citationGuda, 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
First citationJoly, Y. (1992). Phys. Rev. Lett. 68, 950–953.Google Scholar
First citationJoly, Y. (1994). Phys. Rev. Lett. 72, 392–395.Google Scholar
First citationJoly, Y. (1997). J. Phys. IV Fr. 7, C2-111–C2-115.Google Scholar
First citationJoly, Y. (2001). Phys. Rev. B, 63, 125120.Google Scholar
First citationJoly, Y., Cabaret, D., Renevier, H. & Natoli, C. R. (1999). Phys. Rev. Lett. 82, 2398–2401.Google Scholar
First citationKaspar, 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
First citationKimball, G. E. & Shortley, G. H. (1934). Phys. Rev. 45, 815–820.Google Scholar
First citationManceau, 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
First citationMathon, O., Baudelet, F., Itié, J.-P., Pasternak, S., Polian, A. & Pascarelli, S. (2004). J. Synchrotron Rad. 11, 423–427.Google Scholar
First citationNieminen, R. M. & Puska, M. J. (1983). Phys. Rev. Lett. 50, 281–284.Google Scholar
First citationSchwitalla, J. & Ebert, H. (1998). Phys. Rev. Lett. 80, 4586–4589.Google Scholar
First citationWood, J. H. & Boring, A. M. (1978). Phys. Rev. B, 18, 2701–2711.Google Scholar

to end of page
to top of page