Effective spin-orbit models using correlated first-principles wave functions

Diffusion Monte Carlo is one of the most accurate scalable many-body methods for solid state systems. However, to date, spin-orbit interactions have not been incorporated into these calcualtions at a first-principles level; only having been applied to small systems. In this technique, we use explicitly correlated first-principles quantum Monte Carlo calculations to derive an effective spinorbit model Hamiltonian. The simplified model Hamiltonian is then solved to obtain the energetics of the system. To demonstrate this method, benchmark studies are performed in main-group atoms and monolayer tungsten disulfide, where high accuracy is obtained.

The interplay between the electron correlation and relativistic effects may give rise to a plethora of intriguing quantum phases in condensed matter systems, including the proposed axion insulating state, topological Mott insulating state, Weyl semimetal and quantum spin liquid, etc. [1,2]. Spin-orbit interactions in materials have been widely investigated in the weakly interacting regime using mean-field theory, which has resulted in the prediction of topological phases [3][4][5]. The spin-orbit effect can remove orbital degeneracies and modify the electronic and magnetic structures of materials; it thus finds applications in spintronic devices [6]. In the strong-tointermediate correlation regime, spin-orbit interaction is proposed to play important roles in unconventional Mott insulating and topological semimetal states [7,8], quantum spin liquid [9,10], axion insulating state [11] and multipolar orders [12,13].
It is a major challenge to describe both electron correlation and spin-orbit effects in materials. Many-body wave function techniques are the most direct way of accessing electron correlations with high accuracy. In particular, quantum Monte Carlo techniques have seen success [14][15][16][17] in computing the properties of materials with strong correlation. However, to our knowledge, there are no first-principles quantum Monte Carlo calculations on bulk materials that also include spin-orbit effects, previous to this research.
Spin-orbit effects combined with electron correlation are commonly considered for small molecular systems. Expansion of the spin-orbit interaction operator in terms of configuration state functions suffer from a slow convergence of the electron correlation contribution [18]. As for methods based on stochastic sampling procedures, a recent work employs both an auxiliary continuous parameter to represent the electron spin in the wave function and a modified sampling process in the quantum Monte Carlo (QMC) calculation [19][20][21]. This method is shown to obtain accurate results for atoms and few-atom molecules, two-dimensional homogeneous electron gas, quantum wells, and circular quantum dots with Rashba interaction [22][23][24]. However, since the spin-orbit energy and the total Coulomb energy are treated on the same footing (a positive feature for some reasons), the compu-tational cost needed to resolve the tiny energy differences among different spin configurations can be very high for solid state systems.
In this paper, we demonstrate a method for treating spin-orbit interactions using quantum Monte Carlo in a scalable and accurate way through rigorous model Hamiltonians. We implement spin-orbit interaction calculations using samples of static-spin wave functions, then downfold this set of first-principles data to a model Hamiltonian. On this manifold, the computational cost needed to resolve different spin configurations is greatly reduced compared to dynamical spins in the calculation, very similar to that needed in a standard fixed-node diffusion Monte Carlo (FN-DMC) calculation. To demonstrate our method, we perform benchmark studies in small atomic systems and a solid state system. Using this technique, we calculate the ground state fine-structure splittings induced by spin-orbit interaction for maingroup atoms. Our calculations generate values agree with the experiments for most atoms, with errors attributed mainly to the lack of jj coupling in the calculation. To show that this method can be easily generalized to realistic materials, we compute the band splitting at K point for monolayer tungsten disulfide (WS 2 ), yielding a value of 0.39(2) eV, as compared to the photoluminescence measurement of 0.4 eV [25].
Effective spin-orbit Hamiltonian. Our method is quite general and can treat effective interactions as well as multiple bands [26,27]. However, for this paper, we will focus on the case of electronic states that are degenerate in the absence of spin-orbit interaction; for example the six occupations of a single electron on a p manifold. Consider many-body wave function |Ψ N within that degenerate subspace. Then the expectation value of the energy is given by where we evaluate E SOI by subtracting an averagedrelativistic effective core potential (ARECP) from the relativistic effective core potential (RECP) [21]: where operatorsv iI l,j=l± 1 2 are the spin-up/down radial parts of the two-component RECP between the ith electron and the Ith ion, respectively. The spin-orbit energy in Eq. 1 is given by We take the effective Hamiltonian of the degenerate manifold to be where the operator L iI = (r i − r I ) × 1 i ∇ i is the orbital angular momentum of the ith electron relative to Ith ion, and S i is the spin angular momentum operator acting on the ith electron.
Using the results proved in [26,27], the effective λ can be estimated using the first-principles wave functions |Ψ n by writing Using Eq. 1, we thus have This relationship is true for any wave function |Ψ n sampled from the degenerate manifold if the spin-orbit model is valid. Since both E SOI [Ψ n ] and Ψ n | iI L iI · S i |Ψ n can be computed using first-principles wave functions, the coefficient λ can be obtained by fitting to the linear relationship between these two values. For periodic systems, the effective Hamiltonian at momentum k is given by H eff k,SOI = iI λ I (k) L iI · S i (for the monolayer WS 2 , the spin-orbit interaction is mostly contributed by W atoms). This method can be easily extended to the case where the wave functions do not belong to a degenerate manifold by including parameters that represent orbital energies, interaction terms, and so on. This has been performed for non-spin-orbit-coupled systems [27].
Details of the calculations. For a condensed matter system, we first construct single Slater determinants from the single-particle orbitals computed using Hartree-Fock or density functional theory. These Slater determinants are sampled such that they differ only by single-particle orbitals that have the same scalar-relativistic energy. We construct the following Slater-Jastrow trial wave functions, where the spins are collinear, chosen to be along the z axis. Determinants D ↑ and D ↓ are constructed from occupying the lowest energy spin-up and spin-down singleparticle orbitals, respectively. To incorporate the electron correlation explicitly, the trial wave functions are constructed by multiplying these Slater determinants by Jastrow factors U s which include up to three-body interactions (see [28]). The Jastrow factors are then optimized by minimizing the variance of the trial scalarrelativistic energy. Finally, we apply FN-DMC to project out the lowest energy wave function that has the same nodal structure as the corresponding trial wave functions. The scalar-relativistic Hamiltonian we use to evolve the trial wave function in FN-DMC is given by expressing the Coulomb potential using the sum of electron-electron interaction among the valence electrons and the averagedrelativistic effective core potential (ARECP). T -moves are used to evaluate the nonlocal effective core potentials during the imaginary-time evolution of the trial wave functions [29]. In order to correct the time-step error introduced by Trotter product, we extrapolate the result to zero-time-step limit using a method based on the Bayesian probability distribution described in [30]. All VMC and FN-DMC calculations are performed using the open-source quantum Monte Carlo code QWalk [28]. Because the ab initio spin-orbit interaction Hamiltonian (Eq. 3) does not commute with the scalar-relativistic Hamiltonian, we correct the mixed estimator error up to first order by computing the extrapolated estimators given in [31]. We compute the extrapolated estimator values of the spin-orbit interaction energies and the corresponding descriptor iI L iI · S i for several fixed-node wave functions, then find the coefficient λ using linear regression. This procedure maps a set of values calculated from first-principles wave functions to a model Hamiltonian Eq. 3. Such procedures can be generalized to expectation values calculated using dynamic-spin wave functions, but for the brenchmarking systems considered here, static-spin wave functions are sufficient.
For the atomic systems, we use GAMESS [32,33] to generate the single-particle orbitals for the cations and anions (Ga + , In + , Ge 2+ , Sn 2+ , Br − and I − ), such that the np x , np y and np z orbitals have exactly the same energies for each ion. Thus, we can reduce the orbital relaxation effect introduced by performing restricted openshell Hartree-Fock calculations directly on atomic systems, Then, we use these ionic orbitals to construct the Slater determinants for the corresponding atoms. We perform the calculations using the Dolg-Stoll and the Trail-Needs ECPs and the available basis sets provided in [34][35][36][37][38][39]. The Dolg-Stoll ECPs include relativistic effects implicitly by adjusting the orbitals at the multiconfiguration Dirac-Hartree-Fock level and including the Breit interaction [34][35][36][37]. We also compare the large-core and small-core Dolg-Stoll ECPs, which exclude or include respectively the semi-core (n − 1)spd shells in the valence shell. The Trail-Needs ECPs are large-core only for the elements considered here, and are finite at the origin, thus are particularly suitable for QMC calculations [38,39].
For the WS 2 calculations, we use CRYSTAL17 [40,41] to perform mean-field calculations for monolayer WS 2 , with the lattice constant a = 3.23Å, and the spacing between two S atoms in the same unit cell d = 3.19Å. The molecular orbitals are calculated for the primitive cell by solving the Kohn-Sham equation, where the exchangecorrelation interactions are included within the generalized gradient approximation in the Perdew-Burker-Ernzerhof form [42]. We use a converged 12 × 12 shrinking factor for reciprocal lattice vectors. The Dolg-Stoll averaged-relativistic ECP (60 core electrons for W and 10 core electrons for S) is used [43,44]. A truncated Gaussian basis set including up to triple ζ functions is used to expand the molecular orbitals for W [43]. For S, we use the starting basis set provided in reference [44].
Benchmark studies. Taking the Ga atom as an example, the ground-state electron configuration without spinorbit interaction is [Ar]3d 10 4s 2 4p 1 with L = 1, S = 1 2 ( 2 P , 6-fold degenerate). The spin-orbit interaction splits the 6-fold degenerate state into two manifolds, the J = 3 2 quartet and the J = 1 2 doublet (see Fig. 2, panel (a)). In the Ga atom, the sampling of wave functions on a sub-Hilbert space with the lowest scalar-relativistic energy is performed by projecting out the fixed-node wave functions from the trial Slater-Jastrow wave functions, whose Slater parts are linear superpositions of the Slater determinants including 4p x , 4p y and 4p z orbitals. In our calculations, we generate 20 samples of such trial wave functions with different SOI descriptors values iI L iI · S i (a). Ga, In ranging from −0.5 to 0.5 in atomic units, when taking the dimensionless S z for spin-up/down electron to be ±1/2. Within this context, the spin-orbit splitting ∆ of Ga atom ground-state configuration is given by ∆ = 3/2λ.
In Fig. 1, we plot the SOI energies E SOI versus the descriptors iI L iI · S i calculated on 20 Slater determinants using variational Monte Carlo (VMC) (labeled by the red triangles), and on the corresponding fixednode wave functions for Ga atom (labeled by the blue circles). The Dolg-Stoll small-core (10 core electrons) and (12s12p9d)/[6s6p4d] basis set are used [34]. There is a linear relationship between these two quantities, as required by Eq. 4.
To demonstrate the accuracy of the technique, we compute the spin-orbit splittings of the ground state configurations of the following main-group atoms: Ga, In, Ge, Sn, Br, and I. The In, Br and I atom have similar spin-orbit splitting pattern as Ga, while the ground state splittings of Ge and Sn atoms ∆ J=1 and ∆ J=2 satisfy ∆ J=1 = 1/2 |λ| and ∆ J=2 = 3/2 |λ| (see Fig. 2 (a)).
To check if there is dependence on the ECP, we perform the same calculations using the Dolg-Stoll and the Trail-Needs ECPs [34][35][36][37][38][39]. Fig. 2 (b) shows the FN-DMC extrapolated estimator values of the SOI splittings of the ground state configurations, versus the experimentally determined fine-structure splittings for the chosen main-group atoms (extracted from NIST atomic spectra database [45]). The light gray line indicates an equality between the calculated values and the experimentally determined values. From the figure, we can see that our method yields better agreement with the experiments when using the Dolg-Stoll small-core ECP and the Trail-Needs ECP, compared with that when using the Dolg-Stoll large-core ECP.
We attribute the discrepancy between the calculated SOI splitting and the experimentally determined finestructure splitting of the Sn atom J = 1 state to the interplay of the jj-coupling and the LS-coupling (the spin-orbit interaction). Our model Hamiltonian Eq. 3 only includes the spin-orbit interaction, thus the Landé interval rule follows, i.e., ∆ J=2 /∆ J=1 = 3/1. However, due to the non-negligible effect of the jj-coupling, the experimentally observed fine-structure splittings of Sn atom is not entirely induced by spin-orbit interaction of electrons. As a consequence, Landé interval rule does not apply, i.e., ∆ J=2 /∆ J=1 = 3/1.
In order to demonstrate that our model can be easily generalized to larger scale systems, we perform a benchmark study in monolayer WS 2 . Due to the presence of tungsten atoms, the strong spin-orbit interaction splits the valence band maximum at K. We apply our method to compute this band splitting. We construct trial wave functions from superpositions of four Slater determinants multiplied by three-body Jastrow factors. These Slater determinants are constructed by promoting one spin-up/down electron from the two-fold degenerate valence band maximum (A ) to the conduction band minimum (A 1 ) at K, such that they have degenerate scalar-relativistic energy (see Fig. 3). Fig. 4 shows the FN-DMC extrapolated SOI energies versus the SOI descriptors calculated using samples of many-body wave functions of monolayer WS 2 at K point. The slope of the fitted line is the spin-orbit interaction strength λ in our linear model. The calculated spin-orbit splitting ∆ for the valence band maximum at K is 0.39 (2) eV. This is statistically in agreement with the band splitting 0.43 eV calculated using PBE exchange-correlation functional reported in reference [47]. It also agrees with the reported band splitting 0.4 eV, determined by measuring the differential reflectance spectra of exfoliated single layer 2H-WS 2 [25]. This demonstrates that our method can be easily generalized to perform spin-orbit interaction calculations in solid state systems.
Conclusion. We have demonstrated a scalable firstprinciples quantum Monte Carlo technique for extended systems; this will allow future studies to blend treatment of electron correlation and spin-orbit interactions using QMC calculations. The method is based on de- Band structure of monolayer WS2 calculated on the unit cell, using the Perdew-Burke-Ernzerhof exchangecorrelation functional [42] without spin-orbit interaction. The red dashed line represents the Fermi level. The conduction band minimum and valence band maximum at K are labeled by A 1 and A . The corresponding single-particle orbitals are plotted on the right using the Visual Molecular Dynamics software [46], with orbital shape derived from the isosurface of the magnitude, colored by the phase. riving an effective Hamiltonian for the spin-orbit interaction. In this method, the major computational cost is contributed by the evaluation of spin-orbit interaction energy on the FN-DMC wave functions, less than twice that of a standard FN-DMC calculation. One can include also electron-electron interactions in the effective Hamiltonian.
We demonstrated this technique in atomic systems and monolayer WS 2 . For the main-group atoms, we compute the spin-orbit splittings of the ground state configurations and obtain results in agreement with the experimentally determined fine-structure splittings. For the monolayer WS 2 , we compute the band splitting induced by spin-orbit interaction at K. Our first-principles result agrees with previous calculations and experiments, thus demonstrating the ability of this method to be generalized to larger scale systems. A major promise of this technique is that one can treat electron-electron interactions, spin-orbit effects, and onebody terms in effective Hamiltonians all on the same footing. As mentioned before, the cost is similar to a standard FN-DMC calculation, which will allow it to be applied to realistic models of materials. We envision this opening a new frontier for modeling spin-orbit effects in correlated materials.
This work was funded by the grant DOE FG02-12ER46875 (SciDAC). The authors want to thank Cody Melton for precious discussions.