# Meta-generalized-gradient exchange-correlation functionals for spin-orbit
coupled electrons

Jacques K. Desmarais<EMAIL_ADDRESS>Dipartimento di
Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy  Jefferson
Maul  Dipartimento di Chimica, Università di Torino, via Giuria 5, 10125
Torino, Italy  Bartolomeo Civalleri  Dipartimento di Chimica, Università di
Torino, via Giuria 5, 10125 Torino, Italy  Alessandro Erba  Dipartimento di
Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy  Giovanni
Vignale  Institute for Functional Intelligent Materials, National University
of Singapore, 4 Science Drive 2, Singapore 117544  Stefano Pittalis
<EMAIL_ADDRESS>Istituto Nanoscienze, Consiglio Nazionale delle
Ricerche, Via Campi 213A, I-41125 Modena, Italy

###### Abstract

The prominence of density functional theory (DFT) in the field of electronic
structure computation stems from its ability to usefully balance accuracy and
computational effort. At the base of this ability is a functional of the
electron density: the exchange-correlation (xc) energy. This functional
satisfies known exact conditions that guide the derivation of approximations.
The strongly-constrained-appropriately-normed (SCAN) approximation stands out
as a successful, modern, example. In the presence of spin-orbit coupling,
however, the xc energy must depend not only on the electron density but also
on spin currents – a fact which is not captured by SCAN or by any mainstream
functional approximation. In this work, we demonstrate how the SU(2) gauge-
invariance of the exact functional allows us to add the necessary dependence
on spin currents in the SCAN functional (here called JSCAN) — and similar
meta-generalized-gradient functional approximations — in a non-empirical
manner. The prescription further ensures that the conjugate xc-fields of the
spin currents have vanishing local torques as in the exact case. In passing, a
spin-current dependent generalization of the electron localization function
(here called JELF) is also derived. The extended forms are implemented in a
developer’s version of the Crystal23 program. Applications on molecules and
materials confirm the practical relevance of the extensions.

###### pacs:

71.15.Mb, 71.15Rf, 31.15.E-

Introduction. Relativistic effects play an important role in our understanding
of the electronic structure of matter. Greiner (2012) Spin-orbit coupling
(SOC), in particular, emerges naturally in the two-component Pauli equation (a
simplified low-energy version of the Dirac equation) and plays a crucial role
in our understanding of the band structure of materials such as Si, Ge and
GaAs, which powered the semiconductor revolution. Yu and Cardona (2010) More
recently, the interest in spin-orbit interactions in condensed matter has
grown even further following the discovery of a host of materials which
exhibit nontrivial topological properties Bernevig and Hughes (2013);
Vanderbilt (2018) and the explosion of interest in new techniques to
accurately control the spin of the electrons in spintronic and quantum
computing devices.$\bm{Z}$utić et al. (2004)

Just as an electric potential couples to the particle density and a magnetic
field couples to the spin and to the particle (i.e. ordinary) current, SOC can
be described as a tensor field which couples to the spin current (SC), i.e.,
the tensor product of spin and ordinary current. Unlike the ordinary current,
the spin current is invariant under time reversal (being the product of two
time-reversal odd quantities) and therefore can freely flow in the time-
reversal invariant ground state of a non magnetic system.Tokatly (2008);
Droghetti et al. (2022)

It is sometimes adequate (for systems composed of light atoms) to include SOC
as a small perturbation on top of the Kohn-Sham (KS) potential of ordinary
density functional theory (DFT).Kohn (1999) This approach can in principle
find the particle density and total energy of real interacting ground states
by solving a single-particle problem in the presence of a self-consistent
exchange-correlation (xc) potential, which is a functional of the density.
Also, the eigenvalues of the KS single-particle equation provide a useful
representation of the band structure and related properties.Seidl et al.
(1996); Perdew et al. (2017); Nazarov and Vignale (2011); Prokopiou et al.
(2022)

However, there is no reason why a density functional approximation that was
developed for systems lacking SOC would still yield accurate results when
strong SOC are present. To accomplish this more ambitious task, we need an
exchange-correlation energy functional that depends not only on the density
but also — and crucially — on the spin current density, and does so in a
universal (i.e., system-independent) manner, as mandated by the original spin-
current density functional theory (SCDFT)Vignale and Rasolt (1988); Bencheikh
(2003) and its more recent developments. Rohra and Görling (2006); Heaton-
Burgess et al. (2007); Trushin and Görling (2018); Desmarais et al. (2023a) In
practice, this means including in the effective potential, not only the “bare”
SOC, but also its “dressing” due to many-body effects, in analogy to the
dressing of the bare electron-nuclear potential by Hartree and exchange-
correlation (xc) effects.

Crucially, any xc energy functional that depends on the spin current density
must be invariant — not only under U(1) gauge transformationsVignale and
Rasolt (1987); Tao and Perdew (2005) — but also under SU(2)
transformationsVignale and Rasolt (1988); Bencheikh (2003); Pittalis et al.
(2017). The latter can be viewed as local rotations of the spinorial state of
the electronsFröhlich and Studer (1993); Jin et al. (2006). Failing to satisfy
SU(2) invariance would have disastrous consequences for a density functional
approximation: it would mean that the approximation does not distinguish
between actual variations of the spin current density — which can have
physical consequences — and variations arising from gauge transformations,
which must have no consequence! Thus the gauge invariance of the exact xc-
energy functional acts as a most important exact constraint, which guides the
construction of approximate functionals.

While different approximate functionals obviously have different strengths and
weaknesses, _they all must satisfy the constraint of gauge-invariance._ Hence,
we focus on the strongly-constrained-appropriately-normed (SCAN)
functional,Sun et al. (2015) which is currently recognized as the least
empirical general-purpose approximation in the class of Meta-Generalized-
Gradient Approximations (MGGA).Perdew et al. (2009); Sun et al. (2013, 2015,
2016); Tran et al. (2016); Peng et al. (2016); Brandenburg et al. (2016);
Isaacs and Wolverton (2018); Tran et al. (2019) However, our procedure is
applicable to any approximate functional in this class.

The exceptional performance of the SCAN is fundamentally due to its ability to
satisfy altogether a large number of known exact conditions for the xc-energy
functional of DFT. Furthermore, among its crucial ingredients, the KS kinetic-
energy density appears in a combination of variables that was originally
proposed for the so-called electron localization function (ELF) – a quantity
that provides a vivid visualization of molecular bonds and atomic shells (more
below).Becke and Edgecombe (1990a); Silvi and Savin (1994)

But DFT (and thus SCAN) cannot deal with non-perturbative SOC! This brings
forth the necessity to switch to SCDFT both in principleVignale and Rasolt
(1988); Bencheikh (2003) and in practice.Rohra and Görling (2006); Heaton-
Burgess et al. (2007); Trushin and Görling (2018); Desmarais et al. (2023a,
2020a); Bodo et al. (2022); Comaskey et al. (2022) The dependence of the xc
energy functional on spin currents in SCDFT is strongly constrained by the
requirement that its form remains invariant under local SU(2) gauge
transformations (defined below). Taking as a starting point the SCAN and the
ELF forms, in this work we show how to make them functionals of the spin
currents while preserving the fundamental SU(2) gauge invariance of the
theory. The resulting generalizations are here called JELF and JSCAN. A set of
computed results for molecules and materials demonstrates the practical
usefulness of our non-empirical extension of the SCAN functional. An
implementation is presented in a developer’s version of the Crystal23
program.Erba et al. (2023)

SCDFT in a nutshell. We shall consider ground states at geometries that do not
magnetize nor host particle currents spontaneously. Furthermore, magnetic
couplings (Zeeman field, Abelian vector potential) that break time-reversal
symmetry are not included. Accounting for SOC, however, we must consider
together with the particle density, $n$, the spin currents, ${\vec{\bf J}}$,
and corresponding external couplings: a scalar, $v$, and a non-Abelian vector
potential, ${\vec{\bf A}}$. The former accounts for the electron-nuclear
potential, while the latter provides a useful formal representation of
SOC.Fröhlich and Studer (1993); Bencheikh (2003); Abedinpour et al. (2010);
Zhang and Vignale (2016); Droghetti et al. (2022) The energy density
functional of interest, thus, reads as follows:

$\displaystyle E_{v,{\vec{\bf A}}}[n,{\vec{\bf J}}]=T_{s}[n,{\vec{\bf
J}}]+E_{{\rm H}}[n]+E_{\rm xc}[n,{\vec{\bf J}}]+\frac{1}{c}\int d^{3}r~{}{\bf
J}^{a}\cdot{\bf A}^{a}+\int d^{3}r~{}{n}\left(v+\frac{1}{2c^{2}}{\bf
A}^{a}\cdot{\bf A}^{a}\right)$ (1)

where $T_{s}$ is the kinetic energy of the KS system, while $E_{{\rm H}}$ and
$E_{\rm xc}$ are Hartree and xc energies. Here and in the following, we denote
with bold characters, $\bm{v}$, quantities with spatial indices (Greek lower
indices, $v_{\mu}$, when written explicitly); and use an arrow, $\vec{v}$, to
denote quantities with spin indices (upper Latin indices, $v^{a}$, when
written explicitly) — $\mu$ and $a$ have values $x$,$y$,$z$. Therefore
${\vec{\bf J}}$ and ${\vec{\bf A}}$ have both spatial and spin indices.
Contractions over spatial indices are denoted with “$\cdot$” and Einstein
convention is used for summing over repeated indices. Unless otherwise stated,
we use Hartree atomic units in which $\hbar=m=1$. The corresponding KS
equations are:

$\displaystyle\left\\{\frac{1}{2}\left[-i\nabla+\frac{1}{c}{\mathbfcal{A}}_{\rm{\tiny\tiny
GKS}}\right]^{2}+\left[\left(v+v_{\rm
Hxc}\right)+\frac{1}{2c^{2}}\left(\mathbfcal{A}^{2}-{\mathbfcal
A}^{2}_{\rm{\tiny\tiny
GKS}}\right)\right]\right\\}\Phi_{\kappa}=\varepsilon_{\kappa}\Phi_{\kappa}\;$
(2)

where ${\mathbfcal{A}_{\rm{\tiny\tiny
GKS}}}={\mathbfcal{A}}+{\mathbfcal{A}_{\rm xc}}={\sigma}^{a}{{\bf
A}}^{a}+{\sigma}^{a}{{\bf A}}_{\rm xc}^{a}$, with $\vec{\sigma}$ being the
vector of Pauli matrices, $v_{\rm Hxc}=\delta E_{\rm Hxc}/\delta n$ is a
Hartree-xc-scalar potential (analogous to the Hxc-scalar potential of DFT) and
${\frac{1}{c}}{{\vec{\bf A}}}_{\rm xc}=\delta E_{\rm xc}/\delta{\vec{\bf J}}$
is a non-Abelian xc-vector potential (which is specific to SCDFT). The
occupied (lowest in energy) two-component KS spinors $\Phi_{k}$ allow us to
compute the particle density,
$n=\sum_{k=1}^{N}\Phi^{\dagger}_{k}\Phi_{k}^{\vphantom{\dagger}}$, and the
spin current ${\vec{\bf
J}}=\frac{1}{2i}\sum_{k=1}^{N}\Phi^{\dagger}_{k}\vec{\sigma}\left[\nabla\Phi_{k}^{\vphantom{\dagger}}\right]-\left[\nabla\Phi^{\dagger}_{k}\right]\vec{\sigma}\Phi_{k}^{\vphantom{\dagger}}$
of the real state and, thus, the energy as well. All this is exact in
principle, while, in practice, we are bound to approximate $E_{\rm xc}$.

Invoking a MGGA, it is convenient to consider the so-called generalized-KS
approach to SCDFT.Desmarais et al. (2023a) This means that for generating the
single-particle equations, we differentiate the energy functional expression
explicitly w.r.t. the orbitals rather than w.r.t. to the densities not (a) —
nowadays, this has become a standard procedure in DFT.Neumann et al. (1996)

A central property of the xc-energy functional in SCDFT is its form invariance
under local SU(2) transformations, i.e., transformations that act on the KS
spinors as follows: $\Phi({\bf r})\rightarrow U_{S}({\bf r})\hat{\Psi}({\bf
r})$, where $U_{S}({\bf
r})=\exp\big{[}\frac{i}{c}{\sigma}^{a}{\lambda}^{a}({\bf r})\big{]}$, i.e.,
subjecting them to different rotations at different points in space. This
invariance is expressed by the equation

$\displaystyle E_{\rm xc}[n^{\prime},{{\vec{\bf J}}^{\prime}}]=E_{\rm
xc}[n,{\vec{\bf J}}]$ (3)

where $n\to n^{\prime}=n$, ${\bf J}^{a}\to{{\bf
J}^{\prime}}^{a}=R^{ab}\big{[}{\bf J}^{b}+\frac{n}{c}{\bf G}^{b}\big{]}$,
$R^{ab}$ is a $3\times 3$ matrix describing a rotation in $\mathbb{R}^{3}$
around the direction $\hat{\lambda}$ by an angle $\varphi=-2\lambda/c$, and
${\vec{\bf G}}=-\frac{ic}{2}{\rm Tr}\left(\vec{\sigma}U^{\dagger}_{\rm
S}\nabla U_{\rm S}^{\vphantom{\dagger}}\right)$. not (b)

It is important to notice that the spin current itself is not a gauge
invariant quantity since it can appear as a consequence of a generic SU(2)
transformation acting on a state in which ${\vec{\bf J}}=0$. The invariance of
$E_{\rm xc}$ stated in Eq. (3), implies that only certain gauge-invariant
combinations of ${\vec{\bf J}}$ with other KS quantities can legitimately
appear in $E_{xc}$. We now tackle the task of generalizing the SCAN functional
(and the closely related ELF) to include a dependence on spin currents that
complies with Eq. (3).

Extension of SCAN and ELF to spin-current carrying states. The SCAN functional
setsSun et al. (2015) $E_{\mathrm{xc}}=\int d^{3}r~{}n({\bf
r})\epsilon^{\mathrm{SCAN}}_{\mathrm{xc}}({\bf r})$ where
$\epsilon_{\mathrm{xc}}^{\mathrm{SCAN}}=\epsilon_{\mathrm{xc}}^{1}+\left(\epsilon_{\mathrm{xc}}^{0}-\epsilon_{\mathrm{xc}}^{1}\right)f_{\mathrm{xc}}(\alpha)$
interpolates between the semi-local energy densities for single-orbital
densities ($\epsilon_{\mathrm{xc}}^{0}$) and for slowly-varying densities
($\epsilon_{\mathrm{xc}}^{1}$), which only depend on $n$ (a gauge-invariant
quantity). The interpolation function $f_{\mathrm{xc}}$ is controlled by the
variable $\alpha=(\tau-\tau_{\rm W})/\tau_{\rm unif}$. The latter quantity
entails three kinetic energy densities: the positive-definite conventional
$\tau=1/2\sum_{i}^{\mathrm{occ.}}|\nabla\phi_{i}|^{2}$ defined in terms of the
occupied (one-component) KS orbitals {$\phi_{i}$}, the von Weizsäcker
$\tau_{W}=|\nabla n|^{2}/(8n)$ — that is the bosonic expression of $\tau$ —
where $n=\sum_{i}^{\mathrm{occ.}}|\phi_{i}|^{2}$, and the kinetic energy of
the noninteracting Fermi gas at uniform density $n$: $\tau_{\rm
unif}=\frac{3}{10}(3\pi^{2})^{2/3}\,n^{5/3}$.

The success of the SCAN is explained not only by the fact that it satisfies a
large set of known exact conditions but also by the fact that it uses,
specifically, $\alpha$ as a key variable. The $\alpha$ tends to zero in the
iso-orbital limit (i.e., when the density is dominated by a single occupied
orbital), tends to one in the uniform-density limit, and tends to infinity in
regions dominated by density overlap between closed shells. In fact, $\alpha$
is a useful building block of density functional approximations,Becke (1996a);
Perdew et al. (2009); Sun et al. (2013) and it is the main ingredient of the
electron localization function,Becke and Edgecombe (1990b); Silvi and Savin
(1994); Savin et al. (1997) ${\rm ELF}=1/(1+\alpha^{2})$ – a simple but
effective descriptor of molecular bonds and atomic shells.

More in general, the kinetic energy density must be expressed in terms of the
occupied two-component KS spinors
$\tau=1/2\sum_{k}^{\mathrm{occ.}}\nabla\Phi^{{\dagger}}_{k}\cdot\nabla\Phi_{k}^{\vphantom{\dagger}}$.
Crucially, let us next consider the transformation of the kinetic energy under
a local SU(2) transformation, $U_{\rm S}$, of the spinors. For the states
considered in this work, one readily finds: Pittalis et al. (2017)
$\tau\to\tau^{\prime}=\tau+\frac{1}{c}{\bf J}^{a}\cdot{\bf
G}^{a}+\frac{n}{2c^{2}}{\bf G}^{a}\cdot{\bf G}^{a}$; for the definition of
${\vec{\bf G}}$ see below Eq. (3). As a consequence of which, it is apparent
that $\alpha$ will change, too. Hence, neither the SCAN functional nor the ELF
are SU(2)-form invariant. Such an invariance, however, is not only an exact
mathematical property of the xc-energy functional in SCDFT but — because a
local SU(2) transformation may be regarded as part of an overall gauge
transformation (i.e., a change in the description which cannot change the
physics) — it is also a proper feature for any quantity that should carry
direct chemical-physical information. Fortunately, in one stroke, the
substitution $\tau\to\widetilde{\tau}=\tau-\frac{{\bf J}^{a}\cdot{\bf
J}^{a}}{2n}$, solves both issues: i.e., it introduces an explicit dependence
on the spin currents and enforces the sought form invariance. In detail,

$\alpha\rightarrow\widetilde{\alpha}=\alpha-\frac{{\bf J}^{a}\cdot{\bf
J}^{a}}{2n\tau_{\rm unif}}\;$ (4)

implies

$\epsilon_{\mathrm{xc}}^{\mathrm{JSCAN}}:=\epsilon_{\mathrm{xc}}^{1}+\left(\epsilon_{\mathrm{xc}}^{0}-\epsilon_{\mathrm{xc}}^{1}\right)f_{\mathrm{xc}}(\widetilde{\alpha}),$
(5)

and

${{\rm JELF}}:=\frac{1}{1+\widetilde{\alpha}^{2}}\;.$ (6)

Note, the above substitution should not be confused with the analogous well-
known substitution that involves the (paramagnetic) particle current:Tao and
Perdew (2005) $\tau\to\widetilde{\tau}=\tau-\frac{{\bf j}\cdot{\bf j}}{2n}.$
The latter derives from the consideration of local U(1) transformations. The
importance of the dependence of density functionals on the particle current
has been largely demonstrated. Dobson (1993); Becke (1996b, 2002); Maximoff et
al. (2004); Burnus et al. (2005); Tao and Perdew (2005); Pittalis et al.
(2007, 2009); Räsänen et al. (2009); Oliveira et al. (2010); Zhu et al.
(2016); James W. Furness and Teale (2016); Holzer et al. (2022)

Next, from Eq. (5), taking a functional derivative with respect to ${\vec{\bf
J}}$, we obtain

$\frac{1}{c}{{\vec{\bf A}}}^{{}_{{\rm JSCAN}}}_{\mathrm{xc}}=-{\vec{\bf
J}}\left[\frac{\left(\epsilon_{\mathrm{xc}}^{0}-\epsilon_{\mathrm{xc}}^{1}\right)f^{\prime}_{\mathrm{xc}}(\widetilde{\alpha})}{n\tau_{\rm
unif}}\right]\;,$ (7)

where $f^{\prime}_{\mathrm{xc}}$ denotes the derivative of $f_{xc}$ wrt its
own argument. This is a non-Abelian effective connection; i.e., an xc-spin-
vector potential. This expression shows that ${{\vec{\bf A}}}^{{}_{{\rm
JSCAN}}}_{\mathrm{xc}}$ is parallel to ${\vec{\bf J}}$ at each point in space
and, therefore,

$\frac{1}{c}{{\vec{A}}}^{{}_{{\rm JSCAN}}}_{{\rm
xc},\mu}\times\vec{J}_{\mu}=0\;.$ (8)

For the non-magnetic ground states considered in this work, the JSCAN
satisfies the non-Abelian torque conditions exactly. not (c)

In summary, Eqs. (4)-(8) are the key equations of this work — Eqs. (4), (7)
and (8) apply to any MGGA which uses $\tau$ through $\alpha$ (for example the
TASK functional).Aschebrock and Kümmel (2019)

Applications. We have implemented a self-consistent treatment of the JSCAN
functional, including effective non-Abelian vector potentials, as well as the
JELF in a developer’s version of the Crystal23 code.Erba et al. (2023)
Computational and implementation details are provided in the supporting
informationESI  (see also Refs. Desmarais et al., 2021a, 2020b, 2019a, 2019b;
Metz et al., 2000; Peterson et al., 2003, 2007; Doll et al., 2001; Doll, 2001;
Doll et al., 2006; Civalleri et al., 2001; Laun et al., 2018; Heyd et al.,
2005; Lebedev, 1976, 1977; Towler et al., 1996; Lehtola et al., 2020 therein).
In the calculations, we employ the restored-regularized r2-SCAN, rather than
the original SCAN functional, to improve numerical stability and
performance.Furness et al. (2022)

Table 1: We report the effect of SOC on the HOMO-LUMO gap at the level of Exchange-only and Exchange-Correlation using J(SCAN) and the SOC effect on the optical gap at the RPAx level.NWC ; Aprà et al. (2020) The results show that these effects are comparable, and the agreement is improved when spin currents are included in the MGGA. All values are reported in eV. | Exx-Only |  | Exx & Corr
---|---|---|---
| J-r2SCAN | HF |  | J-r2SCAN | RPAx
FI | 0.248 (0.225) | 0.322 |  | 0.232 (0.222) | 0.360
HI | 0.337 (0.276) | 0.343 |  | 0.324 (0.281) | 0.343
I2 | 0.333 (0.276) | 0.364 |  | 0.317 (0.278) | 0.396
HAt | 0.927 (0.675) | 0.891 |  | 0.884 (0.213) | 0.887

We apply our SU(2) gauge-invariant J-r2SCAN approach to molecules on
reproducing the SOC-induced modification of the valence electronic structure
in halogen diatoms and hydrides, which have served as model systems in
previous studies.Desmarais et al. (2021b, c, 2023b) Values of the SOC-induced
modification to the HOMO-LUMO gap are reported in Table 1, first at the
exchange-only level of theory, and then including effects of exchange and
correlation. Values are presented with exchange-only (gauge-invariant)
J-r2SCAN functional, and numbers in parentheses represent corresponding
results with the standard (gauge-dependent) r2SCAN approach.

Inclusion of the spin currents in the functional greatly improves the
agreement with Hartree-Fock (HF) results. In the case of the HI molecule, for
example, the HOMO-LUMO gap is increased from 0.276 eV to 0.337 eV; which
reduces the relative error from 20% to 2%. In the case of HAt, the effect of
SOC on the HOMO-LUMO gap is increased from 0.675 eV to 0.927 eV; relative
error decreased from 24% to 4%.

At the xc level of theory, values are benchmarked against the random phase
approximation including the dynamical response of Fock exchange (RPAx).Ullrich
(2011) In detail, this comparison assumes that the effect of SOC on the
fundamental gap is to a good approximation the same as the effect of SOC on
the optical gap. In all cases, inclusion of currents in the xc model improves
agreement against the RPAx. The most extreme case is HAt where the effect of
SOC on the gap is increased from 0.213 eV to 0.884 eV by inclusion of spin-
currents in the MGGA, essentially coinciding with the RPAx value of 0.887 eV.

Next, we discuss electron localization through the JELF [see Eq. (6)]. Let us
begin by recalling that large values of $\tilde{\alpha}$, corresponding to a
(J)ELF close to one, indicate that electrons are highly localized. Here we
look at the effect of SOC on the (J)ELF. For the I2 molecule, after a self-
consistent solution with J-r2SCAN, including both exchange and correlation
effects, in Fig. 1(a) we plot the ELF and JELF indicators as differences
w.r.t. the calculation without SOC. Yellow contours indicate that SOC
localizes the electron. Blue contours indicate that the electron is
delocalized by SOC. The leftmost panel, containing the $\Delta{\rm SOC}$ JELF
shows that electrons are localized along the bond axis by SOC. On the other
hand, we observe two blue “onion rings” around the atomic centers. Thus, SOC
localizes electrons along the bond and delocalizes them in the orthogonal
directions. Although this analysis is consistent with previous
reports,Desmarais et al. (2023c); Gulans and Draxl (2022) the degree of
delocalization is largely overestimated by the gauge-dependent ELF, compared
to the SU(2) form-invariant JELF, see their difference in the rightmost panel
of Fig. 1(a). The JELF displays more localization on the bonding compared to
the (gauge-dependent) ELF.

We proceed to apply our approach to the description of the electronic band
structure of infinite periodic systems. In Table 2, we provide values of the
splitting of the valence band at the $\mathbf{K}$ high-symmetry point in
molybdenum dichalcogenide Rashba systems: the inversion asymmetric hexagonal
single layer MoSe2 and the inversion-symmetric hexagonal bulk $\alpha$-MoTe2.
The valence band structure is provided in panel (b) of Figure 1.

For both systems, about one third of the total effect of SOC on the splitting
is accounted for through many-body effects by inclusion of spin currents in
the functional. Low-temperature experimental values on the transport gap are
available from Refs. Choi et al., 2017; Island et al., 2016; Zelewski and
Kudrawiec, 2017; Conan et al., 1979, while band splittings are provided in
Refs. Oliva et al., 2020; Zhang et al., 2014; Shim et al., 2014; Ross et al.,
2013; Beal et al., 1972; Ruppert et al., 2014. Including spin currents in the
xc form has a very small effect on the band gaps (see Table 2). On the other
hand, the values of valence band splittings are considerably enhanced (from
0.14 eV to 0.20 eV in MoSe2, and from 0.28 eV to 0.32 eV in MoTe2), leading to
an improved agreement with experiment (relative error decreased from 29% to
10% in MoSe2 and leading to an exact coincidence with the mean experimental
value in MoTe2). Overall, the accuracy of the J-r2SCAN approach is comparable
to more involved hybrid-GGA SCDFT treatments (which provide splittings of 0.17
and 0.32 eV, respectively for MoSe2 and MoTe2 using optimal fractions of Fock
exchange) — without requiring the determination of an optimal fraction of Fock
exchange:Desmarais et al. (2023a) a step that involves either empiricism or
extra/external ab initio calculations.Skone et al. (2014)

Table 2: Splitting of the valence band at $\mathbf{K}$ and band gap of molybdenum dichalcogenide Rashba systems. | 2D MoSe2 |  | 3D $\alpha$-MoTe2
---|---|---|---
| split | gap |  | split | gap
r2SCAN | 0.14 | 1.57 |  | 0.28 | 0.82
J-r2SCAN | 0.20 | 1.53 |  | 0.32 | 0.81
Exp. | 0.18 | 1.6-2.3 |  | 0.3-0.34 | 1.03

Finally, we report on the application of our SU(2) form-invariant approach to
topological materials. We consider the case of Weyl semi-metals. We report on
the splitting of Weyl nodes in the orthorhombic TaAs phase. We study the W1
node pair, which lies along the $k_{z}=0$ mirror plane. The pair of nodes is
located at $k_{y}=0.5066$ (in units of $\pi/\mathbf{b}$), in excellent
agreement with the experimental value of $k_{y}=0.5173$.Lv et al. (2015);
Huang et al. (2015) Spin-current dependent terms in the xc functional account
for $21\%$ of the total splitting along $k_{x}$. See Fig. S1 of the
supplementary material.

Figure 1: (a) Effect of SOC on the ELF and JELF (isosurface at 0.005) in the
I2 molecule. $\Delta{\rm SOC}$ (J)ELF is the difference between calculations
with and without SOC on the (J)ELF. The rightmost panel reports the difference
JELF-ELF (isosurface at 0.002). Positive (negative) isosurfaces are in yellow
(blue), while the atoms are represented as purple balls. (b) (top) Valence
bands of single-layer MoSe2 and (bottom) bulk $\alpha$-MoTe2 Rashba materials
without SOC (black) and with SOC, employing the r2SCAN and J-r2SCAN xc models
(dashed and continuous blue lines, respectively).

Conclusions. Focusing on the prominent example of the SCAN functional, we have
shown how to include a dependence on spin currents in a meta-GGA xc functional
form for evaluating SOC effects in such a way as to fulfill the prerequisite
of SU(2) gauge invariance of the xc energy. For the states considered in this
work, the current-dependent form of the approximation is furthermore strongly
suggested by previous studies on the structure of the exact exchange hole (see
Ref. Pittalis et al., 2017), which incidentally, also motivate the inclusion
of the kinetic energy density as a variable on which the functional must
depend. A large body of works have highlighted the usefulness of the kinetic
energy in combination with other quantities as in the ELF (and, thus, as in
the iso-orbital indicator $\alpha$ of SCAN) as a practical means for capturing
relevant local features of the many-electron state in correlated regimes
beyond exchange-only effects. Of course, the deeper question of including the
“right physics” beyond JSCAN remains open. Aside from exact constraints
provided by the SU(2) gauge principle, the success of the benchmark
calculations reported in our paper leads us to believe that we are proceeding
on the right path.

Looking ahead, the approach proposed in this work will help disentangle the
role of the SU(2) invariance in density functional approximations from
questions more directly related to the physics of the electron-electron
interaction. It is particularly appealing to explore the consequence of a full
(compound) U(1)$\times$SU(2) gauge-invariance on general states that may also
break time-reversal symmetry. This may help solve standing difficulties of
present DFT methodologies for magnetized systems.Isaacs and Wolverton (2018);
Jana et al. (2018); Ekholm et al. (2018); Fu and Singh (2018, 2019); Mejía-
Rodríguez and Trickey (2019); Romero and Verstraete (2018); Buchelnikov et al.
(2019); Tran et al. (2020)

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

## References

  * Greiner (2012) W. Greiner, _Relativistic Quantum Mechanics: Wave Equations_ (Springer Berlin Heidelberg, 2012), ISBN 9783642880827\.
  * Yu and Cardona (2010) P. Yu and M. Cardona, _Fundamentals of Semiconductors: Physics and Materials Properties_ , Graduate Texts in Physics (Springer Berlin Heidelberg, 2010), ISBN 9783642007101.
  * Bernevig and Hughes (2013) B. Bernevig and T. Hughes, _Topological Insulators and Topological Superconductors_ (Princeton University Press, 2013), ISBN 9780691151755.
  * Vanderbilt (2018) D. Vanderbilt, _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018), ISBN 9781107157651\.
  * $\bm{Z}$utić et al. (2004) I. $\bm{Z}$utić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).
  * Tokatly (2008) I. V. Tokatly, Phys. Rev. Lett. 101, 106601 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.101.106601.
  * Droghetti et al. (2022) A. Droghetti, I. Rungger, A. Rubio, and I. V. Tokatly, Phys. Rev. B 105, 024409 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.024409.
  * Kohn (1999) W. Kohn, Rev. Mod. Phys. 71, 1253 (1999), URL https://link.aps.org/doi/10.1103/RevModPhys.71.1253.
  * Seidl et al. (1996) A. Seidl, A. Görling, P. Vogl, J. A. Majewski, and M. Levy, Phys. Rev. B 53, 3764 (1996).
  * Perdew et al. (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, et al., PNAS 114, 2801 (2017).
  * Nazarov and Vignale (2011) V. U. Nazarov and G. Vignale, Phys. Rev. Lett. 107, 216402 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.107.216402.
  * Prokopiou et al. (2022) G. Prokopiou, M. Hartstein, N. Govind, and L. Kronik, J. Chem. Theory Comput. 18, 2331 (2022).
  * Vignale and Rasolt (1988) G. Vignale and M. Rasolt, Phys. Rev. B 37, 10685 (1988).
  * Bencheikh (2003) K. Bencheikh, J. Phys. A 36, 11929 (2003).
  * Rohra and Görling (2006) S. Rohra and A. Görling, Phys. Rev. Lett. 97, 013005 (2006).
  * Heaton-Burgess et al. (2007) T. Heaton-Burgess, P. Ayers, and W. Yang, Phys. Rev. Lett. 98, 036403 (2007).
  * Trushin and Görling (2018) E. Trushin and A. Görling, Phys. Rev. B 98, 205137 (2018).
  * Desmarais et al. (2023a) J. K. Desmarais, G. Ambrogio, G. Vignale, A. Erba, and S. Pittalis, _Generalized kohn-sham approach for the electronic band structure of spin-orbit coupled materials_ (2023a), eprint 2309.11158.
  * Vignale and Rasolt (1987) G. Vignale and M. Rasolt, Phys. Rev. Lett. 59, 2360 (1987), URL http://link.aps.org/doi/10.1103/PhysRevLett.59.2360.
  * Tao and Perdew (2005) J. Tao and J. P. Perdew, Phys. Rev. Lett. 95, 196403 (2005).
  * Pittalis et al. (2017) S. Pittalis, G. Vignale, and F. G. Eich, Phys. Rev. B 96, 035141 (2017).
  * Fröhlich and Studer (1993) J. Fröhlich and U. M. Studer, Rev. Mod. Phys. 65, 733 (1993).
  * Jin et al. (2006) P.-Q. Jin, Y.-Q. Li, and F.-C. Zhang, Journal of Physics A: Mathematical and General 39, 7115 (2006).
  * Sun et al. (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.
  * Perdew et al. (2009) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin, and J. Sun, Phys. Rev. Lett. 103, 026403 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.026403.
  * Sun et al. (2013) J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G. I. Csonka, G. E. Scuseria, and J. P. Perdew, Phys. Rev. Lett. 111, 106401 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.106401.
  * Sun et al. (2016) J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, et al., Nat Chem 8, 831 (2016).
  * Tran et al. (2016) F. Tran, J. Stelzl, and P. Blaha, The Journal of Chemical Physics 144, 204120 (2016), ISSN 0021-9606, eprint https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4948636/15511401/204120_1_online.pdf, URL https://doi.org/10.1063/1.4948636.
  * Peng et al. (2016) H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Phys. Rev. X 6, 041005 (2016), URL https://link.aps.org/doi/10.1103/PhysRevX.6.041005.
  * Brandenburg et al. (2016) J. G. Brandenburg, J. E. Bates, J. Sun, and J. P. Perdew, Phys. Rev. B 94, 115144 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.115144.
  * Isaacs and Wolverton (2018) E. B. Isaacs and C. Wolverton, Phys. Rev. Mater. 2, 063801 (2018), URL https://link.aps.org/doi/10.1103/PhysRevMaterials.2.063801.
  * Tran et al. (2019) F. Tran, L. Kalantari, B. Traoré, X. Rocquefelte, and P. Blaha, Phys. Rev. Mater. 3, 063602 (2019), URL https://link.aps.org/doi/10.1103/PhysRevMaterials.3.063602.
  * Becke and Edgecombe (1990a) A. D. Becke and K. E. Edgecombe, The Journal of Chemical Physics 92, 5397 (1990a), eprint http://dx.doi.org/10.1063/1.458517, URL http://dx.doi.org/10.1063/1.458517.
  * Silvi and Savin (1994) B. Silvi and A. Savin, Nature 371, 683 (1994), URL https://www.nature.com/articles/371683a0.
  * 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).
  * Comaskey et al. (2022) W. P. Comaskey, F. Bodo, A. Erba, J. L. Mendoza-Cortes, and J. K. Desmarais, Phys. Rev. B 106, L201109 (2022).
  * Erba et al. (2023) 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. Theor. Comput. 19, 6891 (2023).
  * Abedinpour et al. (2010) S. H. Abedinpour, G. Vignale, and I. V. Tokatly, Phys. Rev. B 81, 125123 (2010).
  * Zhang and Vignale (2016) S. S.-L. Zhang and G. Vignale, Phys. Rev. B 94, 140411 (2016).
  * not (a) At the level of GGAs, this implies no substantial difference because the approximate functional, however, depends on the particle density (and its gradients) explicitly. At the level of a MGGA, this implies that an extra semi-local potential shows up because of the explicit dependence on the the KS kinetic energy. This potential is usually grouped with the scalar potential Lehtola et al. (2020).
  * Neumann et al. (1996) R. Neumann, R. H. Nobes, and N. C. Handy, Mol. Phys. 87, 1 (1996).
  * not (b) ${\vec{\bf G}}$ is pure gauge vector potential, fully originated by the SU(2) transformation.
  * Becke (1996a) A. D. Becke, The Journal of Chemical Physics 104, 1040 (1996a), ISSN 0021-9606, eprint https://pubs.aip.org/aip/jcp/article-pdf/104/3/1040/10779816/1040_1_online.pdf, URL https://doi.org/10.1063/1.470829.
  * Becke and Edgecombe (1990b) A. D. Becke and K. E. Edgecombe, The Journal of Chemical Physics 92, 5397 (1990b), ISSN 0021-9606, eprint https://pubs.aip.org/aip/jcp/article-pdf/92/9/5397/9731640/5397_1_online.pdf, URL https://doi.org/10.1063/1.458517.
  * Savin et al. (1997) A. Savin, R. Nesper, S. Wengert, and T. F. Fässler, Angewandte Chemie International Edition in English 36, 1808 (1997).
  * Dobson (1993) J. F. Dobson, The Journal of Chemical Physics 98, 8870 (1993).
  * Becke (1996b) A. D. Becke, Canadian Journal of Chemistry 74, 995 (1996b).
  * Becke (2002) A. D. Becke, The Journal of Chemical Physics 117, 6935 (2002).
  * Maximoff et al. (2004) S. N. Maximoff, M. Ernzerhof, and G. E. Scuseria, The Journal of Chemical Physics 120, 2105 (2004).
  * Burnus et al. (2005) T. Burnus, M. A. Marques, and E. K. Gross, Physical Review A 71, 010501 (2005).
  * Pittalis et al. (2007) S. Pittalis, S. Kurth, S. Sharma, and E. K. U. Gross, The Journal of Chemical Physics 127, 124103 (2007).
  * Pittalis et al. (2009) S. Pittalis, E. Räsänen, and E. K. U. Gross, Phys. Rev. A 80, 032515 (2009).
  * Räsänen et al. (2009) E. Räsänen, S. Pittalis, C. R. Proetto, and E. K. U. Gross, Phys. Rev. B 79, 121305 (2009).
  * Oliveira et al. (2010) M. J. T. Oliveira, E. Räsänen, S. Pittalis, and M. A. L. Marques, Journal of Chemical Theory and Computation 6, 3664 (2010).
  * Zhu et al. (2016) W. Zhu, L. Zhang, and S. B. Trickey, The Journal of Chemical Physics 145, 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, Ulf Ekström and A. M. Teale, Molecular Physics 114, 1415 (2016), eprint https://doi.org/10.1080/00268976.2015.1133859, URL https://doi.org/10.1080/00268976.2015.1133859.
  * Holzer et al. (2022) C. Holzer, Y. J. Franzke, and A. Pausch, The Journal of Chemical Physics 157, 204102 (2022), ISSN 0021-9606.
  * not (c) SU(2) invariance of the exact xc-energy functional on time-reversal invariant states, implies $\frac{1}{c}\left[\vec{A}_{xc,\mu}\times\vec{J}_{\mu}\right]^{b}=\frac{1}{2c}\partial_{\mu}\left(A_{xc,\mu}^{b}n\right)$ Vignale and Rasolt (1988); Pittalis et al. (2017). The dependence of the exact xc-energy on the spin current via the covariant vorticity, Abedinpour et al. (2010) for the same class of states further implies $\partial_{\mu}\left(A_{xc,\mu}^{b}n\right)=0$. From the two relations, we readily obtain $\frac{1}{c}\left[\vec{A}_{xc,\mu}\times\vec{J}_{\mu}\right]^{b}$ = 0.
  * Aschebrock and Kümmel (2019) T. Aschebrock and S. Kümmel, Phys. Rev. Res. 1, 033082 (2019), URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033082.
  * (61) Computational and implementation details are available at INSERT URL HERE.
  * Desmarais et al. (2021a) J. K. Desmarais, S. Komorovsky, J.-P. Flament, and A. Erba, J. Chem. Phys. 154, 204110 (2021a).
  * Desmarais et al. (2020b) J. K. Desmarais, J.-P. Flament, and A. Erba, Phys. Rev. B 101, 235142 (2020b).
  * 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).
  * Metz et al. (2000) B. Metz, H. Stoll, and M. Dolg, J. Chem. Phys. 113, 2563 (2000).
  * Peterson et al. (2003) K. A. Peterson, D. Figgen, E. Goll, H. Stoll, and M. Dolg, J. Chem. Phys. 119, 11113 (2003).
  * Peterson et al. (2007) K. A. Peterson, D. Figgen, M. Dolg, and H. Stoll, J. Chem. Phys. 126 (2007).
  * Doll et al. (2001) K. Doll, V. Saunders, and N. Harrison, Int. J. Q. Chem. 82, 1 (2001).
  * Doll (2001) K. Doll, Comput. Phys. Commun. 137, 74 (2001).
  * Doll et al. (2006) K. Doll, R. Dovesi, and R. Orlando, Theor. Chem. Acc. 115, 354 (2006).
  * Civalleri et al. (2001) B. Civalleri, P. D’Arco, R. Orlando, V. Saunders, and R. Dovesi, Chem. Phys. Lett. 348, 131 (2001).
  * Laun et al. (2018) J. Laun, D. Vilela Oliveira, and T. Bredow, J. Comp. Chem. 39, 1285 (2018).
  * Heyd et al. (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 et al. (1996) M. D. Towler, A. Zupan, and M. Causà, Comput. Phys. Commun. 98, 181 (1996).
  * Lehtola et al. (2020) S. Lehtola, F. Blockhuys, and C. Van Alsenoy, Molecules 25 (2020).
  * Furness et al. (2022) J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Chem. Phys. 156 (2022).
  * (80) RPAx calculations were performed with NWCHEM version 7.0, see https://nwchemgit.github.io/Excited-State-Calculations.html.
  * Aprà et al. (2020) E. Aprà, E. J. Bylaska, W. A. de Jong, N. Govind, K. Kowalski, T. P. Straatsma, M. Valiev, H. J. J. van Dam, Y. Alexeev, J. Anchell, et al., J. Chem. Phys. 152, 184102 (2020).
  * Desmarais et al. (2021b) J. K. Desmarais, A. Erba, J.-P. Flament, and B. Kirtman, J. Chem. Theor. Comput. 17, 4697 (2021b).
  * Desmarais et al. (2021c) J. K. Desmarais, A. Erba, J.-P. Flament, and B. Kirtman, J. Chem. Theor. Comput. 17, 4712 (2021c).
  * Desmarais et al. (2023b) J. K. Desmarais, A. Boccuni, J.-P. Flament, B. Kirtman, and A. Erba, J. Chem. Theor. Comput. 19, 1853 (2023b).
  * Ullrich (2011) C. A. Ullrich, _Time-dependent density-functional theory: concepts and applications_ (OUP Oxford, 2011).
  * Desmarais et al. (2023c) J. K. Desmarais, A. Erba, and J.-P. Flament, Phys. Rev. B 108, 134108 (2023c).
  * Gulans and Draxl (2022) A. Gulans and C. Draxl, arXiv preprint arXiv:2204.02751 (2022).
  * Choi et al. (2017) B. K. Choi, M. Kim, K.-H. Jung, J. Kim, K.-S. Yu, and Y. J. Chang, Nanoscale Res. Lett. 12, 1 (2017).
  * Island et al. (2016) J. O. Island, A. Kuc, E. H. Diependaal, R. Bratschitsch, H. S. Van Der Zant, T. Heine, and A. Castellanos-Gomez, Nanoscale 8, 2589 (2016).
  * Zelewski and Kudrawiec (2017) S. J. Zelewski and R. Kudrawiec, Sci. Reports 7, 15365 (2017).
  * Conan et al. (1979) A. Conan, D. Delaunay, A. Bonnet, A. Moustafa, and M. Spiesser, physica status solidi (b) 94, 279 (1979).
  * Oliva et al. (2020) R. Oliva, T. Woźniak, F. Dybala, J. Kopaczek, P. Scharoch, and R. Kudrawiec, Mat. Res. Lett. 8, 75 (2020).
  * Zhang et al. (2014) Y. Zhang, T.-R. Chang, B. Zhou, Y.-T. Cui, H. Yan, Z. Liu, F. Schmitt, J. Lee, R. Moore, Y. Chen, et al., Nature nano. 9, 111 (2014).
  * Shim et al. (2014) G. W. Shim, K. Yoo, S.-B. Seo, J. Shin, D. Y. Jung, I.-S. Kang, C. W. Ahn, B. J. Cho, and S.-Y. Choi, ACS nano 8, 6655 (2014).
  * Ross et al. (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, et al., Nature Commun. 4, 1474 (2013).
  * Beal et al. (1972) A. Beal, J. Knights, and W. Liang, J. Phys. C Sol. State Phys. 5, 3540 (1972).
  * Ruppert et al. (2014) C. Ruppert, B. Aslan, and T. F. Heinz, Nano Lett. 14, 6231 (2014).
  * Skone et al. (2014) J. H. Skone, M. Govoni, and G. Galli, Phys. Rev. B 89, 195112 (2014).
  * Lv et al. (2015) B. Lv, H. Weng, B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. Huang, L. Zhao, G. Chen, et al., Phys. Rev. X 5, 031013 (2015).
  * Huang et al. (2015) S.-M. Huang, S.-Y. Xu, I. Belopolski, C.-C. Lee, G. Chang, B. Wang, N. Alidoust, G. Bian, M. Neupane, C. Zhang, et al., Nat. Comms. 6, 7373 (2015).
  * Jana et al. (2018) S. Jana, A. Patra, and P. Samal, The Journal of Chemical Physics 149, 044120 (2018), ISSN 0021-9606, eprint https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5040786/14927796/044120_1_online.pdf, URL https://doi.org/10.1063/1.5040786.
  * Ekholm et al. (2018) M. Ekholm, D. Gambino, H. J. M. Jönsson, F. Tasnádi, B. Alling, and I. A. Abrikosov, Phys. Rev. B 98, 094413 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.094413.
  * Fu and Singh (2018) Y. Fu and D. J. Singh, Phys. Rev. Lett. 121, 207201 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.207201.
  * Fu and Singh (2019) Y. Fu and D. J. Singh, Phys. Rev. B 100, 045126 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.045126.
  * Mejía-Rodríguez and Trickey (2019) D. Mejía-Rodríguez and S. B. Trickey, Phys. Rev. B 100, 041113 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.041113.
  * Romero and Verstraete (2018) A. Romero and M. Verstraete, Eur. Phys. J. B 91, 193 (2018), URL https://link.springer.com/article/10.1140/epjb/e2018-90275-5#citeas.
  * Buchelnikov et al. (2019) V. D. Buchelnikov, V. V. Sokolovskiy, O. N. Miroshkina, M. A. Zagrebin, J. Nokelainen, A. Pulkkinen, B. Barbiellini, and E. Lähderanta, Phys. Rev. B 99, 014426 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.014426.
  * Tran et al. (2020) F. Tran, G. Baudesson, J. Carrete, G. K. H. Madsen, P. Blaha, K. Schwarz, and D. J. Singh, Phys. Rev. B 102, 024407 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.024407.


# Learning Dense Correspondence from Synthetic Environments

###### Abstract

Estimation of human shape and pose from a single image is a challenging task.
It is an even more difficult problem to map the identified human shape onto a
3D human model. Existing methods map manually labelled human pixels in real 2D
images onto the 3D surface, which is prone to human error, and the sparsity of
available annotated data often leads to sub-optimal results. We propose to
solve the problem of data scarcity by training 2D-3D human mapping algorithms
using automatically generated synthetic data for which exact and dense 2D-3D
correspondence is known. Such a learning strategy using synthetic environments
has a high generalisation potential towards real-world data. Using different
camera parameter variations, background and lighting settings, we created
precise ground truth data that constitutes a wider distribution. We evaluate
the performance of models trained on synthetic using the COCO dataset and
validation framework. Results show that training 2D-3D mapping network models
on synthetic data is a viable alternative to using real data.

Index Terms—  Synthetic Data Generation, SMPL, Densepose, Semantic
Segmentation, Human modelling.

##  1 Introduction

Estimating human shape and pose from a single 2D image and establishing
correspondence with a 3D human surface is a hard problem. Challenges from the
lack of 3D information in 2D images range from depth ambiguity to occlusion or
imaging conditions (e.g lighting) as well as the difficulty and tediousness of
creating (sparse) ground truth data and overall natural variability in human
pose and shape. It is a primary task for many applications including virtual
reality, human action recognition, augmented reality and human-robot
interaction [1, 2]. 2D-3D human mapping involves other difficult problems
including object detection, pose estimation and semantic segmentation as
auxiliary tasks. Patel et al. [3] introduced a synthetic human dataset for
regression analysis containing 3D human scans placed in scenes that performed
well on fine tuning the SMPL (Skinned Multi-Person Linear) oPtimization IN the
loop (SPIN)[4] model that fits a human shape model to humans in 2D images. [5]
introduced another synthetic dataset that improves pedestrian detection and
tracking. Thus, the adaptability and reliability of synthetic data over real
data are increasing over time.

There is a substantial body of work trying to solve the problem of mapping a
2D image to a 3D model with dense point correspondence for inanimate objects
[6, 7, 8]. For humans, 3D representations, parametric surface models such as
SMPL [9] or Adam model [10] are often used to create controllable 3D surface
deformations. Pixel aligned implicit functions have been used for high
resolution clothed human digitization by transferring the human geometry onto
a 3D model based on a single image. The method is non-parametric, which
decreases the bias towards training data [11], however, when complex motions
are involved, parametric models like SMPL can be constrained more easily and
provide more temporally consistent and realistic animations.

DenseReg [12] uses Convolutional Neural Networks (CNNs) to establish dense
correspondence between 2D human face images and 3D models. Challenges involved
a high structural complexity and pose variability. A 24 part IUV (Texture
coordinate space) regression to achieve dense correspondence between image
pixels and surface points was introduced in [13]. A cascaded extension of
Densepose-RCNN and an in-painting network were used to interpolate a dense
supervision from sparse supervision signals. Guler et al. [13] published a
method to regress human surface coordinates using a Mask-RCNN architecture.
They also released the COCO (Common Objects In Context) Denspose dataset that
consists of human 2D images annotated with human semantic part segmentation
and UV coordinate mapping within each semantic part [13]. These UV coordinate
mappings establish human correspondence with the 3D human model [9]. To create
the annotations, humans were hired who spent a lot of time to create
annotations that where rather sparse and prone to errors. This is a big
limitation for a task like dense correspondence where the reliability on human
annotators results in the creation of erroneous data. Our work focuses on
resolving this issue with the dataset.

Our contributions are two fold. First, we present a method to automatically
generate synthetic data and extract ground truth annotations for the problem
of dense 2D-3D human correspondence. The data generated consists of human
semantic segmentation, exact 2D and 3D human pose keypoints and ground truth
2D pixels to 3D human voxel mappings. Second, we evaluate the performance of
our labelled synthetic data for the problem of dense correspondence mapping
between 2D human images and 3D human SMPL models [9] using the method in [13].

##  2 Methodology

In this section, we describe the completely automated method used to generate
accurate ground truth UV data of a human mesh using a synthetic 3D environment
where human avatars are made to perform human actions using motion capture
(MOCAP) data.

###  2.1 Avatar animation and rendering

Part of the MOCAP data was collected at the Queensland Centre for Advanced
Technologies. Ten volunteers performed a set of actions (e.g., walking,
running, seating) while wearing reflective markers. The 3D position of these
markers was captured using 24 synchronised infrared cameras and Qualisys Track
Manager was used to reconstruct and animate 3D human skeletons. Additional
publicly available MOCAP data were included to increase diversity [14].

The SMPL parametric model was used to create the human avatar [9]. It is a
compact representation of the human body as an average human surface paired
with modes of variation (statistical deformations driving vertices). It uses
user-defined shape known as blend-shapes to create the general shape of the
human avatar, and pose blend-shapes that are automatically updated by the
model based on the human pose to account for tissue deformation. The model can
thus be used to generate an infinite number of statistically plausible human
body shapes, which animate realistically based on MOCAP data [9].

Publicly available textures from the Surreal dataset[15] were mapped onto the
avatar to generate a realistic human look. The male and female textures were
separated and the semantic parts (torso, arms etc.) mixed and stitched
together to increase variety. To render realistic human images, we used 3456
random background images from the COCO Densepose train set and the synthetic
human model is made to move in 3D space between the cone of the camera and
image plane based on animation.

The Unity Real-Time Development Platform was used to render all video data. 9
different synthetic cameras were used to create synthetic data where avatar
animation is overlayed onto the existing background image that is placed
within the cone of the rendering camera. Synthetic camera noise and camera
lens distortion is added to the rendered scene to impart more realism.
Separate SMPL human models are used for male and female bodies to accommodate
precise physiological representation. We used a single windows desktop with
GPU to render the images.

Fig. 1: Synthetic Data Generation Pipeline

###  2.2 Mesh Vertex visibility test in 3D

The vertex visibility in 3D is dependent on the camera view $\alpha$ and pose
parameters $\theta$. We compute the concave hull of the SMPL human mesh in 3D
from the mesh vertices to get precise visible vertex locations of the mesh
surface. We used the method of Ray casting from each vertex locations to the
corresponding camera center to solve the vertex visibility for a given view.
All visible vertices in the world coordinate 3D space are transformed into
image space in-order to retrieve visible vertices of the human 3D mesh in the
rendered 2D image space.

$Visibility(v_{i})=F(\alpha,\theta),$ (1)

where $v_{i}$ denotes 3D vertex coordinate and $\alpha$, $\theta$ are camera
parameters, pose parameters respectively.

###  2.3 Semantic region mapping and IUV segmentation

We divide the human mesh in SMPL UV space into 24 part atlas UV space as
discussed in [13]. The whole human mesh is divided into 14 semantic regions
like Head, Torso, Lower/Upper Arms, Lower/Upper Legs, Hands and Feet etc. To
extract ground truth segmentations, all image space coordinates are projected
on to the world coordinate space to see if these points collide with the
concave hull of the human mesh. All such coordinates of collision points in 2D
are obtained as a dense mapping of the human mesh in the image space. To
extract semantic region mapping, we project all image space coordinates from
the corresponding camera to check whether the rays hit the concave hull field.
If a point in image space when casted, hits the field, then the collision
point is associated with the geometrically closest vertex coordinate in world
coordinate space such that a correlation is established between all 2D points
in image space on the human surface and mesh vertices. The association created
between image space coordinates and UV space coordinates helps to gather
precise semantic region information in the image domain as shown in Figure 2.

Fig. 2: Human Dense Correspondence Mapping on SMPL Human model

###  2.4 Occlusion generation and Alpha blending

To create occlusions, we used 227 segmented 2D mundane objects which were
extracted using a background removal tool. An object is randomly sampled and
re-scaled according to the human Region of Interest (ROI) in the scene. Then,
the object is positioned randomly over the the rendered human in the image.
The 2D objects that are used to create occlusion have very sharp edges when
randomly positioned on top of the human surface. To provide smooth transition
between the edge of the 2D object and the rendered image , a random Gaussian
filter and alpha blending are applied over a thick inner and outer boundary of
the object. This allows for seamless transition between object edge and
background.

###  2.5 Scene harmonisation for visual data generation

To better match the “style” of the human avatar and occluding object with the
background, a scene harmonisation method was used [16]. It uses a generative
adversarial network conditioned on the ROI (the human mask and the occluded
object mask) that transforms the target human pixel features to adapt more
with the background style of the image as depicted in Figure 3.

Fig. 3: Synthetic Human Dataset Post Processing

##  3 EXPERIMENTS & RESULTS

###  3.1 COCO Densepose Training

Our goal is to evaluate the performance of synthetic data generated for the
problem of dense correspondence estimation. We used the detectron2 for
training the dense correspondence model. For all experiments, we used
Resnet-50 backbone as the primary feature extractor and the Panoptic head from
[17] and DeepLabV3 head from [18]. A hyperparameter free training strategy is
followed according to the linear scaling rule from [19]. Comparison between
the same model setting trained using real data and synthetic data are
discussed in the coming sub-sections. We discovered that training the model
from scratch using our synthetic data caused the model to over-fit the
synthetic data in fewer iterations, thus losing generalisation capability on
real data. To avoid this, we initialised the model with weights trained on
real data and only the last layers i.e, the part of the model where the IUV
estimation occurs, is fine-tuned. We also used an additional mask-based
supervision for training the coarse segmentation parts (14 parts) whose
weights are randomly initialised. Thus, the human mask estimation and the IUV
dense correspondence estimation is completely learned from our synthetic data.
We trained the model for 35,000 iterations with batch size 2, which is almost
equivalent to an epoch at a learning rate 0.001 on a single Nvidia Tesla T4
GPU to avoid over-fitting. To train the model, we used 73,671 synthetic human
instances containing directly rendered images without occlusion and images
post-processed by adding occlusion along with scene harmonisation.

Task | AP | $\textbf{AP}_{50}$ | $\textbf{AP}_{75}$ | AR | $\textbf{AR}_{50}$ | $\textbf{AR}_{75}$
---|---|---|---|---|---|---
bbox | 58.35 | 84.89 | 63.99 |  |  |
GPS | 52.18 | 85.69 | 57.96 | 61.07 | 91.48 | 69.19
GPSm | 54.38 | 88.02 | 63.79 | 62.08 | 93.40 | 72.89
Segm | 57.30 | 91.69 | 66.62 | 64.91 | 94.07 | 76.19

Table 1: Model Fine-tuned with Synthetic data

Task | AP | $\textbf{AP}_{50}$ | $\textbf{AP}_{75}$ | AR | $\textbf{AR}_{50}$ | $\textbf{AR}_{75}$
---|---|---|---|---|---|---
bbox | 61.08 | 88.04 | 67.10 |  |  |
GPS | 65.54 | 91.73 | 75.36 | 73.02 | 95.09 | 81.32
GPSm | 66.69 | 92.74 | 79.62 | 72.61 | 96.03 | 84.31
Segm | 68.33 | 94.54 | 81.51 | 73.74 | 96.70 | 86.22

Table 2: Model trained on real data (Pre-trained model available in detectron2
repo)

Task | AP | $\textbf{AP}_{50}$ | $\textbf{AP}_{75}$ | AR | $\textbf{AR}_{50}$ | $\textbf{AR}_{75}$
---|---|---|---|---|---|---
bbox | 61.08 | 88.04 | 67.10 |  |  |
GPS | 53.87 | 88.18 | 59.54 | 61.25 | 91.84 | 68.83
GPSm | 56.56 | 90.52 | 66.90 | 62.44 | 93.89 | 73.25
Segm | 59.07 | 92.44 | 69.67 | 65.38 | 94.78 | 76.59

Table 3: Model Fine-tuned with Synthetic data inferred with real bbox Fig. 4:
Semantic part segmentation comparison  Fig. 5: Texture extraction comparison

###  3.2 Evaluation Metrics

We report the standard COCO densepose evaluation metrics described in [13]. We
consider the human bounding box estimation (bbox), human segmentation mask
estimation (Segm), Geodesic Point Similarity(GPS) and Geodesic Point
Similarity within the correctly predicted human part mask (GPSm) to evaluate
the dense correspondence. For all tasks, mIoU (mean Intersection over Union)
is considered with thresholds and AP (average precision averaged over IoU
thresholds) is the primary metric followed by $\textbf{AP}_{50}$ and
$\textbf{AP}_{75}$ being the selected supplementary metrics.

###  3.3 Quantitative Evaluation

The COCO densepose test set containing 1508 images with 5581 person instances
are used for performing quantitative evaluation. Table 1 shows the results of
the dense correspondence model fine-tuned with synthetic data. Table 2 shows
the state-of-the-art results[13]. Table 3 shows another experiment where the
bounding box estimation is performed by the state-of-the-art model while all
the remaining metrics are evaluated by the model fine-tuned with synthetic
data. The results depict that with better known human ROIs, all other metrics
showcase improvements.

###  3.4 Qualitative Evaluation

We provide some comparative results between the human semantic region
estimation by the model fine-tuned using synthetic data and estimation
performed by the pre-trained model trained on real data in Figure 4. The human
pixels inferred within the ROI are then transformed to the SMPL UV space for
further comparison of the 2D human image to 3D surface mapping. Figure 5 shows
comparison of the human textures inferred by the model fine-tuned with
synthetic data and the pre-trained model trained using real data.

##  4 CONCLUSION

We studied the possibility of training deep-learning models for human 2D-3D
mapping using synthetic data, which has the advantage of being dense, precise,
unbiased, labour free and solves issues of data volume and privacy. Initial
experiments show promising quantitative results and qualitative results that
the synthetic model has the potential to improve fine edge accuracy. The
current gap between synthetic and real data can cause an over-fitting of the
model to the synthetic data, which can be solved using advanced domain
adaptation techniques. This presents a clear avenue for work on synthetic data
for 2D-3D mapping.

## References

  * [1] Gül Varol, Ivan Laptev, Cordelia Schmid, and Andrew Zisserman,  “Synthetic humans for action recognition from unseen viewpoints,”  International Journal of Computer Vision, vol. 129, no. 7, pp. 2264–2287, 2021.
  * [2] Denis Zherdev, Larisa Zherdeva, Sergey Agapov, Anton Sapozhnikov, Artem Nikonorov, and Sergej Chaplygin,  “Producing synthetic dataset for human fall detection in ar/vr environments,”  Applied Sciences, vol. 11, no. 24, pp. 11938, 2021.
  * [3] Priyanka Patel, Chun-Hao P Huang, Joachim Tesch, David T Hoffmann, Shashank Tripathi, and Michael J Black,  “Agora: Avatars in geography optimized for regression analysis,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 13468–13478.
  * [4] Nikos Kolotouros, Georgios Pavlakos, Michael J Black, and Kostas Daniilidis,  “Learning to reconstruct 3d human pose and shape via model-fitting in the loop,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 2252–2261.
  * [5] Matteo Fabbri, Guillem Brasó, Gianluca Maugeri, Orcun Cetintas, Riccardo Gasparini, Aljoša Ošep, Simone Calderara, Laura Leal-Taixé, and Rita Cucchiara,  “Motsynth: How can synthetic data help pedestrian detection and tracking?,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 10849–10859.
  * [6] Carl Martin Grewe and Stefan Zachow,  “Fully automated and highly accurate dense correspondence for facial surfaces,”  in European conference on computer vision. Springer, 2016, pp. 552–568.
  * [7] Jonathan Taylor, Jamie Shotton, Toby Sharp, and Andrew Fitzgibbon,  “The vitruvian manifold: Inferring dense correspondences for one-shot human pose estimation,”  in 2012 IEEE Conference on Computer Vision and Pattern Recognition. IEEE, 2012, pp. 103–110.
  * [8] Tinghui Zhou, Philipp Krahenbuhl, Mathieu Aubry, Qixing Huang, and Alexei A Efros,  “Learning dense correspondence via 3d-guided cycle consistency,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 117–126.
  * [9] Matthew Loper, Naureen Mahmood, Javier Romero, Gerard Pons-Moll, and Michael J Black,  “Smpl: A skinned multi-person linear model,”  ACM transactions on graphics (TOG), vol. 34, no. 6, pp. 1–16, 2015\.
  * [10] Hanbyul Joo, Tomas Simon, and Yaser Sheikh,  “Total capture: A 3d deformation model for tracking faces, hands, and bodies,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 8320–8329.
  * [11] Shunsuke Saito, Zeng Huang, Ryota Natsume, Shigeo Morishima, Angjoo Kanazawa, and Hao Li,  “Pifu: Pixel-aligned implicit function for high-resolution clothed human digitization,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 2304–2314.
  * [12] Riza Alp Guler, George Trigeorgis, Epameinondas Antonakos, Patrick Snape, Stefanos Zafeiriou, and Iasonas Kokkinos,  “Densereg: Fully convolutional dense shape regression in-the-wild,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 6799–6808.
  * [13] Rıza Alp Güler, Natalia Neverova, and Iasonas Kokkinos,  “Densepose: Dense human pose estimation in the wild,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 7297–7306.
  * [14] Yang Yang, Chunping Hou, Yue Lang, Dai Guan, Danyang Huang, and Jinchen Xu,  “Open-set human activity recognition based on micro-doppler signatures,”  Pattern Recognition, vol. 85, pp. 60–69, 2019.
  * [15] Linxi Fan, Yuke Zhu, Jiren Zhu, Zihua Liu, Orien Zeng, Anchit Gupta, Joan Creus-Costa, Silvio Savarese, and Li Fei-Fei,  “Surreal: Open-source reinforcement learning framework and robot manipulation benchmark,”  in Conference on Robot Learning, 2018.
  * [16] Jun Ling, Han Xue, Li Song, Rong Xie, and Xiao Gu,  “Region-aware adaptive instance normalization for image harmonization,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 9361–9370.
  * [17] Alexander Kirillov, Ross Girshick, Kaiming He, and Piotr Dollár,  “Panoptic feature pyramid networks,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019, pp. 6399–6408.
  * [18] Liang-Chieh Chen, George Papandreou, Florian Schroff, and Hartwig Adam,  “Rethinking atrous convolution for semantic image segmentation,”  arXiv preprint arXiv:1706.05587, 2017.
  * [19] Priya Goyal, Piotr Dollár, Ross Girshick, Pieter Noordhuis, Lukasz Wesolowski, Aapo Kyrola, Andrew Tulloch, Yangqing Jia, and Kaiming He,  “Accurate, large minibatch sgd: Training imagenet in 1 hour,”  arXiv preprint arXiv:1706.02677, 2017.


# Kittel law and domain formation mechanism in PbTiO3/SrTiO3 superlattices.

Fernando Gómez-Ortiz  Departamento de Ciencias de la Tierra y Física de la
Materia Condensada, Universidad de Cantabria, Avenida de los Castros s/n 39005
Santander. Spain.  Hugo Aramberri  Materials Research and Technology
Department, Luxembourg Institute of Science and Technology (LIST), Avenue des
Hauts-Fourneaux 5, L-4362 Esch/Alzette, Luxembourg  Juan M. López  Instituto
de Física de Cantabria (IFCA), Universidad de Cantabria-CSIC, 39005 Santander,
Spain  Pablo García-Fernández  Departamento de Ciencias de la Tierra y Física
de la Materia Condensada, Universidad de Cantabria, Avenida de los Castros s/n
39005 Santander. Spain.  Jorge Íñiguez  Materials Research and Technology
Department, Luxembourg Institute of Science and Technology (LIST), Avenue des
Hauts-Fourneaux 5, L-4362 Esch/Alzette, Luxembourg  Department of Physics and
Materials Science, University of Luxembourg, Rue du Brill 41, L-4422 Belvaux,
Luxembourg  Javier Junquera<EMAIL_ADDRESS>Departamento de
Ciencias de la Tierra y Física de la Materia Condensada, Universidad de
Cantabria, Avenida de los Castros s/n 39005 Santander. Spain.

###### Abstract

We report second-principles simulations on the structural and energetic
properties of domains in (PbTiO3)n/(SrTiO3)n superlattices. For the explored
layer thickness ($n$ ranging between 8 and 16 unit cells) and lateral sizes of
the domains, the most stable configuration corresponds to polar domains
separated by a sequence of counter-rotating vortices
(clockwise/counterclockwise) perpendicular to the stacking direction and
acting as domain walls. The balance between the domain wall energy and the
electrostatic energy yields to an optimal domain period $\omega$ that is
proportional to the square-root of the thickness of the PbTiO3 layer,
following the Kittel law. For a given lateral size of the simulation box,
suboptimal domain structures (with a width larger than the one predicted by
the Kittel law) can be obtained in a metastable form. However, at finite
temperature, molecular dynamics simulations show the spontaneous change of
periodicity, which implies the formation of new domains whose generation is
initiated by the nucleation of vortices and antivortices at the interface
between the SrTiO3 and the PbTiO3 layers. The vortices progressively elongate
and eventually annihilate with the antivortices yielding the formation of new
domains to comply the Kittel law via a topological phase transition.

††preprint: APS/123-QED

##  I Introduction

A common feature among the family of ferroelectric materials is the formation
of domain structures, regions of space with different polarization separated
by boundaries called domain walls [1, 2]. Domains of opposite polarization
lead to an overall charge neutrality at the surfaces reducing the
depolarization field and the associated electrostatic energy.

The structure and energetics of domains in ferroic materials were first
addressed by Landau and Lifshitz [3], and one decade later by Kittel [4, 5] in
his studies on ferromagnetic domains. The delicate balance between the energy
of the boundary between domains, the magnetic field energy of the
configuration, and the anisotropy energy of the spin orientation determine the
relationship between the width of the domains, $\omega$, and the thickness of
the material, $d$ [4]. Adding up all the energy costs and minimizing this with
respect to the domain size leads to a square-root dependence of $\omega$ as a
function of $d$. This is the so-called Landau-Kittel law,
$\frac{w^{2}}{\delta}=A\cdot d$, where $A$ is an adimensional proportionality
constant and $\delta$ is the thickness of the domain wall. This law was
extended to ferroelectric materials by Mitsui and Furuichi [6] studying the
domain structure of the Rochelle salt, where the electrostatic, elastic, and
gradient energies determine the number and width of the domains for a given
thickness of the material. The square-root dependence was further generalised
under specific periodicities and screening conditions to the case of ultrathin
ferroelectric layers [7] and to the case of superlattices with paraelectric
materials [8]. Moreover, Roitburd expanded it also for ferroelastic thin films
under epitaxial strain [9]. Therefore, it seems that the Landau-Kittel law of
stripe domain width on film thickness is a general property of all ferroics
[2].

Beyond the analytical derivations, the validity of the law has been confirmed
by first-principles-based studies in ferroelectric [10] and multiferroic [11]
thin films with thicknesses down to three unit cells. Following the spirit of
these two works, we widen the applicability of the Kittel law to the case of
ferroelectric/dielectric superlattices characterized by a ground state
consisting of polar domains separated by a sequence of counter-rotating
vortices (clockwise/counterclockwise) acting as domain walls [12, 13, 14].
Using second-principles simulations we validate the law for this complex
polarization texture. Very interestingly we show how when the system is
initialized from a metastable state, where the density of domains is smaller
than the one predicted by the Kittel law, it evolves upon heating to the
ground state via a topological phase transition. The driving mechanism for the
generation and closure of new domains is the recombination of vortex and
antivortex defects generated at the interface between PbTiO3 and SrTiO3.

##  II Methodology

The second-principles simulations were performed using the same methodology
presented in previous works [15, 16], as implemented in the scale-up package
[15, 17]. The second-principles parameters of both materials were fitted from
density functional theory imposing a hydrostatic pressure of $-11.2$ GPa to
counter the underestimation obtained by the local density approximation of the
cubic lattice constant that was taken as the reference structure. We imposed
an epitaxial constraint assuming in-plane lattice constants of $a=b=3.901$ Å,
forming an angle $\gamma=90^{\circ}$ mimicking the conditions of a SrTiO3
substrate. The interatomic potentials, and the approach to simulate the
interface, are the ones first introduced in Ref. [16]. For a given value of
supercell periodicity $n$, (PbTiO3)n/(SrTiO3)n several values of lateral size
$L$ were relaxed making them commensurate with the number of simulated
domains. We solved the models by running Monte Carlo-simulated annealing from
$60$ K down to very low temperatures, typically comprising 20,000 relaxation
sweeps. Regular Langevin molecular dynamics methods at constant temperature of
$T=90$ K were also used to solve the models in order to follow the dynamics of
the emergent domains. For computational feasibility, we have focused on a
simulation supercell made of a periodic repetition of $L\times 1\times 2n$
elemental perovskite unit cells for the Monte Carlo-simulated annealings and
of $L\times 10\times 2n$ for following the dynamics of the domains. As proven
in [18], at low temperatures ($T<73$ K), the vortices do not vary along the
axial $y$-direction. Therefore, the simplification of taking one unit cell
along this direction does not affect the validity of the model while it speeds
up the calculations. However, when $T=90$ K a sufficiently high number of unit
cells along the $y$-direction must be considered in order to account for the
variation along the axial direction.

We obtained the force-constant band calculations using the direct supercell
approach as implemented in the phonopy package [19]. To this end, we
considered the high-symmetry unit cell of the superlattices (in which the
atoms in the PbTiO3 and SrTiO3 layers occupy the cubic-like perovskite
positions), and repeated it 4$\times$4 times in the $xy$-plane to build the
supercell for the calculations, which we found to be large enough to yield
well-converged results. We included the non-analytical contribution to the
bands which accounts for the splitting between the longitudinal and transverse
polar bands as implemented in [20].

Local polarizations are computed within a linear approximation of the product
of the Born effective charge tensor times the atomic displacements from the
reference structure positions divided by the volume of the unit cell.

##  III Results

###  III.1 Validation of Kittel law

We have checked the validity of the Kittel law in our superlattices following
a similar recipe as in Refs. [10, 11], For different layer thicknesses $n$
with constant and equal dielectric/ferroelectric ratio, ranging between 8 and
16 unit cells, the lateral size of the supercell to host two domains was
optimized [see Fig. 1(a)]. In order to achieve this goal for every value of
$n$ different lengths of the supercell along the $x$-direction, $L$, were
simulated. Once $n$ and $L$ were fixed, the initial atomic positions were
chosen to mimic a couple of pure Ising domains, where the polarization changes
abruptly from pointing upwards (up-domain) to downwards (down-domain) along
the $z$-direction in just one unit cell, as shown in Fig. 1(b). This
configuration was taken as the starting point of the Monte Carlo annealing.
The resulting typical dipole configuration, a local minimum at low
temperature, is shown in Fig. 1(c). The spontaneous formation of alternating
pairs of clockwise/counterclockwise vortices along the $x$-direction is
clearly visible, together with the development of an axial component of the
polarization along the $y$-direction. These vortices were already
theoretically predicted from phenomenological theories [21, 22], first-
principles-based effective Hamiltonian [23], second-principles simulations
[16, 14, 18] or full first-principles calculations [12], and experimentally
demonstrated [13] in PbTiO3/SrTiO3 superlattices.

Figure 1: Structural relaxation of a two-domain structure in
(PbTiO3)n/(SrTiO3)n superlattices. (a) Schematic view of the simulation cell
periodically repeated in the three directions of the space. Red (blue) regions
indicate the positive (negative) polarization domains along the $z$ direction.
(b) Local dipoles of the initial Ising-like configuration presenting a clear
two-domain structure for a layer thickness of $n$ = 12, and a lateral size
$L=12$ u.c. (c) Corresponding pattern of local dipoles after relaxation,
showing an alternating clockwise/counterclockwise polar vortex configuration
[14] along $x$. The domain wall thickness $\delta$ extends over one unit cell.
Colors represent the axial component of the polarization, perpendicular to the
plane defined by the vortices.

For a given layer thickness, the energy per five-atom unit cell as a function
of the lateral size of the supercell, always assuming the presence of two
domains in the simulation box, is shown in Fig. 2(a). The first observation
that can be drawn is that the larger the layer thickness, $n$, the smaller the
energy per unit cell. This fact stems from two different causes. On the one
hand, the polarization in the PbTiO3 layers increases with $n$, tending to the
bulk polarization value and approaching the ground state of the domain. On the
other hand, the larger the layer thickness the smaller the polarization within
the SrTiO3 layer; the system undergoes a transition from an electrostatically
“coupled” regime for small $n$ to a “decoupled” regime for large $n$ [24, 25,
26]. Since the SrTiO3 layers will be closer to the ground-state unpolarized
configuration, the energy is also reduced. The second observation that can be
drawn is that the energy curve of the two-domain structure as a function of
the lateral size $L$ presents a minimum corresponding to the most stable
geometry. The larger the layer thickness the shallower the minimum, that is
localized for longer values of $L$. From these minima, we can infer the
optimal width of the domain, $\omega$. Assuming that the domain wall is one-
unit-cell thick [in good agreement with our simulations; see Fig. 1(c) for a
typical case], $\omega=L/2-1$ unit cells. Plotting the square of these optimal
widths against the thickness of the PbTiO3 layer, red dots in Fig. 2(b), we
recover a linear behavior as predicted by the Kittel law, a tendency also
shown in BiFeO3 [11] and Pb$($Zr,Ti$)$O3 ultrathin films [10].

Figure 2: Optimization of the lateral size of a two-domain structure in
(PbTiO3)n/(SrTiO3)n superlattices. (a) Energy per five atoms unit cell as a
function of the lateral size of the supercell for different periodicities: red
dots ($n=8$), blue squares ($n=10$), green up-triangles ($n=12$), black
diamonds ($n=14$), and magenta down-triangles ($n=14$). (b) Linear fit of the
squared optimized width of the domain as a function of the layer thickness by
two procedures. Red dots correspond to minimizing the total energy of the
supercell while blue squares correspond to the minimum of the force-constant
bands along $\Gamma-X$.

Alternatively, one can try to predict the optimal width of the multidomain
structure by studying force-constant bands like the ones presented in Fig.
3(a). Following the unstable modes along $\Gamma-X$ we can identify $q_{\rm
min}$ as the wave-vector associated with the strongest instability describing
vortex structures as the ones presented in Fig. 3(b). From this value we can
infer the optimal width as $\omega=\frac{1}{2q}-1$ [see blue squares in Fig.
2(b)]. Interestingly, while the square-root dependence of the domain width as
a function of the layer periodicity is nicely reproduced, there exists a
discrepancy with the results obtained via the first method. The (harmonic)
force-constant analysis predicts narrower domains than the ones obtained from
a full energy minimization. This difference must be due to anharmonic effects
and/or the fact that the fully relaxed structures feature a combination of
phonon mode distortions to optimize the energy. Thus, for example, the
development of a slight offset coupled with an in-plane component of the
polarization as the one shown in Fig. 1 (c) and experimentally attained [13],
although not captured in the force-constant analysis [see how vortices are
centered in Fig. 3(b)], reduces the normal component of the polarization to
the surface. The tilt of the polarization reduces the depolarization charges
at the surface and allows the widening of the domains, resulting in the
underestimate of the domain width predicted by the harmonic analysis.

Figure 3: (a) Force-constant bands along $\Gamma-X$ obtained after
diagonalization of the force-constant matrix in a(PbTiO3)9/(SrTiO3)9
supercell. Dashed line indicates the energy of the centrosymetric
configuration which is taken as the reference. Red dot indicates the position
of the strongest instability located at $q_{\rm min}=0.123$ in fractional
units.(b) Relaxed structure following the strongest instability of the force-
constant bands.

###  III.2 Domain formation

Up to now, the existence of two domains in the simulation box has been imposed
in the calculations. In the following we shall study the phase competition
between different configurations under suboptimal domain widths. This is shown
in Fig. 4 where we compare the energy per unit cell considering two or four-
domains for a given layer thickness, $n$, and lateral size of the supercell,
$L$. Again, for each pair of $(n,L)$ a Monte-Carlo annealing was carried out
starting from purely 180∘ Ising-like structures containing two- or four-domain
walls in the simulation box. In every case, the final relaxed structures
display the typical polar vortex configuration, similar to the one shown in
Fig. 1(c). In Fig. 4 we plot the evolution of the energy profile per five-
atoms unit cell as a function of the lateral size of the supercell. As already
discussed in Fig. 2(a), for a given amount of domain walls in the supercell,
the larger the layer thickness $n$ the smaller the total energy. But the most
important conclusion that can be drawn from Fig. 4 is that, for a given $n$, a
crossover between the two- and four- domain configurations is found for a
critical length $L_{\rm c}$, whose value increases with $n$, as marked by the
filled squares in Fig. 4.

Figure 4: Energy profile per five atom unit cell as a function of the lateral
size of the supercell for the two- (filled dots) and four- (open dots) domain
structures. Different layer thicknesses are indicated by colors: red ($n=8$),
blue ($n=10$), green ($n=12$), black ($n=14$) and magenta ($n=16$). Filled
black squares indicate the crossing point where the four-domain structure
becomes more stable. Numbers 1-4 label the different dipole patterns plotted
in Fig. 5.

Studying lateral sizes far from the coexistence regions (below or above
$L_{\rm c}$) different characteristic features can be observed for the
evolution of the unstable phases. In Fig. 5 we show the behaviour for a layer
thickness of $n=12$. Below $L_{\rm c}$ (short lateral sizes) the four-domain
structure is only metastable since the large penalty coming from the gradient
energy term is not compensated by the saving in electrostatic energy. Indeed,
in this metastable regime the system tends to reduce the energy gradient
contribution by forming polarization waves [27] [Fig. 5(a)] along the
$[100]_{\rm pc}$ direction, with the concomitant onset of a net in-plane
polarization, and the displacement of the center of the vortices in the PbTiO3
layer towards the interface with SrTiO3. The smaller the lateral size of the
supercell $L$, the larger the offset between neighboring vortex cores that do
not fit in the center of the PbTiO3 layer as in the case of larger lateral
sizes [see Fig. 5(d)]. This phenomenon has also been observed in BiFeO3
ultrathin films [11].

The energetically most favourable configuration presents only two domains
[Fig. 5(c)], with larger domain sizes and a smaller number of domain walls.
However, the energy of this configuration increases with $L$. Above $L_{\rm
c}$ (long lateral sizes), the two-domain structure [Fig. 5(d)] becomes
metastable since the electrostatic energy penalty starts to grow and become
dominant. Therefore, we observe new patterns containing two vortices and two
antivortices at the interface with SrTiO3, as will be further discussed in
Fig. 6. Indeed, increasing the temperature we observe how the system is able
to escape from this metastable configuration and transit to a state with four
domains, as shown in Fig. 5(b).

Figure 5: Polarization map of the relaxed structures labeled from 1 to 4 in
Fig. 4 for different lateral sizes $L$ and a layer thickness of $n=12$. Black
arrows indicate the local dipoles, projected onto the $(x,y)$ plane. The axial
component of the polarization along the $[010]_{\rm pc}$ direction is
represented by the green and magenta color map. Units of the dipoles in
$e\times$ Bohr, where $e$ is the electron charge.

The new vortices formed at the interface between PbTiO3 and SrTiO3 layers,
shown in Fig. 5(d), serve as nucleation points of a new down (respectively,
up) domains dividing the already existing up (respectively, down) polarization
regions. The formation of these new domains reduces the polarization charges
generated at the interface.

Assuming in our simulations the in-plane lattice constant of SrTiO3 [14], and
at low-enough temperatures ($T<$ 50 K), this state is a long-lived metastable
phase. The new topological defects formed at the interface are not able to
propagate and close the new domain. Increasing the temperature (beyond 90 K),
or inducing compressive strain (beyond -0.5 %) on the sample, the energy
barrier can be overcome and we observe the formation of new domains.
Interestingly, the transition to the optimal domain configuration is not
observed at high domain-densities. Below $L_{\rm c}$ the system is trapped in
a polarization wave and cannot transit to a lower-density domain configuration
by means of increasing temperature. This asymmetry comes as a consequence of
the different nature between the electrostatic and gradient energy penalties.

Molecular dynamics simulations at constant temperature show how the
recombination of vortex-antivortex pairs is the driving mechanism for the
domain propagation through the sample until the new domain is completely
formed.

Figure 6: Polarization map for a $n=14$, $L=28$ supercell at finite
temperature ($T=90$ K) and compressive epitaxial strain of $-0.5$%. (a)
Initial two-domain structure configuration. Dashed square delimits the upwards
domain studied for the dynamics in (b), dashed arrows within the SrTiO3 are a
guide to the eye to locate one of the different antivortex structures. (b)
Temporal evolution obtained by molecular dynamics simulations at finite
temperature. Red dots indicate the location of vortex while light- and dark-
blue dots indicate the location of antivortices of vorticity $-1$ and $-2$,
respectively. Meaning of the arrows and colors as in Fig. 5.

In Fig. 6(a) we show in detail the case of a $n=14$, $L=28$ at constant
temperature $T=90$ K and slight compressive strain of $-0.5$%. There we can
notice the balance of vortex and antivortex defects resulting in a zero net
vorticity on the supercell as stated by the Poincaré-Hopf theorem for our
specific periodic boundary conditions. The antivortex textures are mostly
formed at the SrTiO3 layers where the magnitude of the polarization and the
concomitant electrostatic energy of head-to-head and tail-to-tail domains is
smaller. This is in accordance with first-principles calculations [12].

In Fig. 6(b) we analyze the time evolution of a portion within the up domain
of the same superlattice [see dashed square in Fig. 6(a)], in a region where
new polarization vortices have been formed at the interfaces between the
PbTiO3 and the SrTiO3 layers. The presence of two vortices (red circles) and
an antivortex (light-blue circle) is clearly observed both at the top and the
bottom interface. Starting from this configuration these vortices and
antivortices change their shapes in order to reduce the total energy of the
system. First, the vortices elongate, while keeping their centers essentially
at the same positions. This comes with two main consequences. First, locally,
the number of unit cells with down polarization increases with time [from
three at the initial configuration to four at 300 or 400 fs, or even five at
550 fs, see dashed ovals at Fig. 6(b)]. Second, the region where the local
polarization points in-plane to close the vortices moves towards the center of
the PbTiO3 layer, and so does the center of the antivortices. At 550 fs, the
two antivortices merge to form an antivortex with topological charge -2 at the
center of the PbTiO3 layer [dark-blue point in Fig. 6(b)]. The field
disturbance doubles its charge with a high energetic cost [in accordance with
Kosterlitz-Thouless analysis in the sample XY model [28, 29] where the energy
of the vortices increases (quadratically) with the vorticity]. This is the
reason why this state is very short-lived in the molecular dynamic
simulations. Only 50 fs later, it annihilates with two vortices. In this
process a new domain with down polarization is formed, together with two new
elongated clockwise/counterclockwise pair. Finally, the new domain widens till
recovering the optimal lateral size determined by the Kittel law.

##  IV Conclusions

In summary, we theoretically extend the application of the Kittel law to the
polar vortex phase in (PbTiO3)n/(SrTiO3)n superlattices. For the explored
layer thicknesses, the square-root dependence of the domain period with the
thickness of PbTiO3 is restored by two different procedures: ($i$) full
minimization of the energy where all possible interactions are considered, and
($ii$) analyzing the harmonic force-constant bands. We find that the harmonic
approach predicts narrower domains, which is consistent with the fact that
anharmonic effects, like the development of an offset, tend to reduce the
depolarizing fields on the structure.

Moreover, studying the phase competition under suboptimal domain widths we
showed how at low-domain density new domains can be formed to relax
electrostatic constraint. These domains nucleate as vortex/antivortex pair
defects at the interfaces with SrTiO3 and propagate through the lattice by
means of recombination until the new domains are completely formed. This
recombination of vortex/antivortex is driven by the high energy costs of
polarization patterns containing vortex/antivortex pairs.

###### Acknowledgements.

F.G.-O., P.G.-F., and J.J. acknowledge financial support from Grant No.
PGC2018-096955-B-C41 funded by MCIN/AEI/10.13039/501100011033 and by ERDF “A
way of making Europe” by the European Union. F.G.-O acknowledge financial
support from grant FPU18/04661 funded by MCIN/AEI/ 10.13039/501100011033.
J.M.L. was supported by Grant No. PID2021-125543NB-I00 funded by MCIN/
AEI/10.13039/501100011033/ and by ERDF “A way of making Europe” by the
European Union. H. A. and J. I. were funded by the Luxembourg National
Research Fund through Grant C21/MS/15799044/FERRODYNAMICS. The authors
thankfully acknowledge computing time at Altamira supercomputer and the
technical support provided by the Instituto de Física de Cantabria (IFCA) and
Universidad de Cantabria (UC). The authors would also like to thank Jose Ángel
Herrero for his valuable assistance with the supercomputing environment
HPC/HTC cluster “Calderon”, supported by datacenter 3Mares, from Universidad
de Cantabria.

## References

  * Tagantsev  _et al._ [2010] A. K. Tagantsev, L. E. Cross, and J. Fousek,  _Domains in Ferroic Crystals and Thin Films_ (Springer, New York, 2010).
  * Catalan  _et al._ [2012] G. Catalan, J. Seidel, R. Ramesh, and J. F. Scott, Domain wall nanoelectronics, Rev. Mod. Phys. 84, 119 (2012).
  * Landau and Lifshitz [1935] L. Landau and E. Lifshitz, On the theory of the dispersion magnetic permeability in ferromagnetic bodies, Phys. Z. Sowjetunion 8, 153 (1935).
  * Kittel [1946] C. Kittel, Theory of the structure of ferromagnetic domains in films and small particles, Phys. Rev. 70, 965 (1946).
  * Kittel [1949] C. Kittel, Physical theory of ferromagnetic domains, Rev. Mod. Phys. 21, 541 (1949).
  * Mitsui and Furuichi [1953] T. Mitsui and J. Furuichi, Domain structure of rochelle salt and $\mathrm{K}$$\mathrm{H}_{2}$$\mathrm{P}$$\mathrm{O}_{4}$, Phys. Rev. 90, 193 (1953).
  * Catalan  _et al._ [2008] G. Catalan, H. Béa, S. Fusil, M. Bibes, P. Paruch, A. Barthélémy, and J. F. Scott, Fractal dimension and size scaling of domains in thin films of multiferroic $\mathrm{Bi}\mathrm{Fe}\mathrm{O}_{3}$, Phys. Rev. Lett. 100, 027602 (2008).
  * Bennett  _et al._ [2020] D. Bennett, M. Muñoz Basagoiti, and E. Artacho, Electrostatics and domains in ferroelectric superlattices, R. Soc. Open Sci. 7, 201270 (2020).
  * Roitburd [1976] A. L. Roitburd, Equilibrium structure of epitaxial layers, Phys. Status Solidi (a) 37, 329 (1976).
  * Lai  _et al._ [2007] B.-K. Lai, I. Ponomareva, I. Kornev, L. Bellaiche, and G. Salamo, Thickness dependency of 180∘ stripe domains in ferroelectric ultrathin films: A first-principles-based study, Appl. Phys. Lett. 91, 152909 (2007).
  * Prosandeev  _et al._ [2010] S. Prosandeev, S. Lisenkov, and L. Bellaiche, Kittel law in BiFeO3 ultrathin films: A first-principles-based study, Phys. Rev. Lett. 105, 147603 (2010).
  * Aguado-Puente and Junquera [2012] P. Aguado-Puente and J. Junquera, Structural and energetic properties of domains in PbTiO3/SrTiO3 superlattices from first principles, Phys. Rev. B 85, 184105 (2012).
  * Yadav  _et al._ [2016] A. K. Yadav, C. T. Nelson, S. L. Hsu, Z. Hong, J. D. Clarkson, C. M. Schlepütz, A. R. Damodaran, P. Shafer, E. Arenholz, L. R. Dedon, D. Chen, A. Vishwanath, A. M. Minor, L. Q. Chen, J. F. Scott, L. W. Martin, and R. Ramesh, Observation of polar vortices in oxide superlattices, Nature 530, 198 (2016).
  * Shafer  _et al._ [2018] P. Shafer, P. García-Fernández, P. Aguado-Puente, A. R. Damodaran, A. K. Yadav, C. T. Nelson, S.-L. Hsu, J. C. Wojdeł, J. Íñiguez, L. W. Martin, E. Arenholz, J. Junquera, and R. Ramesh, Emergent chirality in the electric polarization texture of titanate superlattices, Proc. Natl. Acad. Sci. U.S.A. 115, 915 (2018).
  * Wojdeł  _et al._ [2013] J. C. Wojdeł, P. Hermet, M. Ljunberg, P. Ghosez, and J. Íñiguez, First-principles model potentials for lattice-dynamical studies: general methodology and example of application to ferroic perovskite oxides, J. Phys.:Condens. Matter 25, 305401 (2013).
  * Zubko  _et al._ [2016] P. Zubko, J. C. Wojdeł, M. Hadjimichael, S. Fernandez-Pena, A. Sené, I. Luk’yanchuk, J.-M. Triscone, and J. Íñiguez, Negative capacitance in multidomain ferroelectric superlattices, Nature 534, 524 (2016).
  * García-Fernández  _et al._ [2016] P. García-Fernández, J. C. Wojdeł, J. Íñiguez, and J. Junquera, Second-principles method for materials simulations including electron and lattice degrees of freedom, Phys. Rev. B 93, 195137 (2016).
  * Gómez-Ortiz  _et al._ [2022] F. Gómez-Ortiz, P. García-Fernández, J. M. López, and J. Junquera, Melting of crystals of polarization vortices and chiral phase transitions in oxide superlattices, Phys. Rev. B 105, L220103 (2022).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1 (2015).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * Stephenson and Elder [2006] G. B. Stephenson and K. R. Elder, Theory for equilibrium 180∘ stripe domains in PbTiO3 films, J. Appl. Phys. 100, 051601 (2006).
  * Bratkovsky and Levanyuk [2009] A. M. Bratkovsky and A. P. Levanyuk, Continuous theory of ferroelectric states in ultrathin films with real electrodes, J. Comput. Theor. Nanosci. 6, 465 (2009).
  * Kornev  _et al._ [2004] I. Kornev, H. Fu, and L. Bellaiche, Ultrathin films of ferroelectric solid solutions under a residual depolarizing field, Phys. Rev. Lett. 93, 196104 (2004).
  * Stephanovich  _et al._ [2003] V. A. Stephanovich, I. A. Luk’yanchuk, and M. G. Karkut, Domain proximity and ferroelectric transition in ferro-paraelectric superlattices, Ferroelectrics 291, 169 (2003).
  * Stephanovich  _et al._ [2005] V. A. Stephanovich, I. A. Luk’yanchuk, and M. G. Karkut, Domain-enhanced interlayer coupling in ferroelectric/paraelectric superlattices, Phys. Rev. Lett. 94, 047601 (2005).
  * Zubko  _et al._ [2012] P. Zubko, N. Jecklin, A. Torres-Pardo, P. Aguado-Puente, A. Gloter, C. Lichtensteiger, J. Junquera, O. Stéphan, and J.-M. Triscone, Electrostatic coupling and local structural distortions at interfaces in ferroelectric/paraelectric superlattices, Nano Lett. 12, 2846 (2012).
  * Lu  _et al._ [2018] L. Lu, Y. Nahas, M. Liu, H. Du, Z. Jiang, S. Ren, D. Wang, L. Jin, S. Prokhorenko, C.-L. Jia, and L. Bellaiche, Topological defects with distinct dipole configurations in $\mathrm{Pb}\mathrm{Ti}\mathrm{O}_{3}$/$\mathrm{Sr}\mathrm{Ti}\mathrm{O}_{3}$ multilayer films, Phys. Rev. Lett. 120, 177601 (2018).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, J. Phys. C: Solid State Phys. 6, 1181 (1973).
  * Kosterlitz [1974] J. M. Kosterlitz, The critical properties of the two-dimensional $\mathrm{XY}$ model, J. Phys. C: Solid State Phys. 7, 1046 (1974).