Chapter 2.11. Finite-difference method for the calculation of X-ray spectroscopies

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² and Δ_ii = −6/d². 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/3d² and −1/12d² when j is the first and second neighbour, respectively, and Δ_ii = = −5/2d².

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/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,

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).

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 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³ 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 non­cubic 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/9d² and Δ_ij = −1/18d², 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)² 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)

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,

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_lY_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)² 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.

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₂ 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.

Figure 3 Linear dichroism at the O K edge in rutile TiO2. Top: experiment from Kaspar et al. (2012). 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.

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.

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)

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 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 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 XMCD at the Fe K edge in a simple metallic iron foil. The data (dotted line) are from Mathon et al. (2004); the full line shows the FDM simulation and the dashed line shows the MST simulation.

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₃ 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).

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₂ and L₃ 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.

Figure 6 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.