Full Breit Hamiltonian in the Multiwavelets Framework

New techniques in core–electron spectroscopy are necessary to resolve the structures of oxides of f-elements and other strongly correlated materials that are present only as powders and not as single crystals. Thus, accurate quantum chemical methods must be developed to calculate core spectroscopic properties in such materials. In this contribution, we present an important development in this direction, extending our fully adaptive real-space multiwavelet basis framework to tackle the four-component Dirac-Coulomb-Breit Hamiltonian. We show that multiwavelets can reproduce one-dimensional grid-based approaches. They are however a fully three-dimensional approach which can later be extended to molecules and materials. Our multiwavelet implementation attained precise results irrespective of the chosen nuclear model, provided that the error threshold is tight enough and that the chosen polynomial basis is sufficiently large. Furthermore, our results confirmed that in two-electron species, the magnetic and Gauge contributions from s-orbitals are identical in magnitude and can account for the experimental evidence from K and L edges.


Introduction
Core-electron spectroscopies like X-ray photoelectron spectroscopy, X-ray absorption spectroscopy and electron energy loss spectroscopy are powerful tools to investigate the electronic structure of transition-metal and rare-earth materials.For example, multi-layered transition metal carbides and carbonitrides M n+1 AX n , where M is an early transition metal, A is an Agroup element (mostly groups 13 and 14), X is C or/and N and n is 1 to 3. 1 These materials can be employed for energy storage systems, such as lithium-ion batteries, [1][2][3][4][5] lithium-ion capacitors, 6 aqueous pseudocapacitors, 7,8 and transparent conductive films. 9[12] Unfortunately, their spectra are not straightforwardly interpretable due to relativistic effects.4][15][16][17][18][19][20][21][22] A computational approach based on first-principle calculations that will take into account both relativity and electron correlation could help the interpretation of such spectra.A recent, promising approach in quantum chemistry is based on multiresolution analysis (MRA), by making use of multiwavelets (MWs). 235][26][27][28] A variational treatment of relativistic effects into MRA will allow modelling the spectra of transition metal and rare earth materials.An important step in this direction was presented to tackle the mean-field atomic and molecular Dirac-Coulomb problem in an adaptive, 4-component multiwavelets basis. 29,30 such a model the electrons are considered static charges where the average interaction between electrons is modelled with the Coulomb-like term only.This is the lowest-order relativistic approximation for the two-electron interaction, which disregards the magnetic interactions, such as spin-other-orbit, and the retardation effects due to the finite speed of light.These effects are important and must be taken into account for a realistic modelling of core-electron spectroscopies.6][37][38][39] The Breit Hamiltonian adds two negative terms, called Gaunt and Gauge, respectively: The first term in Eq. ( 1) is the non-relativistic Coulomb interaction.The second term, called Gaunt, can be seen, in the non-relativistic limit, as the scalar product between the curl of two spin orbitals: ⃗ α i ∼ ∇ × ϕ i . 37⃗ α denotes a Cartesian vector collecting the 4 × 4 Dirac matrices α x , α y and α z (vide infra).When ⃗ α acts on a 4-component orbital, it couples its components, as detailed later on in this contribution.This means that the spin rotation of one electron on its axis generates a vector potential that will interact with the vector potentials generated by all other electrons present in the system, 37 resulting in a scalar potential.Finally, the third term, called Gauge, describes the retardation effects due to the reciprocal interaction between the rotational vector fields (α i • ∇ i ) of two electrons. 37[33][34] In this contribution, we will present the adaptive MRA multiwavelet implementation of the full Breit interaction as a perturbative correction on top of a 4-component Dirac-Coulomb-Hartree-Fock wavefunction.We will demonstrate the precision of our implementation by comparing ground-state energies of highly-charged helium-like ions with increasing Z, X (Z -2)+ , performed with our Python code, VAMPyR (Very Accurate Multiresolution Python Routines) 40 with numerical radial integration in GRASP 41 and Gaussian basis set calculations with the DIRAC 42 software.
2 Theory and Implementation

Multiresolution Analysis and Multiwavelets
Multiresolution analysis 43 is constructed by considering a set of orthonormal functions called scaling functions ϕ i (x) supported on the interval [0, 1].They can be dilated and translated to obtain a corresponding basis in subintervals of [0, 1].The most common procedure is a dyadic subdivision, such that at scale n there will be 2 n intervals defined by a translation index l = 0, 2 n − 1 such that the scaling functions in the l-th interval [l/2 n , (l + 1)/2 n ] are obtained as: Additionally, functions at subsequent scales are connected by the two-scale relationships which allow to obtain the scaling function at scale n as a linear combination of scaling functions at scale n − 1.
This construction leads to a ladder of scaling spaces in a telescopic sequence which is dense in L 2 : The multiwavelet functions are then obtained as the orthogonal complement of the scaling functions at scale n + 1 with respect to the ones at scale n.
In the construction of Alpert, 44 the scaling functions are a simple set of polynomials, and the wavelet functions are then piecewise polynomial functions.The possibility to construct efficient algorithms, with precise error control, relies on the combination of several properties of such a construction.Here, it will suffice to say that the most important aspects concern the disjoint support of the basis, which enables function-based adaptivity, the vanishing moments of the wavelet functions, which guarantees fast decay of the representation coefficients, the non-standard (NS) form of operators, 45 which uncouples scales during operator application thus preserving adaptivity, the separated representation of integral kernels, which leads to low-scaling algorithms.The interested reader is referred to the available literature for details about those aspects. 23,44,46,472 Mean-field two-electron operators in a Multiwavelets Basis We will summarise the main methodological developments enabling the results in this contribution.We first recall that in a relativistic framework, molecular orbitals are vectors with four complex components.We will use indices: • u, w ∈ {x, y, z} for Cartesian components, • p, q, . . .for occupied 4-component orbitals, • A, B, . . .∈ {1, 2, 3, 4} for orbital components.
Furthermore, Greek capital letters will be used for the 4-component orbitals and their lowercase counterparts will be used for the corresponding components: The corresponding Hermitian conjugate (transposed and complex conjugate) orbital is: with † denoting Hermitian conjugation and overline complex conjugation of a component.
To avoid confusion we will also refer to the instantaneous electron interaction (first term in Eq. 1 as the Coulomb term, whereas we will use the terms direct and exchange to refer to the two parts of each term, arising from the fermionic nature of the electrons. For the Coulomb operator g Coulomb (⃗ r 1 , ⃗ r 2 ) = I 1 •I 2 r 12 , the direct and exchange operators are straightforward and shown in Eqs.(7a) and (7b) in the Supporting Information, respectively.
In practice, these operators are applied as convolutions.Efficient and accurate convolution with an integral operator is implemented in a separated representation (see Ref. 47 for details).We underline that the Coulomb part of the two-electron interaction is in this framework diagonal, in the sense that it is not coupling the four components of the spinor.
In a Gaussian Type Orbital (GTO) framework, the exchange part would instead couple the four components of the spinor, because the formalism is tied to the atomic orbital (AO) densities, thus generating an artificial coupling once the exchange operation is performed. 48 proceed similarly for the Gaunt operator g appearing in the numerator are Cartesian vectors whose components are 4 × 4 anti-diagonal block matrices: with σ u , u ∈ {x, y, z}, the Pauli matrices.Applying α u on a 4-component orbital, in practice reorders the components, possibly multiplied by a phase factor.
The two-electron energy for the Gaunt operator is thus: where we have introduced the current density Cartesian vector, with components: to rewrite the expression more compactly.The corresponding mean-field, effective oneelectron, direct and exchange operators are: ⃗ j is the trace of the matrix collecting the orbital-pair current densities j pq;u .
The Gaunt direct and exchange operators use the same primitive as the Coulomb operators for the convolution with the inverse-distance kernel.Thus: 1.Although the expressions for the Gaunt mean-field operators appear more complicated than those stemming from the Coulomb interaction, their computational load is only three times higher, because each component of the ⃗ α vector only has four non-zero elements.
2. For each Cartesian component, one can compute a "Gaunt potential" which is then multiplied by the ⃗ α-transformed orbital, exactly as for the Coulomb interaction.
Turning our attention to the gauge two-electron potential, we follow the suggestion of Sun et al. 49 −∇ 1 1 r 12 , and rewrite it as: where the sign/index pairs (−∇ 1 or +∇ 2 ) can be chosen independently for each of the two terms, giving rise to four equivalent expression.
The energy expressions corresponding to each of the above forms can be considerably simplified using integration by parts, thus avoiding the need for differentiating the inversedistance kernel.However, of the four forms presented above, the energy expression obtained by choosing +∇ 2 in both terms of Eq. ( 13) is the most compact and computationally parsimonious: The former three terms are the direct contributions and the latter three the exchange contributions.The use of the inverse-distance kernel is the most significant advantage of this formulation, since that is already an efficient and robust computational primitive in a multiwavelet basis.Note that the calculation of the divergence of the orbital current densities ∇ • ⃗ j pq ≡ ∂j pq;x ∂x + ∂j pq;y ∂y + ∂j pq;z ∂z is both efficient and precise in a multiwavelet basis. 50nally, we present the expressions for the direct and exchange Gauge mean-field operators: All terms in both the direct and exchange operators are applied using the inverse-distance integral operator only.
For completeness, we report also the expressions for the Gauge term when using the inverse-cube-distance form for the operator: The two-electron energy reads: While this is arguably more compact than the sum of all six terms in the previous equation (Eqs.( 14)-( 19)), it has two main disadvantages.First, it is harder to appreciate the physical content of the expression at a glance.Second, it requires the application of a different convolution operator.The latter point is apparent when looking at the expressions for the direct and exchange operators: The new convolution operator, G, is a matrix convolution operator with 6 unique elements, each of which must be implemented by approximating the integral representation of the inverse-cube-distance kernel 51 as a finite exponential sum: 52 Each term, though anisotropic, can be applied in each Cartesian direction separately.Coefficients and exponents in the sum are obtained similarly to those for the inverse-distance convolution operator, see Ref. 47 for details.This form has been tested in our code, but it turned out to be less stable numerically and significantly more demanding computationally.

Computational Details
DIRAC calculations were performed using a nuclear point-charge model and a threshold of 10 −7 on the norm of the error vector (electronic gradient) was chosen as the convergence criterion for the SCF procedure.The chosen basis set for He, Ne 8+ , Ar 16+ , Kr 34+ , Xe 52+ and Rn 84+ was dyall-aug-cvqz. 38,53,54Furthermore, the calculations were performed using default settings for 4-center integral screening and replacing (SS|SS) integrals by a simple Coulombic correction.In our MW implementation it is not possible to perform such a correction, because 4-center integrals do not appear in the formalism.We investigated whether this could impact our perturbative/variational comparisons: with the full two-electron integral tensors the total energy computed with DIRAC changes slightly and computational cost increases significantly.However, the relative error with respect to both our implementation in VAMPyR and in GRASP was practically unaffected.This shows that the error is dominated by the intrinsic limitation of the basis set.

Results and Discussion
We present results for closed-shell, helium-like species: the core 1s-orbitals are doubly occupied and our code explicitly enforces Kramers' time-reversal symmetry (TRS), 55,56 such that the 4-component 1s α is related to 1s β by a quaternionic unitary transformation. 57 a mean-field treatment -e.g.HF and Kohn-Sham DFT -the Coulomb two-electron operator is replaced by the corresponding Direct and Exchange terms, indicated with J and K, respectively.Further inclusion of the Gaunt and Gauge interactions in Eq. ( 1) will result into additional J-and K-like terms.Making use of Kramers TRS has a significant impact on the computational cost: the Coulomb interaction will only encompass the direct term, whereas exchange one will be equal to zero.The Gaunt and Gauge interactions will give rise to both direct and exchange terms but several contributions will either vanish or be identical to each other.
Previous work by Anderson et al. 30 on full 4-component Dirac-Coulomb relativistic calculations used smeared nuclear charge models. 58In particular for the isolated atoms they used the Fermi nuclear model. 58This was done to mitigate numerical issues treating core orbitals with a point-charge model and improve precision.The Fermi model represents the nuclear charge using the Fermi-Dirac distribution for the nuclear charge density, introducing two parameters: the skin thickness and the half-charge radius.The former is set to 2.30 fm (2.30×10 −5 Å) for all nuclei. 58The latter is the radius of a sphere containing half of the total nuclear charge.This parameter depends on the atomic mass of the nucleus M N , with one expression used when M N ≤ 5 atomic mass units and another for M N > 5. 58 The Fermi model for the nuclear charge is smooth and is thus more physically meaningful.Furthermore, it avoids singularities at the nuclei, in contrast to a point-like model.However, the results of Anderson et al. 30 showed that the achieved precision of multiwavelet methods with respect to the grid-based approach available in GRASP decreases with increasing Z, even though a more physically motivated nuclear model was used.
Our multiwavelet implementation in VAMPyR uses two parameters to tune the precision of the calculation: the tolerance, ε, and the polynomial order, k.Furthermore, both pointcharge and Fermi models are available for the nuclei.In order to validate our Dirac-Coulomb Hartree-Fock (DCHF) implementation and reassess the impact of the nuclear model, we performed DCHF calculations with a point-charge model and increasingly tighter precision settings.We report the comparison of our results with GRASP in Figure 1.The relative errors obtained at looser precision settings, as shown in Fig. 1, are not consistent with the user-requested ϵ for heavy elements.The desired precision is user-selected through the settings for ϵ and k and should, in principle, be achieved irrespective of the nuclear model.nel, as shown in by Eqs. ( 14)-( 19), from which it is evident how the magnetic energy term arises as half of the Gaunt term, since both the direct and exchange Gauge contributions (third and sixth terms) contain half of the Gaunt term.For the specific case of 1s 2 systems, the terms involving a gradient in the Gauge energy (i.e., 1st Eq. ( 14), 2nd Eq. ( 15), 4th Eq. ( 17) and 5th Eq. ( 18)) are either zero or cancel each other out, up to the chosen numerical precision ε.Thus, the ratio between the Gauge term (E Gauge ) and the magnetic interaction energy, which corresponds to half of the Gaunt term (E M ag = 1 2 E Gaunt ), should be one (i.e., identical Magnetic and Gauge terms).This was verified comparing the Breit energy corrections from VAMPyR and GRASP results, see Figure 3 and Table 5 in the SI.
The E Gauge /E M ag ratio was calculated previously using Gaussian atomic orbital basis sets for several atoms from Z = 9 to Z = 79. 49It was shown to range between 0.90 (Fluorine) and 0.80 for Z > 56, converging asymptotically.In Table 5 of the SI, where we have considered 1s 2 systems exclusively, we have obtained a unitary ratio between Gauge and magnetic term.
Furthermore, the magnitude of the Gauge term from our results in Table 5 in SI confirms what was previously found by Halbert et al. 62 that in core-electron spectroscopy the Gauge term remains quite significant for the K and L edges, and it must be accounted for, especially for 1s to 2s transitions. 63

Conclusions
We have shown that the 4-component Dirac-Coulomb-Hartree-Fock equations can be solved self-consistently with a fully adaptive MW basis irrespective of the chosen nuclear model, as required with Gaussian basis sets. 64,65e use of multiresolution analysis (MRA) with a MW basis to solve the KS-DFT equations allows to separate model errors from discretization (i.e.basis set) errors, with the latter precisely quantifiable.Thus, the use of a MW basis provides fundamental insight to understand the range of applicability of KS-DFT with localized basis sets.This issue is especially relevant for 4-component relativistic calculations on heavy elements where the description of the core electrons is challenging due the nature of the Dirac equation combined with the extremely high nuclear charge and a reduced availability of GTO bases.
We have shown that the DCHF ground state combined with the Breit Hamiltonian as a perturbative correction can reproduce grid-based calculations performed with GRASP .Albeit not performed in this work, the fully variational inclusion of the Gaunt and Gauge terms can be obtained by making use of the corresponding operator expressions (Equations (11a) and (11b) and Equations (20a) and (20b) for Gaunt and Gauge, respectively).This has not been done for the current work both to simplify the comparison with the GRASP code and because of excessive memory demands of the current pilot implementation.The latter is indeed the main challenge for future extensions to general molecular systems where the simplifications that enabled our results (time-reversal symmetry, spherical symmetry of the 1s orbital) will no longer hold.Work is in progress in our group to overcome these hurdles.
The unitary E Gauge /E M ag ratio for s-orbitals explains how neither Gaunt nor Gauge terms can be neglected for core-electron spectroscopy and explains the importance of considering both these terms when x-ray photoelectron spectra are calculated to fit the experimental ones. 62,63,66Our results confirm the validity of the MW approach for future development of core-electron spectroscopy to resolve the structures of oxides of f -elements and other strongly correlated systems.

Figure 2 :
Figure 2: Comparison between the spinorbit energies (a) and Gaunt terms (b) coming from VAMPyR and DIRAC for selected systems in electronic configuration 1s 2 .The y axis shows the logarithm of the unsigned relative difference between VAMPyR and DIRAC results.The VAMPyR calculations were done with Legendre polynomial order k = 10 and tolerance ϵ = 10 −8 .All codes have used a nuclear point charge model as described in Ref. 58.

Figure 3 :
Figure 3: Comparison between the Breit perturbative corrections computed VAMPyR and GRASP for noble gases and actinides in electronic configuration 1s 2 .The y axis shows the logarithm of the unsigned relative difference between VAMPyR and GRASP results.The VAMPyR calculations were done with Legendre polynomial order k = 10 and tolerance ϵ = 10 −8 .All codes have used a nuclear point charge model as described in Ref. 58.