Theory of orbital magnetoelectric response

We extend the recently-developed theory of bulk orbital magnetization to finite electric fields, and use it to calculate the orbital magnetoelectric response of periodic insulators. Working in the independent-particle framework, we find that the finite-field orbital magnetization can be written as a sum of three gauge-invariant contributions, one of which has no counterpart at zero field. The extra contribution is collinear with and explicitly dependent on the electric field. The expression for the orbital magnetization is suitable for first-principles implementations, allowing to calculate the magnetoelectric response coefficients by numerical differentiation. Alternatively, perturbation-theory techniques may be used, and for that purpose we derive an expression directly for the linear magnetoelectric tensor by taking the first field-derivative analytically. Two types of terms are obtained. One, the `Chern-Simons' term, depends only on the unperturbed occupied orbitals and is purely isotropic. The other, `Kubo' term, involves the first-order change in the orbitals and gives isotropic as well as anisotropic contributions to the response. In ordinary magnetoelectric insulators both terms are generally present, while in strong Z2 topological insulators only the former type is allowed, and is quantized. In order to validate the theory we have calculated under periodic boundary conditions the linear magnetoelectric coupling for a 3-D tight-binding model of an ordinary magnetoelectric insulator, using both the finite-field and perturbation-theory expressions. The results are in excellent agreement with calculations on bounded samples.


Introduction
In insulating materials in which both spatial inversion and time-reversal symmetries are broken, a magnetic field B can induce a first-order electric polarization P , and conversely an electric field E can induce a first-order magnetization M [1,2]. This linear magnetoelectric (ME) effect is described by the susceptibility tensor where indices label spatial directions. This tensor can be divided into a "frozen-ion" contribution that occurs even when the ionic coordinates are fixed, and a "latticemediated" contribution corresponding to the remainder. Each of these two contributions can be decomposed further according to whether the magnetic interaction is associated with spins or orbital currents, giving four contributions to α in total. All of those contributions, except the frozen-ion orbital one, are relatively straightforward to evaluate, at least in principle, and ab initio calculations have started to appear. For example, the lattice-mediated spin-magnetization response was calculated in [3] for Cr 2 O 3 and in [4] for BiFeO 3 (including the strain deformation effects that are present in the latter), and calculations based on the converse approach (polarization response to a Zeeman field) were recently reported [5]. One generally expects the lattice-mediated couplings to be larger than the frozen-ion ones, and insofar as the spin-orbit interaction can be treated perturbatively, interactions involving spin magnetization are typically larger than the orbital ones. However, we shall see that there are situations in which the spin-orbit interaction cannot be treated perturbatively, and in which the frozen-ion orbital contribution is expected to be dominant. Therefore, it is desirable to have a complete description which accounts for all four contributions.
The frozen-ion orbital contribution is, in fact, the one part of the ME susceptibility for which there is at present no satisfactory theoretical or computational framework, although some progress towards that goal was made in two recent works [6,7]. Following Essin et al. [7] we refer to it as the "orbital magnetoelectric polarizability" (OMP). For the remainder of this paper, we will focus exclusively on this contribution to (1), and shall denote it simply by α. Accordingly, the symbol M will be used henceforth for the orbital component of the magnetization.
The question we pose to ourselves is the following: what is the quantum-mechanical expression for the tensor α of a generic three-dimensional band insulator? We note that the conventional perturbation-theory expression for α [8,9] does not apply to Bloch electrons, as it involves matrix elements of unbounded operators. The proper expressions for P [10] and M [11,12,13,14] in periodic crystals have been derived, but so far only at B = 0 and E = 0 respectively. The evaluation of equation (1) remains therefore an open problem.
Phenomenologically, the most general form of α is a 3 × 3 matrix where all nine components are independent. Dividing it into traceless and isotropic parts, the latter is conveniently expressed in terms of a single dimensionless parameter θ as The presence of an isotropic ME coupling is equivalent to the addition of a term proportional to θE · B to the electromagnetic Lagrangian. Such a term describes "axion electrodynamics" [15] and (2) may therefore also be referred to as the "axion OMP." The electrodynamic effects of the axion field are elusive (in fact, the very existence of α θ was debated until recently: see [16,17] and references therein). For example, in a finite, static sample cut from a uniform ME medium those effects are only felt at the surface [15,18]. In particular, α θ gives rise to a surface Hall effect [19]. An essential feature of the axion theory is that a change of θ by 2π leaves the electrodynamics invariant [15]. The profound implications for the ME response of materials were recognized by Qi et al. [6], and discussed further by Essin et al. [7]. These authors showed that there is a part of the isotropic OMP which remains ambiguous up to integer multiples of 2π in the corresponding θ until the surface termination of the sample is specified. For example, a change by 2πn occurs if the surface is modified by adsorbing a quantum anomalous Hall layer. Hence this particular contribution to θ can be formulated as a bulk quantity only modulo a quantum of indeterminacy, in much the same way as the electric polarization P [10,20]. A microscopic expression for it was derived in the framework of single-particle band theory by the above authors. It is given by the Brillouin-zone integral of the Chern-Simons form [21] in k-space, which is a multivalued global geometric invariant reminiscent of the Berry-phase expression for P [10]. We denote henceforth this "geometric" contribution to the OMP as the Chern-Simons OMP (CSOMP).
A remarkable outcome of this analysis is the prediction [6] of a purely isotropic "topological ME effect," associated with the CSOMP, in a newly-discovered class of time-reversal invariant insulators known as Z 2 topological insulators [22,23,24]. As a result of the multivaluedness of θ, the presence of time-reversal symmetry in the bulk, which takes θ into −θ, is consistent with two solutions: θ = 0, corresponding to ordinary insulators, and θ = π, corresponding to strong Z 2 topological insulators. ‡ The latter case is non-perturbative in the spin-orbit interaction, and θ = π amounts to a rather large ME susceptibility (in Gaussian units it is 1/4π times the fine structure constant, or ∼6×10 −4 , to be compared with ∼1×10 −4 for the total ME response of Cr 2 O 3 at low temperature [25]).
It is not clear from these recent works, however, whether the isotropic CSOMP constitutes the full OMP response of a generic insulator. It does appear to do so for the tight-binding model studied in [7], whose ME response was correctly reproduced by the Chern-Simons expression even when the parameters were tuned to break time-reversal and inversion symmetries (i.e., for generic θ not equal to 0 or π). On the other hand, other considerations seem to demand additional contributions. For example, it is not difficult to construct tight-binding models of molecular crystals in which it is clear that the OMP cannot be purely isotropic.
In this work we derive, using rigorous quantum-mechanical arguments, an expression for the OMP tensor α of band insulators, written solely in terms of bulk quantities (the periodic Hamiltonian and ground state Bloch wavefunctions, and their first-order change in an electric field). We restrict our derivation to non-interacting Hamiltonians, as the essential physics we wish to describe occurs already at the single-particle level. We find that in crystals with broken time-reversal and inversion symmetries there are, in addition to the CSOMP term discussed in [6,7], extra terms which generally contribute to both the trace and the traceless parts of α.
Our theoretical approach closely mimics one type of ME response experiment: a finite electric field E is applied to a bounded sample, and the (orbital) magnetization is calculated in the presence of the field. Then the thermodynamic limit is taken at fixed field. This key step in the derivation must be done carefully, so that crucial surface contributions are not lost in the process, and here we follow the Wannier-based approach of references [12,13], adapted to E = 0. Finally the linear response coefficient α da = ∂M a /∂E d is extracted in the limit that E goes to zero.
In a concurrent work by Essin, Turner, Moore, and one of us [26] an alternative approach was taken, which is closer in spirit to the calculation in [10] of the change in polarization as an integrated current: the adiabatic current induced in an infinite crystal by a change in its Hamiltonian in the presence of a magnetic field is computed, and then expressed as a total time derivative. The two approaches are complementary and lead to the same expression for α, illuminating it from different angles.
The paper is organized as follows. In section 2 we derive the bulk expression for M(E), and reorganize it into three gauge-invariant contributions, one of which yields directly the CSOMP response. The gauge-invariant decomposition of M (E) is done at first in k-space for periodic crystals, and then also for bounded samples working in real space. In section 3 we derive a k-space formula for the OMP tensor α by taking analytically the field-derivative of M (E). Numerical tests on a tight-binding model of a ME insulator are presented at appropriate places throughout the paper in order to validate the bulk expressions for M (E) and α. In Appendix A we describe the tightbinding model, as well as technical details on how the various formulas are implemented on a k-point grid. Appendix B and Appendix C contain derivations of certain results given in the main text.

Preliminaries
The orbital magnetization M is defined as the orbital moment per unit volume, Here e > 0 is the magnitude of the electron charge, V is the sample volume, and |ψ i are the occupied eigenstates. While this expression can be directly implemented when using open boundary conditions, the electronic structure of crystals is more conveniently calculated and interpreted using periodic boundary conditions, in order to take advantage of Bloch's theorem. This poses however serious difficulties in dealing with the circulation operator r × v, because of the unbounded and nonperiodic nature of the position operator r. These subtle issues were fully resolved only recently, with the derivation of a bulk expression for M directly in terms of the extended Bloch states [11,12,13,14].
In previous derivations the crystal was taken to be under shorted electrical boundary conditions. We shall extend the derivation given in [12,13] to the case where a static homogeneous electric field E is present, so that the full Hamiltonian reads The derivation, carried out for an insulator with N valence bands within the independent-particle approximation, involves transforming the set of occupied eigenstates |ψ i of H into a set of Wannier-type (i.e., localized and orthonormal) orbitals |w i and expressing M (E) in the Wannier representation. This is done at first for a finite sample cut from a periodic crystal, and eventually the thermodynamic limit is taken at fixed field. Before continuing, two remarks are in order. First, the assumption that it is possible to construct well-localized Wannier functions (WFs) spanning the valence bands is only valid if the Chern invariants of the valence bandstructure vanish identically [27]. This requirement is satisfied by normal band insulators as well as by Z 2 topological insulators, but not by quantum anomalous Hall insulators [28], which thus far remain hypothetical. Second, because of Zener tunnelling, an insulating crystal does not have a well-defined ground state in a finite electric field. Nonetheless, upon slowly ramping up the field to the desired value, the electron system remains in a quasistationary state which is, for all practical purposes, indistinguishable from a truly stationary state. This is the state we shall consider in the ensuing derivation. As discussed in [29,30], it is Wannier-and Bloch-representable, even though the Hamiltonian (4) is not lattice-periodic.

k-space expression
Our derivation of a k-space (bulk) expression for M (E) is carried out mostly in real space, using a Wannier representation. It is only in the last step that we switch to reciprocal space, by expressing the crystalline WFs |Rn in terms of the cell-periodic Bloch functions |u nk via [31] where R is a lattice vector, V c is the unit-cell volume, [dk] ≡ d 3 k/(2π) 3 , and the integral is over the first Brillouin zone. We begin with a finite sample immersed in a field E, divide it up into an interior region and a surface region, and assign each WF to either one. The boundary between the two regions is chosen in such a way that the fractional volume of the surface region goes to zero as V → ∞, but deep enough that WFs near the boundary are bulk-like. Following [12,13], equation (3) for the orbital magnetization can then be rewritten as an interior contribution plus a surface contribution, denoted respectively as the "local circulation" (LC) and the "itinerant circulation" (IC). Remarkably, in the thermodynamic limit both can be expressed solely in terms of the interior-region crystalline WFs, or equivalently, in terms of the bulk Bloch functions, as shown in the above references at E = 0 and below for E = 0. Specifically, we shall show that where is the contribution from the interior WFs, is the part of the surface contribution coming from the zero-field Hamiltonian, and is the part of the surface contribution coming from the electric field term in the Hamiltonian (4). In the above expressions γ = −e/(2 c), and A mnk is the Berry connection matrix defined in equation (14) below.
Having stated the result we now present the derivation, starting with the interior contribution M LC . Using [r i , r j ] = 0, the velocity operator v = (i/ )[H, r] becomes (i/ )[H 0 , r], so that the circulation operator r × v is unaffected by the electric field. It immediately follows that the local circulation part M LC is given in terms of the fieldpolarized states |u nk by the same expression, equation (6b), as was derived in [13] for the zero-field case.
Consider now the contribution M IC = M IC,0 + M IC,E from the surface WFs |w s . For large samples it takes the form [13] where N c is the number of crystal cells of volume V c , r s = w s |r|w s , and v s = w s |v|w s = 2 Im w s |rH|w s .
Note that H|w s already belongs to the occupied manifold spanned by P = occ j |w j w j |, since we assume a (quasi)stationary state. Thus we can insert a P between r and H above, and using (4) we obtain where v 0 js = (2/ )Im[r sj H 0 js ] is the same as in [12,13] and The reasoning [12,13] by which M IC can be recast in terms of the bulk WFs |Rn relies on the exponential localization of the WFs and on certain properties of v 0 js (antisymmetry under j ↔ s and invariance under lattice translations deep inside the crystallite) which are shared by v E js . Hence we can follow similar steps as in those works, arriving at and similarly for M IC,0 a with v 0 substituting for v E . The latter is identical to the expression for M IC a valid at E = 0 [12,13], and upon converting to k-space becomes (6c).
Let us now turn to M IC,E a and write (12) as In order to recast this expression as a k-space integral it is useful to introduce the N ×N Berry connection matrix where ∂ b ≡ ∂/∂k b . It satisfies the relation [31,32] We also need which follows from (15). Using these two relations, (13) becomes and we arrive at (6d). The sum of equations (6b)-(6d) gives the desired k-space expression for M (E). In the limit E → 0 the term M IC,E vanishes, and equation (31) of [13] is recovered.
We have implemented (6b)-(6d) for the tight-binding model of Appendix A. Since for small electric fields M (E) differs only slightly from M (0), in order to observe the effect of the electric field we consider differences in magnetization rather than the absolute magnetization. Therefore, in all our numerical tests we evaluated the OMP tensor α da . With the help of (6b)-(6d) we calculated it as ∆M a /∆E d , using small fields E d = ±0.01. We then repeated the calculation on finite samples cut from the bulk crystal, using (3) in place of (6b)-(6d). Figure 1 shows the value of the zz and zy components of α plotted as a function of the parameter ϕ, the phase of one of the complex hopping amplitudes (see Appendix A for details). The very precise agreement between the solid and dashed lines confirms the correctness of the k-space formula. The same level of agreement was found for the other components of α.

Gauge-invariant decomposition
2.3.1. Periodic crystals Equation (6a) for M (E) is valid in an arbitrary gauge, that is, the sum of its three terms given by (6b)-(6d) -but not each term individuallyremains invariant under a unitary transformation among the valence-band states at each k. In order to make the gauge invariance of (6a) manifest, it is convenient to first manipulate it into a different form, given in terms of certain canonical objects which we now define. We begin by introducing the covariant k-derivative of a valence state [30], where Q k = 1 − P k and The covariant and ordinary derivatives are related by The generalized metric-curvature tensor is [31] Viewed as an N × N matrix over the band indices, F is gauge-covariant, changing as under the transformation (18). We also note the relation We shall make use of two more gauge-covariant objects, and which enter the relation Coming back to equations (6a)-(6d), for M LC a we use (27) and for M IC a we use (24), leading to where "tr" denotes the electronic trace over the occupied valence bands and we have dropped the subscript k. The second term can be rewritten using (To obtain this relation start from the generalized Schrödinger equation satisfied by the valence states at E = 0 [33], and multiply through by ∂ c u m |.) Let us define the quantities and The total magnetization is given by their sum Referring to (22) and (26) the first two terms read, in a more conventional notation, and These are the only terms that remain in the limit E → 0, in agreement with equation (43) of [13]. At finite field they depend on E implicitly via the wavefunctions. We now show that the term M CS , which gathers all the contributions with an explicit dependence on E, can be recast as To do so it is convenient to introduce the Berry curvature tensor Ω nm,ab = iF nm,ab − iF nm,ba = −Ω nm,ba , where F nm,ab was defined in (22). A few lines of algebra show that In order to go from (33)  where Ω nm,a = 1 2 ǫ abc Ω nm,bc is the Berry curvature tensor written in axial-vector form. This leads to The first term is parallel to the field, and can be rewritten with the help of (36): While not immediately apparent, the second term in (37) also points along the field. To see this, write where we suspended momentarily the implied summation convention. The last term vanishes because the factor ǫ abc forces d = a to equal either b or c, producing terms such as Im tr and restoring the summation convention, we arrive at (34d). Equations (34b)-(34d), which constitute the main result of this section, are separately gauge-invariant. For M LC and M IC this is apparent already from (31) and (32), whose integrands are gauge-invariant, being traces over gauge-covariant matrices. In contrast, equation (34d) for M CS only becomes invariant after taking the integral on the right-hand-side over the entire Brillouin zone (the integrand being familiar from differential geometry as the Chern-Simons 3-form [34,21]).
The Chern-Simons contribution (34d) has several remarkable features: (i) as already noted, it is perfectly isotropic, remaining parallel to E for arbitrary orientations of E relative to the crystal axes; (ii) being isotropic, it vanishes in less than three dimensions, which intuitively can be understood because already in two dimensions polarization must be in the plane of the system and magnetization must be out of the plane; (iii) for N > 1 valence bands it is a multivalued bulk quantity with a quantum of arbitrariness (e 2 /hc)E a , a fact that is connected with the possibility of a cyclic adiabatic evolution that would change (47a) below for θ by 2π [6].
We have repeated the calculation of the OMP presented in figure 1 using (34b)-(34d) instead of (6b)-(6d), finding excellent agreement between them. The electric field derivative of the decomposition (34a) gives the corresponding decomposition of the OMP tensor (1), where each term is also gauge-invariant. The zz components of these terms are plotted separately in figure 2.

Finite samples
It is natural to ask whether the gauge-invariant decomposition of the orbital magnetization given in equation (34a) can be made already for finite samples, before taking the thermodynamic limit and switching to periodic boundary conditions. Here P and Q = 1 − P are the projection operators onto the occupied and empty subspaces, respectively, and "Tr" denotes the electronic trace over the entire Hilbert space. These two expressions, which are manifestly gauge-invariant, remain valid at finite field, reducing to (34b) and (34c) in the thermodynamic limit. We now complete this picture for E = 0 by showing that the remaining contribution M CS = M − M LC − M IC can also be written in trace form, as We first recast the orbital magnetization (3) as and then subtract (41a) and (41b) from it to find, after some manipulations, Replacing H 0 with H − eE d r d and using QHP = 0, The imaginary part of the trace vanishes if any two of the indices b, c, or d are the same, and therefore d must be equal to a. Using the cyclic property we conclude that all non-vanishing terms in the sum over b and c are identical, leading to (41c). This part of the field-induced magnetization is clearly isotropic, with a coupling strength (see equation (2)) given by This expression can assume nonzero values because the Cartesian components of the projected position operator P rP do not commute [31].
We have used (41a)-(41c) to evaluate the OMP contributions α LC , α IC , and α CS for finite samples, finding excellent agreement with the k-space calculations using (34b)-(34d). As an example, the finite-sample results for the zz component are plotted as the symbols in figure 2.

Linear-response expression for the OMP tensor
In sections 2.2 and 2.3.1 expressions were given for evaluating M (E) under periodic boundary conditions. Used in conjunction with finite-field ab-initio methods for periodic insulators [36,29], they allow to calculate the OMP tensor by finite differences. Alternatively, the electric field may be treated perturbatively [33]. With this approach in mind, we shall now take the E-field derivative in (1) analytically and obtain an expression for the OMP tensor which is amenable to density-functional perturbationtheory implementation [37]. It should be kept in mind that in the context of selfconsistent-field (SCF) calculations the "zero-field" part of the Hamiltonian (4), does depend on E implicitly, through the charge density. As we will see, this dependence gives rise to additional local-field screening terms in the expression for the OMP. We shall only consider the case where the OMP is calculated for a reference state at zero field, which we indicate by a superscript "0." Upon inserting (34a) into (1) we obtain the three gauge-invariant OMP terms in (40). The term α CS is clearly of the isotropic form (2), with This is the same expression as obtained previously by heuristic methods [6,7] §. The other two terms were not considered in the previous works. They are where ∂ D denotes the field-derivative ∂/∂E d and are the first-order field-polarized states projected onto the unoccupied manifold. The terms containing ∂ D H 0 k describe the screening by local fields. They vanish for tightbinding models such as the one in this work, but should be included in self-consistent calculations, in the way described in [37]. We shall sometimes refer to α LC and α IC as 'Kubo' contributions because, unlike the Chern-Simons term, they involve firstorder changes in the occupied orbitals and Hamiltonian, in a manner reminiscent of conventional linear-response theory.
Equations (47a)-(47c) are the main result of this section. The derivation of (47b) and (47c) is somewhat laborious and is sketched in Appendix B. We emphasize that the Kubo-like terms, besides endowing the tensor α with off-diagonal elements, also generally contribute to its trace, which therefore is not purely geometric. Writing the isotropic part of the OMP response in the form (2), we then have The terminology 'Kubo terms' for α LC and α IC is only meant to be suggestive. A Kubo-type linearresponse calculation of the OMP should produce all three terms, including α CS .   The two contributions are plotted for our model in figure 3. Moreover, the open circles in figure 1 show the zz and zy components of the OMP tensor computed from (47a)-(47c), confirming that the analytic field derivative of the magnetization was taken correctly.
In the case of an insulator with a single valence band, the partition (40) of the OMP tensor acquires some interesting features. The terms α IC and α CS become purely itinerant, i.e., they vanish in the limit of a crystal composed of non-overlapping molecular units, with one electron per molecule. Also, the first term in the expression (47c) for α IC -the only term for tight-binding models -becomes traceless, as can be readily verified in a Hamiltonian gauge (where H 0 k |u 0 nk = E 0 nk |u 0 nk ) with the help of the perturbation theory formula [33] In order to verify these features numerically, we calculated the various contributions treating only the lowest band of our tight-binding model as occupied. The molecular limit was taken by setting to zero the hoppings between neighbouring eight-site cubic "molecules." It could have been anticipated from the outset that the Chern-Simons term (47a) could not be the entire expression for the OMP, based on the following argument [26]. Consider an insulator with N > 1 valence bands, all of which are isolated from one another. By looking at α da as ∂P d /∂B a one can argue that, since each band carries a certain amount of polarization P (n) , the total OMP should satisfy the relation where α (n) is the OMP one would obtain by filling band n while keeping all other bands empty. We shall refer to this property as the "band-sum-consistency" of the OMP. It only holds exactly for models without charge self-consistency (see the analytic proof in Appendix C), but that suffices for the purpose of the argument. We note that the Chern-Simons contribution (47a) alone is not band-sum-consistent, because the second term therein vanishes for a single occupied band. Hence an additional contribution, also band-sum-inconsistent, must necessarily exist. Indeed, both α LC and α IC are bandsum-inconsistent, in such a way that the total OMP satisfies (51). This is illustrated in figure 4 for our tight-binding model.

Summary and outlook
In summary, we have developed a theoretical framework for calculating the frozen-ion orbital-magnetization response (OMP) to a static electric field. This development fills an important gap in the microscopic theory of the magnetoelectric effect, paving the way to first-principles calculations of the full response. While the OMP is often assumed to be small compared to the lattice-mediated and spin-magnetization parts of the ME response, there is no a priori reason why it should always be so. In fact, in strong Z 2 topological insulators it is the only contribution that survives, and the predicted value is large compared to that of prototypical magnetoelectrics. Although the measurement of the θ = π ME effect in topological insulators is challenging, as time-reversal symmetry must be broken to gap the surfaces [6,7,23], there may be other related materials where those symmetries are broken already in the bulk. The present formalism should be helpful in the ongoing computational search for such materials with a large and robust OMP.
A key result of this work is a k-space expression for the orbital magnetization of a periodic insulator under a finite electric field E (equations (6a)-(6d), or equivalently, (34a)-(34d)). In addition to the terms (34b)-(34c) already present at zero field [13], in three dimensions the field-dependent magnetization comprises an additional purely isotropic 'Chern-Simons' term, given by (34d). This new term depends explicitly on E and only implicitly on H 0 k , while the converse is true for the other terms. Moreover, it is a multivalued quantity, with a quantum of arbitrariness M 0 = Ee 2 /hc along E. Thus, the analogy with the Berry-phase theory of electric polarization [10,20], where a similar quantum arises, becomes even more profound at finite electric field.
The Chern-Simons term M CS is responsible for the geometric part of the OMP response discussed in [6,7] in connection with topological insulators. We have clarified that in materials with broken time-reversal and inversion symmetries in the bulk the CSOMP does not generally constitute the full response, as the remaining orbital magnetization terms, M LC and M IC , can also depend linearly on E. Their contribution to the OMP, given by (47b)-(47c), is twofold: (i) to modify the isotropic coupling strength θ; and (ii) to introduce an anisotropic component of the response.
Another noteworthy result is equation (45) for the Chern-Simons OMP of finite systems. One appealing feature of this expression is that it allows one to calculate the CSOMP without the need to choose a particular gauge. Instead, its k-space counterpart, equation (47a), requires for its numerical evaluation a smoothly varying gauge for the Bloch states across the Brillouin zone. Equation (45) is also the more general of the two, as it can be applied to noncrystalline or otherwise disordered systems.
We conclude by enumerating a few questions that are raised by the present work. Do the individual gauge-invariant OMP terms identified here in a one-electron picture remain meaningful for interacting systems, and can they be separated experimentally? (This appears to be the case for M LC and M IC at E = 0 [35].) How does the expression (47a)-(47c) for the linear OMP response change when the reference state is under a finite electric field E? Finally, we note that equation (41c) for the CSOMP of finite systems has a striking resemblance to a formula given by Kitaev [39] for the 2-D Chern invariant characterizing the integer quantum Hall effect. Can this connection be made more precise, in view of the fact that the quantum of indeterminacy in θ CS is associated with the possibility of changing the Chern invariant of the surface layers? These questions are left for future studies.

Acknowledgments
This work was supported by NSF Grant Nos.
DMR 0706493 and 0549198. Computational resources have been provided by NERSC. The authors gratefully acknowledge illuminating discussions with Andrew Essin, Joel Moore, and Ari Turner.

Tight-binding model
We have chosen for our tests a model of an ordinary (that is, non-topological) insulator. The prerequisites were the following. It should break both time-reversal and inversion symmetries, as the OMP tensor otherwise vanishes identically. It should be threedimensional, as the geometric part of the response vanishes otherwise. Its symmetry should be sufficiently low to render all nine components of the OMP tensor nonzero. Finally, it should have multiple valence bands, for generality.
We opted for a spinless model on a cubic lattice. It can be obtained starting with a one-site simple cubic model, doubling the cell in each direction, and assigning random on-site energies E i and complex first-neighbour hoppings t j→i = te iφ j→i of fixed magnitude t = 1. The Hamiltonian reads Table A1. The parameters of the tight-binding model. Columns I-III give the site coordinates i = (x, y, z), in units of the lattice constant a = 1 of the 2 × 2 × 2 primitive cubic cell. Column IV contains the on-site energies E i , and the last three columns contain the phases of the complex nearest-neighbour hopping amplitudes along bonds in the negativex,ŷ, andẑ directions. where i = (x, y, z) labels the sites and ij denotes pairs of nearest-neighbour sites. The values of E i on two of the eight sites were adjusted to ensure a finite gap everywhere in the Brillouin zone between the two lowest bands (chosen as the valence bands) and the remaining six. We also made sure that nonzero phases φ j→i were not restricted to twodimensional square-lattice planes, otherwise those are mirror symmetry planes, whose existence is sufficient to make the diagonal elements of the OMP tensor vanish. In our calculations all the model parameters were kept fixed except for one phase, which was scanned over the range [0, 2π], and the results are plotted as a function of this phase ϕ. For reference, the on-site energies and the phases of the hopping amplitudes are listed in table A1. The energy bands are shown in figure A1 for ϕ = 0. In order to couple the system to the electric field and to be able to define its orbital magnetization, the position operator r must be specified along with H 0 . We have chosen the simplest representation where r is diagonal in the tight-binding basis.

Technical details
The calculations employing periodic boundary conditions were carried out on an 80 × 80 × 80 k-point mesh, and the k-space implementation of finite electric fields was done using the method discussed in section V of [30]. The open boundary condition calculations used cubic samples containing L × L × L eight-site unit cells, that is, 2L + 1 sites along each edge. For large L, we expect the magnetization to scale as where a, b, and c account for face, edge, and corner corrections, respectively [13].
Calculations of M (L) under small fields were done using L = 4, 5, 6, 7, and then fitted to (A.2) in order to extract the value M of the magnetization in the L → ∞ limit. The  47c)] on a grid, they need to be properly discretized. The presence of the gauge-dependent Berry connection in (6d) demands the use of a "smooth gauge" for its evaluation, where the valence Bloch states given by (18) are smoothly varying functions of k. This is achieved by projecting a set of trial orbitals onto the set of occupied Bloch eigenstates according to the prescription in equations (62-64) of [31]. (For the tightbinding model discussed below, when treating the two lowest bands as occupied, the two trial orbitals are chosen as delta functions located at the two sites with lowest on-site energy.) If needed, this one-shot projection procedure can be improved upon by finding an optimally smooth gauge using methods based on minimizing the real-space spread of the WFs [31], but we found our results to change negligibly when performing this extra step. In a smooth gauge the needed k-derivatives of the Bloch states and of the Berry connection matrix are then evaluated by straightforward numerical differentiation. Note that (6b) and (6c) should be evaluated in the same smooth gauge as (6d), as these three equations are not separately gauge-invariant. A smooth gauge must also be used for (34d) and (47a), because, as discussed in section 2.3.1, the Chern-Simons 3-form is locally gauge-dependent.
The same strategy can be used to discretize (34b) and (34c). However, since the k-derivatives appearing in those equations are covariant, the discretized form of the covariant derivative (19) given in [30,13] may be used instead, circumventing the need to work in a smooth gauge. We have implemented both approaches, finding excellent agreement between them.
Finally we come to equations (47b) and (47c). In addition to the k-derivative of the valence Bloch states, we need their (covariant) field-derivative (48), as well as the k-derivative of H 0 k . The latter quantity is easily calculated within the tight-binding method, and for the former we used the linear-response expression (50). Note that this requires choosing the unperturbed states to be in the Hamiltonian gauge. This choice precludes calculating the k-derivative on the right-hand-side of (50) by straightforward finite differences, which can only be done in a smooth gauge. But because u 0 mk |∂ d u 0 nk equals u 0 mk | ∂ d u 0 nk for m > N, the discretized covariant derivative approach may be used instead. Alternatively, one can evaluate the ordinary k-derivative by summation over states as We note that this formula may not be used to calculate the geometric term (47a), because it induces locally a parallel transport gauge (A 0 nn = 0), which cannot be enforced globally since the Brillouin zone is a closed space.

Appendix B. Derivation of equations (47b) and (47c)
For notational simplicity we drop the crystal momentum index k. So, for example, |u nk shall be denoted by |u n . In order to calculate the OMP terms and starting from (31) and (32), we shall first examine the field-and k-derivatives of certain basic quantities. We begin by noting that the field-derivative ∂ D P = −∂ D Q of the projection operator (20) can be written as in terms of the covariant field-derivative (48) (a similar expression holds for the kderivative). This follows from a relation analogous to (21): is the Berry connection matrix along the parametric direction E d . With the help of (B.3) the field-derivative of (22) becomes where F bD is obtained from F bd by replacing ∂ d with ∂ D . We shall also need the fieldand k-derivatives of the matrix H 0 nm defined by (8): where we introduced the notation (∂ D,c H 0 op ) nm ≡ u n |∂ D,c H 0 |u m , where 'op' indicates that the derivative is taken on the operator itself, not its matrix representation. These two relations follow directly from (B.4) and (21). We will also make use of identities such as Equations (B.14) and (B.13) are respectively equivalent to (47b) and (47c) in the main text. The gauge invariance of these equations follows from the fact that they are written as traces over gauge-covariant objects. (We also note that the covariant derivative transforms according to (18) regardless of the parameter with respect to which the differentiation is carried out.)

Appendix C. Band-sum consistency of the OMP
Here we show analytically that the OMP tensor α satisfies the band-additivity relation (51) in models without charge self-consistency. In order to isolate the contribution α (n) coming from valence band n (assumed to be well-separated in energy from all other bands), we choose the Hamiltonian matrix to be diagonal at zero field, i.e., H 0 mnk (E = 0) = E 0 nk δ mn . If in addition we use a parallel-transport gauge for the linear electric field perturbation [33] (this is achieved by setting to zero the matrix A D defined in (B.5)) we find, using (B.7), ∂ D H 0 mnk | E=0 = 0. With the help of these two relations, the field-derivative ∂ D M a | E=0 of (6a) is easily taken. From the first two terms therein we obtain (dropping the index k) (C.1) In the parallel-transport gauge |∂ D u n is given by (50), and combining the resulting expression with the field-derivative of the third term in (6a) yields The LHS is proportional to the difference between α da and N n α (n) da , and vanishes as a result of an exact cancellation between the terms (n, m) and (m, n) in the double sum. The integrand of the (n, m) term is