Optimally Localized Wannier Functions for 2D Chern Insulators

The construction of optimally localized Wannier functions (and Wannier functions in general) for a Chern insulator has been considered to be impossible owing to the fact that the second moment of such functions is generally infinite. In this manuscript, we propose a solution to this problem in the case of a single band. We accomplish this by drawing an analogy between the minimization of the variance and the minimization of the electrostatic energy of a periodic array of point charges in a smooth neutralizing background. In doing so, we obtain a natural regularization of the diverging variance and this leads to an analytical solution to the minimization problem. We demonstrate our results numerically for a particular model system. Furthermore, we show how the optimally localized Wannier functions provide a natural way of evaluating the electric polarization for a Chern insulator.


I. INTRODUCTION
The quantum mechanics of particles in a periodic potential forms a cornerstone of condensed matter physics.A standard approach in solid-state physics is to work with states that are eigenstates of the underlying Hamiltonian that have crystal momentum as a good quantum number [1].However, since such Bloch states are extended in real space, they can be ill-suited for understanding local phenomena like covalent bonding in a crystal.Therefore, the construction of Wannier functions [2], which are spatially localized states describing a given band or collection of bands, is often desirable.Roughly speaking, these Wannier states bridge the gap between solid-state physics and quantum chemistry, where a local description is preferred.
Since early days, quantum chemists have used so-called Localised Molecular Orbitals (LMO's) to study chemical bonding in molecular systems [3,4].Wannier functions have provided a powerful means to extend this to studying bonding and other structural properties such as defects in general solid state systems [5][6][7][8].Furthermore, these states are connected to the Berry phase theory of polarization, providing an alternative and illuminating viewpoint for studying things such as dielectric properties [9].They are also useful as a set of basis states in studying quantum ballistic transport via the Landauer formalism [10,11].More recently, Wannier functions have taken a central role in classifying topological band systems [12,13].
To construct Wannier functions for a single band, one integrates the Bloch functions over the Brillouin zone as [2] w nR (r) = 1 V BZ BZ dk e −ik•R ψ nk (r).
(Details of notation from this section are described more systematically in the following sections.)Such states are typically strongly localized in real space and possess desirable properties: they are orthogonal and the Wannier states for a given Bloch band are simply translations of each other by real space lattice vectors.The above definition has an important ambiguity, namely that the phases of the Bloch states are not specified.One can make the replacement ψ nk (r) → e iθ(k) ψ nk (r), where e iθ(k) has kspace periodicity, and obtain an equally good definition of the Wannier function.For the multi-band case, there is a larger ambiguity corresponding to a unitary matrix.
In general, changing the gauge in this fashion can have a non-trivial effect on the properties of the Wannier functions.Therefore, it is natural to seek a recipe to the fix the gauge such that we obtain the most desirable Wannier functions.One approach which is well-suited for conventional (i.e.non-topological) systems is to minimize the spatial variance, which directly reflects localization, of the Wannier function with respect to the gauge freedom.This approach, and the resulting Maximally Localized Wannier Functions (MLWFs) [14], has achieved wide-scale adoption in the electronic structure community and has numerous applications.
A major difficulty emerges [14] when attempting to construct Wannier states for topological insulators such as Chern insulators.Chern insulators are canonical models exhibiting topological bands [15,16], with the most famous being the Haldane model [17].The difficulty arises due to the presence of a topological obstruction preventing the construction of Bloch states that smoothly depend on k while having the correct periodicity.Indeed, if one could find such a smooth dependence, then a straightforward application of Stokes' theorem would reveal that the Chern number will always vanish.We remark that other types of topological obstructions do not necessarily cause issues.For example, it has been shown that globally smooth Bloch states can be constructed in the case of an Z 2 obstruction as long as one does not demand that the Wannier functions preserve time-reversal symmetry [18].
To deal with the difficulty in Chern insulators, one approach is to use different gauges to cover two or more patches of the Brillouin zone so that the Bloch states have smooth k-dependence in each patch.However, such Bloch states will not generally be smooth along lines where the patches meet.Another approach, which will be taken in this work, is to allow the Bloch states to have vortices in reciprocal space.This approach is sensible because we need to use a single choice of gauge for the Bloch states in (1), and we would like the Bloch states to be as smooth as possible.Generally, the sum of the winding numbers of all the vortices in the Brillouin zone will correspond to the Chern number of the band under consideration.For example, for a band with Chern number one, the simplest possibility is to allow for one singlyquantized vortex.Now, the gauge ambiguity allows one a choice not just of local variations of the phase but also of the positions of vortices, which can be moved by singular gauge transformations, which will be discussed in detail later in this article.Topology dictates the number of defects but not their location.This phenomenon is similar to the Hairy Ball theorem from algebraic topology [19].
The aim of this article is to construct optimally localized Wannier functions for Chern insulating bands.Besides the ambiguity of the vortex positions we just discussed, it is well-known that it is not possible to construct Wannier functions for such systems that decay exponentially at large distances [20,21].This is because we are effectively evaluating the Fourier transform of a function with vortex singularities.The importance of this fact has taken a central role in giving a defining relation for topological band systems.In particular, in [12], a system is defined to be topological if there is no Wannier representation composed of exponentially localized states respecting symmetries of the underlying lattice.In [13], it was argued conversely that states without such an atomic limit are either topological or exhibit so-called fragile topology.The theory of Wannier functions, originally a concept in electronic structure theory, has thus taken a prominent role in the newer field of topological states of matter.Although we cannot construct exponentially localized Wannier functions for a Chern insulator, many of the well-known applications of Wannier functions (e.g.Wannier interpolation and, as will be discussed in this article, the computation of polarization) do not rely on exponential decay.In this article, we construct optimal algebraically localized Wannier functions which are wellsuited for these applications.
Wannier functions for a Chern insulating band in 2D will generically decay asymptotically as 1/r 2 and this is an issue because it leads to a logarithmic divergence of the variance [21].Additionally, the location of the previously discussed Brillouin zone vortices, which depend on the gauge chosen, will influence quantities which are gauge invariant for the non-topological case such as the so-called Wannier center, which is the average position of a Wannier state.As an example, one can imagine nucleating a vortex and anti-vortex pair infinitesimally close together in the Brillouin zone and then transporting the vortex around the Brillouin zone torus until it returns to the anti-vortex and annihilates.Such a process can be achieved by a (singular) gauge transformation alone.This process will continuously shift the Wannier center by one lattice constant.The importance of a system-atic way of fixing this singular part of the gauge field is therefore clear.
We solve the problem of minimizing the variance over different vortex positions by identifying a dual picture in which the vortices can be thought of as periodic arrays of point charges in a smooth neutralising background charge distribution.In this setting, minimizing the variance corresponds to minimizing the electrostatic energy per unit cell.This approach provides a natural way to regularize the diverging variance.The source of its divergence is identified with the self-energy of the point charges, which can be neglected (as in ordinary electrostatics) from the total energy.Crucially, this self-energy does not depend on where the point charges are placed.
To test this theoretical development, we next turn to a concrete example of a Chern insulator lattice model.We construct Wannier functions where the non-singular portions of the gauges are optimised (i.e. are in the Coulomb gauge) but with vortex locations at various places in the Brillouin zone.We demonstrate that the localization of these Wannier functions are in quantitative agreement with the theoretical predictions, with the most localized being the one that minimizes the electrostatic energy in the dual picture.Finally, as a first application of the developed scheme of optimal gauge fixing, we demonstrate how it can be used to compute the polarization of a Chern insulator, once again illustrating on a particular model.Even though we have conducted our numerical tests on a particular model, we emphasize that the theoretical results presented in this article are applicable to any Chern insulator.

II. BACKGROUND AND PRELIMINARY CONCEPTS
In this section, we give some background to the development of localized Wannier functions and set notation.Readers familiar with this may wish to only skim this section.We begin by considering a periodic quantum system with Hamiltonian H, whose eigenstates are given by Bloch's theorem as ψ nk (r) = e ik•r u nk (r) where the u nk (r) possesses the real-space periodicity of the Hamiltonian.As usual, n labels the band index and k labels reciprocal space wave vectors.We use the normalization convention ⟨ψ nk |ψ mk ′ ⟩ = V BZ δ nm δ(k − k ′ ) for the Bloch states, where V BZ is the area (or volume in higher dimensions) of the Brillouin zone.
The corresponding Wannier functions [2] are then defined by (1) where R denotes a lattice vector.These states are orthonormal and span the same space as the corresponding Bloch states with which they are constructed.Furthermore, the Wannier functions for a given band are all translations of each other, which can be observed by applying position translation operators to the Wannier functions in the home unit cell.

A. Maximally Localized Wannier Functions
The primary motivation for Wannier functions comes from the fact they are generally very localized compared to the highly delocalized Bloch states.However, Wannier functions as defined by ( 1) are not unique.The Bloch states have an inherent gauge degree of freedom which usually does not have any physical significance.Yet, this choice of gauge has a non-trivial effect on the localization properties of Wannier functions, giving rise to different shapes and decay properties.
To choose between these different Wannier functions, one needs a concrete and systematic way of fixing the gauge of the Bloch states.The canonical way of achieving this for non-topological systems is to minimize the total spatial variance of the set of Wannier functions corresponding to a particular set of Bloch bands as originally described in [22].Since the Wannier functions for a given band are translations of the Wannier function in the home unit cell, it is sufficient to minimize the total variance of the Wannier functions in the home unit cell, which is given by the localization functional where n ranges over the set of bands of interest.
It is often more useful to express this variance in momentum space.For the rest of this article, we will focus on a single isolated band and so we drop band indices in what follows.Following work by Blount [23], we can write F in reciprocal space as [14,22] where A(k) = i⟨u k |∇ k u k ⟩ is the Berry connection of the band under consideration (see Appendix A for a detailed derivation).Moreover, we note that where the Einstein summation convention is used and α = x, y, z.
We now wish to minimize the above functional with respect to smooth gauge degrees of freedom of the Bloch states.To this end, consider the following family of Bloch states where we assume that θ(k) is a smooth and periodic function and |ũ k ⟩ is a fixed reference state.With respect to this reference state, the terms in (3) are Then, taking the functional derivative of F with respect to θ gives Thus the maximally localized state will satisfy the Coulomb gauge condition ∇ k • A = 0.This quite remarkable and perhaps under-appreciated result was first found by Blount [23].Using this, MLWFs can be constructed numerically via a steepest descent algorithm in such a way that one maintains the unitarity of the gauge transformation.This approach was formulated in [22] and is applicable to the many-band generalization as well and it led to the development of the Wannier90 package for computing MLWFs for real materials [24].

B. MLWFs for 1D systems
Let us first consider the 1D case, for completeness and also to demonstrate trouble that arises when extending the treatment to higher dimensions.Here, it turns out that one can always find exponentially localized Wannier functions [25].Furthermore, it turns out that in 1D, the MLWFs defined previously are simply the eigenstates of the projected position operator P xP , where P projects onto the band under consideration.This can be seen quite easily by splitting the variance (2) into a gauge invariant part F I and a gauge dependent part F G as [14] and where the operator Q = 1 − P projects onto the complement of the band under consideration.Thus, if we suppose that the |w R ⟩s are eigenstates of P xP with corresponding eigenvalues λ R , then we find and so we see that F G = 0 for these Wannier functions.Since F G is clearly positive definite and is the only gauge dependent piece, eigenstates of P xP are clearly maximally localized Wannier states.We note that this argument generalizes straightforwardly to the 1D multi-band case.

C. MLWFs for 2D systems
Moving to 2D systems, the situation becomes more complicated.We no longer consider eigenstates of the projected position operators because P xP and P yP do not commute in general and are therefore not always simultaneously diagonalizable.In fact, the commutator of P xP and P yP is directly related to the Chern number [26][27][28].It was shown in [20] that a necessary and sufficient condition for the existence of exponentially localized Wannier functions in 2D is the vanishing of the Chern number of the bands considered.Additionally, the variance becomes infinite if the system has non-zero Chern number [21].As discussed in the introduction, this is due to the presence of irremovable vortices in the Bloch states.Due to these reasons, one may be led to believe that MLWFs cannot be constructed in a sensible way for a 2D Chern insulator [29].

III. THEORY OF OPTIMALLY LOCALIZED WANNIER FUNCTIONS
In this section, we present the main result of this article.We provide a solution to the problem of constructing optimally localized Wannier functions for a 2D Chern insulator in the case of a single isolated band.We emphasize that the method described applies even in the case of a complex band structure as long as the band considered is isolated from the rest of the bands.

A. Berry connection and Bloch state vortices
A Chern insulator is characterized by having a nonzero Chern number, which is defined by [30] where Ω(k) is the Berry curvature which is gauge invariant.The Chern number is a topological invariant and is always an integer.In what follows, we restrict ourselves to the case of C = ±1 to avoid unnecessary complexity.However, it is straightforward to generalize our work to higher Chern numbers and we will remark on this later.One usually defines the Berry curvature via Ω = ∇ k × A but care must be taken when A is not smooth.A standard approach in computing the Chern number is to use different gauges for different patches of the Brillouin zone to ensure A is smooth within each region.Then one simply integrates the curl of the smooth connections over each region to evaluate C. On the other hand, in the standard construction of Wannier functions (1) we require using a single gauge to cover the entire Brillouin zone.If one insists on such a gauge, a straightforward application of Stokes' theorem gives due to the Brillouin zone periodicity of A.
There is no contradiction of the above with the possibility of non-zero Chern numbers if we allow for a vortex in the Bloch states.If there is a vortex at k v , then the Berry connection will also be singular at that point leading to a delta function in its curl.So, we have where δ P (k) = G δ(k+G) is the periodic delta function with summation over reciprocal lattice vectors.Note that the curl of a 2D vector field is shorthand for the scalar quantity ∂A y /∂k x − ∂A x /∂k y .Now, ( 11) is satisfied because of the cancellation between the delta function and the net Berry curvature.(See Appendix E for a related discussion in the context of polarization which will be studied later in this article.)It is crucial to note that for a Chern insulator, having a vortex is inevitable in this construction.Furthermore, the position of the vortex is a gauge dependent property because it can be moved by a singular gauge transformation i.e. a gauge transformation by a phase with a vortex in it.If the gauge transformation has a vortex of winding number −C at the initial position and a vortex of winding number C at some other position, then it moves the position of the vortex.Therefore, the vortex position presents us with an extra degree of freedom to minimize over in addition to the smooth gauge degree of freedom.The latter leads to the Coulomb gauge condition ∇ k • A = 0, which determines the gauge once the position of the vortex is given.Note that for a non-topological system, the Chern number is zero and so there are no vortices to begin with.There, singular gauge transformations should be completely avoided as they would introduce vortices that worsen the decay of the Wannier functions.Thus, one only considers the Coulomb gauge condition in such cases.

B. Electrostatics analogy
We next proceed to develop the correspondence of the minimization problem with two-dimensional electrostatics.We begin with the variance functional for a single band in the momentum space representation, as given by (3).Note that this can be broken down into gauge invariant and gauge dependent parts.In section II.B, we wrote these down in real space form.In momentum space, the gauge dependent part, which is all we need to minimize, takes the form The second term above is gauge invariant under smooth (and periodic) gauge transformations, but we cannot ignore it for the present scenario because it is not invariant under singular gauge transformations.While the minimization of ( 13) over smooth gauge transformations yields the Couloumb gauge condition as usual, this process cannot move the vortex.Therefore, we must solve the Coulomb gauge condition at each vortex position k v and then use the localization functional to identify the optimal vortex position.In other words, we must solve for A satisfying the equations These equations uniquely determine A up to a constant.Furthermore, the integral of A over the Brillouin zone with specified k v is a gauge invariant quantity which fixes this constant.
In order to solve ( 14), we introduce a vector field E such that where c = 1 V BZ BZ dk A and ẑ is a unit vector perpendicular to the 2D plane on which A lies.Then, in this dual picture, E satisfies the electrostatic equations with periodic boundary conditions, where ρ(k; is the corresponding charge distribution.In the above, E is interpreted as the electric field of a charge distribution composed of a periodic array of point charges, one per reciprocal space unit cell, together with a smooth background Ω.Note that the net charge per unit cell is zero, which means that this periodic problem is solvable.We emphasize that the Berry connection A should not be confused with an electromagnetic vector potential in the dual picture.Further note that (15) implies that the average electric field is zero from which it follows that it can be expressed in terms of a periodic potential.Fig. 1 provides an example of a Berry connection for a set of Bloch states with a vortex and the corresponding zero-average electric field.Using (15) to express the localization functional in terms of E, we arrive at the simple relation The electrostatic analogy therefore goes even further.In the dual picture, the localization functional becomes the electrostatic energy per unit volume of an infinite periodic charge-neutral system.Minimizing over the singular gauge transformations corresponds to finding the positions of the point charges which minimize this electrostatic energy.

C. Analytical solution
We now proceed to find analytical expressions for quantities entering the localization functional.We split the charge distribution into two parts as ρ(k; corresponds to a periodic lattice of point charges in a uniform neutralizing background given by ρ 0 = C/V BZ and corresponds to a smooth periodic charge distribution in the same uniform neutralizing background.We will solve ( 16) for these two charge distributions separately and thus obtain a solution to the entire problem by linearity.For the charge distribution ρ v (k; k v ), one typically uses the Ewald summation technique to find the electric field [31][32][33].In the present problem, we are in 2D so we can instead leverage results from complex analysis.For each wave vector k, we introduce a corresponding complex representation given by z = k x + ik y .Furthermore, we represent the electric field E v (k; k v ) for this problem in complex form as E v (z; z v ) = E x − iE y , where z v is the complex representation of the vortex position k v .Then, following ideas from work by Tkachenko on superfluid vortex lattices [34], we find that the solution to the electric field is given by where α is a constant chosen to make E v periodic (see Appendix C for more details) and ζ(z) is the Weierstrass zeta function given by Here, G = G x + iG y is the complex representation of reciprocal lattice vector G.
For this electric field, the corresponding electrostatic potential is given by We note that this electrostatic potential is both real and periodic (see Appendix C).The solution is much simpler for the smooth charge distribution ρ s (k).We can write the Fourier series representation where the Ω R s are the Fourier components of the Berry curvature.Note that there is no R = 0 term in the above equation due to charge neutrality of ρ s (r).By taking a Fourier transform of Poisson's equation, we therefore obtain the solutions FIG.
1.An example of a Berry connection for a set of Bloch states with a vortex (left) and the corresponding zero-average electric field (right).The colors represent the magnitudes of the vectors.Note that there is just one vortex and this is located at (π, π) (the Berry connection magnitude vanishes at other points where its direction has a defect). and for the electric field and electrostatic potential respectively.
The solution to the full problem ( 16) is obtained by summing together the two solutions E v and E s we found above.Note that the total potential ϕ = ϕ v + ϕ s is periodic in reciprocal space so the electric field obtained indeed has zero-average.

D. Minimizing F
We now proceed to describe the minimization of the localization functional.Using the solution found in Section III.C, we can solve the problem of minimizing the localization functional F .In the electrostatics picture, we have where the constant term does not depend on the vortex position.The non-constant term follows from integrating by parts the term involving the interaction between the point particle and the smooth charge distribution.Note that this immediately gives us a sensible way to regularize the blowing up of the variance mentioned previously.In the electrostatics picture, this infinity simply corresponds to the self energy of the point charges, which is contained in the term involving E 2 v .This does not depend on where the charges are placed and so it can be neglected in the usual fashion.In the variance picture, this term corresponds to a logarithmic divergence which is independent of the vortex position.Therefore, we can regularize the localization functional by introducing a long distance cutoff and thus effectively minimize the finite short-range contribution to the localization functional over the vortex position.This will be discussed further in the next section.Moreover, (27) tells us that to find optimal Wannier functions for C = +1(−1) in our setting, the vortex needs to be placed at the maximum (minimum) of the electrostatic potential ϕ s for the smooth charge distribution in a uniform neutralizing background.
As stated in the beginning of this section, the generalization to higher Chern numbers is relatively straightforward.For a Chern number |C|= 1, there is only one vortex per unit cell.The vortex-vortex interaction contribution to the electrostatic energy will not depend on the vortex position because changing the vortex position corresponds to uniformly translating the vortex lattice.Thus we have lumped the vortex-vortex interaction into the constant term in (27).This, however, does not remain true for |C|≥ 2 where there will be two or more vortices in the Brillouin zone.Here, one must explicitly consider the vortex-vortex interaction terms involving (22) in the minimization which makes the problem richer.For this case, there will generally be two competing terms in (27): one describing vortex-vortex interaction and the other describing the vortices interacting with the smooth potential.

IV. NUMERICAL RESULTS
In this section, we present numerical results to test the developments of the previous section.We use a minimal Chern insulator model that is related to the Qi-Wu-Zhang model [35] via a spin rotation.The lattice Hamiltonian considered is This is a Hamiltonian on a lattice with sublattices A and B where the creation operators a † i,j and b † i,j create a particle at the A and B sites of the (i, j) unit cell respectively.We set the lattice constants to one so that the distance between the A and B sites in a given unit cell is 1/2.The Hamiltonian involves nearest-neighbor hoppings in the x and y directions and also cross hoppings between A and B sites of different unit cells, as shown in Fig. 2. One of the hoppings is imaginary in order to break time-reversal symmetry (similar to the Haldane model's complex next-nearest-neighbor hoppings).
Since this Hamiltonian is periodic in real space, one can Fourier transform to obtain the following Bloch Hamil-tonian where σ x,y,z are the Pauli matrices.Note that we have chosen to Fourier transform the Hamiltonian in such a way that it is not periodic in reciprocal space.This has been done to incorporate the different positions of the A and B sublattice sites within the unit cell into the real space position operator.This choice of position operator is much more intuitive for our discussion since the Wannier representation is an inherently real space quantity.The Hamiltonian (29) has lowest-band Chern number C = 1 for 0 < u < 2, C = −1 for −2 < u < 0 and C = 0 for |u| > 2. For the numerical work in this section, we set u = 1 so that the Chern number is always +1.Throughout we will consider Wannier functions of the lower band.
In order to compute optimal Wannier functions for this model, we begin by discretizing the Brillouin zone into N equally spaced points in each direction.Moreover, we need to make an initial choice of gauge for the Bloch states.Let us begin by fixing the first entry of the eigenvector of H(k) to always be real.It can then be seen that this fixes the vortex in the position (k x , k y ) = (π, 3π/2).
We now require a general singular gauge transformation capable of shifting this vortex to an arbitrary location in the Brillouin zone.Such a gauge transformation is not unique but we can fix it by also demanding that it leaves the value of ∇ k • A unchanged.This is a useful choice because it allows us to carry out minimizations over smooth and singular gauge transformations in any order.It follows directly from the solution to the electrostatic problem in the previous section (see also Appendix D) that the singular gauge transformation shifting a vortex from position a to b (where both are written in complex representation), is given by e iθs(z;a,b) , where Here, z = k x + ik y is once again the complex representation of the reciprocal space vector k and α is the same constant as in (20).This gauge transformation satisfies ∇ 2 θ s = 0 by construction (thus preserving ∇ k • A) and e iθs(z;a,b) is Brillouin zone periodic.Fig. 3 (a) shows a vector plot of the singular gauge transformation for N = 29 and a = π + 3πi/2 and b = π + πi/2.The arrows represent the direction of the phase e iθs on the unit circle.We observe that this gauge transformation creates an anti-vortex at a and a vortex at b, thereby completely shifting the vortex from a to b without leaving any trace behind.Next, we compute the optimal Wannier functions for our lattice model using the analytical result from the previous section.Namely, noting that we have C = +1, we first move the vortex to the maximum of the smooth potential ϕ s (k) and then minimize the variance functional with respect to smooth gauge transformations.The latter is done via a steepest descent algorithm.In our single band case, this amounts to solving the differential equation where U (k, τ ) = e iθ(k,τ ) is the smooth unitary gauge transformation acting on the Bloch states and τ is a suitably scaled relaxation time variable.Note that we always time-evolve U in such a way that unitarity is preserved at each step.It is crucial to note that the steepest descent process cannot alter the vortex position by definition.Fig. 3 (c) shows the density of the resulting optimal Wannier function in real space.In contrast, Fig. 3 (d) shows the Wannier function for a gauge satisfying ∇ • A = 0 but with the vortex at (π/4, π/4), which is sub-optimal.It is clear that this Wannier function overlaps non-trivially with a neighboring unit cell and so the importance of choosing the right vortex position is evident.
We now move on to consider the variance of the optimally localized Wannier functions we have found.For an unbiased check of the analytical arguments, we will endeavor to compute the variance working entirely within real space.This involves finding a real-space method of regularizing the divergence of the variance.As we have seen, in k-space the divergence is regularized simply by dropping the vortex self-interaction term.In real space, we expect that the variance should be regularized by long-distance modifications.The variance in real space can be written in the form where w(r) is the Wannier function and V (r) = r 2 /2.
(The idea of expressing localization functionals in this form is not new -two such examples are the Edmiston-Ruedenberg criterion [3] and the von Niessen criterion [36].)The regularization procedure then is to alter the "potential" V (r) at large distances to achieve a convergent result.We work with the potential given by where R is a cutoff parameter.We observe that Ṽ (r) ∼ r 2 /2 for r ≪ R while Ṽ (r) ∼ Rr for r ≫ R. Thus, this choice of potential has the desired properties, given the 1/r 2 decay of the Wannier functions.For our calculations, we require a value of R that is large enough to capture the finite part of the variance.This can be done by choosing R such that R 2 is large compared to the ratio of the band curvature and the band gap, at k = 0. Thus, we now have a sensible choice of cut-off for the localization functional which can be applied to arbitrary system sizes (which are larger than the cut-off).We will refer to this version of the localization functional as the modified variance.
With this, we can verify that the optimal Wannier functions are indeed found by placing the vortex in the maximum of the potential and that the dependence of the localization functional on the vortex position is of the form given by (27).We carry out the steepest descent procedure while placing the vortex in different positions across the Brillouin zone.For each vortex position, we evaluate the modified variance of the resulting Wannier function using a cutoff radius of R = 30.The result is shown in Fig. 3(f).We compare this with Fig. 3(e), which is a plot of the negative of the modified smooth potential φs (see (27)), which is defined to be We see that these two plots differ only by an additive constant as predicted by ( 27), thus validating our analytical result.

V. APPLICATION: POLARIZATION
Having developed and tested optimally localized Wannier functions for Chern insulators, we now move on to describe an application.In particular, we will show how they can be used to compute the polarization of a Chern insulator in an unambiguous and physically relevant way.
For topologically trivial systems, Wannier functions provide an elegant means of computing the electric polarization, following the modern theory of polarization [9,37].The fundamental idea is that an absolute polarization cannot be defined uniquely for an extended system, but the change in polarization of a system subject to some deformation is well-defined and experimentally measurable.It can also be shown that this change in polarization is given by the change in the position of the Wannier center in real space, providing a very intuitive picture [14].
Polarization for a Chern insulator, on the other hand, has been a topic under some amount of debate, with recent works studying this extensively (see, for example, [38][39][40]).Indeed, the polarization is conventionally determined by integrating the Berry connection over the Brillouin zone.However, as we have discussed, this quantity is not invariant with respect to singular gauge transformations.Therefore, a consistent mechanism of fixing the singular gauge is crucial.We will show that the construction we have for optimally localized Wannier functions provides a very natural way of fixing the gauge and thereby obtain an expression for the polarization in terms of the Wannier center shift plus a correction.Our results are fully consistent with those of [41] where a k-space expression for the polarization of a Chern insulator is developed.
Polarization is encapsulated by the fundamental relation [9] ∆P n = 1 0 dλ J n (λ), (35) which describes the change in polarization as the integral of the current during the deformation process.Here, ∆P n is the change in polarization of the n th band and λ parametrizes the adiabatic deformation of the Hamiltonian (λ = 0 is usually the centrosymmetric case).Furthermore, (36) is the cell-averaged adiabatic current where −e is the charge quantum, and we note that the Bloch states are now functions of λ.This adiabatic current is a measurable and therefore a gauge invariant quantity, thus providing the correct starting point for considering polarization in a Chern insulator.
We now recall the Berry connection A nk = i⟨u nk |∇ k u nk ⟩ and introduce A nλ = i⟨u nk |∂ λ u nk ⟩.In terms of these, we may write the adiabatic current as The important term here is the final one.Mixed partial derivatives commute when acting on a smooth function but for a Chern insulator, they would not commute at vortices and so we cannot ignore this term.Therefore, for a band with C = ±1, we find that (see Appendix E) where ∆r n is the change in the Wannier center, ∆k vn is the change in the vortex position of the Bloch states due to the adiabatic deformation, and V c is the area of the unit cell.This relation ( 38) is another main result of the article.For conventional (non-topological) systems, the correct expression for the polarization is given by the first term of (38) only.However, this term is not invariant with respect to singular gauge transformations and the additional term in (38) serves to cancel off this gauge dependence, making the full expression of the polarization gauge invariant as it must be.We now wish to use the developed formalism to compute the polarization of our model system.In order to induce a change in polarization, we add a staggered onsite potential which breaks centrosymmetry and so we consider the modified Hamiltonian where δ is the onsite potential.We set u = 0.5 so that the hopping amplitudes alternate in strength in the xdirection.Starting at the centrosymmetric case δ = 0, we compute the change in polarization as the staggering is adiabatically increased up to δ = 0.75.We note that we cannot cross δ = √ 3/2 ≈ 0.866 since the gap would close, breaking adiabaticity.Throughout we fix the average particle number to be one particle per every two lattice sites.By considering the symmetry of Fig. 2, we see that turning on δ would only induce a polarization in the x-direction.We use optimally localized Wannier functions in each case to compute the polarization via (38), where we use an 80 × 80 grid of points in reciprocal space.The results, which can be computed fairly efficiently, are shown by the blue curve in Fig. 4.
To check the validity of the expression for the polarization, we now turn to numerically computing the dipole moment of our model for large finite system sizes, mostly following the method described in [41].We take the system to have open boundary conditions in the x-direction (which is the direction of the induced polarization) and periodic boundary conditions in the y-direction.For a system with M unit cells along the x-direction, the dipole moment per unit length is given by where V c is the area (volume in higher dimensions) of the real space unit cell, n is the band index and k is the wave number.Note that k is now a scalar because we only have periodicity in the y-direction, and the number of bands involved scales with M .There is then a question of how to fill the edge states when computing the dipole moment as δ is varied.Here we appeal to a sensible physical scenario.We start with the zero-temperature ground state when δ = 0. We take the rate of variation to be slow with respect to the bulk energy gap, but fast with respect to edge-to-edge tunneling which exponentially decays as a function of system length.This means that, during the course of the evolution, the population of edge states may differ from those of instantaneous ground states.Specifically, the final state will acquire a different chemical potential on each side through the quasi-adiabatic evolution.This is how occupied states in (40) are determined.
The change in dipole moment with respect to δ = 0, calculated as prescribed above, is shown by red dots in Fig. 4. We compute the dipole moment for two widths M = 20 and M = 30, and two k-space discretizations with 500 and 1000 k-points.We then apply Richardson extrapolation in both M and number of k-points to eliminate the leading order errors in both quantities.This gives us an excellent approximation for the dipole moment in the thermodynamic limit M → ∞, so that it may be compared with the polarization.Indeed, we observe very good agreement as seen in the plot, where the largest relative error is only 0.17% and occurs for δ = 0.75, which is in the vicinity of where the gap closes.

VI. SUMMARY AND OUTLOOK
In this article, we have proposed a means of computing optimally localized Wannier functions for a single band in Chern insulators.We showed that the problem of minimizing the variance functional for such systems has a dual interpretation involving minimizing the electrostatic energy of a periodic array of point charges in a smooth neutralising background.In this dual picture, the electrostatic energy naturally involved an infinite self-energy term.We observed that if this term were neglected in the usual fashion, then we obtained a finite expression which depended non-trivially on the positions of the vortices.Minimizing this finite part corresponds to minimizing the short range spread of the Wannier function around its center so that the overlap with neighboring atomic sites is minimized.This makes our optimal Wannier functions useful in practical calculations.We used Chern numbers of ±1, where there is always only one vortex, for clarity and showed that the vortex position minimizing the electrostatic energy was an extremum of the electrostatic potential corresponding to the (neutralized) smooth background.We then remarked that everything works similarly for higher Chern numbers in the case of a single band -the electrostatic energy simply involves extra terms due to pairwise vortex interactions (since there can be more than one vortex) and it must be minimized over all over vortex arrangements.
As an immediate application of our work, we presented a Wannier function approach to computing polarization in a Chern insulator.In doing so, we showed that the fact that the Wannier center is no longer gauge invariant in the usual sense is not problematic given that we start with a physical gauge invariant quantity (namely the adiabatic current).Moreover, we verified our result by comparing it with the dipole moment per unit length which is also a physical quantity.
Our work leaves many avenues open for future research.Foremost among them is the generalization to multi-band settings.There, one must minimize the localization functional over all gauge transformations within that subspace of bands.This includes singular non-Abelian gauge transformations which dictate the distribution of the vortices across the multiple bands.Furthermore, multi-band problems often involve band crossings which give rise to further discontinuities.In the nontopological case, such discontinuities can be removed using smooth non-abelian gauge transformations [14].We speculate that it may be possible to do the same for band crossings when the total Chern number is non-zero.If this is possible, then we can pick a gauge in which one of the new bands (no longer eigenstates of the Hamiltonian) has all the vortices while the rest of the new bands are smooth.In that case, all but one of the Wannier functions would be exponentially localized and so it may be possible that such a scenario leads to the optimal choice.An additional generalization is the case of a smaller subspace of states disentangled from a larger set of entangled bands (for example by restricting to states within a specific energy window).This is discussed for conventional systems in [14].
Another issue to consider is that of Wannier functions displaying spatial symmetries of the underlying Hamiltonian.This is a topic that has been studied in depth in many different settings since the early days [42][43][44][45].In particular, there have been examples such as the s-like Wannier function of Cu [46], where the optimal Wannier function does not exhibit the right symmetries of the system.As a result, methods have been introduced to obtain maximally localized symmetric Wannier functions for non-topological systems by extending the methods of Marzari and Vanderbilt [47].It remains to be studied whether the optimally localized Wannier functions we propose preserve key symmetries or whether our methods can similarly be extended to obtain symmetric Wannier functions by trading off some amount of localization.Indeed, to study this fully, one must first extend our work to multiple bands.Moreover, the fact the Wannier center in our case can be shifted arbitrarily via singular gauge transformations may be crucial to this discussion.
Other potential applications of our work remain to be investigated.For example, our optimal Wannier functions may be used to construct minimal tight-binding models for Chern insulating bands in a similar fashion to techniques applied in topologically trivial settings.Finally, how the Wannier construction described in this article applies to so-called fragile topological phases [13] warrants further consideration.
where we used the fact that ⟨u nk |∇ k u lk ⟩ = −⟨∇ k u nk |u lk ⟩.
Summing together (A4) and (A5) and noting that l |u lk ⟩ ⟨u lk | = 1, we obtain from which it follows that Appendix B: Summary of special functions used In this appendix to be self-contained, we provide a short summary of the Weierstrass functions, which are special functions central to our solution of the electrostatics problem.Tailoring to the problem of interest, we cast everything in reciprocal space.The classic text Ref. [48] is recommended for additional details.
Elliptic functions are doubly periodic functions in the complex plane and their only singularities are a finite number of poles per unit cell.In order to connect this directly to our problem, let us cast the 2D reciprocal space in complex form, as discussed in the main text.That is, we define z = k x + ik y to be the complex number representing the reciprocal wave vector k = (k x , k y ).Then, we can identify Brillouin zone periodicity with double periodicity in the complex representation.
With this setup, the Weierstrass ℘-function in complex reciprocal space is defined as where G = G x + iG y is the complex representation of reciprocal lattice vector G.It can be seen that this function is even and it has poles at reciprocal lattice vectors.Moreover, let the complex primitive reciprocal lattice vectors be denoted 2ω 1 and 2ω 2 .Then, it can be shown that ℘(z) is doubly periodic with half-periods ω 1 and ω 2 .(Here, we choose ω 1 and ω 2 such that ω 2 /ω 1 has positive imaginary part.)In other words, we have Therefore, ℘(z) is an elliptic function and it has the required Brillouin zone periodicity.This function, however, is not immediately helpful in solving the electrostatics problem from the main text (16) as we require simple, not double poles.Instead, we use the Weierstrass zeta function which is defined by together with the condition By term-by-term integration, it can be shown that the ζ-function has the series representation which exhibits the desired simple poles.
In contrast to the ℘-function, the ζ-function is not elliptic because it is only quasi-periodic in the Brillouin zone.To see this, one can first integrate (B3) to find where the last term is a constant of integration.Next, note that ζ(z) is an odd function as can be observed from its series expansion.Substituting z = −ω i into B6 now allows us to determine the integration constants η i = ζ(ω i ) which are generally non-zero.We now state an important relation between the ω i 's and the η i 's, namely This can be derived by taking a contour integral of ζ(z) around the boundary of a unit cell.Using the residue theorem together with the fact the integrals along opposite edges are related via periodicity, we can obtain (B7).
(This result relies on our assumption that ω 2 /ω 1 has positive imaginary part since this determines the handedness of the integration contour.) The final special function we require, for the purpose of constructing the electromagnetic potential, is the Weierstrass sigma function, which is best defined via its relationship to the ζ-function as In this appendix, we provide more details for the solution of the electrostatic equations for a lattice of point charges in a uniform neutralizing background.We recall that the electrostatic equations to be solved are given by where k v is the position of the vortex in the Brillouin zone.
The electric field for a single point charge in an infinite 2D plane is given according to Coulomb's law as In the complex representation, this electric field is given by E(z) = E x − iE y = 1/(z − z v ) where z = k x + ik y and z v is the point charge position as usual.Moving to the periodic case, we can try to extend the point charge solution in the most obvious way as where G represents reciprocal lattice vectors.This sum, however, does not converge which is problematic.
The solution is to use the Weierstrass ζ-function introduced in the previous section.Informally, we can understand the ζ-function as a regularization of (C2) -that is, it introduces additional terms which make the new sum convergent outside of the poles.Despite fixing the divergence, and having poles in the correct locations, the ζ-function is not doubly periodic.This can be dealt with by adding terms linear in k x and k y .In particular, a direct calculation using relations from the previous Appendix shows that the expression has full Brillouin-zone periodicity.A similar expression was obtained in early work on vortex lattices in superfluid helium [34].The quantity α depends solely on the lattice geometry and for a square lattice it can be found that α = 0.This can be observed by noting that ω 2 = iω 1 and using (B5) to show that η We now show that the electric field found indeed satisfies the required equations (C1).Using the identity d/dz * ((z − z v ) −1 ) = πδ(k − k v ), and noting the ζ function has simple poles, one finds Noting that E is real, the equations (C1) are satisfied.Finally, as stated in the main text, the electrostatic potential for this problem is given by This is also periodic as claimed, which can be seen by using the quasi-periodic conditions (B10) for the σfunction.

Appendix D: Derivation of Singular Gauge Transformation
In this appendix, we derive Equation ( 30) from the main text, which is the formula for the singular gauge transformation which shifts the vortex position from one point in the Brillouin zone to another while preserving the Coulomb gauge condition.
We begin with the fact that under a gauge transformation of the form |u nk ⟩ → e iθ(k) |u nk ⟩, the Berry connection transforms as A → A − ∇ k θ.If we move the vortex from position a to b, then the Berry connection initially and finally must satisfy respectively.If we demand the additional condition that this gauge transformation keeps ∇ k • A unchanged (so that we can shift the vortex after minimizing over continuous gauge transformations in one position without having to carry out steepest descent again), then it follows from the above that v = ∇ k θ must satisfy where v = v x − iv y and S is a different constant.
Our aim is to obtain an expression for θ from (D4).To achieve this, we note that (D2) are fluid equations for a periodic incompressible potential flow around a vortex and an anti-vortex.The complex velocity is a meromorphic function of z.So using methods familiar from fluid dynamics in 2D, we proceed by introducing a complex potential f = θ + iψ where θ is the velocity potential, corresponding precisely to our gauge function, and ψ is the stream function.This complex potential is defined to satisfy Note that we have thrown out the constant of integration in (D6) because this would simply correspond to a global phase, which we can ignore.Now, our singular gauge θ(z; a, b) is given by the real part of (D6) and we require that e iθ is periodic in the Brillouin zone.We now use this condition to obtain the constant S. It follows from (B10) that

Appendix E: Polarization
In this appendix, we provide a detailed derivation of equation (38) for the change in polarization in a Chern insulator with C = ±1.We recall once again that the change in polarization is defined by [9] ∆P n = This can be written in terms of the Berry connection A nk = i⟨u nk |∇ k u nk ⟩ and the new quantity A nλ = i⟨u nk |∂ λ u nk ⟩ as (E3) It is easy to see that the first term in the current gives the change in Wannier center when integrated over λ, while integrating the second term gives zero due to periodicity of the Berry connection.Let us now focus on the last term in the current.This clearly vanishes except at the vortex itself.Therefore, let us consider a small disk D ϵ of radius ϵ around the vortex.In such a region, we may write where |ũ nk ⟩ is smooth and k vn (λ) is the vortex position.Thus, in our small region, we have extracted the vortex out of the Bloch state and into a phase factor.(It is not possible to factor out a phase from |ũ nk ⟩ over the whole Brillouin zone without creating a new vortex somewhere else, as vortices cannot be eliminated from a band with a Chern number, but we will use this expression only within D ϵ .) In D ϵ , it is true that and therefore, where ∂ x , ∂ y are taken with respect to k x and k y respectively, and y vn is the y-component of k vn .In the first step, we have used the fact that the derivatives commute except when applied to the phase.Note that we used (E8) for the third equality.A similar expression can be found for the other component.Then substituting the expression found for the adiabatic current density back into (E1) we obtain the desired result for the polarization (38).

FIG. 2 .
FIG.2.Diagram of the real space lattice connected to the Hamiltonian given by(29).Sites on sublattices A and B are colored red and blue respectively.The arrows indicate the hopping terms.

FIG. 3 .
FIG. 3. Numerical Results for the Lattice Model.N = 87 is used everywhere except for the vector plot.(a) Vector plot of the singular gauge transformation for N = 29 which moves the vortex from (π, 3π/2) to (π, π/2).The arrows are the vectors formed by the real and imaginary parts of this gauge transformation.(b) Magnitude of the Berry connection after applying the singular gauge transformation from (a).(c) Zoomed-in plot of the optimal Wannier function in the home unit cell.This has a modified variance (see text) of 0.8461.(d) Zoomed-in plot of the Wannier function obtained by optimising over smooth gauges but placing the vortex at the sub-optimal position (π/4, π/4).This has a modified variance (see text) of 0.9758.(e) Plot of minus the modified smooth potential φs(k) = (4π/VBZ)ϕs(k).(f) Plot of the modified variance (see text) of the optimal Wannier function obtained by placing the vortex in different points across the Brillouin zone.The number of different vortex positions used is 441.

FIG. 4 .
FIG.4.Plots of change in polarization (blue) for the lattice model with respect to the centrosymmetric case and change in dipole moment (red dots) across an infinite strip of the lattice model.Both quantities are given in units of −e/a where a is the lattice constant.A reciprocal space grid with N = 80 was used to compute the polarization.To compute the dipole moment, Richardson extrapolation was applied with two strip widths M = 20 and M = 30 and two reciprocal space discretisations with 500 and 1000 k-points respectively.