Geometric orbital magnetization in adiabatic processes

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.


I. INTRODUCTION
We consider periodic adiabatic processes of spinless shortrange entangled phases with period T . 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 1,2 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 φ = (φ x , φ y ) as illustrated in Fig. 1. For concreteness we work in two spatial dimensions throughout this work. Then the current operator can be expressed aŝ (For brevity we show the dependence on time t, flux φ, and etc., in the subscript.) Further taking an average over all values of φ, 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 φ where |Φ tφ is the instantaneous ground state of the Hamil-tonianĤ tφ . This expression assumes the periodicity in φ (see Eqs. (17), (18) below). It is not possible to further impose the periodicity in time simultaneously. Instead we have |Φ T φ = e −iφ·Q |Φ 0φ , where Q = T 0 dt d 2 φ (2π) 2 ĵ tφ ∈ Z 2 is the pumped charge during in one cycle.
The result (2) is formally similar to the constitutive relation for Maxwell's equations where p and m are the bulk polarization and the bulk magnetization. Later it was shown 3-6 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 |Φ T φ = |Φ 0φ instead. In this setting it is useful to integrate over time rather than the solenoidal flux. Then the Thouless result reads 2 where ϕ φ is the many-body Berry phase associated with the adiabatic time-evolution. In the thermodynamic limit, ∂ φ ϕ φ is independent of φ 2,7 and one can set φ = 0 for instance. There is also a contribution from ∂ φ E tφ 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 m 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 [8][9][10][11] goes under the name of "the modern theory" of the orbital magnetization [see Eq. (12)]. Unlike the bulk polarization, m 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 ∇ × m pers (t, x). We find that, after coarse-graining in time, this contribution can be expressed as the curl of an additional term m to the orbital magnetization so that m in Eq. (3) is given by m = m pers + m. (6) Our main result is that 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 where V represents the system size and ϕ Bz is defined by Eq. (5) upon substitution φ → B z . This expression is welldefined 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 A for the detailed discussion. As comparison, in the presence of an external uniform magnetic field B = (0, 0, B z ) T , the instantaneous orbital magnetization m pers gives an energy shift T 0 dt(E tBz − E t0 )/T = −V m pers · B + O(B 2 ) of the many-body ground state. (Throughout this work we set = 1.) Accordingly, after the period T , the ground state acquires an additional phase proportional T V m pers · B. On the other hand, a non-zero value of m shows up as Berry phase T V m · B + O(B 2 ) acquired by the many-body ground state. For this reason, we name m geometric orbital magnetization. The bulk quantity T 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. 11 and 12 and find that geometric orbital magnetization consists of two contributions the topological contribution m 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 m non-top is written in terms of instantaneous Bloch states and energies [Eq. (52)]. The obtained expression for m of band insulators has a formal similarity with the expression for the magnetoelectric polarizability of three-dimensional band insulators 12 upon identification t/T ↔ k z /2π. To gain intuitive understanding of the two contributions (8), one can think of the topological piece m top to be originating from the Aharanov-Bohm contribution to the many-body Berry phase in the magnetic field. Thus m top describes the magnetization from electrons, whose positions are moving during the adiabatic process, as depicted with dashed lines and a) FIG. 2. Different contributions to orbital magnetization of twodimensional periodic adiabatic process with period T . a): Persistent current jpers within each unit cell produces the instantaneous orbital magnetization m. b): An adiabatic process where an electron trapped in a potential well whose center x = r(t) is moving along dashed curve. In the presence of an externally applied magnetic field, the many-body Berry phase ϕB z is given by the Aharonov-Bohm flux (the hatched area). c): Periodic boundary conditions are necessary when each of two potential wells comes back to its initial position after time T by passing though the seam. Two possible areas to define Aharonov-Bohm flux (the hatched one and the non-hatches one) differ by an integer flux quanta. d): Unit cell consists of a single anisotropic potential well that is spinning during adiabatic process. e): Two identical potential wells that exchange their positions after single adiabatic cycle. All adiabatic processes shown here have vanishing integrated current (4). arrows in Fig. 2b-c and e. Although Fig. 2a may look similar, j 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 m 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. [13][14][15][16][17] In Refs. 14 and 15 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. 16 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. 13 and 17 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 m, which for the case of band insulators is captured by non-abelian third Chern-Simons form, see Sec. IV C.
Analogous to the bulk polarization, crystalline symmetries can quantize topological geometric orbital magnetization m top . We show that, under certain crystalline symmetries, m top is related to recently discussed higher-order topological phases. [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] Among them, the topological insulators that exhibit quantized corner charges in the presence of certain crystalline symmetries attracted recently a lot of theoretical 34,35 and experimental 36,37 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. [38][39][40] Recent work in Ref. 41 revealed a connection between higher-order topological insulators 18-33 protected by roto-inversion symmetries and adiabatic processes that involve topological insulators with quantized corner charges. In Sec. III E we show that the adiabatic processes discussed by van Miert and Ortix 41 are characterized by quantized geometric orbital magnetization, and we relate the value of T m top z 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 E 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 E calculation of magnetization induced by rotation of an insulator. A long time ago, 42,43 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 42 argued that rotation creates a torque that acts to align local magnetic moments with rotation axis. This torque gives rise to magnetization M = χΩ/γ, where χ is paramagnetic susceptibility, Ω is rotation frequency and γ is electron gyromagnetic ratio. Barnett's measurement of this effect 42 provided first accurate measurement of electron gyromagnetic ratio. We calculate m for this model, which, as seen from Eq. (7), is also proportional to rotational frequency Ω = 2π/T and estimate quantum correction to Barnett effect. Electron contribution to m has both topological and non-topological piece, but since the system is uncharged, we find that electron contribution to m top is canceled by corresponding ionic contribution. Thus resulting 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" 13,17 or "shaking" 14-16 of the whole system. 44 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. 3-6, 8-11. To simplify notations, we assume primitive lattice vectors of the square lattice type, but this general framework is not restricted to this special choice.

A. Modern theory
Let us denote by ψ tkn (x) = (a/ √ V )e ik·x u tkn (x) the instantaneous Bloch function of n-th occupied band, satisfying h t |ψ tkn = ε tkn |ψ tkn . 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 R and reciprocal lattice vector G 45 According to the modern theory, the bulk polarization density p(t) is given by where e (< 0) is the electric charge. The sum over k can be replaced with the integral V d 2 k (2π) 2 over the first Brillouin zone. Similarly, the orbital magnetization density m pers (t) is given by where h tk ≡ e −ik·x h t e ik·x . The ambiguity in Eq. (11) can be seen by a smooth gauge transformation |u tkn = m |u tkm (U 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 R.
In terms of the Wannier function, p(t) is the deviation of the Wannier center from R, i.e., p(t) = e a 2 When the origin of unit cell is changed by δ, we find m pers (t) = m pers (t).
Namely, p(t) depends on the specific choice of the origin, while m pers (t) does not. Therefore, it is not p(t) itself but rather the change ∆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 φ i in this formulation is encoded in the relation whereP i is the polarization operator (see Ref. 46 for example).

B. Topological response
To discuss the physical consequence of ∆p(t), let us recall first the topological linear response in (1 + 1) dimension 1 that holds at a mesoscopic scale after coarse-graining Here, x µ (µ = 0, 1) represents (t, x), and j µ corresponds to (n, j). A θ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 (θ(t)) is the 1D version of p(t) in Eq. (11). This response is derived starting from the Chern-Simons theory j µ = C1 2π νλ ε µνλ ∂ ν A ex λ in (2 + 1) dimensions that describes the response toward an external field A ex and reducing the dimension to (1 + 1) dimensions. The parameter θ in Eq. (19) is a slowly varying field interpolating between two different systems. For example, an adiabatic time evolution θ(t) induces the bulk current j(t, x) = ∂ t P 1 (θ(t)). The bulk charge transfer from t = 0 to t = T is thus given by Similarly, a transition of one 1D system to another can be described by θ(x), giving rise to a charge density n(t, x) = −∂ x P 1 (θ(x)). Therefore, the total charge Q edge accumulated to the bound- . For a given θ that specifies P 1 (θ) as a continuous function of t and x, even the integer part of Q edge is well-defined. However, only the fractional part of Q 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 n is The analog of Eq. (19) in (3 + 1) dimensions reads 47 Here, µ, ν, ρ, λ = 0, 1, 2, 3. Again, A θk is defined by occupied Bloch states. This response is derived from the Chern- dimensions by a dimensional reduction. This topological response implies, for example, the magnetoelectric effect 12,47,48

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).

A. Berry phases in adiabatic process
Suppose we are interested in the expectation value of the quantityX, given byX for some parameter in the Hamiltonian. For example, in the case of the averaged current operator, can be identified with the solenoidal flux φ [see Eq. (1)]. Likewise, for the orbital magnetization we use the external magnetic field B z . Now suppose that the HamiltonianĤ t has a periodic adiabatic dependence on t, and let |Φ t be the instantaneous ground state with the energy eigenvalue E t . We assume an excitation gap ∆ t and the time-dependence of the Hamiltonian must be slow enough so that ∆ t T 1. Using the density matrixρ t = |Ψ t Ψ t | obeying the time-dependent Schrödinger equation ∂ tρt = −i[Ĥ t ,ρ t ], we express the time-average of the expectation value ofX t as In the absence of the time-evolution, the density matrix is identical to |Φ t Φ t |. It acquires contributions from excited statesĤ t |Φ M t = E M t |Φ M t due to the time evolution. To the lowest-order perturbation theory with respect to (∆ t T ) −1 , the relevant matrix elements are given by 1,46 Now we plugX t ≡ ∂ Ĥ t | =0 and make use of the Sternheimer identity which follows by differentiatingĤ t |Φ t = E t |Φ t with respect to at = 0. In our notation,Ĥ t | =0 =Ĥ t and |Φ t | =0 = |Φ t . Combining these equations, we find where is the Berry curvature in (t, ) space. Further assuming the periodicity in time we arrive at our general expression where X inst is the time average of the expectation value using the instantaneous ground state and X 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 whereP t = |Φ t Φ t | is the projector onto the many-body ground state,Ĝ t = (z −Ĥ 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 = B z . Then X geom in Eq. (35) gives the geometric orbital magnetization m in Eq. (7), while X inst is the persistent current contribution. Note the additional minus sign because ofM z = −∂ BzĤBz . Previously, Refs. 13 and 17 considered the Berry curvature in (t, B z ) space to describe orbital magnetization induced by rotation of molecules. See also examples in Sec. IV D and V A 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 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π. 49,50 However, each eigenstate ofĤ Bz may not be analytic as a function of B z despite the fact that the magnetic field B z = 2π/L x L y itself can be made small for a large systems. The expression Eq. (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 12,51 . However, to our knowledge there is no general proof for gapped interacting systems.
Before closing this subsection, let us make a few remarks on the (first) Chern number in the (t, φ x ), (t, φ y ), and (φ x , φ y ) space. Those who are happy to assume that they are zero can skip to Sec. III B.
When periodic boundary conditions are applied in at least one direction, corresponding solenoidal flux in the presence of an external magnetic field can be defined by specifying a big loop on torus, for example, by the dashed curves on Fig. 1. Solenoidal fluxes, defined in this way, need to be fixed when evaluating the derivative in Eq. (7), unless Chern numbers in (t, φ x ) and (t, φ y ) spaces vanish. To better understand the origin of these requirements, consider a modification of the adiabatic process in Fig. 2c, with only one instead of two potential wells. In this modified example, the Chern number in (t, φ y ) space does not vanish. It is easy to check that m z from Eq. (7) is proportional to Aharonov-Bohm flux defined by the area between the big loop, that defines φ y (that is being kept constant), and the trajectory of adiabatic process.
Lastly, if the system has a non-zero quantized Hall conductance, periodic boundary conditions need to be assumed since otherwise the ground state is not gapped. Aside from above mentioned remarks regarding periodic boundary conditions in the presence of magnetic field, we do not see any issue in applying Eq. (7) to the systems with non-zero Chern number, although we have not constructed any concrete examples to test this statement. Therefore we leave this question for future works. We remark that recent work in Ref. 52 found interesting bulk charge pumping as applied magnetic flux is varied in the systems with non-zero Chern number under periodic boundary conditions in only one direction.

B. Noninteracting systems
Let us apply this general expression to noninteracting electrons described by the quadratic Hamiltonian We label single-particle states in such a way that ε t n+1 ≥ ε t n for all n = 1, 2, · · · . We also assume a finite gap ∆ = ε t N +1 − ε t N between N -th and (N + 1)-th levels. Then the N -particle ground state can be written as For a later purpose, we allow for a unitary transformation among the occupied levelŝ 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 appropriately. (See Sec. III C for an example). After all we find that the many-body Berry phase ϕ is given by the sum of singleparticle Berry phases ϕ of occupied levels Therefore, we get the following expressions for single-particle problems. The latter two expressions are basis-independent where P t = N =1 |ψ t ψ t | is projector onto occupied single-particle states, h t = n ε tn |γ tn γ tn | is singleparticle Hamiltonian, g t = (z −h t ) −1 is single-particle Green function, and the integration contour encloses all the occupied states at z = ε tn (n = 1, 2, · · · , N ).
For the orbital magnetization, we again set = 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. 49,50 For single-particle systems, assuming symmetric gauge A ex (x) = B × x/2, the magnetic periodic boundary conditions can be taken into account explicitly by restricting the form of the projectorP tBz to 12,51 x 1 |P tBz |x 2 =P tBz (x 2 , x 1 )e ieBzx1×x2·ẑ/2 , where P tBz (x 1 , x 2 ) is an arbitrary N × N matrix function (not necessarily projector) that satisfies P tBz (x 1 + R, x 2 + R) = P tBz (x 1 , x 2 ), where R is an element of Bravais lattice. The expression forP tBz (x) can be found perturbatively in B z , 12,51 which, after substituting back to Eq. (41), yields an expression for the Berry phase and m. The second option is to apply a spatially modulating magnetic field as we discuss below.

C. 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 11,12 in order to avoid changing of the boundary condition as discussed in the previous subsection. We use the vector potential with e z ≡ (0, 0, 1) T , q ≡ 2π/L, and f (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 and |∂ w t nR | =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 → 0. The single-particle Berry phase in this basis, summed over occupied bands, takes the following form If we further sum over R, or equivalently if we work in the Bloch basis |ψ tnk , we get 0 reflecting the fact that for bulk systems Fourier component m z (q) vanishes for q = 0. Thus, care must be taken to correctly read off local contribution to m z -an unintentional integration over R of a term proportional to f (R) makes it impossible to find the correct value of m z . The rest calculation follows the appendix in Ref. 12. Upon taking the limit q → 0, we find This last expression can precisely be expressed as the sum of two terms, m top + m non-top . The topological piece m top reads The Berry connection (A K ) n,m ≡ −i u Kn |∇ K u Km is defined using occupied Bloch states as a function of K ≡ (t, k). The smoothness and the periodicity of A K are assumed in the integral in Eq. (51). Such a choice is possible only when both the pumped charge through the bulk 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 Here, h tk is Bloch Hamiltonian, P tk = n∈occ |u tnk u tnk | is the projector onto occupied bands at k, g tk = (z − h tk ) −1 is Bloch's Green function, the curly brackets denote sym-

D. Topological contribution from response theory
Here we give an alternative, easier derivation of m 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 47 : where is the Berry curvature in the (k x , k y , θ, φ) space and θ and φ are two slowly varying fields: θ(t) denotes an adiabatic and periodic time dependence and φ(x) describes a smooth interface of domains (Fig. 3). In this setting, we find so that This reproduces Eq. (50). In the derivation we used the relation It is important to note that j(t, x) itself cannot be written as a curl of a vector field -Eq. (56) holds only after the time convolution (or equivalently the time average).
The above derivation relies on the connection (3) between m top and topological edge current in adiabatically driven twodimensional systems. To see this more concretely, let us consider the boundary of two regions with φ 0 ≡ φ(x 0 ) and φ 1 ≡ φ(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 φ 0 and φ 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. 1 To capture the fractional bulk contribution to the edge current, one can separately compute m z (x 0 ) and m z (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 m z (x 0 ) and m z (x 1 ) = m z (x 1 ) (Fig. 3) is given by (57) Notice that I edge of two adjacent edges may differ by an integer. To see this formally, let us consider a charge flow ∆Q corner into a corner surrounded by a closed curve x α with x 1 = x 0 (see Fig. 3). The net charge flow in the process is given by the second Chern number For example, when the corner is formed by two edges along x and y directions, we have meaning that the charge transfer along two intersecting edges can only differ by an integer multiple of e. Clearly, ∆Q 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 (Fig. 3).

E. Symmetry constraints and corner charge
Here we consider adiabatic process of two-dimensional systems constrained by certain symmetries that quantize T m z . We show that, if the symmetry allows one to define the bulk contribution to quadrupole moment, the quantized quadrupole moment is equal to T m z . 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. 53 For concreteness, let us consider the four-fold rotation C 4 mapping x = (x, y, 0) to C 4 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 ∆Q corner , see Fig. 4a. We consider an arbitrary interpolation between the system of interest h 0k and the reference system h T /2 k that has no corner charge. The second half of the cycle is performed in a C 4 -symmetric manner 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 m z : This does not mean that T m z vanishes since it is defined only mod 1. Thus in the presence of C 4 symmetry T m z is quantized either 0 or e/2 mod e. Note that the inversion symmetry, for example, also quantizes T m z 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. 4b). Now we show that the parity of the corner charge accumulation ∆Q corner is actually a bulk topological invariant for symmetric adiabatic processes satisfying Eq. (60). 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 edge y = −I edge x and that Furthermore, Fig. 4b tells us that the corner charge accumulation during the symmetric process is ∆Q corner = 2q corner . Therefore, We will discuss an example of quadrupole insulators with P 3 = e/2 in Sec. IV C using this result. On the other hand, a C 4 -symmetric phase that hosts a corner charge of q corner = e/4 were recently reported. 35 The fact that ∆Q corner ∈ 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. 53, the C 4symmetric adiabatic process h tk considered above, can be viewed as a 3D topological insulator protected by the rotoinversion symmetry IC 4 upon identification k z = 2π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 xyplane 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 0 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, m pers and m. Additionally, we want to consider examples where there is only topological geometric magnetization, Sec. IV C, only non-topological geometric magnetization, Sec. IV D, and both topological and non-topological contributions, Sec. IV E. To keep the discussions simple while capturing the relevant physics, we focus on isolated orbitals without any overlap between them.
A. Bloch functions in the localized limit Let us consider a time-dependent deep potential v 0 t (x) centering at x = r(t) that accommodates at least one bound state. Let φ 0 t (x) be the wavefunction of the instantaneous lowestenergy bound state, satisfying h 0 Here A ex t (x) describes an external field. In these expressions, the superscript 0 implies the quantities for an isolated orbit. When the potential v 0 t (x) is deep enough, φ 0 t (x) should be well-localized around x = r(t) with the localization length ξ a. Hence, we assume that and that both φ 0 t (x) and v 0 t (x) decays fast enough, i.e., |v 0 t (x)|, |φ 0 t (x)| → 0 as |x − r(t)| ξ. With these building blocks, we construct a periodic potential and the cell-periodic Bloch state.
We assume that A ex (x) respects the periodicity, i.e., A ex (x − R) = A ex (x). Then, as far as φ 0 is an eigenstate of the periodic Hamiltonian with a completely flat band dispersion ε tk = ε 0 t .

B. Polarization and instantaneous magnetization
Let us first demonstrate the modern theory formula for the polarization and the orbital magnetization by deriving p and 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 change and current distribution from a single orbit φ 0 t (x) can be written as We introduce vector fields p 0 t (x) and m 0 t (x) such that The existence of such m 0 t (x) is guaranteed by the divergencefree nature of the instantaneous current density j 0 t (x). The current density induced by the adiabatic motion of r(t) is captured by j 0 t (x) in Eq. (81) whose divergence may not vanish. We assume both p 0 t (x) and m 0 t (x) decay rapidly for |x − r(t)| > ξ, 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 and analogously for n, p, and 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. 54). Here we choose the Gaussian g(x) = (πR 2 ) −1 e −|x| 2 /R 2 (R a) We do the same for other quantities. Relations such as j(t, x) = ∇ × m(t, x) are preserved by the convolution. Because the convolution is identical to the average for the periodic distribution, we find n(t, x) =n = e a 2 , j(t, x) = 0, 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 nonoverlapping assumption of φ 0 t (x), it can be readily shown that the formula in Eqs. (11), (12) for the Bloch function (67) can be simplified to

C. Topological geometric magnetization
Next, we discuss the topological geometric contribution m top for this model. To this end, suppose that the position of the potential minimum r(t) adiabatically moves as a function of t ∈ [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 m top in Eq. (50) to the Bloch function (67). Thanks to the non-overlapping assumption, the vector potential A (t,k) is k-independent: . Plugging this into Eq. (51), we find where S r represents the area enclosed by the orbit of r(t) in one cycle. Therefore, Observe the analogy to 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 x = r(t) induces a local current distribution It becomes divergence-free if averaged over one period As we have seen above, the sum of such microscopic currents from each unit cell produces the bulk magnetization This agrees with Eq. (80) because the integral in the parenthesis is precisely 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 x = r n (t) (n = 1, 2, · · · ) in each unit cell, which are adiabatically varied as a function of t ∈ [0, T ]. This time, each orbit is allowed to form an open curve, as far as the total polarization p(t) = (e/a 2 ) n r n (t) satisfies p(T ) = p(0). Under such an assumption, we find that 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. 34. For the wallpaper group p4, there exist three spacial Wyckoff positions: the unit cell origin at x a = (0, 0), the center of the plaquette at x b = (a/2, a/2), and the center of bonds at x c = (a/2, 0), (0, a/2). 55 In the nontrivial phase, the two occupied Wannier orbitals locate at x b , while in the trivial phase they are at 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 tk 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. 34 A variant of this adiabatic process was also discussed in Ref. 53.

D. Non-topological geometric magnetization
In this example we first consider a single electron in an anisotropic and rotating two-dimensional well, 44 see Fig. 2d. We assume a harmonic confining potential, i.e., Hamilto- where x t ≡ (x cos Ωt + y sin Ωt, −x sin Ωt + y cos Ωt) and Ω = 2π/T . To obtain the geometric orbital magnetization for this model, we consider external magnetic field B = B z e z described by the vector potential A ex (x) = B × x/2.
where N is the normalization factor and ω c ≡ eB z /m is the cyclotron frequency. The Berry phase ϕ Bz during the adiabatic process t ∈ [0, T ] is From Eq. (7) it follows that the adiabatic process (87) has nonzero geometric orbital magnetic moment m 0 Now we construct the Bloch function (67) using φ 0 t (x) as the building block and compute the geometric orbital magnetization 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 tk is given by Eq. (68) with v 0 t (x) given by Eq. (87) and A ex = 0. We assume that there is no overlap between wavefunctions belonging to different unit cells as before. (When ω x , ω y are large enough, such an assumption is valid at least for relevant low-energy states.) Bloch wavefunctions read where n ≡ (n x , n y ) labels energy levels of two-dimensional the anisotropic harmonic oscillator and n = (0, 0) corresponds to the ground state in Eq. (88) with ω c = 0. Substituting above expressions to Eqs. (50) and (52), we find

E. 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. 42,43 The Barnett effect predicts magnetization χ/γΩ, where χ is the paramagnetic susceptibility, γ is the electron gyromagnetic ratio, and Ω is the rotation frequency. Since the rotation axis does not necessarily coincide with potential well minima we have v 0 t (x − r(t)) with v 0 t from Eq. (87) and r(t) = (R cos Ωt, R sin Ωt, 0) T , where R is the distance of the potential well minima to the rotation axis. The lowest-energy instantaneous wavefunction φ t (x) can be obtained from Eq. (88) by performing gauge transformation ×e − m √ (ωx +ωy ) 2 +ω 2 c (ωx (x t −R) 2 +ωy y 2 t ) 2(ωx+ωy ) .
As compared to Eq. (89), the Berry phase ϕ Bz during the adiabatic process t ∈ [0, T ] acquires an additional contribution eB z πR 2 from the Aharonov-Bohm phase. Therefore electrons contribute to the following geometric orbital magnetization The first term can be identified with m top z 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 m non-top z = m 0 z remains since ions are much more localized compared to electrons. Assuming anisotropy ω x /ω 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 χ ∼ 10 −5 . For comparison, paramagnets have typically magnetic susceptibility χ ∼ 10 −3 − 10 −5 . 57

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 13,17 and translating the whole system. 14,16,44 We consider many-body systems under the open boundary condition in two spatial dimensions. We start with a timeindependent HamiltonianĤ that can contain arbitrary interactions. The total charge, current, polarization, and orbital magnetization operator for this Hamiltonian can be written aŝ N = V d 2 xn(x),Ĵ = d 2 xĵ(x),X = V d 2 xxn(x), M = (1/2) V d 2 xx×ĵ(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Ĥ and its energy by |Φ and E, respectively.
To compute the many-body Berry phase, letĤ B be the Hamiltonian with the vector potential in the symmetric gauge A ex (x) = (1/2)B × x with B = B z e z . Expanding to the linear order in B z and usingĵ(x) = −∂ A(x)Ĥ , we get Therefore, the ground state ofĤ Bz to the leading order in B z can be expressed as HereQ ≡ 1 − |Φ Φ| is the projector onto excited states.

A. Rotation
Let us first introduce the time-dependence to the problem by the rotationĤ where Ω = e z Ω is the rotation frequency andL is the angular momentum operator. For the time-dependent HamiltonianĤ t , the orbital magnetization We evaluate the instantaneous contribution m pers and the geometric contribution m to the orbital magnetization via the formulae in Eqs. (34) and (35). The instantaneous contribution is given by the instantaneous ground state |Φ t = e −iLzΩt |Φ The geometric contribution is given by the many-body Berry phase. Since the instantaneous ground state of the Hamil-tonianĤ tBz ≡ e −iLzΩtĤ Bz e iLzΩt is given by |Φ tBz = e −iLzΩt |Φ Bz , we have This is the expectation value ofL z Ω in the presence of the perturbation −m z B z in Eq. (96). Using Eq. (97), we get m = Φ|L z ΩQ 1 We verify these results by solving time-dependent problem. The solution to the time-dependent Schrödinger equation i∂ t |Ψ t =Ĥ t |Ψ t can be readily constructed using the ground state |Φ Ω of the time-independent Hamiltonian The solution that is smoothly connected to the ground state in the static limit Ω → 0 reads The time-average of the orbital magnetization is thus given by This is the expectation value ofM in the presence of the perturbation −L z Ω as in Eq. (103). The first-order perturbation theory with respect to Ω gives m = m pers + Φ|L z ΩQ 1 This is precisely m pers + 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 energyF 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 timedependent Hamiltonian bŷ whereT t = e −iP ·r(t) is the translation by amount r(t) andP is the momentum operator. ForĤ t the orbital magnetization operator becomeŝ where the second term in the parenthesis is due to the change of the origin. The instantaneous ground state |Φ t =T t |Φ gives m pers as in Eq. (100), where we used Φ|Ĵ |Φ = 0.
Next, we compute the geometric contribution via the manybody Berry phase. In the presence of magnetic field translation operatorT tBz ≡T Bz (r(t)) becomes translation followed by gauge transformation 49,50 ∂ tTtBz ≡ −i(P + e 2 B ×X) · ∂ t r(t)T tBz .
The instantaneous ground state ofĤ tBz ≡T tBzĤBzT † tBz is |Φ tBz =T tBz |Φ Bz , thus the many-body Berry phase reads Here, S r ≡ 1 2 r × dr represents the area swept by r(t) in one cycle. In the derivation, we used Therefore, when the whole system is translated, the geometric contribution to the orbital magnetization captures the Aharonov-Bohm phase To verify the above results, we consider the time-dependent Schrödinger equation i∂ t |Ψ t =Ĥ t |Ψ t . An approximate solution is given by |Ψ t =T t |Φ tr , where |Φ tr is the instantaneous ground state of the HamiltonianĤ tr ≡Ĥ −P · ∂ t r(t). Therefore, the time-average of the orbital magnetization is dt r(t) × ∂ t r(t).
In the adiabatic limit, this reproduces m pers + m in Eqs. (100) and (113). In the derivation we used Φ tr |Ĵ |Φ tr = eN ∂ t r(t) for the ground state ofĤ tr .

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 3-6 and orbital magnetization [8][9][10][11] 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 m. This new contribution is defined only after performing the timeaverage 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 m in terms of instantaneous Bloch states and energies. Interestingly, we find that for band insulators m = m top + m non-top consists of two pieces, where topological piece m 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. 54 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, 58 it would be interesting to see to which extent our findings can be generalized to such systems.
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. 34,35 In fact, it is a topic of current research whether higher moments can be established as bulk quantities in general. [38][39][40]59,60 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 m top z in systems with proper symmetries that allow bulk definition of the quadrupole moment. 40 In this work we succeeded in separating 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, 61,62 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, 12 thus the analogous "separation question" arises also in that context. Additionally, defining quantized quadrupole moment for interacting systems is one of the open questions. [38][39][40] Since for band insulators we find connection between m top z and quantized quadrupole moment, separating m top z 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. [13][14][15][16][17]