Generalized Kohn-Sham Approach for the Electronic Band Structure of Spin-Orbit Coupled Materials (2024)

Jacques K. Desmaraisjacqueskontak.desmarais@unito.itDipartimento di Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy Giacomo AmbrogioDipartimento di Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy Giovanni VignaleInstitute for Functional Intelligent Materials, National University of Singapore, 4 Science Drive 2, Singapore 117544 Alessandro ErbaDipartimento di Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy Stefano Pittalisstefano.pittalis@nano.cnr.itIstituto Nanoscienze, Consiglio Nazionale delle Ricerche, Via Campi 213A, I-41125 Modena, Italy

(January 1, 2024)

Abstract

Spin-current density functional theory (SCDFT) is a formally exact framework designed to handle the treatment of interacting many-electron systems including spin-orbit coupling (SOC) at the level of the Pauli equation. In practice, robust and accurate calculations of the electronic structure of these systems call for functional approximations that depend not only on the densities and currentsbut also on spinors explicitly.Here we extend the generalized Kohn-Sham (GKS) approach of [Seidl et all. “Generalized Kohn-Sham schemes and the band-gap problem”, Phys. Rev. B 53, 3764 (1996)] to SCDFT.This framework entails the prominent cases of hybrid forms and meta-generalized-gradient-approximations.We clarify that the exchange-correlation potentials conjugate to the currentsneed to be computed within the GKS approach only when the spin currents are included in the functional form explicitly.We analyze the consequence of this fact for various approximations and numerical procedures for the evaluation of SOC effects.The practical power of the extended approach is demonstrated by calculating the spin-orbit induced/enhanced band-splittings of inversion-asymmetric single-layer MoSe ${}_{2}$ and inversion-symmetric bulk $\alpha$ -MoTe ${}_{2}$ . Key to these results is the capacity to account for SOC self-consistently while employing energy functionals and effective potentials that depend (implicitly or explicitly) on spin currents.

pacs:

71.15.Mb, 71.15Rf, 31.15.E-

I Introduction

Since the early days of quantum mechanics, spin-orbit interactions have played a central role in our understanding of the electronic properties of atoms, molecules, and solids. The Dirac equation and its simplified version,Greiner (2012) the two-component Pauli equation, were pinnacle achievements of that era, leading to a unified description of fine structure and Zeeman splittings in practically all systems known at the time, including those which would later turn out to becrucial to the semiconductor revolution (e.g., Ge, Si, and GaAs). Yu and Cardona (2010)

In this century, the discovery of nontrivial topological properties of the band structure of periodic solids, such as topological insulators and Weyl semimetals, whose extraordinary properties include quantized transport coefficients, magnetoelectric response, chiral anomalies, non-reciprocity, etc. has led to an explosion of interest in spin-orbit interactions. Indeed, by making the electronic wave functions complex even in the absence of a magnetic field, the spin-orbit interaction sets the scenario for the nontrivial response properties, band inversions and topological quantum numbers that underlie the above mentioned effects. Bernevig and Hughes (2013); Vanderbilt (2018) In a parallel development, the emergence of spintronics has raised the interest in non-collinear spin textures both in real and in momentum space.Žutićetal. (2004) For instance, the phenomenon of spin-momentum locking – the emergence of a spin texture in momentum space – is responsible for remarkable magneto-transport effects, such as the unidirectional magnetoresistance.Zhang and Vignale (2016)

In this context, it has become more pressing than ever to develop the tools of computational electronic structure so that they can be trusted to quantitatively predict the impact of spin-orbit interactions on properties such as spin-orbit splitting of bands, closing and re-openings of gaps at topological phase transitions, the positions of conical intersection (Dirac and Weyl points) in the Brillouin zone, the dispersion of Fermi arcs and the shape of non-collinear spin textures.

Given the dominance of density functional theory (DFT) on the landscape of computational electronic structure, it seems natural to seek to include spin-orbit interaction effects through an extension of the DFT framework. What we mean by this is much more than simply including the spin-orbit interaction as an additional one-body term in the Kohn-Sham equation of DFT – an option that is already incorporated and widely available in existing electronic structure packages. Rather, in order to achieve quantitative accuracy and predictive power, we believe it is essential to include the effect of the spin-orbit interaction in the many-body potentials that appear in a (suitably generalized) Kohn-Sham theory.

The formal frameworkfor implementing this program has been known for a long time: it is the U(1) $\times$ SU(2)-invariant Spin Current DFT (SCDFT), see Refs. Vignale and Rasolt, 1988; Bencheikh, 2003.This theory includes 16 external fields coupling to 16 densities, i.e., the scalar potential coupling to the particle density, the Zeeman magnetic field (three components) coupling only to the spin density, the charge vector potential (3 components) coupling only to the orbital current density, and, lastly, the SU(2) vector potential (a $3\times 3$ tensor) coupling to the spin current densities. Depending on which functional form is invoked in the calculations, there are several convenient ways of organizing this extended set of densities and potentials; see for example Refs. Rohra and Görling, 2006; Heaton-Burgessetal., 2007; Abedinpour etal., 2010.Because it deals in a unified fashion withthe magnetic interactions and the spin-orbit coupling, SCDFT appears to be the ideal framework to simulate a multitude of materials useful for magnetism, spintronics, orbitronics, valleytronics, and topologically non-trivial states as described above.

In spite of its great promise, SCDFT has so far lagged behind other DFT and non-DFT methods in its application to real material. The reason for this delay can be traced to the lack of good and transferable approximations for the exchange-correlation (xc) energy functional in terms of spin-current densities.In the last two decades, it has become increasingly evident that accurate calculations of the electronic structure, including in particular band gaps and band splittings require functionals that depend on the densities not only explicitly (as in the traditional formulation of SCDFT) but also implicitly, through single-particle spin orbitals.Kurth and Pittalis (2006); Pittalis etal. (2006); Rohra and Görling (2006); Heaton-Burgessetal. (2007); Sharma etal. (2007a, b); Trushin and Görling (2018); Ullrich (2018); Pluhar and Ullrich (2019)The emergence of orbital-dependent functionals began with the widespread practice of including exact exchange, or a fraction thereof, in the energy functional (the so-called “hybrid” functionals), and gained momentum with the development of “meta-GGA” functionals, in which the traditional set of densities is augmented by the inclusion of the (spin-)kinetic-energy density. Most importantly for SCDFT, it was realized that spin-orbital-dependent functionals are explicitly required in any non-trivial gauge-invariant formulation.Pittalis etal. (2017)

The problem with orbital functionals is that, because they are regarded as implicit nonlocal functionals of the density, they must be differentiated with respect to the densities in order to yield the Kohn-Sham potentials. This differentiation is difficult, as it involves the functional derivative of the orbitals with respect to the densities. The procedure is usually referred to asthe “Effective Potential Method” (OPM), and, whilethe resulting “Optimized Effective Potential” (OEP) is a legitimate local Kohn-Sham potential, the benefits of locality arewiped out by the complexity and costliness of the numerical treatment.WU and YANG (2003); Gidopoulos and Lathiotakis (2012); Eich and Hellgren (2014); Gidopoulos and Lathiotakis (2013); Oueis etal. (2023)

Experience in regular (Spin-)DFT has demonstrated that the cost of implementing spin-orbital functionals can be lowered, and the numerical treatment simplified by switching to a generalized Kohn-Sham (GKS) framework,Becke (1993); Seidl etal. (1996); Perdew etal. (2017)which admits the use of non-local effective potentials, as naturally appear in hybrid functional forms.In fact, the key ideas of this approach are also used in calculations involving meta-GGA functionals. Neumann etal. (1996); Eich and Hellgren (2014); Lehtola etal. (2020)In this method, as in the OPM, the xc potential is expressed as the sum of two parts: a functional of the spin-orbitals – typically, but not necessarily, a fraction of the exact exchange – plus a regular explicit functional of local densities. This simple shift in perspective has far-reaching consequences. The functional derivative of the explicitly orbital-dependent part of the functional yields a nonlocal, but simple potential – the Fock potential in the case of exact exchange – while the functional derivative of the regular part yields a local potential as in the standard Kohn-Sham formalism. The resulting GKS equation combines the accuracy of exact nonlocal exchange with the flexibility of semilocal density functional approximations for the correlation energy. Crucially, band gaps calculated in this manner become more closely related to the KS gaps, Perdew etal. (2017); Wing etal. (2021)since part of the derivative discontinuity of the exact functional is captured by the discontinuous dependence of the orbitals on band occupation. The rigorous theoretical foundation of the GKS approach is presented in Ref.Seidl etal., 1996 for DFT. Here, we extend this approach to SCDFT.

In doing so, we lay down the framework that allows us to merge two previous works from some of the same authors of the present work:the implementationDesmaraisetal. (2020a, b, 2021a); Desmarais etal. (2019) and applicationComaskey etal. (2022); Bodo etal. (2022) of the (regular) global hybrids in SCDFT via the Crystal codeand the (formal) proposal of Meta-GGAs for SCDFT.Pittalis etal. (2017)This allows us to demonstrate that at the heart of the success of the method is its ability to include the dependence of the effective many-body potentials on spin currents. Even when this dependence is only implicit (as discussed below) its inclusion is essential to obtain agreement with experimental results — whereas conventional treatments of spin-orbit coupling fail. When the inclusion of spin currents is explicit, it becomes necessary to include the feedback on the effective-vector potentials explicitly as well.

The paper is organized as follows: We start with an introductory section on the (regular) Kohn-Sham approach to SCDFT. We then describe the GKS approach and proceed to its application.We discuss first the prominent case of exact exchange followed by an in-depth discussion of global hybrid forms.We also discuss the case of spin-current dependent Meta-GGAs. In this way, we are able to highlight several crucial features which are peculiar to the GKS approach of SCDFT.We apply the approach to the calculations of valence-band splittings induced/enhanced by SOC in inversion-asymmetric single-layer, 2D, MoSe ${}_{2}$ with spin-splitting (Rashba-I effect) and inversion-symmetric bulk $\alpha$ -MoTe ${}_{2}$ with spin-valley locking (Rashba-II effect).After touching on near-future developments, we conclude.

II Formal Aspects

II.1 Spin-Current Density Functional Theory

In order to appreciate the key difference between SCDFT and Spin-DFT (SDFT, the most popular flavor of DFT), it is useful to start with the SDFT Hamiltonian:

\displaystyle\hat{H}_{\rm{\tiny\tiny SDFT}}=\frac{1}{2}\int d^{3}r~{}\hat{\Psi%}^{\dagger}({\bf r})\left(-i\nabla\right)^{2}\hat{\Psi}({\bf r})+\int d^{3}r~{%}\left[\hat{n}({\bf r})v({\bf r})+\hat{m}^{a}({\bf r}){B}^{a}({\bf r})\right]+%\hat{W}\;,

(1)

where $\hat{\Psi}^{\dagger}({\bf r})=(\hat{\Psi}^{\dagger}_{\uparrow}({\bf r}),\hat{%\Psi}^{\dagger}_{\downarrow}({\bf r}))$ denotes a two-component creation field operator ( $\uparrow$ and $\downarrow$ refer to the spin “up” and “down”); $v$ represent an external scalar-multiplicative potentialthat couples to the electrons via the particle-density operator $\hat{n}=\hat{\Psi}^{\dagger}\hat{\Psi}$ ; ${B}^{a}$ is the $a$ compoenent of a magnetic field that couples to the electrons via the spin-density operator $\hat{m}^{a}=\hat{\Psi}^{\dagger}{\sigma}^{a}\hat{\Psi}$ ( ${\sigma}^{a}$ denotes the Pauli matrices $\sigma^{x},\sigma^{y},\sigma^{z}$ ).The last term of the right hand side denotes the electron-electron interaction

\displaystyle\hat{W}

\displaystyle=\int d^{3}r\int d^{3}r^{\prime}\frac{\hat{\Psi}^{\dagger}({\bf r%})\hat{\Psi}^{\dagger}({\bf r}^{\prime})\hat{\Psi}({\bf r}^{\prime})\hat{\Psi}%({\bf r})}{2|{\bf r}-{\bf r}^{\prime}|}\;.

(2)

Unless otherwise stated, we use Hartree atomic units in which $\hbar=m=1$ .Also note that, for notational convenience, the Bohr magneton $\mu_{B}$ factor is absorbed in the symbol ${B}^{a}$ .

The SDFT Hamiltonian does not include either the vector potential corresponding to ${B}^{a}$ nor spin-orbit couplings. For taking into account the latter interactions in DFT fashhion,we consider the SCDFT Hamiltonian, which can be obtained via the“minimal” substitution $-i\nabla\rightarrow-i\nabla+\frac{1}{c}{\bf A}({\bf r})+\frac{1}{c}\sigma^{a}{%\bf A}^{a}({\bf r})$ :

\displaystyle\hat{H}_{\rm{\tiny\tiny SCDFT}}=\frac{1}{2}\int d^{3}r~{}\hat{%\Psi}^{\dagger}({\bf r})\left[-i\nabla+\frac{1}{c}{\bf A}({\bf r})+\frac{1}{c}%\sigma^{a}{\bf A}^{a}({\bf r})\right]^{2}\hat{\Psi}({\bf r})+\int d^{3}r~{}%\left[\hat{n}({\bf r})v({\bf r})+\hat{m}^{a}({\bf r}){B}^{a}({\bf r})\right]+%\hat{W}\;.

(3)

Above and in the following, we denote with bold characters, $\boldsymbol{A}$ , quantities with spatial indices (Greek lower indices, $A_{\mu}$ , when written explicitly); and use an arrow, $\vec{A}$ , to denote quantities with spin indices (upper Latin indices, $A^{a}$ , when written explicitly).Thus ${\vec{\bf A}}$ denotes a tensor with two indices $A^{a}_{\mu}$ . Both $\mu$ and $a$ have values $x$ , $y$ , $z$ . Contractions over spatial indices are denoted with “ $\cdot$ ”, Einstein convention is applied to intend summation over repeated indices.

Eq. (3) includes besides the terms of the SDFT Hamiltonian also a(charge-) vector potential ${\bf A}({\bf r})$ and a spin-vector potential ${\bf A}^{a}({\bf r})$ .While ${\bf A}({\bf r})$ is useful to represent an external magnetic field in the usual fashion $B^{a}=\epsilon^{a\mu\nu}\partial_{\mu}A_{\nu}$ , ${\bf A}^{a}({\bf r})$ is useful to represent the (one-body) spin-orbit couplings in the system.Fröhlich and Studer (1993); Bencheikh (2003); Abedinpour etal. (2010)Note, in our notation we absorb the prefractor $\frac{\mu_{B}}{2}$ in the symbol ${\bf A}^{a}$ .These vector potentials may be viewed as some “induction” fields, in the sense that they induceparticle and spin currents in the systems on which they act.To make the latter fact apparent and towards a proper density functionalization, let us expand the first term in Eq. (3)

\displaystyle\hat{H}_{\rm{\tiny\tiny SCDFT}}=\hat{T}+\hat{W}+\int d^{3}r~{}%\hat{n}({\bf r})\tilde{v}({\bf r})+\int d^{3}r~{}\hat{m}^{a}({\bf r})\tilde{B}%^{a}({\bf r})+\frac{1}{c}\int d^{3}r~{}{\hat{\bf j}}({\bf r})\cdot{\bf A}({\bfr%})+\frac{1}{c}\int d^{3}r~{}{\hat{\bf J}^{a}}({\bf r})\cdot{{\bf A}^{a}}({\bf r%})\;,

(4)

where

\displaystyle\hat{T}=\int d^{3}r\;\hat{\Psi}^{\dagger}({\bf r})\left(-\frac{%\nabla^{2}}{2}\right)\hat{\Psi}({\bf r})\;,

(5)

\tilde{v}=v+\frac{1}{2c^{2}}\left[{\bf A}\cdot{\bf A}+{\bf A}^{a}\cdot{\bf A}^%{a}\right]\;,

(6)

and

\tilde{B}^{a}=B^{a}+\frac{1}{2c^{2}}{\bf A}\cdot{\bf A}^{a}\;.

(7)

Crucially, in Eq. (4) we also find the (paramagnetic) particle-current operator $\hat{{\bf j}}=\frac{1}{2i}\left[\hat{\Psi}^{\dagger}\nabla\hat{\Psi}-\left(%\nabla\hat{\Psi}^{\dagger}\right)\hat{\Psi}\right]$ and the (paramagnetic) spin-current operator $\hat{{\bf J}}^{a}=\frac{1}{2i}\left[\hat{\Psi}^{\dagger}{\sigma}^{a}\nabla\hat%{\Psi}-\left(\nabla\hat{\Psi}^{\dagger}\right){\sigma}^{a}\hat{\Psi}\right]\;.$ Therefore we may anticipate that while SDFT may only accounts for particle- and spin-density self-consistently,SCDFT must account also for particle- and spin-current self-consistently (below).

In fact, given the external fields $v$ , ${B}^{a}$ , ${\bf A}$ , and ${\bf A}^{a}$ , the ground-state energy may then be determined by means of a constrained-search minimization:Levy (1982); Lieb (1983)

\displaystyle E=\min_{(n,~{}{\vec{m}},~{}{\bf j},~{}{\vec{\bf J}})}\Big{\{}F[n%,{\vec{m}},{\bf j},{\vec{\bf J}}]+\int d^{3}r~{}{n}({\bf r})\tilde{v}({\bf r})%+\int d^{3}r~{}{m}^{a}({\bf r})\tilde{B}^{a}({\bf r})+\frac{1}{c}\int d^{3}r~{%}{{\bf j}}({\bf r})\cdot{\bf A}({\bf r})+\frac{1}{c}\int d^{3}r~{}{{\bf J}^{a}%}({\bf r})\cdot{{\bf A}^{a}}({\bf r})\Big{\}}\;,

(8)

with

II.2 From Regular to Generalized-KS Systems in SCDFT

GKS systems can be introduced in SCDFT in a way that is similar to (S)DFT by noting thatthe minimization in Eq. (8) can equivalently be performed by invoking differentsplittings of $F[n,{m}^{a},{\bf j},{\bf J}^{a}]$ and a different minimization procedure. In detail, let us consider:

\displaystyle F[n,{\vec{m}},{\bf j},{\vec{\bf J}}]=F_{\rm{\tiny\tiny GKS}}[n,{%\vec{m}},{\bf j},{\vec{\bf J}}]+E^{\rm{\tiny\tiny GKS}}_{\rm Hxc}[n,{\vec{m}},%{\bf j},{\vec{\bf J}}]\;,

(21)

where

\displaystyle F_{\rm{\tiny\tiny GKS}}[n,{\vec{m}},{\bf j},{\vec{\bf J}}]=\min_%{\Phi\rightarrow(n,{\vec{m}},{\bf j},{\vec{\bf J}})}\langle\Phi|\hat{O}_{\rm{%\tiny\tiny GKS}}|\Phi\rangle

(22)

is the analogous of Eq. (10) but here $\hat{O}_{\rm{\tiny\tiny GKS}}$ may differ from $\hat{T}$ by including some interaction (below).Next, note that

$\displaystyle E$	$\displaystyle=\min_{(n,{\vec{m}},{\bf j},{\vec{\bf J}})}~{}\Big{\{}\min_{\Phi%\rightarrow(n,{\vec{m}},{\bf j},{\vec{\bf J}})}\langle\Phi\|\hat{O}_{\rm{\tiny%\tiny GKS}}\|\Phi\rangle+E^{\rm{\tiny\tiny GKS}}_{\rm Hxc}\Big{[}n,{\vec{m}},{%\bf j},{\vec{\bf J}}\Big{]}$
	$\displaystyle+\int d^{3}r~{}{n}({\bf r})\tilde{v}({\bf r})+\int d^{3}r~{}{m}^{%a}({\bf r})\tilde{B}^{a}({\bf r})+\frac{1}{c}\int d^{3}r~{}{{\bf j}}({\bf r})%\cdot{\bf A}({\bf r})+\frac{1}{c}\int d^{3}r~{}{{\bf J}^{a}}({\bf r})\cdot{{%\bf A}^{a}}({\bf r})\Big{\}}$
	$\displaystyle=\min_{\Phi}~{}\Big{\{}\langle\Phi\|\hat{O}_{\rm{\tiny\tiny GKS}}\|%\Phi\rangle+E^{\rm{\tiny\tiny GKS}}_{\rm Hxc}\Big{[}n[\Phi],{{\vec{m}}}[\Phi],%{\bf j}[\Phi],{\vec{\bf J}}[\Phi]\Big{]}$
	$\displaystyle+\int d^{3}r~{}{n}[\Phi]({\bf r})\tilde{v}({\bf r})+\int d^{3}r~{%}{m}^{a}[\Phi]({\bf r})\tilde{B}^{a}({\bf r})+\frac{1}{c}\int d^{3}r~{}{{\bf j%}}[\Phi]({\bf r})\cdot{\bf A}({\bf r})+\frac{1}{c}\int d^{3}r~{}{{\bf J}^{a}}[%\Phi]({\bf r})\cdot{{\bf A}^{a}}({\bf r})\Big{\}}$	(23)

may be admissible, provided interacting and non-interacting $N$ - and $v$ -representabilityhold true.

In practice, the form of GKS schemes depends upon the detail of $\hat{O}_{\rm{\tiny\tiny GKS}}$ and $E^{\rm{\tiny\tiny GKS}}_{\rm Hxc}$ .A prominent example is $\hat{O}_{\rm{\tiny\tiny GKS}}\equiv\hat{T}+\alpha\hat{W}$ where $\alpha\in(0,1]$ turns on the interaction in the GKS reference system — yet the minimization is restricted to single Slater determinants $\Phi$ , only —and $E^{{\rm{\tiny\tiny GKS}}}_{\rm Hxc}\equiv(1-\alpha)E^{\rm{\tiny\tiny DFA}}_{%\rm Hx}+E^{\rm{\tiny\tiny DFA}}_{\rm c}\equiv E^{\rm{\tiny\tiny DFA}}_{\rm Hxc}$ ; i.e.,we mix Fock exchange with (standard) LDAs or GGAs, in the form of typical global hybrid approximations.

Approximations of this kind can fix, at least partially, the self-interaction error of DFAs.They have been justified by the necessity of mimicking an exact (almost) semi-local xc-hole by the combination of a non-local exact-exchange hole with an approximate (semi-)local correlation hole. Global hybrids and refinements thereof, have been guided by the (so-called) adiabatic connection integration.Harris (1984); Becke (1993); Perdew etal. (1996a); Koch and Holthausen (2nd edition,2001); Garrick etal. (2022)But there is no prescription for choosing an optimal value of the hybridization parameter $\alpha$ that works for all systems.One needs to devise ways to find optimal values, driven by first principle calculations,Marques etal. (2011); Skone etal. (2014); Erba (2017) or consider more sophisticated forms and procedure of hybridizations.Kronik etal. (2012); Nguyen etal. (2018); Wing etal. (2021); Prokopiou etal. (2022)Furthermore, the optimization task always targets specific observables; most commonly, energy gaps.These aspects have received a huge attention and application in (S)DFT. In this work we show that switching from Spin-DFT to SCDFT is crucial for determining optimal mixing that can work both for band gaps and band splittings of spin-orbit coupled materials.

Hence, let us start with the functional form

	$\displaystyle E[\Phi]$	$\displaystyle=$	$\displaystyle T_{\rm{\tiny\tiny GKS}}[\Phi]+\alpha E^{\rm Fock}_{\rm x}[\Phi]+%E^{\rm{\tiny\tiny DFA}}_{\rm Hxc}\Big{[}n[\Phi],{{\vec{m}}}[\Phi],{\bf j}[\Phi%],{\vec{\bf J}}[\Phi]\Big{]}+\int d^{3}r~{}{n}[\Phi]\tilde{v}({\bf r})+\int d^%{3}r~{}v[\Phi]({\bf r})\tilde{B}^{a}({\bf r})$
		$\displaystyle+$	$\displaystyle\frac{1}{c}\int d^{3}r~{}{{\bf j}}[\Phi]({\bf r})\cdot{\bf A}({%\bf r})+\frac{1}{c}\int d^{3}r~{}{{\bf J}^{a}}[\Phi]({\bf r})\cdot{{\bf A}^{a}%}({\bf r})\;,$

where

\displaystyle E^{\rm Fock}_{\rm x}[\Phi]\equiv-\frac{1}{2}\int d^{3}r\int d^{3%}r^{\prime}~{}\frac{{\rm Tr}\left\{\Gamma({\bf r},{\bf r}^{\prime})\Gamma({\bfr%}^{\prime},{\bf r})\right\}}{|{{\bf r}-{\bf r}^{\prime}}|}\;,

(24)

is the Fock exchange, here evaluated with GKS spinors, and

\displaystyle\Gamma({\bf r},{\bf r}^{\prime})\equiv\sum_{k=1}^{N}\Phi_{k}({\bfr%})\Phi^{\dagger}_{k}({\bf r}^{\prime})\;

(25)

is the one-electron reduced density matrix (1RDM). In Eq. (24), Tr denotes the trace over spin.

As announced, we shall consider “typical” global hybrids forms. With “typical”, we intend forms that mix a fraction of Fock exchange, $E^{\rm Fock}_{x}[\Phi]$ , with GGAs (or lower rung approximations).

Note, GGAs in SCDFT (as in Spin-DFT) may dependent on all the basic variables and their gradients, but not on other quantities (e.g. kinetic energy densities).

The corresponding (generalized) KS equation, then, reads as follows

	$\displaystyle\hat{H}_{{\rm{\tiny\tiny GKS}}}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(-i\nabla+\frac{1}{c}{\mathbfcal{A}}_{\rm{\tiny%\tiny GKS}}\right)^{2}+\alpha\hat{{\cal V}}_{\rm x}^{\rm NL}$		(26)
		$\displaystyle+$	$\displaystyle{\cal V}_{{\rm{\tiny\tiny GKS}}}\;$		(26)

where

\displaystyle{\mathbfcal{A}_{\rm{\tiny\tiny GKS}}}=\mathbfcal{A}+(1-\alpha){%\mathbfcal{A}}^{{\rm{\tiny\tiny DFA}}}_{\rm x}+{\mathbfcal{A}}^{{\rm{\tiny%\tiny DFA}}}_{\rm c}

(27)

with

\displaystyle\mathbfcal{A}^{\rm{\tiny\tiny DFA}}_{\rm x/c}={\bf A}^{\rm{\tiny%\tiny DFA}}_{\rm x/c}+\sigma^{a}{\bf A}^{{\rm{\tiny\tiny DFA}},a}_{\rm x/c}\;,

(28)

\displaystyle\hat{{\cal V}}^{\rm NL}_{\rm x}\Phi_{k}

\displaystyle=\frac{\delta E^{\rm Fock}_{\rm x}[\Phi]}{\delta{\Phi_{k}}^{%\dagger}}=-\int d^{3}r^{\prime}~{}\frac{\Gamma({\bf r},{\bf r}^{\prime})\Phi_{%k}({\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}\;;

(29)

$\hat{{\cal V}}_{\rm x}^{\rm NL}$ is the Non-Local Fock potential — here evaluated with GKS spinors.Next,

	$\displaystyle{\cal{V}_{\rm{\tiny\tiny GKS}}}$	$\displaystyle=$	$\displaystyle{\cal{V}}+v_{\rm H}+(1-\alpha){\cal{V}}^{{\rm{\tiny\tiny DFA}}}_{%\rm x}+{\cal{V}}^{{\rm{\tiny\tiny DFA}}}_{\rm c}$		(30)
		$\displaystyle+$	$\displaystyle\frac{1}{2c^{2}}\left[\mathbfcal{A}^{2}-{\mathbfcal A}^{2}_{\rm{%\tiny\tiny GKS}}\right]$		(30)

with

\displaystyle v_{\rm H}

\displaystyle=\int d^{3}r^{\prime}~{}\frac{{\rm Tr}\;\Gamma({\bf r}^{\prime},{%\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}\;

(31)

and

\displaystyle{\cal V}^{\rm{\tiny\tiny DFA}}_{\rm x/c}=\left(v^{{\rm{\tiny\tinyDFA%}}}_{\rm x/c}+\sigma^{a}B^{{\rm{\tiny\tiny DFA}},a}_{\rm x/c}\right)\;.

(32)

The GKS equations reduce to the regular KS equations for $\alpha=0$ .For $\alpha\neq 0$ , the xc-scalar, xc-magnetic, and xc-vector potentials produced by $E^{\rm{\tiny\tiny DFA}}_{\rm Hxc}$ get, partially,replaced by a fraction of the non-local potential $\hat{{\cal V}}_{\rm x}^{\rm NL}$ . At $\alpha=1$ , the DFA gives no contribution: we end up with the HF approximation for the Pauli equation.In passing, also note that, had we invoked a meta-GGA instead of a GGA, the differentiation w.r.t. the spin-orbitals would have generated additional terms from the explicit dependence on the (spin-)kinetic-energy density — yielding to terms like the ones already accounted for in (non-collinear) SDFT.Peralta etal. (2007)

It is expedient to contrast the GKS equations including exact-exchange, Eqs. (26)-(32), against the exact-exchange approximation of the regular KS approach. In the present GKS scheme, the non-local Fock potential is directly given in terms of $\Gamma({\bf r},{\bf r}^{\prime})$ [see Eq. (29)]. On the other hand, in the regular KS approach to SCDFT, exact-exchange leads to the 16 integro-differential OEP equations that produce 16 local exchange potentialsin response to variations in 16 basic density components.Rohra and Görling (2006); Heaton-Burgessetal. (2007)At the present stage of development,the determination of local exact-exchange potentials are both numerically more involved and more costly than the evaluation of the non-local Fock potential.

III SCDFT versus SDFT+SOC

The GKS-SCDFT framework allows us to include SOC non-perturbatively and self-consistently in a density functional calculation. It is interesting to see how this general framework adapts to various approximations. For sake of simplicity, let us restrict to systems with vanishing magnetization ( ${\vec{m}}=0$ ) and vanishing particle currents ( ${\bf j}=0$ ) but with non-vanishing spin currents ( ${\vec{\bf J}}\neq 0$ , i.e. a typical time-reversal symmetry preserving system with SOC). Additionally, we consider systems for whichspin-currents are vanishing when SOC is turned off.Thus, there are two main types of approximations to be considered: approximations depending on all the basic densitiesbut currents; and those also including currents. Of the latter class, there is also the case of those approximations that includecurrents but only implicitly. In view of the complexity of this scenario, and in preparation of the calculations we shall perform in the next section,we first review analogous Gedanken calculations.

Standard GGAs. As a first example, we may perform a calculation by using a GGA that, like any standard GGA, does not include a dependence on spin currents. In detail,

\displaystyle E^{\rm GKS}_{\rm xc}[\Phi]\approx E^{\rm GGA}_{\rm xc}[n[\Phi]]\;.

(33)

Correspondingly, Eq. (15) implies ${\vec{\bf A}}^{\rm GGA}_{\rm xc}=0$ .Then, the GKS equations we must solve look like the GKS equations of SDFT calculations, if it were not that SOC is also added directly.Therefore, more appropriately, we shouldregard the resulting equations as some approximate GKS SCDFT equations. This interpretation is forced on us by the fact that, by construction, the SDFT energy functional is derived for systems described by the Hamiltonian reported in Eq. (1): i.e., for systems thatdo not include SOC. Via SDFT we may, however, include SOC as a perturbation.This state of affair must be borne in mind also whendiscussing the more involved cases reported below.An evaluation of SOC can be performed in a “second variational” step after the convergence of a “first variational” SDFT calculations.Huhn and Blum (2017); Koelling and Harmon (1977) Here, the “first variational” step is the diagonalization of the SDFT Hamiltonian (without SOC) in a set of basis functions, and the “second variational” step is the diagonalization of the SDFT Hamiltonian plus SOC in the basis of SDFT single particle orbitals (a one-shot calculation, that is ordinary first-order quasi-degenerate perturbation theory).It may be tempting to iterate the procedure self-consistently until convergence but, at the level of a regular GGA, the effects of self-consistency are, usually, irrelevant.Desmaraisetal. (2021b); Desmarais etal. (2023)Our formalism makes it apparent that the degree of the aforementioned self-consistency cannot make up for the missing dependence of an xc-approximation on the spin currents in the selected GGA.

Global Hybrid. As a more subtle and advanced example, let us perform a calculation by employing a regular hybrid functional; i.e. consisting in a fraction of Fock exchange plus a complementary fraction of a regular GGA. In detail,

\displaystyle E^{\rm GKS}_{\rm xc}[\Phi]\approx E^{\rm Hybrid}_{\rm xc}[\Phi]%\equiv\alpha E^{\rm Fock}_{\rm x}[\Phi]+(1-\alpha)E^{\rm GGA}_{\rm x}[n[\Phi]]%+E^{\rm GGA}_{\rm c}[n[\Phi]]\;.

(34)

In the absence of SOC, the Fock exchange can be restricted to the usual one component (i.e., globally collinear) form. When SOC is included, the Fock exchange must be upgraded to a two-component (i.e. non collinear) form — as we have described in the previous sections.Next, to appreciate the difference of more or less self-consistent calculations, it is instrumental to scrutinize the short-range behaviour of the Fock-energy density.Pittalis etal. (2017)

For this purpose, coming back to Eq. (24), let us employ the shorthand notation

Q_{\rm x}({\bf r},{\bf r}^{\prime})={\rm Tr}\left\{\Gamma({\bf r},{\bf r}^{%\prime})\Gamma({\bf r}^{\prime},{\bf r})\right\}\;,

(35)

and change integration variables by introducing the inter-particle separation $\mathbf{u}$ :

\displaystyle E^{\rm Fock}_{\rm x}[\Phi]

\displaystyle=

\displaystyle-\frac{1}{2}\int d^{3}r\int d^{3}u~{}\frac{Q_{\rm x}({\bf r}+{\bfu%}/2,{\bf r}-{\bf u}/2)}{u}\;.

We recall that $Q_{\rm x}$ from Eq. (35) is the trace of a $2\times 2$ matrix, and thus may be decomposed in the basis $I$ , $\sigma^{x},\sigma^{y},\sigma^{z}$ . Next, a Taylor expansion of the spherical average around $\mathbf{r}$ , $\langle Q_{\rm x}\rangle$ , to second-order in $u$ gives:

\displaystyle\langle Q_{\rm x}({\bf r},u)\rangle=\frac{n^{2}({\bf r})}{2}+%\frac{u^{2}}{6}\left[2n({\bf r})\tau({\bf r})-{\vec{\bf J}}({\bf r})\odot{\vec%{\bf J}}({\bf r})-\frac{n({\bf r})\;\nabla^{2}n({\bf r})}{4}\right]+\mathcal{O%}\left(u^{4}\right)\;.

(36)

Here, “ $\circ$ ” denotes contraction w.r.t. spin indices, and “ $\odot$ ” denotes a double contraction w.r.t. spin and real-space indices,

\tau({\bf r})=\frac{1}{2}\sum_{k=1}^{N}\Big{(}\nabla\Phi_{k}^{\dagger}({\bf r}%)\Big{)}\cdot\Big{(}\nabla\Phi_{k}({\bf r})\Big{)}\;,

(37)

is the kinetic energy density of the occupied GKS spinors.Eq. (36) shows that both the particle density and the spin currents contribute — via the spinors — to the Fock exchange energy-density at short range and, thus, they also contribute to the corresponding non-local potential $\hat{{\cal V}}_{\rm x}^{\rm NL}$ [see Eq. (29)]In passing also note that because the kinetic energy density does not belong to the basic variables of SCDFT, it must be regraded as a purely spinor-dependent term.

From Eq. (36), we can appreciate the important fact that a perturbative evaluation of SOC done in GKS-SDFT at the level of a hybrid approximation unavoidably misses the feedback from the implicit dependence on the spin currents. For going beyond such a perturbative evaluation, it is necessary to include SOC self-consistently. Adding SOC will turn on the spin-currents. Therefore,unconstrained two-component spin-orbital non-perturbative GKS-SDFT+SOC calculations are nothing else but GKS-SCDFT calculations.Because the dependence of $E^{\rm Fock}_{\rm x}[\Phi]$ on the spin currents is implicit, however,such an equivalence remains subtle.The crucial point is that, at variance with the previous GGA-only case,the self-consistency in this latter case gains on the feedback from the spin currents:both via the energy functional $E^{\rm Fock}_{\rm x}[\Phi]$ and in the non-local potential $\hat{{\cal V}}_{\rm x}^{\rm NL}$ .The effect of such a self-consistency may not be expected to be irrelevant.

Spin-current-dependent MGGAs.Let us lastly perform a calculation that invokes a functional form with an explicit dependence on the spin currents.This type of form can be generated, for example, by “localizing” the Fock exchange functional.Pittalis etal. (2017) A most direct way of doing this is to perform a Gaussian re-summation of Eq. (36).Straightforwardly, one gets

E^{\rm JLP}_{\rm x}[\Phi]=-\frac{3\pi}{4}\int d^{3}r~{}\frac{n^{3}({\bf r})}{%\left[\tau({\bf r})-\frac{{\vec{\bf J}}({\bf r})\odot{\vec{\bf J}}({\bf r})}{2%n({\bf r})}-\frac{\nabla^{2}n({\bf r})}{8}\right]}\;,

(38)

Eq. (38) is a spin-current dependent generalization of the Lee-Parr functional for exchange (JLP),Pittalis etal. (2017); Lee and Parr (1987) which, through Eq. (16), gives rise to the non-Abelian x-vector potential

{\vec{\bf A}}^{\rm JLP}_{\rm x}=-\frac{3\pi c}{16}\frac{n^{2}({\bf r})}{\left[%\tau({\bf r})-\frac{{\vec{\bf J}}({\bf r})\odot{\vec{\bf J}}({\bf r})}{2n({\bfr%})}-\frac{\nabla^{2}n({\bf r})}{8}\right]^{2}}{\vec{\bf J}}({\bf r})\;.

(39)

Therefore, Eq. (39) makes the switch from SDFT+SOC to SCDFT explicit also at the level of the GKS equations.The contribution of ${\vec{\bf A}}^{\rm JLP}_{\rm x}$ may not be expected to be minor.

In the next section, we shall make use of JLP to reduce the cost of the numerical calculationof hybrid calculation by using

E^{\rm GKS}_{\rm xc}[\Phi]\approx E^{\rm JMGGA}_{\rm xc}[\Phi]\equiv\alpha E^{%\rm JLP}_{\rm x}[\Phi]+(1-\alpha)E^{\rm GGA}_{\rm x}[n]+E^{\rm GGA}_{\rm c}[n]\;,

(40)

where $E^{\rm Fock}_{\rm x}[\Phi]$ in Eq. (34) has been replaced with $E^{\rm JLP}_{\rm x}[\Phi]$ .

In passing, we recall that the importance of the dependence of density functionals on the particle current has been largely demonstrated.Dobson (1993); Becke (1996, 2002); Maximoff etal. (2004); Burnus etal. (2005); Tao and Perdew (2005); Pittalis etal. (2007, 2009); Räsänen etal. (2009); Oliveira etal. (2010); Zhu etal. (2016); James W.Furness and Teale (2016); Holzer etal. (2022)The point stressed above concerns instead the spin current which, to the best of our knowledge, remains to be explored.

Last but not least, the derivation just above may be extended to generate MGGA-correlation forms by reducing the non-locality of higher-level xc models. This is left to future efforts together with other appealing possibilities as outlined in the road map of developments given below.

IV Applications

The aim of this section is to demonstrate the practical importance and flexibility of the GKS approach of SCDFT for going beyond the perturbative treatment of SOC.We shall consider the case of the SOC-induced/enhanced band splittings that occur near the top of the valence band of layered molibdenum dichalcogenides. These systems aretime-reversal symmetric and the ground states have vanishing magnetization ( ${\vec{m}}=0$ ), vanishing particle current ( ${\bf j}=0$ ),and, because of the presence of SOC, non-vanishing spin currents ( ${\vec{\bf J}}\neq 0$ ). The calculations we consider here showcase numerically the different SOC evaluations discussed formally in Sec. III.

Generalized Kohn-Sham Approach for the Electronic Band Structure of Spin-Orbit Coupled Materials (1)

Let us discuss formal aspects of valence band splittings in the presence of SOC and different symmetry constraints. The systems here considered preserve time-reversal symmetry (TRS):

\varepsilon_{k}^{\uparrow}({\bf k})=\varepsilon_{k}^{\downarrow}({\bf-k})\;,

(41)

where $\varepsilon_{k}$ are the energy values of band $k$ at different points of the first Brillouin zone (FBZ). Let us recall that space inversion symmetry (SIS) results in the following constraint on the band structure:

\varepsilon_{k}^{\sigma}({\bf k})=\varepsilon_{k}^{\sigma}({\bf-k})\;.

(42)

TRS and SIS together imply bands which are doubly degenerate in spin.The inclusion of SOC makes the Hamiltonian spin-dependent, correspondingly spin-up and spin-down states feel a different potential and split, if allowed by symmetry.

The single-layer MoSe ${}_{2}$ system preserves TRS but breaks SIS, which, in the presence of SOC, leads to possible spin-splittings of bands that would otherwise be doubly-degenerate at the SDFT level:

\varepsilon_{k}^{\uparrow}({\bf k})\neq\varepsilon_{k}^{\downarrow}({\bf k})\;.

(43)

In uniaxial (or low-dimensional) systems, such as 2D hexagonal MoSe ${}_{2}$ , the spin-splittings are embodied by the Rashba Hamiltonian (Rashba-I effect).Rashba (1960) Figure 1 A) shows such spin-splitting at the high-symmetry point K of the FBZ and along K- $\Gamma$ and K-M paths. At the SDFT level (black line), the top valence band is doubly degenerate. The spin degeneracy is lifted by SOC according to Eq. (43); for instance, see the SCDFT description (blue lines).

The $\alpha$ -MoTe ${}_{2}$ hexagonal crystal is characterized by stacked MoTe ${}_{2}$ layers along the c crystallographic axis, separated by van-der-Waals gaps. As both TRS and SIS are preserved, the combination of Eqs. (41) and (42) leads to:

\varepsilon_{k}^{\uparrow}({\bf k})=\varepsilon_{k}^{\downarrow}({\bf k})\;,

(44)

so that all bands are necessarily spin degenerate.Therefore, in the case of $\alpha$ -MoTe ${}_{2}$ , SOC enhanced band splittings are related to the dipole field of the locally asymmetric Mo crystallographic sites. This so-called Rashba-II effect results in the appearance of spatially localized “hidden spin valleys” associated with the band splittings.Zhang etal. (2014a); Oliva etal. (2020)Figure 1 B) shows such enhanced band splitting around K in the top of the valence band. At the SDFT level (black lines), the two top valence bands are spin doubly degenerate and are already split. With SOC, the bands are still doubly degenerate according to Eq. (44) but get further split (see the SCDFT blue lines).

Generalized Kohn-Sham Approach for the Electronic Band Structure of Spin-Orbit Coupled Materials (2)

We perform calculations with the Crystal23 package.Erba etal. (2022) We first report on calculations employing a global-hybrid functionalas in Eq. (34) , making sure to allow for unrestricted two-component single particle spinors, with the PBE generalized-gradient approximation.Perdew etal. (1996b) Computational details are reported in the supplementary material.ESI (see also Refs. Desmarais etal., 2019; Doll etal., 2001; Doll, 2001; Doll etal., 2006; Civalleri etal., 2001; Metz etal., 2000; Peterson etal., 2003, 2007; Laun etal., 2018; Heyd etal., 2005; Lebedev, 1976, 1977; Towler etal., 1996 therein).Below calculations are also performed and results reported for a regular GGA (PBE) andfor the spin current-dependent MGGA (JLP) as defined in Eq.(40).

We start by discussing Rashba-I type SOC-induced spin-splitting in 2D single-layer MoSe ${}_{2}$ . A graphical representation (side and top views) of the atomic structure of this system is given in Figure 2. The splitting occurs at the K point of the FBZ of the system and has been measured by angle-resolved photo-electron spectroscopy (ARPES) experiments (0.180-0.185 eV).Zhang etal. (2014b); Shim etal. (2014); Ross etal. (2013)

Figure 2 reports the computed spin-splittings as a function of $\alpha$ (i.e. fraction of Fock exchange).All geometries were fully optimized for each value of $\alpha$ .

The following is observed:(i) The black line reports the SDFT results: no spin-splitting is observed, as expected;(ii) The yellow lines describes the results from one-shot second-variational treatment of SOC (SDFT+SOC@SV). A value of 0.14 eV is obtained that significantly underestimates the experimental values, (almost) independently of $\alpha$ .;iii) The blue line shows the results from SCDFT calculations. The experimental band splittings are reproduced at values of $\alpha$ in the range 14-17%. Correspondingly, the effect on the band splitting is found to amount to 22% of the total SOC-induced splitting.

Generalized Kohn-Sham Approach for the Electronic Band Structure of Spin-Orbit Coupled Materials (3)

	GGA	MGGA	Hybrid	JMGGA
SDFT	0	0	0	0
SDFT+SOC@SV	0.14	0.14	0.14	0.14
SCDFT	0.14	0.14	0.18	0.18
Exp.			0.180 - 0.185

	GGA	MGGA	Hybrid	JMGGA
SDFT	0.22	0.25	0.23	0.25
SDFT+SOC@SV	0.22	0.31	0.29	0.31
SCDFT	0.22	0.31	0.32	0.33
Exp.			0.30 - 0.34

	MoSe ${}_{2}$	MoTe ${}_{2}$
GGA	1.47	0.71
Hybrid	1.99	1.05
JMGGA	1.65	0.71
Exp.	1.6-2.3	1.03

Next, we discuss Rashba-II type SOC-enhanced band-splitting in bulk $\alpha$ -MoTe ${}_{2}$ . A graphical representation of the atomic structure of this system is given in Figure 3.Here, two spin doubly-degenerate bands near the top of the valence are already split at the SDFT level (black lines) at the K point of the FBZ. The energy gap further widens upon inclusion of SOC, by an extent that depends on how spin-currents are treated.

The splitting has been measured by optical experiments (0.30-0.34 eV).Beal etal. (1972); Oliva etal. (2020); Ruppertetal. (2014a)Figure 3 compares the experimental values with computed band-splittings from different treatments of SOC as a function of $\alpha$ (i.e. fraction of Fock exchange). We note that, for this system, experimental values are more significantly spread, which results in a more difficult quantitative assessment of the different theoretical approaches. However, the following is observed: i) SDFT values visibly underestimate experimental results; ii) The SDFT+SOC@SVresults are better than the SDFT results, as expected;iii) The slope of the SDFT(+SOC@SV) results, however, is significantly different from the slope of the SCDFT results. As a consequence of which, agreement for the band splittings is obtained at an $\alpha$ which does not yield a band gap in agreement with experiments. Indeed, the SV calculation at a fraction $\alpha=0.2$ provides a splitting of 0.30 eV, but the band gap is much too large at 1.47 eV.;iv) The experimental band splitting is reproduced via SCDFT calculations at values of $\alpha$ in the range 0.06-0.15. It amounts to about 20% of the total band-splitting at a fraction $\alpha=0.10$ .

Tables 1 – 3 summarize our results on splittings and fundamental band gaps for MoSe ${}_{2}$ and MoTe ${}_{2}$ , employing, respectively, fractions $\alpha=0.15$ and $\alpha=0.10$ of Fock exchange. We reiterate that by employing one and the same hybrid form,we can reproduce the experimental band gaps and splittings with one and the same value of $\alpha$ (splittings of 0.18 and 0.32 eV, gaps of 1.99 and 1.05 eV, respectively, on MoSe ${}_{2}$ and MoTe ${}_{2}$ ).

We also stress that the difference in the slopes between the yellow and blue lines in Fig. 2 and Fig. 3 can unambiguously be attributed to the different treatment of the spin-currents.Such a dependence is implicitly encoded in the Fock exchange and it can be exploited by taking as an input spinors derived under the action of SOC.SDFT, however, neglects SOC from the outset; so spin currents vanish in the corresponding solutions.SDFT+SOC@SV accounts for SOC butevaluates Fock exchange at the level of SDFT spinors; thus, after the second variation step, spin-currents do not vanish but are not used as a feedback in the calculation.SCDFT, by construction, evaluates Fock exchange under the action of SOC; thus spin currents can drive the convergence toward more accurateself-consistent results.

Finally, we pass to the GGA and the JMGGA cases.Not surprisingly, of course, neither the experimental gaps nor the splittings can be reproduced with a regular GGA functional.Eq. (40) via Eqs. (38) and (39) makes explicit the dependence of the exchange energy on spin-currents and brings forth the corresponding non-Abelian exchange potential. JMGGA lowers the cost of the analogous global hybrid calculations — by reducing Eq. (34) to (40) — yet maintaining accurate band splittings (0.18 and 0.33 eV). In doing so, however, the fundamental band gaps are decreased (i.e. from 1.99 to 1.65 eV on MoSe ${}_{2}$ , from 1.05 to 0.71 eV on MoTe ${}_{2}$ ), which worsens the agreement against the experiment, when compared to the results of the full hybrid approximation. Nonetheless, the JMGGA gap of 1.65 eV on MoSe ${}_{2}$ is still an improvement over the pure GGA value of 1.47 eV. SDFT+SOC@SV calculations underestimate splittings (0.14 eV on MoSe ${}_{2}$ , 0.22-0.31 eV on MoTe ${}_{2}$ ).

In conclusion, among the cases here considered, only the non-local Fock potential allows for a simultaneous agreement against the experiment on both fundamental band gaps and band splittings. If only band splittings are required, comparably good results can be obtained by replacing the non-local Fock operator by a computationally cheaper and formally simpler semi-local spin-current dependent approximation.

IV.1 Near-future road-map for GKS-SCDFT

The results illustrated above show that GKS-SCDFT is readily useful. Two questions can be posed, however: (i) Will it be possible to get rid of the empiricism involved in the determination of $\alpha$ , the “optimal” fraction of exchange; or – we may ask — can the fraction be determined self-consistently in SCDFT without having to resort toother (computationally more demanding) methodologies? (ii) Will it be possible to derive more accurate functional approximations with an explicit dependence on the spin-current?We foresee that the answer to both questions is likely to be positive.

Question (i) may be resolved by upgrading a very recent development foroptimally-tuned range separated hybrids,Wing etal. (2021) which has been shown to work both for molecular and periodic materials.

Question (ii) may be answered by invoking the U(1) $\times$ SU(2)-gauge invariant extensionPittalis etal. (2017) of more evolved meta-GGAs than the case reported here for illustrative purposes.Work is in progress on the SCANSun etal. (2015) and TASKAschebrock and Kümmel (2019) energy functionals and the like.

V Outlooks and conclusions

We have put forward a generalization of the Kohn-Sham formalism (GKS), which admits the use of non-local effective potentials firmly rooted in SCDFT. This formulation is the analogous of the popular GKS formulation of (Spin-)DFT.Seidl etal. (1996)Here, we have spelled out and analyzed the novel and subtle aspects that are uniquely brought forth by the SCDFT framework.We have demonstrated via applications that GKS-SCDFT readily allows us to obtain results beyond the state-of-the-art in electronic structure calculations for spin-orbit coupled materials. By considering time-reversal symmetric spin-orbit coupled states, we have demonstrated that the dependence of the energy functional on spin currents is important even when it is only implicit, as in the prominent case of Fock exchange. Global hybrid approximations can yield significantly more accurate results when used in GKS-SCDFT calculations rather than in perturbative SDFT+SOC calculations.

In particular, we have applied GKS-SCDFT to the evaluationof band gaps and SOC-induced band splittings in materials of great interest in spintronics and valleytronics, with Rashba-I, Rashba-II effects.At the level of the global hybrid approximations, we have shown that by applying the self-consistent SCDFT treatment of spin-orbit interactions one can find an optimal fraction of Fock exchange which works well for both the fundamental band gaps and the SOC-induced or -enhanced band splittings.We have shown that the widely used method for refining Spin-DFT results via a second-variational treatment of SOC canfail to reproduce the experimental results – superior agreement can be achieved by switching from SDFT to full-fledged SCDFT calculations.

Efforts, in the near future, will be devoted to reduce empiricism in finding the “optimal” fraction of Fock exchange. We believe that the optimally-tuned range separated hybrids offer, presently, a valid and very promising option.Wing etal. (2021)Furthermore, the illustrative case reported here of a simple spin-current dependent meta-GGA suggests that it is appealing to develop theextension of more recent and more evolved meta-GGA forms to SCDFT as well.

Most importantly, already at this stage of the development, the GKS approach of Spin-Current DFT can offer significant improvements in the calculation of the electronic structure of challenging spin-orbit coupled materials.

Acknowledgements.

This research has received funding from the Project CH4.0 under the MUR program “Dipartimenti di Eccellenza 2023-2027” (CUP: D13C22003520001). GV was supported by the Ministry of Education, Singapore, under its Research Centre of Excellence award to the Institute for Functional Intelligent Materials (I-FIM, project No. EDUNC-33-18-279-V12). We are grateful to Stephen Dale for a reading of the manuscript.

References

Greiner (2012)W.Greiner,Relativistic Quantum Mechanics: Wave Equations(Springer Berlin Heidelberg, 2012), ISBN9783642880827.
Yu and Cardona (2010)P.Yu andM.Cardona,Fundamentals of Semiconductors: Physics and MaterialsProperties, Graduate Texts in Physics (Springer BerlinHeidelberg, 2010), ISBN 9783642007101.
Bernevig and Hughes (2013)B.Bernevig andT.Hughes,Topological Insulators and TopologicalSuperconductors (Princeton University Press,2013), ISBN 9780691151755.
Vanderbilt (2018)D.Vanderbilt,Berry Phases in Electronic Structure Theory: ElectricPolarization, Orbital Magnetization and Topological Insulators(Cambridge University Press, 2018), ISBN9781107157651.
Žutićetal. (2004)I.Žutić,J.Fabian, andS.DasSarma,Rev. Mod. Phys. 76,323 (2004).
Zhang and Vignale (2016)S.S.-L. Zhang andG.Vignale,Phys. Rev. B 94,140411 (2016).
Vignale and Rasolt (1988)G.Vignale andM.Rasolt,Phys. Rev. B 37,10685 (1988).
Bencheikh (2003)K.Bencheikh,J. Phys. A 36,11929 (2003).
Rohra and Görling (2006)S.Rohra andA.Görling,Phys. Rev. Lett. 97,013005 (2006).
Heaton-Burgessetal. (2007)T.Heaton-Burgess,P.Ayers, andW.Yang,Phys. Rev. Lett. 98,036403 (2007).
Abedinpour etal. (2010)S.H. Abedinpour,G.Vignale, andI.V. Tokatly,Phys. Rev. B 81,125123 (2010).
Kurth and Pittalis (2006)S.Kurth andS.Pittalis,The Optimized Effective Potential Method andLDA + U (John von Neumann Institute for Computing,Jülich, 2006), vol.31of NIC series, pp. 299–334.
Pittalis etal. (2006)S.Pittalis,S.Kurth,N.Helbig, andE.K.U. Gross,Phys. Rev. A 74,062511 (2006).
Sharma etal. (2007a)S.Sharma,J.K. Dewhurst,C.Ambrosch-Draxl,S.Kurth,N.Helbig,S.Pittalis,S.Shallcross,L.Nordström,and E.K.U.Gross, Phys. Rev. Lett.98, 196405(2007a).
Sharma etal. (2007b)S.Sharma,S.Pittalis,S.Kurth,S.Shallcross,J.K. Dewhurst,and E.K.U.Gross, Phys. Rev. B76, 100401(2007b).
Trushin and Görling (2018)E.Trushin andA.Görling,Phys. Rev. B 98,205137 (2018).
Ullrich (2018)C.A. Ullrich,Phys. Rev. B 98,035140 (2018).
Pluhar and Ullrich (2019)E.A. Pluhar andC.A. Ullrich,Phys. Rev. B 100,125135 (2019).
Pittalis etal. (2017)S.Pittalis,G.Vignale, andF.G. Eich,Phys. Rev. B 96,035141 (2017).
WU and YANG (2003)Q.WU andW.YANG, J.Theor. Comput. Chem. 02, 627(2003).
Gidopoulos and Lathiotakis (2012)N.I. Gidopoulosand N.N.Lathiotakis, Phys. Rev. A85, 052508(2012).
Eich and Hellgren (2014)F.G. Eich andM.Hellgren,J. Chem. Phys. 141(2014).
Gidopoulos and Lathiotakis (2013)N.I. Gidopoulosand N.N.Lathiotakis, Phys. Rev. A88, 046502(2013).
Oueis etal. (2023)Y.Oueis,G.N. Sizov, andV.N. Staroverov,J. Phys. Chem. A 127,2664 (2023).
Becke (1993)A.D. Becke,J. Chem. Phys. 98,1372 (1993).
Seidl etal. (1996)A.Seidl,A.Görling,P.Vogl,J.A. Majewski,and M.Levy,Phys. Rev. B 53,3764 (1996).
Perdew etal. (2017)J.P. Perdew,W.Yang,K.Burke,Z.Yang,E.K. Gross,M.Scheffler,G.E. Scuseria,T.M. Henderson,I.Y. Zhang,A.Ruzsinszky,etal., PNAS 114,2801 (2017).
Neumann etal. (1996)R.Neumann,R.H. Nobes, andN.C. Handy,Mol. Phys. 87,1 (1996).
Lehtola etal. (2020)S.Lehtola,F.Blockhuys,andC.VanAlsenoy,Mol. 25 (2020).
Wing etal. (2021)D.Wing,G.Ohad,J.B. Haber,M.R. Filip,S.E. Gant,J.B. Neaton,and L.Kronik,Proceedings of the National Academy of Sciences118, e2104556118(2021).
Desmaraisetal. (2020a)J.K. Desmarais,J.-P. Flament,and A.Erba,Phys. Rev. B 101,235142 (2020a).
Desmaraisetal. (2020b)J.K. Desmarais,J.-P. Flament,and A.Erba,Phys. Rev. B 102,235118 (2020b).
Desmaraisetal. (2021a)J.K. Desmarais,S.Komorovsky,J.-P. Flament,and A.Erba,J. Chem. Phys. 154,204110 (2021a).
Desmarais etal. (2019)J.K. Desmarais,J.-P. Flament,and A.Erba,J. Chem. Phys. 151,074107 (2019).
Comaskey etal. (2022)W.P. Comaskey,F.Bodo,A.Erba,J.L. Mendoza-Cortes,and J.K.Desmarais, Phys. Rev. B106, L201109(2022).
Bodo etal. (2022)F.Bodo,J.K. Desmarais,and A.Erba,Phys. Rev. B 105,125108 (2022).
Fröhlich and Studer (1993)J.Fröhlich andU.M. Studer,Rev. Mod. Phys. 65,733 (1993).
Levy (1982)M.Levy, Phys.Rev. A 26, 1200(1982).
Lieb (1983)E.H. Lieb,Int. J. Quantum Chem. 24,243 (1983).
(40)We note that the mathematical conditions for $N$ -representability and $v$ -representability in SCDFT are presently unknown.Non-trivial notable differences with the case of standard DFT have recentlybeen pointed out in SDFT, which make an ensemble formulation necessary evenfor $N$ -representabilityGontier (2013). Similar or bigger differences areexpected in SCDFT. Thus, the domain of the $F$ functional in Eq.(9) isnot known, and we have informally assumed that it is sufficiently dense toapproximate most density sets that may occur in physical systems. Likewise,when introducing the Kohn-Sham system, we will informally assumenon-interacting $v$ -representability based on single Slater determinants. Itis natural to expect that, the $N$ -representability and $v$ -representabilityissues which are prone to arise in SCDFT will be ameliorated (if notresolved) within an ensemble formulation, in which the admissiblenon-interacting pure states would also comprise linear combination ofdeterminants.
Garrick etal. (2022)R.Garrick,T.Gould, andL.Kronik,Adv. Theory Simul. 10(2022).
Harris (1984)J.Harris,Phys. Rev. A 29,1648 (1984).
Perdew etal. (1996a)J.P. Perdew,M.Ernzerhof,and K.Burke,J. Chem. Phys. 105,9982 (1996a).
Koch and Holthausen (2nd edition,2001)W.Koch andM.C. Holthausen,A Chemist’s Guide to Density FunctionalTheory (Wiley - VCH, Weinheim -New York, 2nd edition, 2001).
Marques etal. (2011)M.A.L. Marques,J.Vidal,M.J.T. Oliveira,L.Reining, andS.Botti,Phys. Rev. B 83,035119 (2011).
Skone etal. (2014)J.H. Skone,M.Govoni, andG.Galli,Phys. Rev. B 89,195112 (2014).
Erba (2017)A.Erba, J.Phys.: Condens. Matter 29,314001 (2017).
Kronik etal. (2012)L.Kronik,T.Stein,S.Refaely-Abramson,and R.Baer,J. Chem. Theory Comput. 8,1515 (2012).
Nguyen etal. (2018)N.L. Nguyen,N.Colonna,A.Ferretti, andN.Marzari,Phys. Rev. X 8,021051 (2018).
Prokopiou etal. (2022)G.Prokopiou,M.Hartstein,N.Govind, andL.Kronik,J. Chem. Theory Comput. 18,2331 (2022).
Peralta etal. (2007)J.E. Peralta,G.E. Scuseria,and M.J.Frisch, Phys. Rev. B75, 125119(2007),URL https://link.aps.org/doi/10.1103/PhysRevB.75.125119.
Huhn and Blum (2017)W.P. Huhn andV.Blum,Phys. Rev. Mater. 1,033803 (2017),URL https://link.aps.org/doi/10.1103/PhysRevMaterials.1.033803.
Koelling and Harmon (1977)D.Koelling andB.Harmon,J. Phys. C Sol. State Phys. 10,3107 (1977).
Desmaraisetal. (2021b)J.K. Desmarais,A.Erba,J.-P. Flament,and B.Kirtman,J. Chem. Theor. Comput. 17,4712 (2021b).
Desmarais etal. (2023)J.K. Desmarais,A.Boccuni,J.-P. Flament,B.Kirtman, andA.Erba, J.Chem. Theor. Comput. 19, 1853(2023).
Lee and Parr (1987)C.Lee andR.G. Parr,Phys. Rev. A 35,2377 (1987).
Dobson (1993)J.F. Dobson,The Journal of Chemical Physics98, 8870 (1993).
Becke (1996)A.D. Becke,Canadian Journal of Chemistry74, 995 (1996).
Becke (2002)A.D. Becke,The Journal of Chemical Physics117, 6935 (2002).
Maximoff etal. (2004)S.N. Maximoff,M.Ernzerhof,and G.E.Scuseria, The Journal of Chemical Physics120, 2105 (2004).
Burnus etal. (2005)T.Burnus,M.A. Marques,and E.K. Gross,Physical Review A 71,010501 (2005).
Tao and Perdew (2005)J.Tao andJ.P. Perdew,Phys. Rev. Lett. 95,196403 (2005).
Pittalis etal. (2007)S.Pittalis,S.Kurth,S.Sharma, andE.K.U. Gross,The Journal of Chemical Physics127, 124103(2007).
Pittalis etal. (2009)S.Pittalis,E.Räsänen,and E.K.U.Gross, Phys. Rev. A80, 032515(2009).
Räsänen etal. (2009)E.Räsänen,S.Pittalis,C.R. Proetto,and E.K.U.Gross, Phys. Rev. B79, 121305(2009).
Oliveira etal. (2010)M.J.T. Oliveira,E.Räsänen,S.Pittalis, andM.A.L. Marques,Journal of Chemical Theory and Computation6, 3664 (2010).
Zhu etal. (2016)W.Zhu,L.Zhang, andS.B. Trickey,The Journal of Chemical Physics145, 224106(2016),eprint http://aip.scitation.org/doi/pdf/10.1063/1.4971377,URL http://aip.scitation.org/doi/abs/10.1063/1.4971377.
James W.Furness and Teale (2016)T.H. James W.Furness,UlfEkström andA.M. Teale,Molecular Physics 114,1415 (2016).
Holzer etal. (2022)C.Holzer,Y.J. Franzke,and A.Pausch,The Journal of Chemical Physics157, 204102(2022), ISSN 0021-9606.
Camino etal. (2023)B.Camino,H.Zhou,E.Ascrizzi,A.Boccuni,F.Bodo,A.Cossard,D.Mitoli,A.M. Ferrari,A.Erba, andN.M. Harrison,Comput. Phys. Commun. 292,108853 (2023).
Rashba (1960)E.Rashba,Sov. Phys.-Solid State 2,1109 (1960).
Zhang etal. (2014a)X.Zhang,Q.Liu,J.-W. Luo,A.J. Freeman,and A.Zunger,Nat. Phys. 10,387 (2014a).
Oliva etal. (2020)R.Oliva,T.Woźniak,F.Dybala,J.Kopaczek,P.Scharoch, andR.Kudrawiec,Mat. Res. Lett. 8,75 (2020).
Zhang etal. (2014b)Y.Zhang,T.-R. Chang,B.Zhou,Y.-T. Cui,H.Yan,Z.Liu,F.Schmitt,J.Lee,R.Moore,Y.Chen, etal.,Nature nano. 9,111 (2014b).
Shim etal. (2014)G.W. Shim,K.Yoo,S.-B. Seo,J.Shin,D.Y. Jung,I.-S. Kang,C.W. Ahn,B.J. Cho, andS.-Y. Choi,ACS nano 8,6655 (2014).
Ross etal. (2013)J.S. Ross,S.Wu,H.Yu,N.J. Ghimire,A.M. Jones,G.Aivazian,J.Yan,D.G. Mandrus,D.Xiao,W.Yao, etal.,Nature Commun. 4,1474 (2013).
Erba etal. (2022)A.Erba,J.K. Desmarais,S.Casassa,B.Civalleri,L.Donà,I.J. Bush,B.Searle,L.Maschio,L.Edith-Daga,A.Cossard,etal., J. Chem. Theory Comput.(2022), URL 10.1021/acs.jctc.2c00958.
Perdew etal. (1996b)J.P. Perdew,K.Burke, andM.Ernzerhof,Physical review letters 77,3865 (1996b).
(79)Computational details are available at INSERT URL HERE.
Doll etal. (2001)K.Doll,V.Saunders, andN.Harrison,Int. J. Q. Chem. 82,1 (2001).
Doll (2001)K.Doll,Comput. Phys. Commun. 137,74 (2001).
Doll etal. (2006)K.Doll,R.Dovesi, andR.Orlando,Theor. Chem. Acc. 115,354 (2006).
Civalleri etal. (2001)B.Civalleri,P.D’Arco,R.Orlando,V.Saunders, andR.Dovesi,Chem. Phys. Lett. 348,131 (2001).
Metz etal. (2000)B.Metz,H.Stoll, andM.Dolg, J.Chem. Phys. 113, 2563(2000).
Peterson etal. (2003)K.A. Peterson,D.Figgen,E.Goll,H.Stoll, andM.Dolg, J.Chem. Phys. 119, 11113(2003).
Peterson etal. (2007)K.A. Peterson,D.Figgen,M.Dolg, andH.Stoll, J.Chem. Phys. 126 (2007).
Laun etal. (2018)J.Laun,D.VilelaOliveira,and T.Bredow,J. Comp. Chem. 39,1285 (2018).
Heyd etal. (2005)J.Heyd,J.E. Peralta,G.E. Scuseria,and R.L.Martin, J. Chem. Phys.123 (2005).
Lebedev (1976)V.I. Lebedev,USSR Comput. Math. Math. Phys.16, 10 (1976).
Lebedev (1977)V.I. Lebedev,Sib. Math. J. 18,99 (1977).
Towler etal. (1996)M.D. Towler,A.Zupan, andM.Causà,Comput. Phys. Commun. 98,181 (1996).
Beal etal. (1972)A.Beal,J.Knights, andW.Liang, J.Phys. C Sol. State Phys. 5, 3540(1972).
Ruppertetal. (2014a)C.Ruppert,B.Aslan, andT.F. Heinz,Nano Lett. 14,6231 (2014a).
Choi etal. (2017)B.K. Choi,M.Kim,K.-H. Jung,J.Kim,K.-S. Yu, andY.J. Chang,Nanoscale Res. Lett. 12,1 (2017).
Island etal. (2016)J.O. Island,A.Kuc,E.H. Diependaal,R.Bratschitsch,H.S. Van DerZant,T.Heine, andA.Castellanos-Gomez,Nanoscale 8,2589 (2016).
Zelewski andKudrawiec (2017)S.J. Zelewski andR.Kudrawiec,Sci. Reports 7,15365 (2017).
Ruppertetal. (2014b)C.Ruppert,B.Aslan, andT.F. Heinz,Nano Lett. 14,6231 (2014b).
Conan etal. (1979)A.Conan,D.Delaunay,A.Bonnet,A.Moustafa, andM.Spiesser,physica status solidi (b) 94,279 (1979).
Sun etal. (2015)J.Sun,A.Ruzsinszky,and J.P.Perdew, Phys. Rev. Lett.115, 036402(2015),URL https://link.aps.org/doi/10.1103/PhysRevLett.115.036402.
Aschebrock and Kümmel (2019)T.Aschebrock andS.Kümmel,Phys. Rev. Res. 1,033082 (2019),URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033082.
Gontier (2013)D.Gontier,Phys. Rev. Lett. 111,153001 (2013).

	$\displaystyle{\mathbfcal{A}}_{\rm{\tiny\tiny KS}}$	$\displaystyle=\left({\bf A}+\sigma^{a}{\bf A}^{a}\right)+\left({\bf A}_{{\rm xc%}}+{\bf A}^{a}_{{\rm xc}}\right)$
		$\displaystyle={\mathbfcal{A}}+{\mathbfcal{A}}_{\rm xc}\;,$		(13)

	$\displaystyle{\cal V}_{\rm{\tiny\tiny KS}}$	$\displaystyle=v_{\rm H}+\left(v+v_{\rm xc}\right)+\sigma^{a}\left(B^{a}+B^{a}_%{\rm xc}\right)$
		$\displaystyle+\frac{1}{2c^{2}}\left[\mathbfcal{A}^{2}-{\mathbfcal A}^{2}_{\rm{%\tiny\tiny KS}}\right]\;,$		(14)