Many-Body Electric Multipole Operators in Extended Systems

TL;DR

This paper develops a formalism for defining and calculating higher-order multipole moments in extended quantum systems with periodic boundary conditions, enabling better understanding of topological phases and electric multipole responses.

Contribution

It introduces many-body operators for higher multipole moments and demonstrates their effectiveness in characterizing topological quadrupole and octupole phases.

Findings

Operators correctly identify phase diagrams of topological multipole models

Capture adiabatic multipole pumping phenomena

Distinguish bulk multipole moments from corner charge effects

Abstract

The quantum mechanical position operators, and their products, are not well-defined in systems obeying periodic boundary conditions. Here we extend the work of Resta who developed a formalism to calculate the electronic polarization as an expectation value of a many-body operator, to include higher multipole moments, e.g., quadrupole and octupole. We define $n$-th order multipole operators whose expectation values can be used to calculate the $n$-th multipole moment when all of the lower moments are vanishing (modulo a quantum). We show that changes in our operators are tied to flows of $n-1$-st multipole currents, and encode the adiabatic evolution of the system in the presence of an $n-1$-st gradient of the electric field. Finally, we test our operators on a set of tightbinding models to show that they correctly determine the phase diagrams of topological quadrupole and octupole models, capture an adiabatic quadrupole pump, and distinguish a bulk quadrupole moment from other mechanisms that generate corner charges.

Equations

$${\hat{\mathcal{U}}}_j=\exp \left[\frac{2\pi i{\hat{X}}^j}{L_j}\right]$$

$${\hat{\mathcal{U}}}_{ab}^Q$$

$${\hat{\mathcal{U}}}_{abc}^O$$

$$Q^{ab}$$

$$O^{abc}$$

$$H=\sum _{n=1}^{N_{elec}}\sum _{i=x,y,z}\frac{({\hat{p}}_n^i{)}^2}{2m}+V(\hat{\mathbf{X}}).$$

$${\hat{U}}_{\lambda }(\alpha ){\hat{p}}_n^i{\hat{U}}_{\lambda }^{†}(\alpha )={\hat{p}}_n^i-\hslash \alpha \frac{\partial \lambda (\hat{\mathbf{X}})}{\partial {\hat{x}}_n^i}.$$

$$U_{{\hat{X}}^1}(2\pi /L_1)|{\Psi }_0⟩\approx$$

$$e^{i{\gamma }_{p1}}\left(|{\Psi }_0⟩-\frac{2\pi \hslash }{m{L}_1}\sum _{j\ne 0}|{\Psi }_j⟩\frac{⟨{\Psi }_j|{\hat{\mathcal{P}}}^1|{\Psi }_0⟩}{E_0-E_j}\right)$$

$$\frac{d{P}^1}{dt}\approx \frac{ie\hslash }{m{L}_1}\sum _{j\ne 0}⟨{\dot{\Psi }}_0|{\Psi }_j⟩\frac{⟨{\Psi }_j|\widehat{{\mathcal{P}}^1}|{\Psi }_0⟩}{E_0-E_j}+\text{c.c},$$

$$U_{{\hat{\mathcal{Q}}}^{12}}(2\pi /L_1{L}_2)|{\Psi }_0⟩\approx$$

$$e^{i{\gamma }_{q12}}\left(|{\Psi }_0⟩-\frac{2\pi \hslash }{L_1{L}_2}\sum _{j\ne 0}|{\Psi }_j⟩\frac{⟨{\Psi }_j|{\hat{\mathcal{J}}}_D^{12}|{\Psi }_0⟩}{E_0-E_j}\right)$$

$$\frac{d{Q}^{12}}{dt}\approx \frac{ie\hslash }{L_1{L}_2}\sum _{j\ne 0}⟨{\dot{\Psi }}_0|{\Psi }_j⟩\frac{⟨{\Psi }_j|{\hat{\mathcal{J}}}_D^{12}|{\Psi }_0⟩}{E_0-E_j}+\text{c.c}$$

$${\gamma }_{p1}$$

$$A^i=-\frac{ht{\sigma }^{ij}{x}^j}{e{L}_1{L}_{2}T}i,j=1,2,$$

$${\gamma }_{q12}$$

$$⟨{\Phi }_0|{\hat{\mathcal{U}}}_a^P|{\Phi }_0⟩=⟨{\Phi }_0|{\hat{T}}_{{\mathbf{a}}_c}^{-1}{\hat{T}}_{{\mathbf{a}}_c}{\hat{\mathcal{U}}}_a^P{\hat{T}}_{{\mathbf{a}}_c}^{-1}{\hat{T}}_{{\mathbf{a}}_c}|{\Phi }_0⟩$$

$$⟨{\Phi }_0|{\hat{\mathcal{U}}}_{ab}^Q|{\Phi }_0⟩=\exp \left[2\pi i{n}_{0ab}{\delta }_{ac}{\delta }_{bc}\right]$$

$$H(\mathbf{k},\Phi )=$$

$$\displaystyle\left(\begin{array}[]{cccc}\delta&\Gamma_{x}(k_{x})&\Gamma_{y}(k_{y})&0\\ \Gamma_{x}(-k_{x})&-\delta&0&\Gamma_{y}(k_{y})\\ \Gamma_{y}(-k_{y})&0&-\delta&e^{i\Phi}\Gamma_{x}(k_{x})\\ 0&\Gamma_{y}(-k_{y})&e^{-i\Phi}\Gamma_{x}(-k_{x})&\delta\end{array}\right)$$

$$H(\mathbf{k})$$

$$⟨{\hat{\mathcal{U}}}_{xy}^{\prime }⟩=⟨{\hat{\mathcal{U}}}_{xy}⟩⟨{\hat{\mathcal{U}}}_y⟩+\mathcal{O}\left(\frac{1}{E_{\mathrm{gap}}}\right),$$

$${\hat{\mathcal{U}}}_{xy}^Q=\exp \left(\frac{2\pi i{\sum }_j{x}_j{y}_j}{L_x{L}_y}\right).$$

$$T_{\mathbf{R}}{\hat{\mathcal{U}}}_{xy}^Q{T}_{\mathbf{R}}^{-1}={\hat{\mathcal{U}}}_{xy}^Q\exp \left[2\pi i\sum _j{h}_1{h}_2/N_x{N}_y\right]\exp \left[2\pi i{h}_1\sum _j{y}_j/N_x{L}_y\right]\exp \left[2\pi i{h}_2\sum _j{x}_j/L_x{N}_y\right]$$

$$={\hat{\mathcal{U}}}_{xy}^Q\exp \left[2\pi i\sum _{\mathbf{R}}\sum _{\alpha =1}^{\nu }{h}_1{h}_2/N_x{N}_y\right]\exp \left[2\pi i{h}_1\sum _{\mathbf{R}}\sum _{\alpha =1}^{\nu }{y}_{\mathbf{R},\alpha }/N_x{L}_y\right]\exp \left[2\pi i{h}_2\sum _{\mathbf{R}}\sum _{\alpha =1}^{\nu }{x}_{\mathbf{R},\alpha }/L_x{N}_y\right]$$

$$={\hat{\mathcal{U}}}_{xy}^Q\exp \left[2\pi i\nu h_1{h}_2\right]\exp \left[2\pi i{h}_1\sum _{n_2}\sum _{\alpha =1}^{\nu }({\bar{y}}_{\alpha }+n_2{a}_y)/L_y\right]\exp \left[2\pi i{h}_2\sum _{n_1}\sum _{\alpha =1}^{\nu }({\bar{x}}_{\alpha }+n_1{a}_x)/L_x\right]$$

$$={\hat{\mathcal{U}}}_{xy}^Q\exp \left[2\pi i\nu h_1{h}_2\right]\exp \left[\frac{2\pi i{h}_1}{a_y}\sum _{\alpha =1}^{\nu }{\bar{y}}_{\alpha }\right]\exp \left[\frac{2\pi i{h}_2}{a_x}\sum _{\alpha =1}^{\nu }{\bar{x}}_{\alpha }\right]\exp \left[\pi i{h}_1\nu (N_y+1)\right]\exp \left[\pi i{h}_2\nu (N_x+1)\right],$$

$$\exp \left[2\pi i\nu \right]\exp \left[\frac{2\pi i}{a_y}\sum _{\alpha =1}^{\nu }{\bar{y}}_{\alpha }\right]\exp \left[\frac{2\pi i}{a_x}\sum _{\alpha =1}^{\nu }{\bar{x}}_{\alpha }\right]\exp \left[\pi i\nu (N_y+N_y+2)\right]=1.$$

$$\exp \left[\frac{2\pi i\nu (N_x+1)(N_y+1)(N_z+1)}{4}\right]=1.$$

$$P_x=\frac{e{L}_x}{2\pi \mathcal{V}}\mathrm{Im}\log \left[⟨{\hat{\mathcal{U}}}_x⟩\right]$$

$$P_x=e{x}_0+e\frac{N_x+1}{2},$$

$${\psi }_{x_0}(x)=\frac{1}{\sqrt{{\pi }^{1/2}{\sigma }_x}}\exp \left[-\frac{(x-x_0{)}^2}{2{\sigma }_x^2}\right].$$

$$⟨{\hat{\mathcal{U}}}_x⟩=\exp \left[\frac{2\pi i{x}_0}{N_x}\right]\exp \left[-\frac{{\pi }^2{\sigma }_x^2}{N_x^2}\right].$$

$$⟨{\hat{\mathcal{U}}}_x⟩=\exp \left[\frac{2\pi i{x}_0}{N_x}\right]\exp \left[-\frac{{\pi }^2{\sigma }_x^2}{N_x^2}\right]$$

Full text

The modern theory of electric polarization in crystals has had a wide-ranging impact over the past 25 yearsKing-Smith and Vanderbilt (1993); Vanderbilt and King-Smith (1993); Resta (1994); Ortiz and Martin (1994); Resta and Vanderbilt (2007); Vanderbilt (2018). Aspects of the theory are not only useful for predictions of observable properties of real materials, but are tightly connected to a variety of topological insulator phenomena. For example, insulating crystals having a symmetry that would naively forbid the existence of a dipole moment, may instead have non-vanishing, but quantized, polarization if they are in a topological insulator phaseZak (1989); Qi et al. (2008); Hughes et al. (2011); Turner et al. (2012). Furthermore, by tracking changes in polarization in an adiabatic cycle, one can observe quantized charge pumping characterized by a Chern number topological invariantThouless (1983). More recently, the theory of quantized electronic polarization in topological insulators was extended to account for higher electric multipole momentsBenalcazar et al. (2017a, b). Therein it was shown that two and three dimensional, non-interacting topological crystalline phases of matter can support quantized electric quadrupole or octupole moments.

The electric polarization is calculated by the Berry phase of the electronic energy bands over non-contractible cycles in the Brillouin zoneKing-Smith and Vanderbilt (1993); Vanderbilt and King-Smith (1993). Similarly, Ref. Benalcazar et al., 2017a adapts this Berry phase formalism to describe the hierarchy of electric multipole moments. While the Berry phase calculations for higher multipole moments can also be used on non-topological insulators with moments outside the quantized regime, the formulae are not directly applicable to many-body interacting systems. For polarization, one can consider the Berry phase of a many-body ground state parameterized by twisted periodic boundary conditionsOrtiz and Martin (1994); Souza et al. (2000), but it is not obvious how this approach can be practically extended to the nested Berry phase (Wilson loop) approach to calculating higher multipole moments.

In this letter we instead provide a formulation of the quadrupole and octupole moments in terms of ground-state expectation values of many-body operators, analogous to Resta’s formulation of the charge polarizationResta (1998). We define many-body multipole operators, and discuss the connection between these operators and the flow of adiabatic multipole currents. In addition, we test our operators for several explicit model systems. We show that our formulae precisely capture the phase diagrams of the quadrupole and octupole models from Ref. Benalcazar et al., 2017a. Furthermore, we use our formula to track the changes in quadrupole moment during the higher order pumping process in Ref. Benalcazar et al., 2017b, and we show that the results match the expected phenomena exactly. As an important comparison, we go on to test our formula on a model with corner charges but no bulk quadrupole moment. We show that our operator does not detect any difference between the phases of this model with and without corner charges, as we expect since there is no change in the bulk quadrupole moment, despite the presence of corner charges.

I Motivation

Let us begin by recounting the calculation of the electric polarization using a many-body expectation valueResta (1998). Consider a translationally invariant, crystalline insulator with periodic boundary conditions. It is well known that the many-body position operator $\hat{{\bf{X}}}=\sum_{n}\hat{{\bf{x}}}_{n},\;(n=1,\ldots N_{e})$ cannot be naively used in the calculation of the electric polarization in extended (and/or periodic) systems. In such systems the operator is not well-defined because it can transform a state in the Hilbert space to a state outside the Hilbert space, e.g., it can take a normalizable state to a non-normalizable state, or one obeying periodic boundary conditions to one that violates the boundary conditions. Instead, the operator

[TABLE]

can be employed in calculations of the polarization of a many-body ground state $|\Phi_{0}\rangle$, i.e., $P^{j}=\frac{eL_{j}}{2\pi\mathcal{V}}{\rm{Im}}\log\langle\Phi_{0}|\hat{\mathcal{U}}_{j}|\Phi_{0}\rangle\equiv\frac{eL_{j}}{2\pi\mathcal{V}}{\rm{Im}}\log z^{(P)}_{j},$ where $\mathcal{V}$ is the volume of the system. This formula for $P^{j}$ approximates a derivative with respect to the many-body momentum by a finite difference Souza et al. (2000), so it is only strictly true in the thermodynamic limit $L_{j}\rightarrow\infty$. For a non-degenerate, insulating ground state one finds $\hat{\mathcal{U}}_{j}|\Phi_{0}\rangle=e^{i\gamma_{j}}|\Phi_{0}\rangle+O(1/L_{j})$ as the thermodynamic limit is approachedResta (1998), i.e., $|\Phi_{0}\rangle$ becomes an eigenstate of $\hat{\mathcal{U}}_{j}$ in the thermodynamic limit111In taking the thermodynamic limit one must take $L_{j}\to\infty$ first before the transverse directions. (we recount this result briefly below) with a polarization given by $P^{j}=\frac{e\gamma_{j}}{2\pi}\frac{L_{j}}{\mathcal{V}}.$ The polarization has an ambiguity from the choice of the branch of the $\log,$ often referred to as the quantum of polarization, i.e., $\gamma_{j}\equiv\gamma_{j}+2\pi n$.

Following this line of reasoning we can consider higher electric multipole moments such as the quadrupole $\hat{q}^{ij}=\sum_{n}\hat{x}^{i}_{n}\hat{x}^{j}_{n}$ or octopole $\hat{o}^{ijk}=\sum_{n}\hat{x}^{i}_{n}\hat{x}^{j}_{n}\hat{x}^{k}_{n}.$ These operators are problematic in periodic systems for the same reasons as the many-body position operator, but we can also consider exponentiated versions. Let ${\bf{b}}_{a}$ be the set of reciprocal lattice vectors satisfying ${\bf{b}}_{a}\cdot{\bf{a}}_{b}=\delta_{ab}$ for the primitive unit vectors ${\bf{a}}_{b}.$ Then we can consider the operators

[TABLE]

where $N_{a}$ is the number of unit cells in the $a$-th lattice direction. The bulk moments are then defined as

[TABLE]

where we have left the crucial step of taking the thermodynamic limit implicit in these formulas. For the majority of this article we will direct our focus to the quadrupole operator, and the octupole case can be studied essentially mutatis mutandis.

II Moments as integrated current

To confirm that our definition of the quadrupole moment matches the expected physical observables we will generalize the arguments from Ref. Resta, 1998. Consider a generic many-body Hamiltonian

[TABLE]

Let us focus on the unitary transformation $\hat{U}_{\lambda}(\alpha)=\exp\left(i\alpha\lambda({\bf{\hat{X}}})\right)$, for an arbitrary real differentiable function $\lambda(\bf{\hat{X}})$, that acts on each momentum operator as

[TABLE]

Now we define $H(\alpha)=\hat{U}_{\lambda}(\alpha)H\hat{U}_{\lambda}^{\dagger}(\alpha),$ noting that the difference between $H(0)$ and $H(\alpha)$ is the presence of an electromagnetic vector potential acting on each electron, e.g., $A^{i}({\bf{\hat{x}}}_{n})=-\frac{\hbar\alpha}{e}\frac{\partial\lambda({\bf{\hat{X}}})}{\partial\hat{x}^{i}_{n}}$ for the $n$-th electron. Assuming that the ground state of $H(0)$, $|\Psi_{0}\rangle$, is always non-degenerate, then the relation $H(\alpha)U_{\lambda}(\alpha)|\Psi_{0}\rangle=E_{0}U_{\lambda}(\alpha)|\Psi_{0}\rangle$ holds.

Let us take the particular case where $\lambda({\bf{\hat{X}}})=\hat{X}^{1},$ which will allow us to reproduce Resta’s argument. This generates the constant vector potential $A^{1}=-\frac{\hbar\alpha}{e},A^{2}=A^{3}=0$. To have $U_{\hat{X}^{1}}(\alpha)|\Psi_{0}\rangle$ satisfy periodic boundary conditions, we can choose $\alpha=2\pi/L_{1}.$ When approaching the thermodynamic limit, $L_{1}$ is large and we can treat $\alpha$ as a small parameter. Thus we can treat $U_{\hat{X}^{1}}(2\pi/L_{1})|\Psi_{0}\rangle,$ which is an eigenstate of $H(2\pi/L_{1}),$ as a perturbed version of the initial ground state $|\Psi_{0}\rangle,$ and hence we can expand $U_{\hat{X}^{1}}(2\pi/L_{1})|\Psi_{0}\rangle$ in terms of the eigenstates of $H(0).$ The perturbation term in the Hamiltonian is $H^{\prime}=-\frac{2\pi\hbar}{mL_{x}}\sum_{n=1}^{N_{elec}}\hat{p}_{n}^{1},$ and the leading order correction to the unperturbed ground state yields

[TABLE]

where $|\Psi_{j}\rangle$ are the excited states of $H(0)$ with energies $E_{j},$ and $\hat{\mathcal{P}}^{i}$ is the many-body momentum operator of all the electrons. In addition to the usual form, we have allowed for a phase factor $\gamma_{p1}$ which determines the $1$-component of the electric polarizationResta (1998). As we stated earlier, we now see explicitly that in the thermodynamic limit $U_{\hat{X}^{1}}(2\pi/L_{1})|\Psi_{0}\rangle=e^{i\gamma_{p1}}|\Psi_{0}\rangle.$ To connect to the physical polarization we can show that the time-derivative of this polarization definition matches the physically expected result of an electric current. Indeed Ref. Resta, 1998 shows from the perturbation theory calculation above that in 1D

[TABLE]

which is the adiabatic electric charge currentThouless (1983).

Now consider the quadrupole moment $\hat{Q}^{12}$. We will consider a unitary transformation $\hat{U}_{\lambda}(\alpha)$ with $\lambda({\bf{\hat{X}}})=\sum_{n=1}^{N_{elec}}\hat{x}^{1}_{n}\hat{x}^{2}_{n}$ and with $\alpha=2\pi/L_{1}L_{2}$. Using the same arguments as the case for polarization, we can use perturbation theory to calculate

[TABLE]

where $\hat{\mathcal{J}}_{D}^{12}=\tfrac{1}{m}\sum_{n=1}^{N_{elec}}(\hat{p}_{n}^{1}\hat{x}_{n}^{2}+\hat{p}_{n}^{2}\hat{x}_{n}^{1})$ is the dipole-current operator. Note that the velocity $d\hat{x}^{i}_{n}/dt=(i/\hbar)[H,\hat{x}^{i}_{n}]=p^{i}_{n}/m$ is like the time derivative of the dipole operator, while the time derivative of the quadrupole operator is $d(\hat{x}^{1}_{n}\hat{x}^{2}_{n})/dt=(i/\hbar)[H,\hat{x}^{1}_{n}\hat{x}^{2}_{n}]=(1/m)(\hat{p}^{1}_{n}\hat{x}^{2}_{n}+\hat{p}^{2}_{n}\hat{x}^{1}_{n}),$ hence our notion of a dipole current operator. From this result we can further calculate

[TABLE]

which relates the time derivative of our definition of the quadrupole to a dipole-current as we would physically expect. We note that the above expression for $\hat{\mathcal{J}}_{D}^{12}$ (and thus Eqs. II and II) cannot be evaluated in a periodic system in its current form because it contains position operators. Nonetheless, the notion of a dipole current is useful for understanding the origin of a bulk quadrupole moment as the result of a process that rearranges charge moments, even in extended systems. The problem with periodicity is reminiscent of the polarization operator ${\bf\hat{X}}$, which cannot be evaluated directly in a periodic system, but can be computed from the operator in Eq. 1. Finding an alternate expression for $\hat{\mathcal{J}}_{D}^{12}$ that is compatible with periodic boundary conditions is an interesting direction for future work.

III Adiabatic evolution

For an intuitive understanding we can characterize the action of our operators in the language of adiabatic evolution. The action of $\hat{U}_{\hat{X}^{1}}(2\pi/L_{1})$ on a state can be treated as adiabatic evolution from a system with a vanishing vector potential to one with $A^{1}=-h/eL_{1}$. The process can be accomplished through a time-dependent vector potential of the form $A^{1}=-\frac{ht}{eL_{1}T},$ and which is [math] for $t<0,$ and $-h/eL_{1}$ for $t>T.$ The process inserts one magnetic flux quantum into the periodic cycle in the $1$-direction. From the Faraday effect, this is equivalent to turning on a uniform electric field in the $1$-direction during the time interval $t\in[0,T].$ We are interested in the evolution of the non-degenerate ground-state of a gapped, neutral insulator, and so the presence of the uniform electric field will only activate the additional phase factor $\gamma_{p1}=-\hbar^{-1}\int_{0}^{T}{\bf{d}}\cdot{\bf{E}}\,dt$ in the adiabatic limit, where ${\bf{d}}$ is the total dipole moment of the material. Besides the weak, time-varying vector potential, we assume the insulator is otherwise static so that ${\bf{d}}$ is time-independent. Hence, we find the phase (c.f. Eq. II)

[TABLE]

which exactly produces the right value of $P^{1}=\frac{e}{2\pi}{\rm{Im}}\log\langle\Phi_{0}|\hat{U}_{\hat{X}^{1}}(2\pi/L_{1})|\Phi_{0}\rangle$ in the thermodynamic limit.

The essential idea in the previous paragraph is that, to detect a polarization in the material, we adiabatically turn on a small electric field and track how the phase of the ground state responds during this process. For the quadrupole, instead of turning on a uniform electric field, we turn on a uniform electric field gradient. The quadrupolar operator $U_{\hat{\mathcal{Q}}^{12}}(2\pi/L_{1}L_{2})$ can be interpreted as adiabatically evolving our ground state from a vanishing vector potential at $t=0$ to the vector potential $A^{i}=-\frac{h\sigma^{ij}x^{j}}{eL_{1}L_{2}}$ at $t=T.$ This process can be carried out using a time-dependent vector potential

[TABLE]

where $\sigma^{12}=\sigma^{21}=1,\sigma^{11}=\sigma^{22}=0,$ and $T$ is large. This vector potential represents a constant electric field gradient. Thus, a natural way to interpret the proposed quadrupole operator is as an evolution process where a small, uniform electric field gradient is turned on for a finite amount of time. During this time, the ground state will develop the phase factor $\gamma_{q12}$ shown in Eq. II, in the thermodynamic limit. Since we are applying the electric field gradient in a neutral, unpolarized insulator, and the $Q^{12}$ quadrupole moment couples to an electric field gradient $\tfrac{1}{2}(\partial_{1}E^{2}+\partial_{2}E^{1})$, its contribution to the phase will be $\gamma_{q12}=-\tfrac{1}{2\hbar}\int_{0}^{T}(q^{12}\partial_{1}E^{2}+q^{21}\partial_{2}E^{1}).$ We assume that the quadrupole moment of our system is static so we can simplify this expression to

[TABLE]

which confirms our operator definition of the quadrupole moment.

IV Practical evaluation

After having shown that $\hat{\mathcal{U}}^{Q}_{ab}$ can be used to determine dipole currents in the bulk, and that it has a physical interpretation via adiabatic evolution in an electric field gradient, we now turn to the practical evaluation of the expectation value of this operator. To determine if $\hat{\mathcal{U}}^{Q}_{ab}$ can have a well-defined, non-vanishing expectation value $z^{(Q)}_{ab}$ in periodic systems we need to evaluate its transformation properties under translations. Let us recall the argument from Ref. Aligia and Ortiz, 1999 that shows that $z_{j}^{(P)}$ is non-vanishing only when the expectation value of $\hat{\mathcal{U}}_{j}$ transforms in the trivial representation of the translation group. Assuming translation invariance, the ground state $|\Phi_{0}\rangle$ is an eigenstate of the translation operator carrying many-body momentum ${\bf{K}}.$ We find

[TABLE]

where we translated all the electrons in the $c$-th primitive direction using $\hat{T}_{{\bf{a}}_{c}},$ and assumed $|\Phi_{0}\rangle$ is an eigenstate of the total number operator. We see that if $n_{0a}=N_{e}/N_{a}$ is an integer for $N_{a}$ unit cells in the $a$-th direction, then the expectation value does not have to vanish. We can relax $n_{0a}$ to be rational with a suitable modification of $\hat{\mathcal{U}}_{a}$Aligia and Ortiz (1999).

Repeating for the quadrupole we find

[TABLE]

where $n_{0ab}=N_{e}/N_{a}N_{b}$, is assumed to be an integer, and we have applied the results of Ref. Resta, 1998 to evaluate $e^{2\pi i\hat{X}^{j}(b_{a}^{j}\delta_{bc}+\delta_{ac}b_{b}^{j})/N_{a}N_{b}}|\Phi_{0}\rangle=e^{i\Gamma}|\Phi_{0}\rangle+O(1/L).$ Hence, in the thermodynamic limit we find the relation $\langle\Phi_{0}|\hat{\mathcal{U}}^{Q}_{ab}|\Phi_{0}\rangle=e^{i\Gamma}\langle\Phi_{0}|\hat{\mathcal{U}}^{Q}_{ab}|\Phi_{0}\rangle,$ and $z^{(Q)}_{ab}$ must therefore vanish unless $\Gamma=2\pi p$ for some integer $p.$ The phase factor $\Gamma$ is different than the phase generated by translations of the polarization operator $\hat{\mathcal{U}}^{P}_{a}$, which depends on the particle number $n_{0a}$. Instead, $\Gamma$ depends on the polarization $\hat{X}^{j}$. Specifying that the planar filling factor $n_{0ab}$ is an integer (or rational fraction for suitably generalized $\hat{\mathcal{U}}^{Q}_{ab}$Aligia and Ortiz (1999)) will not be enough in this case, and one must also specify that the polarization vanishes (up to a quantum).

However, even after satisfying this constraint, the expectation value of $\hat{\mathcal{U}}^{Q}_{ab}$ is not invariant under translations in general. Constraining the polarization to vanish (which should ensure $\Gamma=2\pi p$) will not fix the issue in finite-sized systems due to the non-vanishing fluctuations of the dipole moment. The dipole fluctuations are a measure of the Wannier function localization length, or alternatively, the non-vanishing correlation length controlled by the insulating gap. When the dipole fluctuations are non-vanishing, the state $|\Phi_{0}\rangle$ is not an eigenstate of the $\hat{X}^{j}$ operator, so even with no average dipole, evaluating the expectation value of the quadrupole operator can be problematic. If we have single particle orbitals in a zero-correlation length limit where Wannier functions are point-localized, then the dipole fluctuations vanish and the many-body ground state will be an exact eigenstate of the $\hat{\mathcal{U}}^{P}_{a}$ operators. For this situation the expectation value of $\hat{\mathcal{U}}^{Q}_{ab}$ is well-defined even in finite-size as long as the polarization vanishes, and the particle filling is integer. Indeed, for many-body ground states that are an exact eigenstate of $\hat{\mathcal{U}}^{P}_{a}$ the magnitude of the expectation value of $\hat{\mathcal{U}}^{Q}_{ab}$ will tend to unity in the thermodynamic limit. To evaluate the expectation value of the quadrupole operator more generally we must be more careful.222Similar problems would occur for the polarization operator if the ground state was not an eigenstate of the number operator, but instead had fluctuations in particle number. The effect of dipole fluctuations on $\hat{\mathcal{U}}^{Q}_{ab}$ is discussed further in Appendix C.

In this work we replace the quadrupole operator by an approximation: we evaluate the expectation value of $\hat{\mathcal{U}}_{ab}^{Q}$ on a finite supercell of size $L_{x}\times L_{y}$, then extend the truncated operator periodically outside the supercell. We then evaluate the expectation value of this periodic operator. The error introduced by this approximation depends on the size of the supercell we use, but improves as the ratio of the correlation/localization length and the characteristic length of the supercell goes to zero. From our analytic evaluation of the expectation value of $\hat{\mathcal{U}}_{ab}^{Q}$ for Gaussian charge configurations in Appendix C, we expect that the magnitude of the expectation value will approach a constant in a thermodynamic limit where both $N_{x}$ and $N_{y}$ are taken to infinity together. This constant will not generically be unity, and it will depend on the size of the dipole fluctuations in the $x$ and $y$ directions as well as the aspect ratio $N_{x}/N_{y}$. In this limit the phase factor that determines the quadrupole moment converges to the correct value as well. Discussions of finite size error and more details about the effects of dipole fluctuations are shown in Appendices B and C. We will leave a systematic study of this approximation and alternatives to future work.

V Example Calculations

V.1 Multipole Models with Mirror or $C_{4}$ Symmetry

Now we provide some examples in which we calculate the multipole moments using our operators in tightbinding models. Our first example will be focused on the insulating phases of a particular tightbinding modelBenalcazar et al. (2017a, b) on a square lattice with four degrees of freedom per unit cell, which, for simplicity, we will treat as spinless orbitals (we also provide calculations for a related octupole modelBenalcazar et al. (2017a), and in the subsequent subsection we provide calculations for another tightbinding model with a non-vanishing quadrupole moment, but with $C_{4}T$ symmetrySchindler et al. (2018)). The unit cell basis and tunneling terms are illustrated in Fig. 1, and the Bloch Hamiltonian for a system with periodic boundary conditions is given by

[TABLE]

where $\Gamma_{x}(k_{x})=\gamma_{x}+\lambda_{x}e^{ik_{x}},\Gamma_{y}(k_{y})=\gamma_{y}+\lambda_{y}e^{ik_{y}},$ $\delta,\gamma_{x/y},\lambda_{x/y}$ are real parameters representing onsite energies, intra-cell tunneling, and inter-cell tunneling respectively, and $\Phi$ is the flux in each (intra- and inter-cell) plaquette in our chosen gauge.

This model can be tuned to a variety of insulator phases, and we will consider several test cases. For case (i) we explore this model in the quadrupolar phase protected by $M_{x}$ and $M_{y}$ mirror symmetries where $\Phi=\pi,$ $\delta=0,$ and we fix $\lambda_{x}=\lambda_{y}\equiv\lambda$ to plot a phase diagram as a function of $\gamma_{x}/\lambda$ and $\gamma_{y}/\lambda$ in the interval $[-2,2].$ We find that the calculation of the quadrupole moment using our operator (see Fig. 2a), matches the expected phase diagram from Ref. Benalcazar et al., 2017a where the system has a quantized quadrupole moment with magnitude $Q^{xy}=e/2$ when $\gamma_{x}/\lambda,\gamma_{y}/\lambda$ are both in the interval $[-1,1].$

Let us make some technical remarks about this calculation. Although difficult to see by eye, there is some deviation from the expected phase diagram very close to the phase boundaries where the correlation length is increasing without bound, but this deviation is just a finite-size effect as it decreases rapidly with the size of the supercell. Additionally, while the phase of the operator $\hat{\mathcal{U}}_{ab}^{Q}$ yields the quadrupole moment of the ground state in the thermodynamic limit, it can only be extracted if the magnitude of the expectation value is non-vanishing. According to Eq. II, one would expect its magnitude to approach one in the thermodynamic limit; however, we find that, for our choice of boundary conditions and our choice for the approach to the thermodynamic limit, the magnitude does not go to unity in the presence of dipole fluctuations. In Appendix B, we show the scaling of the magnitude of the operator as a function of system size for two points in the phase diagram shown in Fig. 2a. Furthermore, we discuss the complications associated with dipole fluctuations in more detail in Appendix C.

For case (ii) we consider a variation of the quadrupole phase protected by $C_{4}$ symmetry where we tune $\gamma_{x}=\gamma_{y}\equiv\gamma$, $\lambda_{x}=\lambda_{y}\equiv\lambda$ and allow for $\Phi\in[0,\pi].$ At $\Phi=0$ the model is gapless at half filling, but for any finite $\Phi$ in this interval we expect to find a quantized quadrupole of $Q^{xy}=e/2$ when $|\lambda|>|\gamma|.$ We show this calculation in the extended planes along the diagonals of the phase diagram in Fig. 2a, and find that our operator correctly reproduces the phase diagram.

For case (iii) we consider the adiabatic dipole pumping process introduced in Refs. Benalcazar et al., 2017a, b where the quadrupole moment is continuously tuned from the topological phase where $Q^{xy}=e/2$ to the trivial phase where $Q^{xy}=0.$ We can parameterize this pumping process by $\lambda_{x}=\lambda_{y}=1,$ $\delta\to-\sin(\theta(t)),$ and $\gamma=1/2(1-\cos(\theta(t))).$ We show the result in Fig. 2b, in which we plot the result of our operator compared to calculations of the corner charge (for an open, square geometry) and the edge polarization on a cylinder, and we find they all match exactly as expected for a continuously varying bulk quadrupole moment $Q^{xy}.$

For case (iv) we tune the system to $\Phi=\delta=0.$ This model has $M_{x}$ and $M_{y}$ mirror symmetries so we expect the quadrupole moment to be quantized, though that does not mean it takes the non-trivial topological value. Indeed this model is not expected to have a bulk quadrupole moment despite having insulating phases with fractional $e/2$ corner charge, since this charge is associated with edge polarization rather than a bulk quadrupole moment. Here we show a phase diagram for fixed $\lambda_{y}=1,\gamma_{x}=0,$ where we vary $\lambda_{x}$ and $\gamma_{y}.$ We show that our formula, while yielding a quantized value for $Q^{xy},$ does not distinguish between the insulating phases with and without fractional corner charge. This is a clear and important indication that our invariant is only sensitive to the bulk quadrupole moment, and is not just testing for the possible existence of corner charge alone.

Finally, in Fig. 2d we show a calculation for the octupole moment in the phase diagram of the mirror-symmetric octupole model of Ref. Benalcazar et al., 2017a where we have fixed the inter-cell hopping $\lambda$ and varied the intra-cell hoppings $\gamma_{x,y,z}$ in the interval $[0,2].$ We see that the phase diagram matches the expected result with a topological octupole moment $O^{xyz}=e/2$ when $|\gamma_{i}/\lambda|<1$ for all $i.$

V.2 Quadrupole model with $C_{4}T$ Symmetry

In this subsection we perform a calculation of the quadrupole moment in a 2D version of the chiral hinge-insulator model presented in Ref. Schindler et al., 2018. This model has a quantized quadrupole moment protected by the product of rotation and time-reversal symmetry $C_{4}T.$ The Bloch Hamiltonian we consider is

[TABLE]

where $\Gamma^{0}=\tau^{z}\otimes\mathbb{I},\Gamma^{1}=\tau^{x}\otimes\sigma^{x},\Gamma^{2}=\tau^{x}\otimes\sigma^{y},\Gamma^{3}=\tau^{x}\otimes\sigma^{z}.$ When $\Delta=0$ this model has both $C_{4}$ and time-reversal symmetry $T,$ and is a model representing a 2D time-reversal invariant topological insulator with helical edge states when, e.g., $0<m<2$. When $\Delta\neq 0$ the system breaks both $C_{4}$ and $T$ but preserves the product $C_{4}T.$ The non-vanishing $\Delta$ acts to gap the helical edge states and generate corner modes for an open square geometry that respects global $C_{4}$ symmetry. Thus we expect to find an ill-defined quadrupole moment when $\Delta=0,$ and a non-vanishing topological quadrupole moment $Q^{xy}=e/2$ when $0<m<2,\Delta\neq 0.$ When $m<0,\Delta\neq 0$ we expect to find a well-defined, but vanishing quadrupole moment. We show the calculated phase diagram as a function of $m$ and $\Delta$ in Fig. 3. We note that we look at a range of $\Delta$ that is small enough not to destroy the bulk topology and we see that our operator correctly reproduces the phase diagram.

VI Connection to recent work

There have been several recent papers discussing the evaluation of the quadrupole operators on which we can commentKang et al. (2018); Ono et al. (2019). Kang et. al.Kang et al. (2018) simultaneously proposed the same definitions for the quadrupole and octupole operators as we do in Eqs. 2-5. Similar to our argument in Eq. 16, they investigate the invariance of the operator $\hat{\mathcal{U}}_{xy}$ under translations, and conclude that the polarization must be zero for the operator to be well-defined. We have both pointed out that the ground state is not an eigenstate of polarization $\hat{\mathcal{U}}_{j}$, which is what complicates the evaluation of $\hat{\mathcal{U}}_{xy}.$ They point out a violation of translation invariance of $\mathcal{O}(1/E_{\rm gap})$, inversely proportional to the excitation gap:

[TABLE]

where $\hat{\mathcal{U}}_{xy}^{\prime}$ is the operator $\hat{\mathcal{U}}_{xy}$ translated in the $x$ direction by $x\rightarrow x+L_{x}$. We describe this same error as being a result of dipole fluctuations, which, for non-interacting systems, measure the Wannier localization length. For single-particle orbitals, the error disappears in the limit that the Wannier functions are $\delta$-function localized, equivalent to their limit of $E_{\rm gap}\rightarrow\infty$. We drew a connection between dipole fluctuations and decreasing magnitude of $\langle\hat{\mathcal{U}}_{xy}\rangle,$ which is a signature of the problem with translation invariance. In contrast to their discussion, we point out that zero polarization is not sufficient to have a well-defined quadrupole in finite-sized periodic systems. As such, for systems with dipole fluctuations, evaluating this operator in PBCs necessarily uses an approximate form; we described one approximation that we used in our calculation, which is similar to that used in Ref. Kang et al., 2018.

To confirm the definitions of the quadrupole operator $\hat{\mathcal{U}}_{xy}$ as a reliable estimator of the response to electric field gradients, they also give an argument similar to our adiabatic evolution discussion, but in a more field-theoretical approach.

In response to these proposals, Ono et. al.Ono et al. (2019) comment on some problems that arise when applying the quadrupole operator $\hat{\mathcal{U}}_{xy}$ to certain tight-binding models, that were not considered in our workOno et al. (2019). They find that their evaluation of the quadrupole operator has more serious issues than those discussed here, namely that the result seems to depend non-trivially on the correlation length in an insulator phase when the quadrupole moment should nominally be quantized. They also comment that their results can depend on the parity of the number of lattice sites, but as we show in Appendix A, C this latter issue has a simple origin, and similar effects can appear even for the polarization, and even for classical point-charge configurations.

In the models in which we evaluated the quadrupole moment we did not find any dependence on the correlation length besides typical finite-size effects, so it might be useful to comment on the distinctions between the models we considered and one of their models, though fully resolving the discrepancy will require further research. One of their models is a $C_{4}$ invariant system with four occupied bands that form a Wannier representable obstructed atomic limit. The electrons in the model sit at maximal Wyckoff positions $1a,1b,2c$, i.e., one at the center, one at the corner, and two in the middle of the edges of the unit cell. One possibly important distinction is that this model does not have gapped Wannier bands. While the need for gapped Wannier bands was emphasized in Ref. Benalcazar et al., 2017b, it has not been proven that this is a generic requirement to define the quadrupole moment (though a gapped, neutral edge is required). A second distinction is that in the models we consider, one can always relate the configurations with a non-vanishing quadrupole to a configuration representing a trivial onsite limit through a continuous deformation that preserves a symmetry under which the dipole moment does not change. Thus, in our models one can clearly calculate a change in quadrupole moment when going from the trivial atomic limit to a non-vanishing quadrupole configuration. In their model, such an interpolation does not seem to exist, so it is not clear how to compare the quadrupole moment in their obstructed atomic limit to the trivial atomic limit. While again, this may not be a requirement to define a quadrupole model, it does seem like a natural consideration.

Ref. Ono et al., 2019 also makes a note regarding second-order contributions to the perturbative expansion of the magnitude $|\langle\Phi_{0}|\hat{\mathcal{U}}_{x}|\Phi_{0}\rangle|$, and comments that Resta’s formula for polarization is not guaranteed to work, i.e., the magnitude may go to zero in higher dimensions. This is a well-known fact for polarization. Indeed to obtain a nonzero magnitude of the expectation value of $\hat{\mathcal{U}}_{x}$ in higher dimensions one must take the correct thermodynamic limit, namely, where the direction parallel to the polarization component of interest goes to infinity first. We also find that, for periodic boundary conditions, the quadrupole magnitude $|\langle\hat{\mathcal{U}}_{xy}\rangle|$ goes to zero in the thermodynamic limit if dipole fluctuations are finite, seemingly no matter what thermodynamic limit one uses. Despite this we find that this issue does not affect the consistency of the phase in the tight-binding calculations we performed in finite size systems. Based on the calculations in Appendix C, we expect that if we truncate the system, as indicated by our approximation scheme, then the magnitude $|\langle\hat{\mathcal{U}}_{xy}\rangle|$ can approach a finite (though likely not unity) value if we take the limit $L_{x}=L_{y}\to\infty,$ even in the presence of dipole fluctuations.

In summary, Ono et. al. raise some interesting questions about the evaluation of the many-body quadrupole operator and about what constraints are necessary to enforce in order to have a well-defined quadrupole moment in a crystalline system. The issues they have pointed out will be valuable in specifying careful definitions and treating boundary effects correctly for detecting higher order multipole moments in bulk systems.

VII Conclusion

In conclusion, we have proposed definitions for many-body operators whose expectation values determine the quadrupole and octupole moments of insulators. We showed that the change in quadrupole moment corresponds to dipole currents in the material. We proposed a method to evaluate the quadrupole; however, the method is not entirely satisfactory because it suffers from significant finite size errors related to the fluctuations of the dipole moment, even when the total dipole moment is vanishing. However, we showed that on several nontrivial tightbinding models the operators do indeed capture the bulk properties, and can be used to evaluate topological indices. The operator can be evaluated using many-body wave functions as well as using tightbinding wave functions, so it can be used to generate new many-body topological indices for systems with interactions and/or disorder.

Note: During the preparation of this manuscript we became aware of an independent overlapping work by B. Kang, K. Shiozaki, and G. Y. ChoKang et al. (2018). We thank them for discussions and for coordinating submission.

Acknowledgements.

We thank Jahan Claes, Matthew Foulkes, and Gerardo Ortiz for useful conversations. LKW would like to thank the Simons Foundation Collaboration on the many electron problem for support. TLH thanks the US National Science Foundation under grant DMR 1351895-CAR, and the MRSEC program under NSF Award Number DMR-1720633 (SuperSEED) for support.

Appendix A Constraints on Charge Configurations for the Quadrupole and Octupole Operators

In this Appendix we will evaluate the conditions under which one will generically find a non-vanishing expectation value of multipole operators if the electron orbitals are point-charge localized. For example, in a limit where all dipole fluctuations vanish, we expect that the expectation values of the quadrupole operators are well-behaved in extended systems. More precisely, to evaluate the quadrupole moment we want to consider a many-body state that is an exact eigenstate of the $\mathcal{U}_{x}$ and $\mathcal{U}_{y}$ operators. This is not necessarily an unphysical situtation, as we know it will be approximately true as we approach the thermodynamic limit. That is, if we consider our insulating system to be large compared to the characteristic size of individual Wannier functions of the occupied bands, then we can treat the charge distribution as a collection of point charges located at the Wannier centers that are periodically repeated throughout the crystal latticeVanderbilt and King-Smith (1993). For a many-body state represented by a product of single-particle orbitals, vanishing dipole fluctuations imply that the orbitals are $\delta$-function localized, so we will consider the expectation values of the multipole operators in such a state. Even in this situation, there are constraints under which the expectation values of the quadrupole or octupole operators are well-defined in the thermodynamic limit, i.e, when they effectively transform in the trivial representation of the translation group, for such point-charge product states. We will determine those conditions now.

For simplicity, let us consider our electron configuration to be a rectangular lattice of point charges and focus on the quadrupole moment $Q^{xy}$ associated to the operator

[TABLE]

Now let us transform this operator by a translation ${\bf{R}}=h_{1}a_{x}\hat{x}+h_{2}a_{y}\hat{y}$ where $h_{1},h_{2}$ are integers. We find that

[TABLE]

where the sum over $\alpha$ runs over all the electrons in a single unit cell and $\nu$ is the electron filling.

Ultimately we want this operator to be invariant under any lattice translation of all the electrons, and the strongest constraints arise from taking $h_{1},h_{2}$ to be the smallest non-zero integers, e.g., $h_{1}=h_{2}=1.$ For this choice we find the constraint

[TABLE]

To satisfy this constraint in a manner that is independent of the eveness/oddness of $N_{x},N_{y}$ it is sufficient to choose $\nu=2\mathbb{Z}$ and the $x$ and $y$ components of the dipole moment per unit cell must be integer multiples of $a_{x},a_{y}$ respectively. We can relax the constraint of $\nu$ being integer valued if we modify our starting operator Eq. 24 according to a prescription analogous to that presented in Ref. Aligia and Ortiz, 1999, but in that case we will still require $\nu=2\mathbb{Q}.$ While we have derived these constraints assuming the charge distribution is point-like we expect this result to apply to more complicated distributions in the thermodynamic limit where the system size is much larger than the spatial extent of the Wannier functions that compose our insulating ground state. The coordinates $\bar{x}_{\alpha},\bar{y}_{\alpha}$ will then be associated to the centers of the Wannier functions in each unit cell.

We can repeat this calculation for the octupole moment, e.g, $O^{xyz},$ and translating by a vector ${\bf{R}}=a_{x}\hat{x}+a_{y}\hat{y}+a_{z}\hat{z}.$ We find the constraint that the dipole and quadrupole moments per unit cell must be integer valued as well as the constraint

[TABLE]

In order for this constraint to be satisfied independently of the system size we thus need $\nu=4\mathbb{Z}$ for the octupole moment.

Appendix B Finite size effects of the quadrupole operator

The evaluation of the expectation value of the quadrupole operator $\hat{\mathcal{U}}^{Q}_{xy}$ (Eq. 24) exhibits some challenges. We show how this manifests in our calculation as a function of system size. Fig. 4 shows the $\gamma_{x}$-$\gamma_{y}$ phase diagram of the quantized, mirror-symmetric quadrupole model Eq. 22 from Ref. Benalcazar et al., 2017b for $L_{x}=L_{y}=10,\;20,\;40$. There is an error in a thin region near the phase boundary, where only one of $\gamma_{x}$ and $\gamma_{y}$ is close to (but less than) one. The expectation value predicts the wrong topological phase in this region; however, this region decreases rapidly with system size. In the thermodynamic limit, the error disappears. This type of error is due to the large correlation length that exceeds the system size for parameters near the phase boundaries.

The expectation value of the operator $\hat{\mathcal{U}}^{Q}_{xy}$ has a phase and a magnitude. While the phase gives the quadrupole moment, the magnitude reflects its fluctuations. Fig. 5 shows how the magnitude depends on our finite simulation size for selected points in the $\gamma_{x}$-$\gamma_{y}$ phase diagram shown in Fig. 4. The magnitudes are close to zero and decrease as the system gets larger. This does not mean that the physical quantity is not meaningful; in comparison, for a 3D system, Resta’s polarization operator also decreases to zero magnitude if $N_{x}=N_{y}=N_{z}\rightarrow\infty$. Exploring efficient ways to evaluate the expectation value for large systems is still an open challenge. In Appendix C we will show how this magnitude decrease can be traced back to fluctuations in the dipole moments.

Appendix C Scaling of the Magnitude of the Quadrupole Expectation Value

In this Appendix we discuss the scaling properties of the expectation values of the dipole and quadrupole operators using a test configuration of continuum charge distributions arranged on a crystalline lattice. We will consider the cases of product states of $\delta$-function point charges as well as Gaussian charges, both product and determinant states. In these systems, the dipole fluctuations cause the magnitude of the quadrupole operator to go to a limiting value which is not equal to one.

C.1 Polarization

For the calculation of the polarization let us consider a 1D system with a filling of one electron per unit cell at a position $x_{0}$ with respect to the origin of each unit cell. The electrons in our system are thus located at the positions $x_{0}+n$ in units of the lattice constant $a$ which we set to $a=1$, and where $n\in\mathbb{Z}.$ Now we can calculate the polarization via the expectation value of the operator $\hat{\mathcal{U}}_{x}=\exp\left[\frac{2\pi i\hat{X}}{L_{x}}\right]$ in the many-body ground state consisting of a tensor product of localized $\delta$-function charges. The polarization is given by

[TABLE]

where $\mathcal{V}$ is the volume of the system, i.e., $L_{x}=N_{x}$ in 1D. After integrating over the real-space coordinates of each electron to generate the expectation value, we find the result

[TABLE]

which is only well-defined mod $e,$ so the second term either contributes [math] or $1/2$ depending on the parity of $N_{x}.$

Now let us take localized charges obeying a Gaussian distribution

[TABLE]

We find that, for a single charge at $x_{0}$

[TABLE]

We see that this result includes a phase factor encoding the position $x_{0}$ as well as an exponential factor with magnitude $\leq 1.$ The phase factor is identical to what the $\delta$-function charge distribution produces, and when all electrons across the lattice are considered it will exactly reproduce Eq. 29 for the value of the polarization. On the other hand, if we have a number of electrons $N_{e}=N_{x}$ then the magnitude of the expectation value will be $e^{-\pi^{2}\sigma_{x}^{2}/N_{x}}$ which will tend to unity as $N_{x}\to\infty$ as expected. Thus, this operator recovers unitarity in the thermodynamic limit.

Let us now illustrate an important subtlety by reconsidering the case of Gaussian charges in a 2D system. Suppose that we have one electron per unit cell on a square lattice, and that the electrons are centered at $(x_{0},0)+{\bf{R}}$ where ${\bf{R}}$ is the set of 2D lattice vectors. For a single-electron, which is Gaussian distributed in $x$ and $y$ we have

[TABLE]

which is the same as the 1D case above. The distinction between dimensions appears when we use the fact that $N_{e}=N_{x}N_{y}.$ After including all of the electrons we have

[TABLE]

The important feature is that the magnitude of this operator approaches unity in a specific thermodynamic limit where we take $N_{x}\to\infty$ before $N_{y},$ but it vanishes if we use the opposite order, and approaches a constant $<1$ if we take the limit while keeping a fixed aspect ratio $N_{x}/N_{y}$. Under the same conditions, but in 3D, we would find a factor $e^{-\pi^{2}\sigma_{x}^{2}N_{y}N_{z}/N_{x}}.$ Thus, generically to keep the operator expectation value non-vanishing we must always take the thermodynamic limit in the direction of the polarization component of interest before letting the transverse directions approach infinity.

C.2 Quadrupole Moment

Now let us consider the quadrupole moment for a 2D system on a square lattice. We can begin with the periodic arrangement of $\delta$-function point charges and calculate the expectation value of $\hat{\mathcal{U}}_{xy}=\exp\left[\frac{2\pi i\hat{Q}_{xy}}{L_{x}L_{y}}\right].$ We must be careful in this case to specify that the total dipole moment vanishes (up to an integer). We will satisfy this constraint by taking the system with two electrons per cell with coordinates $(x_{0},y_{0})$ and $(-x_{0},-y_{0})$ with respect to the origin of the cell, and we repeat these positions across the entire lattice. For this configuration, and for a lattice indexed from $1\dots N_{x},1\ldots N_{y},$ we find

[TABLE]

which depends on the parity of $N_{x}$ and $N_{y}$ similarly to the polarization.

Now let us move on to charges with a Gaussian distribution described by

[TABLE]

In this case, the expectation value generated by a single electron with a Gaussian wave function centered at $(x_{0},y_{0})$ is

[TABLE]

where $\alpha=2\pi/L_{x}L_{y}.$ We can see from this result that both the magnitude and the phase of this expectation value depend on the system size when $\sigma_{x}$ and/or $\sigma_{y}$ are non-vanishing. This is not the case for the polarization, where the imaginary phase is independent of the system size for a Gaussian distribution. Here the difference arises because we allow for fluctuations of the dipole moment, which is equivalent to having non-vanishing $\sigma_{i}.$ The dipole moment is vanishing on average, but has non-zero fluctuations that lead to this ambiguity in the phase factor. The same would likely occur for the polarization if we let the charge in each unit cell fluctuate, but we have implicitly assumed that it does not fluctuate. We can calculate the quadrupole moment density from the imaginary part of the $\log$ of this expression for the expectation value to find

[TABLE]

for a single electron. We find that this value converges to the correct result as $L_{x},L_{y}\to\infty,$ i.e., when the system size is much larger than the dipole fluctuation lengths. When all of the electrons in the lattice are taken into account this will reproduce the value calculated for the $\delta$-function distribution.

Finally, let us consider the magnitude $\Lambda$ of this expectation value. We find a magnitude for a single electron to be

[TABLE]

We want to determine the value of this magnitude for a collection of electrons on the lattice. The magnitude for the full ground state is harder to calculate since it depends on the positions of the electrons, but for a product state including two electrons per unit cell, i.e., $N_{e}=2N_{x}N_{y},$ and ignoring contributions of the charge density that spreads from the last cells back into the first cells due to periodic boundary conditions, we find

[TABLE]

where in the first approximation we have dropped the terms that shift $N_{x,y}$ by $1$ coming from a discrete sum over unit cells, and in the second approximation we have dropped the terms proportional to $x_{0},y_{0}$ since they are coordinates within a unit cell and are less than $1$ in units of the lattice constant. We find that as $N_{x},N_{y}\to\infty$ that the first factor tends to unity:

[TABLE]

If we make an analogy with the polarization, this factor has a dependence on the quadrupole fluctuations $\sigma_{x}^{2}\sigma_{y}^{2}/4$ and tends to unity as long as those fluctuations are finite.

Now let us consider the second factor

[TABLE]

This factor has interesting features: (i) the numerator has a dependence on the dipole fluctuations in the $x$ and $y$ directions, (ii) the denominator has a dependence on the quadrupole fluctuations. If dipole fluctuations in either direction vanish, then there is a consistent way to take the thermodynamic limit such that this term tends to unity. If the dipole fluctuations are non-vanishing in both directions then the best that one can do is to take $N_{x}=N_{y}\to\infty.$ In this case the factor tends to the finite value $e^{-2\pi^{2}(\sigma_{y}^{2}+\sigma_{x}^{2})/3},$ which is non-vanishing, but not unity. We also note that if one takes $N_{x}$ or $N_{y}$ to infinity first then this factor generically tends to zero. This result seems to imply that even though $\hat{\mathcal{U}}_{xy}$ may not be strictly a well-defined operator in periodic systems with dipole fluctuations, its expectation value is still meaningful and can recover the correct results for the quadrupole moment.

The above discussion considers a product state of Gaussian orbitals; however, the tightbinding calculations presented in the paper use a Slater determinant, so we also show a comparison with a determinant of the Gaussian orbitals. An analytic calculation of the determinant state is impractical, so we compared the numerical results of the determinant and product states.

We confirm that the determinant state yields similar results to the product state except when the Gaussian spreads $\sigma_{x},\sigma_{y}$ grow so large that determinant state exhibits a transition to a conducting state (i.e., the magnitude of the polarization operator goes to zero in the infinite limit), while the product state remains insulating (Fig. 6). For smaller values of $\sigma_{x},\sigma_{y}$, the phases and magnitudes of the product state agree with the determinant state. We confirm that the phase of the polarization operator expectation value in our states is always zero, so the quadrupole is well-defined in principle. We also confirm that the magnitude of the quadrupole operator remains fixed when we vary $\sigma=\sigma_{x}=\sigma_{y}$ and even values of $N_{x}=N_{y}$, and is the same for the determinant and product states. When $\sigma$ increases, even when the magnitude of the polarization goes to one, the limit of the quadrupole magnitude decreases (it does not always approach one, even when finite). Beyond a certain value of $\sigma$, the quadrupole magnitude goes to zero. Based on our numerics, the insulating-conducting transition of the determinant state seems to coincide with the transition where the quadrupole magnitude goes to zero.

Our simple model of Gaussian orbitals can only go so far in explaining the results from our tightbinding calculations; however, we have shown that even with very simple wave functions, the quadrupole magnitude is suppressed by dipole fluctuations, and does not generally tend to one in the thermodynamic limit. We expect that our results would closely match the tightbinding results if we accounted for dipole fluctuations across the periodic boundary conditions.