Geometric orbital magnetization in adiabatic processes

TL;DR

This paper introduces a geometric contribution to orbital magnetization in adiabatic processes of gapped electrons, linking it to the many-body Berry phase and topological invariants.

Contribution

It identifies and formulates the geometric orbital magnetization as a derivative of the many-body Berry phase, including topological and non-topological components.

Findings

The geometric orbital magnetization arises from adiabatic time evolution.

It decomposes into topological and non-topological parts in 2D insulators.

The topological part is expressed as a third Chern-Simons form.

Abstract

We consider periodic adiabatic processes of gapped many-body spinless electrons. We find an additional contribution to the orbital magnetization due to the adiabatic time evolution, dubbed \textit{geometric} orbital magnetization, which can be expressed as derivative of the many-body Berry phase with respect to an external magnetic field. For two-dimensional band insulators, we show that the geometric orbital magnetization generally consists of two pieces, the topological piece that is expressed as third Chern-Simons form in $(t,k_x,k_y)$ space, and the non-topological piece that depends on Bloch states and energies of both occupied and unoccupied bands.

Figures

Equations

$${\hat{j}}_{t\mathbf{\varphi }}^i\equiv \frac{1}{L_i}\int d^{2}x{\hat{j}}_{t\mathbf{\varphi }}^i(\mathbf{x})={\partial }_{{\varphi }_i}{\hat{H}}_{t\mathbf{\varphi }},i=x,y.$$

$$\int \frac{d^2\varphi }{(2\pi {)}^2}⟨{\hat{\mathbf{j}}}_{t\mathbf{\varphi }}⟩={\partial }_t\left(\int \frac{d^2\varphi }{(2\pi {)}^2}⟨{\Phi }_{t\mathbf{\varphi }}|i{\partial }_{\mathbf{\varphi }}|{\Phi }_{t\mathbf{\varphi }}⟩\right),$$

$$\mathbf{j}(t,\mathbf{x})$$

$${\int }_0^{T}dt⟨{\hat{\mathbf{j}}}_{t\mathbf{\varphi }}⟩=-{\partial }_{\mathbf{\varphi }}{\phi }_{\mathbf{\varphi }},$$

$${\phi }_{\mathbf{\varphi }}\equiv {\int }_0^{T}dt⟨{\Phi }_{t\mathbf{\varphi }}|i{\partial }_t|{\Phi }_{t\mathbf{\varphi }}⟩,$$

$$\mathbf{m}={\mathbf{m}}_{\text{pers}}+\text{mathbcal}m.$$

$$TV{m}_z$$

$$\text{mathbcal}m=\text{mathbcal}{m}^{\text{top}}+\text{mathbcal}{m}^{\text{non-top}},$$

$$u_{t\mathbf{k}n}(\mathbf{x}+\mathbf{R})=u_{t\mathbf{k}n}(\mathbf{x}),$$

$$u_{t\mathbf{k}+\mathbf{G}n}(\mathbf{x})=e^{-i\mathbf{G}\cdot \mathbf{x}}{u}_{t\mathbf{k}n}(\mathbf{x}).$$

$$\mathbf{p}(t)=\frac{ei}{V}\sum _{\mathbf{k}n\in \text{occ}}⟨u_{t\mathbf{k}n}|{\nabla }_{\mathbf{k}}{u}_{t\mathbf{k}n}⟩\text{ mod }\frac{e}{a}.$$

$${\mathbf{m}}_{\text{pers}}(t)=\frac{ei}{2V}\sum _{\mathbf{k}n\in \text{occ}}⟨{\nabla }_{\mathbf{k}}{u}_{t\mathbf{k}n}|\times (h_{t\mathbf{k}}+{\epsilon }_{t\mathbf{k}n})|{\nabla }_{\mathbf{k}}{u}_{t\mathbf{k}n}⟩,$$

$$|w_{tn\mathbf{R}}⟩\equiv \frac{a}{\sqrt{V}}\sum _{\mathbf{k}}{e}^{-i\mathbf{k}\cdot \mathbf{R}}|{\psi }_{tn\mathbf{k}}⟩.$$

$$|u_{t\mathbf{k}n}{⟩}^{\prime }=e^{i\mathbf{k}\cdot \mathbf{\delta }}|u_{t\mathbf{k}n}⟩,$$

$${\mathbf{p}}^{\prime }(t)=\mathbf{p}(t)-\frac{e\mathbf{\delta }}{a^2},$$

$${\mathbf{m}}_{\text{pers}}^{\prime }(t)={\mathbf{m}}_{\text{pers}}(t).$$

$$|{\Phi }_{t2\pi {\varphi }_y}⟩=e^{-2\pi i{\hat{P}}_x}|{\Phi }_{t0{\varphi }_y}⟩,$$

$$|{\Phi }_{t{\varphi }_{x}2\pi }⟩=e^{-2\pi i{\hat{P}}_y}|{\Phi }_{t{\varphi }_{x}0}⟩,$$

$$j^{\mu }(t,x)=-\sum _{\nu }{\epsilon }^{\mu \nu }{\partial }_{\theta }{P}_1(\theta ){\partial }_{\nu }\theta ,$$

$$P_1(\theta )\equiv -e\int \frac{dk}{2\pi }\text{tr}{A}_{\theta k},$$

$$(A_{\theta k}{)}_{n,m}\equiv -i⟨u_{\theta kn}|{\partial }_k{u}_{\theta km}⟩.$$

$$\mathbf{Q}\equiv {\int }_0^{T}dt{\partial }_t\mathbf{p}(t)=\frac{1}{T}[\mathbf{p}(T)-\mathbf{p}(0)].$$

$$j^{\mu }(t,\mathbf{x})=-\frac{1}{2\pi }\sum _{\nu ,\lambda ,\rho }{\epsilon }^{\mu \nu \lambda \rho }{\partial }_{\theta }{P}_3(\theta ){\partial }_{\nu }\theta {\partial }_{\lambda }{A}_{\rho }^{\text{ex}},$$

$$P_3(\theta )\equiv -e\int \frac{d^{3}k}{8{\pi }^2}\text{tr}{\mathbf{A}}_{\theta \mathbf{k}}\cdot ({\nabla }_{\mathbf{k}}+\frac{2i}{3}{\mathbf{A}}_{\theta \mathbf{k}})\times {\mathbf{A}}_{\theta \mathbf{k}},$$

$$({\mathbf{A}}_{\theta \mathbf{k}}{)}_{n,m}\equiv -i⟨u_{\theta \mathbf{k}n}|{\nabla }_{\mathbf{k}}{u}_{\theta \mathbf{k}m}⟩.$$

$$\hat{X}={\partial }_{ϵ}{\hat{H}}_{ϵ}{|}_{ϵ=0}$$

$$X\equiv {\int }_0^T\frac{dt}{T}\mathrm{tr}[{\hat{\rho }}_t{\hat{X}}_t].$$

$$⟨{\Phi }_t^M|{\hat{\rho }}_t|{\Phi }_t⟩=⟨{\Phi }_t|{\hat{\rho }}_t|{\Phi }_t^M{⟩}^{\ast }=\frac{i⟨{\Phi }_t^M|{\partial }_t|{\Phi }_t⟩}{E_t^M-E_t}.$$

$$(E_t-E_t^M)⟨{\Phi }_t^M|{\partial }_{ϵ}|{\Phi }_{tϵ}⟩{|}_{ϵ=0}=⟨{\Phi }_t^M|{\partial }_{ϵ}{\hat{H}}_{tϵ}{|}_{ϵ=0}|{\Phi }_t⟩,$$

$$X$$

$${\mathcal{F}}_{tϵ}$$

$$|{\Phi }_{Tϵ}⟩=|{\Phi }_{0ϵ}⟩,$$

$$X$$

$$X_{\text{inst}}$$

$$X_{\text{geom}}$$

Full text

I Introduction

We consider periodic adiabatic processes of spinless short-range entangled phases with period $T$ at zero temperature. The ultimate goal when attacking this type of time-dependent problems on the general ground would be to obtain an expression for the physical observables induced by the time evolution in terms of instantaneous eigenstates and eigenenergies of the system.

In their pioneering works, Niu and Thouless Thouless (1983); Niu and Thouless (1984) found such an expression for the current operator uniformly averaged over the entire space. In the formulation, they assumed the periodic boundary conditions with the period $L_{i}$ for $i=x,y$ and introduced the solenoidal flux $\bm{\phi}=(\phi_{x},\phi_{y})$ as illustrated in Fig. 1. For concreteness we work in two spatial dimensions throughout this work. Then the current operator can be expressed as

[TABLE]

(For brevity we show the dependence on time $t$, flux $\bm{\phi}$, and etc., in the subscript.) Further taking an average over all values of $\bm{\phi}$, the expectation value of the current operator induced by the adiabatic time-evolution can be expressed as the time derivative of the many-body Berry phase for varying $\bm{\phi}$

[TABLE]

where $|\Phi_{t\bm{\phi}}\rangle$ is the instantaneous ground state of the Hamiltonian $\hat{H}_{t\bm{\phi}}$. This expression assumes the periodicity in $\bm{\phi}$ (see Eqs. (17), (18) below). It is not possible to further impose the periodicity in time simultaneously. Instead we have $|\Phi_{T\bm{\phi}}\rangle=e^{-i\bm{\phi}\cdot\bm{Q}}|\Phi_{0\bm{\phi}}\rangle$, where $\bm{Q}=\int_{0}^{T}dt\int\frac{d^{2}\phi}{(2\pi)^{2}}\langle\hat{\bm{j}}_{t\bm{\phi}}\rangle\in\mathbb{Z}^{2}$ is the pumped charge during in one cycle.

The result (2) is formally similar to the constitutive relation for Maxwell’s equations

[TABLE]

where $\bm{p}$ and $\bm{m}$ are the bulk polarization and the bulk magnetization. Later it was shown King-Smith and Vanderbilt (1993); Vanderbilt and King-Smith (1993); Resta (1994); Resta and Vanderbilt (2007) that Thouless result (2) combined with constitutive relation (3) gives a useful formula for the bulk polarization — this development marked the birth of “the modern theory” of electric polarization. The bulk polarization is given by the integral of the Berry connection (the first Chern-Simons form) $P_{1}$ [see Eq. (11) for the formula for band insulators].

Alternatively, let us impose the periodicity in time $|\Phi_{T\bm{\phi}}\rangle=|\Phi_{0\bm{\phi}}\rangle$ instead. In this setting it is useful to integrate over time rather than the solenoidal flux. Then the Thouless result reads Niu and Thouless (1984)

[TABLE]

where $\varphi_{\bm{\phi}}$ is the many-body Berry phase associated with the adiabatic time-evolution. In the thermodynamic limit, $\partial_{\bm{\phi}}\varphi_{\bm{\phi}}$ is independent of $\bm{\phi}$ Niu and Thouless (1984); Watanabe (2018) and one can set $\bm{\phi}=\bm{0}$ for instance. There is also a contribution from $\partial_{\bm{\phi}}E_{t\bm{\phi}}$ in Eqs. (2), (4) but it is negligibly small for the same reason. We find this formulation of the Thouless pump more useful because it can be generalized to wider class of physical observables as we discuss below.

The persistent current associated with a part of the orbital magnetization can also be expressed using the instantaneous eigenstates and eigenenergies of the Hamiltonian. For band insulators, it can be written as the curl of a vector (see Fig. 2a), which together with the constitutive relation (3), allows one to define the orbital magnetization $\bm{m}_{\text{pers}}$. (The subscript refers to the contribution associated with the persistent current.) Alternatively, one can evaluate the change of the instantaneous ground state energy with respect to the external magnetic field. This recent development Xiao et al. (2005); Thonhauser et al. (2005); Ceresoli et al. (2006); Shi et al. (2007) goes under the name of “the modern theory” of the orbital magnetization [see Eq. (12)]. Unlike the bulk polarization, $\bm{m}_{\text{pers}}$ is not related to topological response.

In this work we develop a general formulation of the remaining contribution to the electric current in the constitutive relation  (3) that are neither captured by the averaged current in Eqs. (2), (4) nor by the persistent current $\bm{\nabla}\times\bm{m}_{\text{pers}}(t,\bm{x})$. We find that, after coarse-graining in time, this contribution can be expressed as the curl of an additional term $\mathbcal{m}$ to the orbital magnetization so that $\bm{m}$ in Eq. (3) is given by

[TABLE]

Our main result is that $\mathbcal{m}$ can be obtained as a derivative of the many-body Berry phase with respect to an external magnetic field $B_{z}$ applied in $z$ direction

[TABLE]

where $V$ represents the system size and $\varphi_{B_{z}}$ is defined by Eq. (5) upon substitution $\bm{\phi}\rightarrow B_{z}$. This expression is well-defined in two-dimensional systems with the open boundary condition at least in one direction. There are known subtleties when applying uniform magnetic field to periodic systems. See Sec. III.1 for the detailed discussion. In the following we assume vanishing Chern numbers in $(\phi_{x},\phi_{y})$, $(t,\phi_{x})$ and $(t,\phi_{y})$ spaces.

As comparison, in the presence of an external uniform magnetic field $\bm{B}=(0,0,B_{z})^{\rm T}$, the instantaneous orbital magnetization $\bm{m}_{\text{pers}}$ gives an energy shift $\int_{0}^{T}dt(E_{tB_{z}}-E_{t0})/T=-V\bm{m}_{\text{pers}}\cdot\bm{B}+O(\bm{B}^{2})$ of the many-body ground state. (Throughout this work we set $\hbar=1$.) Accordingly, after the period $T$, the ground state acquires an additional phase proportional $TV\bm{m}_{\text{pers}}\cdot\bm{B}$. On the other hand, a non-zero value of $\mathbcal{m}$ shows up as Berry phase $TV\mathbcal{m}\cdot\bm{B}+O(\bm{B}^{2})$ acquired by the many-body ground state. For this reason, we name $\mathbcal{m}$ geometric orbital magnetization. The bulk quantity $T\mathbcal{m}$ is independent of the period $T$ and is defined “mod $e$”. This ambiguity is reflecting the possibility of decorating the boundary by one-dimensional Thouless pump.

For band insulators, we perform the perturbation theory with respect to the applied magnetic field following Refs. Shi et al., 2007; Essin et al., 2010 and find that geometric orbital magnetization consists of two contributions

[TABLE]

the topological contribution $\mathbcal{m}^{\text{top}}$ is expressed as integral of the third Chern-Simons form $P_{3}$ in $(t,k_{x},k_{y})$ space [see Eq. (51)], while the non-topological contribution $\mathbcal{m}^{\text{non-top}}$ is written in terms of instantaneous Bloch states and energies [Eq. (52)]. The obtained expression for $\mathbcal{m}$ of band insulators has a formal similarity with the expression for the magnetoelectric polarizability of three-dimensional band insulators Essin et al. (2010); Malashevich et al. (2010); Chen and Lee (2011) upon identification $t/T\leftrightarrow k_{z}/2\pi$. It is worth mentioning that due to relatively large gap (order of electronvolts) of band insulators, the adiabaticity conditions is not particularly restrictive, the period $T$ can be as small as several femtoseconds.

To gain intuitive understanding of the two contributions (8), one can think of the topological piece $\mathbcal{m}^{\text{top}}$ to be originating from the Aharanov-Bohm contribution to the many-body Berry phase in the magnetic field. Thus $\mathbcal{m}^{\text{top}}$ describes the magnetization from electrons, whose positions are moving during the adiabatic process, as depicted with dashed lines and arrows in Fig. 2b-c and e. Although Fig. 2a may look similar, $\bm{j}_{\text{pers}}$ in Fig. 2a represents a static persistent current that is uniformly distributed on the ring. In contrast, the current density in Fig. 2b at each time is localized to the position of the potential well and it becomes divergence-free only after averaging over the period $T$. Similarly, the non-topological piece $\mathbcal{m}^{\text{non-top}}$ can be understood to be originating from “spinning” of anisotropic crystalline potentials, see Fig 2d.

Let us mention at this point several related works. Adiabatic dynamics can be induced by time-dependent lattice deformations (phonons), which is the subject of studies on dynamical deformations of crystals. Ceresoli and Tosatti (2002); Juraschek et al. (2017); Juraschek and Spaldin (2018); Dong and Niu (2018); Stengel and Vanderbilt (2018) In Refs. Juraschek et al., 2017; Juraschek and Spaldin, 2018 it was shown that a time-varying polarization gives rise to a contribution to the orbital magnetization, and the semi-classical description developed in Ref. Dong and Niu, 2018 found the same effect within their framework. The time-varying polarizations in these works correspond to the situation depicted in Fig. 2b. In the case of band insulators, we find that they are captured by the abelian third Chern-Simons form. Furthermore, Refs. Ceresoli and Tosatti, 2002; Stengel and Vanderbilt, 2018 showed that rotation of molecules, as in Fig. 2d, gives rise to an orbital magnetization contribution that can be captured by relation (7). The present approach gives unified description of the above-mentioned effects. More importantly, it properly describes the orbital magnetization in adiabatic processes that have not been previously considered: the process in Fig. 2e has inversion symmetry at all times, thus polarization is time-independent, yet it gives rise to non-zero $\mathbcal{m}$, which for the case of band insulators is captured by non-abelian third Chern-Simons form, see Sec. IV.3.

Analogous to the bulk polarization, crystalline symmetries can quantize topological geometric orbital magnetization $\mathbcal{m}^{\text{top}}$. We show that, under certain crystalline symmetries, $\mathbcal{m}^{\text{top}}$ is related to recently discussed higher-order topological phases. Parameswaran and Wan (2017); Schindler et al. (2018); Peng et al. (2017); Langbehn et al. (2017); Song et al. (2017); 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) Among them, the topological insulators that exhibit quantized corner charges in the presence of certain crystalline symmetries attracted recently a lot of theoretical Benalcazar et al. (2017, 2017, 2019) and experimental Serra-Garcia et al. (2018); Peterson et al. (2018) attention. Although, due to crystalline symmetries, the bulk quadrupole moment is well defined in these systems (see Fig. 4a), it is still disputed in the literature whether such definition is possible in the absence of any quantizing crystalline symmetries. Kang et al. (2018); Lapa and Hughes (2018); Ono et al. (2019) Recent work in Ref. van Miert and Ortix, 2018a revealed a connection between higher-order topological insulators Parameswaran and Wan (2017); Schindler et al. (2018); Peng et al. (2017); Langbehn et al. (2017); Song et al. (2017); 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) protected by roto-inversion symmetries and adiabatic processes that involve topological insulators with quantized corner charges. In Sec. III.5 we show that the adiabatic processes discussed by van Miert and Ortix van Miert and Ortix (2018a) are characterized by quantized geometric orbital magnetization, and we relate the value of $T\mathcal{m}_{z}^{\text{top}}$ to the quantized corner charge.

The remaining of this article is organized as follows: in Sec. II we review the modern theory of the polarization and the orbital magnetization, Sec. III contains derivation of our main results, Sec. III.5 discusses the role of symmetries in the adiabatic process, and Sec. IV presents various non-interacting examples that illustrate difference between instantaneous orbital magnetization, topological and non-topological geometric orbital magnetization. More precisely, we consider toy models illustrating systems depicted in Fig. 2. As a more realistic application of physics considered in this work, we present in Sec. IV.5 calculation of magnetization induced by rotation of an insulator. A long time ago, Barnett (1915, 1935) Barnett considered magnetization of an uncharged paramagnetic material when spun on its axis. Modeling paramagnetic material as collection of local magnetic moments that are randomly oriented, Barnett Barnett (1915) argued that rotation creates a torque that acts to align local magnetic moments with rotation axis. This torque gives rise to magnetization $M=\chi\Omega/\gamma$, where $\chi$ is paramagnetic susceptibility, $\Omega$ is rotation frequency and $\gamma$ is electron gyromagnetic ratio. Barnett’s measurement of this effect Barnett (1915) provided first accurate measurement of electron gyromagnetic ratio. We calculate $\mathbcal{m}$ for this model, which, as seen from Eq. (7), is also proportional to rotational frequency $\Omega=2\pi/T$ and estimate quantum correction to Barnett effect. Electron contribution to $\mathbcal{m}$ has both topological and non-topological piece, but since the system is uncharged, we find that electron contribution to $\mathbcal{m}^{\text{top}}$ is canceled by corresponding ionic contribution. Thus resulting $\mathbcal{m}$ is solely due to anisotropy of crystalline potential, analogous to toy model in Fig. 2d. In Sec. V we consider examples of general interacting systems where periodic adiabatic process consists of “spinning” Ceresoli and Tosatti (2002); Stengel and Vanderbilt (2018) or “shaking” Juraschek et al. (2017); Juraschek and Spaldin (2018); Dong and Niu (2018) of the whole system. Goldman et al. (2014) Our conclusions and outlook can be found in Sec. VI.

II Preliminaries

Here we review the formulation of the polarization and the orbital magnetization for band insulators in $2+1$ dimensions developed in Refs. King-Smith and Vanderbilt, 1993; Vanderbilt and King-Smith, 1993; Resta, 1994; Resta and Vanderbilt, 2007; Xiao et al., 2005; Thonhauser et al., 2005; Ceresoli et al., 2006; Shi et al., 2007. To simplify notations, we assume primitive lattice vectors of the square lattice type, but this general framework is not restricted to this special choice.

II.1 Modern theory

Let us denote by $\psi_{t\bm{k}n}(\bm{x})=(a/\sqrt{V})e^{i\bm{k}\cdot\bm{x}}u_{t\bm{k}n}(\bm{x})$ the instantaneous Bloch function of $n$-th occupied band, satisfying $h_{t}|\psi_{t\bm{k}n}\rangle=\varepsilon_{t\bm{k}n}|\psi_{t\bm{k}n}\rangle$. Here $h_{t}$ is the single-particle Hamiltonian with a periodic potential, $V=L_{x}L_{y}$ is the system size and $a$ is the lattice constant. We choose the cell-periodic gauge so that they obey the following conditions for any lattice vector $\bm{R}$ and reciprocal lattice vector $\bm{G}$ Vanderbilt (2018)

[TABLE]

According to the modern theory, the bulk polarization density $\bm{p}(t)$ is given by

[TABLE]

where $e$ $(<0)$ is the electric charge. The sum over $\bm{k}$ can be replaced with the integral $V\int\frac{d^{2}k}{(2\pi)^{2}}$ over the first Brillouin zone. Similarly, the orbital magnetization density $\bm{m}_{\text{pers}}(t)$ is given by

[TABLE]

where $h_{t\bm{k}}\equiv e^{-i\bm{k}\cdot\bm{x}}h_{t}e^{i\bm{k}\cdot\bm{x}}$. The ambiguity in Eq. (11) can be seen by a smooth gauge transformation $|u_{t\bm{k}n}\rangle^{\prime}=\sum_{m}|u_{t\bm{k}m}\rangle(U_{\bm{k}})_{m,n}$ that changes the integral in Eq. (11) by an integer multiple of $e/a$, while the integral in (12) remains unchanged.

In addition to the derivation via the Thouless pump as we described in the introduction, the formula (11) was also verified in terms of the Wannier state localized around the unit cell $\bm{R}$.

[TABLE]

In terms of the Wannier function, $\bm{p}(t)$ is the deviation of the Wannier center from $\bm{R}$, i.e., $\bm{p}(t)=\frac{e}{a^{2}}\int d^{2}x(\bm{x}-\bm{R})|w_{tn\bm{R}}(\bm{x})|^{2}$.

When the origin of unit cell is changed by $\bm{\delta}$, we find

[TABLE]

Namely, $\bm{p}(t)$ depends on the specific choice of the origin, while $\bm{m}_{\text{pers}}(t)$ does not. Therefore, it is not $\bm{p}(t)$ itself but rather the change $\Delta\bm{p}(t)$ that is of physical interest. It also follows that for an periodic adiabatic process, where the system is translated by certain number of unit cells during the period $T$, the orbital magnetization is periodic in time while the polarization is not.

For interacting systems under the periodic boundary condition, the combination in the parenthesis in Eq. (2) replaces Eq. (11). The periodicity in $\phi_{i}$ in this formulation is encoded in the relation

[TABLE]

where $\hat{P}_{i}$ is the polarization operator (see Ref. Watanabe and Oshikawa, 2018 for example).

II.2 Topological response

To discuss the physical consequence of $\Delta\bm{p}(t)$, let us recall first the topological linear response in $(1+1)$ dimension Thouless (1983) that holds at a mesoscopic scale after coarse-graining

[TABLE]

Here, $x^{\mu}$ ($\mu=0,1$) represents $(t,x)$, and $j^{\mu}$ corresponds to $(n,j)$. $A_{\theta k}$ is a (finite dimensional) matrix constructed by occpied Bloch states and the trace in Eq. (20) is the matrix trace. Comparing with above equations, we see that $P_{1}(\theta(t))$ is the 1D version of $\bm{p}(t)$ in Eq. (11). This response is derived starting from the Chern-Simons theory $j^{\mu}=\frac{C_{1}}{2\pi}\sum_{\nu\lambda}\varepsilon^{\mu\nu\lambda}\partial_{\nu}A_{\lambda}^{\text{ex}}$ in $(2+1)$ dimensions that describes the response toward an external field $A^{\text{ex}}$ and reducing the dimension to $(1+1)$ dimensions.

The parameter $\theta$ in Eq. (19) is a slowly varying field interpolating between two different systems. For example, an adiabatic time evolution $\theta(t)$ induces the bulk current $j(t,x)=\partial_{t}P_{1}(\theta(t))$. The bulk charge transfer from $t=0$ to $t=T$ is thus given by $\int_{0}^{T}dtj(t,x)=P_{1}(\theta(T))-P_{1}(\theta(0))$. Similarly, a transition of one 1D system to another can be described by $\theta(x)$, giving rise to a charge density $n(t,x)=-\partial_{x}P_{1}(\theta(x))$. Therefore, the total charge $Q^{\text{edge}}$ accumulated to the boundary is $Q^{\text{edge}}=\int_{x_{0}}^{x_{1}}dxn(t,x)=P_{1}(\theta(x_{0}))-P_{1}(\theta(x_{1}))$. For a given $\theta$ that specifies $P_{1}(\theta)$ as a continuous function of $t$ and $x$, even the integer part of $Q^{\text{edge}}$ is well-defined. However, only the fractional part of $Q^{\text{edge}}$ is independent of the detailed choice of the interpolation — the fractional part depends only on the initial and the final values of $P_{1}$ that can be individually computed by Eq. (20). What we described here can be straightforwardly translated to 2D systems. The pumped charge through the bulk per unit length along $\bm{n}$ is given by $\bm{Q}\cdot\bm{n}$, where

[TABLE]

The analog of Eq. (19) in $(3+1)$ dimensions reads Qi et al. (2008)

[TABLE]

Here, $\mu,\nu,\rho,\lambda=0,1,2,3$. Again, $\bm{A}_{\theta\bm{k}}$ is defined by occupied Bloch states. This response is derived from the Chern-Simons theory $j^{\mu}=\frac{C_{2}}{8\pi^{2}}\varepsilon^{\mu\nu\lambda\rho\sigma}\partial_{\nu}A_{\lambda}^{\text{ex}}\partial_{\rho}A_{\sigma}^{\text{ex}}$ in $(4+1)$ dimensions by a dimensional reduction. This topological response implies, for example, the magnetoelectric effect Qi et al. (2008); Essin et al. (2009, 2010) $\rho(z)=-\frac{1}{2\pi}\partial_{z}P_{3}(\theta(z))B_{z}^{\text{ex}}$.

III Geometric orbital magnetization

In this section we present the derivations of our main results. We start with verifying the most general expression  (7) for the geometric orbital magnetization. Then we derive the formula for the topological and non-topological contributions in Eq. (8).

III.1 Berry phases in adiabatic process

Suppose we are interested in the expectation value of the quantity $\hat{X}$, given by

[TABLE]

for some parameter $\epsilon$ in the Hamiltonian. For example, in the case of the averaged current operator, $\epsilon$ can be identified with the solenoidal flux $\bm{\phi}$ [see Eq. (1)]. Likewise, for the orbital magnetization we use the external magnetic field $B_{z}$.

Now suppose that the Hamiltonian $\hat{H}_{t}$ has a periodic adiabatic dependence on $t$, and let $|\Phi_{t}\rangle$ be the instantaneous ground state with the energy eigenvalue $E_{t}$. We assume an excitation gap $\Delta_{t}$ and the time-dependence of the Hamiltonian must be slow enough so that $\Delta_{t}T\gg 1$. Using the density matrix $\hat{\rho}_{t}=|\Psi_{t}\rangle\langle\Psi_{t}|$ obeying the time-dependent Schrödinger equation $\partial_{t}\hat{\rho}_{t}=-i[\hat{H}_{t},\hat{\rho}_{t}]$, we express the time-average of the expectation value of $\hat{X}_{t}$ as

[TABLE]

In the absence of the time-evolution, the density matrix is identical to $|\Phi_{t}\rangle\langle\Phi_{t}|$. It acquires contributions from excited states $\hat{H}_{t}|\Phi_{t}^{M}\rangle=E_{t}^{M}|\Phi_{t}^{M}\rangle$ due to the time evolution. To the lowest-order perturbation theory with respect to $(\Delta_{t}T)^{-1}$, the relevant matrix elements are given by Thouless (1983); Watanabe and Oshikawa (2018)

[TABLE]

Now we plug $\hat{X}_{t}\equiv\partial_{\epsilon}\hat{H}_{t\epsilon}|_{\epsilon=0}$ and make use of the Sternheimer identity

[TABLE]

which follows by differentiating $\hat{H}_{t\epsilon}|\Phi_{t\epsilon}\rangle=E_{t\epsilon}|\Phi_{t\epsilon}\rangle$ with respect to $\epsilon$ at $\epsilon=0$. In our notation, $\hat{H}_{t\epsilon}|_{\epsilon=0}=\hat{H}_{t}$ and $|\Phi_{t\epsilon}\rangle|_{\epsilon=0}=|\Phi_{t}\rangle$. Combining these equations, we find

[TABLE]

where

[TABLE]

is the Berry curvature in $(t,\epsilon)$ space. Further assuming the periodicity in time

[TABLE]

we arrive at our general expression

[TABLE]

where $X_{\text{inst}}$ is the time average of the expectation value using the instantaneous ground state and $X_{\text{geom}}$ is the geometric contribution originating from the adiabatic time dependence. This is the generalization of Eq. (4) for the electric current to physical observables written as the derivative of Hamiltonian as in Eq. (26). The following basis-independent expressions may also be useful

[TABLE]

where $\hat{P}_{t\epsilon}=|\Phi_{t\epsilon}\rangle\langle\Phi_{t\epsilon}|$ is the projector onto the many-body ground state, $\hat{G}_{t}=(z-\hat{H}_{t})^{-1}$ is the many-body Green function, and the integration contour encloses only the ground state at $z=E_{t}$.

Let us now specialize to the case $\epsilon=B_{z}$. Then $X_{\text{geom}}$ in Eq. (35) gives the geometric orbital magnetization $\mathbcal{m}$ in Eq. (7), while $X_{\text{inst}}$ is the persistent current contribution. Note the additional minus sign because of $\hat{M}_{z}=-\partial_{B_{z}}\hat{H}_{B_{z}}$. Previously, Refs. Ceresoli and Tosatti, 2002; Stengel and Vanderbilt, 2018 considered the Berry curvature in $(t,B_{z})$ space to describe orbital magnetization induced by rotation of molecules. See also examples in Sec. IV.4 and V.1 below.

When applying these formulae, one has to be careful about boundary conditions. If open boundary conditions in at least one direction are imposed, the result (36) is directly applicable. However, a process that is periodic in time under periodic boundary condition may loose its periodicity in time under open boundary conditions. For example, the system in Fig. 2c is not periodic in time if the open boundary condition in $y$ direction is imposed. Similarly, the periodicity in time requires periodic boundary conditions in both directions for the $C_{4}$-symmetric system in Fig. 5. Keeping (original) periodic boundary conditions in both directions in the presence of magnetic field, implies that the net flux through the system has to vanish. If the system is homogeneous, local contributions to $\mathbcal{m}$ cancel out and we cannot obtain a useful information about the system. (For single-particle problems there is a resolution as we discuss below.) Finally, one can impose magnetic periodic boundary conditions assuming that the total magnetic flux applied to the system $B_{z}L_{x}L_{y}$ is an integer multiple of $2\pi$. Brown (1964); Zak (1964) However, each eigenstate of $\hat{H}_{B_{z}}$ may not be analytic as a function of $B_{z}$ despite the fact that the magnetic field $B_{z}=2\pi/L_{x}L_{y}$ itself can be made small for a large systems. The expression (36) is still applicable if the projector onto the instantaneous ground state is analytic function of $B_{z}$, which is the case for band insulators with vanishing Chern number. Essin et al. (2010); Malashevich et al. (2010); Gonze and Zwanziger (2011); Chen and Lee (2011) However, to our knowledge there is no general proof for gapped interacting systems.

III.2 Noninteracting systems

Let us apply this general expression to noninteracting electrons described by the quadratic Hamiltonian

[TABLE]

We label single-particle states in such a way that $\varepsilon_{t\epsilon n+1}\geq\varepsilon_{t\epsilon n}$ for all $n=1,2,\cdots$. We also assume a finite gap $\Delta=\varepsilon_{t\epsilon N+1}-\varepsilon_{t\epsilon N}$ between $N$-th and $(N+1)$-th levels. We write the single particle state $|\gamma_{t\epsilon n}\rangle\equiv\hat{\gamma}_{t\epsilon n}^{\dagger}|0\rangle$. Then the $N$-particle ground state can be written as

[TABLE]

For a later purpose, we allow for a unitary transformation among the occupied levels

[TABLE]

Although such a basis change may sound unnecessary, in the actual application of this framework it is sometimes important to work in the proper basis by choosing $U_{n\ell}$ appropriately. (See Sec. III.3 for an example). After all we find that the many-body Berry phase $\varphi_{\epsilon}$ is given by the sum of single-particle Berry phases $\varphi_{\epsilon\ell}$ of occupied levels

[TABLE]

Therefore, we get the following expressions for single-particle problems. The latter two expressions are basis-independent

[TABLE]

where $P_{t\epsilon}=\sum_{\ell=1}^{N}|\psi_{t\epsilon\ell}\rangle\langle\psi_{t\epsilon\ell}|=\sum_{\ell=1}^{N}|\gamma_{t\epsilon n}\rangle\langle\gamma_{t\epsilon n}|$ is the projector onto occupied single-particle states, $h_{t\epsilon}=\sum_{n}\varepsilon_{tn\epsilon}|\gamma_{t\epsilon n}\rangle\langle\gamma_{t\epsilon n}|$ is the single-particle Hamiltonian, $g_{t}=(z-h_{t})^{-1}$ is single-particle Green function, and the integration contour encloses all the occupied states at $z=\varepsilon_{tn}$ ($n=1,2,\cdots,N$).

For the orbital magnetization, we again set $\epsilon=B_{z}$. The same remarks as in the previous section apply here. In the case of band insulators, one may want to impose periodic boundary conditions to preserve the translation symmetry. As discussed in the previous section there are two possibilities to achieve this. One can change the the boundary condition to magnetic periodic boundary conditions. Brown (1964); Zak (1964) For single-particle systems, assuming symmetric gauge $\bm{A}^{\rm ex}(\bm{x})=\bm{B}\times\bm{x}/2$, the magnetic periodic boundary conditions can be taken into account explicitly by restricting the form of the projector $\hat{P}_{tB_{z}}$ to Essin et al. (2010); Gonze and Zwanziger (2011) $\langle\bm{x}_{1}|\hat{P}_{tB_{z}}|\bm{x}_{2}\rangle=\hat{P}^{\prime}_{tB_{z}}(\bm{x}_{2},\bm{x}_{1})e^{ieB_{z}\bm{x}_{1}\times\bm{x}_{2}\cdot\hat{\bm{z}}/2}$, where $P^{\prime}_{tB_{z}}(\bm{x}_{1},\bm{x}_{2})$ is an arbitrary $N\times N$ matrix function (not necessarily projector) that satisfies $P^{\prime}_{tB_{z}}(\bm{x}_{1}+\bm{R},\bm{x}_{2}+\bm{R})=P^{\prime}_{tB_{z}}(\bm{x}_{1},\bm{x}_{2})$, where $\bm{R}$ is an element of Bravais lattice. The expression for $\hat{P}^{\prime}_{tB_{z}}(\bm{x})$ can be found perturbatively in $B_{z}$, Essin et al. (2010); Gonze and Zwanziger (2011) which, after substituting back to Eq. (41), yields an expression for the Berry phase and $\mathbcal{m}$. The second option is to apply a spatially modulating magnetic field as we discuss below.

III.3 Geometric orbital magnetization for band insulators

Below we consider band insulators and show that geometric orbital magnetization has two contributions as in Eq. (8). Since we assume the periodic boundary condition both in $x$ and $y$, we apply a slowly modulating magnetic field Shi et al. (2007); Essin et al. (2010) in order to avoid changing of the boundary condition as discussed in the previous subsection. We use the vector potential

[TABLE]

with $\bm{e}_{z}\equiv(0,0,1)^{\rm T}$, $q\equiv 2\pi/L$, and $f(\bm{x})=(\cos qx+\cos qy)/2$. (To simplify the notation, we assume $L=L_{x}=L_{y}$ in this subsection). Such a magnetic field induces the change of the Bloch function

[TABLE]

and $|\partial_{\epsilon}w_{t\epsilon n\bm{R}}\rangle|_{\epsilon=0}$ is given via Eq. (13).

We compute the Berry phase using the formula (41) derived above. It is important to work in the Wannier basis for which the magnetic field effectively becomes uniform in the limit $q\rightarrow 0$. The single-particle Berry phase in this basis, summed over occupied bands, takes the following form

[TABLE]

If we further sum over $\bm{R}$, or equivalently if we work in the Bloch basis $|\psi_{tn\bm{k}}\rangle$, we get [math] reflecting the fact that for bulk systems Fourier component $\mathcal{m}_{z}(\bm{q})$ vanishes for $\bm{q}\neq 0$. Thus, care must be taken to correctly read off local contribution to $\mathcal{m}_{z}$— an unintentional integration over $\bm{R}$ of a term proportional to $f(\bm{R})$ makes it impossible to find the correct value of $\mathcal{m}_{z}$. The rest calculation follows the appendix in Ref. Essin et al., 2010. Upon taking the limit $q\rightarrow 0$, we find

[TABLE]

This last expression can precisely be expressed as the sum of two terms, $\mathbcal{m}^{\text{top}}+\mathbcal{m}^{\text{non-top}}$. The topological piece $\mathbcal{m}^{\text{top}}$ reads

[TABLE]

The Berry connection $(\bm{A}_{\bm{K}})_{n,m}\equiv-i\langle u_{\bm{K}n}|\bm{\nabla}_{\bm{K}}u_{\bm{K}m}\rangle$ is defined using occupied Bloch states as a function of $\bm{K}\equiv(t,\bm{k})$. The smoothness and the periodicity of $\bm{A}_{\bm{K}}$ are assumed in the integral in Eq. (51). Such a choice is possible only when both the pumped charge through the bulk $\bm{Q}$ in Eq. (22) and the 2D Chern number for $(k_{x},k_{y})$ vanish.

The non-topological contribution depends also on instantaneous eigenenergies of the Bloch Hamiltonian

[TABLE]

Here, $h_{t\bm{k}}$ is Bloch Hamiltonian, $P_{t\bm{k}}=\sum_{n\in\text{occ}}|u_{tn\bm{k}}\rangle\langle u_{tn\bm{k}}|$ is the projector onto occupied bands at $\bm{k}$, $g_{t\bm{k}}=(z-h_{t\bm{k}})^{-1}$ is Bloch’s Green function, the curly brackets denote symmetrization $\{\bm{A}\times\bm{B}\}=\bm{A}\times\bm{B}+\bm{B}\times\bm{A}$, and the integration contour encloses all the filled Bloch states at $z=\varepsilon_{tn\bm{k}}$. See the appendix of Ref. Essin et al., 2010 for the details. Note that both $\mathbcal{m}^{\text{top}}$ and $\mathbcal{m}^{\text{non-top}}$ are not affected by the shift of the origin in Eq. (14).

III.4 Topological contribution from response theory

Here we give an alternative, easier derivation of $\mathbcal{m}^{\text{top}}$ in Eq. (50) from the topological response theory. To this end let us further reduce one spatial dimension in Eq. (23) to achieve the topological quadratic response in $(2+1)$d Qi et al. (2008):

[TABLE]

where $F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}+i[A_{\mu},A_{\nu}]$ is the Berry curvature in the $(k_{x},k_{y},\theta,\phi)$ space and $\theta$ and $\phi$ are two slowly varying fields: $\theta(t)$ denotes an adiabatic and periodic time dependence and $\phi(\bm{x})$ describes a smooth interface of domains (Fig. 3). In this setting, we find

[TABLE]

so that

[TABLE]

This reproduces Eq. (50). In the derivation we used the relation $\int_{0}^{2\pi}\frac{d\theta}{2\pi}G_{2}(\theta,\phi)=\partial_{\phi}P_{3}(\phi)$. It is important to note that $\mathbcal{j}(t,\bm{x})$ itself cannot be written as a curl of a vector field — Eq. (III.4) holds only after the time convolution (or equivalently the time average).

The above derivation relies on the connection (3) between $\mathbcal{m}^{\text{top}}$ and topological edge current in adiabatically driven two-dimensional systems. To see this more concretely, let us consider the boundary of two regions with $\phi_{0}\equiv\phi(\bm{x}_{0})$ and $\phi_{1}\equiv\phi(\bm{x}_{1})$ (see Fig. 3). Just like in the case of polarization, only the fractional part of the edge current is the bulk contribution that depends only on $\phi_{0}$ and $\phi_{1}$. This can be understood by noticing that decorating the boundary with a 1D chain leads to an integer charge transfer through the Thouless pump. Thouless (1983) To capture the fractional bulk contribution to the edge current, one can separately compute $\mathcal{m}_{z}(\bm{x}_{0})$ and $\mathcal{m}_{z}(\bm{x}_{1})$ without paying attention to their continuity. The geometric contribution to the charge transfer along $i$ direction, $i=x,y$ between two bulk systems with $\mathcal{m}_{z}(\bm{x}_{0})$ and $\mathcal{m}_{z}(\bm{x}_{1})=\mathcal{m}_{z}(\bm{x}_{1}^{\prime})$ (Fig. 3) is given by

[TABLE]

Notice that $I^{\text{edge}}$ of two adjacent edges may differ by an integer. To see this formally, let us consider a charge flow $\Delta Q^{\text{corner}}$ into a corner surrounded by a closed curve $\bm{x}_{\alpha}$ with $\bm{x}_{1}=\bm{x}_{0}$ (see Fig. 3). The net charge flow in the process is given by the second Chern number

[TABLE]

For example, when the corner is formed by two edges along $x$ and $y$ directions, we have

[TABLE]

meaning that the charge transfer along two intersecting edges can only differ by an integer multiple of $e$. Clearly, $\Delta Q^{\text{corner}}$ is not a bulk topological invariant in general, since its value can be changed by closing the boundary gap, i.e., attaching 1D Thouless pump at certain boundaries (Figs. 3 and 4b).

III.5 Symmetry constraints and corner charge

Here we consider adiabatic process of two-dimensional systems constrained by certain symmetries that quantize $T\mathcal{m}_{z}^{\text{top}}$. We show that, if the symmetry allows one to define the bulk contribution to quadrupole moment, the quantized quadrupole moment is equal to $T\mathcal{m}_{z}^{\text{top}}$. Such adiabatic processes were recently discussed by van Miert and Ortix, who found the connection between the quantized corner charge and higher-order topological invariant. van Miert and Ortix (2018b)

For concreteness, let us consider the four-fold rotation $C_{4}$ mapping $\bm{x}=(x,y,0)$ to $C_{4}\bm{x}=(-y,x,0)$. It is easy to see that boundary decorations by polarized one-dimensional chains do not affect the fractional part of the corner charge $\Delta Q^{\text{corner}}$, see Fig. 4a. We consider an arbitrary interpolation between the system of interest $h_{0\bm{k}}$ and the reference system $h_{T/2\,\bm{k}}$ that has no corner charge. The second half of the cycle is performed in a $C_{4}$-symmetric manner

[TABLE]

The $C_{4}$ symmetry defined above behaves as the roto-inversion $IC_{4}$ in $(t,k_{x},k_{y})$-space, resulting in the following transformation law for $\mathcal{m}_{z}^{\text{top}}$:

[TABLE]

This does not mean that $T\mathcal{m}_{z}^{\text{top}}$ vanishes since it is defined only ${\rm mod}\,1$. Thus in the presence of $C_{4}$ symmetry $T\mathcal{m}_{z}$ is quantized either [math] or $e/2\mod e$. When $T\mathcal{m}_{z}=e/2\mod e$, the circulating edge current as in Fig. 3 violates $C_{4}$ symmetry constraint (60)—the only allowed edge current distribution is shown by black arrows in Fig. 4b. Note that the inversion symmetry, for example, also quantizes $T\mathcal{m}_{z}^{\text{top}}$ but the total corner charge accumulation during inversion-symmetric cycles need to vanish since the charge distribution of quadrupole moment is invariant under the inversion (see Fig. 4c).

Now we show that the parity of the corner charge accumulation $\Delta Q^{\text{corner}}$ is actually a bulk topological invariant for symmetric adiabatic processes satisfying constraint (60), see also Fig. 4b. To this end, consider two perpendicular edges along $x$ and $y$ direction, related to each other by $C_{4}$ symmetry. The relations (59) and (61) suggest that $I^{\text{edge}}_{y}=-I^{\text{edge}}_{x}$ and that

[TABLE]

Furthermore, Fig. 4c tells us that the corner charge accumulation during the symmetric process is $\Delta Q^{\text{corner}}=2q^{\text{corner}}$. Therefore,

[TABLE]

We will discuss an example of quadrupole insulators with $P_{3}=e/2$ in Sec. IV.3 using this result. On the other hand, a $C_{4}$-symmetric phase that hosts a corner charge of $q^{\text{ corner}}=e/4$ were recently reported. Benalcazar et al. (2019) The fact that $\Delta Q^{\text{ corner}}\in\mathbb{Z}$ forces us to conclude that $C_{4}$-symmetric adiabatic process cannot be constructed for such a phase.

Alternatively, as discussed in detail in Ref. van Miert and Ortix, 2018b, the $C_{4}$-symmetric adiabatic process $h_{t\bm{k}}$ considered above, can be viewed as a 3D topological insulator protected by the roto-inversion symmetry $IC_{4}$ upon identification $k_{z}=2\pi t/T$. In fact, the 3D topological insulator with $P_{3}=e/2$ obtained this way is a second-order topological insulator. If we consider a geometry with the open boundary conditions in $xy$-plane and the periodic boundary conditions in $z$-direction, such a second-order phase can be translationally invariant in $z$-direction both in the bulk and on the boundary. The boundary hosts an odd number of chiral modes running along each of four hinges in $IC_{4}$-symmetric manner. Going back to the picture of an adiabatic process, it becomes clear that the corner charge accumulation is an odd integer as $t$ is varied from [math] to $T$, which is consistent with the above result (62).

IV Examples: noninteracting systems

In this section, we discuss a simple model of noninteracting spinless electrons in a periodic potential, which highlights the distinction of two contributions to the bulk orbital magnetization, $\bm{m}_{\text{pers}}$ and $\mathbcal{m}$. Additionally, we want to consider examples where there is only topological geometric magnetization, Sec. IV.3, only non-topological geometric magnetization, Sec. IV.4, and both topological and non-topological contributions, Sec. IV.5. To keep the discussions simple while capturing the relevant physics, we focus on isolated orbitals without any overlap between them.

IV.1 Bloch functions in the localized limit

Let us consider a time-dependent deep potential $v_{t}^{0}(\bm{x})$ centering at $\bm{x}=\bm{r}(t)$ that accommodates at least one bound state. Let $\phi_{t}^{0}(\bm{x})$ be the wavefunction of the instantaneous lowest-energy bound state, satisfying $h_{t}^{0}\phi_{t}^{0}(\bm{x})=\varepsilon_{t}^{0}\phi_{t}^{0}(\bm{x})$ with

[TABLE]

Here $\bm{A}_{t}^{\text{ex}}(\bm{x})$ describes an external field. In these expressions, the superscript [math] implies the quantities for an isolated orbit. When the potential $v_{t}^{0}(\bm{x})$ is deep enough, $\phi_{t}^{0}(\bm{x})$ should be well-localized around $\bm{x}=\bm{r}(t)$ with the localization length $\xi\ll a$. Hence, we assume that

[TABLE]

and that both $\phi_{t}^{0}(\bm{x})$ and $v_{t}^{0}(\bm{x})$ decays fast enough, i.e., $|v_{t}^{0}(\bm{x})|,|\phi_{t}^{0}(\bm{x})|\rightarrow 0$ as $|\bm{x}-\bm{r}(t)|\gg\xi$.

With these building blocks, we construct a periodic potential and the cell-periodic Bloch state.

[TABLE]

We assume that $\bm{A}^{\text{ex}}(\bm{x})$ respects the periodicity, i.e., $\bm{A}^{\text{ex}}(\bm{x}-\bm{R})=\bm{A}^{\text{ex}}(\bm{x})$. Then, as far as $\phi_{t}^{0}(\bm{x}-\bm{R})^{*}\phi_{t}^{0}(\bm{x})$ and $v_{t}^{0}(\bm{x}-\bm{R})\phi_{t}^{0}(\bm{x})$ ($\bm{R}\neq\bm{0}$) are entirely neglected, $u_{t\bm{k}}(\bm{x})$ is an eigenstate of the periodic Hamiltonian

[TABLE]

with a completely flat band dispersion $\varepsilon_{t\bm{k}}=\varepsilon_{t}^{0}$.

IV.2 Polarization and instantaneous magnetization

Let us first demonstrate the modern theory formula for the polarization and the orbital magnetization by deriving $\bm{p}$ and $\bm{m}$ in two different ways.

First, we present direct calculation of the polarization and the orbital magnetization from the microscopic charge distribution and the persistent current densities in this insulator. The instantaneous contribution to the local charge and current distribution from a single orbit $\phi_{t}^{0}(\bm{x})$ can be written as

[TABLE]

We introduce vector fields $\bm{p}_{t}^{0}(\bm{x})$ and $\bm{m}_{t}^{0}(\bm{x})$ such that

[TABLE]

The existence of such $\bm{m}_{t}^{0}(\bm{x})$ is guaranteed by the divergence-free nature of the instantaneous current density $\bm{j}_{t}^{0}(\bm{x})$. The current density induced by the adiabatic motion of $\bm{r}(t)$ is captured by $\mathbcal{j}_{t}^{0}(\bm{x})$ in Eq. (81) whose divergence may not vanish. We assume both $\bm{p}_{t}^{0}(\bm{x})$ and $\bm{m}_{t}^{0}(\bm{x})$ decay rapidly for $|\bm{x}-\bm{r}(t)|>\xi$, which specifies the boundary condition for differential equations (71).

Physical quantities of the insulator composed of periodically arranged localized orbits can be written as the sum of the contributions from each orbit. For example microscopic current is given by

[TABLE]

and analogously for $n$, $\bm{p}$, and $\bm{m}$. These microscopic expressions have a strong spatial dependence, periodically oscillating at the scale of $a$. To derive to a smooth mesoscopic description, we need to perform a convolution in space (Sec. 6.6 of Ref. Jackson, 1999). Here we choose the Gaussian $g(\bm{x})=(\pi R^{2})^{-1}e^{-|\bm{x}|^{2}/R^{2}}$ ($R\gg a$)

[TABLE]

We do the same for other quantities. Relations such as $\bm{j}(t,\bm{x})=\bm{\nabla}\times\bm{m}(t,\bm{x})$ are preserved by the convolution. Because the convolution is identical to the average for the periodic distribution, we find $n(t,\bm{x})=\bar{n}=\frac{e}{a^{2}}$, $\bm{j}(t,\bm{x})=\bm{0}$,

[TABLE]

This is the part of the orbital magnetization produced by the persistent current as illustrated in Fig. 2a.

Let us check that we get the same results using the general formulae of the modern theory. Because of the non-overlapping assumption of $\phi_{t}^{0}(\bm{x})$, it can be readily shown that the formula in Eqs. (11), (12) for the Bloch function (67) can be simplified to

[TABLE]

where we used Eqs. (69) and (70). These are well-known expressions in classical electrodynamics for the charge and current distributions in a confined region (see Sec. 4.1 and 5.6 of Ref. Jackson, 1999). The equivalence of Eqs. (74), (75) and (76), (77) can be easily checked by using the definition of $\bm{p}_{t}^{0}$ and $\bm{m}_{t}^{0}$ in Eq. (71) and integrating by parts. The second equality of Eq. (76) follows from Eqs. (65), (69), and (76).

IV.3 Topological geometric magnetization

Next, we discuss the topological geometric contribution $\mathbcal{m}^{\text{top}}$ for this model. To this end, suppose that the position of the potential minimum $\bm{r}(t)$ adiabatically moves as a function of $t\in[0,T]$ and forms a closed curve as illustrated in Fig. 2b. We assume the form of the potential, and thus the localization length, remains unchanged during the adiabatic process.

We first apply our general expression for $\mathbcal{m}^{\text{top}}$ in Eq. (50) to the Bloch function (67). Thanks to the non-overlapping assumption, the vector potential $\bm{A}_{(t,\bm{k})}$ is $\bm{k}$-independent:

[TABLE]

with $A_{t}\equiv-i\int d^{2}x\phi_{t}^{0}(\bm{x})^{*}\partial_{t}\phi_{t}^{0}(\bm{x})$. Plugging this into Eq. (51), we find

[TABLE]

where $S_{\bm{r}}$ represents the area enclosed by the orbit of $\bm{r}(t)$ in one cycle. Therefore,

[TABLE]

Observe the analogy to $\bm{m}$ in Eq. (77). This expression does not have integer ambiguity because it is given by abelian third Chern-Simons form.

Let us verify this result from a direct calculation. The adiabatic motion of the single orbit following the potential minimum at $\bm{x}=\bm{r}(t)$ induces a local current distribution

[TABLE]

It becomes divergence-free if averaged over one period

[TABLE]

As we have seen above, the sum of such microscopic currents from each unit cell produces the bulk magnetization

[TABLE]

This agrees with Eq. (80) because the integral in the parenthesis is precisely $\bm{r}(t)$ due to Eq. (65).

The result in Eq. (80) can be readily generalized to the case with multi-orbitals, such as examples in Fig. 2c and e. Let us introduce potential minima $\bm{x}=\bm{r}_{n}(t)$ ($n=1,2,\cdots$) in each unit cell, which are adiabatically varied as a function of $t\in[0,T]$. This time, each orbit is allowed to form an open curve, as far as the total polarization $\bm{p}(t)=(e/a^{2})\sum_{n}\bm{r}_{n}(t)$ satisfies $\bm{p}(T)=\bm{p}(0)$. Under such an assumption, we find that

[TABLE]

We present the proof in the Appendix A. Although the above expression appears to be the sum of single-band contributions, the “would-be” contribution from each band depends on the specific choice of the origin when it does not form a closed loop. Only after performing the summation over all occupied bands, or in other words, only after fully taking into account the non-abelian nature of the third Chern-Simons form, the result restores the independence from the origin choice.

As the application of the formula (85), let us discuss the corner charge of the $C_{4}$-symmetric quadrupole insulator introduced in Refs. Benalcazar et al., 2017, 2017. For the wallpaper group $p4$, there exist three spacial Wyckoff positions: the unit cell origin at $\bm{x}_{a}=(0,0)$, the center of the plaquette at $\bm{x}_{b}=(a/2,a/2)$, and the center of bonds at $\bm{x}_{c}=(a/2,0)$, $(0,a/2)$. Hahn (2006) In the nontrivial phase, the two occupied Wannier orbitals locate at $\bm{x}_{b}$, while in the trivial phase they are at $\bm{x}_{a}$. We consider a periodic adiabatic process illustrated in Fig. 5 starting with the nontrivial phase at $t=0$ and passing the trivial phase at $t=T/2$. The instantaneous Hamiltonian $h_{t\bm{k}}$ itself breaks the $C_{4}$-symmetry down to $C_{2}$ symmetry except when $t=0$ and $T/2$, while the adiabatic process as a whole implements the full $C_{4}$ in the sense of Eq. (60). We can readily compute $P_{3}$ of this process using Eq. (85) which turns out to be $e/2$. This is the corner charge of the quadrupole insulator as predicted by Eq. (63), which agrees with the original study. Benalcazar et al. (2017, 2017) A variant of this adiabatic process was also discussed in Ref. van Miert and Ortix, 2018b.

IV.4 Non-topological geometric magnetization

In this example we first consider a single electron in an anisotropic and rotating two-dimensional well, Goldman et al. (2014) see Fig. 2d. We assume a harmonic confining potential, i.e., Hamiltonian (64) with

[TABLE]

where $\bm{x}_{t}\equiv(x\cos\Omega t+y\sin\Omega t,-x\sin\Omega t+y\cos\Omega t)$ and $\Omega=2\pi/T$. To obtain the geometric orbital magnetization for this model, we consider external magnetic field $\bm{B}=B_{z}\bm{e}_{z}$ described by the vector potential $\bm{A}^{\rm ex}(\bm{x})=\bm{B}\times\bm{x}/2$. (Strictly speaking, this form is valid only around the origin as it lacks the required periodicity.) The wave function of the instantaneous ground state of this model can be obtained based on Ref. Rebane, 2012:

[TABLE]

where $\cal N$ is the normalization factor and $\omega_{c}\equiv eB_{z}/m$ is the cyclotron frequency. The Berry phase $\varphi_{B_{z}}$ during the adiabatic process $t\in[0,T]$ is

[TABLE]

From Eq. (7) it follows that the adiabatic process (87) has non-zero geometric orbital magnetic moment $\mathcal{m}_{z}^{0}$

[TABLE]

Now we construct the Bloch function (67) using $\phi_{t}^{0}(\bm{x})$ as the building block and compute the geometric orbital magnetization $\mathcal{m}_{z}$ based on Eq. (52) for the corresponding band insulator. To this end we need instantaneous eigenstates and eigenenergies in absence of external magnetic field including unoccupied bands. The Hamiltonian $h_{t\bm{k}}$ is given by Eq. (68) with $v_{t}^{0}(\bm{x})$ given by Eq. (87) and $\bm{A}^{\rm ex}=0$. We assume that there is no overlap between wavefunctions belonging to different unit cells as before. (When $\omega_{x},\omega_{y}$ are large enough, such an assumption is valid at least for relevant low-energy states.) Bloch wavefunctions read

[TABLE]

where $\bm{n}\equiv(n_{x},n_{y})$ labels energy levels of two-dimensional the anisotropic harmonic oscillator and $\bm{n}=(0,0)$ corresponds to the ground state in Eq. (88) with $\omega_{c}=0$. Substituting above expressions to Eqs. (50) and (52), we find

[TABLE]

IV.5 Geometric magnetization by rotation

Here we calculate the contribution to the geometric orbital magnetization of a rotating uncharged body and compare it to the classical Barnett effect. Barnett (1915, 1935) The Barnett effect predicts magnetization $\chi/\gamma\Omega$, where $\chi$ is the paramagnetic susceptibility, $\gamma$ is the electron gyromagnetic ratio, and $\Omega$ is the rotation frequency. Since the rotation axis does not necessarily coincide with potential well minima we have $v_{t}^{0}(\bm{x}-\bm{r}(t))$ with $v_{t}^{0}$ from Eq. (87) and $\bm{r}(t)=(R\cos\Omega t,R\sin\Omega t,0)^{\rm T}$, where $R$ is the distance of the potential well minima to the rotation axis. The lowest-energy instantaneous wavefunction $\phi_{t}(\bm{x})$ can be obtained from Eq. (88) by performing gauge transformation

[TABLE]

As compared to Eq. (89), the Berry phase $\varphi_{B_{z}}$ during the adiabatic process $t\in[0,T]$ acquires an additional contribution $eB_{z}\pi R^{2}$ from the Aharonov-Bohm phase. Therefore electrons contribute to the following geometric orbital magnetization

[TABLE]

The first term can be identified with $\mathcal{m}_{z}^{\text{top}}$ in Eq. (80) and the second term is the contribution in Eq. (90). Since the body is uncharged, the contribution from ions cancels the topological contribution, while $\mathcal{m}_{z}^{\text{non-top}}=\mathcal{m}_{z}^{0}$ remains since ions are much more localized compared to electrons. Assuming anisotropy $\omega_{x}/\omega_{y}=2$, and confinement of electrons on the scale of angstroms, we find that contribution (93) is on the same order as Barnett effect for paramagnets with paramagnetic susceptibility $\chi\sim 10^{-5}$. For comparison, paramagnets have typically magnetic susceptibility $\chi\sim 10^{-3}-10^{-5}$. Wikipedia contributors (2019)

V Examples: finite interacting systems

In this section we demonstrate the validity of Eq. (7) for finite interacting systems. We consider two canonical ways of introducing the time-dependence to the Hamiltonian: rotating Ceresoli and Tosatti (2002); Stengel and Vanderbilt (2018) and translating the whole system. Goldman et al. (2014); Juraschek et al. (2017); Dong and Niu (2018)

We consider many-body systems under the open boundary condition in two spatial dimensions. We start with a time-independent Hamiltonian $\hat{H}$ that can contain arbitrary interactions. The total charge, current, polarization, and orbital magnetization operator for this Hamiltonian can be written as $\hat{\bm{N}}=\int_{V}d^{2}x\hat{\bm{n}}(\bm{x})$, $\hat{\bm{J}}=\int d^{2}x\hat{\bm{j}}(\bm{x})$, $\hat{\bm{X}}=\int_{V}d^{2}x\bm{x}\hat{\bm{n}}(\bm{x})$, $\hat{\bm{M}}=(1/2)\int_{V}d^{2}x\bm{x}\times\hat{\bm{j}}(\bm{x})$. We stress that these expressions are valid only when the system is confined in a finite region; they need to be modified in extended systems under the periodic boundary conditions as done by the modern theory. We denote the many-body ground state of $\hat{H}$ and its energy by $|\Phi\rangle$ and $E$, respectively.

To compute the many-body Berry phase, let $\hat{H}_{\bm{B}}$ be the Hamiltonian with the vector potential in the symmetric gauge $\bm{A}^{\text{ex}}(\bm{x})=(1/2)\bm{B}\times\bm{x}$ with $\bm{B}=B_{z}\bm{e}_{z}$. Expanding to the linear order in $B_{z}$ and using $\hat{\bm{j}}(\bm{x})=-\partial_{\bm{A}(\bm{x})}\hat{H}$, we get

[TABLE]

Therefore, the ground state of $\hat{H}_{B_{z}}$ to the leading order in $B_{z}$ can be expressed as

[TABLE]

Here $\hat{Q}\equiv 1-|\Phi\rangle\langle\Phi|$ is the projector onto excited states.

V.1 Rotation

Here we consider the time-dependence of the problem induced by the rotation of the whole system

[TABLE]

where $\Omega=\bm{e}_{z}\Omega$ is the rotation frequency and $\hat{\bm{L}}$ is the angular momentum operator. For the time-dependent Hamiltonian $\hat{H}_{t}$, the orbital magnetization operator $\hat{\bm{M}}_{t}\equiv(1/2)\int_{V}d^{2}x\bm{x}\times\hat{\bm{j}}_{t}(\bm{x})$

[TABLE]

We evaluate the instantaneous contribution $\bm{m}_{\text{pers}}$ and the geometric contribution $\mathbcal{m}$ to the orbital magnetization via the formulae in Eqs. (34) and (35). The instantaneous contribution is given by the instantaneous ground state $|\Phi_{t}\rangle=e^{-i\hat{L}_{z}\Omega t}|\Phi\rangle$

[TABLE]

The geometric contribution is given by the many-body Berry phase. Since the instantaneous ground state of the Hamiltonian $\hat{H}_{tB_{z}}\equiv e^{-i\hat{L}_{z}\Omega t}\hat{H}_{B_{z}}e^{i\hat{L}_{z}\Omega t}$ is given by $|\Phi_{tB_{z}}\rangle=e^{-i\hat{L}_{z}\Omega t}|\Phi_{B_{z}}\rangle$, we have

[TABLE]

This is the expectation value of $\hat{L}_{z}\Omega$ in the presence of the perturbation $-\hat{m}_{z}B_{z}$ in Eq. (96). Using Eq. (97), we get

[TABLE]

We verify these results by solving time-dependent problem. The solution to the time-dependent Schrödinger equation $i\partial_{t}|\Psi_{t}\rangle=\hat{H}_{t}|\Psi_{t}\rangle$ can be readily constructed using the ground state $|\Phi_{\Omega}\rangle$ of the time-independent Hamiltonian

[TABLE]

The solution that is smoothly connected to the ground state in the static limit $\Omega\rightarrow 0$ reads

[TABLE]

The time-average of the orbital magnetization is thus given by

[TABLE]

This is the expectation value of $\hat{\bm{M}}$ in the presence of the perturbation $-\hat{L}_{z}\Omega$ as in Eq. (103). The first-order perturbation theory with respect to $\Omega$ gives

[TABLE]

This is precisely $\bm{m}_{\text{pers}}+\mathbcal{m}$ predicted above in Eqs. (100) and  (102). As it is clear from the derivation, the agreement of the two independent approaches is guaranteed by the Maxwell relation for the free energy $\hat{F}\equiv\hat{H}-\hat{L}_{z}\Omega-\hat{M}_{z}B_{z}$

[TABLE]

V.2 Translation

Next let us introduce the time-dependence by the translation. All discussions proceed in essentially the same way, while there are still a few differences. First we define the time-dependent Hamiltonian by

[TABLE]

where $\hat{T}_{t}=e^{-i\hat{\bm{P}}\cdot\bm{r}(t)}$ is the translation by amount $\bm{r}(t)$ and $\hat{\bm{P}}$ is the momentum operator. For $\hat{H}_{t}$ the orbital magnetization operator becomes

[TABLE]

where the second term in the parenthesis is due to the change of the origin. The instantaneous ground state $|\Phi_{t}^{\prime}\rangle=\hat{T}_{t}|\Phi\rangle$ gives $\bm{m}_{\text{pers}}$ as in Eq. (100), where we used $\langle\Phi|\hat{\bm{J}}|\Phi\rangle=\bm{0}$.

Next, we compute the geometric contribution via the many-body Berry phase. In the presence of magnetic field translation operator $\hat{T}_{tB_{z}}\equiv\hat{T}_{B_{z}}(\bm{r}(t))$ becomes translation followed by gauge transformation Brown (1964); Zak (1964)

[TABLE]

The instantaneous ground state of $\hat{H}_{tB_{z}}\equiv\hat{T}_{tB_{z}}\hat{H}_{B_{z}}\hat{T}_{tB_{z}}^{\dagger}$ is $|\Phi_{tB_{z}}\rangle=\hat{T}_{tB_{z}}|\Phi_{B_{z}}\rangle$, thus the many-body Berry phase reads

[TABLE]

Here, $\bm{S}_{\bm{r}}\equiv\frac{1}{2}\oint\hat{\bm{r}}\times d\bm{r}$ represents the area swept by $\bm{r}(t)$ in one cycle. In the derivation, we used

[TABLE]

Therefore, when the whole system is translated, the geometric contribution to the orbital magnetization captures the Aharonov-Bohm phase

[TABLE]

To verify the above results, we consider the time-dependent Schrödinger equation $i\partial_{t}|\Psi_{t}^{\prime}\rangle=\hat{H}_{t}^{\prime}|\Psi_{t}^{\prime}\rangle$. An approximate solution is given by $|\Psi_{t}^{\prime}\rangle=\hat{T}_{t}|\Phi_{t\bm{r}}\rangle$, where $|\Phi_{t\bm{r}}\rangle$ is the instantaneous ground state of the Hamiltonian $\hat{H}_{t\bm{r}}\equiv\hat{H}-\hat{\bm{P}}\cdot\partial_{t}\bm{r}(t)$. Therefore, the time-average of the orbital magnetization is

[TABLE]

In the adiabatic limit, this reproduces $\bm{m}_{\text{pers}}+\mathbcal{m}$ in Eqs. (100) and (113). In the derivation we used $\langle\Phi_{t\bm{r}}|\hat{\bm{J}}|\Phi_{t\bm{r}}\rangle=eN\partial_{t}\bm{r}(t)$ for the ground state of $\hat{H}_{t\bm{r}}$.

VI Conclusion

In order to obtain current and charge distribution in a medium, one needs to solve Maxwell’s equation together with two constitutive relations [see Eq. (3)] that fully characterize the medium at the mesoscopic scale. The modern theories, developed in the last 30 years, provide handy formulae to calculate electric polarization King-Smith and Vanderbilt (1993); Vanderbilt and King-Smith (1993); Resta (1994); Resta and Vanderbilt (2007) and orbital magnetization Xiao et al. (2005); Thonhauser et al. (2005); Ceresoli et al. (2006); Shi et al. (2007) for realistic materials.

The focus of this work is on spinless short-range entangled systems under periodic adiabatic evolution. Our main result is to identify an additional contribution to the orbital magnetization that we name geometric orbital magnetization $\mathbcal{m}$. This new contribution is defined only after performing the time-average over the period of the adiabatic process, which makes the current density divergence-free. We find that the geometric orbital magnetization can be expressed compactly as derivative of the many-body Berry phase with respect to an externally applied magnetic field. For band insulators, we obtain handy formulae for the bulk geometric orbital magnetization $\mathbcal{m}$ in terms of instantaneous Bloch states and energies. Interestingly, we find that for band insulators $\mathbcal{m}=\mathbcal{m}^{\text{top}}+\mathbcal{m}^{\text{non-top}}$ consists of two pieces, where topological piece $\mathbcal{m}^{\text{top}}$ depends only on the Bloch states of occupied bands. For spinless systems only electric polarization and orbital magnetization enter constitutive relations, since the contributions from higher moments are typically negligible. Jackson (1999) In this sense, our results together with “the modern theories” provide a complete mesoscopic description of a medium under periodic adiabatic time evolution. In this work we have not considered adiabatic processes with ground state degeneracy, Meidan et al. (2011) it would be interesting to see to which extent our findings can be generalized to such systems.

In the present work, the adiabaticity assumption is crucial for validity of the obtained results. In practice, for band insulators with band gaps on the order of electronvolt, this conditions requires that the period $T$ is larger than couple of femtoseconds. Nevertheless, shorter period $T$ results in a larger geometric orbital magnetization. It would be therefore interesting to extend our results to the case of strong drive that excites unoccupied bands. In the case of Thouless pumps such extension was very fruitful and resulted in recent discovery of shift currents. von Baltz and Kraut (1981); Morimoto and Nagaosa (2016)

Although higher (than dipole) electric and magnetic multiple moments typically do not enter constitutive relations, the knowledge of these quantities may be useful for certain systems. Benalcazar et al. (2017, 2019) In fact, it is a topic of current research whether higher moments can be established as bulk quantities in general. Gao and Xiao (2018); Shitade et al. (2018); Kang et al. (2018); Lapa and Hughes (2018); Ono et al. (2019) In the presence of certain crystalline symmetries, both electric polarization and topological geometric orbital magnetization can be quantized, in which case they can serve as a topological invariants. In this context, we showed that the quantized quadrupole moment is related to $\mathcal{m}^{\text{top}}_{z}$ in systems with proper symmetries that allow bulk definition of the quadrupole moment. Ono et al. (2019)

In this work we succeeded in separating $\mathbcal{m}$ into the topological and the non-topological piece only for band insulators. There, we found that the topological contribution is expressed as the third Chern-Simons form ($P_{3}$), in $(t,k_{x},k_{y})$ space. For interacting systems, based on examples considered in Sec. V, we conclude that it is possible to separate Aharonov-Bohm contribution originating from the center of mass motion. In fact, this contribution can be captured by calculating $P_{3}$ formally defined for the many-body ground state as a function of time and two solenoidal fluxes. Clearly, the many-body $P_{3}$ defined in such manner is abelian and does not capture all possible topological contributions. For example, it vanishes for the model in Fig. 4. As a future direction, it would be interesting to see if separation achieved for band insulators is possible for general single-particle, or even many-body systems. The affirmative answer to this question would provide a way to define $P_{3}$ in two dimensional systems with adiabatic time-dependence lacking the translational invariance or the single-particle description. The formula for $P_{3}$ in many-body three-dimensional systems already exist in the literature, Wang and Zhang (2014); Shiozaki et al. (2018) where it was argued that $P_{3}$ is related to the magnetoelectric polarizability. The magnetoelectric polarizability of three-dimensional materials contains, at least for the case of band insulators, not only topological but also non-topological contribution, Essin et al. (2010) thus the analogous “separation question” arises also in that context. Additionally, defining quantized quadrupole moment for interacting systems is one of the open questions. Kang et al. (2018); Lapa and Hughes (2018); Ono et al. (2019) Since for band insulators we find connection between $\mathcal{m}_{z}^{\text{top}}$ and quantized quadrupole moment, separating $\mathcal{m}_{z}^{\text{top}}$ contribution in interacting systems might provide useful many-body definition of quantized quadrupole moment.

We hope that our work will also have practical implication as it contributes to emerging field of “dynamical material design” by providing a way to calculate additional orbital magnetization contribution that appears in these systems. Ceresoli and Tosatti (2002); Juraschek et al. (2017); Juraschek and Spaldin (2018); Dong and Niu (2018); Stengel and Vanderbilt (2018)

Acknowledgements.

We would like to thank David Vanderbilt for drawing Refs. Ceresoli and Tosatti, 2002; Juraschek et al., 2017; Juraschek and Spaldin, 2018; Dong and Niu, 2018; Stengel and Vanderbilt, 2018 to our attention. The work of S.O. is supported by 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 Topological geometric orbital magnetization in simple models

Here we present the derivation of Eq. (85) in Sec. IV.3 for multi potential minima at $\bm{x}=\bm{r}_{n}(t)$ ($n=1,2,\cdots,N_{\rm occ}$). To calculate the geometric orbital magnetization, we use the cell-periodic Bloch state (67) for each orbital. Substitution of $(\bm{A}_{\bm{K}})_{nm}=(A_{tn},-\bm{r}_{n}(t))^{\rm T}\delta_{nm}$ with $\bm{K}\equiv(t,\bm{k})$ [c.f. Eq. (78)] into Eq. (51), we find

[TABLE]

However, this is not the complete expression of $P_{3}$ because Eq. (78) implies that the above Berry connection violates the periodicity in time when $\bm{r}_{n}(T)$ differs from $\bm{r}_{n}(0)$. To quantify the violation, let us introduce the mismatch matrix

[TABLE]

To repair periodicity of $\bm{A}_{\bm{K}}$ we need to find a family of unitary matrices $\Theta_{\bm{K}}$ that interpolates $\Theta_{0\bm{k}}$ and $\Theta_{T\bm{k}}\equiv\mathds{1}_{N_{\rm occ}\times N_{\rm occ}}$. Such an interpolation exists if and only if the polarization $\bm{p}(t)=(e/a^{2})\sum_{n}\bm{r}_{n}(t)$ satisfy $\bm{p}(T)=\bm{p}(0)$. Accordingly, below we focus on the case $\bm{Q}=0$. The gauge transformation with this $\Theta_{\bm{K}}$ gives a periodic Berry connection

[TABLE]

as far as

[TABLE]

The full expression of $P_{3}$ is given by computing Eq. (51) with $\tilde{\bm{A}}_{\bm{K}}$. It is given by Ryu et al. (2007)

[TABLE]

where the sum over $i,j=x,y$ and $\mu,\nu,\lambda=t,x,y$ are assumed. $\Delta_{\text{PI}}$ is the Pontryagin index of the matrix $\Theta_{\bm{K}}$, which requires an explicit expression for $\Theta_{\bm{K}}$, whereas the knowledge of mismatch matrix $\Theta_{0\bm{k}}$ in Eq. (116) suffices for the boundary term $\Delta_{\text{boundary}}$. It reads

[TABLE]

Thus, in order to prove (85), it remains to show that $\Delta_{\text{PI}}=0$. Below we show this by explicitly constructing $\Theta_{\bm{K}}$.

Consider a family of $2\times 2$ unitary matrices $\Xi_{\bm{K}}(\bm{r})$

[TABLE]

with $f(0)=0$ and $f(T)=T$. The above matrix interpolates between unitary matrix $\Xi_{0\bm{k}}(\bm{r})=e^{i\sigma_{3}\bm{k}\cdot\bm{r}}$ and the identity matrix $\Xi_{T\bm{k}}(\bm{r})=\mathds{1}_{2\times 2}$. We also require $f^{\prime}(0)=f^{\prime}(T)=0$ so that the requirement (118) is satisfied. One such function $f(t)$ is given by

[TABLE]

Since $\Xi_{\bm{K}}(\bm{r})$ depends only on $\bm{k}\cdot\bm{r}$, it is easy to check that its Pontryagin index vanishes. We then define $\Theta_{\bm{K}}$ by a composition

[TABLE]

where each interpolation matrix $\Theta^{(n)}_{\bm{K}}$ is defined as,

[TABLE]

with

[TABLE]

It is straightforward to check that $\Theta^{(N_{\rm occ}-1)}_{T\bm{k}}$ equals to the identity matrix and that the Pontryagin index for each $\Theta^{(n)}_{\bm{K}}$ vanishes.