{'fraction_non_alphanumeric': 0.039455563716489166, 'fraction_numerical': 0.012578038927653324, 'mean_word_length': 4.05675487465181, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wannier representation of Z_2 topological insulators

We consider the problem of constructing Wannier functions for Z_2 topological insulators in two dimensions. It is well known that there is a topological obstruction to the construction of Wannier functions for Chern insulators, but it has been unclear whether this is also true for the Z_2 case. We consider the Kane-Mele tight-binding model, which exhibits both normal (Z_2-even) and topological (Z_2-odd) phases as a function of the model parameters. In the Z_2-even phase, the usual projection-based scheme can be used to build the Wannier representation. In the Z_2-odd phase, we do find a topological obstruction, but only if one insists on choosing a gauge that respects the time-reversal symmetry, corresponding to Wannier functions that come in time-reversal pairs. If instead we are willing to violate this gauge condition, a Wannier representation becomes possible. We present an explicit construction of Wannier functions for the Z_2-odd phase of the Kane-Mele model via a modified projection scheme followed by maximal localization, and confirm that these Wannier functions correctly represent the electric polarization and other electronic properties of the insulator.


I. INTRODUCTION
In the past several years there has been a surge of interest in topological insulators. These are materials that are gapped in the bulk, just like ordinary insulators, but that cannot be adiabatically connected to ordinary insulators without closing the gap or breaking some specified symmetries. They also exhibit chiral metallic edge states that are topologically protected from disorder. [1][2][3] Topological insulators can be distinguished from normal ones based on the manner in which the Bloch eigenfunctions are topologically twisted in k-space.
Two types of topological insulators have received the most attention. First, Thouless et al. 4 pointed out long ago that a two-dimensional (2D) insulator is characterized in general by a topological integer known as the "Chern number" or "TKNN index." A prospective insulator having a non-zero value of this integer would be known as a "Chern" or "quantum anomalous Hall" insulator. The latter name arises because such a crystal would exhibit a quantum Hall effect (QHE) even in the absence of a macroscopic magnetic field, and would have chiral edge states just like the ordinary field-induced QHE. Haldane devised an explicit tight-binding model realizing such a case. 5 Since the Hall conductance is odd under the time-reversal (T ) operator, Chern insulators can only be realized in systems with broken T symmetry, e.g., insulating ferromagnets. Despite the fact that these possibilities have been appreciated now for almost three decades, no known experimental realizations of a Chern insulator are yet known.
Second, a great deal of interest has surrounded the recent discovery of a different class of topological insulators known as Z 2 insulators that realize the quantum spin Hall effect (QSH). 6 Subsequent theoretical 7-9 and experimental [10][11][12][13][14] work has succeeded in identifying several materials systems that realize the case of a Z 2 topological insulator. Unlike the Chern index, which van-ishes unless T is broken, the Z 2 index (which takes values of 0 and 1, or equivalently, "even" and "odd") is only well defined when T is conserved. Z 2 insulators are thus non-magnetic, although a spin-orbit or similar interaction is needed to mix the spins in a non-trivial way. Because T is preserved, the occupied states at k and −k form Kramers pairs, and one can associate a Z 2 invariant with the way in which these Kramers pairs are connected across the Brillouin zone. 15 Since the Z 2 index cannot change along an adiabatic path that is everywhere gapped and T -symmetric, a Z 2 -even (normal) insulator cannot be connected to a Z 2 -odd (topological) one by such a path. In 2D there is a single Z 2 invariant, and T -invariant insulators are classified as "even" or "odd," while in 3D there are four Z 2 invariants and the classification is more complicated. 16 Wannier functions (WFs) have proven to be a valuable tool when working with semiconductors and insulators, providing a real-space description that can be used to understand bonding, construct model Hamiltonians, and directly compute certain physical properties such as the electric polarization. 17,18 Thus, it is desirable to understand the construction of the Wannier representation for topological insulators so that this useful set of techniques can be applied to these novel materials.
For Chern insulators it has been shown that a nonzero Chern number presents a topological obstruction that prevents the construction of exponentially localized WFs. 19,20 Conversely, a general proof has been given that exponentially localized WFs should exist in any 2D or 3D insulator having a vanishing Chern index. 21 In principle this applies to Z 2 -odd as well as Z 2 -even T -invariant insulators, suggesting that a Wannier representation should be possible in both cases. However, it is unclear whether the nontrivial topology of the Z 2 -odd case has any effect on the Wannier representation. In particular, one may wonder whether the procedure for obtaining WFs would be the same as for ordinary insulators, and if not, how it should be modified in order to get well localized WFs in the Z 2 -odd regime.
In this paper we address this question using the model of Kane and Mele 6 as a paradigmatic system that exhibits both Z 2 -odd and Z 2 -even phases. We demonstrate that the usual projection scheme used for constructing the Wannier representation is still applicable to the Z 2odd insulators, but only for gauge choices that do not allow WFs to come in time-reversal pairs. We present an explicit projection procedure for constructing welllocalized WFs in the topologically non-trivial phase, and show that the WFs can be made even more localized using the standard maximal-localization procedure. 17 We also discuss the electric polarization from both Berryphase and Wannier points of view, showing the relations between the viewpoints and confirming that both give identical results.
The paper is organized as follows. In Sec. II we define the Z 2 topological invariant in 2D and briefly discuss methods for determining it numerically. We review the model of Kane and Mele in Sec. III, and describe its spectrum and phase diagram. In Sec. IV we present the projection scheme used to construct WFs and explain how the application of this scheme to Z 2 -odd insulators is different than for ordinary insulators. The localization properties of the constructed WFs are described in Sec. V. The electric polarization properties and locations of the Wannier charge centers are considered in Sec. VI. Finally, we make concluding remarks in Sec. VII.

II.
Z2 INVARIANT Here we briefly review some of the equivalent ways of determining the Z 2 invariant in 2D insulators.
In the work of Ref. 22 the definition of the Z 2 invariant was given in terms of a function P (k) defined as i.e., the Pfaffian of a certain k-dependent antisymmetric N ×N matrix, where N is the number of occupied bands.
Here |u j (k) = e −ik·r |ψ j (k) is the periodic part of the Bloch function of the j'th occupied band andθ = is yĈ is the time-reversal operator (Ĉ is complex conjugation and s y is the second Pauli matrix). If the zeros of P (k) are discrete, then the Z 2 invariant is odd if the number of zeros of the Pfaffian within one half of the Brillouin zone (BZ) (see Fig. 1) is odd, and even otherwise. If the zeros of the Pfaffian occur along lines in the BZ, then the Z 2 invariant depends similarly on whether half the number of sign changes of P (k) along the boundary of the half BZ is odd or even. Using ∆ = 0 and 1 to represent evenness and oddness respectively, the Z 2 invariant can equivalently be determined as 6 where the loop integral runs along the boundary ∂τ of the half BZ, and the δ term is included for convergence. Another approach to the problem of defining ∆ results from considerations of "time-reversal polarization." 23 Here a spin-pumping cycle is considered and it is shown that the Z 2 index is given by the difference between the time-reversal polarizations at the beginning and the midpoint of the cycle. This approach leads to the formula where w mn (k) = u m (−k)|θ|u n (k) and Γ i are the four time-reversal invariant points of the BZ (i.e., those for which −Γ i = Γ i + G with G a reciprocal vector). Note that the matrix w mn is not the same as that in Eq. (1). The definition in Eq.
(3) appears to require a knowledge of the occupied wavefunctions at only four points in the BZ, unlike Eq. (2), for which the wavefunctions must be known at all points along the boundary of the half BZ. However, Eq. (3) is usually not suitable for numerical implementation in practice, since the sign of the Pfaffian at any one of the four points can be flipped by a relabeling of the Kramers-degenerate states at that point. To be more explicit, there is a "gauge freedom" in the choice of states |u m (k) , corresponding to a k-dependent N × N unitary rotation among the occupied states. Eq. (3) is only meaningful when a globally smooth gauge choice enforces a relation between the labels at the four special k-points. 23 This problem may be avoided in the presence of some additional symmetry that can be used to establish the labels of the bands at these points. For example, in Ref. 9 it is shown how the presence of inversion symmetry allows for a simplified calculation of ∆ from Eq. (3).
In the absence of inversion symmetry, one can use yet another definition of the Z 2 index taking the form 23 where A = i N n=1 u n |∇ k |u n is the Berry connection of N occupied states and F = ∇ k × A is the corresponding Berry curvature. 24 Of course, if A and F are both constructed from a common gauge that is smooth over τ , the result would vanish by Stokes' theorem. Thus, Eq. (4) is only made meaningful by the additional specification 23 that the boundary integral of A must be calculated using a gauge that respects time-reversal symmetry, i.e., For the case of the nontrivial Z 2 state, it turns out to be impossible to choose a gauge that satisfies both smoothness over τ and the constraint (5) over ∂τ . In other words, ∆=1 signals the existence of the topological obstruction.
To see how this works more explicitly, the contributions to the integral of A over ∂τ are illustrated in Fig. 1. We choose a gauge that is periodic, |u j (k) = |u j (k+G) , in addition to satisfying Eq. (5). The contributions of the top and bottom segments (solid blue arrows in Fig. 1) then cancel because they are connected by a reciprocal lattice vector G. Thus, the gauge needs to be fixed only along the left and right boundaries (composed of red dashed and gray dotted arrows in Fig. 1), which are separated by a half reciprocal lattice vector. At each of the special points Γ i , one state from each Kramersdegenerate pair is arbitrarily identified as |u 2n−1 (Γ i ) , and the other is constructed via Then we can make an arbitrary gauge choice along the remaining portions of the gray dotted arrows in Fig. 1 e.g., accepting the output of some numerical diagonalization procedure. Finally, the gauge should be transferred to the dashed-arrow segments using Eq. (5), where k and −k belong to the dotted and dashed segments respectively. Eq. (4) can now be evaluated using a uniform discretized mesh K covering the region τ , with the timereversal constraint applied to the boundary ∂τ as described above. To do so, define the link matrices M µ,nm (k) = u n (k)|u m (k + s µ ) and the unimodular link variables L µ (k) = det M µ /| det M µ |, where k ∈ K and s 1 (s 2 ) is the step of the mesh in the direction of the reciprocal lattice vector G 1 (G 2 ). By defining A 1 (k) = log L 1 (k) and one can write the lattice definition of the Z 2 invariant as For a sufficiently fine mesh there will be no ambiguity in the branch choice for the complex log in Eq. (7), since the argument of the log must approach unity as the mesh becomes dense. Moreover, a change in the branch choice determining one of the boundary links A s (k) has no effect (mod 2) on Eq. (7), since each A s (k) appears twice as a result of the gauge-fixing on the boundary. Thus, once the mesh is fine enough so that the branch choices in Eq. (7) are all unambiguous, Eq. (8) gives ∆ exactly. 25

III. THE KANE-MELE MODEL
In their remarkable paper introducing a Z 2 topological classification to distinguish a QSH (Z 2 -odd) insulator from an ordinary (Z 2 -even) insulator, Kane and Mele (KM) 6 also introduced a model tight-binding Hamiltonian that describes a 2D Z 2 -odd insulator in some of its parameter space. In this section we will describe some of the properties of the model suggested therein.
The KM model is a tight-binding model on a honeycomb lattice with one spinor orbital per site. The primitive hexagonal lattice vectors are a 1,2 = a/2( √ 3ŷ ±x) and sites A and B are located at t A = aŷ/ √ 3 and t B = 2aŷ/ √ 3 respectively. The KM Hamiltonian is where the spin indices have been suppressed on the raising and lowering operators, and t is the nearest-neighbor hopping amplitude. In the second term, λ SO is the strength of the spin-orbit interaction acting between second neighbors, with ν ij = (2/ √ 3)[d 1 ×d 2 ] = ±1 depending on the relative orientation of the first-neighbor bond vectorsd 1 andd 2 encountered by an electron hopping from site j to site i, and s z is the z Pauli spin matrix. Next, λ R describes the Rashba interaction 26 that couples differently oriented first-neighbor spins, with s being the vector of Pauli matrices. Finally, λ v is the strength of the staggered on-site potential, for which ξ i is +1 and −1 on A and B sites respectively. Note that the symmetry of the problem is lowered significantly compared to an ideal honeycomb lattice, since the on-site staggered potential makes the A and B sites inequivalent, while the Rashba term breaks s z conservation.
To proceed, we choose the tight-binding basis wavefunctions to be  where σ is a spin index, j = {A, B} denotes the atom type, t j is a vector that specifies the position of the atom in the unit cell, 27 and R is a lattice vector built from the primitive lattice vectors a 1 and a 2 . This allows the Hamiltonian to be written as a 4×4 matrix H jσ,j ′ σ ′ (k) = χ jσk |H|χ j ′ σ ′ k , which can be cast in terms of five Dirac matrices Γ α and their ten commutators where the Dirac matrices are chosen to be Γ 1,2,3,4,5 = (I ⊗ σ x , I ⊗ σ z , s x ⊗ σ y , s y ⊗ σ y , s z ⊗ σ y ) with the Pauli matrices σ k and s k acting in sublattice and spin space respectively. The dependence of the d α and d αβ coefficients on wavevector is detailed in Table I using the notation x = k x a/2 and y = √ 3k y a/2, with the relationship of these variables to the BZ being sketched in Fig. 2.
Since,θΓ αθ−1 = Γ α andθΓ αβθ−1 = −Γ αβ , while d α (k) = d α (−k) and d αβ (k) = −d αβ (−k), the Hamil- tonian (9) is time-reversal invariant, i.e.,θH(k)θ −1 = H(−k). However, it lacks particle-hole symmetry in the sense of Refs. (1-3), because of the action of the on-site and spin-orbit coupling terms. In the general classification of topological insulators and superconductors, 1-3 therefore, the Kane-Mele model falls into the AII symplectic symmetry class, which in two dimensions has a Z 2 classification. This means that by varying parameters of the Hamiltonian of Eq. (9) one can switch between Z 2odd and Z 2 -even phases, with the system experiencing a gap closure and becoming metallic at the transition from one phase to the other.
For the present purposes we assume λ SO > 0 without loss of generality. We also fix λ v > 0. For this case, the transition between Z 2 -odd and Z 2 -even phases is accompanied by a gap closure at the K and K ′ points (the zone-boundary points of three-fold symmetry) in the BZ. The energy is independent of t at these points, and λ SO can be used as the energy scale. The energy gap is then given by |6 Fig. 3.
Note that when λ R = 0 the model reduces to two independent copies of the Haldane model 5 the Z 2 invariant is odd when the Chern numbers are odd, and even otherwise. 28 In what follows we use t as the energy scale and fix the values of the other parameters to be λ SO /t = 0.6 and λ R /t = 0.5. Varying the third parameter λ v /t allows us to switch from the Z 2 -even to the Z 2 -odd phase. The phase transition occurs at |λ v |/t ≃ 2.93, with the system in the Z 2 -odd phase for −2.93 < λ v /t < 2.93. As discussed above, the energy gap closes at the phase transition, and remains open in both the Z 2 -odd and Z 2 -even phases.

A. General considerations
We now consider the problem of constructing Wannier functions (WF) for the Kane-Mele model. We emphasize that we mean by this a set of localized functions spanning the same space as the occupied Bloch bands. Several recent papers have discussed the construction of WFs for an enlarged subspace including also some unoccupied bands for 3D topological insulators such as Bi 2 Se 3 , 8,29 in which case there is typically no topological obstruction, but this is not the context of the present work.
We start with the general definition of the WF in cell R and with band index n in 2D, where A is the unit cell area and Bloch wavefunctions ψ nk are assumed to be normalized within the unit cell. This definition is not unique; not only is there the usual U(1) gauge freedom associated with a k-dependent phase twist of each band n, there is more generally a U(N ) gauge freedom coming from the fact that the N occupied Bloch bands can be mixed with each other by a k-dependent U(N ) transformation. In fact, it is generally necessary to premix the Bloch states using this U(N ) gauge freedom in order that the resulting Bloch-like states (and their phases) will be smooth functions of k. However, having done so, there is still a large gauge freedom associated with the application of a subsequent U(N ) gauge rotation that is smooth in k. This ambiguity in the gauge choice can be removed by applying some criterion to the selection of the WFs. Since electrons are expected to be localized in insulators, 30 a sensible criterion is that of Ref. 17, which specifies maximal localization of the WFs in real space. In this approach, which we adopt here, one chooses some localized trial functions in order to provide a starting guess about where the electrons are localized in the unit cell, and obtains a fairly well-localized set of WFs by a projection procedure to be described shortly. If desired, one can follow this with an iterative procedure to make the resulting WFs optimally localized. 17 Consider an insulator with N occupied bands. We start with a set of N trial states |τ i located in the home unit cell, and at each k we project them onto the occupied subspace at k to get a set of Bloch-like states Since this set of states will not generally be orthonormal, we make use of a Löwdin orthonormalization procedure which consists of constructing the overlap matrix S mn (k) = Υ mk |Υ nk (15) and obtaining the orthonormal set of Bloch-like orbitals Note that theψ nk are not eigenstates of the Hamiltonian, but they span the same space, and have the same form, as the usual Bloch eigenstates. For an insulator whose gap is not too small, and for a set of trial functions embodying a reasonable assumption about character of the localized electrons, theψ nk will be smooth functions of k. In that case, by the usual properties of Fourier transforms, the WFs constructed in analogy with Eq. (12), should be well localized. Such a construction will break down if the determinant of S(k) vanishes at any k. This is guaranteed to occur in a Chern insulator, where time-reversal symmetry is broken and the Chern index of the occupied manifold is non-zero; in this case, construction of exponentially localized WFs becomes impossible. [19][20][21] For a Z 2 insulator, however, the presence of time-reversal symmetry guarantees a zero Chern index, so that exponentially localized WFs must exist. 21 In this case, we should be able to find a set of trial functions such that det S(k) = 0 throughout the BZ.

B. Z2-even phase
Let us first apply the method described above to the case of the Z 2 -even phase of the Kane-Mele model. This phase is topologically equivalent to the ordinary insulator, so we anticipate a picture in which the two electrons per cell are opposite-spin ones approximately localized on the lower-energy (B) site. One way to see this is to look at the weights of the basis states in the occupied subspace. Figure 4(a) shows the distribution of these weights along a high-symmetry line in the BZ for the Kane-Mele model in its Z 2 -even phase. From the figure it is obvious that the two basis states on the B site dominate in the occupied subspace over the whole BZ. It is then natural to choose the two trial functions to be opposite-spin spatial δ-functions localized on the B site in the home unit cell. We choose these to be spin-aligned along z, i.e., where |σ z 1 = | ↑ z and |σ z 2 = | ↓ z . Transforming to k-space we get The two occupied Bloch bands may be written as where ℓ is a combined index for sublattice and spin, ℓ = The overlap matrix S is constructed from these functions, and for the determinant one finds Recall that for the Löwdin orthonormalization procedure to succeed, this determinant must remain non-zero everywhere in the BZ. This is indeed the case for the Z 2even phase, as illustrated in Fig. 5(a), where the solid black curve shows the dependence of the determinant on k along the high-symmetry line in the BZ. In contrast, the dashed red curve in Fig. 5(a) shows the behavior of det [S(k)] in the Z 2 -odd regime. The determinant can be seen to vanish at the K and K ′ points in the BZ. Clearly, this choice of trial functions is not appropriate for building the Wannier representation in the Z 2 -odd phase. Indeed, as we shall see in the next subsection, any choice of trial functions that come in Kramers pairs is guaranteed to fail in the Z 2 -odd case. There we shall also investigate alternative choices of trial functions that allow for a successful construction of WFs.

C. Z2-odd phase
To gain some insight into the appropriate choice of trial functions in the Z 2 -odd regime, consider the weights of the basis functions in the occupied space shown for this case in Fig. 4(b). Unlike the normal insulator, the Z 2 -odd phase does not favor any particular basis states. Instead, different basis states dominate in different portions of the BZ. For example, at points K and K ′ the occupied space is represented by only two of the four basis states; at each of these points the two participating basis states have opposite spin and sublattice indices, and none appear in common at both points. (The states at K are, of course, Kramers pairs of those at K ′ .) It follows that if any of the trial states is simply set equal to one of the four basis states, then at least one of the |Υ would vanish either at K or K ′ , and the determinant would vanish there too. This explains the failure of the naive Wannier construction procedure for the Z 2 -odd phase; with the naive choice of trial functions as in Eq. (18), the determinant vanishes at both K and K ′ , as shown by the red dashed curve in Fig. 5(a). 31 In fact, this failure can be understood from a general point of view. If the two trial functions form a Kramers pair, then the projection procedure of Eqs. (14-16) will result in Bloch-like functions obeying The WFs obtained from Eq. (17) will then also form a Kramers pair. But Eq. (24) is nothing other than the constraint of Eq. (5) defining a gauge that respects timereversal symmetry, and it has been shown 23,32,33 that an odd value of the Z 2 invariant presents an obstruction against constructing such a gauge. In other words, in the Z 2 -odd phase a smooth gauge cannot be fixed by choosing trial functions that are time-reversal pairs of each other, and a choice of WFs as time-reversal pairs is not possible. Hence, in order to construct the Wannier representation in the Z 2 -odd regime, one should choose trial functions that do not transform into one another under time reversal. Following these arguments, we choose the two trial functions to be localized on different sites in the home unit cell. Moreover, in order that they will have components on states with spins both up and down along z, we choose the spins of the trial states so that one is along +x and the other along −x. 34 In k-space this becomes where t 1 = t A and t 2 = t B , leading to The determinant takes the form The dependence det[S(k)] is shown along the diagonal of the Brillouin zone for this choice of trial functions in Fig. 5(b). In the Z 2 -odd phase (dashed line) the determinant remains non-zero everywhere in the BZ. 35 Not surprisingly, the same trial functions are very poorly suited to the normal-insulator phase, as can be seen from solid line in the same panel. In this case det[S(k)] almost vanishes at K and K ′ and remains quite small throughout the rest of the BZ, so that one should clearly revert to the time-reversed pair of trial functions of Eq. (18) and Fig. 5(a) in order to get well-localized WFs.
We made an arbitrary choice above in selecting the two trial functions to be up and down along x. In fact, if we repeat the entire procedure using trial functions that are spin-up and spin-down along any unit vectorn lying in the xy-plane, we find that det[S(k)] changes very little, with only small changes in the size of the dip near the Γ point. Thus, we find that the choice of trial functions in Eq. (25) is not unique. Instead, there is a large degree of arbitrariness in the choice of WFs in the Z 2 -odd case.
To conclude, we have established that the choice of a time-reversal pair of trial functions, Eq. (18), that allows for the construction of well-localized WFs in the ordinary-insulator phase cannot be used in the Z 2 -odd phase. In order for the usual projection method for constructing the Wannier representation to work in this topologically nontrivial phase, the trial functions should explicitly break time-reversal symmetry, i.e., they should not come in time-reversal pairs.

V. LOCALIZATION OF WANNIER FUNCTIONS IN THE Z2-ODD INSULATOR
Now that we know how to construct WFs for the Z 2odd insulator, we discuss their localization properties. As we have noted in the preceding section, the choice of the trial functions, Eq. (25), is not unique; there are other gauge choices arising from different trial functions that also produce well-defined sets of WFs. Since different gauge choices lead to different degrees of localization of the resulting WFs, it is natural to fix the gauge by the condition of maximal possible localization of the WFs.
The problem of constructing maximally-localized WFs was studied by Marzari and Vanderbilt. 17 They considered the total quadratic spread as a measure of the delocalization of WFs in real space, and developed methods for iteratively reducing the spread via a series of unitary transformations, Eq. (13), applied prior to WF construction. The spread functional was decomposed into two parts, Ω = Ω I +Ω, with being the gauge-invariant part and the gauge-dependent part of the spread. Discretized kspace formulas for Eqs. (30) and (31)  for the case that the BZ is represented by a uniform k mesh. The resulting expression for the gauge-invariant spread is, for example, where 33) are overlap matrices and b are "mesh vectors" connecting each k-point to its nearest neighbors. The latter are chosen, together with a set of weights ω b , in such a way as to satisfy the condition A corresponding expression forΩ, and a description of steepest-descent methods capable of minimizing Ω, were also given in Ref. 17. Note that, in order to avoid getting trapped in false local minima, the iterative procedure is normally initialized using the trial-function projection procedure described in Sec. IV above. We now apply this method to the Kane-Mele model. The lattice is hexagonal, and in this case six b j vectors are needed to satisfy the condition (34), namely b 1 = −b 4 = G 1 /q, b 2 = −b 5 = (G 1 + G 2 )/q, and b 3 = −b 6 = G 2 /q. All six have the same length b and weight ω b = 1/(3b 2 ). We start with the WFs obtained with the projection method using the trial functions of Eq. (25), appropriate for the Z 2 -odd phase.
The resulting spreads, both before and after the iterative minimization, are shown in Fig. 6. (Ω I , being gauge-invariant, is the same before and after.) The left part of the figure shows the behavior in the Z 2 -odd phase, where the trail functions are the appropriate ones. The results in this region were not strongly sensitive to the k-point mesh density. The fact thatΩ is similar in magnitude to the unminimized Ω I , and that the localization procedure reducesΩ by only 20−30%, provide additional evidence that the choice of trial functions was a good one. The Wannier charge centers were almost unchanged by the minimization procedure; the x-coordinates were zero, whiler 1y ≃ a/ √ 3 andr 2y ≃ 2a/ √ 3 (see Sec. VI for details), in good agreement with our initial assumption about the WFs being localized on A and B sites.
The right part of Fig. 6, for λ v /t 2.93, shows what happens when we attempt to use the same trial functions in the normal phase. Ω I is of course unaffected by the choice of trial functions, and the fact that it has a smaller value in this region indicates, not surprisingly, that the insulating state is simpler and more localized in the normal state. (For large λ v /t the WFs approach spatial delta functions, explaining the fact that Ω I asymptotes to zero in that limit.) Not surprisingly, however, using the trial functions appropriate to the Z 2 -odd phase in the Z 2 -even regime results in very poor localization of the WFs as measured byΩ. Our data also suggests that in the Z 2 -odd phase MLWFs are less localized than MLWS in the Z 2 -even phase. For example, the use of trial functions (18) with λ v /t = 5 and a 60 × 60 k-mesh results in Ω I = 0.02770 andΩ = 0.00025. We also find that the results are more sensitive to the choice of k-mesh in the Z 2 -odd regime.
To summarize the results of this section, we studied the construction of maximally localized WFs in the Z 2 -odd phase using the Kane-Mele model as an example. We have seen that our initial guess of Sec. IV about the localization of WFs in this topological regime is very good, and that the maximal localization procedure does not greatly reduce the spread.

VI. HYBRID WANNIER CHARGE CENTERS AND POLARIZATION
In this section we discuss the polarization in Z 2 -odd insulators using the example of the Kane-Mele model, and see what insights about the topological insulating phase can be obtained by inspecting this property.
The electronic polarization in a 2D system can be defined either in terms of the Berry phase 36 or via the summation of Wannier charge centers 18 where e is the electronic charge and A is the area of the unit cell. The two definitions are identical and define electronic polarization modulo a polarization quantum |e|R/A, R being a lattice vector. This ambiguity can be understood as a freedom in the choice of branch in Eq. (35) or in the choice of unit cell in Eq. (36). The definition via Wannier charge centers makes the dependence of P on the choice of origin obvious. As described in Sec. III, the origin of the Kane-Mele model is chosen such that atoms are located along the y-axis at t A = ξŷ/3 and t B = 2ξŷ/3, where ξ = |a 1 +a 2 | = a √ 3. Because the Hamiltonian has 3-fold symmetry, we expect the rescaled polarization (A/|e|)P to lie at the origin, at t A , or at t B . To distinguish between these possibilities it is sufficient to compute P y , which is well-defined modulo |e|/a.

A. Total polarization
A direct computation of electronic polarization via Eq. (35) in the Z 2 -even phase results in P y = |e|/3a mod |e|/a, consistent with the fact that both Wannier centers in Eq. (36) lie at t B (since −4|e|ξ/3A = −8|e|/3a = |e|/3a mod |e|/a.) In the Z 2 -odd phase, on the other hand, Eqs. (35) and (36) lead to P y = 0 mod |e|/a. Again, this is consistent with the locations of the WFs. As indicated in Sec. V, the Wannier centersr n in this phase lie approximately at t A and t B . More precisely, we find that they are located atr 1 = (1 − δ)ξŷ/3 and r 2 = (2 + δ)ξŷ/3, where δ is a small correction (e.g., δ = 0.0018 at λ v /t = 1). Thus, the sum of the Wannier centers is just ξŷ, or zero modulo a lattice vector.
It is interesting to note that, in retrospect, the computation of the polarization via Eq. (35) would have given a strong hint about the appropriate choice of trial functions in the Z 2 -odd insulator. That is, knowing only that P y = 0, one might have guessed that both WFs should be centered halfway between t A and t B , or both at the center of the honeycomb ring, or one at t A and the other at t B . The latter possibility becomes the most likely when we also take into account that in the Z 2 -odd phase the two WFs cannot form a Kramers pair.

B. Hybrid Wannier decomposition
In order to obtain a deeper understanding of the origin of the polarization and expose some qualitative differences in the behavior of its k-dependent decomposition in Z 2 -even and odd phases, it is useful to use a hybrid representation in which the Wannier transformation is carried out in one direction only. As indicated above, we know from symmetry considerations that we can set P x = 0 and characterize the polarization by P y mod ξ|e|/A. To compute P y , it is convenient to choose the BZ to be a rectangle extending over k x ∈ [0, 2π/a] and k y ∈ [0, 4π/ξ] (corresponding to the region ζ in Fig. 2). We can then define hybrid WFs 4π/ξ 0 dk y e −ikyly |ψ nk (37) in terms of which the usual WFs are The hybrid Wannier centers are defined as and the total electronic polarization is In practice the k x integral is discretized by a sum over a mesh of k x values, and at each k x theȳ n (k x ) are calculated by considering the corresponding string of k-points along k y . In the case that the gauge has been specified by a particular set of 2D WFs |Rn , or, equivalently, by the corresponding Bloch-like functions |ψ nk , this is done straightforwardly using the discretized Berry-phase formulaȳ where M (j) is a shorthand for the overlap matrix M (kj,kj+1) of Eq. (33) connecting k y -points j and j + 1 along the string. As was emphasized in Sec. IV, theψ nk carry the information about the gauge choice. Thus, different gauge choices -i.e., different choices of WFs -will result in different hybrid WFs and differentȳ n (k x ). However, the sum nȳ n (k x ) at a given k x is gauge-invariant, and as a result P y of Eq. (40) must remain the same in any gauge.
Of special interest is a gauge choice in which, at each k x , the hybrid WFs |nk x l y are maximally localized in the y direction. It was shown in Ref. 17 that in 1D the Wannier charge centers can be obtained by a paralleltransport construction using the overlap matrices M (j) . Specifically, the "unitary part" M (j) of each overlap matrix is obtained by carrying out the singular-value decomposition M = V ΣW † , where V and W are unitary and Σ is real-positive and diagonal, and then setting M = V W † . This is reasonable because, for a sufficiently fine mesh spacing, Σ is almost the unit matrix. Then, the unitary matrix Λ = j M (j) describes the transport of states along the string. The eigenvalues λ n of this matrix are all of unit modulus, and their phases define Wannier centers via 37ȳ Im log λ n .
Note that no iterative procedure is needed. Inserting this equation into Eq. (40), one gets a discretized formula for P y that is consistent with Eq. (35).

C. Results
We illustrate these ideas now for the KM model in its normal and Z 2 -odd phases. In each case we present results forȳ n (k x ) for two choices of gauge: the maximallylocalized one alongŷ as discussed in the previous paragraph, and the one corresponding to the WFs constructed from the trial functions of Eq. (18) for the Z 2 -even phase or those of Eq. (25) for the Z 2 -odd phase. In what follows, we refer to these as the "maxloc" and "WF-based" gauges respectively.
In the ordinary insulating regime, the maxloc and WF-basedȳ n (k x ) curves look very similar to each other. Fig. 7(a) and (b) show the calculated results for the case of λ v /t = 3, very close to the transition on the insulating side (recall the critical value is at λ v /t = 2.93). Three of the infinite number of periodic images along y are shown. The "bumps" in the curves near the K and K ′ points in the BZ are the result of the proximity to the transition; as one goes deeper into the insulating phase, the curves flatten out and become smooth functions of k x . The solid and dashed curves are mirror images of each other; in the maxloc construction of Fig. 7(a) this just reflects the time-reversal invariance of the Hamiltonian, while in Fig. 7(b) it follows from the fact that the WFs form a Kramers pair.
When averaged over k x , each curve is found to have a meanȳ value of 2ξ/3 to numerical precision, or ξ/6 modulo ξ/2, consistent with the discussion in Sec. VI A.
The corresponding results for the Z 2 -odd phase are shown in Fig. 7(c) and (d) for λ v /t = 1. As expected, there is again a mirror symmetry visible in the curves for the maxloc construction in Fig. 7(c), but the connectivity of the curves is qualitatively different: in going from k x = 0 to π/a we see that the n'th solid curve goes up to cross the (n + 1)'th dashed curve, while the n'th dashed curve goes down to cross the (n − 1)'th solid curve. This is exactly the kind of behavior that was exhibited in Fig. 3(a) of Ref. 23 as a signal of the Z 2 -odd phase. Moreover, if we follow the n'th dashed curve all the way across the BZ, we find that it wraps to become the (n + 1)'th one when k x = 2π/a wraps back to k x = 0. This is precisely the kind of behavior that is characteristic of a Chern (or quantum anomalous Hall) insulator, 38 which implies that we can assign a Chern number of +1 to the Bloch subspace spanned by the eigenvectors corresponding to the dashed bands. However, since we are studying here a system with time-reversal symmetry, we find also a partner subspace corresponding to the full curve in Fig. 7(c) having Chern number −1. As a result, of course, the overall occupied space has a total vanishing Chern number, as it must due to the time-reversal symmetry. The evaluation of the polarization P y through Eq. (40) again yields P y = 0 mod |e|/a, consistent with the direct calculation of Sec. VI A.
Finally, Fig. 7(d) shows theȳ n (k x ) curves for the same Z 2 -odd parameters as in Fig. 7(c), but using the WFbased gauge determined by the trial functions of Eq. (25). At any given k x , we confirm thatȳ 1 +ȳ 2 is the same in Fig. 7(d) as in Fig. 7(c), and the total polarization is therefore the same. However, because the two WFs do not form a Kramers pair in this case, the dashed and solid curves do not map into each other under time-reversal symmetry, and there is no degeneracy at k x = π/a. Moreover, the Chern number of each band is individually zero, consistent with the fact that each one is derived from a WF. The averageȳ values for the solid and dashed curves are 0.978ξ/3 and 2.022ξ/3 mod ξ/2, very close to the nominal locations of the trial functions at t A and t B , respectively.
To recap, in both the Z 2 -even and Z 2 -odd cases, we find that the occupied Bloch space can be cast as the direct sum of two subspaces that map into one another under the time-reversal operation, corresponding to the solid and dashed curves of Figs. 7(a-c). These subspaces are not built from Hamiltonian eigenstates, but from suitable k-dependent U(2) rotations among the Hamiltonian eigenstates. In the Z 2 -even case the Chern index of each of these subspaces is separately zero, so that we can also provide a Wannier representation for each subspace separately. This is essentially the case of Fig. 7(b), and since the spaces form a time-reversal pair, the WFs form a time-reversal pair as well. In contrast, for the Z 2 -odd phase, the decomposition into two subspaces that are time-reversal images of each other necessarily results in subspaces having individual Chern numbers of ±1, and these are not individually Wannier-representable. Only by violating the condition that the two spaces be timereversal partners, as was done in Fig. 7(d), can we decompose the space into two subspaces having zero Chern indices individually. By doing so, we can find a Wannier representation of the entire space, but only on condition that the two WFs do not form a Kramers pair.

VII. CONCLUSIONS
In this paper we have considered the question of how to construct a Wannier representation for Z 2 -odd topological insulators in 2D. We have shown that the usual method based on projection onto trial functions fails because of a topological obstruction if one imposes the condition that the trial functions should come in timereversal pairs. On the other hand, the projection method can be made to work if this condition is not imposed, resulting in WFs that do not transform into one another under time reversal.
Such a Wannier representation may have some formal disadvantages. For example, if one writes the Hamiltonian as a matrix in this Wannier representation, its time-reversal invariance is no longer transparent, and the presence of other symmetries may become less obvious as well. On the other hand, it does satisfy all the usual properties of a Wannier representation, as for example the ability to express the electric polarization in terms of the locations of the Wannier centers, and there is every reason to expect that the maximally localized WFs are still exponentially localized. 21 The generalization of our findings to the 3D case should be relatively straightforward. Certainly the topological obstruction to the construction of Kramers-pair WFs re-mains for both weak and strong Z 2 topological insulators in 3D. To see this, consider in turn each of the six symmetry planes in k-space (k 1 = 0, k 2 = 0, k 3 = 0, k 1 = π/a, etc.) on which H k behaves like a 2D time-reversal invariant system. For both weak and strong topological insulators, at least one of these six planes must have a Z 2 -odd 2D invariant. But if a gauge exists obeying the time-reversal condition of Eq. (5) in the 3D k-space, then it does so in particular on the 2D plane, in contradiction with the 2D arguments about a topological construction.
Thus, the general strategy for constructing WFs for 3D topological insulators should be very similar to the one presented here in 2D. Namely, one has to choose pairs of trial functions that do not transform into one another by time-reversal symmetry, and to do it in such a way that the projection of these trial functions onto the Bloch states does not become singular anywhere in the 3D BZ. While it may be interesting to explore how this might best be done in practice for real 3D topological insulators, e.g., in the density-functional context, the choice is likely to depend sensitively on details of the particular system of interest. Thus, an investigation of these issues falls beyond the scope of the present work.

VIII. ACKNOWLEDGMENTS
The work was supported by NSF Grant DMR-0549198.