Electric polarization and magnetization in metals

A feature of the"modern theory"is that electric polarization is not well-defined in a metallic ground state. A different approach invokes the general existence of a complete set of exponentially localized Wannier functions, with respect to which general definitions of microscopic electronic polarization and magnetization fields, and free charge and current densities are always admitted. These definitions assume no particular initial electronic state of the crystal, and the set of microscopic fields satisfy the usual relations of classical electrodynamics. Notably, when applied to a trivial insulator initially occupying its $T=0$ ground state, the expressions for the unperturbed polarization and orbital magnetization, and for the orbital magnetoelectric polarizability tensor obtained from these different approaches can agree. However, the"modern theory of magnetization"has been extended via thermodynamic arguments to include metals and Chern insulators. We here compare with that generalization and find disagreement; the manner in which the expressions differ elucidates the distinct philosophies of these approaches. Our approach leads to the usual electrical conductivity tensor in the long-wavelength limit; in the absence of any scattering mechanisms, the dc divergence of that tensor is due to the free current density and the finite-frequency generalization of the anomalous Hall contribution arises from a combination of bound and free current densities. As well, in the limit that the electronic ground state is that of a trivial insulator, our expressions reduce to those expected for the unperturbed polarization and magnetization, and the electric susceptibility.


I. INTRODUCTION
In elementary classical electrodynamics, the macroscopic charge and current densities in material media are written in terms of (electric) polarization P (x, t) and magnetization M (x, t) fields, and "free" charge and current densities F (x, t) and J F (x, t), (x, t) = −∇ • P (x, t) + F (x, t), J (x, t) = ∂P (x, t) ∂t + c∇ × M (x, t) + J F (x, t).(1) Going back to the time of Lorentz, P (x, t) and M (x, t) have typically been taken to involve those charges that remain "bound" within individual atoms and molecules, while F (x, t) and J F (x, t) are associated with other charges that are "free" to move through the medium.
In more modern treatments of metallic crystals and doped semiconductors, if the motion of the ion cores is neglected then F (x, t) and J F (x, t) are associated with intraband electronic transitions within partly occupied energy bands.In the "long-wavelength limit," where the wavelength of an applied electric field is much larger than the lattice constant, the response of those carriers is calculated as if that field were uniform.For example, in a too-simplistic model in which scattering is neglected and the relevant carriers are all assumed to have the same effective mass m 0 and carry the electric charge e = −|e|, for a uniform electric field oscillating at frequency ω with amplitude E(ω) the amplitude of the uniform current density driven in linear response is where the superscript (1) indicates the linear response of the quantity, and N is the density of relevant carriers.Turning then to the other terms in the second of (1), in the long-wavelength limit the macroscopic magnetization is uniform and the "bound" current density ∂P (t)/∂t is associated with interband electronic transitions involving occupied or partly occupied energy bands.The simplest procedure, even more elementary than a Kubo approach, is to calculate the interband absorption rate using Fermi's Golden Rule, and associate that absorption with the absorption that would result from a model in which the polarization responded to the electric field through a dielectric tensor δ il + il inter (ω), which would give where superscript indices indicate Cartesian components and are summed over if repeated.This association identifies the imaginary part of il inter (ω), and the real part of il inter (ω) can then be found using the Kramers-Kronig relation [1].Using both (2) (or a less simplistic version) and (3), with il inter (ω) so determined, a calculation of the linear response in the long-wavelength limit is complete.Sometimes one even introduces an "effective" dielectric constant il eff (ω), formally writing the linear response J (1) (t) of the full current density from the second of (1) as just J (1) (t) = ∂P Then, in terms of the calculated il inter (ω) and within the simple model (2) for the intraband response, we have This strategy is somewhat indirect.One might suppose that the polarization would be defined, and then its response to the electric field calculated.But such a definition is bypassed by calculating il inter (ω), that is, the contribution that a purported polarization would make to the optically induced current density.And there is no definition of either a purported polarization or magnetization that would exist before the electric field is applied.
Of course, the use of an approach that bypasses such definitions is not surprising.Any consideration of the response of "bound" charges and their currents, and the polarization and magnetization to be associated with them, is at least initially suspect from the perspective of the quantum theory of solids.In fact, problems in defining a polarization and magnetization arise even for the ground state of a crystal [2].In recent years the "modern theories of polarization and magnetization" have been developed to clarify these concepts, primarily focused on insulators [3][4][5], and have provided many physical insights, including the "quantum of ambiguity" inherent to the unperturbed polarization of the T = 0 ground state [6], the existence of two distinct contributions to the orbital magnetization [4,5], and that a static and uniform magnetic (electric) field can induce a polarization (magnetization) [7,8].However, the "modern theories" are based on static or adiabatically varying uniform fields, and are not immediately applicable to treat the optical properties of materials, especially at wavelengths so small -beyond the "long-wavelength limit" -that one has to take into account the variation of the optical fields over a unit cell.
As well, the main focus of the "modern theories" has been "topologically trivial" insulators, a class of band insulators that we define below.Indeed, among the contributions of the "modern theory of polarization" is that there may be a relationship between a certain "localization" of the electronic ground state and the polarization of an unperturbed crystal [9][10][11], and it has been argued that a "localized" ground state is necessary for the polarization to be well-defined; the ground state of a metallic crystal is found to violate that condition.In contrast, the "modern theory of magnetization" has been extended to include metals and Chern insulators [5,12].These extensions are based on thermodynamic arguments, and thus again are not applicable to optical fields.Indeed, there seems no straightforward roadmap for extending the approach of the "modern theories" to frequency dependent polarizations, magnetizations, and free currents.
A set {P (x, t), M (x, t), F (x, t), J F (x, t)} that satisfies (1) is far from unique.An underlying pillar of the "modern theories" is that, in finite-sized media, P and M , when taken as the usual charge and current density dipole moments, are experimentally accessible; it is implicitly assumed that the numerical value of the bulk quantities should coincide with those [4,5,13].In this way, strange properties of the bulk expressions, for example, that the ground state magnetization in a Chern insulator involves a chemical potential or that polarization is not a well-defined bulk quantity in a metal while magnetization is, are justified.However, the relation between these quantities in bulk and finite-sized systems is not straightforward and considerations at the boundary are often important, even in insulators; for example, the bulk topological magnetoelectric coefficient does not generically determine that of a thin film [14,15].
In recent work [16] we have taken a different approach, which is related more directly to the classical strategy of Lorentz, and our focus has been on bulk crystals with static ions.Here, polarization and magnetization fields, and free charge and current densities, serve as intermediary quantities that aid calculation and provide physical insight, but in general only the appropriate combinations that lead to the charge and current densities have direct physical significance [17].To identify the electronic component of these quantities, we employ a complete set of exponentially localized Wannier functions (ELWFs) [18] (or "modified" versions thereof) with respect to which we decompose the charge and current density expectation values [19] as sums of spatially-localized contributions, one associated with each lattice site R [20].From these we define microscopic "site" polarization p el R (x, t) and magnetization mR (x, t) fields in a manner similar to that of atomic and molecular physics.Since charge continuity does not generally hold site-wise [21], a "corrective" contribution mR (x, t) to mR (x, t) arises, and in all m R (x, t) = mR (x, t) + mR (x, t).Natural definitions for site charges and currents that link the lattice sites emerge as well, which are used to define microscopic free "site" charge and current densities.After identifying the ionic contribution to those site quantities, we take their lattice sums to be the microscopic polarization and magnetization fields, and free charge and current densities.
For an unperturbed "trivial" insulator occupying its T = 0 ground state, the spatial integral of mR (x) ( mR (x)) coincides with that site's "atomic-like" ("itinerant") contribution to M of the modern theory [5].Here M is unique (i.e.not "gauge dependent" in that it does not depend on the choice of smooth frame of the bundle of occupied electronic Hilbert spaces over the first Brillouin zone (BZ), or equivalently on how the ELWFs are chosen, assuming they are taken "occupied.");from this perspective, this is but a special case.In contrast, the spatial integral of p el R (x), which coincides with that site's contribution to P el of the modern theory [6], is gauge dependent and only unique modulo a "quantum of ambiguity."And in the optical response of such an insulator, wherein the linearly induced charge and current densities arise entirely from induced electric and magnetic multipole moments [16], it is only the combinations of such moments corresponding to those densities that are generally gauge invariant and of direct physical significance [22]; we there find agreement with the usual approach involving a q-expansion of the conductivity tensor [23].
In this paper we implement this approach to treat the optical response of a metal.In this initial treatment we restrict ourselves to the long-wavelength limit and consider the independent particle approximation, in which the interaction between electrons is approximately treated through an effective potential energy characterizing the lattice and by taking the "applied" electric field in our calculations to be the macroscopic Maxwell electric field, having frequency components E(ω).
Our identification of the electronic component of the "site quantities" requires the existence of a complete set of ELWFs, which depends entirely on topological considerations [24].In Sec.II we briefly discuss such issues, but the result is that if the Bloch energy eigenfunctions associated with all of the energy bands are employed, a complete set of ELWFs can always be constructed.Thus, the microscopic polarization and magnetization fields we introduce, and the corresponding macroscopic fields, are always defined, as are the macroscopic free charge and current densities.In this paper we will restrict our study to those crystalline solids for which ELWFs can be constructed from the energy eigenfunctions associated with any set of isolated bands -including the completely occupied and partly occupied bands in a p-doped semiconductor in its T = 0 ground state, which is the model of a metal we adopt in this first communication.
In Sec.II we also introduce the basic equations of our approach, relying heavily on earlier work [16].In Sec.III we calculate the polarization and magnetization in the unperturbed ground state, and discuss their form; in Sec.IV we calculate the linear response in the longwavelength limit.If the crystal is assumed to initially possess time-reversal symmetry and its energy bands are isolated, then our results follow the pattern sketched in the first three paragraphs of this section: The induced free current can be associated with intraband transitions, and the polarization current with interband transitions.Here, of course, we have an explicit expression for the polarization and can calculate the polarization current by directly taking ∂P (t)/∂t; we can thus construct an expression for ij eff (ω) by direct calculation, which is gaugeinvariant as expected.
The situation is more complicated in the absence of time-reversal symmetry.There we find that the first order response of each frequency component of the microscopic charge density, ρ(x, ω) (1) , contains a term proportional to ω −1 and thus diverges as ω → 0. It is then not surprising that the contribution associated with each lattice site R, ρ R (x, ω), also diverge as ω → 0, and thus that P (1) (ω) -associated with the electric dipole moments of those localized charge densities -does as well.Such a result is inevitable in the approach we adopt, where polarization and magnetization are associated with quantities localized about each lattice site.This divergent term originates from an "intraband contribution" to P (1) (ω), and leads to a finite contribution to the induced macroscopic current density −iωP (1) (ω) as ω → 0. In addition to a contribution to J (1) F (ω) that is divergent as ω → 0, which arises as an expected generalization of (2), we also find a contribution that is finite as ω → 0. When this is combined with the contribution to −iωP (1) (ω) that is finite as ω → 0 we find a gauge-invariant contribution to J (1) (ω) that is finite as ω → 0, and can be identified as giving rise to the anomalous Hall current.The other contributions to the full J (1) (ω), which are also gauge-invariant, correspond to a generalization of (2) and to a pure interband response of P (1) (ω).Here we can also introduce an il eff (ω), but it is perhaps more natural to write where σ il (ω) = −iω( il eff (ω) − δ il )/(4π), and at the end of Sec.IV we give the general expression for σ il (ω).
A reader might ask, "why bother?"Expressions for σ il (ω) can always be derived using Kubo's approach [25], and there is an intrinsic ambiguity in how polarization and magnetization fields are defined.And isn't the idea of a polarization -and certainly an induced polarization -in a metal suspect if it goes beyond merely the "formal" role played by, for example, the effective polarization (4)?
One reason is that usual calculations made in minimal coupling can require the identification of sum rules to show properly behaved results at low frequencies [26][27][28][29] -especially if nonlinear optical response is calculated, which is a future direction for this work -and that is not a difficulty with the calculations presented here, since the response is calculated as due to electric and magnetic fields directly.A second reason is that in an insulator there is a clear physical significance to the response of the polarization to applied fields, as has been demonstrated within the "modern theory" [6], and we feel it is interesting to see how that response can be seen to follow from the response of a p-doped semiconductor in the limit of vanishing doping.After all, since one can move from metal-like behavior to insulator-like behavior in this limit, it would seem physically reasonable to expect a polarization that would continuously evolve from that of a metal to that of an insulator.A third reason is that with this approach we can establish a connection to an earlier generation of calculations based on strategies introduced by Blount and co-workers [30].A fourth reason, we feel, is the interesting way a calculation based on polarization and free currents highlights the way broken time-reversal symmetry leads to a response qualitatively different than usual.And a fifth reason is that, with its emphasis on ELWFs and the interest in those functions for electronic structure and response calculations in general, we can hope that the approach here will be useful in numerical calculations.
Our conclusions and perspectives on future work are presented in Sec.V. Ultimately, when implemented in a metallic crystal, that our general definitions agree with past work of Blount and co-workers in a simple limit, and that the well-known σ il (ω) results, provides positive support for our approach.

II. SINGLE-PARTICLE DENSITY MATRIX
We consider a simple instance of a metal, a p-doped semiconductor, perturbed by a uniform electric field.We restrict our study to bulk crystalline solids of spatial dimension two or three (d = 2, 3), and implement the frozen-ion and independent-particle approximations; the spin degree of freedom is neglected.Thus, in the Heisenberg picture the only dynamical degree of freedom of the crystal is the electron field operator.In an unperturbed crystal we denote that field by ψ0 (x, t) with dynamics governed by the equation of motion i d dt ψ0 (x, t) = [ ψ0 (x, t), Ĥ0 ] for one-body Hamiltonian operator Ĥ0 = ψ † 0 (x, t)H 0 (x, p(x)) ψ0 (x, t)dx, where with V (x) = V (x + R) for any Bravais lattice vector R characterizing H 0 x, p(x) , and with to allow for the presence of an "internal," static, cellperiodic magnetic field described by the vector potential A static (x), where A static (x) = A static (x+R), that generally breaks time-reversal symmetry.We assume that the set of energy eigenvalues E nk of the cell-periodic Hamiltonian (6) admits a band gap, below which we take the Fermi energy E F to lie, unless otherwise stated.Thus, a distinction can be made between the set of Bloch energy eigenvectors that are associated with partly occupied energy bands and those that are associated with completely unoccupied energy bands in the T = 0 ground state, which we take to be the initial state of the crystal.
In recent times, it has become clear that the spectral data of the relevant Hamiltonian does not entirely characterize a bulk crystal and that some topological data must also be identified.In particular, the existence of a complete set of ELWFs [32] is equivalent to the existence of a global smooth frame of the Hilbert bundle over BZ with fibres constructed point-wise for each k ∈ BZ as the linear span of the cell-periodic parts u nk (x) ≡ x|nk of the Bloch energy eigenfunctions ψ nk (x) ≡ x|ψ nk associated with all of the energy bands, which we term the Bloch bundle [33].In fact, such a frame always exists [34] and its components |αk -which, for each k ∈ BZ, constitute an orthonormal basis of the fibre of the Bloch bundle at that k -can generally be written [24,35,36] where the U nα (k) constitute a unitary matrix U (k) at each k; in what follows, sums are generally taken over all band indices n or all "type" indices α unless otherwise indicated.It is then each of the |αk , which are smooth over the BZ, that can be mapped to an ELWF W αR (x) ≡ x|αR via the (inverse) Bloch-Floquet-Zak transform [37,38], where Ω uc is the volume of the real space unit cell; each ELWF is identified by a type index α and the Bravais lattice vector R with which it is associated.Additionally, the existence of such a global smooth frame of the Bloch bundle is equivalent to the existence of a global trivialization thereof, thus any (collection of) Chern number(s) characterizing it vanish [39]; this is often understood implicitly in the physics literature [40].The construction outlined above corresponds to using the eigenvectors associated with all of the energy bands to construct a complete set of ELWFs; in general, each U nα (k) is nonvanishing.However, often times a number of Hilbert subbundles (of the Bloch bundle) that are associated with sets of isolated energy bands are trivial, in which case the corresponding subsets of energy eigenfunctions can be separately used to construct subsets of a complete set of ELWFs.
We here restrict our study to crystals for which the Hilbert bundle associated with any set of isolated energy bands is globally trivial.Taking there to be an energy gap between bands N and N + 1, U (k) can always be taken of block diagonal form with the "upper left" block being N ×N dimensional.If the Fermi energy lies in that gap -for example, if E F coincides with the upper (red) dashed line of Fig. 1 -then we will classify the material as a "trivial" insulator, and at each k the U (k) acts on the occupied and unoccupied states |nk separately.On the other hand, if the Fermi energy lies below that gapfor example, if E F coincides with the lower (blue) dashed line of Fig. 1 -then at each k the U (k) acts on the states |nk associated with the N partly occupied energy bands separate from the remaining unoccupied states.
Such considerations generally apply to any crystal whose electronic spectrum has a band gap.A more general approach to generate such smooth frames in metallic crystals where it is not necessary to have isolated sets of energy bands has been formulated [41]; in future work we plan to implement this construction.However, even within the simplified scheme that we implement, an important distinction between metals and trivial insulators arises: For a trivial insulator with N occupied energy bands there exists a global smooth frame of the occupied subbundle with components |αk labelled by integers α ∈ {1, 2, . . ., N } that satisfies ∀k ∈ BZ : span C ({|nk |E nk < E F }) = span C ({|αk |α ∈ {1, 2, . . ., N }}), while for metals with N partly occupied bands there exists only such a frame that satisfies ∀k ∈ BZ : span C ({|nk |E nk < E F }) ⊆ span C ({|αk |α ∈ {1, 2, . . ., N }}) [35,41].This is because for metallic systems the construction of a vector bundle over BZ whose fibre at each k ∈ BZ is the Hilbert space spanned by the occupied |nk fails [42].Instead, one can construct a vector bundle whose fibre at each k contains, as a subspace, the occupied

A. FLAPW method
In order to determine the optical response functions in a full band structure approach, one requires the eigenvalues and velocity matrix elements at many k points in the Brillouin zone ~BZ!.The velocity matrix elements, in turn, require a knowledge of the electronic wave functions.For this purpose, we employ a first-principles approach in the form of the FLAPW method.As this method has been previously discussed, 31,32 we highlight only a few of its pertinent features.
The spirit of the FLAPW approach is the partitioning of FIG. 1. FLAPW electronic band structures for GaAs and GaP.The fundamental band gap has been adjusted within the scissors approximation.
1: Schematic of the energy bands whose associated eigenvectors would be used in the construction of ELWFs in a hypothetical d > 1 crystalline solid (the bandstructure of GaAs, which we import from a past publication [31], is used only for illustrative purposes).Upper (red) and lower (blue) horizontal dashed lines indicate possible Fermi energies for a trivial insulator and for a p-doped semiconductor (that we here take as a simple instance of a metal), respectively.
Hilbert space at that k, which yields the subset relation.Thus, while for trivial insulators the subspace span C ({|αk |α ∈ {1, 2, . . ., N }}) contains only "ground state data," this is not so for metals.This does not pose an issue since we do not assert that the electronic polarization and magnetization fields involve only the initially occupied energy eigenvectors.Rather, we introduce p el (x, t) and m(x, t) using any set of functions that are sufficiently localized spatially -ELWFs are the most natural and convenient choice -and, by construction, from those fields the ground state expectation values of the charge and current density operators can be found [16].Since topological notions underlie the existence of EL-WFs, the appearance of related geometric objects in many identities involving ELWFs is less opaque than it might otherwise be.One such identity that will be useful in this work is [35] where, denoting the inner product on the Hilbert space spanned by the set of cell-periodic functions u nk (x) (which are here taken normalized over the real-space unit cell are components of the non-Abelian Berry connection that is induced by a global smooth frame with components |αk .Here u αk (x) ≡ x|αk and we adopt the shorthand ∂ a ≡ ∂/∂k a .The components (11) are related to the components of the non-Abelian Berry connection that is induced by a local smooth frame with components |nk , via the gauge transformation In a periodic gauge choice (8), which we always employ, where G is a reciprocal lattice vector, all the objects appearing here, including the Hermitian matrix W a (k) [2] populated by elements are periodic over BZ.In what follows, the k-dependence of the preceding objects is usually kept implicit.
The matrix elements of W a (k) that are nonvanishing depend on the structure of U (k), which for the materials we consider can take the block-diagonal form discussed above.Then W a mn (k) = 0 only if m and n lie in the same block, for if they are associated with different blocks then the values of α for which U nα (k) = 0 differ from the values of α for which U mα (k) = 0.In a trivial insulator the first, "upper left" block acts only on the occupied |nk and the second, "bottom right" block only on the unoccupied |nk .If we introduce Fermi filling factors f nk = 1 (0) if the state |nk is initially occupied (unoccupied), then for the trivial insulator f nk = f n , depending only on the band, and W a mn (k) = 0 only if f m = f n .But in a p-doped semiconductor, which is our simple instantiation of a metal in this paper, the first block also acts on some unoccupied |nk , and so in general we can have We account for the interaction between the electron field and the "applied" electromagnetic field via the usual minimal coupling prescription.From the resulting minimal coupling Hamiltonian, the field-theoretic charge and current density operators constituting the Noether current can be found in the usual way [43].Associated with each spatial component of this current density operator is the differential operator where m is the electron mass, and where the vector and scalar potentials A(x, t) and φ(x, t) describe the classical applied electromagnetic field.As a consequence, another useful identity will be where the matrix elements are found to be [44] Under the frozen-ion approximation implemented here, we take the positively charged ion cores that compose the underlying crystal structure of the material to be fixed, even in the presence of an applied electromagnetic field, and introduce the charge density ρ ion (x) to describe the periodic distribution of these static charges.We take ρ ion (x) such that the crystal as a whole is electrically neutral.
In this work, we implement a previously developed formalism [16] restricted to the "long-wavelength limit," wherein we take the applied electric field to be uniform and the magnetic field to vanish.To simplify these initial considerations we here neglect local field corrections, which can be important [45], and take the applied electric field be the macroscopic Maxwell field denoted E(t) [22].Consideration of phenomena related to spatially-varying electromagnetic fields is left for future work.A quantity central to this formalism is the so-called (electronic) single-particle density matrix η αR ;βR (t), the definition of which (see Eq. (15,27,30,33,36) of Mahon et al. [16]) involves a generalized Peierls phase, the fermionic operators that generate ELWFs W αR (x) and "modified" versions WαR (x, t) thereof, and the initial electronic state of the unperturbed crystal, here taken to be the T = 0 ground state.The operators ânk and â † nk generating the eigenvectors |ψ nk of the unperturbed Hamiltonian Ĥ0 are also relevant in perturbative calculations and we take |ψ nk ≡ â † nk |vac , where Ĥ0 |ψ nk = E nk |ψ nk .Still in the Heisenberg picture, the electronic field operators ψ(x, t) here evolve as In what follows, we account for the effect of the applied electric field perturbatively.Thus we assume the existence of a valid expansion of all electronic quantities in powers of E(t).In particular, we take η αR ;βR (t) of the form where the superscript (0) denotes the contribution to a quantity that is independent of E(t), the superscript (1) denotes the contribution that is linear in E(t), and ". .." denotes non-linear contributions, which we here neglect.In Appendices A and B we find and implementing the usual Fourier series analysis, we also find that the first-order perturbative modification to η αR ;βR (t) due to E(ω) is where f nm,k ≡ f nk − f mk .The first term of (20) is the straight-forward generalization of the previously found perturbative modification for trivial insulators, and can be understood in the context of time-dependent pertur-bation theory as arising from the interaction term [46] −eE a (t) Due to the form of this interaction term we will later describe any first-order modifications that involve the first term of (20) as being "interband."The second term of (20) is a new contribution, here related to the presence of a Fermi surface.Notably this term diverges in the dc limit, and indeed it is this term that will lead to the expected dc divergence of the induced free current density, as we later show.This term can here be understood as arising from the interaction term The first contribution (that involving ξ a nn ) to this interaction term gives a vanishing contribution to η (1) αR ;βR (ω) for both metals and trivial insulators initially occupying their T = 0 electronic ground state (see Appendix B).In contrast, although the second and third contributions as well give vanishing contributions to η (1) αR ;βR (ω) for trivial insulators, they give rise to finite contributions if the crystal is metallic; these finite contributions involve only those occupied |nk with energies "near" the Fermi energy.Due to the form of this interaction term we will later describe the first-order modifications that involve the second term of (20) as being "intraband."Such an identification of interaction terms that give rise to interand intraband contributions at linear response is implicit in the earlier works of Blount [47] and others [48].The primary difference between our approach and those is that this investigation is a limiting case of a more general framework within which spatial and temporal variation of electric and magnetic fields can be taken into account; that is not the case in earlier works.
The limit of a trivial insulator can be reached from (18) by taking f nk → f n and requiring each |αk to be an element of either the initially occupied or unoccupied Hilbert subspace; the second condition implies that, in general, W i nm (k) = 0 only if f n = f m (see discussion below ( 14)).In this limit, one can define an analogous filling factor f α associated with |αk ; we set the f α associated with |αk to equal the f n associated with the set of |nk used in its construction.The sum over n in (18) then corresponds to the matrix multiplication of U (k) and its inverse giving the unit matrix at each k, which in components is δ αβ .It then follows that, in this limit, as expected [16].Implementing this limit in (20) we find again, as expected [16].

III. DIPOLE MOMENTS
In general, the presence of an applied electromagnetic field will break the discrete translational symmetry of the unperturbed crystal.Consequently, in the minimallycoupled system, contributions to a given electric or magnetic multipole moment that are associated with distinct lattice sites of the crystal will generally differ.However, in the long-wavelength limit, all of the site contributions to a given multipole moment are equivalent.In particular, in this limit the site electric and magnetic dipole moments satisfy for any Bravais lattice vectors R and R of H 0 x, p(x) .It follows that the macroscopic polarization and magnetization fields are uniform [22] and can be written as for any such R.Because we consider the ionic cores within the crystal to be fixed, these charges do not contribute to the magnetization; there will however be a static contribution to the polarization that is found from the "site" polarization fields that are defined from the constituents of a decomposition of ρ ion (x) into "site" contributions [16,22].Irrespective of the initial electronic state of the crystal, the electric dipole moment associated with lattice site R is defined to be (see Eq. (42,44,55,57) of Mahon et al. [16]) and the magnetic dipole moment associated with R to be (see also Eq. (64, 66, 67) of Mahon et al. [16]), where ρ βR ;αR (x, R; t), j βR ;αR (x, R; t), and jβR ;αR (x, R; t) are termed generalized (electronic) "site-quantity matrix elements" and were introduced previously [16].Again, with the assumption of a valid perturbative expansion (22,23) can be written and the same can be done for (21).The first term in each such expansion is identified as the "unperturbed contribution," and these are the focus of the rest of this section.
Although the matrix elements appearing in (22,23) are generally of quantities arising from a minimally coupled Hamiltonian and written in a basis of "modified" Wannier functions [16], we explicitly show in the following subsections (and in Sec.IV) that, as usual, terms appearing at each order in a perturbative expansion can be written in terms of energy eigenvectors (or equivalently in terms of ELWFs in a crystalline solid) of the unperturbed system.When practical, we include the expression for a quantity written as a product of a BZ integral and an integral involving ELWFs, and as a single BZ integral.Although both forms are equivalent, the numerical implementation of the expressions might favor a particular form.For example, when written as a single BZ integral, some quantities involve diagonal matrix elements of the Berry connection ξ a nn (k), which are typically not easy to evaluate numerically.For such quantities, evaluating integrals involving ELWFs may be more tractable.

A. Electric polarization
From (22), the electronic contribution µ and implementing (10,13,18) we find µ el(0) R to be independent of R, as expected.The resulting unperturbed polarization is which is formally similar to that of a trivial insulator [6,44], for which f nk → f n .As is the case there, apart from a gauge dependent contribution, P (0) vanishes if the unperturbed Hamiltonian is inversion symmetric; as discussed previously [44], we take the gauge dependence of the electronic quantities to be contained entirely within the U (k) and consider any terms that involve this object, including the W i (k), to be "gauge dependent."While it appears that the gauge dependent term appearing in (24) no longer generally evaluates to an element of a set of discrete values, at least not following from the same argument that is presented for trivial insulators [6], P (0) maintains the physically sensible characteristic that upon shifting the origin of all ELWFs by a constant Bravais lattice vector R s , the polarization is altered by an additive constant that is proportional to R s .That is, although there is no longer a "quantum of ambiguity" associated with P (0) for a general change in U nα (k), as occurs for a trivial insulator [6], taking |αR → |αR + R s , or equivalently U nα (k) → e −ik•Rs U nα (k) and thus W a nm (k) → W a nm (k) + δ nm R a s , yields where n f nk is the number of electrons per unit volume.Thus, with respect to simple shifts in the positions of the Wannier function a discrete ambiguty does arise.We do note however that an expression formally similar to (24) arises in the case of a Chern insulator [13], and while in that case the gauge dependent contribution again would not be discretely valued by the original argument of Resta [6], it indeed has this property when treated carefully.We here consider a more general notion of gauge dependence than in the "modern theories" and generalizations thereof, therefore those results are not directly applicable.Nevertheless it may still be the case that the gauge dependent contribution to (24) always evaluates to an element of a set of discrete values and we postpone such an investigation for a later work.Finally, (24) is manifestly invariant under a translation of the energy zero, as one would expect.

B. Orbital magnetization
We first identify the "atomic-like" contribution [4,16] to the unperturbed magnetization, which arises from the term involving j βR ;αR (x, R; t) in (23).We find The "itinerant contribution" [4,16], which arises from the term involving jβR ;αR (x, R; t) in (23), is found to be In the trivial insulator limit described in Sec.II, ( 25) and ( 26) separately reduce to the usual expressions; taking f nk → f n and W a nm (k) = 0 only if f n = f m , and using i∂u nk (x)/∂k a = m ξ a mn (k)u mk (x) to find (∂ , the expressions for M LC and M IC of the "modern theory of magnetization" in a trivial insulator [5] are recovered from ( 25) and (26), respectively.More generally, combining (25) and (26) we have where the second equality follows under the assumption that the integrand is sufficiently well-behaved such that all surface terms can be taken to vanish upon an integration by parts; at first order in the perturbative analysis we also employ such an assumption.Of course, in the trivial insulator limit the usual expression [5] is again recovered; that is, in this limit ( 27) reduces to In particular, in that limit ( 27) is gauge invariant.Moreover, in Appendix C we show that (27) generally vanishes if the unperturbed Hamiltonian is time-reversal symmetric, as expected.
If we again consider the effect of shifting the origin of each ELWF by a Bravais lattice vector R s , we find where we have used p a nn (k) = m ∂ a E nk from (17).The term involving R s vanishes as the net current that flows in an unperturbed crystal occupying its T = 0 ground state is zero; that is, Thus, ( 27) is unaffected by shifting ELWFs, as physically expected.Moreover, it is manifest that ( 27) is unchanged by a translation of the energy zero.

IV. FIRST ORDER MODIFICATIONS
We here consider the linearly induced macroscopic charge and current densities, which can be understood to arise from the induced macroscopic polarization and the induced free charge and current densities; in the longwavelength limit considered here, the induced macroscopic magnetization would be uniform [22] and thus not contribute to (1).Under the frozen-ion approximation that we implement, there are only electronic contributions to such quantities.

A. Electric polarization
We first consider the contribution to (22) that is first order in E(ω).Making contact with past work [22,44], we mention that in general the ρ βR ;αR (x, R; ω) do not involve the electric field and so the only contribution to µ (1) R (ω) is "dynamical," arising from the modification of the single-particle density matrix due to E(ω).Implementing (20) we find This expression has two notable features; it is gauge dependent, and it diverges in the dc limit.The gauge dependence is not troubling because induced free charges and currents are also involved here; ultimately it is only the net induced charge and current densities that need be gauge invariant.Also, in the limit of a trivial insulator the second term (that involving ∂ l (ξ i nn + W i nn )) vanishes and the expected gauge invariant result is recovered [49].It is notable however that the distinct terms of ( 28) are sensitive to different aspects of the gauge transformation; the first term, the interband term, involves only offdiagonal elements of W i (k), while the second term, the intraband term, involves only diagonal elements.This is to be expected because of the way in which the Lie algebra components of the Berry connection appear.Second, it is notable that a diverging linearly induced polarization in the dc limit is not unprecedented.For example, if one considers a hydrogen atom initially occupying its 2s state, dc divergences occur as a result of non-vanishing matrix elements between 2s and 2p states facilitated by an electric dipole interaction term.Such a divergence could arise from the first term of ( 28), but does not occur here as we take the crystal to initially occupy its unique electronic ground state, in contrast to this example for the hydrogen atom.So although such a divergence is not entirely novel in principle, the mechanism underlying the divergence of ( 28) is distinct from that of atomic and molecular physics.We return to this issue in Sec.V.

B. Macroscopic bound and free currents
Like the macroscopic polarization and magnetization, the spatial uniformity of the electric field renders the macroscopic bound and free current densities uniform [22].Thus we do not indicate any spatial dependence of such quantities.Moreover, both the macroscopic bound and free current densities can found from any one of the "site quantities" used in their construction [16].
Implementing (28) we find the linearly induced macroscopic bound current density [16,22], This is non-diverging in the ω → 0 limit, as would be expected physically.Furthermore, in Appendix C we show that (29) vanishes in the ω → 0 limit if the unperturbed Hamiltonian is time-reversal symmetric.
We now consider the linearly induced macroscopic free current density.The corresponding microscopic density is defined as [16] From the definitions presented in that past work, we find the first-order modification to the link currents I(R, R ; ω) to be of the form αR;λR (ω)η The first term of ( 31) is termed a "compositional" modification, arising due to a dependence of the generalized site quantity matrix elements on the electromagnetic field, and the second a "dynamical" modification.An expression for ( 31) is given in Appendix D, which can explicitly be shown to satisfy as required, as well as as one would physically expect for a translationally invariant system subject to a uniform electric field.The latter can be understood by noting that the electronic "site charges" evolve according to [16] dQ and in this case we expect there to be no build up of charge at any particular lattice site; we therefore expect dQ R (t)/dt to vanish.In fact, from the definition of Q R (t) and using (20) it can be show that αR;αR (ω) = 0.
Then, in Appendix D we show where the first term in the second equality results from the compositional modification of (31), while the second and third terms from the dynamical modification.Notably J F (ω) diverges in the dc limit, which is as one would physically expect given that we have not accounted for any scattering mechanisms.In fact, it is the third term in the second equality of (32) that will lead to the dc divergence of the electrical conductivity tensor; this term involves the second term of (20).Moreover, J F (ω) is gauge dependent, akin to J (1) B (ω), and it is only through this gauge dependence that (32) involves "interband" contributions.Notably if all of the energy bands of the unperturbed crystal were isolated from one another and the corresponding Hilbert bundles assumed trivial, then one can take U nα (k) proportional to δ nα such that W a nm (k) is proportional to δ nm and thus the interband contributions to J (1) F (ω) vanish; in this limit-ing case the induced free current density involves only "intraband" contributions, which is as one would expect for simple models.

C. Time-reversal symmetry
The general expression (28) that we derive for P (1) (ω) has the feature that it contains an intraband contribution and this contribution diverges in the dc limit.From the simple picture of polarization presented in Sec.I, the presence of such a contribution is unexpected.While in general our description thus asserts that this simple picture is not complete, in Appendix C we show that such an intraband contribution vanishes if the unperturbed Hamiltonian is time-reversal symmetric and P (1) (ω) takes the more expected form Here T = will denote an equality that holds in the presence of time-reversal symmetry.Adopting the approach of (3), we find il inter (ω) which, apart from the gauge dependence, is consistent with the insight from analogies with molecular response and the more simple approaches mentioned in Sec.I.That is, for crystalline solids in which time-reversal symmetry holds, it is only interband contributions that are involved in P (1) (ω).However, even in this simple case P (1) (ω) and il inter (ω) remain gauge dependent, and thus the introduction of ELWFs and the ambiguity in their choice need be involved in any discussion of such quantities.In this case the induced macroscopic free current density (32) reduces to still having both interband and intraband contributions.Notably the interband contribution is gauge dependent and cancels with the gauge dependent term appearing in the induced bound current density −iωP (1) (ω) and thus the net induced current density is gauge independent, as one would expect.In the special case of isolated bands, U nα (k) can be chosen proportional to δ nα and the gauge dependent contributions to P (1) (ω) and J (1) F (ω) separately vanish.
In addition, if the"parabolic band approximation" is implemented, that is, if one takes each energy eigenvalue of an occupied state to be E nk = 2 |k| 2 /2m, J F (ω) agrees with (2) after the identification m 0 = m and N = N el .In fact, we find under the parabolic band approximation that il eff (ω) In moving from the first to the second equality in this expression for σ il (ω), relations that hold only in the presence of time-reversal symmetry are implemented.However, we will find that this latter form of σ il (ω) holds even in the absence of time-reversal symmetry.Moreover, in the absence of that symmetry P (1) (ω) takes the more complicated form (28), which generally involves intraband contributions.This results in P (1) (ω) having a more general gauge dependence and as well J F (ω) having a more general gauge dependence.However, as was the case here, when these more general expressions are combined, for instance when constructing il eff (ω) or σ il (ω), the gauge dependent terms again cancel.

D. Induced macroscopic current density
Returning to the more general investigation, and thus allowing the possible breaking of time-reversal symmetry, we again find that although (29) and (32) are not individually gauge invariant and thus are not themselves directly physically observable, their sum is.Indeed, combining ( 29) and ( 32) we find The first term comes from taking −iω = (E nk − E mk ) + (E mk −E nk −iω) in the interband contribution of P (1) (ω) to J (1) B (ω) (( 28) to ( 29)).Then it is only the term that is gauge invariant and explicitly energy dependent in this particular contribution to J (1) B (ω) that is not cancelled when combined with J (1) F (ω), (32).In particular, it is part of the first term in the second equality of (32) (all but that involving f nm,k Im[ξ l nm W i mn ]), which is a compositional modification, that combines with the second term of the contribution of ( 28) to ( 29) when we calculate J (1) (ω), and ultimately it is the combination of these terms that cancel with the interband contribution of ( 28) to ( 29) that is gauge invariant and does not explicitly depend on energy.The remaining gauge dependent terms all involve products of the form f nm,k ξ a nm W b mn and cancel one another.The second term arises from the induced free current density alone and is the only term in (32) that does not cancel with terms from (29).While the "origin" of each of the terms can most easily be seen in the above form of the expression, it can be rewritten in the more familiar form This is in agreement with usual perturbative calculations that implement the minimal coupling Hamiltonian.In particular, using (17) to rewrite the integrands of (33) to involve velocity matrix elements v nn (k) = p nn (k)/m, for example, Eq. (25,26) of Allen [50] are reproduced.The final term of ( 33) can be understood as a "Drude" contribution.This term follows from the final term of (20), and enters here via the induced free current density (32).Notably, such a term can lead to an induced current density that is orthogonal to the applied electric field.This is not to be confused with the well-understood anomalous Hall conductivity however, because in this case since the Cartesian components i and l are symmetric there exists a basis in which this contribution to the conductivity tensor is diagonal.Physically this means that, were the applied electric field characterized by a single non-vanishing component with respect to such a basis, the induced current density arising from this term would be parallel to that field.Thus, we understand the possibility of such an induced orthogonal current density to be entirely a consequence of crystalline anisotropy.
In contrast, the first and second terms of ( 33) are related to both the induced bound and free current densi-ties.Notably, the second term can be understood as a finite-frequency generalization of the "anomalous Hall" current density [51].This portion of the induced current density is unique because, unlike the contribution from final term of (33), the spatial components i and l are asymmetric and consequently there does not exist a basis in which this contribution is diagonal; there does not exist a basis in which the induced current associated with this term is parallel to the applied electric field.

E. Microscopic charge and current densities
The divergence of (28) in the dc limit may raise concerns about our identification of the polarization.We are thus motivated to consider the first-order modifications of the expectation values of the electronic charge and current density operators due to E(ω) -quantities that could be found from traditional perturbation theory with the minimal coupling Hamiltonian [52] -with the hope that further insight might be gained.Implementing (20) into previously developed expressions [16], we find ρ(x, ω) and ĵi (x, ω) The electronic charge and current density operators that we implement are those that arise via Noether's theorem and thus satisfy the continuity equation Assuming an expansion of these operators in powers of the electric field exists, continuity must then hold at each order in E(ω).The same must then be true of the expectation values of such operators.This can explicitly be shown to be the case at first order in E(ω); implementing (34,35), we find = 0, given that in principle charge continuity holds in the unperturbed system in a perturbative scheme.
Notably, (34) has a dc divergence taking a form similar to that of (28).Like that second term of ( 28), the second term of (34) vanishes if the unperturbed system is time-reversal symmetric, although this symmetry does not cause the second term of (35) to vanish.Thus, it appears that if one insists on defining electric multipole moments by way of partitioning the electronic charge density into portions that are used to define "site" polarization fields from which "site" multipole moments are extracted and summed to give the full electric multipole moments of the crystal, whether that be via the approach we implement here or some other method, it is unavoidable that one will find a such a dc divergence.In a sense, this unexpected dc divergence is not arising as a consequence of our identification of the polarization, but rather it is inherent to the induced charge density at low frequencies.

V. CONCLUSION
In this work we have considered how polarization and magnetization fields can be defined for metallic systems.In contrast to the approach of the "modern theories of polarization and magnetization," we employ a previously developed strategy [16] for defining microscopic polarization and magnetization fields in general crystalline solids, the macroscopic analogues of which are defined by spatial averaging.Exponentially localized Wannier functions play a central role in how the electronic components of such quantities are defined.In a trivial insulator the macroscopic charge and current densities can be obtained from the macroscopic polarization and magnetization fields alone, both for the ground state and in linear response, while for a metal one would naturally expect contributions from the macroscopic free charge and free current densities, and we have identified them here.
We implemented this approach for a simple instance of a metal, a p-doped semiconductor, initially occupying its T = 0 ground state, and we assume that the Hilbert bundle over the first Brillouin zone associated with any set of isolated energy bands is globally trivial.With this, and because we assume the existence of a band gap above the Fermi energy, contact with expressions for a trivial insulator can readily be reached as a limiting case of the more general expressions we obtain.Indeed, in Sec.III we employ the general definitions in this setting to obtain expressions for P (0) and M (0) , and in the limit of vanishing doping our expressions reduced to those of the "modern theories."While in that limit P (0) is unique modulo a "quantum of ambiguity" and M (0) is gaugeinvariant, this is not so for a metal.Nonetheless, P (0) exhibits the expected property that under translation of the origin of all ELWFs by a Bravais lattice vector R, P (0) is changed by an additive constant proportional to R; M (0) is unaffected by such a translation, and both quantities are unchanged by a shift of the energy zero.
Although the expressions we obtain for P (0) and M (0) agree with the "modern theories" in the limit of a trivial insulator, the two approaches disagree more generally.In the "modern theory of polarization" it has been argued that P (0) is not well-defined in metallic systems [9,11].In the approach implemented here, a definition is always admitted and we obtain a P (0) that is formally similar to that of a trivial insulator.Meanwhile, the "modern theory of magnetization" has been generalized using thermodynamic arguments to obtain an expression for M (0) valid for metals and Chern insulators [11,12], but even so the expression we derive does not agree.This disagreement is not surprising; there is an inherent ambiguity in what one might identify as a magnetization, and the underlying philosophies of these approaches differ.We consider polarization and magnetization to fundamentally arise as microscopic quantities from which macroscopic analogues are obtained, while the "modern theories" view such quantities as being fundamentally macroscopic.These differences are elucidated in the way the expressions for M (0) differ; we find M (0) to be gauge dependent, owing to the central role played by a set of ELWFs in its identification, while in the "modern theory" it is found to explicitly involve a chemical potential, even in the case of a Chern insulator, emphasizing the inherent thermodynamic considerations and the assumed relation to finite-sized systems.In bulk crystals both approaches are valid, each with positive features particularly evident in the domain of considerations that motivate them.Some advantages of the approach implemented here is that the polarization and magnetization are on the same footing, both being defined for all media, that definitions for free charge and current densities are admitted, and that the charge and current densities (1) arise directly from an analysis of the underlying microscopic theory.
In Sec.IV we investigated the linear response of a metallic crystal to an optical field at finite frequency ω, a more general response than is typically considered in the "modern theories."We considered the "long-wavelength limit," within the independent particle and frozen-ion approximations, where the applied electric field is taken to be the macroscopic Maxwell field.Here only P (t) and J F (t) make a contribution to the linearly induced macroscopic current density, J (1) (t) = ∂P (1) (t)/∂t + J While in elementary models of the optical response of metals ∂P (1) (t)/∂t is associated with interband response and J (1) F (t) with intraband response, here we find a more general scenario; in general, that simple association is no longer the case and both contributions are gauge dependent.However, we do find that if all of the energy bands of the unperturbed crystal were isolated from one another then J (1) F (t) would have only intraband contributions and would be gauge invariant, in agreement with those more simple models.Nevertheless, the general σ il (ω) we obtain is gauge invariant and reproduces the usual conductivity tensor of a metal, consisting of a finite-frequency generalization of the "anomalous Hall" and a "Drude" contribution; the latter is entirely due to J F (ω).We also found that if an unperturbed metallic crystal violates time-reversal symmetry, then there is a term in the linear response of the microscopic charge density proportional to ω −1 ; in an approach such as ours that relates the macroscopic polarization to electric dipole moments associated with "site" contributions to the microscopic charge density, this leads to a term in P (1) (ω) proportional to ω −1 .It is the same mechanism that gives rise to the dc divergences of both P (1) (ω) and J (1) F (ω).That is, both divergences involve the second term of (20), which we show in Appendix B is a consequence of an interaction term that gives rise to the intraband response; in this way of identifying inter-and intraband contributions to the linear response we can make contact with earlier work by Blount and others, although the formalism in which we work is indeed much more general.
Given that the association of J F (t) with intraband response and ∂P (1) (t)/∂t with interband response does not hold, that both contributions are gauge dependent, and that P (1) (ω) involves a term proportional to ω −1 , one could argue that a different definition of polarization would be more appropriate.However, such a purported new polarization could not be associated with the dipole moment of microscopic charge densities localized about individual lattice sites.In fact, a more general argument could be made against the philosophy of our investigations.Our goal, a critic might assert, should be to seek what could be taken as "unique" definitions of P , M , F , and J F , and for a metal we do not even demonstrate that for P and M in the ground state.We would reply that such uniqueness is not a reasonable goal.After all, even in the ground state of a trivial insulator the value of P is subject to a "quantum of ambiguity."And once one moves to a general temporal and spatial dependence there are clearly a host of fields P (x, t), M (x, t), F (x, t), and J F (x, t) that could be used to de-scribe the physical quantities (x, t) and J (x, t) via (1).Our perspective is that the focus should be on exploring what might be useful ways of introducing such quantities, for the purpose of both physical insight and calculation.Within that framework this paper can be taken as one such contribution.

VI. ACKNOWLEDGMENTS
We thank Jason Kattan for useful discussions.This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).P. T. M. acknowledges an Ontario Graduate Scholarship.

VII. APPENDICES
Appendix A: Perturbation theory -A strict approach Recall Eq. ( 37) of [16], the general equations of motion for the single-particle density matrix η αR;βR .We here consider the long-wavelength limit, E(x, t) → E(t) and B(x, t) → 0, of those general expressions.
Assuming valid power series expansions for all quantities with respect to the applied electric field E(t), we have αR;βR (t) evolves as the unperturbed single-particle density matrix gs| e i Ĥ0t/ â † αR âβR e −i Ĥ0t/ |gs under the unperturbed Hamiltonian Ĥ0 , as we expect.In particular, starting from the equation of motion for the unperturbed electron Green function i gs| ψ † 0 (y, t) ψ0 (x, t) |gs , the related single-particle density matrix evolves as (A2); the argument is analogous to that which yields (A1) from the "global" Green function (which is related to the minimal coupling Green function by a generalized Peierls phase).Now, via (8,9) the relation between the operators generating ELWFs and those generating the |ψ nk is found to be which we then implement to find which is independent of time.
We now consider (A3) and will closely follow the procedure of Appendix B of [44], however we will not introduce filling factors associated with the ELWFs.It is useful to define the intermediate quantity and from (A3) it follows that i ∂η Then, implementing the usual Fourier analysis via (19), we find where 0 + entering in the denominator describes the "turning on" of the electric field at t > −∞.Finally, using the inverse of (A6), we find Now, using the identity (10) and the result (A5), we find where the first term in brackets results from the two terms in Q µR1;νR2 (ω) that involve dipole moments of the ELWFs, and the second term in brackets results from the term µR1;νR2 (ω).In recasting this second term as a single BZ integral, an integration by parts is performed and all surface terms are taken to vanish.While the integrand is periodic over BZ, it may not be smooth.Thus, this result is valid only if the ground state projector n f nk |ψ nk ψ nk | and therefore n f nk U † αn (k)U nβ (k) for any α, β, which appears in the integrand, is smooth over BZ.While this is always true for insulators -topologically trivial or not -it is here an assumption.However, in the case of p-doped semiconductors considered here, we believe this to be valid if there are no degeneracies at the Fermi energy.With this we arrive at the result (20).Although Eq. ( 20) can be found as the extension of our earlier work presented in Appendix A, we believe some insight can be gained by looking at its derivation using a more traditional perturbation theory approach.
Consider first a molecule, where nuclei are considered fixed and the dynamics of the electron field operator ψ(x, t) follows from the usual minimal coupling Hamiltonian, where p(x) is given previously (7), the applied electromagnetic field is described by the scalar φ(x, t) and vector A(x, t) potentials, and V (x) is the potential energy that confines the electrons to the nuclei.If the wavelength of light is much larger than the molecule, then the electric field E(x, t) can be taken as uniform over the molecule, E(x, t) = E(t), and the magnetic field can be neglected.Via usual strategies [16], it can be shown that the dynamics of the electron field follows from the dipole Hamiltonian where H 0 (x, p(x)) = 1 2m (p(x)) 2 + V (x).The use of (B1) to describe instead the response of the electrons in an infinite crystal to long-wavelength radiation, where H 0 (x) is now taken to be the Bloch Hamiltonian, is a strategy followed by Blount and others [47]; it has even been used to describe the nonlinear optical response of metals [53].The appearance of a position operator in the interaction Hamiltonian requires calculations to be done cautiously, for giving meaning to matrix elements of the position with respect to the Bloch functions of the infinite crystal in a careful way is obviously problematic.In fact, the usual position operator is generally ill-defined to act on the Hilbert space containing such Bloch functions [9].Moreover, it does not seem possible to implement a generalization of this kind of approach to treat instances where the electromagnetic field cannot be approximated as uniform.Indeed, that is one of the reasons the approach applied in this paper was developed.Nonetheless, this strategy does allow for the interaction Hamiltonian to be written as the sum of two terms, which can be identified as "interband" and "intraband."This permits the identification of the interband and intraband contributions to (20), at least within this perspective, and allows us to make contact with earlier work.And so we here present a derivation of (20) using this approach.Although most derivations [47] work in the electronic Hilbert space spanned by Bloch functions from the onset, some issues related to the position operator can be avoided if one works, at least initially, in an isomorphic Hilbert space spanned by a set of exponentially localized Wannier functions; the latter is a subspace of the space of square-integrable functions, where the usual position operator is well-defined.This is the approach we follow here.Moreover, we believe that this approach elucidates the physics of the two terms.Yet we ask the reader to forgive the mathematically questionable steps that are part of the derivation and that are not characteristic of the rest of this paper.We feel that the cavalier approach we take in this Appendix is justified by the insight that the resulting expressions provide.
Working in the Heisenberg picture, the one-body operator on the electronic Fock space related to (B1) is where The primary quantities of interest, the expectation values of the electronic charge and current density operators for a crystal initially occupying its T = 0 ground state |gs , can be extracted from the single-particle electron Green function G(x, y; t) ≡ i gs| ψ † (y, t) ψ(x, t) |gs .

(B3)
We now move from the Heisenberg picture to the interaction picture, wherein operators on Fock space evolve under Ĥ0 and the effect of the perturbation is accounted for in the evolution of the electronic state |ψ(t) = Û(t) |gs , where the time-evolution operator Û(t) is given by [54] Û for VI (t) ≡ −eE a (t) ψ † 0 (x, t)x a ψ0 (x, t)dx.The electron Green function (B3) is then rewritten as Noting that a (complete) set of ELWFs spans the singleparticle electronic Hilbert space, the related operators can be used as a basis with respect to which the electron field operator ψ0 (x, t) can be expanded [55] Implementing (B7) in (B4) and using that result in (B5), we find νR1;µR2 (t) + . . .W * µR1 (y).(B8) Note that in Eq. ( 36) of past work [16] we introduced the single-particle density matrix η αR;βR such that it involved operators generating "adjusted Wannier functions" WαR (x, t) (see Eq. (27,30,33) of Mahon et al. [16]) as well as a generalized Peierls phase Φ(x, y; t) (see Eq. ( 15) there).Thus, in general, it is not the minimal coupling Green function G(x, y; t) to which η αR;βR is "naturally" related, but rather the "global" Green function.However, in the case of a uniform electric field considered here, the corresponding vector potential A is necessarily uniform and Φ(x, y; t) = e c (x − y) • A(t) for a choice of straight-line path in the relators.Then, in this case, Eq. ( 32) of that work simplifies as In this Appendix we employ the gauge choice φ(x, t) = −x • E(t), A(t) = 0, such that the phase Φ(x, y; t) on the RHS of (B9) vanishes.Thus, the identification of the single-particle density matrix in (B8) is consistent with past work.Now, and we thus identify yielding (18).Next consider η νR2;µR1 (t), which from (B8) we identify as We group the second and final lines of (B11) into η νR2;µR1 (t) arising from the second term of (B7).After some algebra we find where we have integrated by parts and taken any surface terms to vanish; this again demands smoothness of the integrand over BZ and requires the same assumption described in Appendix A. We have also taken the electric field E(t) to be adiabatically applied at t = −∞ resulting in the "i0 + " in the denominator and in the phase of e −i(ω+i0 + )t .We now consider η where Notably terms resulting from the first term of the completeness relation (B13), which would involve diagonal matrix elements, cancel one another.Then combining (B12) with (B14) and implementing (19), we find where we have again used an integration by parts.Notably the term in (B15) that diverges in the dc limit arises from the interaction term −eE a (t) αR R a â † αR (t)â αR (t), the second term of (B7).At first one might suspect that it is the sum over Bravais lattice vectors R that leads to the dc divergence, or if not, some other divergence.But, in fact, this is not the case because in the linear response calculation the relevant objects are of the form R R a gs| â † µR1 (t)â νR2 (t)â † αR (t )â αR (t ) |gs (see the first equality of (B12)); thus not all R's contribute equally and the result of such a sum appears to be finite.
To gain further insight into origin of the terms appearing in (B15), it is useful to rewrite VI (t) in terms of the operators that generate the single-particle Bloch energy eigenvectors.The second term of (B7) involves (B16) When implemented in the linear response calculation, the first two terms of (B16) give non-zero contributions only for those k "near" the Fermi surface and indeed gives vanishing contribution if |gs is the ground state of a trivial insulator.That such an interaction term leads to a diverging induced free current density is in-line with physical expectation.The first term of (B7) can also be rewritten,  which is gauge independent, as expected.As described above, due to the relative negative sign between terms involving VI (t) and V † I (t) in the perturbative expansion of the electron Green function, the interaction term involving E a (t) BZ dk nm â † nk (t)ξ a nm (k)â mk (t) gives rise only to terms for which n = m, which we refer to as being related to the "interband response" (see, for example, the cancellation of "intraband" terms in (B14)).In contrast, we refer to the terms resulting from the interaction term involving iE a (t) BZ dk n â † nk (t) ∂ a ânk (t) − ∂ a â † nk (t) ânk (t) as being related to the "intraband response." B (ω = 0) T = 0, or equivalently that the term in (28) that diverges in the dc limit vanishes.Moreover from the relations (C1,C2) it follows that M (0) T = 0.
Appendix D: Link currents and the related free current density Recall from past work [16] that in the "long-wavelength limit" αR ;λR (ω), with all higher order contributions vanishing in this case, we identify The first line of both equalities of (D4) is the result of a "compositional" modification, while the remainder is the result of a "dynamical" modification; the second term of (31) is the result of the first term of (20) and the final term of ( 31) is the result of the second term of (20).Notably the first line of ( 31) is independent of energy and involves frequency only through E(ω), while this is generally not the case for the other terms.In Sec.IV we are interested, among other things, in the macroscopic free current density, J F (x, ω), related to the microscopic free current density j F (x, ω).In past work [22] we have described this averaging procedure in some detail, in particular for the microscopic polarization and magnetization fields.In the limit of a uniform applied electric field, the expressions Eq. ( 7), ( 9), (B4)-(B6), and (B8) presented there result in the macroscopic polarization and magnetization fields being uniform, and the only contributions being the dipole moments, (21).We here focus on the macroscopic free current density found by implement a spatial averaging function w(x) to relate the microscopic and macroscopic quantities.That is, J F (x, ω) ≡ w(x − x )j F (x , ω)dx . (D5) Implementing the definition (30), the relator expansion [22] s i (w; x, y) (x i − y i )δ(w − y) − 1 2 (x i − y i )(x j − y j ) ∂δ(w − y) ∂w j + . . ., (D6) and noting that the first-order modification to the link currents here takes the form I (1) (R, R ; ω) = I (1) (R − R , ω), we find where in going to the final line we have used the special case of a uniform applied electric field in Eq. (B8) of [22].Thus, we arrive at (32).