# Perturbation Theory Treatment of Spin-Orbit Coupling. III: Coupled Perturbed
Method for Solids

Jacques K. Desmarais<EMAIL_ADDRESS>Dipartimento di
Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy  Alberto
Boccuni  Dipartimento di Chimica, Università di Torino, via Giuria 5, 10125
Torino, Italy  Jean-Pierre Flament  Université de Lille, CNRS, UMR 8523 —
PhLAM — Physique des Lasers, Atomes et Molécules, 59000 Lille, France  Bernard
Kirtman  Department of Chemistry and Biochemistry, University of California,
Santa Barbara, California 93106, USA  Alessandro Erba
<EMAIL_ADDRESS>Dipartimento di Chimica, Università di Torino, via
Giuria 5, 10125 Torino, Italy

###### Abstract

A previously proposed non-canonical coupled-perturbed Kohn-Sham density
functional theory (KS-DFT)/Hartree-Fock (HF) treatment for spin-orbit coupling
is here generalized to infinite periodic systems. The scalar-relativistic
periodic KS-DFT/HF solution, obtained with a relativistic effective core
potential, is taken as the zeroth-order approximation. Explicit expressions
are given for the total energy through 3rd-order, which satisfy the 2N + 1
rule (i.e. requiring only the 1st-order perturbed wave function for
determining the energy through 3rd-order). Expressions for additional 2nd-
order corrections to the perturbed wave function (as well as related one-
electron properties) are worked out at the uncoupled-perturbed level of
theory. The approach is implemented in the Crystal program and validated with
calculations of the total energy, electronic band structure, and density
variables of spin-current DFT on the tungsten dichalcogenide hexagonal bilayer
series (i.e. WSe2, WTe2, WPo2, WLv2), including 6p and 7p elements as a stress
test. The computed properties through second- or third-order match well with
those from reference two-component self-consistent field (2c-SCF)
calculations. For total energies, $E^{(3)}$ was found to consistently improve
the agreement against the 2c-SCF reference values. For electronic band
structures, visible differences w.r.t. 2c-SCF remained through second-order in
only the single-most difficult case of WLv2. As for density variables of spin-
current DFT, the perturbed electron density, being vanishing in first-order,
is the most challenging for the perturbation theory approach. The visible
differences in the electron densities are, however, largest close to the core
region of atoms and smaller in the valence region. Perturbed spin-current
densities, on the other hand, are well reproduced in all tested cases.

##  I Introduction

In modern electronic structure programs, relativistic effects are typically
accounted for through self-consistent field (SCF) treatments, either in a two-
or four-component spinor basis (2c-SCF or 4c-SCF),Chang et al. (1986); Buenker
et al. (1984); Almlöf et al. (1985); Douglas and Kroll (1974); Dyall (1997);
Iliaš et al. (2005); Kutzelnigg and Liu* (2006); Liu and Kutzelnigg (2007);
Liu and Peng (2006); Liu and Daoling (2009); Peng et al. (2007); Sikkema et
al. (2009); Seino and Hada (2008); Seino and Nakai (2012); Saue and Jensen
(1999); Jorgen Aa. Jensen et al. (1996) or by perturbation methods.Rutkowski
(1986a, b, c); Kutzelnigg (1996); van Wüllen (1995); Stopkowicz and Gauss
(2011, 2019); Cheng et al. (2013, 2018); Desmarais et al. (2021a, b) Within
these approaches, a particularly convenient representation of the Dirac
equation from a computational perspective is provided through the relativistic
effective core potential (RECP), including both scalar-relativistic and spin-
orbit coupling (SOC), or non-scalar, effects.Dolg and Cao (2012); Ermler et
al. (1981) In Part II of this series, we developed a molecular pure-state
coupled perturbed Kohn-Sham (CPKS) density functional theory (DFT) treatment
for including SOC effects, which was implemented within the Crystal
program.Desmarais et al. (2021b); Dovesi et al. (2018); Erba et al. (2022)
This approach has the potential for over an order of magnitude savings in
computation times w.r.t. 2c-SCF, and, moreover, provides a convenient starting
point for improvements to treat strongly-correlated, multi-reference, systems,
where an ensemble treatment would be necessary.Lieb (1983); Gross et al.
(1988a, b); Desmarais et al. (2021a, b); Desmarais and Kirtman (2022);
Helgaker and Teale (2022)

Our CPKS treatment is formulated within so-called spin-current DFT (SCDFT), in
which the exchange-correlation (xc) functional depends not only on the
particle-number density $\rho$ and z-component of the magnetization $m_{z}$,
(as in the usual spin-DFT of von Barth and Hedin),Von Barth and Hedin (1972)
but also on the other Cartesian components of the magnetization $m_{x}$ and
$m_{y}$, as well as the particle-current $\mathbf{j}$ and spin-current
$\mathbf{J}^{x}$, $\mathbf{J}^{y}$ and $\mathbf{J}^{z}$ densities.Vignale and
Rasolt (1988); Bencheikh (2003); Vignale and Rasolt (1987); Desmarais et al.
(2019a, b, 2021c, 2020a); Bodo et al. (2022); Desmarais et al. (2020b);
Comaskey et al. (2022)

In this paper, we generalize the molecular CPKS treatment to periodic systems.
After developing the formalism, our perturbation theory (PT) treatment is
implemented in a developmental version of the Crystal23 program,Erba et al.
(2022) and validated on a family of 2D tungsten dichalcogenide layered
compounds, where excellent agreement with reference 2c-SCF calculations is
demonstrated, not only for the total energy, but also for band structures and
SCDFT density variables.

##  II Formalism

###  II.1 Statement of the Problem

In the case of periodic systems, the spinors
$|\psi_{i,\mathbf{k}}\rangle=|\psi_{i,\mathbf{k}}^{\uparrow}\rangle\otimes|\uparrow\rangle+|\psi_{i,\mathbf{k}}^{\downarrow}\rangle\otimes|\downarrow\rangle$
are 2c crystalline orbitals (COs), with components
$|\psi_{i,\mathbf{k}}^{\sigma}\rangle$, expanded in a set of Bloch functions
(BFs) $\phi_{\mu,\mathbf{k}}^{\sigma}$:

$\psi_{i,\mathbf{k}}^{\sigma}\left(\mathbf{r}\right)=\sum_{\mu}^{N_{\mathcal{B}}}C^{\sigma}_{\mu,i}\left(\mathbf{k}\right)\phi_{\mu,\mathbf{k}}^{\sigma}\left(\mathbf{r}\right)$
(1)

where $\mathbf{k}$ is a point in the first Brillouin zone (FBZ),
$N_{\mathcal{B}}$ is the number of basis functions in a given cell of the
infinite-periodic system and $\sigma=\uparrow,\downarrow$ is a spin index.

In Crystal, the BFs are conveniently represented as a linear combination of
pure real atomic orbitals (LCAO), through the inverse-Fourier relation:

$\phi_{\mu,\mathbf{k}}^{\sigma}\left(\mathbf{r}\right)=\frac{1}{\sqrt{\Omega}}\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\chi_{\mu,\mathbf{g}}^{\sigma}\left(\mathbf{r}-\mathbf{A}_{\mu}\right)$ (2)

Here $\Omega$ is the volume of the FBZ, $\mathbf{g}$ is a direct-lattice
vector and $\mathbf{A}_{\mu}$ is the position in cell $\mathbf{g}$ at which
the AO $\chi_{\mu,\mathbf{g}}^{\sigma}$ is centered. In Eq. (2) we have
introduced the shorthand notation
$\chi_{\mu,\mathbf{g}}^{\sigma}\left(\mathbf{r}-\mathbf{A}_{\mu}\right)=\chi_{\mu}^{\sigma}\left(\mathbf{r}-\mathbf{A}_{\mu}-\mathbf{g}\right)$.
Variation of the orbitals $\psi_{i,\mathbf{k}}^{\sigma}$, under the constraint
of orthonormality:

$\langle\psi_{i,\mathbf{k}}^{\sigma}|\psi_{j,\mathbf{k}^{\prime}}^{\sigma^{\prime}}\rangle=\delta_{i,j}\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{\sigma,\sigma^{\prime}}\Rightarrow\mathbf{C}^{\dagger}\left(\mathbf{k}\right)\mathbf{S}\left(\mathbf{k}\right)\mathbf{C}\left(\mathbf{k}\right)=\mathbf{1}$
(3a) leads to the generalized Kohn-Sham (GKS) equation:
$\mathbf{H}\left(\mathbf{k}\right)\mathbf{C}\left(\mathbf{k}\right)=\mathbf{S}\left(\mathbf{k}\right)\mathbf{C}\left(\mathbf{k}\right)\mathbf{E}\left(\mathbf{k}\right)$
(3b)

where all matrices have size $2N_{\mathcal{B}}\times 2N_{\mathcal{B}}$,
$\mathbf{C}\left(\mathbf{k}\right)$ is the matrix of CO coefficients of Eq.
(1), $\mathbf{S}\left(\mathbf{k}\right)$ is the BF overlap matrix,
$\mathbf{E}\left(\mathbf{k}\right)$ is the matrix of Lagrange multipliers
(i.e. for canonical orbitals, corresponding to the diagonal matrix of band-
structure energy levels $\epsilon_{i,\mathbf{k}}$) and
$\mathbf{H}\left(\mathbf{k}\right)$ is the BF Hamiltonian matrix. Eq. (3b) can
be written more explicitly to highlight the structure in spin space:

$\displaystyle\begin{pmatrix}\mathbf{H}^{\uparrow\uparrow}\left(\mathbf{k}\right)&\mathbf{H}^{\uparrow\downarrow}\left(\mathbf{k}\right)\\\
\mathbf{H}^{\downarrow\uparrow}\left(\mathbf{k}\right)&\mathbf{H}^{\downarrow\downarrow}\left(\mathbf{k}\right)\end{pmatrix}\begin{pmatrix}\mathbf{C}^{\uparrow}\left(\mathbf{k}\right)\\\
\mathbf{C}^{\downarrow}\left(\mathbf{k}\right)\end{pmatrix}=\begin{pmatrix}\mathbf{S}^{\uparrow\uparrow}\left(\mathbf{k}\right)&\mathbf{0}\\\
\mathbf{0}&\mathbf{S}^{\downarrow\downarrow}\left(\mathbf{k}\right)\end{pmatrix}\begin{pmatrix}\mathbf{C}^{\uparrow}\left(\mathbf{k}\right)\\\
\mathbf{C}^{\downarrow}\left(\mathbf{k}\right)\end{pmatrix}\mathbf{E}\left(\mathbf{k}\right)$
(4)

In Eq. (4) and elsewhere, matrices with double and single spin indices have
size $N_{\mathcal{B}}\times N_{\mathcal{B}}$ and $N_{\mathcal{B}}\times
2N_{\mathcal{B}}$, respectively.
$\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)$, for instance, has
elements:

$H^{\sigma\sigma^{\prime}}_{\mu\nu}\left(\mathbf{k}\right)=\Omega\langle\phi_{\mu,\mathbf{k}}^{\sigma}|\hat{H}|\phi_{\nu,\mathbf{k}}^{\sigma^{\prime}}\rangle$
(5)

and:

$\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)=\mathbf{h}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)+\mathbf{J}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)-a\mathbf{K}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)+\mathbf{V}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)\;,$
(6)

in which $\mathbf{h}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)$ contains
the matrix elements that can be built from mono-electronic integrals:

$\mathbf{h}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)=\delta_{\sigma,\sigma^{\prime}}\left[\mathbf{v}\left(\mathbf{k}\right)+\mathbf{u}_{AR}\left(\mathbf{k}\right)\right]+\mathbf{u}_{SO}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)\;.$
(7)

Here, $\mathbf{v}$ consists of the electronic kinetic energy and electron-
nuclear interaction terms, $\mathbf{u}_{AR}$ and
$\mathbf{u}_{SO}^{\sigma\sigma^{\prime}}$ are, respectively, the averaged and
spin-orbit relativistic effective potential (AREP and SOREP) matrices; and
$\mathbf{J}^{\sigma\sigma^{\prime}}=\delta_{\sigma\sigma^{\prime}}\boldsymbol{\mathcal{J}}^{\sigma\sigma}$
and $\mathbf{K}^{\sigma\sigma^{\prime}}$ are the usual Coulomb and exact-
exchange terms (with $a$ being the included fraction of the latter).
$\mathbf{V}^{\sigma\sigma^{\prime}}$ is the matrix of DFT correlation and
exchange potentials (in either collinear or non-collinear
treatments).Desmarais et al. (2021b, c)

Inserting Eq. (2) into Eq. (5) (or the equivalent equation with $\hat{H}$
being replaced by any other operator), we are able to relate the BF matrix
$\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)$, for instance, to
the AO one $\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)$ through
the inverse-Fourier relation:

$\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)=\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)$ (8a) where AO
matrix elements of $\hat{H}$ or any other operator read:
$\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)=\langle\chi_{\mu,\mathbf{0}}^{\sigma}|\hat{H}|\chi_{\nu,\mathbf{g}}^{\sigma^{\prime}}\rangle$
(8b)

Following the arguments of Part I,Desmarais et al. (2021a) we recall (assuming
real AOs) that the diagonal spin-blocks of
$\underline{\mathbf{u}}_{SO}\left(\mathbf{g}\right)$ in Eq. (8b) are pure
imaginary:

$\mathcal{R}\left[\mathbf{u}^{\sigma\sigma}_{SO}\left(\mathbf{g}\right)\right]=\mathbf{0}\;,$
(9)

whereas the off-diagonal spin-blocks are complex.

In order to develop a computationally convenient PT approach that satisfies
Wigner’s 2N+1 rule, it turns out to be important to introduce a relation
similar to Eq. (9) in the BF basis. This can be achieved by returning to Eq.
(8a) and carrying out the inverse-Fourier transform individually on the real
and imaginary parts of
$\mathbf{u}^{\sigma\sigma^{\prime}}_{SO}\left(\mathbf{g}\right)$ and all other
matrices entering the Hamiltonian, for instance:

$\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{R}}\left(\mathbf{k}\right)=\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\Re\left[\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)\right]$
(10a) and:
$\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{I}}\left(\mathbf{k}\right)=\imath\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\Im\left[\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)\right]$
(10b) so that:
$\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)=\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{R}}\left(\mathbf{k}\right)+\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{I}}\left(\mathbf{k}\right)$
(10c)

We note at this point that both
$\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{R}}\left(\mathbf{k}\right)$ and
$\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{I}}\left(\mathbf{k}\right)$ are
complex quantities. In the case of the SOC operator, using Eq. (9) and
proceeding as in Eq. (10a), we obtain:

$\mathbf{u}^{\sigma\sigma\mathcal{R}}_{SO}\left(\mathbf{k}\right)=\mathbf{0}\;.$
(11)

###  II.2 Order by Order Expressions for the Orthonormality Conditions

For the expansion of Eq. (3a) in orders of PT, it proves useful to write the
KS coefficients in terms of $N_{\mathcal{B}}\times N_{\mathcal{B}}$ blocks
with double spin indices:

$\begin{pmatrix}\mathbf{C}^{\uparrow\uparrow}\left(\mathbf{k}\right)&\mathbf{C}^{\downarrow\uparrow}\left(\mathbf{k}\right)\\\
\mathbf{C}^{\uparrow\downarrow}\left(\mathbf{k}\right)&\mathbf{C}^{\downarrow\downarrow}\left(\mathbf{k}\right)\end{pmatrix}^{\dagger}\begin{pmatrix}\mathbf{S}^{\uparrow\uparrow}\left(\mathbf{k}\right)&\mathbf{0}\\\
\mathbf{0}&\mathbf{S}^{\downarrow\downarrow}\left(\mathbf{k}\right)\end{pmatrix}\begin{pmatrix}\mathbf{C}^{\uparrow\uparrow}\left(\mathbf{k}\right)&\mathbf{C}^{\downarrow\uparrow}\left(\mathbf{k}\right)\\\
\mathbf{C}^{\uparrow\downarrow}\left(\mathbf{k}\right)&\mathbf{C}^{\downarrow\downarrow}\left(\mathbf{k}\right)\end{pmatrix}=\mathbf{1}$
(12)

As in Part II, we take the scalar-relativistic (SR) unrestricted GKS solution
as the zeroth-order problem, which yields:

$\left[\mathbf{C}^{\sigma\sigma(0)}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{S}^{\sigma\sigma}\left(\mathbf{k}\right)\mathbf{C}^{\sigma\sigma(0)}\left(\mathbf{k}\right)=\mathbf{1}$
(13)

Then, introducing the $N$th-order orbital rotation matrices $\mathbf{U}^{(N)}$
defined through:

$\mathbf{C}^{\sigma^{\prime}\sigma(N)}\left(\mathbf{k}\right)=\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}\left(\mathbf{k}\right)\mathbf{U}^{\sigma^{\prime}\sigma(N)}\left(\mathbf{k}\right)$
(14)

leads to:

$\mathbf{U}^{\sigma^{\prime}\sigma(1)}\left(\mathbf{k}\right)=-\left[\mathbf{U}^{\sigma^{\prime}\sigma(1)}\left(\mathbf{k}\right)\right]^{\dagger}$
(15)

in first order and:

$\mathbf{U}^{\sigma^{\prime}\sigma(2)}\left(\mathbf{k}\right)+\left[\mathbf{U}^{\sigma^{\prime}\sigma(2)}\left(\mathbf{k}\right)\right]^{\dagger}=-\sum_{\sigma^{\prime\prime}}\left[\mathbf{U}^{\sigma^{\prime}\sigma^{\prime\prime}(1)}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{U}^{\sigma^{\prime\prime}\sigma(1)}\left(\mathbf{k}\right)$
(16)

in second order.

###  II.3 Order by Order Expressions for the Hamiltonian, as well as the
Perturbation Equations and their Solution

For Eq. (3b) we take the SR unrestricted GKS Hamiltonian as the zeroth-order
approximation:

$\displaystyle\mathbf{H}^{\sigma\sigma^{\prime}(0)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\delta_{\sigma\sigma^{\prime}}\Big{[}\mathbf{v}\left(\mathbf{k}\right)+\mathbf{u}_{AR}\left(\mathbf{k}\right)+\mathbf{J}^{\sigma\sigma(0)}\left(\mathbf{k}\right)$
(17) $\displaystyle-$ $\displaystyle
a\mathbf{K}^{\sigma\sigma(0)}\left(\mathbf{k}\right)+\mathbf{V}^{\sigma\sigma(0)}\left(\mathbf{k}\right)\Big{]}$

In first order:

$\mathbf{H}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)=\mathbf{u}_{SO}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)+\mathbf{J}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)-a\mathbf{K}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)+\mathbf{V}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)$
(18)

and for all orders $N$ greater than one:

$\mathbf{H}^{\sigma\sigma^{\prime}(N)}\left(\mathbf{k}\right)=\mathbf{J}^{\sigma\sigma^{\prime}(N)}\left(\mathbf{k}\right)-a\mathbf{K}^{\sigma\sigma^{\prime}(N)}\left(\mathbf{k}\right)+\mathbf{V}^{\sigma\sigma^{\prime}(N)}\left(\mathbf{k}\right)$
(19)

Then, following Part II, expansion of Eq. (3b) in orders of PT leads, in the
first order, to:

$\mathbf{G}^{(1)}\left(\mathbf{k}\right)+\mathbf{E}^{(0)}\left(\mathbf{k}\right)\
\mathbf{U}^{(1)}\left(\mathbf{k}\right)=\mathbf{U}^{(1)}\left(\mathbf{k}\right)\
\mathbf{E}^{(0)}\left(\mathbf{k}\right)+\mathbf{E}^{(1)}\left(\mathbf{k}\right)\;,$
(20)

and, in second order, to:

$\displaystyle\mathbf{G}^{(2)}\left(\mathbf{k}\right)+\mathbf{E}^{(0)}\left(\mathbf{k}\right)\
\mathbf{U}^{(2)}\left(\mathbf{k}\right)+\mathbf{G}^{(1)}\left(\mathbf{k}\right)\
\mathbf{U}^{(1)}\left(\mathbf{k}\right)=\mathbf{U}^{(2)}\left(\mathbf{k}\right)\
\mathbf{E}^{(0)}\left(\mathbf{k}\right)+\mathbf{U}^{(1)}\left(\mathbf{k}\right)\
\mathbf{E}^{(1)}\left(\mathbf{k}\right)+\mathbf{E}^{(2)}\left(\mathbf{k}\right)\;,$
(21)

where we have made use of the $N$th order Hamiltonian matrix in the CO basis:

$\mathbf{G}^{(N)}\left(\mathbf{k}\right)=\left[\mathbf{C}^{(0)}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{H}^{(N)}\left(\mathbf{k}\right)\mathbf{C}^{(0)}\left(\mathbf{k}\right)$
(22)

For further development, and following Eq. (10c), it is useful to define:

$\mathbf{G}^{(N)}\left(\mathbf{k}\right)=\mathbf{G}^{\mathcal{R}(N)}\left(\mathbf{k}\right)+\mathbf{G}^{\mathcal{I}(N)}\left(\mathbf{k}\right)$
(23a) as well as:
$\mathbf{U}^{(N)}\left(\mathbf{k}\right)=\mathbf{U}^{\mathcal{R}(N)}\left(\mathbf{k}\right)+\mathbf{U}^{\mathcal{I}(N)}\left(\mathbf{k}\right)$
(23b) where:
$\mathbf{G}^{\mathcal{C}(N)}\left(\mathbf{k}\right)=\left[\mathbf{C}^{(0)}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{H}^{\mathcal{C}(N)}\left(\mathbf{k}\right)\mathbf{C}^{(0)}\left(\mathbf{k}\right)$
(23c) and $\mathcal{C}=\mathcal{R}\text{ or }\mathcal{I}$.

Recalling from Eq. (11) that
$\mathbf{u}^{\sigma\sigma\mathcal{R}}_{SO}\left(\mathbf{k}\right)$ is zero,
the Coulomb and exchange response to a vanishing perturbation is also
vanishing. Thus, inserting Eqs. (11) and (18) in Eq. (23c), we obtain, for
diagonal spin blocks:

$\mathbf{G}^{\sigma\sigma\mathcal{R}(1)}\left(\mathbf{k}\right)=\mathbf{0}$
(24)

As will be discussed below in Section II.4, partitioning of matrices in
$\mathcal{R}$ and $\mathcal{I}$ contributions, as in Eq. (23a) and (23b)
allows us to relate elements upon inversion in reciprocal space.

As far as the GKS perturbed wave functions and band structures are concerned,
again following Part II, we solve Eqs. (15), (16), (20) and (21). In doing so,
we take advantage of the fact that the off-diagonal occupied-virtual blocks of
the Lagrange-multiplier matrix vanish at all orders $N$:

$\mathbf{E}_{VO}^{\sigma\sigma^{\prime}(N)}\left(\mathbf{k}\right)=\mathbf{E}_{OV}^{\sigma^{\prime}\sigma(N)}\left(\mathbf{k}\right)=\mathbf{0}$
(25)

In first order, the solution of Eqs. (15) and (20), together with Eq. (24)
provides, for the virtual-occupied (VO) block:

$\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{VO}\left(\mathbf{k}\right)=\frac{\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{VO}\left(\mathbf{k}\right)}{\epsilon_{O}^{\sigma^{\prime}(0)}\left(\mathbf{k}\right)-\epsilon_{V}^{\sigma(0)}\left(\mathbf{k}\right)}=-\left[\mathbf{U}^{\sigma^{\prime}\sigma\mathcal{C}(1)}_{OV}\left(\mathbf{k}\right)\right]^{\ast}$
(26)

while for the occupied-occupied (OO′) and virtual-virtual (VV′) blocks:

$\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{OO^{\prime}}\left(\mathbf{k}\right)=\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{VV^{\prime}}\left(\mathbf{k}\right)=\mathbf{0}$
(27)

in the non-canonical treatment. It follows that, for the occupied-occupied
block:

$\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)=\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{R}(1)}_{OO^{\prime}}\left(\mathbf{k}\right)+\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{I}(1)}_{OO^{\prime}}\left(\mathbf{k}\right)$
(28)

with an exactly analogous expression also for the virtual-virtual block
$\mathbf{E}_{VV^{\prime}}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)$.
A single diagonalization of, for instance,
$\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(0)}\left(\mathbf{k}\right)+\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)$
at the end of the calculation (i.e. after a coupled solution of Eq. 20) yields
the occupied band structure correct through first-order in PT.

At second order, we proceed similarly and write $\mathbf{U}^{(2)}$ as in Eq.
(23b) in terms of contributions $\mathbf{U}^{\mathcal{R}(2)}$, and
$\mathbf{U}^{\mathcal{I}(2)}$ that end up having different behaviour upon
inversion in reciprocal space - cf. Eq. (35). Then, the simultaneous non-
canonical solution of Eqs. (21) and (16) yields for the occupied-occupied
blocks:

$\displaystyle\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(2)}_{OO^{\prime}}\left(\mathbf{k}\right)=\left[\mathbf{U}^{\sigma^{\prime}\sigma\mathcal{R}(2)}_{O^{\prime}O}\left(\mathbf{k}\right)\right]^{\ast}=$
$\displaystyle-\frac{1}{2}\sum_{\mathcal{C}}\sum_{\sigma^{\prime\prime}}\sum_{V}\left[\mathbf{U}^{\sigma^{\prime\prime}\sigma\mathcal{C}(1)}_{VO}\left(\mathbf{k}\right)\right]^{\ast}\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}(1)}_{VO^{\prime}}\left(\mathbf{k}\right)$
(29a) and:
$\displaystyle\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(2)}_{OO^{\prime}}\left(\mathbf{k}\right)=\left[\mathbf{U}^{\sigma^{\prime}\sigma\mathcal{I}(2)}_{O^{\prime}O}\left(\mathbf{k}\right)\right]^{\ast}=$
$\displaystyle-\frac{1}{2}\sum_{\mathcal{C}}\sum_{\mathcal{C}^{\prime}\neq\mathcal{C}}\sum_{\sigma^{\prime\prime}}\sum_{V}\left[\mathbf{U}^{\sigma^{\prime\prime}\sigma\mathcal{C}(1)}_{VO}\left(\mathbf{k}\right)\right]^{\ast}\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}^{\prime}(1)}_{VO^{\prime}}\left(\mathbf{k}\right)$

Exactly analogous expressions can also be determined for
$\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(2)}_{VV^{\prime}}\left(\mathbf{k}\right)$
and
$\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(2)}_{VV^{\prime}}\left(\mathbf{k}\right)$.
To obtain Eqs. (29) and (29), we insert Eqs. (26) and (27) into Eq. (16). The
expressions for the VO blocks, on the other hand are provided in Appendix A.
In Eqs. (29), (29) and elsewhere the products inside sums over $O$ or $V$ are
performed element-wise.

Eqs. (29), (29), (45a) and (45b) are consistent with a set of non-canonical
second-order Lagrange multipliers that are provided in Appendix B. A single
diagonalization of, for instance,
$\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(0)}\left(\mathbf{k}\right)+\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(1)}\left(\mathbf{k}\right)+\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(2)}\left(\mathbf{k}\right)$
at the end of the calculation (i.e. after a coupled solution of Eqs. 20 and
21) provides occupied band structure energy levels that are correct through
second-order in PT.

As will be shown in Section IV, it turns out that an excellent approximation
to the full solution of Eq. (3b) is given by the coupled solution of Eq. (20),
followed by an uncoupled solution of Eq. (21). Such a procedure, which we call
the PT2′ approximation is mathematically achieved by simply setting:

$\text{PT2${}^{\prime}$
::}\quad\mathbf{G}^{\mathcal{R}(2)}\left(\mathbf{k}\right),\mathbf{G}^{\mathcal{I}(2)}\left(\mathbf{k}\right)\Rightarrow\mathbf{0}$
(30)

in Eqs. (45a), (45b), as well as Eqs. (46) and (46).

###  II.4 Some $\mathbf{k}$ to $-\mathbf{k}$ Relations

Given that the zeroth-order direct-space SR Hamiltonian, being the inverse-
Fourier transform of Eq. (17), is pure-real, we obtain:

$\mathbf{H}^{\sigma\sigma^{\prime}(0)}\left(\mathbf{k}\right)=\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\mathbf{H}^{\sigma\sigma^{\prime}(0)}\left(\mathbf{g}\right)=\left[\mathbf{H}^{\sigma\sigma^{\prime}(0)}\left(-\mathbf{k}\right)\right]^{\ast}$
(31)

This leads to the following well-known result for the CO coefficient and
Lagrange multiplier matrices at $\mathbf{k}$ and $-\mathbf{k}$:Desmarais et
al. (2020b)

$\mathbf{C}^{(0)}\left(\mathbf{k}\right)=\left[\mathbf{C}^{(0)}\left(-\mathbf{k}\right)\right]^{\ast}$
(32a) and:
$\mathbf{E}^{(0)}\left(\mathbf{k}\right)=\mathbf{E}^{(0)}\left(-\mathbf{k}\right)$
(32b)

Proceeding as in Eq. (31) at order $N>0$, we obtain:

$\displaystyle\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{R}(N)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\Re\left[\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)\right]=\left[\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{R}(N)}\left(-\mathbf{k}\right)\right]^{\ast}$
(33a)
$\displaystyle\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{I}(N)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\imath\sum_{\mathbf{g}}e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\Im\left[\mathbf{H}^{\sigma\sigma^{\prime}}\left(\mathbf{g}\right)\right]=-\left[\mathbf{H}^{\sigma\sigma^{\prime}\mathcal{I}(N)}\left(-\mathbf{k}\right)\right]^{\ast}$
(33b)

Then, inserting Eqs. (33) and (32a) into Eq. (23c) gives:

$\displaystyle\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{R}(N)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\left[\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{R}(N)}\left(-\mathbf{k}\right)\right]^{\ast}$
(34a)
$\displaystyle\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{I}(N)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle-\left[\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{I}(N)}\left(-\mathbf{k}\right)\right]^{\ast}$
(34b)

and then, inserting Eqs. (32b) and (34) into Eqs. (26)-(27) as well as Eqs.
(29), (29), (45a) and (45b) gives:

$\displaystyle\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(N)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(N)}\left(-\mathbf{k}\right)\right]^{\ast}$
(35a)
$\displaystyle\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(N)}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle-\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(N)}\left(-\mathbf{k}\right)\right]^{\ast}$
(35b)

###  II.5 Order by Order Expressions for the Direct-Space Density Matrix

The GKS density matrix, obtained as a solution of Eq. (3b) reads:

$\mathbf{P}^{\sigma\sigma^{\prime}}\left(\mathbf{k}\right)=\sum_{\sigma^{\prime\prime}}\mathbf{C}^{\sigma\sigma^{\prime\prime}}\left(\mathbf{k}\right)\mathbf{f}_{\sigma^{\prime\prime}}\left(\mathbf{k}\right)\left[\mathbf{C}^{\sigma^{\prime\prime}\sigma^{\prime}}\left(\mathbf{k}\right)\right]^{\dagger}$
(36)

where $\mathbf{f}_{\sigma}\left(\mathbf{k}\right)$ is the diagonal matrix of
band occupations (i.e. for insulators the matrix elements are $1$ for occupied
bands and $0$ for virtual bands, or, in general fractional values for
partially occupied bands in metals). Expanding the CO coefficients in orders
of perturbation theory, using Eqs. (14), we obtain, for instance, in zeroth
order:

$\mathbf{P}^{\sigma\sigma^{\prime}(0)}\left(\mathbf{k}\right)=\delta_{\sigma\sigma^{\prime}}\mathbf{C}^{\sigma\sigma(0)}\left(\mathbf{k}\right)\mathbf{f}_{\sigma}(\mathbf{k})\left[\mathbf{C}^{\sigma\sigma(0)}\left(\mathbf{k}\right)\right]^{\dagger}$
(37a) in first order:
$\displaystyle\mathbf{P}^{\sigma\sigma^{\prime}(1)}(\mathbf{k})$
$\displaystyle=$
$\displaystyle\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{f}_{\sigma}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}(1)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}$
(37b) $\displaystyle+$
$\displaystyle\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}(1)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}}(\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}$
and, in second order, an expression that is provided in Appendix C.

The matrices in direct space are obtained from an inverse-Fourier transform:

$\displaystyle\mathbf{P}^{\sigma\sigma^{\prime}(N)}(\mathbf{g})=\frac{1}{\Omega}\int
d\mathbf{k}\ e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\mathbf{P}^{\sigma\sigma^{\prime}(N)}(\mathbf{k})=\frac{1}{\Omega}\int^{\prime}d\mathbf{k}\Big{\\{}e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{P}^{\sigma\sigma^{\prime}(N)}(\mathbf{k})+e^{-\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{P}^{\sigma\sigma^{\prime}(N)}(-\mathbf{k})\Big{\\}}$
(38)

where the prime over the integral indicates that integration is limited to the
portion of the FBZ with positive coordinates. Inserting Eqs. (26), (27) and
(38) into Eq. (37b) and making use of Eq. (32a) and (35) to relate elements at
$\mathbf{k}$ and $-\mathbf{k}$, gives, in first order, for the case of
diagonal spin-blocks of $\mathbf{P}^{(1)}$ (see Appendix C for details):

$\Re\left[\mathbf{P}^{\sigma\sigma(1)}(\mathbf{g})\right]=\mathbf{0}\\\ $ (39)

and

$\displaystyle\Im\left[\mathbf{P}^{\sigma\sigma(1)}(\mathbf{g})\right]=\frac{2}{\Omega}\Im\int^{\prime}\\!\\!d\mathbf{k}\Big{\\{}e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}_{O}(\mathbf{k})\mathbf{f}_{\sigma
O}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma\mathcal{I}(1)}_{OV}(\mathbf{k})\right]^{\dagger}\\!\\!\left[\mathbf{C}^{\sigma\sigma(0)}_{V}(\mathbf{k})\right]^{\dagger}\\!\\!+e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}_{V}(\mathbf{k})\mathbf{U}^{\sigma\sigma\mathcal{I}(1)}_{VO}(\mathbf{k})\mathbf{f}_{\sigma
O}(\mathbf{k})\left[\mathbf{C}^{\sigma\sigma(0)}_{O}(\mathbf{k})\right]^{\dagger}\\!\\!\Big{\\}}$
(40)

The expressions for the full $\mathbf{P}^{(1)}(\mathbf{g})$ and
$\mathbf{P}^{(2)}(\mathbf{g})$ are provided in Appendix C. To obtain Eqs. (39)
and (40), we insert Eqs. (24) and (26) into Eq. (50). As was the case in Part
II of this series,Desmarais et al. (2021b) the vanishing
$\Re\left[\mathbf{P}^{\sigma\sigma(1)}(\mathbf{g})\right]$ leads to a
vanishing first-order electron-density response, and thus a vanishing Coulomb
and exchange-correlation response (for functionals not depending on the
particle- and spin-current densities). This leads to an uncoupled-perturbed
procedure for pure LDA and GGA functionals, and a coupled-perturbed procedure
for hybrid functionals, including a fraction of exact Fock exchange.

###  II.6 Order by Order Expressions for the Total Energy

Having set up all pieces of the PT, we are ready to derive computationally
convenient expressions for the total energy that are consistent with Wigner’s
2N+1 rule. The full derivation, which also resembles that given in Part II of
this series,Desmarais et al. (2021b) is provided in the electronic supporting
information (ESI).ESI  Here we report the final expressions. In first order,
the energy corrections due to SOC are vanishing:

$E^{(1)}=0$ (41)

For second and third order, it proves convenient to introduce the matrix of
SOC integrals in the CO basis:

$\boldsymbol{\Xi}^{\sigma\sigma^{\prime}\mathcal{C}}(\mathbf{k})=\left[\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\right]^{\dagger}\mathbf{u}^{\sigma\sigma^{\prime}\mathcal{C}}_{SO}\left(\mathbf{k}\right)\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})$
(42)

in terms of which we obtain:

$E^{(2)}=\frac{2}{\Omega}\int^{\prime}d\mathbf{k}\sum_{\mathcal{C}}\sum_{\sigma\sigma^{\prime}}\Re\text{Tr}\\{\mathbf{f}_{\sigma
O}(\mathbf{k})\boldsymbol{\Xi}^{\sigma\sigma^{\prime}\mathcal{C}}_{OV}(\mathbf{k})\mathbf{U}^{\sigma^{\prime}\sigma\mathcal{C}(1)}_{VO}(\mathbf{k})\\}$
(43)

and

$\displaystyle E^{(3)}$ $\displaystyle=$
$\displaystyle\frac{2}{\Omega}\int^{\prime}d\mathbf{k}\sum_{\mathcal{C},\mathcal{C}^{\prime}}\sum\limits_{\mathcal{C}^{\prime\prime}}{\phantom{\sum}}^{\prime}\sum_{\sigma\sigma^{\prime}\sigma^{\prime\prime}}\Re\text{Tr}\Big{\\{}\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{VV^{\prime}}(\mathbf{k})\mathbf{U}^{\sigma^{\prime}\sigma^{\prime\prime}\mathcal{C}^{\prime}(1)}_{V^{\prime}O}(\mathbf{k})\mathbf{f}_{\sigma^{\prime\prime}O}(\mathbf{k})\left[\mathbf{U}^{\sigma^{\prime\prime}\sigma\mathcal{C}^{\prime\prime}(1)}_{OV}(\mathbf{k})\right]^{\dagger}$
(44) $\displaystyle-$ $\displaystyle\frac{1}{2}\left[\mathbf{f}_{\sigma
O}(\mathbf{k})\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{OO^{\prime}}(\mathbf{k})+\mathbf{G}^{\sigma\sigma^{\prime}\mathcal{C}(1)}_{OO^{\prime}}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}O^{\prime}}(\mathbf{k})\right]\left[\mathbf{U}^{\sigma^{\prime}\sigma^{\prime\prime}\mathcal{C}^{\prime}(1)}_{O^{\prime}V}(\mathbf{k})\right]^{\dagger}\mathbf{U}^{\sigma^{\prime\prime}\sigma\mathcal{C}^{\prime\prime}(1)}_{VO}(\mathbf{k})\Big{\\}}$

where the prime over the sum on $\mathcal{C}^{\prime\prime}$ means that
$\mathcal{C}^{\prime\prime}=\mathcal{R}$ if $\mathcal{C}=\mathcal{C}^{\prime}$
and, otherwise $\mathcal{C}^{\prime\prime}=\mathcal{I}$. Eqs. (41) to (44)
show that the energy may be determined up to third order, using only the
first-order perturbed wavefunction.

Figure 1: Graphical representation of the atomic structure of the series of
tungsten dichalcogenide 2D layered systems considered in this study to test
the proposed perturbation theory treatment of SOC: WX2 with X = Se, Te, Po,
Lv.

##  III Computational Details

Calculations are performed on tungsten dichalcogenide 2D layers, shown in
Figure 1, which have been chosen based on their previously reported “giant
spin-orbit induced spin-splitting”.Zhu et al. (2011) All calculations are
performed with the PBE0 global hybrid functional.Adamo and Barone (1999) We
did not test other functionals because the previous study on molecular systems
has already shown that the convergence behaviour of the perturbation series is
very similar with five different functionals.Desmarais et al. (2021b) The
present calculations employed the small-core STUTSC potential (for W and Lv)
and both the large- and small-core STUTSC, as well as STUTLC potentials for
Se, Te, Po. For W, the valence basis set was of the form
(6$s$6$p$4$d$2$f$)/[5$s$3$p$4$d$2$f$], being modified starting from the
ecp-60-dhf-SVP set available from the Turbomole package.TUR  For the STUTLC
calculations, the valence basis sets for Se, Te and Po of the form
(5$s$5$p$2$d$)/[3$s$3$p$2$d$] were modified from the ones originally presented
in Ref. Stoll et al., 2002. For the STUTSC calculations, the valence basis
sets for Se, Te and Po are of the form (18$s$15$p$8$d$)/[4$s$3$p$2$d$] or
(18$s$13$p$7$d$)/[4$s$3$p$2$d$]. The valence basis set for Lv is an
uncontracted one of the form (10$s$8$p$7$d$1$f$)/[10$s$8$p$7$d$1$f$]. The full
input decks are available in Crystal format in the ESI.ESI  Reciprocal space
was sampled in a $24\times 24$ Monkhorst-Pack net, with Fermi smearing of
0.001 $E_{h}$. A tolerance of 10-8 $E_{h}$ on the total energy was used as a
convergence criterion for the SCF procedure. The five TOLINTEG parameters that
control truncation of the Coulomb and exact-exchange infinite series were set
to 8 8 8 8 30. The exchange-correlation functional and potential (in their
collinear spin-DFT formulation) were sampled on a direct-space pruned grid
over the unit-cell volume with Lebedev angular and Gauss-Legendre radial
quadratures, employing 99 radial and 1454 angular points (keyword XXLGRID).
The geometries of the layers were initially obtained by cleaving three-atom
thick slabs along the (001) surface of the bulk P63/mmc crystal
structures.Schutte et al. (1987) Then, both the atomic fractional coordinates
and lattice parameters of the layers were fully optimized with analytical
gradients of the total energy for systems periodic in two dimensions, and a
quasi-Newton scheme, using, respectively, the PBE and PBE0 functionals at the
scalar-relativistic 1c-SCF level.Doll (2001); Doll et al. (2001, 2006);
Civalleri et al. (2001) Finally, single-point 2c-SCF and CPKS calculations,
including SOC, were performed on the previously optimized scalar-relativistic
geometries.

##  IV Results and Discussion

In the following, we validate the present pure-state periodic CPKS approach
for SOC against reference periodic 2c-SCF calculations.Desmarais et al.
(2020b) Comparisons are reported on i) total energies, ii) electronic band
structures, and iii) spatial distribution of density variables of spin-current
DFT.

Table 1: Percentage contribution of SOC to the energy with large-core (LC) and small-core (SC) RECPs. Relative percentage differences w.r.t. 2c-SCF values are also reported. |  | $\%^{\text{PT2}}$ | $\%^{\text{PT3}}$ | $\%^{\text{2c}}$ |  | $|\Delta\%_{2}|$ | $|\Delta\%_{3}|$
---|---|---|---|---|---|---|---
LC | $\text{WSe}_{2}$ | $2.26\text{E}^{-2}$ | $2.20\text{E}^{-2}$ | $2.21\text{E}^{-2}$ |  | $2.17$ | $0.50$
$\text{WTe}_{2}$ | $2.97\text{E}^{-2}$ | $2.92\text{E}^{-2}$ | $2.93\text{E}^{-2}$ |  | $1.23$ | $0.51$
$\text{WPo}_{2}$ | $7.97\text{E}^{-2}$ | $8.12\text{E}^{-2}$ | $8.13\text{E}^{-2}$ |  | $1.93$ | $0.13$
SC | $\text{WSe}_{2}$ | $2.69\text{E}^{-3}$ | $2.62\text{E}^{-3}$ | $2.63\text{E}^{-3}$ |  | $2.14$ | $0.46$
$\text{WTe}_{2}$ | $6.18\text{E}^{-3}$ | $6.01\text{E}^{-3}$ | $6.03\text{E}^{-3}$ |  | $2.43$ | $0.34$
$\text{WPo}_{2}$ | $3.37\text{E}^{-2}$ | $3.20\text{E}^{-2}$ | $3.20\text{E}^{-2}$ |  | $5.16$ | $0.05$
$\text{WLv}_{2}$ | $6.57\text{E}^{-1}$ | $7.22\text{E}^{-1}$ | $7.10\text{E}^{-1}$ |  | $7.43$ | $1.77$

###  IV.1 Total Energies

Table 1 reports percentage contributions of SOC to total energies through
second and third order in PT ($\%^{\text{PT2}}$ and $\%^{\text{PT3}}$), as
well as with the 2c-SCF reference method ($\%^{\text{2c}}$). Relative
differences w.r.t. 2c-SCF are also provided ($|\Delta\%_{2}|$ and
$|\Delta\%_{3}|$ through second and third order, respectively). It is found
that third order contributions to the total energy consistently improve the
agreement w.r.t 2c-SCF, as $|\Delta\%_{3}|$ is always smaller than
$|\Delta\%_{2}|$. For $\text{WSe}_{2}$ the error in predicting the SOC
contribution to the energy through second order in PT is
$|\Delta\%_{2}|=2.17\%$ with the LC potential or $2.14\%$ with the SC
potential. These figures are improved to $|\Delta\%_{3}|=0.5\%$ and
$|\Delta\%_{3}|=0.46\%$ through third-order. In fact, an error larger than
$1\%$ through third order, is only found on $\text{WLv}_{2}$, with
$|\Delta\%_{3}|=1.77\%$, being a large improvement over the second-order value
of $|\Delta\%_{2}|=7.43\%$.

###  IV.2 Comparison of PT2′ vs. PT2 Approximations

Before discussing the comparisons of the perturbation theories for one-
electron properties on periodic systems, we first validate the PT2′
approximation for GKS eigenvalues and density variables of SCDFT of Eq. (30)
on the well-studied molecular systems of Part II.Desmarais et al. (2021b)
These are, namely, the halogen diatomic and hydride molecules I2,
I${}_{2}^{-}$, At${}_{2}^{-}$, HI, HAt, HTs. We recall that the PT2′
calculation differs from the full PT2 one by an uncoupled, rather than
coupled, solution of Eq. (21). We refer to Part II for computational details
on the molecular systems.

The GKS eigenvalue spectra are provided in Figs. S1 and S2 of the ESI.
Noticeable differences are only found on the single most difficult case of
At${}_{2}^{-}$, which, as previously noted, would ideally require an ensemble,
rather than pure-state treatment.Desmarais et al. (2021b) For At${}_{2}^{-}$,
in Fig. S2, the 1c-SCF and PT2′ HOMO-LUMO gaps are 2.20 and 2.45 eV, while the
PT2 gap is 2.04 eV and the 2c-SCF one is 2.00 eV. Another difficult case is
represented by the 7$p$ superheavy-element Tennessine hydride HTs system of
Fig. S1. In this case, the 1c-SCF, PT2′, PT2 and 2c-SCF gaps are 7.82, 5.76,
5.66 and 6.12 eV. In all other cases, of Figs. S1 and S2, the PT2′, PT2 and
2c-SCF eigenvalue spectra are essentially matching.

As for density variables of SCDFT, Fig. S3 reports contour maps for the
I${}_{2}^{-}$ molecule. Here small difference between PT2′ and 2c-SCF are only
noticeable on the particle-current density $\mathbf{j}$, and these differences
are resolved with the full PT2 calculation.

###  IV.3 Electronic Band Structure

We here validate the present PT treatment for calculation of GKS band
structures. Predicted band gaps of the tungsten dichalcogenide series are
reported in Table 2, as well as differences w.r.t. 2c-SCF values through first
and second orders in PT ($|\Delta_{1}|$ and $|\Delta_{2^{\prime}}|$). Here the
PT2′ gaps are always found to improve upon PT1, with regards to the comparison
against 2c-SCF. Errors larger than $0.01$ eV through second order are only
found in the single most difficult case of $\text{WLv}_{2}$, where
$|\Delta_{2^{\prime}}|=0.09$ eV. For $\text{WLv}_{2}$, the scalar-relativistic
1c-SCF gap is $E_{g}^{\text{1c}}=1.90$ eV, while the 2c-SCF one is
$E_{g}^{\text{2c}}=0.82$ eV. In this particularly challenging case, first-
order PT wrongly predicts a metallic system, with a gap of $0.00$ eV.
Remarkably, the correct description is recovered with the PT2′ treatment,
where the predicted gap becomes $E_{g}^{\text{PT2}^{\prime}}=0.91$ eV, in good
agreement with the 2c-SCF value.

Table 2: Calculated energy band gaps and relative differences w.r.t. 2c-SCF (eV). |  | $E_{g}^{\text{1c}}$ | $E_{g}^{\text{PT1}}$ | $E_{g}^{\text{PT2}^{\prime}}$ | $E_{g}^{\text{2c}}$ |  | $|\Delta_{1}|$ | $|\Delta_{2^{\prime}}|$
---|---|---|---|---|---|---|---|---
LC | $\text{WSe}_{2}$ | $2.68$ | $2.31$ | $2.34$ | $2.33$ |  | $0.02$ | $0.01$
$\text{WTe}_{2}$ | $2.24$ | $1.80$ | $1.84$ | $1.83$ |  | $0.03$ | $0.01$
$\text{WPo}_{2}$ | $2.05$ | $1.35$ | $1.45$ | $1.44$ |  | $0.09$ | $0.01$
SC | $\text{WSe}_{2}$ | $2.68$ | $2.30$ | $2.34$ | $2.33$ |  | $0.02$ | $0.01$
$\text{WTe}_{2}$ | $2.08$ | $1.61$ | $1.66$ | $1.65$ |  | $0.04$ | $0.01$
$\text{WPo}_{2}$ | $1.83$ | $1.08$ | $1.19$ | $1.18$ |  | $0.11$ | $0.01$
$\text{WLv}_{2}$ | $1.90$ | $0.00$ | $0.91$ | $0.82$ |  | $0.82$ | $0.09$

The full electronic band structures of the tungsten dichalcogenide series, as
calculated with the SC potentials, are provided in Fig. 2. Because these
layers lack an inversion center, the scalar-relativistic bands (in the black
dashed line) are spin-split by inclusion of SOC. As expected, the magnitude of
the spin-splitting increases progressively along the series. In the case of
the two lightest-element systems, being WSe2 and WTe2 in the upper panels, the
PT2′ bands (in magenta) are essentially superimposed on the 2c-SCF ones (in
blue), and visible improvements over the PT1 treatment (in turquoise) can be
noticed, especially for the virtual levels. Going down the series, to the
bottom panels, for WPo2 and WLv2, where SOC is enhanced, differences between
PT1 and PT2′/2c-SCF become more evident. In the case of WPo2, small
differences between PT2′ and 2c-SCF begin to appear, and are again more
apparent for virtual bands. In the single most challenging case of WLv2,
first-order PT appears to have completely failed, as the turquoise PT1 bands
display very large differences w.r.t. 2c-SCF. Here, for WLv2, the improvement
provided by the approximate second-order corrections is remarkable, as the
magenta PT2′ bands show instead a good agreement with the blue 2c-SCF ones.
The occupied WLv2 PT2′ bands nearly match the 2c-SCF ones, but visible
differences remain on virtual bands.

Figure 2: Electronic band structures of the tungsten dichalcogenide 2D
systems, at the scalar-relativistic (black dashed lines) level, and including
SOC, both through first- or approximate second-order PT (PT1 in turquoise and
PT2′ in magenta), or the reference 2c-SCF (in blue). Close-up views on the
occupied and virtual levels are also provided for each system.

###  IV.4 Density Variables of Spin-Current DFT

In the following, we report on calculations that quantify the accuracy of the
CPKS approach to reproduce the effect of SOC on the two-component complex
density matrix, through density variables of SCDFT: i) the electron density
$\rho$, and ii) the spin-current density $\mathbf{J}^{z}$. The density
variables are for calculations employing the LC, rather than SC, potential so
as to present results that are more strongly influenced by contributions
originating from valence orbitals (therefore for WSe2, WTe2 and WPo2).

####  IV.4.1 The Electron Density

Fig. 3 reports on calculations of the effect of SOC on the electron density in
WSe2, WTe2 and WPo2, as determined by different approaches (2c-SCF, and CPKS).
In first order, the PT1 density $\rho_{\textup{PT1}}$ coincides with the
scalar-relativistic density $\rho_{\textup{1c}}$, because the real part of
diagonal blocks of the first-order perturbed density are vanishing, from Eq.
(39). On the other hand, these same blocks of the perturbed density do not
vanish, with the PT2′ approach, so that the $\rho_{\textup{PT2}^{\prime}}$
provides a description of the effect of SOC on the electron density.

We report the quantity
$\Delta\rho_{\textup{SOC}}=\rho_{\textup{2c}}-\rho_{\textup{1c}}$ providing
the total effect of SOC on $\rho$, as well as the difference of CPKS and
2c-SCF values
$\Delta\rho_{\textup{PT2}^{\prime}-2c}=\rho_{\textup{PT2}^{\prime}}-\rho_{\textup{2c}}$.
The figure provides electron density data on two distinct planes, being i) the
plane containing the chalcogen atoms (upper panels) and ii) the plane
containing the tungsten atoms (bottom panels). Generally speaking, the plots
show that $\Delta\rho_{\textup{PT2}^{\prime}-2c}$ is always smaller than
$\Delta\rho_{\textup{SOC}}$, and $\Delta\rho_{\textup{SOC}}$ as well as
$\Delta\rho_{\textup{PT2}^{\prime}-2c}$ are always largest on the plane
containing the tunsten atoms. The maxima of
$\Delta\rho_{\textup{PT2}^{\prime}-2c}$ is smaller than the one for
$\Delta\rho_{\textup{SOC}}$, by a factor of about 3 (chalcogen plane for WSe2,
tunsten plane for WTe2 and WPo2) to about 9 (chalcogen plane for WTe2 and
WPo2) with larger values of $\Delta\rho_{\textup{PT2}^{\prime}-2c}$ in the
core, rather than bonding, region of the atoms. To further quantify the
differences, density profiles in WTe2 are provided along the W-W and W-Te
directions in bottom right panels of Fig. 3. The profiles confirm that
$\Delta\rho_{\textup{PT2}^{\prime}-2c}$ has a peak close to the core of the W
atom (at distance of about 0.5 Å from the W nucleus) and smaller values in the
bonding regions.

Figure 3: Effect of SOC on the electron density $\rho$ of WSe2, WTe2 and WPo2.
For each of the three systems, 2D maps are reported on two planes: the plane
of the chalcogen atoms (top panels) and the plane of the W atoms (bottom
panels). For each system and each plane, two quantities are mapped: the effect
of SOC on the electron density
$\Delta\rho_{\textup{SOC}}=\rho_{\textup{2c}}-\rho_{\textup{1c}}$ (left
panels) and the difference between the PT2′ density and the reference $2c$ one
$\Delta\rho_{\textup{PT2}^{\prime}-2c}=\rho_{\textup{PT2}^{\prime}}-\rho_{\textup{2c}}$
(right panels). The bottom right panels show electron density profiles of WTe2
along W-W (left) and W-Te (right) directions: upper panels report the profile
of the total density $\rho_{\textup{SOC}}$ obtained from the reference $2c$
calculation while bottom panels report profiles of $\Delta\rho_{\textup{SOC}}$
(blue) and $\Delta\rho_{\textup{PT2}^{\prime}-2c}$ (yellow) defined above,
along with profiles of
$\Delta\rho_{\textup{SOC}}^{\textup{PT2}^{\prime}}=\rho_{\textup{PT2}^{\prime}}-\rho_{\textup{1c}}$
(pink).

####  IV.4.2 The Spin Current Density

We now discuss the accuracy of the PT approach to reproduce spin-current
densities of the tungsten dichalcogenide series. In the case of
$\mathbf{J}^{z}=\mathbf{J}^{\uparrow}-\mathbf{J}^{\downarrow}$ the
$z$-component spin-current density represents the local velocity field for
transport of spin-magnetization $m_{z}$ (and similar for $\mathbf{J}^{x}$ and
$\mathbf{J}^{y}$ being the velocity fields for transport of $m_{x}$ and
$m_{y}$). Our calculations reproduce all three $\mathbf{J}^{x}$,
$\mathbf{J}^{y}$ and $\mathbf{J}^{z}$ with similar accuracy, and we therefore
report on results of only $\mathbf{J}^{z}$ for sake of brevity.

The results of the calculations on WSe2, WTe2 and WPo2 with the LC potentials
are reported in Fig. 4. At variance with the electron density $\rho$, the
spin-currents $\mathbf{J}^{x}$, $\mathbf{J}^{y}$ and $\mathbf{J}^{z}$ are not
vanishing in first-order, but are vanishing for scalar-relativistic
calculations (i.e. in zeroth-order), and we correspondingly report three maps
(total 2c-SCF $\mathbf{J}^{z}$ on the left panels,
$\Delta\mathbf{J}^{z}_{\textup{PT1}-2c}=\mathbf{J}^{z}_{\textup{PT1}}-\mathbf{J}^{z}_{\textup{2c}}$
on the middle panels and
$\Delta\mathbf{J}^{z}_{\textup{PT2${}^{\prime}$}-2c}=\mathbf{J}^{z}_{\textup{PT2${}^{\prime}$}}-\mathbf{J}^{z}_{\textup{2c}}$
on the right panels). The maps of $\Delta\mathbf{J}^{z}_{\textup{PT1}-2c}$ and
$\Delta\mathbf{J}^{z}_{\textup{PT2${}^{\prime}$}-2c}$ share a common color
scale. With $\mathbf{J}^{z}$ being a vector field, here the color intensities
represent the vector magnitude $|\mathbf{J}^{z}|$, while the arrow lengths and
directions reflect the vector projection of $\mathbf{J}^{z}$ onto the planes
that contain the chalcogen atoms (upper panels) and the tungsten atoms (lower
panels). The $z$ Cartesian direction is perpendicular to the plane, and the
$\mathbf{J}^{z}$ vectors are found in the figures to form circular orbits in
the planes. The arrows are not visible in the difference maps, indicating that
the CPKS approach provides very accurate orientations for the spin-current
densities. Colors, on the other hand are visible, but are about two (in the
case of WSe2 and WTe2) to about one (in the case of WPo2) orders of magnitude
smaller than the total $\mathbf{J}^{z}$, indicating that also the vector
magnitude $|\mathbf{J}^{z}|$ is well reproduced by all CPKS calculations. The
$\textup{PT2}^{\prime}$ treatment usually improves over the PT1 treatment, as
the maps of $\Delta\mathbf{J}^{z}_{\textup{PT2${}^{\prime}$}-2c}$ usually have
color intensities that are smaller than those of
$\Delta\mathbf{J}^{z}_{\textup{PT1}-2c}$. The only exception is for the map of
$\mathbf{J}^{z}$ in the plane containing the tungsten atoms for the lightest
element system WSe2, where more intense colors are visible in the core region
of tungsten atoms. In this case (W-plane of WSe2) even though the
$\textup{PT2}^{\prime}$ treatment worsens the description of spin-current
densities in inner shells, the differences remain very small (around two
orders of magnitude smaller than the total $\mathbf{J}^{z}$), and the
agreement is improved for the description of $\mathbf{J}^{z}$ in the valence
region, because darker colors are observed in the region separating the
tungsten atoms.

Figure 4: Effect of SOC on the spin-current density $\mathbf{J}^{z}$ of WSe2,
WTe2 and WPo2. The color intensities represent the magnitude
$|\mathbf{J}^{z}|$, while the arrow lengths and directions reflect the vector
projection of $\mathbf{J}^{z}$ onto the plane. For each system and each plane,
three quantities are mapped: the effect of SOC on the spin-current density
$\mathbf{J}^{z}$ (left panels), the difference between the PT1 spin-current
density and the reference $2c$ one
$\Delta\mathbf{J}^{z}_{\textup{PT1}-2c}=\mathbf{J}^{z}_{\textup{PT1}}-\mathbf{J}^{z}_{\textup{2c}}$
(middle panels), and finally
$\Delta\mathbf{J}^{z}_{\textup{PT2${}^{\prime}$}-2c}=\mathbf{J}^{z}_{\textup{PT2${}^{\prime}$}}-\mathbf{J}^{z}_{\textup{2c}}$
(right panels).

##  V Conclusions

A previously proposed molecular non-canonical coupled-perturbed Kohn-Sham
density functional theory (KS-DFT)/Hartree-Fock (HF) treatment for spin-orbit
coupling has been generalized to infinite periodic systems. Explicit
expressions have been provided for the total energy through 3rd-order, which
satisfy the 2N + 1 rule (i.e. requiring only the 1st-order perturbed wave
function for its computation). Satisfaction of the 2N + 1 rule has been
achieved for periodic systems by partitioning the key matrices in terms of
components that transform with even and odd parity upon inversion in
reciprocal space ($\mathbf{k}\to\mathbf{-k}$). Second-order corrections to the
perturbed wave function (and one-electron properties) are calculated at the
uncoupled-perturbed level of theory (and this approximation has been justified
on the well-characterized diatomic halogen series of molecules).

The perturbation-theory approach has been validated for calculating the total
energy, electronic band structure and density variables of spin-current DFT on
the tungsten dichalcogenide hexagonal 2D series (i.e. WSe2, WTe2, WPo2 and
WLv2), including challenging 6p and 7p elements. The computed properties
through second- or third-order match well with those from reference two-
component self-consistent field (2c-SCF) calculations. For total energies,
$E^{(3)}$ is found to consistently improve the agreement against the 2c-SCF
values. For electronic band structures, visible differences w.r.t. 2c-SCF
remain through second-order in only the single-most difficult case of WLv2. As
for density variables of spin-current DFT, the perturbed electron density,
being vanishing in first-order, is the most challenging for the perturbation
theory. Visible differences in the electron densities are, however, found to
be largest close to the core region of atoms and smaller in the valence region
relevant to chemical bonding. Perturbed spin-current densities, on the other
hand, are well-reproduced in all tested cases. Our coupled-perturbed approach
thus provides an accurate alternative to 2c-SCF, with potential for over an
order of magnitude savings in computation times. Moreover, it constitutes a
convenient starting point for future improvements to treat multi-reference
systems, where an ensemble treatment would be necessary.

##  Appendix A Virtual-Occupied Blocks of the Second-Order Matrix of Orbital
Rotations

The simultaneous non-canonical solution of Eqs. (16) and (21) yields, for the
virtual-occupied (VO) blocks:

$\displaystyle\mathbf{U}^{\sigma^{\prime}\sigma\mathcal{R}(2)}_{VO}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{V}^{\sigma^{\prime}(0)}\left(\mathbf{k}\right)-\epsilon_{O}^{\sigma(0)}\left(\mathbf{k}\right)}\Bigg{\\{}\sum_{\mathcal{C}}\sum_{\sigma^{\prime\prime}}\Bigg{[}\sum_{O^{\prime}}\mathbf{U}^{\sigma^{\prime}\sigma^{\prime\prime}\mathcal{C}(1)}_{VO^{\prime}}\left(\mathbf{k}\right)\mathbf{G}^{\sigma^{\prime\prime}\sigma\mathcal{C}(1)}_{O^{\prime}O}\left(\mathbf{k}\right)$
(45a) $\displaystyle-$
$\displaystyle\sum_{V^{\prime}}\mathbf{G}^{\sigma^{\prime}\sigma^{\prime\prime}\mathcal{C}(1)}_{VV^{\prime}}\left(\mathbf{k}\right)\mathbf{U}^{\sigma^{\prime\prime}\sigma\mathcal{C}(1)}_{V^{\prime}O}\left(\mathbf{k}\right)\Bigg{]}-\mathbf{G}^{\sigma^{\prime}\sigma\mathcal{R}(2)}_{VO}\left(\mathbf{k}\right)\Bigg{\\}}=-\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(2)}_{OV}\left(\mathbf{k}\right)\right]^{\ast}$
and:
$\displaystyle\mathbf{U}^{\sigma^{\prime}\sigma\mathcal{I}(2)}_{VO}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{V}^{\sigma^{\prime}(0)}\left(\mathbf{k}\right)-\epsilon_{O}^{\sigma(0)}\left(\mathbf{k}\right)}\Bigg{\\{}\sum_{\mathcal{C}}\sum_{\mathcal{C}^{\prime}\neq\mathcal{C}}\sum_{\sigma^{\prime\prime}}\Bigg{[}\sum_{O^{\prime}}\mathbf{U}^{\sigma^{\prime}\sigma^{\prime\prime}\mathcal{C}(1)}_{VO^{\prime}}\left(\mathbf{k}\right)\mathbf{G}^{\sigma^{\prime\prime}\sigma\mathcal{C}^{\prime}(1)}_{O^{\prime}O}\left(\mathbf{k}\right)$
(45b) $\displaystyle-$
$\displaystyle\sum_{V^{\prime}}\mathbf{G}^{\sigma^{\prime}\sigma^{\prime\prime}\mathcal{C}(1)}_{VV^{\prime}}\left(\mathbf{k}\right)\mathbf{U}^{\sigma^{\prime\prime}\sigma\mathcal{C}^{\prime}(1)}_{V^{\prime}O}\left(\mathbf{k}\right)\Bigg{]}-\mathbf{G}^{\sigma^{\prime}\sigma\mathcal{I}(2)}_{VO}\left(\mathbf{k}\right)\Bigg{\\}}=-\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(2)}_{OV}\left(\mathbf{k}\right)\right]^{\ast}$

To obtain Eqs. (45a) and (45b), we insert Eqs. (25)-(28) into Eq. (21).

##  Appendix B Lagrange-Multipliers in Second Order

Inserting Eqs. (26) (29), (29), (45a) and (45b) into Eq. (21) provides the
following set of non-canonical second-order Lagrange multipliers:

$\displaystyle\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}\mathcal{R}(2)}(\mathbf{k})=\mathbf{G}_{OO^{\prime}}^{\sigma\sigma^{\prime}\mathcal{R}(2)}(\mathbf{k})+\sum_{\mathcal{C}}\sum_{\sigma^{\prime\prime}}\sum_{V}$
$\displaystyle\times\Big{[}\frac{1}{2}\left(\epsilon_{O}^{\sigma(0)}(\mathbf{k})-\epsilon_{O^{\prime}}^{\sigma^{\prime}(0)}(\mathbf{k})\right)\mathbf{U}^{\sigma\sigma^{\prime\prime}\mathcal{C}(1)}_{OV}(\mathbf{k})\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}(1)}_{VO^{\prime}}(\mathbf{k})$
$\displaystyle+\left(\epsilon_{V}^{\sigma^{\prime\prime}(0)}(\mathbf{k})-\epsilon_{O}^{\sigma(0)}(\mathbf{k})\right)\mathbf{U}_{OV}^{\sigma\sigma^{\prime\prime}\mathcal{C}(1)}(\mathbf{k})\mathbf{U}_{VO^{\prime}}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}(1)}(\mathbf{k})\Big{]}\;,$
(46a) and:
$\displaystyle\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}\mathcal{I}(2)}(\mathbf{k})=\mathbf{G}_{OO^{\prime}}^{\sigma\sigma^{\prime}\mathcal{I}(2)}(\mathbf{k})+\sum_{\mathcal{C}}\sum_{\mathcal{C}^{\prime}\neq\mathcal{C}}\sum_{\sigma^{\prime\prime}}\sum_{V}$
$\displaystyle\times\Big{[}\frac{1}{2}\left(\epsilon_{O}^{\sigma(0)}(\mathbf{k})-\epsilon_{O^{\prime}}^{\sigma^{\prime}(0)}(\mathbf{k})\right)\mathbf{U}^{\sigma\sigma^{\prime\prime}\mathcal{C}(1)}_{OV}(\mathbf{k})\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}^{\prime}(1)}_{VO^{\prime}}(\mathbf{k})$
$\displaystyle+\left(\epsilon_{V}^{\sigma^{\prime\prime}(0)}(\mathbf{k})-\epsilon_{O}^{\sigma(0)}(\mathbf{k})\right)\mathbf{U}_{OV}^{\sigma\sigma^{\prime\prime}\mathcal{C}(1)}(\mathbf{k})\mathbf{U}_{VO^{\prime}}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}^{\prime}(1)}(\mathbf{k})\Big{]}\;,$
(46b)

where we have introduced:

$\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}(2)}(\mathbf{k})=\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}\mathcal{R}(2)}(\mathbf{k})+\mathbf{E}_{OO^{\prime}}^{\sigma\sigma^{\prime}\mathcal{I}(2)}(\mathbf{k})$
(47)

Exactly analogous expressions can also be worked out for
$\mathbf{E}_{VV^{\prime}}^{\sigma\sigma^{\prime}(2)}(\mathbf{k})$.

##  Appendix C Direct-Space Perturbed Density Matrix in First and Second
Orders

Inserting Eq. (37b) into Eq. (38) provides:

$\displaystyle\mathbf{P}^{\sigma\sigma^{\prime}(1)}(\mathbf{g})$
$\displaystyle=$ $\displaystyle\frac{1}{\Omega}\int^{\prime}d\mathbf{k}\
\Big{\\{}e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{f}_{\sigma}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}(1)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}$
(48) $\displaystyle+$ $\displaystyle
e^{-\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(-\mathbf{k})\mathbf{f}_{\sigma}(-\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}(1)}(-\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(-\mathbf{k})\right]^{\dagger}$
$\displaystyle+$ $\displaystyle
e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}(1)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}}(\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}+e^{-\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(-\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}(1)}(-\mathbf{k})\mathbf{f}_{\sigma^{\prime}}(-\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(-\mathbf{k})\right]^{\dagger}\Big{\\}}$

Now, to relate elements at $\mathbf{k}$ and $-\mathbf{k}$, we also require:

$\mathbf{f}_{\sigma}(\mathbf{k})=\mathbf{f}_{\sigma}(-\mathbf{k})$ (49)

which follows directly from Eq. (32b) since, in pure-state GKS-DFT
calculations, degenerate bands must have the same occupation. Then, inserting
Eqs. (32a), (32b) and (35) into Eq. (48) gives, after combining complex-
conjugates:

$\displaystyle\mathbf{P}^{\sigma\sigma^{\prime}(1)}(\mathbf{g})$
$\displaystyle=$ $\displaystyle\frac{2}{\Omega}\int^{\prime}d\mathbf{k}\
\Big{\\{}\Re\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}_{O}(\mathbf{k})\mathbf{f}_{\sigma
O}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(1)}_{OV}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}_{V}(\mathbf{k})\right]^{\dagger}\right)\
$ (50) $\displaystyle+$
$\displaystyle\imath\Im\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\
\mathbf{C}^{\sigma\sigma(0)}_{O}(\mathbf{k})\mathbf{f}_{\sigma
O}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(1)}_{OV}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}_{V}(\mathbf{k})\right]^{\dagger}\right)\
$ $\displaystyle+$
$\displaystyle\Re\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}_{V}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(1)}_{VO}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}O}(\mathbf{k})\
\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}_{O}(\mathbf{k})\right]^{\dagger}\right)+\i\Im\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}_{V}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(1)}_{VO}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}O}(\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}_{O}(\mathbf{k})\right]^{\dagger}\right)\Big{\\}}$

In second order, an expansion of the CO coefficients of Eq. (36) using Eq.
(14) provides:

$\displaystyle\mathbf{P}^{\sigma\sigma^{\prime}(2)}(\mathbf{k})$
$\displaystyle=$
$\displaystyle\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{f}_{\sigma}(\mathbf{k})\
\left[\mathbf{U}^{\sigma\sigma^{\prime}(2)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}+\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}(2)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}}(\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\
$ (51) $\displaystyle+$
$\displaystyle\sum_{\sigma^{\prime\prime}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime\prime}(1)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime\prime}}(\mathbf{k})\left[\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}(1)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}$

Then, proceeding as in Eqs. (48) and (50), using Eqs. (29), (29), as well as
Eqs. (45a) and (45b), provides:

$\displaystyle\mathbf{P}^{\sigma\sigma^{\prime}(2)}(\mathbf{g})$
$\displaystyle=$
$\displaystyle\frac{2}{\Omega}\int^{\prime}d\mathbf{k}\Big{\\{}\Re\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{f}_{\sigma}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(2)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\right)$
(52) $\displaystyle+$
$\displaystyle\imath\Im\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{f}_{\sigma}(\mathbf{k})\left[\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(2)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\right)\
$ $\displaystyle+$
$\displaystyle\Re\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{R}(2)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}}(\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\right)+\imath\Im\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime}\mathcal{I}(2)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime}}(\mathbf{k})\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\right)$
$\displaystyle+$
$\displaystyle\sum_{\mathcal{C}}\sum_{\sigma^{\prime\prime}}\Re\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime\prime}\mathcal{C}(1)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime\prime}}(\mathbf{k})\left[\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}(1)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\right)$
$\displaystyle+$
$\displaystyle\sum_{\mathcal{C}}\sum_{\mathcal{C}^{\prime}\neq\mathcal{C}}\sum_{\sigma^{\prime\prime}}\imath\Im\left(e^{\imath\mathbf{k}\cdot\mathbf{g}}\mathbf{C}^{\sigma\sigma(0)}(\mathbf{k})\mathbf{U}^{\sigma\sigma^{\prime\prime}\mathcal{C}(1)}(\mathbf{k})\mathbf{f}_{\sigma^{\prime\prime}}(\mathbf{k})\left[\mathbf{U}^{\sigma^{\prime\prime}\sigma^{\prime}\mathcal{C}^{\prime}(1)}(\mathbf{k})\right]^{\dagger}\left[\mathbf{C}^{\sigma^{\prime}\sigma^{\prime}(0)}(\mathbf{k})\right]^{\dagger}\right)\Big{\\}}$

## References

  * Chang et al. (1986) C. Chang, M. Pelissier, and P. Durand, Phys. Scr. 34, 394 (1986).
  * Buenker et al. (1984) R. J. Buenker, P. Chandra, and B. A. Hess, Chem. Phys. 84, 1 (1984).
  * Almlöf et al. (1985) J. Almlöf, K. Faegri Jr, and H. Grelland, Chem. Phys. Lett. 114, 53 (1985).
  * Douglas and Kroll (1974) M. Douglas and N. M. Kroll, Ann. Phys. 82, 89 (1974).
  * Dyall (1997) K. G. Dyall, J. Chem. Phys. 106, 9618 (1997).
  * Iliaš et al. (2005) M. Iliaš, H. J. A. Jensen, V. Kellö, B. O. Roos, and M. Urban, Chem. Phys. Lett. 408, 210 (2005).
  * Kutzelnigg and Liu* (2006) W. Kutzelnigg and W. Liu*, Mol. Phys. 104, 2225 (2006).
  * Liu and Kutzelnigg (2007) W. Liu and W. Kutzelnigg, J. Chem. Phys. 126, 114107 (2007).
  * Liu and Peng (2006) W. Liu and D. Peng, J. Chem. Phys. 125, 044102 (2006).
  * Liu and Daoling (2009) W. Liu and P. Daoling, J. Chem. Phys. 131 (2009).
  * Peng et al. (2007) D. Peng, W. Liu, Y. Xiao, and L. Cheng, J. Chem. Phys. 127, 104106 (2007).
  * Sikkema et al. (2009) J. Sikkema, L. Visscher, T. Saue, and M. Iliaš, J. Chem. Phys. 131, 124116 (2009).
  * Seino and Hada (2008) J. Seino and M. Hada, Chem. Phys. Lett. 461, 327 (2008).
  * Seino and Nakai (2012) J. Seino and H. Nakai, J. Chem. Phys. 136, 244102 (2012).
  * Saue and Jensen (1999) T. Saue and H. A. Jensen, J. Chem. Phys. 111, 6211 (1999).
  * Jorgen Aa. Jensen et al. (1996) H. Jorgen Aa. Jensen, K. G. Dyall, T. Saue, and K. Fægri Jr, J. Chem. Phys. 104, 4083 (1996).
  * Rutkowski (1986a) A. Rutkowski, J. Phys. B 19, 149 (1986a).
  * Rutkowski (1986b) A. Rutkowski, J. Phys. B 19, 3431 (1986b).
  * Rutkowski (1986c) A. Rutkowski, J. Phys. B 19, 3443 (1986c).
  * Kutzelnigg (1996) W. Kutzelnigg, Phys. Rev. A 54, 1183 (1996).
  * van Wüllen (1995) C. van Wüllen, J. Chem. Phys. 103, 3589 (1995).
  * Stopkowicz and Gauss (2011) S. Stopkowicz and J. Gauss, J. Chem. Phys. 134, 064114 (2011).
  * Stopkowicz and Gauss (2019) S. Stopkowicz and J. Gauss, Mol. Phys. 117, 1242 (2019).
  * Cheng et al. (2013) L. Cheng, S. Stopkowicz, and J. Gauss, J. Chem. Phys. 139, 214114 (2013).
  * Cheng et al. (2018) L. Cheng, F. Wang, J. F. Stanton, and J. Gauss, J. Chem. Phys. 148, 044108 (2018).
  * Desmarais et al. (2021a) J. K. Desmarais, A. Erba, J.-P. Flament, and B. Kirtman, J. Chem. Theor. Comput. 17, 4697 (2021a).
  * Desmarais et al. (2021b) J. K. Desmarais, A. Erba, J.-P. Flament, and B. Kirtman, J. Chem. Theor. Comput. 17, 4712 (2021b).
  * Dolg and Cao (2012) M. Dolg and X. Cao, Chem. Revs. 112, 403 (2012).
  * Ermler et al. (1981) W. C. Ermler, Y. S. Lee, P. A. Christiansen, and K. S. Pitzer, Chem. Phys. Lett. 81, 70 (1981).
  * Dovesi et al. (2018) R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima, S. Salustro, et al., WIREs Comput. Mol. Sci. 8, e1360 (2018).
  * Erba et al. (2022) A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, et al., J. Chem. Theory Comput. (2022).
  * Lieb (1983) E. H. Lieb, Int. J. Quantum Chem 24, 243 (1983).
  * Gross et al. (1988a) E. K. Gross, L. N. Oliveira, and W. Kohn, Phys. Rev. A 37, 2805 (1988a).
  * Gross et al. (1988b) E. K. Gross, L. N. Oliveira, and W. Kohn, Phys. Rev. A 37, 2809 (1988b).
  * Desmarais and Kirtman (2022) J. K. Desmarais and B. Kirtman, unpublished (2022).
  * Helgaker and Teale (2022) T. Helgaker and A. M. Teale, arXiv preprint arXiv:2204.12216 (2022).
  * Von Barth and Hedin (1972) U. Von Barth and L. Hedin, J. Phys. C Solid State Phys. 5, 1629 (1972).
  * Vignale and Rasolt (1988) G. Vignale and M. Rasolt, Phys. Rev. B 37, 10685 (1988).
  * Bencheikh (2003) K. Bencheikh, J. of Phys. A 36, 11929 (2003).
  * Vignale and Rasolt (1987) G. Vignale and M. Rasolt, Phys. Rev. Lett. 59, 2360 (1987).
  * Desmarais et al. (2019a) J. K. Desmarais, J.-P. Flament, and A. Erba, J. Phys. Chem. Lett. 10, 3580 (2019a).
  * Desmarais et al. (2019b) J. K. Desmarais, J.-P. Flament, and A. Erba, J. Chem. Phys. 151, 074107 (2019b).
  * Desmarais et al. (2021c) J. K. Desmarais, S. Komorovsky, J.-P. Flament, and A. Erba, J. Chem. Phys. 154, 204110 (2021c).
  * Desmarais et al. (2020a) J. K. Desmarais, J.-P. Flament, and A. Erba, Phys. Rev. B 102, 235118 (2020a).
  * Bodo et al. (2022) F. Bodo, J. K. Desmarais, and A. Erba, Phys. Rev. B 105, 125108 (2022).
  * Desmarais et al. (2020b) J. K. Desmarais, J.-P. Flament, and A. Erba, Phys. Rev. B 101, 235142 (2020b).
  * Comaskey et al. (2022) W. P. Comaskey, F. Bodo, A. Erba, J. L. Mendoza-Cortes, and J. K. Desmarais, Phys. Rev. B accepted (2022).
  * (48) _See electronic supporting information at (link) for the full derivation of eqs. (40), (42) and (43). full input decks for the calculations on the tungsten dichalcogenides, figs. s1-s3 on validation of the pt2’ approximation, and table s1 with raw energy data on the tungsten dichalcogenides._
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * (50) _TURBOMOLE V7.2 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from
http://www.turbomole.com._

  * Stoll et al. (2002) H. Stoll, B. Metz, and M. Dolg, J. Comput. Chem. 23, 767 (2002).
  * Schutte et al. (1987) W. Schutte, J. De Boer, and F. Jellinek, Journal of Solid State Chemistry 70, 207 (1987).
  * Doll (2001) K. Doll, Comput. Phys. Commun. 137, 74 (2001).
  * Doll et al. (2001) K. Doll, V. Saunders, and N. Harrison, Int. J. Quantum Chem. 82, 1 (2001).
  * Doll et al. (2006) K. Doll, R. Dovesi, and R. Orlando, Theor. Chem. Acc. 115, 354 (2006).
  * Civalleri et al. (2001) B. Civalleri, Ph. D’Arco, R. Orlando, V. R. Saunders, and R. Dovesi, Chem. Phys. Lett. 348, 131 (2001).
  * Zhu et al. (2011) Z. Zhu, Y. Cheng, and U. Schwingenschlögl, Physical Review B 84, 153402 (2011).