# Theory of polarization textures in crystal supercells

Daniel Bennett<EMAIL_ADDRESS>John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts
02138, USA  Wojciech J. Jankowski  Theory of Condensed Matter Group, Cavendish
Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE,
United Kingdom  Gaurav Chaudhary  Theory of Condensed Matter Group, Cavendish
Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE,
United Kingdom  Efthimios Kaxiras  John A. Paulson School of Engineering and
Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA
Department of Physics, Harvard University, Cambridge, Massachusetts 02138,
United States  Robert-Jan Slager  Theory of Condensed Matter Group, Cavendish
Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE,
United Kingdom

###### Abstract

Recently, topologically non-trivial polarization textures have been predicted
and observed in nanoscale systems. While these polarization textures are
interesting and promising in terms of applications, their topology in general
is yet to be fully understood. For example, the relation between topological
polarization structures and band topology has not been explored, and polar
domain structures are typically considered in topologically trivial systems.
In particular, the local polarization in a crystal supercell is not well-
defined, and typically calculated using approximations which do not satisfy
gauge invariance. Furthermore, local polarization in supercells is typically
approximated using calculations involving smaller unit cells, meaning the
connection to the electronic structure of the supercell is lost. In this work,
we propose a definition of local polarization which is gauge invariant and can
be calculated directly from a supercell without approximations. We show using
first-principles calculations for commensurate bilayer hexagonal boron nitride
that our expressions for local polarization give the correct result at the
unit cell level, which is a first approximation to the local polarization in a
moiré superlattice. We also illustrate using an effective model that the local
polarization can be directly calculated in real space. Finally, we discuss the
relation between polarization and band topology, for which it is essential to
have a correct definition of polarization textures.

##  I Introduction

The formation of complex polar structures such as ferroelectric domains Mitsui
and Furuichi (1953); Kopal  _et al._ (1997), analogous to ferromagnetic
domains Kittel (1946), is a phenomenon intrinsic to ferroelectric materials
with finite boundary conditions, and has been studied for many years. When
going from bulk to a lower dimensional system, interfaces between a
ferroelectric material and a non-polar material or vacuum lead to polar
discontinuities which if not screened will lead to depolarizing fields that
suppress ferroelectricity Junquera and Ghosez (2003). Ferroelectric materials
can form polydomain structures to mitigate depolarization effects, such as
$180^{\circ}$ stripe domains, which for not too thin films can be described
using a Landau-like theory Bennett  _et al._ (2020). These sharp domain
structures eventually become unstable in very thin films, at which point it
becomes more favourable to form softer domain walls, better described by a
Ginzburg-Landau theory Luk’yanchuk  _et al._ (2009); such domains have been
observed in ferroelectric materials down to the monolayer limit Ji  _et al._
(2019). In some cases, even more complex structures involving polar vortices
may form Gómez-Ortiz  _et al._ (2023), for example in
ferroelectric/paraelectric (FE/PE) superlattices Callori  _et al._ (2012);
Zhang  _et al._ (2017); Dawber (2017); Park  _et al._ (2018); Susarla  _et
al._ (2021) such as $\mathrm{PbTiO_{3}/SrTiO_{3}}$ (PTO/STO), where the
properties can be tuned through the use of different materials and by tuning
the relative thicknesses of the layers. In FE/PE superlattices, both stripe
domains and vortices can form Dawber (2017), depending on the strength of the
coupling between the ferroelectric layers Luk’yanchuk  _et al._ (2009), which
can be determined by the ratio of the layers and the dielectric permittivities
Bennett  _et al._ (2020).

Ferroelectric materials have been fabricated in many different geometries,
from 2D thin films and FE/PE superlattices to 1D nanowires Urban  _et al._
(2002); Yun  _et al._ (2002) and nanotubes Luo  _et al._ (2003); Morrison  _et
al._ (2003); Mao  _et al._ (2003), and even 0D quantum dots Chu  _et al._
(2004); Shin  _et al._ (2005). Lower-dimensional ferroelectric systems
typically exhibit size-dependent transitions where the polarization fields
become more complex, before eventually becoming unstable and vanishing
completely when the paraelectric phase is favored Fu and Bellaiche (2003);
Naumov  _et al._ (2004); Geneste  _et al._ (2006); Morozovska  _et al._
(2006). Soon after these complex polarization textures were discovered, they
were proposed to be topologically non-trivial Hong  _et al._ (2010); Junquera
_et al._ (2023); skyrmion-like polarization structures were identified, for
example in $\mathrm{BaTiO_{3}}$ (BTO) nanowires embedded in a matrix of STO
Nahas  _et al._ (2015). It has also been proposed that 3D skyrmions, i.e.
hopfions, may be created by controlling domains and domain walls in
ferroelectrics, where at low temperatures the polarization rotates in-plane in
the domain walls Pereira Gonçalves  _et al._ (2019), resulting in a non-
trivial winding. Ferroelectric skyrmions have been experimentally observed in
PTO/STO superlattices Das  _et al._ (2019); Han  _et al._ (2022), and a
skyrmion $\to$ meron transition with strain has recently been observed Shao
_et al._ (2023). Polar merons have also been observed in PTO under epitaxial
strain from a $\mathrm{SmScO_{3}}$ substrate Wang  _et al._ (2020).

Figure 1: Sketch of a twisted hBN bilayer. The MPDs in a single moiré cell are
highlighted in blue and red, which have equal and opposite polarization. The
normalized polarization field in the MPDs are sketched above, which form
antimerons ($Q=-\frac{1}{2}$) and merons ($Q=+\frac{1}{2}$), respectively.

In recent years, a new type of ferroelectricity in layered van der Waals
materials has been proposed Li and Wu (2017); Bennett and Remez (2022);
Bennett (2020) and experimentally observed Zheng  _et al._ (2020); Yasuda  _et
al._ (2021). In layered systems such as 3R-stacked hexagonal boron nitride
(hBN) or $\mathrm{MoS_{2}}$, an out-of-plane polarization occurs via an
interlayer transfer of electronic charge when the relative stacking between,
which can be switched by a relative sliding between the layers (van der Waals
sliding Stern  _et al._ (2021)), resulting in ferroelectricity. When there is
a relative twist or lattice mismatch between the layers, forming a moiré
superlattice (see Fig. 1), the interlayer charge transfer results in an out-
of-plane polarization texture, and the stacking domains can be identified as
moiré polar domains (MPDs), which have been experimentally shown to result in
ferroelectricity Yasuda  _et al._ (2021) via the growing and shrinking of the
MPDs in response to an applied field Bennett and Remez (2022); Bennett (2020).
Recently, we discovered that the MPDs also have an in-plane component of
polarization Bennett  _et al._ (2023), resulting in possible real space non-
trivial topological textures. For hBN and similar materials the topological
charge $Q$ in each stacking domain,

$Q=\frac{1}{4\pi}\int{\bf P}\cdot\left(\partial_{x}{\bf
P}\times\partial_{y}{\bf P}\right)\mathop{}\\!\mathrm{d}{\bf r}\;,$ (1)

where the local polarization ${\bf P}({\bf r})$ is normalized, integrates to
$\pm\frac{1}{2}$, meaning the MPDs are polar merons and antimerons (see Fig.
1). Moiré superlattices have also recently been fabricated with perovskites
Shen  _et al._ (2022), which for FE/PE interfaces such as BTO/STO also result
in polar vortices Sanchez-Santolino  _et al._ (2023).

While the complex structures and non-trivial topology of polar nanostructures
are very interesting and promising for future applications in nanotechnology,
some of their fundamental properties require a better understanding, namely:

(i) The origin of topological polarization structures: Polarization can be
topological in the real space sense, i.e. there may be non-trivial winding of
the polarization vector field of a system. This winding must arise from both
the geometry of the system and underlying crystal symmetries of the
constituent materials. For example, in perovskites superlattices, the
electrostatic boundary conditions at the interfaces between alternating layers
make a uniform polarization in the ferroelectric layers energetically
unfavourable. This promotes a polydomain structure, but does not guarantee
non-trivial topology; a 180∘ stripe domain structure is preferable to a
monodomain polarization, but is topologically trivial. The polar modes in each
unit cell, related to the off-centering of the Ti atoms within the oxygen
octahedra, lower the crystal symmetry, and allow coupling between polarization
and strain such as piezoelectricity, which is forbidden by symmetry in the
cubic phase. As a result, the strain across the domain walls and interfaces
results in an additional component of polarization Pereira Gonçalves  _et al._
(2019), which causes the polarization to wind and become topologically non-
trivial. Rolling a thin film into a nanotube induces a radial polarization and
hence a polar vortex via flexoelectricity Artyukhov  _et al._ (2020);
Springolo  _et al._ (2021); Bennett (2021), which is a property of all
insulators. In twisted or strained hBN, the MPDs are a result of the geometry
of the moiré superlattice. However, the reason that the polarization winds and
is topologically non-trivial is related to the underlying symmetry of the hBN
bilayer Bennett  _et al._ (2023). It was shown with a space group analysis
that different mirror planes are broken by different local stackings, meaning
that the polarization must point out-of-plane in the domain centers, and in-
plane along the domain walls, leading to a network of merons and antimerons.

When considering the topology of polarization, one major problem is that the
local polarization in a crystal in real space is not a well-defined quantity.
From the modern theory of polarization Vanderbilt and King-Smith (1993); King-
Smith and Vanderbilt (1993); Resta (1994); Nakagawa  _et al._ (2020);
Vanderbilt (2018); Van Mechelen  _et al._ (2022), the change in total
polarization of a system can be measured from Berry phases, but the
decomposition into individual contributions in real space is somewhat
arbitrary. There are two types of methods for estimating the local
polarization in a supercell: ‘configuration space’ methods, where the local
configurations in a each unit cell are emulated in a commensurate system, and
real space methods, where the the local polarization is calculated directly
from calculations involving the entire supercell. It is clear that calculating
local polarization directly in real space would be preferable, although this
would require expensive calculations, and defining local polarization in real
space is difficult. One proposal is to measure the individual Wannier centers
in each unit cell Wu  _et al._ (2006), which has been applied to slab-like
systems to measure the out-of-plane polarization, averaged in the in-plane
directions. However, this approximation does not necessarily satisfy gauge
invariance. The more commonly used approximation relies on obtaining the local
displacements in each cell and multiplying them by the Born effective charges
obtained from the bulk system Meyer and Vanderbilt (2002); Stengel  _et al._
(2011). This was originally proposed for slab-like systems, and is the
standard method for calculating the local polarization in perovskite systems.
However, this relies on the assumption that the Born effective charges are
uniform everywhere in space, which is not the case in twisted bilayers Bennett
_et al._ (2023), for example.

(ii) The relation to electronic band topology: Local polarization is not
typically described at the electronic structure level, since it is defined in
terms of local displacements in real space. Thus, the relation between the
topology of polarization textures and band topology Qi and Zhang (2011); Hasan
and Kane (2010), which has seen the development of various analytical
diagnoses Fu (2011); Slager  _et al._ (2013); Ünal  _et al._ (2020); Po  _et
al._ (2017); Bradlyn  _et al._ (2017); Kruthoff  _et al._ (2017); Shiozaki and
Sato (2014); Bouhon  _et al._ (2020a, b); Song  _et al._ (2018), is not clear.
In the original description of polarization in terms of localized Wannier
functions, the polarization is well-defined only for topologically trivial
bands; for topologically non-trivial bands, the Wannier functions are not
exponentially localized. The modern theory of polarization has since been
generalized to Chern insulators by tracking the hybrid Wannier centers
throughout closed loops in the Brillouin zone (BZ) Coh and Vanderbilt (2009);
Song  _et al._ (2021) (equivalent to calculating Wilson loops Alexandradinata
and Bernevig (2016); Bouhon  _et al._ (2019)). However, the similarities and
interplay between polar and band topology are not known. Recently there has
been evidence that polarization can affect the topological properties of a
system: in layered $\mathrm{MnBi_{2}Te_{4}}$, which is antiferromagnetic,
ferroelectric, and exhibits quantum anomalous Hall (QAH) conductance,
inverting the polarization via van der Waals sliding changes the sign of the
Chern number and hence the QAH conductance Liang  _et al._ (2023). Thus, a
better understanding of the role that polarization can play in the properties
of systems arising from band topology is needed.

(iii) Physical consequences: Perhaps most importantly, the physical
consequences of topological polarization are not well-understood. While it is
has recently been shown that ferroelectric switching can change the Chern
numbers in topologically non-trivial systems Liang  _et al._ (2023), this is a
result of uniform polarization; there may be additional phenomena which are
unique to topologically non-trivial polar textures. Furthermore, finding
physical consequences may also provide new ways to indirectly detect polar
topology experimentally, which currently requires very careful microscopy
measurements.

In this work, we address the problem of defining the local polarization in a
crystal supercell, in the context of studying topological polarization. The
most common approach for estimating local polarization, in twisted bilayers
for example Bennett and Remez (2022); Bennett (2020); Bennett  _et al._
(2023), is using the configuration space mapping: the local polarization in
real space is approximated by the total polarization in configuration space,
i.e. a commensurate system with a global distortion such as a relative shift
between the layers. Estimating local polarization in this way requires the
polarization to vary slowly, with a wavelength similar to the supercell
period. Additionally, information about the electronic structure is lost, as
the local polarization is not calculated from the bands of the supercell.
While approximate expressions for the local polarization which can be
calculated directly in real space using Wannier centers and Born effective
charges have been proposed in the literature Wu  _et al._ (2006); Meyer and
Vanderbilt (2002); Stengel  _et al._ (2011), they are not necessarily well-
defined. In the former case, the local polarization is given by a partial sum
over Wannier centers, which in general is not gauge invariant. In the latter
case, partial sums over the Born effective charges may not satisfy the
acoustic sum rule, or charge neutrality. We discuss the different ways in
which local polarization can be estimated, and propose definitions of local
polarization in terms of Born effective charges or Wannier centers, which are
well-defined (gauge invariant) and can be evaluated directly in real space. As
a stepping stone to real space calculations, we show that these definitions
yield the correct polarization in commensurate 3R-stacked hBN in configuration
space. We illustrate using effective models in 1D (Aubry-André model in the
continuum limit Szabó and Schneider (2018); Borgnia  _et al._ (2022); Borgnia
and Slager (2023); Aubry and André (1980)) and 2D (Bistritzer-MacDonald
Bistritzer and MacDonald (2011); Balents (2019)) that our proposed definitions
can be used to calculate the local polarization in real space, without relying
on the mapping to configuration space. Finally, we discuss the relation
between polarization and band topology.

##  II Local polarization

Our aim is to develop a definition of the local polarization field in a
supercell comprised of a number of repeated unit cells, each with real space
position ${\bf r}_{j}$. Each unit cell may have local distortions with respect
to a more symmetric reference configuration, and we define the local
polarization as a result of these local distortions as a discrete vector field
${\bf P}({\bf r}_{j})$, valued in each unit cell. When studying systems
described by large supercells, in particular moiré superlattices, it is
generally much more convenient to work in terms of the local distortions in
each unit cell, which form a configuration space, rather than the full
supercell in real space.

###  II.1 Configuration Space

Figure 2: (a): Sketch of a $5\times 1\times 1$ supercell of a titanium oxide
perovskite $\mathrm{ATiO_{3}}$. The arrows indicate the displacement of the Ti
atom from the center of the oxygen octahedra in each unit cell. (b): Mapping
of the displacements in each unit cell from real space to configuration space,
which has the dimensions of a single 5-atom unit cell. (c): $3\times 3$
section of a bilayer with relative twist angle $\theta$. (d): Mapping of the
relative displacements between the atoms in each layer from real space to
configuration space, which has the dimensions of a single bilayer unit cell.

The mapping between real space and configuration space is used to estimate
local properties in periodic supercells and incommensurate systems Carr  _et
al._ (2017); Cazeaux  _et al._ (2017); Massatt  _et al._ (2017); Carr  _et
al._ (2018). In the context of polar domains and textures, the two most common
examples are oxide perovskite films (BTO, PTO, etc.), in which the supercell
composed of five unit cells is repeated in the direction(s) normal to the
interface in order to allow for the formation of polydomain structures, see
Fig. 2 (a), and a twisted/strained bilayer, in which two layers are twisted or
strained relative to one another to form a moiré superlattice, see Fig. 2 (c).

In an $\mathrm{ABO_{3}}$ perovskite, the polarization arises from the off-
centering between the B cations and the $\mathrm{O_{6}}$ octahedra. In
titanate oxide perovskites such as PTO and BTO, this has been attributed to
the hybridisation of the Ti $3d$ and O $2p$ orbitals Cohen (1992). For a
perovskite supercell with a polydomain structure, the local polarization is a
result of the polar mode displacements in each unit cell. Assuming for
simplicity that the polar modes are dominated by the off-centering of the Ti
atoms within each oxygen octahedron, the displacement in each unit cell is
approximately ${{\bf x}({\bf r}_{j})=\Delta{\bf r}_{j}^{\mathrm{Ti}}}$, where
${\bf r}_{j}$ is the position of unit cell $j$, and ${\bf
r}_{j}^{\mathrm{Ti}}$ is the position of the Ti atom in the cell. The
displacements in each unit cell can be mapped to a single space known as
configuration space, which is the size of a single 5-atom unit cell, see Fig.
2 (b). Thus, local properties such polarization can be estimated in
configuration space using a bulk 5-atom unit cell, which is much more
efficient than performing calculations using the full supercell. For example,
the polarization can be parameterized for every point in configuration space
by calculating the Berry phases for a 5-atom bulk cell with various
displacements of the Ti atom. Assuming the displacements in real space are
known, the local polarization can then be parameterized in real space by using
an inverse mapping from configuration space. While the mapping of
displacements between real space and configuration space is exact, the
parameterization of local properties relies on the approximation that the
displacements vary slowly in real space and can therefore be neglected; we say
that the local changes around a given cell are small, and we neglect them,
which gives a commensurate cell which is much more easy to simulate. Of
course, this approximation is not good for properties which modulate at the
unit cell level, such as Peierls distortions, defects, dislocations, etc.
However, for quantities which vary smoothly on the supercell scale, the
approximation works well; in moiré materials for example, the separation of
supercell and single-layer unit cell length scales justifies such approximate
approach relying on the continuity of the mapping.

While the local polarization can be parameterized in configuration space by
calculating Berry phases as a function of displacements, it is typically
estimated as follows for perovskite systems: the displacements ${\bf x}({\bf
r}_{j})$ are measured in real space, either experimentally or from molecular
dynamics simulations. The displacements in each cell are contracted with the
Born effective charge tensor of the bulk cubic phase:

$P_{\beta}({\bf r}_{j})\approx\frac{1}{\Omega_{0}}\sum_{\kappa\in{\bf
r}_{j}}\sum_{\alpha}Z^{*}_{\kappa,\alpha\beta}x_{\kappa,\alpha}\;,$ (2)

where $\Omega_{0}$ is the unit cell volume, contrary to the total supercell
volume $\Omega_{\text{sc}}$. The first sum is over the atoms $\kappa$ in unit
cell ${\bf r}_{j}$, and the second is over Cartesian directions $\alpha$. This
is the generally accepted method to measure local polarization in oxide
perovskites, and has been used to predict topologically non-trivial
polarization. One problem with Eq. (2) is that the individual unit cells are
not guaranteed to satisfy the acoustic sum rule. In slab-like systems, the
polarization in each layer is typically averaged over neighbouring layers
Stengel  _et al._ (2011), and is therefore not truly localized to each unit
cell. Eq. (2) also requires the Born effective charges to be uniform in
configuration space, which is reasonable when the polarization varies slowly
in real space, such as in the centers of the polar domains, but not when the
polarization varies sharply, such as across domain walls; in PTO/STO
superlattices for example, the domain walls have been shown to be a single
unit cell wide Gómez-Ortiz  _et al._ (2023). Additionally, it is possible for
ferroelectric domain walls to be conducting, in which case it is not clear if
the Born effective charges are well defined Sluka  _et al._ (2013); Bednyakov
_et al._ (2015).

One situation in which Eq. (2) fails to correctly estimate the local
polarization is in van der Waals materials, where the Born effective charges
vary nonlinearly as one layer slides over the other Bennett  _et al._ (2023).
Although the charge transfer and modulation of the Born effective charges is
small, they give rise to the unique polarization textures in twisted bilayers
and thus cannot be neglected.

For a bilayer with relative twist angle $\theta$ between the layers, see Fig.
2 (c), the configuration space mapping is given by Carr  _et al._ (2018)

${\bf x}({\bf r})=(\mathbb{I}-R_{\theta}^{-1}){\bf r}\;,$ (3)

modulo any lattice vectors, where ${\bf r}$ is the real space position and
$R_{\theta}$ is a rotation matrix. While Eq. (3) is exact, for small twist
angles the local changes in environment around the black unit cells are small,
and the local properties in each cell can be described by a commensurate
bilayer with a relative translation ${{\bf
x}\approx\theta\begin{bmatrix}0&-1\\\ 1&0\end{bmatrix}{\bf r}}$ between the
layers, see Fig. 2 (d). Similarly, for a small homogeneous strain $\eta$, the
equivalent mapping is ${{\bf x}=\eta{\bf r}}$. This allows the local
properties in strained or small-angle twisted bilayers to be parameterized
efficiently with first-principles calculations using a single commensurate
cell of a bilayer, and sliding one layer over the other. Importantly, what we
define as the unit cell ${\bf r}_{j}$ for moiré systems, is a unit cell of one
layer and the atoms of the other layer which are contained in that unit cell.
Such a general projective description is consistent with the configuration
space picture, and provides an arena to define a local quantity, given its
shortsightedness with respect to the neighbouring unit cells. However, within
this picture, the configuration space calculations for simultaneously twisted
and strained moiré bilayers performed on an elastically deformed unit cell
might yield local polarization significantly deviating from the real values,
on applying periodic boundary conditions in configuration space. The main
reason for this is that when both unit cell deformation and the stacking
modulation due to the twist are present, the repetition of such approximate
unit cells in configuration space for computational purposes, implicit in the
imposed boundary conditions, does not account for the appreciable structural
variations in the neighbourhood of the studied unit cell. Therefore, such
possibilities also motivate pursuing the definitions of local polarization in
crystal superlattices, beyond the notion of configuration space.

The local polarization can be calculated in the configuration space by sliding
one layer over the other and directly calculating the Berry phases Bennett and
Remez (2022); Bennett (2020); Bennett  _et al._ (2023). However, the Berry
phase obtained for each point in configuration is not physically meaningful.
In real space, the Berry phase is a global property of the system, which
yields the total polarization. A more natural way to define the local
polarization is using the Born effective charges, since they are locally well-
defined in real space. When the Born effective charges are not constant,
Eq.(2) is not valid. Instead, the Born effective charges must be integrated:

$P_{\beta}({\bf r}_{j})\approx\frac{1}{\Omega_{0}}\int_{0}^{{\bf x}({\bf
r}_{j})}Z^{*}_{\kappa,\alpha\beta}({\bf
x}^{\prime})dx^{\prime}_{\kappa,\alpha}\;,$ (4)

where summation over repeated indices is assumed, and again only the atoms
$\kappa\in{\bf r}_{j}$ are displaced. When the Born effective charges are
constant, Eq. (4) simplifies to Eq. (2). The integration is performed from a
reference state 0 to a general configuration ${\bf x}$, where in configuration
space we can change smoothly from one to the other via a relative translation
between the layers. In order to obtain the polarization, rather than an
arbitrary change of polarization between two configurations, the reference
state is chosen to be non-polar. For hBN and similar materials, the only non-
polar configuration is when the layers are perfectly aligned (AA stacking) and
unstrained, and therefore, while Eq. (4) is naturally defined on a torus, the
normalized polarization, used in Eq. (1), is defined on a punctured torus,
similarly to the strain fields Engelke  _et al._ (2023). Eq. (4) yields a
polarization identical to the one obtained from Berry phases Bennett  _et al._
(2023).

While Eqs. (2) and (4) have been successfully used to estimate polarization
textures, they are both only valid for certain approximations such as large
supercells, smoothly varying quantities and constant Born effective charges in
the former case. Furthermore, we also lose any information about the
electronic structure when using this approximation, because the electronic
bands for each point in configuration space are not meaningful in real space.
This is especially apparent when calculating the Berry phases; the Berry phase
is a global property of a system, which we calculate from the bands for each
point in configuration space. Although we get a good approximation to the
local polarization, we do not have any information about the contribution from
different bands of the supercell.

###  II.2 Defining local polarization

From the modern theory of polarization, the total polarization ${\bf
P}^{\text{tot}}$ in a crystal is given by King-Smith and Vanderbilt (1993);
Vanderbilt and King-Smith (1993); Vanderbilt (2018)

${\bf
P}^{\text{tot}}=\frac{-ief}{(2\pi)^{3}}\sum_{n}^{\text{occ}}\oint_{\text{BZ}}\braket{u_{n,{\bf
k}}}{\nabla_{{\bf k}}u_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\;,$ (5)

where $\ket{u_{n,{\bf k}}}$ are the cell-periodic parts of the Bloch
wavefunctions and $f$ is the occupation number of states in the valence bands
(2 for spin-degenerate systems). Eq. (6) can be rewritten as

${\bf
P}^{\text{tot}}=\frac{-ef}{(2\pi)^{3}}\sum_{n}^{\text{occ}}\oint_{\text{BZ}}{\bf
A}_{n}({\bf k})\mathop{}\\!\mathrm{d}{\bf
k}=\frac{-ef}{\Omega}\sum_{n}^{\text{occ}}\frac{1}{2\pi}\phi_{n,\alpha}{\bf
a}_{\alpha}\;,$ (6)

where ${\bf a}_{\alpha}$ are the lattice vectors, $\Omega$ is the system cell
volume.

${\bf A}_{n}({\bf k})=i\braket{u_{n,{\bf k}}}{\nabla_{{\bf k}}u_{n,{\bf k}}}$
(7)

is the Berry connection of band $n$, and

$\phi_{n,\alpha}=\frac{i\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}\braket{u_{n,{\bf
k}}}{{\bf b}_{\alpha}\cdot\nabla_{{\bf k}}u_{n,{\bf
k}}}\mathop{}\\!\mathrm{d}{\bf k}$ (8)

is the Berry phase of band $n$ in direction $\alpha$, defined up to a factor
of $2\pi$, where ${\bf b}_{\alpha}$ are the reciprocal lattice vectors. For
disordered systems, where ${\bf k}$ is not a good quantum number, a
generalization for the computation of macroscopic total polarization with
single-point Berry phases has been proposed by Resta Resta (1998).

It is well-known that the absolute polarization in a crystal is not well-
defined. Only changes in polarization are well-defined, modulo any quanta of
polarization (integer values of the Berry phases in different directions).
Derivatives of the polarization with respect to perturbations (phonon, strain,
electric field), which can be identified as the dielectric and
electromechanical properties of a system, are well-defined and are routinely
calculated from first-principles calculations, either using finite difference
methods or density functional perturbation theory (DFPT) Wu  _et al._ (2005):
the dielectric response of a system is related to the derivative of the total
polarization with respect to electric field Gonze and Lee (1997), and
electromechanical properties can be measured by calculating the derivatives of
the polarization with respect to strain (piezoelectricity Vanderbilt (2000)),
or strain and electric field (electrostriction Bennett  _et al._ (2022)). For
our purposes, the most relevant property is the derivative of the polarization
with respect to phonon displacements, i.e. the Born effective charge tensor
Gonze and Lee (1997); Ghosez  _et al._ (1998):

$Z^{*}_{\kappa,\alpha\beta}=\Omega\frac{\partial P_{\beta}}{\partial
x_{\kappa,\alpha}}=\frac{\partial
F_{\kappa,\alpha}}{\partial\mathcal{E}_{\beta}}\;,$ (9)

where $x_{\kappa,\alpha}$ is the real space (phonon) displacement of atom
$\kappa$ in direction $\alpha$, $F_{\kappa,\alpha}$ is the force on atom
$\kappa$ in direction $\alpha$, and $\mathcal{E}_{\beta}$ is an electric field
in direction $\beta$. Because the Born effective charge tensor is related to
the mixed derivative of the free energy of the system with respect to phonon
displacement and electric field, it can be interpreted as both the dipole
generated by a phonon displacement, and the force generated on an atom by an
electric field.

As mentioned previously, the change in polarization from one configuration to
another can be obtained by integrating the Born effective charges using Eq.
(4). For twisted/strained bilayers, this was done in configuration space in
order to avoid expensive calculations involving large supercells Bennett  _et
al._ (2023). However, we propose that the local polarization may be calculated
in a more well-defined way by calculating the Born effective charges in real
space. While the polarization in configuration space is simply approximate,
the Born effective charges are well-defined in real space because they are the
derivatives of the well-defined total polarization of the supercell with
respect to the local atomic displacements in each unit cell. Expressed in this
way, the Born effective charges are obtained directly from the electronic
bands of the supercell, rather than the fictitious electronic bands in
configuration space.

We define a unit cell ${\bf r}_{j}$ as the smallest structural unit that can
be mapped to configuration space, which captures the complete set of all
possible configurations, i.e. spanning the entire system. Writing the
dynamical charges in each unit cell as Ghosez and Gonze (2000); Sai  _et al._
(2002)

$Z^{*}_{\kappa,\alpha\beta}({\bf x}({\bf
r}_{j}))=\frac{-2ief\Omega_{\text{sc}}}{(2\pi)^{3}}\sum_{n}^{\text{occ}}\oint_{\text{scBZ}}\braket{\partial_{x_{\kappa,\alpha}}u_{n,{\bf
k}}}{\partial_{k_{\beta}}u_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\;,$ (10)

for all atoms $\kappa$ in cell ${\bf r}_{j}$, the local polarization in each
unit cell is given by

$P_{\beta}({\bf r}_{j})=\frac{-2ief}{(2\pi)^{3}}\int_{0}^{{\bf x}({\bf
r}_{j})}\sum_{n}^{\text{occ}}\oint_{\text{scBZ}}\braket{\partial_{x_{\kappa,\alpha}}u_{n,{\bf
k}}}{\partial_{k_{\beta}}u_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf
k}\mathop{}\\!\mathrm{d}x^{\prime}_{\kappa,\alpha}\;.$ (11)

It is important to stress that the momentum space integral is evaluated over
the supercell Brillouin zone (scBZ) and the summation is performed over the
bands of the supercell, while the integral with respect to ${\bf x}$ is
performed over the relative displacements not in real space, but in
configuration space, where the different configurations are connected by the
phonon displacements ${\bf x}$. The limits of the integral are the nonpolar
reference state 0 and the local configuration in each unit cell ${\bf r}_{j}$.

The key difference between Eq. (11) and Eq. (4) is that the Born effective
charges are calculated correctly: in real space, and using the electronic
bands of the supercell. There are a few subtle details associated with the
definition of local polarization in this way. First, the system must be a
supercell comprised of a number of smaller unit cells, within each a local
polarization is defined. The local polarization in each cell is really defined
as a change in polarization with respect to a reference cell, but taking the
reference cell to be non-polar, we write the polarization as $P$ rather than
$\Delta P$. The integral over configurations in Eq. (11) could be mapped to
positions in real space. However, in real space the polarization is a discrete
vector field, defined in each unit cell of the supercell, whereas
configuration space is generally continuous and simply connected. A
commensurate supercell is mapped to a finite subset of configuration space
containing a discrete set of points, but for an incommensurate supercell,
where the period goes to infinity, there is a one to one mapping between the
two spaces. The integral over ${\bf x}$ in Eq. (11) should be discretized over
the unit cells in real space, but we can take advantage of the fact that
configuration space is continuous and interpolate the Born charges (or any
local quantity), and obtain the local polarization field in configuration
space which is continuous and varies smoothly.

###  II.3 Local polarization from a 2D continuum model

In the configuration space method the polarization is computed as a global
quantity in each configuration and then related to the local polarization via
a mapping between configuration and the real space. This approach has its
benefit as the global polarization is well defined for each point in
configuration space. However, physically its relation to the local
polarization becomes less transparent. The direct real-space picture described
in the previous section is obtained at the unit cell level. However in moiré
superlattices it is often more convenient to work with continuum models. While
superlattices constructed from microscopic unit-cells are only well-defined at
some commensurate twist angles or strains, the continuum model description
accurately captures the low energy physics at all small angles and strains
with smooth moiré periods, irrespective of microscopic periodicity Bistritzer
and MacDonald (2011).

Here, we show that an expression for local polarization can be obtained in
real space using a continuum field approach, derived in the context of
deformation fields Balents (2019). For illustrative purposes we consider a
moiré bilayer formed by a small strain or twist angle, although the
generalization to the more complicated moiré patterns is straightforward.

The continuum model describes the low energy physics near a band extremum that
is located at the momentum ${\bf K}$ of the un-deformed monolayer. The
electron field at ${\bf K}$ is given by

$c({\bf x})=\psi({\bf x})e^{i{\bf K}\cdot{\bf x}}\;.$ (12)

For concreteness, we assume that near this extremum the low energy physics of
the monolayer is described by a 2D massive-Dirac model. This is appropriate
for bilayer hBN, the main example considered in this work and in Ref. Bennett
_et al._ (2023), although generalization to other models and dimensions is
straightforward. The monolayer Hamiltonian is given by

$H_{\text{ML}}=\int\psi^{\dagger}\left(m\tau_{3}-iv\tau^{\mu}\partial_{{\bf
x}_{\mu}}\right)\psi\mathop{}\\!\mathrm{d}^{2}{\bf x}\;,$ (13)

where $m$ is the mass gap, $v$ is the Dirac-velocity and the $\tau$ Pauli-
matrices act on some internal degrees, which for hBN are the two sublattices.
Summation is assumed, with $\mu=1,2$.

The effect of small strain and twist can be captured by a local deformation
$\mathcal{D}({\bf r})$ field as

${\bf r}={\bf x}+\mathcal{D}({\bf r})\;,$ (14)

where ${\bf r}$ is real space position in a “laboratory frame”, ${\bf x}$ is
the position in the monolayer and $\mathcal{D}({\bf r})$ is a deformation
field as a result of strain and twist. The deformation field is assumed to be
locally small, i.e. $\partial_{\mu}\mathcal{D}\ll 1$. It is crucial to define
the “small” deformation $\mathcal{D}$ as a function of the variable ${\bf r}$
instead of the variable ${\bf x}$, since when considering the multilayers,
${\bf x}$ are associated with the individual layers correspond to very
different locations ${\bf r}$ in the real space. Thus one cannot define a
“small local” deformation field in the variable ${\bf x}$. This is equivalent
to the configuration-space consideration where configuration space vector is
simply the local change in configurations of the unit cells from the top and
the bottom layers, instead of overall global shift of the particular unit cell
as one of the layer is twisted or strained.

As a result of a small deformation in the layer, the electron field is locally
modified as

$c({\bf r})=\left|\text{det}\left(\frac{\partial x_{\mu}}{\partial
r_{\nu}}\right)\right|^{1/2}c({\bf x}({\bf r}))\;.$ (15)

The $\psi$ field is correspondingly modified as

$\psi({\bf r})=\left(1-\bm{\nabla}\cdot\mathcal{D}({\bf
r})\right)^{1/2}\psi({\bf x}({\bf r}))e^{-i{\bf K}\cdot\mathcal{D}({\bf
r})}\;,$ (16)

and the integral measure is modified as

$\mathop{}\\!\mathrm{d}^{2}{\bf x}=\text{det}\left(\frac{\partial
x_{\mu}}{\partial r_{\nu}}\right)\mathop{}\\!\mathrm{d}^{2}{\bf
r}\sim(1-\bm{\nabla}\cdot\mathcal{D}({\bf r}))\mathop{}\\!\mathrm{d}^{2}{\bf
r}\;.$ (17)

The continuum Hamiltonian of the decoupled bilayers can be obtained when each
layer experiences an independent deformation field $\mathcal{D}_{l}$, where
$l=\text{t,b}$ is the layer index:

$\begin{split}H_{\text{BL}}&=\sum_{l=\text{t,b}}\int\psi^{\dagger}_{l}\biggl{[}m\tau_{3}-iv\left((\tau^{\mu}+\frac{\partial\mathcal{D}_{l,\mu}}{\partial
r_{\nu}}\tau^{\nu}\right)\partial_{r_{\mu}}\\\ &\hskip 56.9055pt+v({\bf
K}\cdot\partial_{r_{\mu}}\mathcal{D}_{l})\tau^{\mu}\biggr{]}\psi_{l}\mathop{}\\!\mathrm{d}^{2}{\bf
r}\end{split}\;,$ (18)

where we have kept only the terms linear in the deformation field. For
example, if the two layers are twisted rigidly by angles $\pm\theta/2$, the
deformation field is given by

$\mathcal{D}_{\text{t}}({\bf r})=-\mathcal{D}_{\text{b}}({\bf
r})=\frac{\theta}{2}\hat{z}\times{\bf r}\;.$ (19)

The twist deformation is simply the continuum field analog of the real space
to configuration space mapping for a twisted moiré bilayer in Eq. (3). The
inter-layer tunneling Hamiltonian is also modified under the twist deformation
and takes a general form

$H_{\text{tun}}=\int\psi^{\dagger}_{2}T(\mathcal{D}_{\text{t}}-\mathcal{D}_{\text{b}})\psi_{1}\mathop{}\\!\mathrm{d}^{2}{\bf
r}\ +\ \text{h.c.}\;,$ (20)

Here, we have assumed that the interlayer tunneling $T$ is purely local and
does not depend on the gradients of the deformation fields. Furthermore, if
the two layers are deformed identically, the local interlayer tunneling must
remain unchanged. Thus we take the tunneling Hamiltonian to only depend on the
relative deformation of the two layers. The tunneling Hamiltonian can further
be constrained by the relevant symmetries near ${\bf K}$. A number of those
symmetry-related constraints are system dependent, however, the discrete
lattice translation symmetry ${T(\mathcal{D}({\bf r}+\alpha))=T(\mathcal{D})}$
is a common feature. Thus,

$T(\mathcal{D})=\sum_{{\bf G}}T_{{\bf G}}e^{i{\bf G}\cdot\mathcal{D}}\;,$ (21)

where ${\bf G}$ are the reciprocal lattice vectors of the undeformed
monolayer. Finally the Hamiltonian of the moiré bilayer can be represented as

$H_{\text{moir\'{e}}}=H_{\text{BL}}+H_{\text{tun}}\equiv
H_{0}+H[\mathcal{D},\bm{\nabla}\mathcal{D}],$ (22)

where $H_{0}$ is the Hamiltonian of the undeformed bilayer, which we take to
be non-polar, i.e. a commensurate bilayer with AA stacking. Thus, the local
polarization can be characterized by the evolution of the total polarization
as the deformation is turned on adiabatically:

${\bf P}(\mathcal{D}({\bf r}))=\int^{\mathcal{D}({\bf
r})}_{0}\frac{\partial{\bf P}}{\partial\mathcal{D}}d\mathcal{D}\;,$ (23)

where,

$\frac{\partial{\bf
P}}{\partial\mathcal{D}}=\frac{-2ief}{(2\pi)^{2}}\sum^{\text{occ}}_{n}\oint_{\text{mBZ}}\braket{\partial_{\mathcal{D}}u^{n}_{\mathcal{D},{\bf
G}}({\bf k})}{\partial_{{\bf k}}u^{n}_{\mathcal{D},{\bf G}}({\bf k})}\;,$ (24)

denoting mBZ as the moiré BZ, a specific case of the more general scBZ. Here,
$u^{n}_{\mathcal{D},{\bf G}}({\bf k})$ are the Bloch wavefunctions, obtained
by solving the continuum model Hamiltonian in Eq. (22).

###  II.4 Wannier functions and gauge invariance

When considering polarization in a crystal, it is natural to work in terms of
localized states such as Wannier functions Marzari  _et al._ (2012):

$\begin{split}\ket{w_{n,{\bf
R}}}&=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}e^{-i{\bf k}\cdot{\bf
R}}\ket{\psi_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\\\ \ket{\psi_{n,{\bf
k}}}&=\sum_{{\bf R}}e^{i{\bf k}\cdot{\bf R}}\ket{w_{n,{\bf R}}}\end{split}\;,$
(25)

which are the Fourier transforms of the Bloch states. We have one for each
band $n$ and each lattice vector ${\bf R}$. In the case of a supercell, ${\bf
R}$ represents a supercell vector. The Wannier functions are orthonormal, have
translational invariance, and are exponentially localized for a system with
topologically trivial electronic bands. In seminal works by King-Smith and
Vanderbilt Vanderbilt and King-Smith (1993); King-Smith and Vanderbilt (1993),
it was shown that the Wannier centers ${\bar{{\bf
w}}_{n}\equiv\braket{w_{n,0}}{{\bf r}}{w_{n,0}}}$, the expectation values of
the position operator in the Wannier basis, can be identified as the integral
of the Berry phases, with units of length (see Appendix A):

$\bar{{\bf w}}_{n}\equiv\bra{w_{n,0}}{\bf
r}\ket{w_{n,0}}=\frac{i\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}\braket{u_{n,{\bf
k}}}{\nabla_{{\bf k}}u_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\;.$ (26)

A well-known property of the Wannier functions is that their centers are
invariant, modulo a lattice vector, with respect to single-band gauge
transformations of the Bloch states:

$\ket{u_{n,{\bf k}}}\to\ket{u^{\prime}_{n,{\bf k}}}=e^{-i\beta_{n}({\bf
k})}\ket{u_{n,{\bf k}}}\;,$ (27)

where $\beta({\bf k})=\beta({\bf k}+{\bf G})+{\bf G}\cdot{\bf R}$, and ${\bf
G}$ is a reciprocal lattice vector. For a so-called small gauge transformation
defined by ${\bf R}=0$, each Wannier center is invariant: $\bar{{\bf
w}}^{\prime}_{n}=\bar{{\bf w}}_{n}$, whereas for a large transformation with
${\bf R}\neq 0$: $\bar{{\bf w}}^{\prime}_{n}=\bar{{\bf w}}_{n}+{\bf R}$, which
contributes a quantum of polarization to the total polarization, but does not
affect the physical observables.

All such transformations provide automorphisms of isolated bands. However,
when there are crossings between bands, the identity of single bands within
the band subspaces is lost, which can lead to potential problems with the
smoothness of integrands when computing Wannier centers Vanderbilt (2018). In
this case, on isolating an occupied band subspace from a Bloch bundle, more
general multi-band gauge transformations apply

$\ket{u_{n,{\bf k}}}\to\ket{u^{\prime}_{n,{\bf
k}}}=\sum^{\text{occ}}_{m}U_{nm}({\bf k})\ket{u_{m,{\bf k}}}\;,$ (28)

with the single-band gauge transformations ${U_{nm}({\bf
k})=\delta_{n,m}e^{-i\beta_{n}({\bf k})}}$ constituting only a subset of the
non-Abelian matrix transformations. Under such transformations, the trace of
the matrix-valued non-Abelian Berry connection,

${\bf A}_{nm}({\bf k})=i\braket{u_{n,{\bf k}}}{\nabla_{{\bf k}}u_{m,{\bf
k}}}\;,$ (29)

is preserved, rather than the individual components. Therefore in general, the
gauge invariant quantity is not the individual Wannier centers, but the sum of
Wannier centers of occupied bands Gresch  _et al._ (2017):

$\sum^{\text{occ}}_{n}\bar{{\bf
w}}^{\prime}_{n}=\sum^{\text{occ}}_{n}\bar{{\bf w}}_{n}\;.$ (30)

This gauge freedom is typically used to change the representation of the
states to obtain maximally localized Wannier functions Marzari  _et al._
(2012), using the Wannier90 code Pizzi  _et al._ (2020), for example. The
Wannier centers can be calculated in configuration space using first-
principles calculations and Wannier90:

$\bar{w}_{n}({\bf x})=\braket{w_{n,0}({\bf x})}{{\bf r}}{w_{n,0}({\bf x})}\;.$
(31)

As mentioned previously, the Wannier centers have been used to estimate the
local polarization in real space using Wu  _et al._ (2006):

${\bf P}({\bf r}_{j})=-\frac{ef}{\Omega_{0}}\sum^{\text{occ}}_{\bar{{\bf
w}}_{n}\in{\bf r}_{j}}\bar{{\bf w}}_{n}\;,$ (32)

where the local polarization in cell ${\bf r}_{j}$ is related to the sum over
Wannier centers in cell ${\bf r}_{j}$. However, Eq. (32) is not gauge
invariant in general, since individual sums of Wannier centers are not gauge
invariant in the multi-band case. Going from a unit cell to a supercell
shrinks the BZ and introduces significant band folding, in which case it might
not always be possible to disentangle the Wannier functions into isolated unit
cells.

We propose that the correct way to define the local polarization in terms of
Wannier centers is to calculate the changes of all the Wannier centers of
occupied bands in a supercell with respect to local perturbations:

${\bf P}({\bf
r}_{j})=-\frac{ef}{\Omega_{0}}\sum_{n}^{\text{occ}}\int_{0}^{{\bf x}({\bf
r}_{j})}\partial_{x^{\prime}_{\kappa,\alpha}}\bar{{\bf
w}}_{n}\mathop{}\\!\mathrm{d}x^{\prime}_{\kappa,\alpha}\;,$ (33)

where $\kappa\in{\bf r}_{j}$. Manifestly, this shows that the local
polarization is defined mod $e{\bf R}/\Omega_{0}$ under large gauge
transformations of the Bloch bands over the supercell BZ. As with Eq. (11),
Eq. (33) calculates the change in a global property of the system, which is
well-defined, with respect to local perturbations. If the perturbation in a
given cell only affects the Wannier centers attributed to that cell, then Eq.
(33) reduces to Eq. (32).

The expression for local polarization given by Eq. (33) is exact, but
evaluation of the Wannier centers and their variation in a supercell is
computationally and technically demanding. In practice, Wannier functions are
typically obtained in configuration space, which works well in many
circumstances; for example, tight-binding models of twisted bilayers
parametrized using Wannier functions have been shown to give accurate
descriptions of the electronic bands for a range of supercell sizes Carr  _et
al._ (2019a, b).

Nonetheless, it is interesting to assess how large a discrepancy there can be
between the local polarization obtained from configuration space and real
space Wannier centers. Consider a unit cell, repeated to form a supercell, but
with an identical configuration in each cell. Now, we switch on an additional
supercell potential correction, $\lambda\Delta V_{\text{sl}}({\bf x}({\bf
r}))$, to the potential imposed by configuration space, where
$\lambda=0\rightarrow 1$ parameterizes the switching on/off of the correction.
For $\lambda=0$, the Wannier functions should be the same in each unit cell by
translational invariance, and are thus equivalent to the Wannier functions at
a given point in configuration space: ${\ket{w^{\text{cs}}_{n,{\bf
R}}}\equiv\ket{w^{\lambda=0}_{n,{\bf R}}}}$. For $\lambda=1$, the potential in
each unit cell is different, and the Wannier functions are no longer
equivalent by translational invariance. We denote these as the supercell
Wannier functions: ${\ket{w^{\text{sl}}_{n,{\bf
R}}}\equiv\ket{w^{\lambda=1}_{n,{\bf R}}}}$, i.e. those obtained directly in
real space. In both $\ket{w^{\text{cs}}_{n,{\bf R}}}$ and
$\ket{w^{\text{sl}}_{n,{\bf R}}}$, ${\bf R}$ is a supercell vector. The
Wannier functions are given by ${w^{\lambda}_{n,{\bf R}}({\bf r}-{\bf
R})=\braket{{\bf r}}{w^{\lambda}_{n,{\bf R}}}}$ in the position
representation, and for any $\lambda$, and ${\bf R}=0$ we have

$\bar{{\bf w}}^{\lambda}_{n}=\bra{w^{\lambda}_{n,0}}{\bf
r}\ket{w^{\lambda}_{n,0}}\;,$ (34)

which can be used to track the evolution of the Wannier centres after
switching on the supercell potential $\Delta V_{\text{sl}}({\bf x}({\bf r}))$.
The cumulative change on applying the supercell potential correction can be
written as

$\ket{w^{\text{sl}}_{n,0}}=\ket{w^{\text{cs}}_{n,0}}+\ket{\Delta w_{n,0}}\;,$
(35)

where

$\ket{\Delta
w_{n,0}}=\int^{1}_{0}\ket{\partial_{\lambda}w_{n,0}}\mathop{}\\!\mathrm{d}\lambda\;.$
(36)

The change in the Wannier functions can be determined from the corresponding
changes in the Bloch states:

$\ket{\partial_{\lambda}w_{n,0}}=\frac{\Omega_{\text{sc}}}{(2\pi)^{3}}\int_{\text{scBZ}}e^{i{\bf
k}\cdot{\bf r}}\ket{\partial_{\lambda}u_{n{\bf k}}}\mathop{}\\!\mathrm{d}{\bf
k}\;,$ (37)

where the changes in the Bloch states are determined by the Sternheimer
equation Sternheimer (1954); Gonze (1995); Baroni  _et al._ (2001),

$\ket{\partial_{\lambda}u_{n{\bf k}}}=-A_{n}({\bf k})\ket{u_{n{\bf
k}}}+\sum_{m\neq n}\frac{\ket{u_{m{\bf k}}}\bra{u_{m{\bf
k}}}}{E_{n}-E_{m}}\left(\partial_{\lambda}H\right)\ket{u_{n{\bf k}}}\;,$ (38)

which can be evaluated using DFPT. The difference between the local
polarization obtained in real space and configuration space is determined by
the error in the Wannier centers:

$\begin{split}\Delta\bar{{\bf w}}_{n}=\bar{{\bf w}}^{\text{sl}}_{n}-\bar{{\bf
w}}^{\text{cs}}_{n}=\bra{\Delta w_{n,0}}{\bf
r}\ket{w^{\text{cs}}_{n,0}}+\text{c.c.}+O(\Delta\bar{{\bf
w}}_{n}^{2})\end{split}\;,$ (39)

to first order in $\Delta\bar{{\bf w}}_{n}$. The correction to local
polarization is then given by

$\begin{split}{\bf P}^{\text{err}}({\bf r}_{j})&={\bf P}^{\text{sl}}({\bf
r}_{j})-{\bf P}^{\text{cs}}({\bf r}_{j})\\\
&=-\frac{ef}{\Omega_{0}}\sum_{n}^{\text{occ}}\int_{0}^{{\bf x}({\bf
r}_{j})}\partial_{x^{\prime}_{\kappa,\alpha}}\Delta\bar{{\bf
w}}_{n}\mathop{}\\!\mathrm{d}x^{\prime}_{\kappa,\alpha}\end{split}\;,$ (40)

where ${\bf P}^{\text{sl/cs}}({\bf r}_{j})$ represent the local polarization
calculated using Wannier centers in real space / configuration space:

${\bf P}^{\text{sl/cs}}({\bf
r}_{j})=-\frac{ef}{\Omega_{0}}\sum_{n}^{\text{occ}}\int_{0}^{{\bf x}({\bf
r}_{j})}\partial_{x^{\prime}_{\kappa,\alpha}}\bar{{\bf
w}}^{\text{sl/cs}}_{n}\mathop{}\\!\mathrm{d}x^{\prime}_{\kappa,\alpha}.$ (41)

##  III Results

###  III.1 First-principles calculations

Figure 3:  (a) Orbital character of the 8 valence bands and the bottom of the
conduction band in bilayer hBN (AA stacking). (b) electronic band structure
obtained from first-principles calculations using abinit (black) and Wannier
interpolation using Wannier90 (red). (c) Out-of-plane and (d) in-plane
polarization in fractional coordinates along the configuration space diagonal,
calculated from the Berry phases (red), integrating the Born effective charges
(blue) and the Wannier centers (black). The positions high symmetry stackings
AA ($x=0$), AB ($x=\frac{1}{3}$), DW ($x=\frac{1}{2}$) and BA
($x=\frac{2}{3}$) are indicated by the ticks. (e-f): Maximally localized
Wannier functions obtained in bilayer hBN: (e) an $sp$ orbital aligned with
the B-N bond (three per layer) and (f) an $sp$ orbital normal to the layers
(one per layer).

While we propose that the most correct way to calculate local polarization is
from the Wannier centers / Born effective charges in real space, this is a
very computationally heavy and technically demanding task for large
supercells, which deserves its own dedicated study and is beyond the scope of
this work. As a stepping stone, we first show that our proposed expressions
for local polarization, Eqs. (11) and (33), give the correct expression in
configuration space, i.e. they are in agreement with the result obtained from
Berry phases. As an example, we consider 3R-stacked bilayer hBN, following the
methodology in Refs. Bennett and Remez (2022); Bennett (2020); Bennett  _et
al._ (2023).

First-principles density functional theory (DFT) calculations were performed
using the abinit Gonze  _et al._ (2009) code, using psml García  _et al._
(2018) norm-conserving pseudopotentials Hamann (2013), obtained from Pseudo-
Dojo Van Setten  _et al._ (2018). abinit employs a plane wave basis set, which
was determined using a kinetic energy cutoff of $1000$ eV, and a Monkhorst-
Pack ${\bf k}$-point grid Monkhorst and Pack (1976) of $12\times 12\times 1$
was used. Calculations were converged until the relative changes in the total
energy were less than $10^{-10}$ Ha. The revPBE exchange-correlation
functional was used Zhang and Yang (1998), and the vdw-DFT-D3(BJ) Becke and
Johnson (2006) correction was used to treat the long-range interactions
between the layers.

The top hBN layer was translated along the unit cell diagonal over the bottom
layer, which was held fixed. At each point a geometry relaxation was performed
to obtain the equilibrium layer separation, while keeping the in-plane lattice
vectors fixed. The out-of-plane and in-plane polarization were then obtained
by calculating the Berry phases of the Bloch states. The out-of-plane
polarization $P_{\perp}$ was found to be odd with respect to the relative
stacking, and the in-plane polarization, $\Delta{\bf P}_{\parallel}$, was
found to be even, as reported in Ref. Bennett  _et al._ (2023). At each point
along the unit cell diagonal, the Born effective charges were calculated by
calculating the change in the Hamiltonian and Bloch states in response to
phonon and electric field perturbations, using the DFPT routines in abinit. At
each point in configuration space, the Wannier functions were also calculated
using Wannier90, which interfaces with abinit. Maximally localized Wannier
functions Marzari  _et al._ (2012) were obtained by using the gauge freedom to
minimize the spread. The orbital characters of the 8 valence bands and the
lowest lying conduction band, see Fig. 3, were used to determine the initial
projections onto atomic orbitals. The lowest valence bands are comprised of
in-plane $sp$-like states, with three in each layer. The highest valence bands
and lowest conduction bands have $p_{z}$ character, originating from the N and
B atoms, respectively. The spread was minimized until the relative change was
less than $10^{-10}$ Å2 and the bands obtained from Wannier interpolation were
in agreement with the bands obtained from abinit, see Fig. 3 (b). The bands
are well-reproduced with three in-plane $sp$-like states (Fig. 3 (e)) and one
$sp_{z}$-like state (Fig. 3 (f)) for each layer.

The out-of-plane and in-plane polarization of bilayer hBN in configuration
space are shown in Fig. 3 (c) and (d), respectively. The red data show the
polarization obtained from Berry phases, i.e. calculating the polarization
from the electronic bands of a commensurate bilayer plus a relative
translation between the layers. The blue data show the polarization obtained
by integrating the Born effective charges along the unit cell diagonal, i.e.
Eq. (11). As shown in Ref. Bennett  _et al._ (2023), Eq. (2) is not
appropriate in layered ferroelectrics, but integrating the Born effective
charges yields a polarization identical to the one obtained from the Berry
phases. The black data show the polarization obtained by measuring the sum of
Wannier centers in configuration space, which is in exact agreement with the
other two methods. In configuration space, the 8 Wannier functions describe
all of the occupied bands, so the sum of Wannier centers for each
configuration is always gauge invariant, and Eq. (33) reduces to Eq. (32).

###  III.2 Effective model

Figure 4: (a) Sketch of the two-cosine model for a one-dimensional supercell:
a periodic chain of atoms with spacing $a$. The potential from the ions
$V_{\text{c}}(x)$ is shown in gray. The supercell potential
$V_{\text{sc}}(x)$, with period $Na$, is shown in black. (b) Sketch of the
system with a modified potential due to a depolarizing perturbation from a
probe. The action of the perturbation brings the cell $j$ to a non-polar
configuration and the response of the Wannier centers of occupied bands can be
determined to measure the local polarization.

The calculation of polarization from first-principles calculations in
configuration space are a valid approximation to the local polarization in
real space for large supercells. As a proof of concept, we show using a one-
dimensional effective model that our proposed definition of local
polarization, Eq. (33), can be calculated directly in real space without
relying on this approximation. We consider a model inspired by the Aubry-André
(AA) model Aubry and André (1980) which consists of a one-dimensional chain of
atoms with an atomic spacing $a$ and a supercell potential with period $Na$,
where ${N\in\mathbb{N}}$ is not necessarily large, see Fig. 4 (a). The model
is described by the one-body Hamiltonian

$H=-\frac{1}{2}\frac{d^{2}}{dx^{2}}+V_{\text{c}}(x)+V_{\text{sc}}(x)\;,$ (42)

in terms of one-dimensional position $x$, where the first term represents the
kinetic energy of an electron, ${V_{\text{c}}(x)=2V_{0}\cos(Gx)}$ is the ionic
potential with $G=\frac{2\pi}{a}$, and ${V_{\text{sc}}(x)=2\lambda
V_{0}\cos(G_{N}x)}$ is the supercell potential, which can be interpreted as a
contribution due to the local displacement of the cores, where
${G_{N}=\frac{G}{N}}$. As in the previous section, we introduce a parameter
${\lambda}$ which determines the relative strength between the core and
supercell potentials and can be used to switch on the supercell potential.
Contrary to the previous case, $\lambda$ can be larger than $\lambda=1$ if the
supercell potential is stronger than the core potential.

For weak core and superlattice potentials, we can perturbatively obtain the
eigenstates of (42) using the nearly-free electron gas approach via nearly-
degenerate perturbation theory, using a basis of plane-wave states $\ket{k}$
with energies ${\epsilon_{k}=\frac{1}{2}k^{2}}$ from the free-particle
problem. The potential mixes different plane wave states
$\\{\ket{k+mG_{N}}\\}$, for ${\\{m\in\mathbb{Z}:|m|\leq N\\}}$, yielding a set
of secular equations. Evaluating the corresponding matrix elements of the
effective Hamiltonian,

$\begin{split}\bra{k+mG_{N}}H\ket{k+nG_{N}}&=\epsilon_{k+mG_{N}}\delta_{m,n}\\\
&+V_{0}(\delta_{m,n+N}+\delta_{m,n-N})\\\ &+\lambda
V_{0}(\delta_{m,n+1}+\delta_{m,n-1})\end{split}\;,$ (43)

we obtain the following matrix equation:

$\leavevmode\resizebox{184.28981pt}{}{$\begin{pmatrix}\epsilon_{k}&\lambda
V_{0}&\lambda V_{0}&\cdots&V_{0}&V_{0}\\\ \lambda
V^{*}_{0}&\epsilon_{k+G_{N}}&0&\cdots&0&0\\\ \lambda
V^{*}_{0}&0&\epsilon_{k-G_{N}}&\ddots&0&0\\\
\vdots&\vdots&\ddots&\ddots&\vdots&\vdots\\\
V^{*}_{0}&0&0&\cdots&\epsilon_{k+G}&0\\\
V^{*}_{0}&0&0&\cdots&0&\epsilon_{k-G}\par\end{pmatrix}\begin{pmatrix}c_{k}\\\
c_{k+G_{N}}\\\ c_{k-G_{N}}\\\ \vdots\\\ c_{k+G}\\\
c_{k-G}\end{pmatrix}=E_{k}\begin{pmatrix}c_{k}\\\ c_{k+G_{N}}\\\
c_{k-G_{N}}\\\ \vdots\\\ c_{k+G}\\\ c_{k-G}\end{pmatrix}$}\;.$ (44)

The approximate (unnormalized) eigenstates are given by

$\ket{\psi}=\sum^{N}_{m={-N}}c_{k+mG_{N}}\ket{k+mG_{N}}\;,$ (45)

with position representation

$\psi(x)=e^{ikx}\sum^{N}_{m={-N}}c_{k+mG_{N}}e^{imG_{N}x}\;.$ (46)

The approximate eigenstates are Bloch states, ${\psi(x)=e^{ikx}u_{n{\bf
k}}(x)}$, with cell-periodic parts which are periodic over a supercell period:
${u_{n{\bf k}}(x)=u_{n{\bf k}}(x+Na)}$. Solving the secular equations yields
$2N+1$ coefficients $c_{k+mG_{N}}$ for each $k$, and a set of energy bands
$E_{k}$ which are sensitive to the parameters $V_{0}$ and $\lambda$. The basis
states corresponding to larger reciprocal lattice vectors can also be included
to increase the accuracy of the approximate solution to the model. Eq. (44)
was solved numerically for $N=5$, and the resulting bands are shown in Fig. 5
(a). For $\lambda=0$, the folding of the bands into the scBZ introduces many
band crossings. Switching on the supercell potential with $\lambda\neq 0$
opens up several gaps where the bands cross.

Having found the approximate bands and eigenstates, we can proceed to obtain
the Wannier functions $w^{\lambda}_{n,X}$ for a given $\lambda$, where $n$ is
the band index and $X$ is a lattice vector. For $\lambda=0$, the Wannier
functions are given by

$w^{0}_{n,X}(x-ja)\propto\int_{\text{BZ}}e^{ik(x-ja)}\left(c_{k}+c_{k+G}e^{iGx}\right)\mathop{}\\!\mathrm{d}k\;,$
(47)

where $X$ is a unit cell vector and the integral is over the BZ,
${\left[-\frac{\pi}{a},\frac{\pi}{a}\right]}$, because the periodicity reduces
to $a$ in the absence of a supercell potential. For $\lambda\neq 0$, the
Wannier functions are given by

$w^{\lambda}_{n,X}(x-jNa)\propto\int_{\text{scBZ}}e^{ik(x-jNa)}\sum^{N}_{m={-N}}c_{k+mG_{N}}e^{imG_{N}x}\mathop{}\\!\mathrm{d}k\;,$
(48)

where $X$ is a supercell vector and the integral is over the scBZ,
${[-\frac{\pi}{Na},\frac{\pi}{Na}]}$. The Wannier functions were obtained and
their centers are plotted in Fig. 5 (b) for two different gauges. The first
gauge projected the Wannier centers onto the atomic sites, as indicated by the
black lines. For the second gauge, there is a displacement between the Wannier
centers and the atomic sites in each cell, as indicated by the red lines.
However, the sum of Wannier centers is the same in each case, and both sets of
Wannier functions can be used as a localized basis. This illustrates the
arbitrariness of describing local polarization in a supercell using partial
sums over Wannier centers with Eq. (32), as the partial sums are not gauge
invariant.

Figure 5: (a) $N=5$ two-cosine model bands for $\lambda=0$ (dashed) and for
$\lambda=1$ (solid), with $V_{0}=-1$. As the supercell potential strength is
switched on, additional gaps open. (b) Wannier centers for the five occupied
supercell bands, using two different multi-band gauges. The first gauge
results in atom-centered Wannier functions (black), and the second gauge
results in Wannier functions with cell-dependent off-centering (red). The
Wannier centers are indicated by the vertical dashed lines in both cases. The
total sum of Wannier centers is the same for both gauges. As the Wannier
functions in arbitrary gauges are not well-localized, their representations
around different Wannier centers are purely illustrative.

In the spirit of Eq. (33), we introduce an additional local potential to Eq.
(42) in order to calculate the displacements of all the Wannier centers in
response to local perturbations. This can be achieved by switching a cell $j$
to a non-polar configuration using a depolarizing potential
$V_{\text{dep},j}(x)$:

$\leavevmode\resizebox{188.62363pt}{}{$V_{\text{dep},j}(x)=\left[\theta\left(x-ja+\frac{a}{2}\right)-\theta\left(x-ja-\frac{a}{2}\right)\right]\left[V_{\text{sc}}(x-ja)-V_{\text{sc}}(x)\right]$}\;,$
(49)

where $\theta(x)$ is the Heaviside step function. Such a potential could be
achieved using a tip to locally probe the system, see Fig. 4 (a). The
supercell potential in cell $j$ is removed and replaced with the potential of
the non-polar configuration. The matrix elements of $V_{\text{dep},j}$ are
given by

$\leavevmode\resizebox{188.62363pt}{}{$\begin{split}(V_{\text{dep},j})_{nm}&=\Lambda_{j}\int^{ja+\frac{a}{2}}_{ja-\frac{a}{2}}e^{iG_{N}(m-n)x}\sin(G_{N}(x-ja/2))\mathop{}\\!\mathrm{d}x\\\
\Lambda_{j}&=\frac{4\lambda V_{0}}{Na}\sin{\left(\frac{\pi
j}{N}\right)}\end{split}$}\;.$ (50)

Adding these matrix elements to the secular equations, Eq. (44), the change in
the Wannier centers of the occupied bands in response to local perturbations
in each cell can be obtained, see Fig. 6 (a), which allows the local
polarization to be calculated, see Fig. 6 (b). We note that the polarization
is odd, and reminiscent of the out-of-plane polarization in bilayer hBN along
the unit cell diagonal, see Fig. 3 (c). By the symmetry of the potential, the
local polarization is odd about $x=0$, and sums to zero.

Figure 6: (a) Shift of Wannier centers in response to local perturbations
$V_{\text{dep},j}$ in each cell, for $\lambda=2$ and $V_{0}=-2$. (b) Local
polarization calculated from the displacement of all Wannier centers in
response to the perturbation potential $V_{\text{dep},j}$.

##  IV Discussions and Conclusions

In this work, we propose a definition of the local polarization in a crystal
supercell, obtained either by integrating the Born effective charges, Eq.
(11), or the changes in the Wannier centers of all occupied bands with respect
to local perturbations, Eq. (33). In this way, the polarization can be
calculated directly in real space, and is gauge invariant. We show that Eqs.
(11) and (33) are in agreement with the polarization obtained from Berry
phases when evaluated in configuration space. Although we propose that it is
more correct to calculate the polarization directly in real space,
calculations involving DFPT or Wannierization with large supercells deserve
their own dedicated study and we do not pursue this here. However, this is an
important subject of future research: in many systems where topologically non-
trivial polarization textures are observed, it should be checked that the
shape and topology are the same when calculating the polarization in
configuration space and in real space. For example, in twisted bilayers, the
configuration space approximation is only valid for small angles, and while it
is clear there should be a polar–nonpolar transition somewhere between
${0^{\circ}<\theta<60^{\circ}}$, it is not immediately clear where this
transition occurs, or what the polarization field looks like beyond the small
angle regimes Bennett (2020). Such large-scale calculations may not be
unrealistically expensive: fortunately, larger twist angles result in smaller
supercells, meaning the length scales where the configuration space
approximation breaks down naturally contain fewer atoms. Furthermore, because
the Born effective charge tensor is a mixed derivative of the free energy, it
can be obtained using either type of perturbation (phonon or electric field).
Thus, instead of $3N$ phonon perturbations, where $N$ is very large, the
effective charges can be obtained with only 3 electric field perturbations.
This is known as the ‘interchange theorem’ Gonze and Lee (1997).

We also demonstrate as a proof of concept that the local polarization can be
calculated directly in real space in a one-dimensional effective model. The
polarization is calculated in a small supercell, for which the configuration
space approximation is not valid, by calculating the change in the Wannier
centers of all occupied bands in response to a depolarizing potential from a
local probe. The calculated local polarization is consistent with the form of
the potential experienced by electrons.

A better definition and understanding of local polarization is essential for
considering topological polarization. As mentioned previously, because
polarization is topological in the real space sense, the topological
properties are solely determined by the geometry of the supercell and the
underlying crystal symmetry of the lattice. Having a correct definition of the
local polarization such as Eq. (11) shows that topological polarization can
indeed be calculated from the bands of a supercell, and that this real space
topology can be described at the electronic level. This is important when
considering the relation between polarization and band topology. Recently, it
was predicted that ferroelectric switching via van der Waals sliding can lead
to a change in the QAH conductance in a topological insulator Liang  _et al._
(2023). Another recent theory proposes that that for a material in a moiré
potential, which results in electronic bands with nonzero Chern numbers, the
Chern numbers can be altered by applying an electric field Ghorashi  _et al._
(2022). In this example, the moiré superlattice potential was achieved using a
non-polar material. Substituting a ferroelectric material may lead to Chern
bands which can be switched at zero field. Although in both Refs. Liang  _et
al._ (2023) and Ghorashi  _et al._ (2022) the polarization was uniform, they
both suggest that polarization can influence the topological properties of a
system. There may be additional phenomena which are unique to topologically
non-trivial polarization textures, although currently no such phenomena have
been considered or proposed.

Similar to the effect of polarization on the band topology, one may consider
the effect of the band topology on the local polarization textures. When
considering local polarization, we can define a higher-dimensional space
spanned by the BZ and configuration space, which for the case of a bilayer
would be a 4-dimensional torus. In this higher dimensional space, the local
polarization resembles the phase-space Berry curvature (see Eq. 11). Thus the
local polarization can be linked to the evolution of the phase-space Berry
curvature and Chern numbers as a function of relative configurations. However,
the local polarization is defined as an integral from the nonpolar reference
state which does not form a closed loop in configuration space, and thus the
Chern theorem does not apply unless considering a translation by a unit cell.
For example, the most important shift in a bilayer is the one associated to
the ferroelectric switching of the polarization, achieved by a relative
sliding of a third of a unit cell diagonal: ${{\bf x}=\frac{1}{3}({\bf
a}_{1}+{\bf a}_{2})}$. The change in polarization associated to this switching
process does not form a closed loop in configuration space, and thus cannot
result in a nonzero Chern number. Furthermore, in bilayer hBN the electronic
bands are topologically trivial for every point in configuration space, and
sliding one layer over the other does not change this. Because polarization
alone does not break time-reversal symmetry, a polarization texture cannot by
itself lead to any non-trivial Chern band topology. For this to occur, time-
reversal symmetry must be broken by some other means. In order to see the
interplay between topological polarization and band topology, we must consider
a system which is already topologically non-trivial, but is also polar, such
as those considered in Refs. Liang  _et al._ (2023) and Ghorashi  _et al._
(2022). In this case, the description of local polarization developed here
would need to be generalized to the case of Chern insulators, by tracking the
evolution of hybrid Wannier centers Coh and Vanderbilt (2009); Song  _et al._
(2021). Here, bulk-boundary relations Qi and Zhang (2011); Hasan and Kane
(2010); Slager  _et al._ (2015); Rhim  _et al._ (2018); Hatsugai (1993) also
pose an intriguing question for future pursuits.

##  V Acknowledgements

D.B. and E. K. acknowledge funding from the US Army Research Office (ARO) MURI
project under grant No. W911NF-21-0147 and by the National Science Foundation
DMREF program under Award No. DMR-1922172. W. J. J. acknowledges funding from
the Rod Smallwood Studentship at Trinity College, Cambridge. R. J. S. and G.
C. acknowledge funding from a New Investigator Award, EPSRC grant
EP/W00187X/1. R. J. S. also acknowledges funding from Trinity College,
University of Cambridge.

##  VI Appendix A: Relation between Berry phases and Wannier centers

For illustrative purposes, we briefly review the derivation for the important
relation between the Berry phases and the Wannier centers Vanderbilt (2018).
First, we apply the position operator to a Wannier state $\ket{w_{n,{\bf
R}}}$:

$\begin{split}{\bf r}\ket{w_{n,{\bf
R}}}&=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}{\bf r}e^{-i{\bf k}\cdot{\bf
R}}\ket{\psi_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\\\
&=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}{\bf r}e^{i{\bf k}\cdot({\bf
r}-{\bf R})}\ket{u_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\end{split}\;,$
(51)

where we write the Bloch state as $\ket{\psi_{n,{\bf k}}}=e^{i{\bf k}\cdot{\bf
r}}\ket{u_{n,{\bf k}}}$. Next we rewrite the product of the position operator
and the exponential:

${\bf r}e^{i{\bf k}\cdot({\bf r}-{\bf R})}=({\bf R}-i\nabla_{{\bf k}})e^{i{\bf
k}\cdot({\bf r}-{\bf R})}\;,$ (52)

giving

${\bf r}\ket{w_{n,{\bf R}}}=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}({\bf
R}-i\nabla_{{\bf k}})e^{i{\bf k}\cdot({\bf r}-{\bf R})}\ket{u_{n,{\bf
k}}}\mathop{}\\!\mathrm{d}{\bf k}\;.$ (53)

For the term proportional to $\nabla_{{\bf k}}$, we perform integration by
parts. The boundary term vanishes because the integral is around a closed loop
in the BZ and $\ket{u_{n,{\bf k}}}$ has translational invariance. This gives

${\bf r}\ket{w_{n,{\bf R}}}=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}e^{i{\bf
k}\cdot({\bf r}-{\bf R})}({\bf R}+i\nabla_{{\bf k}})\ket{u_{n,{\bf
k}}}\mathop{}\\!\mathrm{d}{\bf k}\;.$ (54)

Now we multiply by

$\bra{w_{n,0}}=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}e^{-i{\bf
k}^{\prime}\cdot{\bf r}}\bra{u_{n,{\bf k}^{\prime}}}\mathop{}\\!\mathrm{d}{\bf
k}^{\prime}$ (55)

which gives

$\bra{w_{n0}}{\bf r}\ket{w_{n{\bf
R}}}=\frac{\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}e^{-i{\bf k}\cdot{\bf
R}}\bra{u_{n{\bf k}}}{\bf R}+i\nabla_{{\bf k}}\ket{u_{n{\bf
k}}}\mathop{}\\!\mathrm{d}{\bf k}\;.$ (56)

The term proportional to ${\bf R}$ gives ${\bf R}\delta_{0{\bf R}}$, which
vanishes on setting ${\bf R}\to 0$, as we get the expected expression for the
Wannier centers:

$\bar{{\bf w}}_{n}\equiv\bra{w_{n,0}}{\bf
r}\ket{w_{n,0}}=\frac{i\Omega}{(2\pi)^{3}}\oint_{\text{BZ}}\braket{u_{n,{\bf
k}}}{\nabla_{{\bf k}}u_{n,{\bf k}}}\mathop{}\\!\mathrm{d}{\bf k}\;.$ (57)

## References

  * Mitsui and Furuichi (1953) Toshito Mitsui and Jiro Furuichi, “Domain structure of Rochelle salt and KH2PO4,” Phys. Rev. 90, 193 (1953).
  * Kopal  _et al._ (1997) A. Kopal, T. Bahnik, and J. Fousek, “Domain formation in thin ferroelectric films: the role of depolarization energy,” Ferroelectrics 202, 267–274 (1997).
  * Kittel (1946) Charles Kittel, “Theory of the structure of ferromagnetic domains in films and small particles,” Phys. Rev. 70, 965–971 (1946).
  * Junquera and Ghosez (2003) Javier Junquera and Philippe Ghosez, “Critical thickness for ferroelectricity in perovskite ultrathin films,” Nature 422, 506–509 (2003).
  * Bennett  _et al._ (2020) Daniel Bennett, Maitane Muñoz Basagoiti, and Emilio Artacho, “Electrostatics and domains in ferroelectric superlattices,” R. Soc. Open Sci. 7, 201270 (2020).
  * Luk’yanchuk  _et al._ (2009) Igor A. Luk’yanchuk, Laurent Lahoche, and Anaïs Sené, “Universal properties of ferroelectric domains,” Phys. Rev. Lett. 102, 147601 (2009).
  * Ji  _et al._ (2019) Dianxiang Ji, Songhua Cai, Tula R Paudel, Haoying Sun, Chunchen Zhang, Lu Han, Yifan Wei, Yipeng Zang, Min Gu, Yi Zhang, _et al._ , “Freestanding crystalline oxide perovskites down to the monolayer limit,” Nature 570, 87–90 (2019).
  * Gómez-Ortiz  _et al._ (2023) Fernando Gómez-Ortiz, Hugo Aramberri, Juan M López, Pablo García-Fernández, Jorge Íñiguez, and Javier Junquera, “Kittel law and domain formation mechanism in PbTiO3/SrTiO3 superlattices,” arXiv:2303.01755 (2023).
  * Callori  _et al._ (2012) S. J. Callori, J. Gabel, D. Su, J. Sinsheimer, M. V. Fernandez-Serra, and M. Dawber, “Ferroelectric PbTiO3/SrRuO3 superlattices with broken inversion symmetry,” Phys. Rev. Lett. 109, 067601 (2012).
  * Zhang  _et al._ (2017) Qingteng Zhang, Eric M. Dufresne, Pice Chen, Joonkyu Park, Margaret P. Cosgriff, Mohammed Yusuf, Yongqi Dong, Dillon D. Fong, Hua Zhou, Zhonghou Cai, Ross J. Harder, Sara J. Callori, Matthew Dawber, Paul G. Evans, and Alec R. Sandy, “Thermal fluctuations of ferroelectric nanodomains in a ferroelectric-dielectric PbTiO3/SrTiO3 superlattice,” Phys. Rev. Lett. 118, 097601 (2017).
  * Dawber (2017) Matthew Dawber, “Balancing polar vortices and stripes,” Nat. Mater. 16, 971–972 (2017).
  * Park  _et al._ (2018) Joonkyu Park, John Mangeri, Qingteng Zhang, M Humed Yusuf, Anastasios Pateras, Matthew Dawber, Martin V Holt, Olle G Heinonen, Serge Nakhmanson, and Paul G Evans, “Domain alignment within ferroelectric/dielectric PbTiO3/SrTiO3 superlattice nanostructures,” Nanoscale 10, 3262–3271 (2018).
  * Susarla  _et al._ (2021) Sandhya Susarla, Pablo García-Fernández, Colin Ophus, Sujit Das, Pablo Aguado-Puente, Margaret McCarter, Peter Ercius, Lane W Martin, Ramamoorthy Ramesh, and Javier Junquera, “Atomic scale crystal field mapping of polar vortices in oxide superlattices,” Nat. Commun. 12, 6273 (2021).
  * Urban  _et al._ (2002) Jeffrey J Urban, Wan Soo Yun, Qian Gu, and Hongkun Park, “Synthesis of single-crystalline perovskite nanorods composed of barium titanate and strontium titanate,” J. Am. Chem. Soc. 124, 1186–1187 (2002).
  * Yun  _et al._ (2002) Wan Soo Yun, Jeffrey J Urban, Qian Gu, and Hongkun Park, “Ferroelectric properties of individual barium titanate nanowires investigated by scanned probe microscopy,” Nano Lett. 2, 447–450 (2002).
  * Luo  _et al._ (2003) Yun Luo, Izabela Szafraniak, Nikolai D Zakharov, Valanoor Nagarajan, Martin Steinhart, Ralf B Wehrspohn, Joachim H Wendorff, Ramamoorthy Ramesh, and Marin Alexe, “Nanoshell tubes of ferroelectric lead zirconate titanate and barium titanate,” Appl. Phys. Lett. 83, 440–442 (2003).
  * Morrison  _et al._ (2003) Finlay D Morrison, Laura Ramsay, and James F Scott, “High aspect ratio piezoelectric strontium–bismuth–tantalate nanotubes,” J. Phys.: Condens. Matter 15, L527 (2003).
  * Mao  _et al._ (2003) Yuanbing Mao, Sarbajit Banerjee, and Stanislaus S Wong, “Hydrothermal synthesis of perovskite nanotubes,” Chem. Commun. , 408–409 (2003).
  * Chu  _et al._ (2004) Ming-Wen Chu, Izabela Szafraniak, Roland Scholz, Catalin Harnagea, Dietrich Hesse, Marin Alexe, and Ulrich Gösele, “Impact of misfit dislocations on the polarization instability of epitaxial nanostructured ferroelectric perovskites,” Nat. Mater. 3, 87–90 (2004).
  * Shin  _et al._ (2005) Hyung-Joon Shin, Je Hyuk Choi, Hee Jun Yang, Young Dae Park, Young Kuk, and Chi-Jung Kang, “Patterning of ferroelectric nanodot arrays using a silicon nitride shadow mask,” Appl. Phys. Lett. 87, 113114 (2005).
  * Fu and Bellaiche (2003) Huaxiang Fu and L Bellaiche, “Ferroelectricity in barium titanate quantum dots and wires,” Phys. Rev. Lett. 91, 257601 (2003).
  * Naumov  _et al._ (2004) Ivan I Naumov, L Bellaiche, and Huaxiang Fu, “Unusual phase transitions in ferroelectric nanodisks and nanorods,” Nature 432, 737–740 (2004).
  * Geneste  _et al._ (2006) Gregory Geneste, Eric Bousquet, Javier Junquera, and Philippe Ghosez, “Finite-size effects in BaTiO3 nanowires,” Appl. Phys. Lett. 88, 112906 (2006).
  * Morozovska  _et al._ (2006) Anna N Morozovska, Eugene A Eliseev, and Maya D Glinchuk, “Ferroelectricity enhancement in confined nanorods: Direct variational method,” Phys. Rev. B 73, 214106 (2006).
  * Hong  _et al._ (2010) Jiawang Hong, G Catalan, DN Fang, Emilio Artacho, and JF Scott, “Topology of the polarization field in ferroelectric nanowires from first principles,” Phys. Rev. B 81, 172101 (2010).
  * Junquera  _et al._ (2023) Javier Junquera, Yousra Nahas, Sergei Prokhorenko, Laurent Bellaiche, Jorge Íñiguez, Darrell G. Schlom, Long-Qing Chen, Sayeef Salahuddin, David A. Muller, Lane W. Martin, and R. Ramesh, “Topological phases in polar oxide nanostructures,” Rev. Mod. Phys. 95, 025001 (2023).
  * Nahas  _et al._ (2015) Y Nahas, S Prokhorenko, L Louis, Z Gui, Igor Kornev, and Laurent Bellaiche, “Discovery of stable skyrmionic state in ferroelectric nanocomposites,” Nat. Commun. 6, 1–6 (2015).
  * Pereira Gonçalves  _et al._ (2019) Mauro António Pereira Gonçalves, Carlos Escorihuela-Sayalero, Pablo Garca-Fernández, Javier Junquera, and Jorge Íñiguez, “Theoretical guidelines to create and tune electric skyrmion bubbles,” Sci. Adv. 5, eaau7023 (2019).
  * Das  _et al._ (2019) S Das, YL Tang, Z Hong, MAP Gonçalves, MR McCarter, C Klewe, KX Nguyen, F Gómez-Ortiz, P Shafer, E Arenholz, _et al._ , “Observation of room-temperature polar skyrmions,” Nature 568, 368–372 (2019).
  * Han  _et al._ (2022) Lu Han, Christopher Addiego, Sergei Prokhorenko, Meiyu Wang, Hanyu Fu, Yousra Nahas, Xingxu Yan, Songhua Cai, Tianqi Wei, Yanhan Fang, _et al._ , “High-density switchable skyrmion-like polar nanodomains integrated on silicon,” Nature 603, 63–67 (2022).
  * Shao  _et al._ (2023) Yu-Tsun Shao, Sujit Das, Zijian Hong, Ruijuan Xu, Swathi Chandrika, Fernando Gómez-Ortiz, Pablo García-Fernández, Long-Qing Chen, Harold Y Hwang, Javier Junquera, _et al._ , “Emergent chirality in a polar meron to skyrmion phase transition,” Nat. Commun. 14, 1355 (2023).
  * Wang  _et al._ (2020) YJ Wang, YP Feng, YL Zhu, YL Tang, LX Yang, MJ Zou, WR Geng, MJ Han, XW Guo, B Wu, _et al._ , “Polar meron lattice in strained oxide ferroelectrics,” Nat. Mater. 19, 881–886 (2020).
  * Li and Wu (2017) Lei Li and Menghao Wu, “Binary compound bilayer and multilayer with vertical polarizations: Two-dimensional ferroelectrics, multiferroics, and nanogenerators,” ACS Nano 11, 6382–6388 (2017).
  * Bennett and Remez (2022) Daniel Bennett and Benjamin Remez, “On electrically tunable stacking domains and ferroelectricity in moiré superlattices,” npj 2D Mater. Appl. 6, 1–6 (2022).
  * Bennett (2020) Daniel Bennett, “Theory of polar domains in moiré heterostructures,” Phys. Rev. B 105, 235445 (2020).
  * Zheng  _et al._ (2020) Zhiren Zheng, Qiong Ma, Zhen Bi, Sergio de la Barrera, Ming-Hao Liu, Nannan Mao, Yang Zhang, Natasha Kiper, Kenji Watanabe, Takashi Taniguchi, _et al._ , “Unconventional ferroelectricity in moiré heterostructures,” Nature 588, 71 (2020).
  * Yasuda  _et al._ (2021) Kenji Yasuda, Xirui Wang, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, “Stacking-engineered ferroelectricity in bilayer boron nitride,” Science 372, 1458 (2021).
  * Stern  _et al._ (2021) M. Vizner Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. Ben Shalom, “Interfacial ferroelectricity by van der Waals sliding,” Science 372, 1462 (2021).
  * Bennett  _et al._ (2023) Daniel Bennett, Gaurav Chaudhary, Robert-Jan Slager, Eric Bousquet, and Philippe Ghosez, “Polar meron-antimeron networks in strained and twisted bilayers,” Nat. Commun. 14, 1629 (2023).
  * Shen  _et al._ (2022) Jiaying Shen, Zhengang Dong, MingQun Qi, Yang Zhang, Chao Zhu, Zhenping Wu, and Danfeng Li, “Observation of moiré patterns in twisted stacks of bilayer perovskite oxide nanomembranes with various lattice symmetries,” ACS Appl. Mater. Interfaces 14, 50386–50392 (2022).
  * Sanchez-Santolino  _et al._ (2023) Gabriel Sanchez-Santolino, Victor Rouco, Sergio Puebla, Hugo Aramberri, Victor Zamora, Fabian A Cuellar, Carmen Munuera, Federico Mompean, Mar Garcia-Hernandez, Aandres Castellanos-Gomez, _et al._ , “A 2D ferroelectric vortex lattice in twisted BaTiO3 freestanding layers,” arXiv:2301.04438 (2023).
  * Artyukhov  _et al._ (2020) Vasilii I Artyukhov, Sunny Gupta, Alex Kutana, and Boris I Yakobson, “Flexoelectricity and charge separation in carbon nanotubes,” Nano Lett. 20, 3240–3246 (2020).
  * Springolo  _et al._ (2021) Matteo Springolo, Miquel Royo, and Massimiliano Stengel, “Direct and converse flexoelectricity in two-dimensional materials,” Phys. Rev. Lett. 127, 216801 (2021).
  * Bennett (2021) Daniel Bennett, “Flexoelectric-like radial polarization of single-walled nanotubes from first-principles,” Electron. Struct. 3, 015001 (2021).
  * Vanderbilt and King-Smith (1993) David Vanderbilt and RD King-Smith, “Electric polarization as a bulk quantity and its relation to surface charge,” Phys. Rev. B 48, 4442 (1993).
  * King-Smith and Vanderbilt (1993) RD King-Smith and David Vanderbilt, “Theory of polarization of crystalline solids,” Phys. Rev. B 47, 1651 (1993).
  * Resta (1994) Raffaele Resta, “Macroscopic polarization in crystalline dielectrics: the geometric phase approach,” Rev. Mod. Phys. 66, 899–915 (1994).
  * Nakagawa  _et al._ (2020) Masaya Nakagawa, Robert-Jan Slager, Sho Higashikawa, and Takashi Oka, “Wannier representation of floquet topological states,” Phys. Rev. B 101, 075108 (2020).
  * Vanderbilt (2018) David Vanderbilt,  _Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators_ (Cambridge University Press, 2018).
  * Van Mechelen  _et al._ (2022) Todd Van Mechelen, Sathwik Bharadwaj, Zubin Jacob, and Robert-Jan Slager, “Optical ${N}$-insulators: Topological obstructions to optical Wannier functions in the atomistic susceptibility tensor,” Phys. Rev. Res. 4, 023011 (2022).
  * Wu  _et al._ (2006) Xifan Wu, Oswaldo Diéguez, Karin M Rabe, and David Vanderbilt, “Wannier-based definition of layer polarizations in perovskite superlattices,” Phys. Rev. Lett. 97, 107602 (2006).
  * Meyer and Vanderbilt (2002) B Meyer and David Vanderbilt, “Ab initio study of ferroelectric domain walls in PbTiO3,” Phys. Rev. B 65, 104111 (2002).
  * Stengel  _et al._ (2011) Massimiliano Stengel, Pablo Aguado-Puente, Nicola A Spaldin, and Javier Junquera, “Band alignment at metal/ferroelectric interfaces: Insights and artifacts from first principles,” Phys. Rev. B 83, 235112 (2011).
  * Qi and Zhang (2011) Xiao-Liang Qi and Shou-Cheng Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, 1057–1110 (2011).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium: Topological insulators,” Rev. Mod. Phys. 82, 3045–3067 (2010).
  * Fu (2011) Liang Fu, “Topological crystalline insulators,” Phys. Rev. Lett. 106, 106802 (2011).
  * Slager  _et al._ (2013) Robert-Jan Slager, Andrej Mesaros, Vladimir Juričić, and Jan Zaanen, “The space group classification of topological band-insulators,” Nat. Phys. 9, 98–102 (2013).
  * Ünal  _et al._ (2020) F. Nur Ünal, Adrien Bouhon, and Robert-Jan Slager, “Topological Euler class as a dynamical observable in optical lattices,” Phys. Rev. Lett. 125, 053601 (2020).
  * Po  _et al._ (2017) Hoi Chun Po, Ashvin Vishwanath, and Haruki Watanabe, “Symmetry-based indicators of band topology in the 230 space groups,” Nat. Commun. 8, 50 (2017).
  * Bradlyn  _et al._ (2017) Barry Bradlyn, Luis Elcoro, Jennifer Cano, Maia G Vergniory, Zhijun Wang, Claudia Felser, Mois I Aroyo, and B Andrei Bernevig, “Topological quantum chemistry,” Nature 547, 298–305 (2017).
  * Kruthoff  _et al._ (2017) Jorrit Kruthoff, Jan De Boer, Jasper Van Wezel, Charles L Kane, and Robert-Jan Slager, “Topological classification of crystalline insulators through band structure combinatorics,” Phys. Rev. X 7, 041069 (2017).
  * Shiozaki and Sato (2014) Ken Shiozaki and Masatoshi Sato, “Topology of crystalline insulators and superconductors,” Phys. Rev. B 90, 165114 (2014).
  * Bouhon  _et al._ (2020a) Adrien Bouhon, QuanSheng Wu, Robert-Jan Slager, Hongming Weng, Oleg V Yazyev, and Tomas Bzdusek, “Non-Abelian reciprocal braiding of Weyl points and its manifestation in ZrTe,” Nat. Phys. 16, 1137–1143 (2020a).
  * Bouhon  _et al._ (2020b) Adrien Bouhon, Tomas Bzdusek, and Robert-Jan Slager, “Geometric approach to fragile topology beyond symmetry indicators,” Phys. Rev. B 102, 115135 (2020b).
  * Song  _et al._ (2018) Zhida Song, Tiantian Zhang, Zhong Fang, and Chen Fang, “Quantitative mappings between symmetry and topology in solids,” Nat. Commun. 9, 3530 (2018).
  * Coh and Vanderbilt (2009) Sinisa Coh and David Vanderbilt, “Electric polarization in a Chern insulator,” Phys. Rev. Lett. 102, 107603 (2009).
  * Song  _et al._ (2021) Xue-Yang Song, Yin-Chen He, Ashvin Vishwanath, and Chong Wang, “Electric polarization as a nonquantized topological response and boundary Luttinger theorem,” Phys. Rev. Research 3, 023011 (2021).
  * Alexandradinata and Bernevig (2016) A Alexandradinata and B. A. Bernevig, “Berry-phase description of topological crystalline insulators,” Phys. Rev. B 93, 205104 (2016).
  * Bouhon  _et al._ (2019) Adrien Bouhon, Annica M Black-Schaffer, and Robert-Jan Slager, “Wilson loop approach to fragile topology of split elementary band representations and topological crystalline insulators with time-reversal symmetry,” Phys. Rev. B 100, 195135 (2019).
  * Liang  _et al._ (2023) Yan Liang, Fulu Zheng, Thomas Frauenheim, and Pei Zhao, “Ferroelectric antiferromagnetic quantum anomalous Hall insulator in two dimensional van der Waals materials,” arXiv:2302.05091 (2023).
  * Szabó and Schneider (2018) Attila Szabó and Ulrich Schneider, “Non-power-law universality in one-dimensional quasicrystals,” Phys. Rev. B 98, 134201 (2018).
  * Borgnia  _et al._ (2022) Dan S. Borgnia, Ashvin Vishwanath, and Robert-Jan Slager, “Rational approximations of quasiperiodicity via projected green’s functions,” Phys. Rev. B 106, 054204 (2022).
  * Borgnia and Slager (2023) Dan S. Borgnia and Robert-Jan Slager, “Localization as a consequence of quasiperiodic bulk-bulk correspondence,” Phys. Rev. B 107, 085111 (2023).
  * Aubry and André (1980) Serge Aubry and Gilles André, “Analyticity breaking and anderson localization in incommensurate lattices,” Ann. Israel Phys. Soc 3, 18 (1980).
  * Bistritzer and MacDonald (2011) Rafi Bistritzer and Allan H. MacDonald, “Moiré bands in twisted double-layer graphene,” Proceedings of the National Academy of Sciences 108, 12233–12237 (2011), https://www.pnas.org/doi/pdf/10.1073/pnas.1108174108 .
  * Balents (2019) L. Balents, “General continuum model for twisted bilayer graphene and arbitrary smooth deformations,” SciPost Phys. 7, 48 (2019).
  * Carr  _et al._ (2017) Stephen Carr, Daniel Massatt, Shiang Fang, Paul Cazeaux, Mitchell Luskin, and Efthimios Kaxiras, “Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle,” Phys. Rev. B 95, 075420 (2017).
  * Cazeaux  _et al._ (2017) Paul Cazeaux, Mitchell Luskin, and Ellad B Tadmor, “Analysis of rippling in incommensurate one-dimensional coupled chains,” Multiscale Model. Simul. 15, 56–73 (2017).
  * Massatt  _et al._ (2017) Daniel Massatt, Mitchell Luskin, and Christoph Ortner, “Electronic density of states for incommensurate layers,” Multiscale Model. Simul. 15, 476–499 (2017).
  * Carr  _et al._ (2018) Stephen Carr, Daniel Massatt, Steven B Torrisi, Paul Cazeaux, Mitchell Luskin, and Efthimios Kaxiras, “Relaxation and domain formation in incommensurate two-dimensional heterostructures,” Phys. Rev. B 98, 224102 (2018).
  * Cohen (1992) Ronald E Cohen, “Origin of ferroelectricity in perovskite oxides,” Nature 358, 136–138 (1992).
  * Sluka  _et al._ (2013) Tomas Sluka, Alexander K Tagantsev, Petr Bednyakov, and Nava Setter, “Free-electron gas at charged domain walls in insulating BaTiO3,” Nat. Commun. 4, 1808 (2013).
  * Bednyakov  _et al._ (2015) Petr S Bednyakov, Tomas Sluka, Alexander K Tagantsev, Dragan Damjanovic, and Nava Setter, “Formation of charged ferroelectric domain walls with controlled periodicity,” Sci. Rep. 5, 15819 (2015).
  * Engelke  _et al._ (2023) Rebecca Engelke, Hyobin Yoo, Stephen Carr, Kevin Xu, Paul Cazeaux, Richard Allen, Andres Mier Valdivia, Mitchell Luskin, Efthimios Kaxiras, Minhyong Kim, Jung Hoon Han, and Philip Kim, “Topological nature of dislocation networks in two-dimensional moiré materials,” Phys. Rev. B 107, 125413 (2023).
  * Resta (1998) Raffaele Resta, “Quantum-mechanical position operator in extended systems,” Phys. Rev. Lett. 80, 1800–1803 (1998).
  * Wu  _et al._ (2005) Xifan Wu, David Vanderbilt, and D. R. Hamann, “Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory,” Phys. Rev. B 72, 035105 (2005).
  * Gonze and Lee (1997) Xavier Gonze and Changyol Lee, “Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory,” Phys. Rev. B 55, 10355 (1997).
  * Vanderbilt (2000) D Vanderbilt, “Berry-phase theory of proper piezoelectric response,” J. Phys. Chem. Solids 61, 147–151 (2000).
  * Bennett  _et al._ (2022) Daniel Bennett, Daniel Tanner, Philippe Ghosez, Pierre-Eymeric Janolin, and Eric Bousquet, “Generalized relation between electromechanical responses at fixed voltage and fixed electric field,” Phys. Rev. B 106, 174105 (2022).
  * Ghosez  _et al._ (1998) Ph Ghosez, J-P Michenaud, and Xavier Gonze, “Dynamical atomic charges: The case of ABO3 compounds,” Phys. Rev. B 58, 6224 (1998).
  * Ghosez and Gonze (2000) Philippe Ghosez and Xavier Gonze, “Band-by-band decompositions of the Born effective charges,” J. Phys.: Condens. Matter 12, 9179 (2000).
  * Sai  _et al._ (2002) Na Sai, Karin M Rabe, and David Vanderbilt, “Theory of structural response to macroscopic electric fields in ferroelectric systems,” Phys. Rev. B 66, 104108 (2002).
  * Marzari  _et al._ (2012) Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and David Vanderbilt, “Maximally localized Wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012).
  * Gresch  _et al._ (2017) Dominik Gresch, Gabriel Autès, Oleg V. Yazyev, Matthias Troyer, David Vanderbilt, B. Andrei Bernevig, and Alexey A. Soluyanov, “Z2pack: Numerical implementation of hybrid Wannier centers for identifying topological materials,” Phys. Rev. B 95, 075146 (2017).
  * Pizzi  _et al._ (2020) Giovanni Pizzi, Valerio Vitale, Ryotaro Arita, Stefan Blügel, Frank Freimuth, Guillaume Géranton, Marco Gibertini, Dominik Gresch, Charles Johnson, Takashi Koretsune, _et al._ , “Wannier90 as a community code: new features and applications,” J. Phys.: Condens. Matter 32, 165902 (2020).
  * Carr  _et al._ (2019a) Stephen Carr, Shiang Fang, Ziyan Zhu, and Efthimios Kaxiras, “Exact continuum model for low-energy electronic states of twisted bilayer graphene,” Phys. Rev. Res. 1, 013001 (2019a).
  * Carr  _et al._ (2019b) Stephen Carr, Shiang Fang, Hoi Chun Po, Ashvin Vishwanath, and Efthimios Kaxiras, “Derivation of Wannier orbitals and minimal-basis tight-binding Hamiltonians for twisted bilayer graphene: First-principles approach,” Phys. Rev. Res. 1, 033072 (2019b).
  * Sternheimer (1954) R. M. Sternheimer, “Electronic polarizabilities of ions from the Hartree-Fock wave functions,” Phys. Rev. 96, 951–968 (1954).
  * Gonze (1995) Xavier Gonze, “Adiabatic density-functional perturbation theory,” Phys. Rev. A 52, 1096–1114 (1995).
  * Baroni  _et al._ (2001) Stefano Baroni, Stefano de Gironcoli, Andrea Dal Corso, and Paolo Giannozzi, “Phonons and related crystal properties from density-functional perturbation theory,” Rev. Mod. Phys. 73, 515–562 (2001).
  * Gonze  _et al._ (2009) Xavier Gonze, Bernard Amadon, P-M Anglade, J-M Beuken, François Bottin, Paul Boulanger, Fabien Bruneval, Damien Caliste, Razvan Caracas, Michel Côté, _et al._ , “Abinit: First-principles approach to material and nanosystem properties,” Comput. Phys. Commun. 180, 2582 (2009).
  * García  _et al._ (2018) Alberto García, Matthieu J Verstraete, Yann Pouillon, and Javier Junquera, “The PSML format and library for norm-conserving pseudopotential data curation and interoperability,” Comput. Phys. Commun. 227, 51 (2018).
  * Hamann (2013) DR Hamann, “Optimized norm-conserving Vanderbilt pseudopotentials,” Phys. Rev. B 88, 085117 (2013).
  * Van Setten  _et al._ (2018) MJ Van Setten, Matteo Giantomassi, Eric Bousquet, Matthieu J Verstraete, Don R Hamann, Xavier Gonze, and G-M Rignanese, “The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table,” Comput. Phys. Commun. 226, 39 (2018).
  * Monkhorst and Pack (1976) Hendrik J Monkhorst and James D Pack, “Special points for Brillouin-zone integrations,” Phys. Rev. B 13, 5188 (1976).
  * Zhang and Yang (1998) Yingkai Zhang and Weitao Yang, “Comment on “Generalized gradient approximation made simple”,” Phys. Rev. Lett. 80, 890 (1998).
  * Becke and Johnson (2006) Axel D Becke and Erin R Johnson, “A simple effective potential for exchange,” J. Chem. Phys. 124, 221101 (2006).
  * Ghorashi  _et al._ (2022) Sayed Ali Akbar Ghorashi, Aaron Dunbrack, Jiacheng Sun, Xu Du, and Jennifer Cano, “Topological and stacked flat bands in bilayer graphene with a superlattice potential,” arXiv:2206.13501 (2022).
  * Slager  _et al._ (2015) Robert-Jan Slager, Louk Rademaker, Jan Zaanen, and Leon Balents, “Impurity-bound states and Green’s function zeros as local signatures of topology,” Phys. Rev. B 92, 085126 (2015).
  * Rhim  _et al._ (2018) Jun-Won Rhim, Jens H. Bardarson, and Robert-Jan Slager, “Unified bulk-boundary correspondence for band insulators,” Phys. Rev. B 97, 115143 (2018).
  * Hatsugai (1993) Yasuhiro Hatsugai, “Chern number and edge states in the integer quantum Hall effect,” Phys. Rev. Lett. 71, 3697–3700 (1993).