Electric polarization as a nonquantized topological response and boundary Luttinger theorem

We develop a nonperturbative approach to the bulk polarization of crystalline electric insulators in $d\geq1$ dimensions. Formally, we define polarization via the response to background fluxes of both charge and lattice translation symmetries. In this approach, the bulk polarization is related to properties of magnetic monopoles under translation symmetries. Specifically, in $2d$ the monopole is a source of $2\pi$-flux, and the polarization is determined by the crystal momentum of the $2\pi$-flux. In $3d$ the polarization is determined by the projective representation of translation symmetries on Dirac monopoles. Our approach also leads to a concrete scheme to calculate polarization in $2d$, which in principle can be applied even to strongly interacting systems. For open boundary condition, the bulk polarization leads to an altered `boundary' Luttinger theorem (constraining the Fermi surface of surface states) and also to modified Lieb-Schultz-Mattis theorems on the boundary, which we derive.


I. INTRODUCTION
The bulk electric polarization of an insulator is a concept of fundamental importance in condensed matter physics [1][2][3][4][5][6][7][8][9][10]. Powerful methods have been established to calculate the polarization, at least within the independent electron approximation of band theory [1,2,[9][10][11][12]. The precise definition and interpretation of polarization, however, is a subtle issue (see, for example, Refs. [3,13]). Intuitively the polarization measures the density of electric dipole moment in the bulk, but precisely defining this notion is known to be nontrivial. The physics is most easily understood in one dimension, where the polarization can be defined by a topological Θ-term in the low energy, long wavelength (IR) electromagnetic response (setting = c = e = 1): which can be motivated by the fact that a dipole moment d couples to electric field as −d · E. Since dA = dxdt (∂ t A x −∂ x A t ) is always an integer (the first Chern number) multiple of 2π on a closed 2-manifold, P is defined mod 1. This periodicity of P can be understood on the lattice by noticing that shifting an integer charge by one lattice unit (we set to be a = 1) in every unit cell is equivalent to a relabeling of lattice coordinates and should not have any physical effect.
The polarization, defined via Eq. (1), has several consequences. For a system with periodic boundary condition, consider an adiabatic flux-threading process where dx A x changes by 2π, followed by a large gauge transform that brings the Hamiltonian back to the original form. This process corresponds to dA = 2π in the space-time path integral, and Eq. (1) leads to a Berry phase Φ = 2πP . This Berry phase is in principle a measurable quantity and is sometimes used as the definition of polarization in one dimension [14,15]. For open boundary condition Eq. (1) becomes a boundary term ±P dtA t , which represents a fractional charge q = ±P (mod 1) at each boundary -the mod 1 condition comes from the fact that one can always deposit an integer charge on the boundary without affecting the bulk. This is consistent with the intuitive connection between polarization and dipole moment.
Raising this to higher dimensions poses some challenges. A simple generalization of the electromagnetic response term Eq. (1) to higher dimensions does not produce a topological term. One can consider it as a Berry phase term, and measure the polarization through the Berry phase of a flux-threading process (say in the x-direction) similar to that in 1D. The Berry phase, however, is given by where V = L x L y ... is the system volume and L x is the length in x-direction. If we assume lattice translation symmetries (which we do for the rest of the paper), the intensive quantity P can be extracted from the Ldependence of Φ (but simply dividing by V /L x will not work since the phase is defined mod 2π). However it raises the conceptual question whether P itself bears any physical meaning. For example, for a two-dimensional crystalline insulator with L y = 2N (N → ∞ in thermodynamic limit) and a polarization density P x = 1/2, the Berry phase from Eq. (2) is always trivial. Is there a formula for the polarization in this case? One can always define polarization by starting begin with a reference state with a known polarization (for example where polarization is constrained by symmetry) and connect this to the Hamiltonian of interest by an adiabatic path, and integrating the currents obtained while connecting initial and final states. However this algorithm requires defining such an adiabatic path and is conceptually different from a direct measure of polarization that we seek.
A similar issue appears with open boundary condition: the density of dipole moment, which is the classical definition of polarization, is given by the surface charge density. Unlike the 1D case where the boundary charge is robustly determined mod 1, the surface charge density in higher dimensions can be continuously tuned by boundary perturbations (for example a boundary chemical potential). It then appears that the boundary does not necessarily reflects the bulk polarization. An exception was observed in Ref. [12]: when the boundary is gapped and non-degenerate, the boundary charge density faithfully represents the bulk polarization mod 1.
In this work we develop a more topological approach to bulk polarization in higher dimensions. Specifically we use a topological term similar to Eq. (1) to define polarization density in any dimension, and explore its consequences for both periodic and open boundary conditions. For periodic boundary condition, the polarization determines the properties of magnetic monopoles under translation symmetries (such as momenta). For open boundary condition, the bulk polarization determines the degree of Luttinger theorem violation on the boundary, and more generally is related to quantum anomalies of the boundary low-energy theory. Our approach applies to systems of interacting fermions as well as bosons (or spins) as long as there is a conserved U (1) charge with a gap to charged excitations (an insulator) and the systems have lattice translation symmetries.
As we shall see, a topological approach is needed because polarization cannot be measured by local probessomething global has to be introduced such as symmetry fluxes or physical boundaries. This justifies the use of the term "topological response", even though the response itself (the polarization) is in general not quantized and hence its value will change in response to symmetric perturbations. The familiar electromagnetic polarizability (the Θ-term) in 3d also falls into this category when time-reversal and mirror symmetries (or more generally symmetries that invert an odd number of spacetime coordinates) are absent.
Before discussing the main results, we note that there is a feature in the definition of polarization when the unit cell has a nontrivial geometric structure. The polarization P consists of two parts: a classical electric dipole moment within each unit cell and an extra part that measures inter-unit-cell entanglement [16,17] -the latter is denoted asP in Ref. [13]. Both P andP have been discussed in the literature depending on context. We review these notions briefly in Appendix A. Our results below can be applied to both the full P and toP as long as we adopt the appropriate calculation scheme as described below.

II. POLARIZATION FROM TOPOLOGICAL TERMS
Interesting low energy and long wavelength (IR) physics can often be probed by the response to background gauge fields. In the study of polarization density in d > 1, the relevant symmetries include charge conservation and lattice translation symmetries, so we shall consider coupling the system to gauge fields associated with these symmetries. For charge conservation the gauge field is simply the electromagnetic field A µ . For each translation symmetry Z, say in the i'th direction, we introduce a Z-gauge field x i . This "translation gauge field" [18] is less familiar so we review below. The gauge field x i is locally flat (dx i = 0) so only its Wilson loops C1 x i ∈ Z over loops (or 1-cycles C 1 ) in space-time is meaningful -formally this means that where M is the space-time manifold. Furthermore, just like the Wilson loops in other gauge theories, the integer C1 x i measures the number of x i -translations one has to go through to travel across C 1 . To be more concrete consider a path integral description, with dynamical degrees of freedom ψ (bosonic or fermionic) defined in continuous time t ∈ [0, T ) and on discrete lattice sites s in space: where we have used locality and translation symmetries to write the Lagrangian as a sum of local terms of identical form, L s [ψ, A], which involves only fields near site s. We take periodic boundary conditions in space and time (so M is a torus). The translation gauge fields enter the partition function by specifying exactly how the periodic boundary contitions are taken: We now explain these equations in more details. The Wilson loop of x i in the x i direction gives the lattice size i x i = L i . For j = i the number i x j measures how much the slice of the lattice at x i = L i is displaced along the x j direction before it is identified with the slice at x i = 0. Similarly the time component t x i measures the displacement of the entire lattice at t = T before identified with t = 0. In another word, while the "longitudinal" parts of the translation gauge fields measure the lattice size, the "transverse" parts measure the quantized shear strains of the lattice in both space and time. We can also consider a (d − 2) dimensional defect in space, around which x i = n = 0: this is simply a lattice dislocation with Burgers vector B = n x i . The translation gauge field x i is closely related to the concept of tetrad in the theory of elasticity [19], which has been used to characterize three dimensional integer quantum Hall effect recently [20]. Consider embedding the lattice into a continuous space, so that each site s can be assigned a continuous coordinate u s . We can treat u as a field, then the tetrad ∇u i will have all the properties of x i discussed above, and can be used as a representation (a gauge choice) of x i . The gauge invariant properties of x i such as the Wilson lines, however, do not depend on how the lattice is embedded into a continuous space. In this sense the x i gauge field measures the topological part of elasticity response. Now recall that the electric polarization in 1d can be defined through the topological term P dA (Eq. (1)). The natural generalization to higher (d + 1) dimensions is the following term: Here ∧ should really mean cup product for discrete cohomology instead of the usual wedge product, but the distinction does not matter for our purpose.
We now give some justifications for Eq. (5) as a definition of bulk polarization. First, it is the only topological term involving dA and x i that is first order in the field strength dA, as we expect for the polarization. Each component of polarization P i is defined mod 1 since the integral always gives integral multiples of 2π on closed manifolds (the term is therefore a topological Θ-term), and is in agreement with the intuition that shifting integer charges by one lattice unit does not have physical effect. When evaluated for a uniform electric field E on a perfect lattice (free of dislocation and shear strain) of size L 1 × L 2 ..., this term becomes which agrees with the expectation that the total polarization is (V /L i )P i when the system is viewed as 1d in x i direction. In addition, if P has a time-dependence P(t), then by taking derivative with respect to A from the above action we obtain the charge current j = ∂P/∂t, which agrees with physical expectations and is sometimes used as a practical way to define polarization.

III. POLARIZATION AND MONOPOLES
We now examine the consequences of the polarization as defined through Eq. (5) on a closed manifold (like periodic boundary condition). Motivated by the 1d case, it is useful to consider instantons of the A field. In 1d the instanton is the familiar adiabatic flux-threading, a smooth configuration in spacetime. In higher dimensions the instantons become operators supported on (d−2) dimensional sub-manifolds in space, with dA = 2π on the two complementary spatial dimensions. For d = 2 this is simply a unit flux insertion in space, and for d = 3 it corresponds to a unit flux tube in space whose open ends become Dirac monopoles. On the (d − 2) + 1 manifold of the instanton, the topological terms reduces to the following Dijkgraaf-Witten [21] type: Let us look at some physically relevant examples. At d = 2 we obtain which means that the 2d monopole -a point operator in space -carries "charge" of 2π(−P 2 , P 1 ) under translation symmetries in x 1 and x 2 , respectively. But "charge" under translation symmetry is simply the crystal momentum. We then conclude that in 2d the monopole carries lattice momentum It may be helpful to have some simple semi-classical picture here. Consider a 2π-flux quanta spread uniformly over a region much larger than the lattice unit. We can then consider the momentum of such a monopole k M , i.e. the Berry phase from moving the unit flux configuration by one lattice unit. Equivalently we can consider the many-body momentum of the fermions k e under the flux configuration, which would be the inverse of the monopole momentum. Now imagine a semi-classical continuum system, a non-uniform electric field is induced during the turning-on of the magnetic flux, which then induces a momentum on a small electric dipole moment Since the dipole density is given by P, we have which is what we obtained from the topological term. We can use the monopole momentum as a practical way to calculate electric polarization in 2d. We outline the calculation scheme here and discuss more details in Appendix B. Consider a 2d crystalline insulator on a torus, and smoothly spread a total magnetic flux of 2π on the lattice -say 2π/L x L y flux per plaquette. The total momentum of the many-electron system k e can be measured from the ground state wave function, as Berry phase factors associated with lattice translations, and from Eq. (10) we have k e = 2πP × z. The virtue of this calculational scheme is that it is well defined (although possibly complicated) even for strongly correlated systems, where band theory techniques cannot be used.
Another consequence of the connection between polarization and monopole momentum is that we can now define polarization in 2d even in the presence of gapless Dirac fermions. The only subtlety is that with a unit flux in space, each Dirac cone contributes a zero-energy mode, leading to multiple degenerate ground states depending on which zero modes are occupied. For each of the degenerate monopoles we can nevertheless define its lattice momentum and interpret it as a bulk polarization, which then also depends on the zero mode fillings. In Fig. 1 we report a numerical calculation of the polarization using the monopole momentum, for a lattice system of gapless Dirac fermions with a specific choice of zero mode filling. The same polarization can also be calculated using the standard method from band theory which we also report in Fig. 1. The two results clearly agree as we vary a parameter t in the Hamiltonian. More details of the calculational recipe and the lattice model can be found in Appendix B. In fact, there is a long history of numerically calculating monopole momenta for lattice Dirac fermions, motivated by the study of monopole operators in Dirac spin liquids [22][23][24][25]. We showed that what this calculation really produces is the polarization of the underlying Dirac fermions. For d = 3 the term in Eq. (7) becomes a twodimensional integral which describes a 2π flux loop decorated with a 1d topological phase enriched by translation symmetries. This leads to nontrivial boundary modes when the flux loop has open ends, which are nothing but Dirac monopoles. The boundary mode is characterized by a projective representation [26][27][28] of translation symmetries, namely translations in different directions commute up to a phase when acting on a Dirac monopole: EM Dual e mB = 2πD The semi-classical picture relates projective representation of monopoles to polarization by electro-magnetic duality that exchanges magnetic monopole and electric charge. The electric displacement field (left) D = P is mapped to magnetic field (right)B = 2πD and the monopole to an electric charge. The Aharonom-Bohm (AB) phase θ AB seen by the electric charge is proportional to magnetic field(right), hence the monopole Berry phase proportional to displacement field(left).
This also has a simple semi-classical picture shown in Fig. 2. Consider a 3d continuum system with polarization density P. The polarization leads to an electric displacement field D = P. A magnetic monopole sees the D field as an effective "dual magnetic field" B = 2πD = 2πP. A particle moving in an effective magnetic field realizes translation symmetries projectively, namely different translation operations commute up to a phase factor according to Eq. (12). Similar to the 2d case, the relation Eq. (12) is relevant for U (1) quantum spin liquids in three dimensions, described at low energy by an emergent Maxwell U (1) gauge theory that is potentially realized, for example, in quantum spin ice materials [29]. Our results indicate that the monopoles in a U (1) spin liquid will carry projective translation quantum numbers if the emergent electric charges form an insulator with nontrivial polarization density.
It was known from earlier approaches that polarization is related to other topological quantities including Hall conductance and magneto-polarizability (the axion Θ-angle). In Appendix C we show that these connections can be understood easily using the polarizationmonopole connection.
We summarize the connection between bulk polarization and monopole (instanton) properties in d = 1, 2, 3 in Table. I.

IV. BOUNDARY LUTTINGER THEOREM AND ANOMALY
We now explore the consequences of bulk polarization for open boundaries. Consider a boundary at separating the vacuum at x 1 > 0 and the polarized bulk at x 1 < 0, which preserves all translation symmetries except the one along x 1 . The Θ-term in Eq. (5) becomes a boundary term The meaning of this term can be seen by taking functional derivative with A 0 : it simply means a (fractional) charge density of ρ ∂ = P 1 on the boundary. If the boundary has trivial dynamics in the IR, namely with a unique gapped ground state, then Eq. (13) is the only nontrivial term in the IR description of the boundary. This is the well known statement that for a trivially insulating boundary, the charge density is given by the bulk polarization density mod 1. [12] In general, depending on details at the boundary, the boundary can also host nontrivial low energy degrees of freedom. Let us first consider the simplest scenario: a Fermi liquid metal on the boundary. In this case the boundary charge density ρ is obviously not fixed by bulk polarization P (in direction perpendicular to the boundary) since it can be continuously tuned by perturbations that live only on the boundary. But we expect the Fermi surface volume V F to be tuned simultaneously with the charge density following ∆V F /(2π) d−1 = ∆ρ from Luttinger theorem. We therefore expect where n is the normal vector of the boundary.From this relation the polarization density P can be viewed as a source of Luttinger theorem violation on the boundary. Alternatively, we can view V F as a "quantum correction" to the classical expectation of P = ρ. The "boundary Luttinger theorem" Eq. (14) can be understood using an anomaly-matching argument. It is useful to first phrase the usual Luttinger theorem in terms of anomaly matching, following ideas similar to those in Ref. [30]. Consider a theory in (d − 1) space dimensions of low-energy fermions near a Fermi surface, where fermion modes far away from the Fermi surface have been integrated out already. We couple the background gauge fields A and x i minimally to these fermions and denote the action as S F S [ψ, A, x i ]. It is known that under a large gauge transform, in which the real space Wilson loop along the x i direction Ci A changes by 2π, the total crystal momentum of these low energy fermions changes by [30] To see Eq. (15) from Eq. (16), simply recall that Ci x i = L i and that the total momentum along x i is the coefficient of dtx i . Eq. (16) is also related to the familiar chiral anomaly in (1 + 1) dimension, which we briefly explain in Appendix D. Now for a purely (d − 1) dimensional system that is not the boundary of another space, we should add a background term so that the full theory is gauge invariant. The meaning of the counter term, as we discussed under Eq. (13), is simply a charge density of ρ = V F /(2π) d−1 -this is nothing but the familiar Luttinger theorem! It is now straightforward to extend to the case of Eq. (14). Consider a Fermi liquid on the (d − 1) dimensional boundary of a d dimensional bulk, with Fermi volume V F and boundary charge density ρ. The "surface" theory reads The We should therefore view the last term as a polarization term in d space dimensions, hence Eq. (14). The lesson is that bulk polarization does not directly give a boundary charge density, rather it leads to a boundary quantum anomaly. In Appendix F we numerically study a free fermion model on square lattice and verify that Eq. (14) is always satisfied across a range of parameters with qualitatively different edge behaviors. Eq. (14) also applies with P replaced byP if we also replace the bound charge density by the excess charge density.
When viewed as an anomaly-matching condition, Eq. (14) can also be applied to surface states other than Fermi liquids -we simply need to replace V F by the appropriate anomaly indicators of the low energy theories.
Namely we demand ρ = n A + P mod 1 where n A is the anomaly indicator of the low energy effective theory. For example, for rational values of ρ − P the anomaly can be matched by a gapped ground state with intrinsic topological order, which typically hosts nontrivial quasiparticles with fractional electric charge. In such states the anomaly is encoded in how the topologically nontrivial excitations (like anyons in 2d and flux-loops in 3d) transform under translation symmetries. These anomalies are closely related to Lieb-Schultz-Mattis (LSM) type of theorems that constrain the possible low energy theories of a given lattice system [18,[30][31][32][33][34]. We briefly describe the LSM-type of anomaly for topological orders in two and three dimensions in Appendix E. These results reduce to the previously obtained boundary Luttinger relations in the absence of polarization, as described in Ref. [35].
The logic we used to study the boundary can also be used to study dislocations. A dislocation has space dimension (d−2), and therefore can preserve at most (d−2) translation symmetries. For simplicity we consider a dislocation with Burgers vector B = x 2 and unbroken translation symmetries Z x3 × ...Z x d . The polarization term Eq. (5) reduces to the following on the dislocation (with space-time dimension d − 1): This has the same form as the boundary term Eq. (13), only in one dimension lower. This term leads to the same Luttinger theorem violation as Eq. (14) on the (d − 2) dimensional dislocation, where V F is again interpreted as the anomaly indicator of the low energy effective theory. A special case is d = 2 which has been discussed in Ref. [36], where the V F term is not needed and the polarization directly determines the fractional electric charge nucleated at the dislocation point.

V. SUMMARY
In this work we proposed an nonperturbative definition of the physically measurable polarization density in a crystalline insulator through translation properties of test magnetic monopoles. Our formalism is applicable in any space dimension to systems of interacting electrons but equally to interacting bosons or spins that enjoy a U (1) symmetry, as long as there is a unique gapped ground state. The central result is a response involving background U (1) fluxes and translation fields, captured by a topological term. This response is topological despite the fact that the coefficient which is identified with the polarization is not quantized. Indeed, to probe this response one needs to implement a global (non-perturbative) change, e.g. a U (1) flux (monopole/instanton), a lattice shear, a dislocation or a physical boundary. This surprising connection seems natural in light of the necessity of charge quantization in order to properly define polarization [12], which is the consequence of a compact U (1) symmetry group, which, by the Dirac quantization argument, is related to the existence of magnetic monopole operators. The subtleties in previous literature associated with defining polarization (P vsP) are neatly accounted for by gauge field configurations that have different distributions within a unit cell, but correspond to the same instanton in the continuum limit. Besides given a recipe to obtain the polarization in numerical calculations, the connection bears conceptual significance to boundary physics. We see that the classical relation between polarization and boundary charge density receives a quantum correction in the form of anomaly associated with the boundary low energy theory. For a boundary Fermi liquid the anomaly is associated with the familiar Luttinger theorem, and for a more general boundary phase, it is associated with a Lieb-Schultz-Mattis like theorem, but for a general filling.
Here we discuss different definitions of polarization to emphasize the lattice point of view. Most of the physics below are discussed in Refs. [13,17].
To illustrate the essential point, it suffices to consider a simple 1d lattice with two sub-lattices and one electron orbital on each -generalizations to more complicated unit cells (or even higher dimensions) will be straightforward. We label the unit cells by i ∈ Z and sub-lattices by {a, b}. We consider a simple insulator with unit charge occupation per unit cell, namely n i,a + n i,b = 1, with an extra unit negative (ion) charge q ion = −1 sitting on site a to make the entire system charge-neutral.
The electric polarization is defined via the response to a smooth gauge field (ϕ, A) defined on the latticebut what does "smooth" mean? Clearly we want the gauge field to be slowly varying when moving from one unit cell to another, which means the lattice momentum of the gauge field is close to zero. For example A i,a;i,b ≈ A i+1,a;i+1,b and ϕ i,a ≈ ϕ i+1,a . However this does not uniquely specify how the gauge field should be distributed within a unit cell, namely how A i,a;i,b should be compared with A i,b;i+1,a , or how ϕ i,a should be compared with ϕ i,b .
Apparently we have a choice to make here. One simple choice is to demand A i,a;i,b = 0 and ϕ i,a = ϕ i,b . The continuum limit of the gauge fields (now a smooth function of the continuum coordinate x with no explicit dependence on the sublattice index) will be A(x = i) = A i,b;i+1,a and ϕ(x = i) = ϕ i,a . This is equivalent to viewing the entire unit cell as a single point in space. The polarization defined via the response to such field configurations effectively measures only the inter-unit cell entanglement, and is calledP in Ref. [13].
We can make a more general choice as follows: we demand for some constant α. The continuum limit is then The physical meaning is also clear: we interpret the two sublattices a, b as separated by a distance α in real space (with lattice unit normalized to unity). If the lattice system originates from a continuum with the two sites separated physically by a distance α, then this choice corresponds to a physically uniform electric field, and therefore produces the polarization P for a uniform electric field. To obtain the polarization difference P −P , consider the lattice-scale action difference S[A µ P ]−S[A μ P ], where A µ P , A μ P are the corresponding lattice gauge field distributions with the same continuum limit A µ (x, t). Using the definition of A µ P in Eq. (A1),(A2) and the fact that j i,a;i,b − j i,b;i+1,a = ∂ t n i,b where j is the lattice current operator that couples to A field, one can see that which simply means that P −P = α n i,b , and can be interpreted as a classical contribution, namely a charge q b sitting at site b contributes a dipole moment αq b . Different choices of the intra-cell gauge field distribution also lead to different definitions of charge density (per unit cell) in the continuum limit. In the continuum limit we define charge density as ρ(x) = δS/δϕ(x) where ϕ(x) is in the continuum limit as discussed above. This is easy for the gauge field probingP, where ϕ is unique within a unit cell, and we simply get ρ i = n i,a + n i,b + q ion . However with nontrivial α (as for the standard P), there is a correction from the nonuniformity of ϕ within the unit cell. A simple calculation gives δρ i = −α∇ n i,b , which is nonzero only at the boundary. The relation ρ = −P (mod 1) at the boundary (an insulating boundary in d > 1)holds for both P andP as long as the corresponding definitions for charge density are used. A similar difference between current operators from j(x) = δS/δA(x) in different continuum limits was also discussed in the literature. [13] The above discussion implies that when one considers small momentum (long wavelength) components of dA, e.g. a uniform electric field induced by monopole translation, Eq. (5) gives a response controlled by the conventional polarization P; while dA around higher momenta n i G i (n i ∈ Z) probes sites inside a unit cell with different weights. In one extreme case, where one concentrates inductive electric field on only inter-cell bond, the response term givesP.
Appendix B: Details of 2d calculation

Recipes for calculation
Consider a square lattice with L x × L y unit cells, and assume a unique gapped ground state. We put a total magnetic flux of 2π uniformly on the entire surface. To be concrete, let us take the following gauge (analogue of Landau gauge on a discrete torus): In this gauge a unit translation in x (denoted T x ) should be followed by a gauge transform that acts nontrivially only on the x = 0 strip: where q i is the charge density operator on site i. The y-translation (denoted T y ), in contrast, does not need an additional gauge transform. Strictly speaking, however, on a finite torus neither G x T x nor T y is a true symmetry since the Wilson loop along the nontrivial y and x cycles cannot be translationally invariant under T x and T y , respectively. As we can see explicitly from Eq. (B1), the Wilson loops changes by dyδA y = −2π/L x on every y-cycle and dxδA x = 2π/L y after T x and T y , respectively. This non-invariance of Wilson loops cannot be cured by a gauge transform. To overcome this issue, we consider modified translationsT x = F y G x T x andT y = F x T y , where F y , F x are adiabatic evolutions that modify the A fields by δA at the end of the evolutions, where F y : δA x = 0, δA y = 2π L x L y , The composite operationsT x ,T y preserve the Hamiltonian, and therefore produce well-defined Berry phases φ x , φ y which we identify with k e = −k M = 2π(P y , −P x ). The connection between the monopole momentum and polarization can also be understood from the structure ofT x ,T y . Consider the operationsT Lx x andT Ly y . Using T Lx x = T Ly y = 1, one can see that the two operations become the familiar adiabatic 2π flux threading in the y and − x directions, respectively. The corresponding Berry phases are (Φ x , Φ y ) = 2π(L x P y , −L y P x ), in agreement with our previous result. Notice that sincẽ T Li i = 1, the translation Berry phase defined above is not quantized on a finite system -this is consistent with the fact that polarization can take continuous value in a finite system.
In practice, since F x , F y only threads a small flux of order O(1/L), one would expect their actual effect to be small, especially at large L. One can then consider the simpler amplitudes Ω|G x T x |Ω and Ω|T y |Ω (|Ω being the ground state in the flux background). These will have magnitudes smaller than one on a finite torus, but as long as it is non-vanishing (in fact we expect it to approach unity in the thermodynamic limit), one can extract the phase of the amplitude, and this phase should give the monopole momentum, which in turn gives the polarization density. More explicitly where ρ x,y are magnitudes that are non-vanishing in the thermodynamic limit (in practice they → 1, see Sec. B 3). Eq. (B4) is in the same spirit with Resta's formula [15] for polarization in one dimension, which is the phase of the (smaller than one) amplitude Ω| exp(ix q x /L)|Ω .
In higher dimensions Resta's amplitude vanishes in the thermodynamic limit and cannot be used to extract polarization [13]. Our prescription using the amplitudes Ω|G x T x |Ω and Ω|T y |Ω can be viewed as a proper generalization of Resta's formula to two dimensions. In fact this prescription has been carried out in previous studies of monopoles in two-dimensional U (1) spin liquids [25]. Strictly speaking our recipe gives the polarization of the ground state in the 2π-flux background |Ω , which is slightly different from the original ground state without the flux |Ω 0 . The two should agree in the thermodynamic limit. To see this let us consider insulators with zero Hall conductance. If there is no symmetry other than charge conservation and translations, the leading order term in the response theory that can cause a magnetic flux to change the polarization is ∆L ∼ α i BE i for some constants α i (i = x, y). This means that a total 2πflux will change polarization by O(B) ∼ O(1/L 2 ). With time-reversal symmetry the leading order term becomes ∼ B 2 E, and the change of polarization in the 2π-flux background becomes O(B 2 ) ∼ O(1/L 4 ). This error will likely be dominated by other finite-size effects such as omitting the flux-threading F x,y in the calculation. This argument is reliable for insulators without Hall conduc-tance since we expect all terms in the response theory to be local and manifestly gauge-invariant.
If the unit cell contains 2 neighboring sites in x direction, i.e. L x even and a unit cell contains (2n, m), (2n + 1, m)(n, m ∈ Z), the above recipe only distributes nonvanishing A field on bonds between unit cells of choice, and A vanishes within each unit cell, which corresponds to calculatingP in Ref. [13]. The actual unit cell structure and geometry do not contribute to monopole translation properties, or polarization, obtained in such ways. In general the polarization and the monopole momentum depend on choice of unit cells.
To obtain the polarization P, with both intra-(classical) and inter-cell effects, we give a recipe to account for the unit cell geometry, applicable to generic systems. To this end we first give a continuum function for gauge field A on the torus which can then be used to determine the discrete gauge fields. Take the distance between neighboring unit cells to be 1 and the Bravais lattice to be square, the continuum gauge field reads Note the function is not single-valued, but is well defined and hence poses no problems for obtaining the gauge fields on the discrete lattice. When put on the lattice, the gauge connection on one bond l is given by l dx · A(x), i.e., the line integral of continuum A along the bond.
Once put on a lattice, the flux close to the "slit" at y = L y , L x − 1 ≤ x ≤ L x should have an O(1) deviation from 2π/(L x L y ) due to the discontinuity in eq (B5) (the total flux threading the unit cell at (L x − 1, L y − 1) is hence 2π/(L x L y ) − 2π). One could compensate for this deviation by altering the gauge connection on bonds inside the slit, such that the deviation is concentrated to a set of elementary plaquettes (i.e., not containing any smaller plaquettes) that contain the point (L x , L y ), whose flux equals 2π/(L x L y )A − 2π (A is the area of the elementary plaquette). This fixes the translation symmetry breaking of flux derived from the continuum recipe. Upon translation T y , one carefully performs a gauge transform on sites in the slit G y to restore the gauge connection as much as possible, the amplitude of Ω|G y T y |Ω is comparable to unity; the phase converges in thermodynamic limit to the conventional polarization P.
The two recipes have the same flux configuration on torus and hence are connected by a gauge transform. However, this gauge transform generally does not commute with G y T y and will change the momentum obtained, consistent with gettingP versus P for the two recipes. For example, in the above square lattice model, we assume a unit cell at (n, m) contains two sites at (n, m), (n + 1/2, m), respectively (in notation of eq (B1) the coordinates read (2n, m), (2n + 1, m)). Then the two recipes built upon eqs (B1) and (B5) differ by a gauge transform on sites with x = L x − 1/2 by the operator e −i y πy Ly ρ(Lx−1/2,y) . From the commutation relation between this gauge transform and the G y T y operation, one can see that the change of momentum from the gauge transform as L → ∞ is precisely δk y = −π ρ(L x − 1/2, y) , which leads to a change in polarization P x −P x = ρ(L x − 1/2, y) /2, in agreement with the intuition that the difference between P andP can be seen as a classical dipole moment within the unit cell.

Review of band theory calculation
For a d-dimensional lattice system with translation symmetries and periodic boundary conditions in all directions, the polarization corresponds intuitively to the dipole moment in each unit cell. For free fermions the polarization contributed by an occupied band is given by the integrated Berry connection (the Wilson loop) in the Brillouin zone [12]: where |u k is the periodic part of the Bloch state at momentum k and the integration is taken over the entire Brillouin zone. (Here u k (r) = e −ik·r ψ(r).) We also dis-cussP [13] if we instead useũ k (r) = e −ik·R ψ(r) where r = R + r i and R is the Bravais lattice vector associated with r. (ũ k (r) =ũ k+Gi (r).) In the presence of gapless Dirac cones, the band theory polarization Eq. (B6) is not uniquely defined. This ambiguity can also be understood from monopole momentum: in a 2π-flux background there are fermion zero modes associated with the Dirac fermions, and filling different zero modes gives different ground states, with different total momenta. We now discuss this within the usual band theory formulation. For concreteness consider a system of spin-1/2 fermions forming two Dirac valleys, say at momenta K, K . The Wilson loop for each spin α in the k 1 direction has a discontinuity of ±π when k 2 passes through K 2 and K 2 . The polarization requires a choice of the jump in P α (π or −π) at each Dirac point. In order for the polarization to be gaugeinvariant, α P α (k 2 ) should be single-valued in the entire Brillouin zone (i.e. no net Chern number). This leads to six different choices of the jumps in P α at the Dirac points. Now from the monopole momentum point of view, in a 2π-flux background there are four zero modes (one from each Dirac cone), and gauge-invariance requires the ground state to fill half of the zero modes, which leads to C 4 2 = 6 different choices -in exact agreement with the band theory consideration.

An example
Our numerical prescription for calculating polarization density through amplitudes like Ω|T y |Ω (|Ω being the many-body ground state in the presence of a uniform 2π flux background) is well-defined for generic many-body systems. In the special case of free fermions we expect our prescription to agree with the band theory results from Eq. (B6). We demonstrate this through an example of a Dirac semimetal (with a specific choice of zero-mode fillings). We consider a square lattice, labeled by two orthogonal unit lattice vectors e 1,2 , with 2 orbitals and where we have relabeled subscripts iα by [l 1 , l 2 ] through l 1 = 2 * i 1 + α, l 2 = i 2 (site i with coordinates (i 1 , i 2 ) in e 1,2 basis) and t is the tuning parameter. Hopping amplitudes on other diagonal or neighboring bonds not covered in eq (B10) vanish. The two limits t = 0, 1 correspond to a square with C 4 rotation and an effective triangular lattice with C 6 rotation, respectively. Diagonalizing this Hamiltonian in momentum space gives gapless dispersion at half-filling. To avoid the ambiguity for Wilson loop operator when crossing Dirac fermions as discussed in sec B 2, we stipulate the two bands for spin up/down has Chern number ±1, respectively, i.e. effectively open an infinitesimal quantum spin hall mass. For the monopole momentum calculation, the gauge connection on links a ij analogous to Landau gauge in eq (B1) gives a total flux of 2π and the quantum spin hall mass indicates that one fills only 2 zero modes of one spin species, giving a monopole carrying spin 1. Fig 1 in the main text shows a comparison of polarization P 1 , calculated numerically using eq (B6) along the direction of reciprocal vector for e 1 and the monopole momentum k 2 along the orthogonal e 2 direction, calculated as in sec B 1 as one tunes t from 0 to 1. The polarization obtained from Eq. (B6) is discretized as summation of Berry phase ln u k |u k+ for 30000 points in Brillouin zone; the momentum is calculated on a lattice of linear size L = 50. For the special cases when t = 0, 1 the momentum k 2 = π, 2π/3 agrees with results in Refs. [25,37].

Appendix C: Polarization and other topological quantities
In 2d the polarization density is not invariant under large gauge transforms if the system has a nonzero Hall conductance σ xy = 0. This is known in band theory, where Eq. (B6) is not invariant under large gauge transforms in real space -in fact this is one way to define integer quantum Hall effect within band theory. Beyond band theory, it is also easy to understand why this is so from the monopole momentum: a 2π-flux induces an extra charge δQ = σ xy in the ground state, which makes the total momentum non-invariant under large gauge transforms. The total polarizations P x L y and P y L x are still well-defined (gauge invariant) mod 1. Similarly, if the system forms a quantum spin Hall insulator, with a nonzero S z spin trapped in a magnetic flux unit, then the polarization is not invariant under a large S z -gauge transform. In all such cases the polarization remains meaningful (unambiguisly defined) for a given gauge if the gauge field remains non-dynamical.
Contrary to the Hall conductance, a nonzero magnetoelectric angle Θ in 3d: does not obstruct the gauge-invariance of polarization density. Within band theory the Θ-angle can be interpreted as the magnetoelectric polarizability [38,39], i.e. a magnetic field induces an extra polarization density The monopole point of view provides a simple understanding of the above relation beyond band theory: when Θ = 0, the monopole traps a fractional charge q = Θ/2π and becomes a dyon [40]. When a magnetic field is turned on, say in z, the monopole also sees the field due to the fractional q. This contributes to the non-commutativity of T x and T y , with the additional phase factor given by qB = (Θ/2π)B. Using Eq. (12) as the definition of polarization we immediately obtain Eq. (C2).

Appendix D: Anomaly from a Fermi surface
In this appendix, we derive the anomaly term Eq. (16) which proves Luttinger theorem in any spatial dimension d. The logic is to partition the Fermi surface into infinitesimally small patches in whose proximity reside "chiral" fermions, that effectively live in (1 + 1)d. The chiral anomaly from each of these fermions adds up to give Eq. (16).
Let us first write down the anomaly for a 1d chiral fermion with the free Hamiltonian ψ † i(±∂ x )ψ, where we set velocity to unity and ± represents right(left)-movers under consideration. Next we couple to theory to both a U (1) electromagnetic field A and an x-translation gauge field (elasticity tetrad) x. The momentum of the chiral fermion k F becomes the coupling constant between the translation gauge field x and the fermion (in analogy to electric charge e as the coupling constant between EM field A and a fermion, here the charge of translation -momentum -mediates the coupling ). Hence the covariant derivative i∂ x,t → i∂ x,t + A x,t + k F x x,t where subscript denotes the space-time component of the 1− form gauge field, omitted hereafter. To obtain the mixed anomaly between A, x, one goes to one higher dimension (2 + 1)d bulk of the chiral fermion -a quantum hall insulator with Chern number C = ±1, with the low-energy topological quantum field theory action Note that since we introduce elasticity tetrads in addition to A, the familiar Chern-Simons term A ∧ dA is modified as such. The mixed term in Eq. (D1) reads ± k F 2π x ∧ dA from which descends a boundary term ∓ k F 2π A ∧ x. Now that we have the desired (1 + 1)d anomaly, consider in d space dimension system, compactify (d − 1) dimensions and derive similar anomaly for the effective (1 + 1)d system along the remaining ith primitive lattice vector direction. We inspect a small patch on the Fermi surface with momentum range (k 1 ±δk 1 /2, · · · k i , · · · k d ± δk d /2) (δk j ≥ 0,the variation of k i is neglected to zeroth order of the anomaly). On such a patch with an "area" ∆s i = j =i δk j , there are ∆si (2π) d−1 j =i L j chiral fermions along the ith direction, each associated with an anomaly ∓ki 2π A ∧ x i (∓ in numerator results from right(left)-movers given by the orientation of the small patch projected onto ith reciprocal vector direction). Adding all patches up, the self Chern-Simons terms vanish and the remaining anomaly reads where F S counts all patches on the Fermi surface, η i = ±1 denotes the orientation of each patch along/against ith reciprocal lattice vector and we use the identity on luttinger volume F S ∆s i η i k i = V F . The second line arises after we introduce translation gauge fields along the other (d−1) directions and the numerator in the first line j =i L j → ∧ j =i x j . The final result puts all x i 's on equal footing and hence it correctly captures the anomaly of Fermi surface under large gauge transforms along any spatial directions. Adding the anomaly term to the Fermi surface theory will make the full theory anomaly-free, as promised in Eq. (17).

Appendix E: LSM anomaly indicators for topological orders
First we review an important notion for a topological order with a global U (1) symmetry in general d dimensions known as the fluxon. Consider an instanton of the A field, which is an operator supported on a (d − 2) dimensional sub-manifolds in space, with dA = 2π on the two complementary spatial dimensions. For d = 2 it is a point flux insertion and for d = 3 it is a unit flux loop. Dirac quantization requires this object to be un-observable from far away. However in a topologically ordered state, there can be nontrivial quasiparticles that carry fractional electric charge, and moving these fractional charges around the instanton will naively produce an observable Aharanov-Bohm phase. The resolution is that the bare instanton is attached with another nontrivial excitation, called the fluxon, from the topological order. The property of the fluxon is such that the combined object becomes unobservable from far away. For example, a fractionally charged quasiparticle will have a nontrivial braiding phase with the fluxon so that it braids trivially with the combination of fluxon and bare instanton. In 2d the fluxon is an anyon excitation and in 3d it is a loop excitation.
In general, anomalies involving a U (1) global symmetry in topological quantum field theories are encoded in the properties of fluxons. Essentially if the fluxon carries a fractional quantum number under other symmetries, in our case lattice translations, then the instanton will also carry the fractional symmetry quantum numbers since it is bound with a fluxon. Since the instanton is supposed to be unobservable, this becomes an anomaly. The fluxon has space-time dimension d − 1, and crystal symmetry fractionalization can be described using a partition function in d space-time dimension: for which the fluxon lives on the boundary of the ddimensional (space-time) manifold, and n A ∈ [0, 1) is the LSM anomaly indicator. At d = 2 the fluxon is an (abelian) anyon particle, and Eq. (E1) means that the fluxon transforms projectively under translation symmetries: This relation has been discussed in Ref. [41]. At d = 3 the fluxon is a loop excitation with finite tension. Eq. (E1) has the following interpretation. First consider a straight fluxon tube, say pointing in z (assuming periodic boundary condition). If we translate the entire loop in the (x, y) plane, T x and T y will commute only up to a phase: where L z is the number of layers of the entire system in z. Another way to describe this property, without relying on having a finite L z , is to consider a closed fluxon loop that links with a dislocation (a line defect in 3D), say with Burgers vector z. Translation symmetries will act on the loop projectively. More generally, for a fluxon loop linked with a dislocation with Burgers vector B, we have Another consequence, following similar reasoning, is a nonabelian three-loop braiding [42] for a fluxon loop and two dislocations.
Appendix F: Square model Numerics to verify boundary Luttinger theorem Our square lattice model consists of spinless fermions with nearest neighbor and diagonal hopping, detailed configuration shown in Fig. 3a. It's modified from πflux square hopping with variation of vertical hopping, time-reversal breaking imaginary diagonal hopping and complex diagonal hopping t further breaks remaining rotation (inversion), reflection symmetries to allow a generic polarization. The only symmetries that remain preserved is translation.
When put on a cylinder geometry with x 2 = 0, L boundaries and periodic along x 1 direction, we can calculate boundary charge densities for e.g. x 2 = 0 boundary as [12] where ρ(x 1 , x 2 ) is the charge density including ions (for a neutral system), x 0 locates deep in the bulk, a is size of unit cell along x 2 direction, A bd denotes any segment covering exactly one unit cell on the boundary and Ω the unit cell area. This amounts to first averaging charge densityρ(x 2 ) = 1 Ω x2+a/2 x2−a/2 dx A bd dx 1 ρ(x 1 , x ) within a window [x 2 − a/2, x 2 + a/2] to smoothen any irrelevant periodic oscillations in the bulk (ρ(x) = 0 for x inside the bulk) while retain the extra charge accumulation [43], then integrating the averaged density. From the field-theoretic point of view, this window functionρ(x) use comes naturally from the application of the longwavelength limit in eq (5) to discrete lattices. In continuum, one identifies each unit cell with a single point x and hence the vector potential A 0 (x) couples to the average density inside the unit cellρ(x). On the other hand, the lattice-scale oscillation of bare ρ(x) renders it incompatible with continuum limit in long wavelength.  (14) while (c) has non-vanishing Chern number and gauge-dependent polarization, boundary k F , so boundary Luttinger theorem applies with an appropriate Berry connection integral rule.
ρ(x) for the boundary charge density, however,matches withP(see last paragraph).
Similarly, we get boundary k F as one varies chemical potential. The bulk polarization is calculated by Eq. (B6). For simplicity we put the positive ions at sites with integral coordinates in units of Bravais lattice vectors, i.e. site (0, 0) and its equivalents by lattice translations. The ions don't contribute to polarization in this way; Eq (B6) gives the entire polarization then.
As we vary the parameters t, , the system enters multiple physical regimes with a gapped bulk. For example, the system hosts non-chiral edge states in Fig. 3b and when chemical potential stays inside the bulk gap, the edge density ρ, Fermi momentum k F and bulk polarization P always obeys Eq. (14). (We note it's important for the bulk to remain insulating with the chemical potential in between the gap.) A relatively trivial scenario in Fig. 3d is in the absence of edge states, the edge density equals bulk polarization in line with previous knowledge of polarization. A tricky case is when the model has a nonzero Chern number C shown in Fig. 3c and edge density ρ, k F , polarization will change upon a large gauge transform along the orthogonal direction, i.e. gaugedependent. We find that under a fixed gauge, Eq. (14) still holds given appropriate recipe for bulk polarization Eq. (B6), i.e., 2πC discontinuity of Wilson loop phase θ k1 = dk 2 u k |i∂ k2 |u k occurs only at k 1 = 0. In all cases, the momenta of monopole on a torus geometry satisfy (k 1 , k 2 ) = 2π(−P 2 , P 1 ) obtained by the method in appendix B 1.
Finally we remark that all calculations above apply also to the inter-cell part of polarizationP , when we calculate the boundary charge density as ρ 0 = na −∞ dx 2 A bd dx 1 ρ(x 1 , x 2 ) (n ∈ Z) instead of eq (F1).
The relation between these two reads [17] ρ bd = ρ 0 − 1 a na+a/2 na−a/2 where we take x 0 = na in eq (F1) and use the neutrality condition. In passing we remark this rewriting makes explicit the equivalence between the window function method and the charge density as derived in Appendix A. There the charge density in 2 spatial dimensions reads where Ω x0 is the unit cell at x 0 and ∇ Ω denotes gradient of ρ w.r.t its value at the same sublattice in neighboring unit cells. When there's a boundary to vacuum, the boundary charge x0 ρ(x 0 ) has two parts: the first term integrates to ρ 0 and second term integrates to give bulk electric dipole moment. Hence it agrees with eq (F2).
It is now clear that ρ 0 extracts the excess charge at the boundary. The boundaries we considered preserve complete unit cells in bulk. This extra charge accumulation part depends solely on the inter unit cell structure. Hence we use the gauge recipe eq (B1) to calculate monopole momentum related toP , whose A fields reside only on bonds between different unit cells and indeed they agree with Berry connection integral eq (B6) using periodic functionũ k (r).