Difficulties in operator-based formulation of the bulk quadrupole moment

TL;DR

This paper discusses the challenges in formulating the bulk quadrupole moment in insulators, highlighting issues with recent operator-based approaches and emphasizing the complexity of extending dipole formulations to higher multipoles.

Contribution

The paper critically analyzes recent proposals for operator-based formulations of the bulk quadrupole moment, identifying key issues and limitations.

Findings

Recent operator-based formulas for quadrupole moments have unresolved issues.

The extension from dipole to quadrupole moments is non-trivial and problematic.

The paper clarifies the fundamental difficulties in defining higher-order multipole moments.

Abstract

Electric multipole moments are the most fundamental properties of insulating materials. However, the general formulation of bulk multipoles has been a long standing problem. The solution for the electric dipole moment was provided decades ago by King-Smith, Vanderbilt, and Resta. Recently, there have been attempts at generalizing Resta's formula to higher-order multipoles. We point out several issues in the recent proposals.

Figures

Tables

Table 1. Table 1: The analytic expression of Δ ​ q ~ i ​ j Δ subscript ~ 𝑞 𝑖 𝑗 \Delta\tilde{q}_{ij} for the tight-binding model in Eq. ( 10 ) in the limit of large L 𝐿 L .

$(x_1,L)$	$(1,\text{even})$	$(1,\text{odd})$	$(-\frac{L-1}{2},\text{odd})$
$\Delta {\tilde{q}}_{xx}$, $\Delta {\tilde{q}}_{yy}$	$\frac{1}{4}$	$\frac{1}{4}$	$\frac{3}{4}$
$\Delta {\tilde{q}}_{xy}$, $\Delta {\tilde{q}}_{yx}$	$\frac{1}{2}+\frac{1}{\pi }\arctan \frac{2}{\pi }$333$=0.680454\cdots$	$\frac{1}{2}$	$\frac{3}{4}$

Table 2. Table 2: Analytic expression of Δ ​ q ~ x ​ y Δ subscript ~ 𝑞 𝑥 𝑦 \Delta\tilde{q}_{xy} for the first model in the limit of large L 𝐿 L .

$(x_1,L)$	$(1,\text{even})$	$(1,\text{odd})$	$(-\frac{L-1}{2},\text{odd})$
$\Delta {\tilde{q}}_{xy}^{\text{(a)}}$	$\frac{1}{4}+\frac{1}{\pi }\arctan \frac{2}{\pi }$555$=0.430454\cdots$	$\frac{3}{4}$	$0$
$\Delta {\tilde{q}}_{xy}^{\text{(b)}}$	$\frac{3}{4}$	$\frac{1}{4}$	$\frac{1}{4}$

Table 3. Table 3: Analytic expression of Δ ​ q ~ x ​ y Δ subscript ~ 𝑞 𝑥 𝑦 \Delta\tilde{q}_{xy} for the second model in the limit of large L 𝐿 L .

$(x_1,L)$	$(1,\text{even})$	$(1,\text{odd})$	$(-\frac{L-1}{2},\text{odd})$
$\Delta {\tilde{q}}_{xy}^{\text{(a)}}{}^{\prime }$	$0$	$0$	$0$
$\Delta {\tilde{q}}_{xy}^{\text{(b)}}{}^{\prime }$	$\frac{1}{2}$	$0$	$0$

Table 4. Table 4: The analytic expression of Δ ​ q ~ x ​ y Δ subscript ~ 𝑞 𝑥 𝑦 \Delta\tilde{q}_{xy} for the stacked model in the limit of large L 𝐿 L .

$(x_1,L)$	$(1,\text{even})$	$(1,\text{odd})$	$(-\frac{L-1}{2},\text{odd})$
$\Delta {\tilde{q}}_{xy}^{\text{(a)}}+\Delta {\tilde{q}}_{xy}^{\text{(a)}}{}^{\prime }$	$\frac{1}{4}+\frac{1}{\pi }\arctan \frac{2}{\pi }$	$\frac{3}{4}$	$0$
$\Delta {\tilde{q}}_{xy}^{\text{(b)}}+\Delta {\tilde{q}}_{xy}^{\text{(b)}}{}^{\prime }$	$\frac{1}{4}$	$\frac{1}{4}$	$\frac{1}{4}$

Equations

$$n^{\text{(OBC)}}\equiv \frac{1}{V}\sum _{\mathbf{x}}⟨{\Phi }_0|{\hat{n}}_{\mathbf{x}}|{\Phi }_0⟩,$$

$$p_i^{\text{(OBC)}}\equiv \frac{1}{V}\sum _{\mathbf{x}}{x}_i⟨{\Phi }_0|{\hat{n}}_{\mathbf{x}}|{\Phi }_0⟩,$$

$$q_{ij}^{\text{(OBC)}}\equiv \frac{1}{V}\sum _{\mathbf{x}}{x}_i{x}_j⟨{\Phi }_0|{\hat{n}}_{\mathbf{x}}|{\Phi }_0⟩.$$

$$p_i^{\text{(Resta)}}\equiv \frac{1}{2\pi }\text{Im}\ln ⟨{\Phi }_0|{\hat{U}}_i|{\Phi }_0⟩\mathrm{mod}1,$$

$${\hat{U}}_i\equiv e^{2\pi i{\hat{P}}_i/L},{\hat{P}}_i\equiv \sum _{\mathbf{x}}{x}_i{\hat{n}}_{\mathbf{x}}.$$

$${\tilde{q}}_{ij}\equiv \frac{1}{2\pi }\text{Im}\ln ⟨{\Phi }_0|{\hat{U}}_{ij}|{\Phi }_0⟩\mathrm{mod}1,$$

$${\hat{U}}_{ij}\equiv e^{2\pi i{\hat{Q}}_{ij}/L^2},{\hat{Q}}_{ij}\equiv \sum _{\mathbf{x}}{x}_i{x}_j{\hat{n}}_{\mathbf{x}}.$$

$$\psi (\mathbf{x}+L{\mathbf{e}}_i)=\psi (\mathbf{x}).$$

$${\hat{U}}_{ij}|\psi ⟩=\sum _{\mathbf{x}}{e}^{2\pi i{x}_i{x}_j/L^2}\psi (\mathbf{x})c_{\mathbf{x}}^{†}|0⟩$$

$$H_{\mathbf{k}}=-\left(\begin{array}{cccc}0& e^{i{k}_x}& 0& e^{i{k}_y}\\ e^{-i{k}_x}& 0& e^{i{k}_y}& 0\\ 0& e^{-i{k}_y}& 0& e^{-i{k}_x}\\ e^{-i{k}_y}& 0& e^{i{k}_x}& 0\end{array}\right).$$

$$U_{C_4}{H}_{k_x,k_y}=H_{-k_y,k_x}{U}_{C_4},U_{C_4}\equiv \left(\begin{array}{cccc}0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right).$$

$${\hat{\gamma }}_{\mathbf{R}}^{†}\equiv \frac{1}{2}({\hat{c}}_{\mathbf{R}a}^{†}+{\hat{c}}_{\mathbf{R}+{\mathbf{e}}_x,b}^{†}+{\hat{c}}_{\mathbf{R}+{\mathbf{e}}_x+{\mathbf{e}}_y,c}^{†}+{\hat{c}}_{\mathbf{R}+{\mathbf{e}}_y,d}^{†}).$$

$${\tilde{q}}_{ij,\text{ions}}=\frac{1}{L^2}\sum _{\mathbf{R}}\left(R_i+\frac{1}{2}\right)\left(R_j+\frac{1}{2}\right).$$

$$\Delta p_i^{\text{(Resta)}}=p_i^{\text{(Resta)}}-p_{i,\text{ions}}^{\text{(Resta)}}=0,$$

$$\Delta {\tilde{q}}_{ij}={\tilde{q}}_{ij}-{\tilde{q}}_{ij,\text{ions}}.$$

$${\theta }_{\mathbf{x}}=\frac{m_{x}x+m_{y}y}{L}+{\vartheta }_{\mathbf{x}}(m_i\in \mathbb{Z}),$$

$$C_{\vec{x}}=$$

$${\int }_{\text{BZ}}\frac{d^{2}k}{(2\pi {)}^2}\sum _{n=1}^{N_{\text{occ}}}i⟨u_{\mathbf{k}n}|{\partial }_{k_i}{u}_{\mathbf{k}n}⟩$$

$${\int }_{\text{BZ}}\frac{d^{2}k}{(2\pi {)}^2}\sum _{n=1}^{N_{\text{occ}}}⟨{\partial }_{k_i}{u}_{\mathbf{k}n}|{\partial }_{k_j}{u}_{\mathbf{k}n}⟩.$$

$$|u_{\mathbf{k}{n}^{\prime }}{⟩}^{\prime }=\sum _{n=1}^{N_{\text{occ}}}|u_{\mathbf{k}n}⟩(w_{\mathbf{k}}{)}_{n,n^{\prime }}.$$

$$T_{ij}\equiv {\int }_{\text{BZ}}\frac{d^{2}k}{(2\pi {)}^2}\sum _{n=1}^{N_{\text{occ}}}⟨{\partial }_{k_i}{u}_{\mathbf{k}n}|(1-P_{\mathbf{k}})|{\partial }_{k_j}{u}_{\mathbf{k}n}⟩.$$

$${\hat{U}}_x|{\Phi }_0⟩\doteq e^{i\gamma }\left(|{\Phi }_0⟩+\frac{2\pi }{L_x}\sum _{N>0}|{\Phi }_N⟩\frac{⟨{\Phi }_N|{\hat{J}}_x|{\Phi }_0⟩}{E_N-E_0}+\cdots \right),$$

$$|⟨{\Phi }_0|{\hat{U}}_x|{\Phi }_0⟩|\doteq 1+O(L_x^{-2}),$$

$$({\hat{U}}_x{)}_{{\mathbf{k}}^{\prime }{n}^{\prime },\mathbf{k}n}={\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}+(2\pi /L_x){\mathbf{e}}_x}(B_{\mathbf{k}}{)}_{n^{\prime }n},$$

$$(B_{\mathbf{k}}{)}_{n^{\prime }n}$$

$$|⟨{\Phi }_0|{\hat{U}}_x|{\Phi }_0⟩|=e^{-2{\pi }^2{G}_{xx}V/L_x^2+O(V/L_x^3)},$$

$$G_{ij}\equiv \text{Re}\int \frac{d^{d}k}{(2\pi {)}^d}(g_{\mathbf{k}}{)}_{ij},$$

$$(g_{\mathbf{k}}{)}_{ij}=\sum _{n=1}^{N_{\text{occ}}}⟨{\partial }_{k_i}{u}_{\mathbf{k}n}|(1-P_{\mathbf{k}})|{\partial }_{k_j}{u}_{\mathbf{k}n}⟩.$$

$${\tilde{q}}_{ij}=\frac{1}{2\pi }\sum _{\mathbf{R}}\text{Im}\ln ({\hat{U}}_{ij}{)}_{\mathbf{R},\mathbf{R}}.$$

$$({\hat{U}}_{xx}{)}_{\mathbf{R},\mathbf{R}}=\frac{1}{2}\times \{\begin{array}{ll}{F}_{R_x^2}+F_{(R_x+1{)}^2}& \text{for }{R}_x<x_L\\ F_{x_L^2}+F_{x_1^2}& \text{for }{R}_x=x_L\end{array},$$

$$({\hat{U}}_{yy}{)}_{\mathbf{R},\mathbf{R}}=\frac{1}{2}\times \{\begin{array}{ll}{F}_{R_y^2}+F_{(R_y+1{)}^2}& \text{for }{R}_y<y_L\\ F_{y_L^2}+F_{y_1^2}& \text{for }{R}_y=y_L\end{array},$$

$$({\hat{U}}_{xy}{)}_{\mathbf{R},\mathbf{R}}=({\hat{U}}_{yx}{)}_{\mathbf{R},\mathbf{R}}=\frac{1}{4}\times \{\begin{array}{ll}{F}_{R_x{R}_y}+F_{(R_x+1)R_y}+F_{(R_x+1)(R_y+1)}+F_{R_x(R_y+1)}& \text{for }{R}_x<x_L\text{ and }{R}_y<y_L\\ F_{x_L{R}_y}+F_{x_1{R}_y}+F_{x_1(R_y+1)}+F_{x_L(R_y+1)}& \text{for }{R}_x=x_L\text{ and }{R}_y<y_L\\ F_{R_x{y}_L}+F_{(R_x+1)y_L}+F_{(R_x+1)y_1}+F_{R_x{y}_1}& \text{for }{R}_x<x_L\text{ and }{R}_y=y_L\\ F_{x_L{y}_L}+F_{x_1{y}_L}+F_{x_1{y}_1}+F_{x_L{y}_1}& \text{for }{R}_x=x_L\text{ and }{R}_y=y_L\end{array}.$$

$${\tilde{q}}_{xy,\text{ions}}$$

$${\tilde{q}}_{xy}^{\text{(a)}}=\frac{1}{2\pi }\sum _{\mathbf{R}}\text{Im}\ln ({\hat{U}}_{xy}^{\text{(a)}}{)}_{\mathbf{R},\mathbf{R}},$$

Full text

I Introduction

Electric multipole moments are of fundamental importance in understanding the property of insulators. Most significantly, they characterize the charge distribution near boundaries of a finite sample. For two dimensional systems, boundaries can be either one-dimensional edges or zero-dimensional corners. A nonzero electric polarization results in nonzero charge density for edges. Recently, so-called quadrupole insulators that feature charged corners instead of edges have attracted extensive research interest because of their relation to higher-order topological insulators. Parameswaran and Wan (2017); Schindler et al. (2018); Peng et al. (2017); Langbehn et al. (2017); Song et al. (2017); Benalcazar et al. (2017a, b); Fang and Fu (2017); Ezawa (2018); Shapourian et al. (2018); Zhu (2018); Yan et al. (2018); Wang et al. (2018a, b); Khalaf et al. (2018); Khalaf (2018); Trifunovic and Brouwer (2019); Okuma et al. (2018) To develop a systematic understanding of quadrupole insulators, we need to establish a general framework that allows us to compute the quadrupole moment of insulators.

Despite their importance, the precise formulation of multipole moments is known to be a difficult task. For finite systems under the open boundary condition (OBC), multipole moments can be simply defined by the classical formula Resta and Vanderbilt (2007); Lines and Glass (2001); Kittel (2004):

[TABLE]

Here, $|\Phi_{0}\rangle$ is the many-body ground state, $\hat{n}_{\bm{x}}$ is the number density operator at the position $\bm{x}$, and $V=L^{2}$ is the volume of the system. (Throughout this work, we consider a two dimensional square-shaped system for brevity.) The polarization $p_{i}^{\text{(OBC)}}$ is independent of the arbitrary choice of the origin only in neutral systems ($n^{\text{(OBC)}}=0$). Similarly, the quadrupole moment $q_{ij}^{\text{(OBC)}}$ is well-defined only when all of $p_{i}^{\text{(OBC)}}$ ($i=x,y$) and $n^{\text{(OBC)}}$ vanish.

The situation gets a lot more complicated when considering extended systems under a periodic boundary condition (PBC), since the position operator in the above expressions becomes ill-defined. Martin (1974); Aligia and Ortiz (1999) For band insulators, King-Smith and Vanderbilt formulated the bulk polarization in terms of the Berry phase of Bloch wavefunctions.Vanderbilt and King-Smith (1993); King-Smith and Vanderbilt (1993) The Berry phase approach was generalized to many-body systems with interactions and/or disorders by replacing the single-particle crystal momentum to the twisted angle of the boundary condition.Ortiz and Martin (1994); Souza et al. (2000); Watanabe and Oshikawa (2018)As an alternative formulation, Resta Resta (1998) proposed the following formula of the bulk electric polarization in many-body setting under a PBC:

[TABLE]

According to Ref. Resta, 1998, this relation holds even in the presence of disorders and many-body interactions as long as the excitation gap is non-vanishing.

There have been several recent proposals on how to compute the bulk quadrupole moment. For band insulators, the nested-Wilson loop approach formulated in Refs. Benalcazar et al., 2017a, b aims at providing a way of computing the bulk contribution to the quantized corner charge, under assumptions of the spatial symmetry and the so-called “Wannier gap”. We remark here that, although the nested-Wilson loop approach gives a topological invariant, this invariant is generally not a bulk topological invariant—according to Ref. Benalcazar et al., 2017b the nested Wilson loop invariant can change its value if the Wannier gap is closed while the bulk band gap and the protecting symmetry are maintained. Yet, the claim of Refs. Benalcazar et al., 2017a, b is that combining the “bulk” contribution to the edge polarization $p_{x}^{\text{edge}}$, $p_{y}^{\text{edge}}$ obtained this way with an independent input on the corner charge $Q^{\text{corner}}$ computed under a certain open boundary condition, one gets the bulk quadrupole moment via $q_{xy}=Q^{\text{corner}}-p_{x}^{\text{edge}}-p_{y}^{\text{edge}}$. Benalcazar et al. (2017a, b) As later pointed out by us, Trifunovic et al. (2019) decorations by polarized one-dimensional chains change the corner charge; the boundary Hamiltonian and thus the boundary polarization are not generally defined quantities. Therefore it is not possible to subtract the contribution from edge polarization.

As a more general definition of the bulk quadrupole moment in many-body systems under a PBC, two independent groups Wheeler et al. (2018); Kang et al. (2018) proposed a possible generalization of Resta’s formula to higher-order multipoles. For example, their formulas for quadrupole moments read

[TABLE]

Arguments supporting Eq. (6) are based on a field theoretical calculation of the bulk response against a non-uniform electric field Kang et al. (2018) and a “perturbation” theory 111We explain why this perturbation theory fails in Appendix A. expanding the effect of the operator $\hat{U}_{ij}$ in the series of $x_{i}x_{j}/L^{2}$.Wheeler et al. (2018) The authors of Refs. Kang et al., 2018; Wheeler et al., 2018 have also provided a numerical proof of Eq. (6) using tight-binding models of quadrupole insulators. More recently, Refs. Agarwala et al. (2019); Lin et al. (2019) employed the formula (6) in their study of quadrupole insulators.

The goal of this work is to point out issues in the formula (6) for the bulk quadrupole moment. A satisfactory formulation of the bulk quadrupole moment would fulfill the following requirements: (i) independence from the choice of origin and the period $L$ and (ii) quantization in the presence of sufficiently large point group symmetries so that it can serve as a topological invariant characterizing the quantized corner charge of quadrupole insulators. However, we discuss that the operator $\hat{U}_{ij}$ in Eq. (7) is inconsistent with the assumed PBC, and, as a consequence, $\tilde{q}_{ij}$ in Eq. (6) meets none of these criteria.

II Problems in the proposed formula

II.1 Violation of the periodicity

We first point out an obvious issue in Eqs. (6) and (7). Let us consider a single-particle state $|\psi\rangle=\sum_{\bm{x}}\psi(\bm{x})c_{\bm{x}}^{\dagger}|0\rangle$. The PBC requires that the wavefunction $\psi(\bm{x})$ has the periodicity in $\bm{x}$ with the period $L$ in both the $x$ and $y$ directions:

[TABLE]

where $\bm{e}_{i}$ represents the unit vector along $i=x,y$. For example, the wavefunction of the Bloch state $\psi_{\bm{k}n}(\bm{x})=(1/\sqrt{V})u_{\bm{k}n}(\bm{x})e^{i\bm{k}\cdot\bm{x}}$ fulfills this condition thanks to the quantization of $\bm{k}$ to the integer multiples of $2\pi/L$.

The operator $\hat{U}_{i}$ in Eq. (5) preserves this periodicity. The wavefunction of the state $\hat{U}_{i}|\psi\rangle$ is $e^{2\pi ix_{i}/L}\psi(\bm{x})$, which satisfies Eq. (8). In contrast, $\hat{U}_{ij}$ defined in Eq. (7) violates the periodicity. The wavefunction of the state

[TABLE]

is not invariant under $\bm{x}\rightarrow\bm{x}+L\bm{e}_{i}$. Thus, $\hat{U}_{ij}|\psi\rangle$ does not belong to the Hilbert space specified by the PBC. Therefore the quantity $\langle\psi|\hat{U}_{ij}|\psi\rangle$ (the inner-product of two states $|\psi\rangle$ and $\hat{U}_{ij}|\psi\rangle$) lacks the physical meaning. Reference Wheeler et al., 2018 proposed to fix this issue by allowing for discontinuities, but then the analyticity of the wavefunction would be lost.

To discuss more general states, let $\hat{T}_{L\bm{e}_{i}}$ be the translation operator that shifts $\bm{x}$ to $\bm{x}+L\bm{e}_{i}$ ($i=x,y$). Imposing the PBC is equivalent to identifying $\hat{T}_{L\bm{e}_{i}}$ with an identity operator. $\hat{U}_{i}$ in Eq. (5) is consistent with this identification because it commutes with $\hat{T}_{L\bm{e}_{i}}$. On the other hand, $\hat{U}_{ij}$ in Eq. (7) does not commute with $\hat{T}_{L\bm{e}_{i}}$ and consequently it violates the boundary condition.

This simple discussion already poses a serious question about the physical meaning of $\tilde{q}_{ij}$ in Eq. (6). Below we discuss the immediate consequence of the lack of periodicity.

II.2 Tight-binding example of quadrupole insulator

Here we consider a simplified version of the four-band tight-binding model introduced in Ref. Benalcazar et al., 2019. See Appendix C for more physical tight-binding models. As illustrated in Fig. 1 (a), the model has four orbitals ($a$, $b$, $c$, and $d$) at each lattice site. The Hamiltonian in Fourier space reads

[TABLE]

We assume $C_{4}$ rotation symmetry to see if $\tilde{q}_{ij}$ is quantized as proposed by Ref. Kang et al., 2018.

[TABLE]

We write the position of lattice sites as $\bm{R}=(x_{n},y_{m})$, where $x_{n}=x_{1}+n-1$ and $y_{m}=x_{1}+m-1$ ($n,m=1,2,\ldots,L$). The lattice constant is set to be unity for simplicity. Here, $x_{1}$ is introduced to investigate the origin dependence of the results. 222See also Appendix of Ref. Kang et al., 2018 where the dependence of $x_{1}$ is discussed using a tight-binding model. Below we consider two familiar choices of $x_{1}$: $x_{1}=1$ ($x_{L}=L$) for every $L$ and $x_{1}=-(L-1)/2$ ($x_{L}=-x_{1}$) for an odd $L$.

The four eigenvalues of $H_{\bm{k}}$ are [math] (doubly degenerate) and $\pm 2$. The Bloch function of the lowest band ($\varepsilon_{\bm{k}}=-2$) reads $u_{\bm{k}}=(1,e^{-ik_{x}},e^{-i(k_{x}+k_{y})},e^{-ik_{y}})^{T}/2$. In the Wannier basis, the ground state that completely occupies the lowest band can be written $|\Phi_{0}\rangle=\prod_{\bm{R}}\hat{\gamma}_{\bm{R}}^{\dagger}|0\rangle$ with

[TABLE]

Using this real space expression, one can analytically evaluate $\tilde{q}_{ij}$. We present the detailed calculation in Appendix B.

Next we consider ionic contributions. To cancel the electric charge and polarization, we place an ion at every plaquette center, i.e. $\bm{x}=\bm{R}+(\bm{e}_{x}+\bm{e}_{y})/2$ as illustrated in Fig. 1 (b). $\tilde{q}_{ij}$ in Eq. (6) for ions reduces to

[TABLE]

We subtract this reference value from the electronic contribution to impose the charge neutrality and vanishing polarization.

[TABLE]

We summarize our results in Table 1. We immediately realize that the value of $\Delta\tilde{q}_{ij}$ depends sensitively on the detailed choice of origin [compare $x_{1}=1$ and $x_{1}=-(L-1)/2$ for an odd $L$] and also on the parity of $L$ [compare the even and odd $L$ for $x_{1}=1$]. Also, the value of $\Delta\tilde{q}_{xy}$ is not necessarily quantized [see the case of $x_{1}=1$ with an even $L$].

III Difficulties in the improvement

III.1 Absence of Resta’s type formula

Let us ask if one can fix the issues in $\tilde{q}_{ij}$ identified above. Here we explore the Resta-type expression $\tilde{q}^{\prime}\equiv\frac{1}{2\pi}\text{Im}\ln\langle\Phi_{0}|\hat{U}^{\prime}|\Phi_{0}\rangle$ with $\hat{U}^{\prime}\equiv e^{2\pi i\sum_{\vec{x}}\theta_{\vec{x}}\hat{n}_{\vec{x}}}$. To be consistent with the PBC, $\theta_{\vec{x}}$ must have the following form

[TABLE]

where $\vartheta_{\bm{x}}\in\mathbb{R}$ is a periodic function of $\bm{x}$. In particular, quadratic terms such as $x_{i}x_{j}/L^{2}$ cannot appear in $\theta_{\bm{x}}$ contrary to Eq. (7).

One may introduce a cut-off function $C_{\vec{x}}\in[0,1]$ that is $1$ (constant) when $\bm{x}$ is far away from the boundary (i.e., $x_{1}\ll x\ll x_{L}$ and $y_{1}\ll y\ll y_{L}$) and smoothly approaches to [math] when $\vec{x}$ is near the boundary with an intermediate length scale $\ell$ ($1\ll\ell\ll L$). For example, we can use

[TABLE]

Then $\vartheta_{\bm{x}}=C_{\bm{x}}\,xy/L^{2}$ becomes effectively a periodic function when terms proportional to $e^{-L/4\ell}$ are neglected. However, we found that such a modification does not resolve the issues, because the contribution to $\tilde{q}^{\prime}$ from the decaying region ($|\partial_{x}C_{\bm{x}}|\sim 1/\ell$) is non-neglegible and spoils the bulk contribution from the interior ($C_{\bm{x}}\simeq 1$).

III.2 Absence of a Berry-phase type formula

Next, let us examine Berry-phase type formulas. According to the modern theory of electric polarization,Vanderbilt and King-Smith (1993); King-Smith and Vanderbilt (1993) the Berry phase

[TABLE]

gives the electric polarization. Here, $|u_{\bm{k}n}\rangle$ is the Bloch function of $n$-th occupied band. This Berry phase can be interpreted as the expectation value of the position operator measured from the origin of the unit cell $\int d^{2}x\sum_{n=1}^{N_{\text{occ}}}w_{n\bm{R}}(\bm{x})^{*}(x_{i}-R_{i})w_{n\bm{R}}(\bm{x})$ of the Wannier state $w_{n\bm{R}}(\bm{x})\equiv L^{-2}\sum_{\bm{k}}e^{i\bm{k}\cdot(\bm{x}-\bm{R})}u_{\bm{k}n}(\bm{x})$. Therefore one may guess that the bulk quadrupole is given by $\int d^{2}x\sum_{n=1}^{N_{\text{occ}}}w_{\bm{R}}(\bm{x})^{*}(x_{i}-R_{i})(x_{j}-R_{j})w_{\bm{R}}(\bm{x})$, i.e.,

[TABLE]

However, this cannot be the case because the quantity in Eq. (19) is not invariant under the gauge transformation $w_{\vec{k}}\in\text{U}(N_{\text{occ}})$ among occupied bands,

[TABLE]

We can fix the gauge-invariance by inserting the projector onto unoccupied bands $1-P_{\bm{k}}$ Marzari and Vanderbilt (1997); Resta and Sorella (1999); Souza et al. (2000)

[TABLE]

However, we then found that this quantity does not produce a useful topological invariant, because $T_{ij}$ is always proportional to an identity matrix in the presence of $n=3$, $4$, or $6$-fold rotation symmetry. [This is because it satisfies $\sum_{i,j}(p)_{i^{\prime}i}T_{ij}(p^{T})_{jj^{\prime}}=T_{i^{\prime}j^{\prime}}$ under a point group symmetry $p\in\text{O}(2)$.] A useful topological invariant would be constructed from an integral with an integer ambiguity, but, to our knowledge, the Berry phase in Eq. (18) is the only combination with that property in two dimensions.

IV Summary

We have analyzed the definition of the bulk quadrupole moment independently proposed by two groups. Kang et al. (2018); Wheeler et al. (2018) We find that the proposed definition of the bulk quadrupole moment fails even for a simple non-interacting example. Our analysis reveals that the issues with the proposed definition are related to violation of periodicity. Possible strategies to fix these issues all seem to fall short.

The obstacles in obtaining a generalization of Resta’s formulation of bulk polarization to higher multipoles are perhaps best illustrated by our findings that even a simpler task, namely, finding a formulation for single-particle systems seems to fail. Although we cannot provide a general proof that such single-particle formulation of bulk quadrupole moment does not exist, we give a strong indication that this task may be a difficult one.

A proper definition of higher multiples in crystals and formulas that allow practical calculations are topics of broad interest, not limited to computational and theoretical solid state physics. Despite the fact that our findings support in some sense a ‘no-go’ statement, we hope that this work will serve as the first step toward the future resolution to defining bulk multiple moments.

Acknowledgements.

The authors would like to thank A. Furusaki, M. Oshikawa, and A. Shitade for useful discussions. We would also like to acknowledge the communications with authors of Refs. Kang et al., 2018; Wheeler et al., 2018 prior to the submission of this manuscript. The work of S.O. is supported by the Materials Education program for the future leaders in Research, Industry, and Technology (MERIT). The work of H. W. is supported by JSPS KAKENHI Grant No. JP17K17678 and by JST PRESTO Grant No. JPMJPR18LA.

Appendix A Subtlety in the perturbative expansion

Resta Resta (1998) used the “first-order perturbation theory” to derive Eq. (4) for interacting systems. Ref. Wheeler et al., 2018 used a similar perturbative argument to verify the formula (7). Here, we review an issue with such perturbative treatment following Ref. Watanabe and Oshikawa, 2018. Only in this appendix do we assume a general $d$ dimensional system with the period $L_{i}$ in $i$-th direction. According to Resta, Resta (1998)

[TABLE]

where $\hat{J}_{x}$ is the sum of the $x$-component of the current operators over the entire space. If the above relation were a controlled expansion in the series of $L_{x}^{-1}$, the expectation value $\langle\Phi_{0}|\hat{U}_{x}|\Phi_{0}\rangle$ would behave as

[TABLE]

which converges to $1$ in the large $L_{x}$ limit. However, this turns out not to be the case because of the volume sum hidden in $\hat{J}_{x}$. Watanabe and Oshikawa (2018) We put the dot over the equality in Eqs. (22) and  (23) as a caution to the reader.

This issue can be readily seen in the case of band insulators. The matrix element of $\hat{U}_{x}$ among Bloch states is given by

[TABLE]

where $n=1,2,\cdots,N_{\text{occ}}$ is the band index of occupied bands and $B_{\bm{k}}$ is an $N_{\text{occ}}$-dimensional matrix representing the discretized Berry connection defined by

[TABLE]

By expanding $B_{\bm{k}}$ to the second order in $L_{x}^{-1}$, we find

[TABLE]

where $V=L_{x}L_{y}\cdots$ is the volume and $G^{ij}$ is the quantum metric tensor defined by Marzari and Vanderbilt (1997); Resta and Sorella (1999); Souza et al. (2000)

[TABLE]

Note that $G_{ij}=O(1)$, i.e., it does not depend on the system size. For example, for the isotropic case in which all $L_{i}$’s are identical to $L$, the right-hand side of Eq. (26) is $e^{-2\pi^{2}G^{xx}L^{{d-2}}+O(L^{{d-3}})}$. In the large $L$ limit, it converges to a number in the range [math] and $1$ in two dimension, and it vanishes in higher dimensions. This is in sharp contrast to the behavior in Eq. (23). Thus the first-order perturbation (22) does not hold in general. 444Even when $d=1$, Eq. (23) is still violated because the exponent of Eq. (26) is $O(L_{x}^{-1})$, not $O(L_{x}^{-2})$.

Quite remarkably, despite the lack of a general proof of Eq. (4) for interacting systems in multi-dimensions, to best of our knowledge there is no counterexample to Eq. (4). Unfortunately, as we discuss in the main text, the circumstances are not so favorable for the validity of the proposed expression (7), where a similar perturbative expansion argument was used. Wheeler et al. (2018)

Appendix B Details on the tight-biding calculation

Here we present the derivation of the result for the model in Eq. (10) summarized in Table 1. To compute $\tilde{q}_{ij}$ in Eq. (6), one has to compute the matrix element of $\hat{U}_{ij}$ in Eq. (7) among occupied single-particle states. To this end, working in the Wannier basis instead of the Bloch basis is advantageous, because all the off-diagonal elements [i.e. $(\hat{U}_{ij})_{\bm{R},\bm{R}^{\prime}}$, $\bm{R}\neq\bm{R}^{\prime}$] vanish for the Wannier state in Eq. (12). Because of this nice property, we have

[TABLE]

The expression for the diagonal matrix element $(\hat{U}_{ij})_{\bm{R},\bm{R}}$ depends on the position of plaquettes, i.e., whether $\bm{R}$ is in the bulk or around the boundary. Introducing the shorthand notation $F_{z}=e^{2\pi iz/L^{2}}$, the diagonal elements can be written as

[TABLE]

and

[TABLE]

To reproduce Table 1, one just has to plug these expressions into Eq. (29) and extract the value in the large $L$ limit.

Appendix C Examples of tight-biding models

In this section, we discuss two tight-binding models that are more physically natural than the simplified one discussed in Sec. II.2. This exercise will support our claim, clarifying the issues in the proposed formula (6).

C.1 Model 1

The first model is illustrated in Fig. 2 (a). As far as the electronic degrees of freedom are concerned, this model is identical to the one in Eq. (10). The only difference is the position of ions. Since lattice sites are located at $\bm{R}=(x_{n},y_{m})$ (see the main text), here we assume that one ion sits at every lattice site $\bm{R}$. This choice gives the reference value

[TABLE]

The lowest band of the model (10) has a Wannier orbital centering at $\bm{x}=\bm{R}+(\frac{1}{2},\frac{1}{2})$ for each $\bm{R}$. The band insulator that completely occupies this band has

[TABLE]

where $(\hat{U}_{xy}^{\text{(a)}})_{\bm{R},\bm{R}}$ is given in Eq. (32). We tabulate the value of $\Delta\tilde{q}_{xy}^{\text{(a)}}=\tilde{q}_{xy}^{\text{(a)}}-\tilde{q}_{xy,\text{ions}}$ in Table 2.

We compare this to the atomic limit illustrated in Fig. 2 (b). In this limit, all electrons are strictly localized at the Wannier center $\bm{x}=\bm{R}+(\frac{1}{2},\frac{1}{2})$ for every $\bm{R}$. The above band insulator can be adiabatically connected to this limit without closing the bulk gap or breaking the $C_{4}$ symmetry. For the atomic limit, we find

[TABLE]

We tabulate the value of $\Delta\tilde{q}_{xy}^{\text{(b)}}=\tilde{q}_{xy}^{\text{(b)}}-\tilde{q}_{xy,\text{ions}}$ in Table 2. Clearly $\Delta\tilde{q}_{xy}^{\text{(b)}}$ does not agree with the atomic limit $\Delta\tilde{q}_{xy}^{\text{(a)}}$, despite the fact that they are smoothly connected to each other. This implies that $\Delta\tilde{q}_{xy}$ cannot, in general, serve as the topological invariant.

To pin down the origin of the mismatch, we compare $(\hat{U}_{xy}^{\text{(a)}})_{\bm{R},\bm{R}}$ and $(\hat{U}_{xy}^{\text{(b)}})_{\bm{R},\bm{R}}$ by plotting

[TABLE]

as a function of $\bm{R}$. Figure 3 shows that the discrepancy originates purely from the boundary. In fact, the sum of $\Delta_{\bm{R}}$ over the boundary (i.e. $R_{x}=x_{L}$ or $R_{y}=y_{L}$) precisely accounts for the difference $\Delta\tilde{q}_{xy}^{\text{(a)}}-\Delta\tilde{q}_{xy}^{\text{(b)}}$. Furthermore, $(\hat{U}_{xy})_{\bm{R},\bm{R}}^{\text{(b)}}$ for the plaquettes near the boundary has a small amplitude and $|\langle\Phi_{0}|\hat{U}_{xy}|\Phi_{0}\rangle|$ decays exponentially with the system size ($e^{-cL}$ with $c\simeq 1.5$ in this particular model). These issues are the manifestation of the violation of the periodicity discussed in Sec. II.1.

C.2 Model 2

The second model is the orthogonal stacking of Su-Schrieffer-Heeger chains (Fig. 4). The tight-binding Hamiltonian reads

[TABLE]

This time we have two occupied bands and we place two ions at every lattice site:

[TABLE]

We analyze this model in the same way as in Appendix B and find $\Delta\tilde{q}_{xy}^{\text{(a)}}{{}^{\prime}}=\tilde{q}_{xy}^{\text{(a)}}{{}^{\prime}}-\tilde{q}_{xy,\text{ions}}{{}^{\prime}}$. For the atomic limit of this insulator where an electron is localized at both $\bm{x}=\bm{R}+(\frac{1}{2},0)$ and $\bm{R}+(0,\frac{1}{2})$ for every unit cell, we find

[TABLE]

We list $\Delta\tilde{q}_{xy}^{\text{(b)}}{{}^{\prime}}=\tilde{q}_{xy}^{\text{(b)}}{{}^{\prime}}-\tilde{q}_{xy,\text{ions}}{{}^{\prime}}$ in Table 3. Although we do not expect any quadrupole moment in this model, we found $\Delta\tilde{q}_{xy}=1/2$ mod 1 for the case (b) when $x_{1}=1$ and $L$ is even.

C.3 Stacking of the two models

Band insulators considered in Appendices C.1 and C.2 possesse a nonzero bulk polarization. We can fix it by stacking them together. We show the values of $\Delta\tilde{q}_{xy}$ in the large $L$ limit in Table 4. As we can see, problems of $\tilde{q}_{xy}$ still persist in the absence of the bulk polarization.

Let us summarize what we learned about $\tilde{q}_{xy}$ through these examples: (i) $\Delta\tilde{q}_{xy}=\tilde{q}_{xy}-\tilde{q}_{xy,\text{ions}}$ depends on $x_{1}$ and the parity of $L$; (ii) it is not always quantized even under the $C_{4}$ symmetry; and (iii) it takes different values for two states (a) and (b) in the same phase (i.e., adiabatically connected).