Theory of polarization textures in crystal supercells

Recently, topologically nontrivial 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\'e 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.


INTRODUCTION
The formation of complex polar structures such as ferroelectric domains [1,2], analogous to ferromagnetic domains [3], 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 nonpolar material or vacuum lead to polar discontinuities which if not screened will lead to depolarizing fields that suppress ferroelectricity [4].Ferroelectric materials can form polydomain structures to mitigate depolarization effects, such as 180 • stripe domains, which for not too thin films can be described using a Landau-like theory [5].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 [6]; such domains have been observed in ferroelectric materials down to the monolayer limit [7].In some cases, even more complex structures involving polar vortices may form [8], for example in ferroelectric/paraelectric (FE/PE) superlattices [9][10][11][12][13] such as 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 [11], depending on the strength of the coupling between the ferroelectric layers [6], which can be determined by the ratio of the layers and the dielectric permittivities [5].
Ferroelectric materials have been fabricated in many different geometries, from 2D thin films and FE/PE superlattices to 1D nanowires [14,15] and nanotubes [16][17][18], and even 0D quantum dots [19,20].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 [21][22][23][24].Soon after these complex polarization textures were discovered, they were proposed to be topologically non-trivial [25,26]; skyrmion-like polarization structures were identified, for example in BaTiO 3 (BTO) nanowires embedded in a matrix of STO [27].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 [28], resulting in a non-trivial winding.Ferroelectric skyrmions have been experimentally observed in PTO/STO superlattices [29,30], and a skyrmion → meron transition with strain has recently been observed [31].Polar merons have also been observed in PTO under epitaxial strain from a SmScO 3 substrate [32].
In recent years, a new type of ferroelectricity in layered van der Waals materials has been proposed [33][34][35] and experimentally observed [36,37].In layered systems such as 3R-stacked hexagonal boron nitride (hBN) or MoS 2 , an outof-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 [38]), 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 [37] via the growing and shrinking of the MPDs in response to an applied field [34,35].Recently, we discovered that the MPDs also have an in-plane component of polarization [39], resulting in possible real space non- 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 = − 1  2 ) and merons (Q = + 1  2 ), respectively.trivial topological textures.For hBN and similar materials the topological charge Q in each stacking domain, where the local polarization P(r) is normalized, integrates to ± 1 2 , meaning the MPDs are polar merons and antimerons (see Fig. 1).Moiré superlattices have also recently been fabricated with perovskites [40], which for FE/PE interfaces such as BTO/STO also result in polar vortices [41].
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 [28], 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 [42][43][44], 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 [39].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 inplane 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 [45][46][47][48][49][50], 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 [51], 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 [52,53].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 [39], 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 [54,55], which has seen the development of various analytical diagnoses [56][57][58][59][60][61][62][63][64][65], 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) [66,67] (equivalent to calculating Wilson loops [68,69]).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 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 [70].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 [70], 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 [34,35,39], 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 [51][52][53], 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 [71][72][73][74]) and 2D (Bistritzer-MacDonald [75,76]) 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.

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 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 P(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.

Configuration Space
The mapping between real space and configuration space is used to estimate local properties in periodic supercells and incommensurate systems [77][78][79][80].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 ABO 3 perovskite, the polarization arises from the offcentering between the B cations and the 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 [81].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 x(r j ) = ∆r Ti j , where r j is the position of unit cell j, and r Ti j 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 x(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: where Ω 0 is the unit cell volume, contrary to the total supercell volume Ω sc .The first sum is over the atoms κ in unit cell r j , and the second is over Cartesian directions α.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 [53], 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 [8].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 [82,83].
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 [39].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 θ between the layers, see Fig. 2 (c), the configuration space mapping is given by [80] x modulo any lattice vectors, where r is the real space position and R θ 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 x ≈ θ 0 −1 1 0 r between the layers, see Fig. 2 (d).Similarly, for a small homogeneous strain η, the equivalent mapping is x = η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 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 [34,35,39].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: where summation over repeated indices is assumed, and again only the atoms κ ∈ 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 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 [84].Eq. ( 4) yields a polarization identical to the one obtained from Berry phases [39].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.

Defining local polarization
From the modern theory of polarization, the total polarization P tot in a crystal is given by [45,46,49] where |u n,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 where a α are the lattice vectors, Ω is the system cell volume.
is the Berry connection of band n, and is the Berry phase of band n in direction α, defined up to a factor of 2π, where b α are the reciprocal lattice vectors.For disordered systems, where 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 [85].
It is well-known that the absolute polarization in a crystal is not well-defined.Only changes in polarization are welldefined, 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) [86]: the dielectric response of a system is related to the derivative of the total polarization with respect to electric field [87], and electromechanical properties can be measured by calculating the derivatives of the polarization with respect to strain (piezoelectricity [88]), or strain and electric field (electrostriction [89]).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 [87,90]: where x κ,α is the real space (phonon) displacement of atom κ in direction α, F κ,α is the force on atom κ in direction α, and E β is an electric field in direction β .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 [39].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 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 [91,92] for all atoms κ in cell r j , the local polarization in each unit cell is given by 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 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 x.The limits of the integral are the nonpolar reference state 0 and the local configuration in each unit cell 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 ∆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 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.

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 [75].
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 [76].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 K of the un-deformed monolayer.The electron field at K is given by 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. [39], although generalization to other models and dimensions is straightforward.The monolayer Hamiltonian is given by where m is the mass gap, v is the Dirac-velocity and the τ Pauli-matrices act on some internal degrees, which for hBN are the two sublattices.Summation is assumed, with µ = 1, 2. The effect of small strain and twist can be captured by a local deformation D(r) field as where r is real space position in a "laboratory frame", x is the position in the monolayer and D(r) is a deformation field as a result of strain and twist.The deformation field is assumed to be locally small, i.e. ∂ µ D 1.It is crucial to define the "small" deformation D as a function of the variable r instead of the variable x, since when considering the multilayers, x are associated with the individual layers correspond to very different locations r in the real space.Thus one cannot define a "small local" deformation field in the variable 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 The ψ field is correspondingly modified as and the integral measure is modified as The continuum Hamiltonian of the decoupled bilayers can be obtained when each layer experiences an independent deformation field D l , where l = t,b is the layer index: where we have kept only the terms linear in the deformation field.For example, if the two layers are twisted rigidly by angles ±θ /2, the deformation field is given by 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 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 K.A number of those symmetry-related constraints are system dependent, however, the discrete lattice translation symmetry T (D(r + α)) = T (D) is a common feature.Thus, where G are the reciprocal lattice vectors of the undeformed monolayer.Finally the Hamiltonian of the moiré bilayer can be represented as 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: where, denoting mBZ as the moiré BZ, a specific case of the more general scBZ.Here, u n D,G (k) are the Bloch wavefunctions, obtained by solving the continuum model Hamiltonian in Eq. ( 22).

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 [93]: which are the Fourier transforms of the Bloch states.We have one for each band n and each lattice vector R.In the case of a supercell, 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 [45,46], it was shown that the Wannier centers wn ≡ w n,0 |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): 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: where and G is a reciprocal lattice vector.For a so-called small gauge transformation defined by R = 0, each Wannier center is invariant: w n = wn , whereas for a large transformation with R = 0: w n = wn + 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 [49].In this case, on isolating an occupied band subspace from a Bloch bundle, more general multi-band gauge transformations apply with the single-band gauge transformations U nm (k) = δ n,m e −iβ n (k) constituting only a subset of the non-Abelian matrix transformations.Under such transformations, the trace of the matrix-valued non-Abelian Berry connection, 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 [94]: This gauge freedom is typically used to change the representation of the states to obtain maximally localized Wannier functions [93], using the WANNIER90 code [95], for example.The Wannier centers can be calculated in configuration space using first-principles calculations and WANNIER90: As mentioned previously, the Wannier centers have been used to estimate the local polarization in real space using [51]: where the local polarization in cell r j is related to the sum over Wannier centers in cell 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: where κ ∈ r j .Manifestly, this shows that the local polarization is defined mod eR/Ω 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 [96,97].
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, λ ∆V sl (x(r)), to the potential imposed by configuration space, where λ = 0 → 1 parameterizes the switching on/off of the correction.For λ = 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: . For λ = 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: |w sl n,R ≡ |w λ =1 n,R , i.e. those obtained directly in real space.In both |w cs n,R and |w sl n,R , R is a supercell vector.The Wannier functions are given by w λ n,R (r − R) = r|w λ n,R in the position representation, and for any λ , and R = 0 we have which can be used to track the evolution of the Wannier centres after switching on the supercell potential ∆V sl (x(r)).The cumulative change on applying the supercell potential correction can be written as where The change in the Wannier functions can be determined from the corresponding changes in the Bloch states: where the changes in the Bloch states are determined by the Sternheimer equation [98][99][100], 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: to first order in ∆ wn .The correction to local polarization is then given by where P sl/cs (r j ) represent the local polarization calculated using Wannier centers in real space / configuration space:

First-principles calculations
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.[34,35,39].
First-principles density functional theory (DFT) calculations were performed using the ABINIT [101] code, using PSML [102] norm-conserving pseudopotentials [103], obtained from Pseudo-Dojo [104].ABINIT employs a plane wave basis set, which was determined using a kinetic energy cutoff of 1000 eV, and a Monkhorst-Pack k-point grid [105] of 12 × 12 × 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 [106], and the vdw-DFT-D3(BJ) [107] 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 ⊥ was found to be odd with respect to the relative stacking, and the in-plane polarization, ∆P , was found to be even, as reported in Ref. [39].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 [93] 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. [39], 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).

Effective model
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 [74] which consists of a one-dimensional chain of atoms with an atomic spacing a and a supercell potential with period Na, where N ∈ N is not necessarily large, see Fig. 4 (a).The model is described by the one-body Hamil- tonian in terms of one-dimensional position x, where the first term represents the kinetic energy of an electron, V c (x) = 2V 0 cos(Gx) is the ionic potential with G = 2π a , and V sc (x) = 2λ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 = G N .As in the previous section, we introduce a parameter λ 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, λ can be larger than λ = 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 |k with energies ε k = 1 2 k 2 from the free-particle problem.The potential mixes different plane wave states {|k + mG N }, for {m ∈ Z : |m| ≤ N}, yielding a set of secular equations.Evaluating the corresponding matrix elements of the effective Hamiltonian, we obtain the following matrix equation: The approximate (unnormalized) eigenstates are given by with position representation The approximate eigenstates are Bloch states, ψ(x) = e ikx u nk (x), with cell-periodic parts which are periodic over a supercell period: u nk (x) = u nk (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 λ .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 λ = 0, the folding of the bands into the scBZ introduces many band crossings.Switching on the supercell potential with λ = 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 λ n,X for a given λ , where n is the band index and X is a lattice vector.For λ = 0, the Wannier functions are given by where X is a unit cell vector and the integral is over the BZ, − π a , π a , because the periodicity reduces to a in the absence of a supercell potential.For λ = 0, the Wannier functions are given by where X is a supercell vector and the integral is over the scBZ, [− π Na , π 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.
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 dep, j (x): where θ (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 dep, j are given by 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.

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 • < θ < 60 • , it is not immediately clear where this transition occurs, or what the polarization field looks like beyond the small angle regimes [35].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.Fur-thermore, 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' [87].
We also demonstrate as a proof of concept that the local polarization can be calculated directly in real space in a onedimensional 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 [70].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 [108].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.[70] and [108] 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: x = 1 3 (a 1 + 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.[70] and [108].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 [66,67].Here, bulk-boundary relations [54,55,[109][110][111] also pose an intriguing question for future pursuits.

FIG. 2 .
FIG. 2. (a): Sketch of a 5 × 1 × 1 supercell of a titanium oxide perovskite 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 × 3 section of a bilayer with relative twist angle θ .(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.

∆P = P • 1 √ 2 (a 1 + a 2 )FIG. 3 .FIG. 4 .
FIG. 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) Outof-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 = 1 3 ), DW (x = 1 2 ) and BA (x = 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).

FIG. 5 .
FIG. 5. (a) N = 5 two-cosine model bands for λ = 0 (dashed) and for λ = 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 celldependent 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.

FIG. 6 .
FIG. 6.(a) Shift of Wannier centers in response to local perturbations V dep, j in each cell, for λ = 2 and V 0 = −2.(b) Local polarization calculated from the displacement of all Wannier centers in response to the perturbation potential V dep, j .