Symmetry, topology, and geometry: The many faces of the topological magnetoelectric effect

A delicate tension complicates the relationship between the topological magnetoelectric effect (TME) in three-dimensional (3D) $\mathbb{Z}_2$ topological insulators (TIs) and time-reversal symmetry (TRS). TRS underlies a particular $\mathbb{Z}_2$ topological classification of the electronic ground state of crystalline band insulators and the associated quantization of the magnetoelectric response coefficient calculated using bulk linear response theory but, according to standard symmetry arguments, simultaneously forbids a nonzero magnetoelectric coefficient in any physical finite-size system. This contrast between theories of magnetoelectric response in formal bulk models and in real finite-sized materials originates from the distinct approaches required to introduce notions of (electronic) polarization and orbital magnetization in these fundamentally different environments. In this work we argue for a modified interpretation of the bulk linear response calculations in non-magnetic $\mathbb{Z}_2$ TIs that is more plainly consistent with TRS, and use this interpretation to discuss the effect's observation - still absent over a decade after its prediction. Motivated by analytical results, we conjecture a type of microscopic bulk-boundary correspondence: a bulk insulator with (generalized) TRS supports a magnetoelectric coefficient that is purely itinerant (which is generically related to the geometry of the ground state) if and only if magnetic surface dopants are required for the TME to manifest in finite samples thereof. We conclude that in non-magnetic $\mathbb{Z}_2$ TIs the TME is activated by magnetic surface dopants, that the charge density response to a uniform dc magnetic field is localized at the surface and specified by the configuration of those dopants, and that the TME is qualitatively less robust against disorder than the integer quantum Hall effect.


I. INTRODUCTION
One of Peierls' Surprises in Theoretical Physics [1] is that the orbital magnetization of metals can be correctly calculated using infinite lattice models that neglect the surfaces of realistic material samples, even though that magnetization can be understood to arise from bound currents at those surfaces.This mysterious success is now often taken for granted.In recent times, strong interest in the multitude of topological insulators (TIs) and the linear response thereof has highlighted similar issues, in particular when disentangling the roles played by the topology ascribed to the bulk electronic ground state and by topologically protected electronic surface states.The integer quantum Hall effect (IQHE) provides a familiar example.In a two-dimensional (2D) crystalline band insulator the IQHE occurs when the vector bundle of occupied electronic Bloch states over the Brillouin zone (BZ) torus is characterized by a nonzero Chern invariant [2][3][4].Meanwhile, experiments are often interpreted in terms of the gapless chiral edge states [5][6][7][8][9] whose presence is * perry.mahon@austin.utexas.edumore immediately related to charge conduction.In the IQHE, it is generally accepted that this bulk topology of a band insulator implies a quantized bulk Hall conductivity but also implies the existence of chiral edge states that yield a consistent conductance in finite-size samples thereof, unifying the two interpretations.Notably, breaking time-reversal symmetry (TRS) is necessary for a nonzero Chern invariant and, by the usual symmetry arguments, for a nonzero Hall conductivity in finite-sized systems.The topological magnetoelectric effect (TME) presents an even more stark puzzle.In three-dimensional (3D) crystalline band insulators, TRS leads to a Z 2 topological classification of the electronic ground state and in Z 2 -odd phases to a magnetoelectric linear response coefficient that is argued [10,11] to equal (n + 1/2)e 2 /hc for n ∈ Z.On the other hand, in any finite-sized system with TRS, the usual symmetry arguments dictate that the magnetoelectric coefficient must vanish.Indeed, the requirement of magnetic surface dopants for realization of the TME has been noted previously [10][11][12][13][14].There is no physical bulk response directly related to the 3D Z 2 invariant, but a non-trivial Z 2 invariant implies the existence of surface states that can then be gapped to activate the TME.This more convoluted connection has implications for the experimental robustness of the effect, making it qualitatively less robust against disorder compared with the IQHE.
It seems therefore that the magnetoelectric response of non-magnetic Z 2 TIs is a rare violation of the Peierls principle referred to above, i.e. that evaluating a quantity using a band theoretic description does not always correctly produce the value of that quantity in a large finite-sized sample thereof (see Fig. 1).In this paper we explicitly address the relationship between the magnetoelectric linear response coefficient of a bulk insulator that exhibits TRS and that of a finite size sample counterpart.As always, the bulk crystal is merely a convenient theoretical construct.Our interest is in understanding when it yields the physically correct magnetoelectric response.In Sec.II we describe this conundrum in more detail and highlight the main results of the calculations that follow.We argue for a slightly different interpretation of the bulk magnetoelectric coefficient derived in previous work [15], which is informed by the relationship (see Fig. 1) between particular linear response coefficients in finite-sized samples and bulk that we explicitly establish in Sec.III and IV.Since the magnetoelectric coefficient lacks physical significance in macroscopically uniform bulk insulators that exhibit TRS, the bulk interpretation that we propose implicitly involves the consideration of surfaces, a seemingly unavoidable feature.In this interpretation the magnetoelectric coefficient adheres to the Peierls principle and we reach a conclusion that is manifestly consistent with TRS; in non-magnetic Z 2 TIs, magnetoelectric response occurs only when the surface states are gapped by magnetic dopants, the magnetoelectric coefficient is quantized only when the surface magnetization configuration satisfies stringent conditions, and currents proportional to the magnetic field can flow through the material when the surface magnetization profile is changed.
In Sec.III we extract the magnetoelectric coefficient for thin film and bulk samples of non-magnetic Z 2 TIs in the V 2 VI 3 family of materials by directly including a uniform dc magnetic field in a coupled-Dirac cone model of the low-energy electronic states, then computing the polarization of the occupied energy eigenfunctions.In the case of thin films, we find that when the Dirac cones associated with the top and bottom surfaces are gapped by magnetic dopants of opposite magnetization, the magnetoelectric coefficient α me is nonzero and approaches a quantized value as the film thickness is increased.The induced polarization is sensitive to the magnetic dopant configuration and α me moves continuously through zero to −α me as the component of magnetization perpendicular to the surface is continuously taken to its opposite value.When there is no surface magnetization, there is no magnetoelectric response.We show that the quantization of α me in the thick film limit can be understood in terms of the topological properties of a Su-Schrieffer-Heeger model that arises in the bulk description.
In Sec.IV we employ the 3D tight-binding model that we obtain from a lattice regularization of the coupled-

TRS preserving
FIG. 1. Relationship between αCS in bulk insulators that exhibit time-reversal symmetry (TRS) and αme in thin films thereof.In the bulk, TRS leads to a Z2 topological classification of the electronic ground state, to which αCS is sensitive.In three-dimensional Z2-odd phases, αCS is quantized at e 2 /2hc up to an integer multiple of e 2 /hc.In non-magnetic thin films, however, the physically realized magnetoelectric response αme is nonzero only when the energy dispersion of the topologically protected surface states is everywhere gapped by magnetic surface dopants which break TRS locally.
Dirac cone model used in Sec.III and calculate α CS , the purported magnetoelectric coefficient in crystalline band insulators that exhibit TRS, semi-analytically using the well-known bulk linear response expression; we reproduce the expected quantization.Calculating α CS is technically challenging because it must be evaluated with respect to a smooth global gauge of the vector bundle of occupied Bloch states over the 3D BZ.Generically, Bloch energy eigenvectors that are smooth over the entire BZ do not exist in a Z 2 TI [16] and a Wannierization-like procedure [17] is required to obtain an adequate gauge choice [18].Fortunately, the model we employ exhibits a fermionic time-reversal symmetry and an inversion symmetry, thus each eneregy band is at least two-fold degenerate over the BZ.In fact, the energy bands are everywhere pairwise isolated.Thus there exists a smooth global Hamiltonian gauge, and with respect to such a gauge we make our calculations.In Sec.V we analytically demonstrate that, in a particular smooth global Hamiltonian gauge, α CS is entirely attributed to the itinerant portion of the electricfield-induced orbital magnetization, and speculate that this is a universal feature unique to non-magnetic Z 2 TIs.Our calculations explicitly demonstrate that for Z 2 TIs, the bulk theory of magnetoelectric response does not always capture the properties of large finite size samples.

II. LINEAR RESPONSE THEORY
The interaction between macroscopic electromagnetic fields and the charged constituents of material media is often described at the level of linear response by intro-ducing phenomenological susceptibility tensors that are nonlocal in space and time, and relate changes in the (macroscopic) charge ϱ(r, t) and current J (r, t) densities of the material to the applied electric E(r, t) and magnetic B(r, t) Maxwell fields that induce those changes.The charge and current densities and the Maxwell fields are understood to be coarse grained [19] over a length scale intermediate between the atomic spacing and the wavelength of light; indeed local field corrections, which can be important, are not of interest here and are neglected.It can often be assumed that a bulk material is spatially uniform at that intermediate length scale, in which case the susceptibility tensors are translationally invariant.Then, for example, the linear response of J (r, t) is of the form J i(1) (q, ω) = σ il (q, ω)E l (q, ω), where σ il (q, ω) is the effective conductivity tensor, q and ω are respectively the wave vectors and frequencies that arise in the Fourier transforms of E(r, t) and B(r, t), and superscript indices here and below identify Cartesian components that are summed over when repeated.
The first term in Eq. ( 3) is the familiar long-wavelength frequency-dependent conductivity tensor.Using Faraday's law, the second term can be shown to include contributions that involve B(r, ω) and symmeterized spatial derivatives of E(r, ω).Any susceptibility in a material that is uniform at the course grained length scale can be analyzed in this way.Physical insight into the distribution and dynamics of the charged (quasi-)particles that constitute a material medium can be gleaned by identifying bound and free contributions to ϱ(r, t) and J (r, t), and associating the former with macroscopic polarization P(r, t) and (orbital) magnetization M (r, t) fields [20] such that ϱ(r, t) = −∇ • P(r, t) + ϱ F (r, t), For a given ϱ(r, t) and J (r, t) there is always ambiguity in defining P(r, t), M (r, t), ϱ F (r, t), J F (r, t) that satisfy Eq. ( 4).Nevertheless, let us assume that definitions for P(r, t) and M (r, t) in a bulk material have been made and that the applied Maxwell fields are in the linear response regime so that expansions of the electric and magnetic dipole moments [21] as sums of spontaneous and field-dependent contributions, are justified [22].At low frequencies ω, α il P (ω) and α il M (ω) are well approximated by their static ω = 0 values, which are in fact related by a thermodynamic Maxwell relation In this case, using Eq. ( 5) in ( 4) and comparing with Eq. (2) yields where ". .." denotes purely electric contributions to σ ilj .Explicit forms for phenomenological susceptibilities can be obtained from quantum linear response theory.The semi-classical Hamiltonian of the coupled lightmatter system consists of electron and Maxwell contributions, and an interaction term that results from minimal coupling.The latter involves the electromagnetic scalar and vector potentials multiplied by charge and current density operators that are obtained from components of the Noether 4-current of the electron theory.In a finite sized material those charge and current densities are nonzero only in a localized region of space, in which case there exists well-motivated definitions for polarization and magnetization fields developed by Power, Zinau and Wooley (PZW) [23], Healy [24], and others.This approach begins with a unitary transformation of the minimal-coupling Hamiltonian (differential operator) to yield the physically equivalent PZW Hamiltonian.After defining the polarization and magnetization fields, the interaction of the electronic degrees of freedom with the electromagnetic field as described in the PZW Hamiltonian involves terms in which the polarization and magnetization fields are multiplied by E(r, t) and B(r, t), respectively [25]; notably this formulation does not require an explicit choice of the electromagnetic gauge.When a multipole expansion of those polarization and magnetization fields is made, the PZW Hamiltonian takes the classically anticipated form of a multipole Hamiltonian [26].Thus, the PZW definitions for polarization and magnetization fields are physically reasonable, and from them the (minimally-coupled) electronic charge and current density expectation values, and thus the finite size sample analog of the susceptibility tensors mentioned above (which are position dependent), can be rigorously obtained.Moreover, that multipole expansion results in Hermitian operators in the electronic Hilbert space; that is, the electric and magnetic multipole moments are genuine physical observables.For example, the ground state expectation value of the electronic contribution to the PZW electric dipole moment per unit volume (which is equal to the interior polarization of the material) is where e > 0 is the elementary charge, Ω is the finite volume of the sample, E F is the Fermi energy, and Ψ E (r) are the electronic energy eigenfunctions of the unperturbed Hamiltonian [27].
Unfortunately the PZW approach cannot be directly applied to crystalline solids since, even in the absence of electromagnetic fields, the electronic charge and current density expectation values are expressed in terms of Bloch energy eigenfunctions, the support of which is all of space [28].Indeed, this is related to the property that the usual position operator is not well-defined in the Hilbert space of Bloch functions [29].There are a variety of different approaches [29][30][31][32][33][34] that can be used to extend the notions of polarization and magnetization to bulk crystals.Inspired by the PZW approach, Sipe et al. have developed a formalism [34][35][36] applicable to general extended systems, crystalline or otherwise, which has been employed to account for spatial variation of electromagnetic fields in crystal insulators [37,38], and more [39][40][41][42][43].However, none of these bulk crystal approaches are able to define electric and magnetic multipole moments as genuine physical observables in the PZW sense.
The most commonly used approach to sidestep difficulties in defining the electric and magnetic dipole moments in crystalline insulators is the so-called modern theories of polarization [30,44] and (orbital) magnetization [31,32].Focusing on the former, one aims to deduce an electric dipole moment from calculation of the current density, rather than to propose a general definition for it.To do so, the macroscopic free current density (appearing in the second of Eq. ( 4)) is assumed to vanish at linear response in typical insulators, which seems to be a sensible assignment, thus the current density is related only to the polarization and magnetization fields.It is less obvious that this assignment is sensible in Chern insulators, but materials of that type are not considered here.Next, focus is restricted to describing the influence of electric and magnetic fields that are spatially uniform, in which case all of the macroscopic densities appearing in Eq. ( 4) are also expected to be uniform.Thus the polarization and magnetization fields are entirely described by their dipole moments, and J (t) = ∂P (t)/∂t.
The magnetoelectric susceptibility α il can generally be written as the sum of an isotropic δ il α θ and a traceless contribution [15].Examining Eq. ( 7) makes it evident that in a static and macroscopically uniform bulk material, the σ ilj tensor is insensitive to the former.That is, α θ cannot be determined by calculating the wave vector dependent optical conductivity.In other words, the physical implications of α θ are not equivalent to those of a wave-vector dependent conductivity.
If the isotropic magnetoelectric response cannot be calculated from the bulk current response to non-uniform electric fields, how can it be found?Two independent strategies that can be used to obtain α θ have been identified in the literature: (i) Essin et al. [15] consider an insulating state of a Bloch Hamiltonian that is adiabatically varied in time, such that α il becomes time-dependent and there is an additional contribution to Eq. ( 7) (originating from ∂P /∂t) proportional to (∂α il /∂t)B l ; (ii) Malashevich et al. [45] consider a finite sized material, in which case α il becomes position-dependent and there is an additional contribution to Eq. ( 7) (originating from ∇ × M ) that is proportional to (ϵ iab ∂α bl /∂r a )E l .Both approaches reach the same conclusion, which is that in a crystalline insulator initially occupying its zerotempetature electronic ground state where the explicit form of α il G resembles that of a conventional linear response tensor, involving cross-gap matrix elements between Bloch states that are initially occupied and unoccupied, and α CS (the Chern-Simons contribution) can be expressed in terms of occupied Bloch states alone.If the bulk insulator exhibits TRS then α il = δ il α CS and is quantized since α CS evaluates to an element of a discrete lattice of values [10], which does not include zero in the case of Z 2 TIs.It is at the step that the conundrum mentioned in the Introduction is introduced, since this is where the conclusion is reached that α ii mod e 2 /hc [46] can be nonzero in insulators with TRS.
The strategy (i) of Essin et al. [15] considers a 3D crystalline insulator described by a Bloch Hamiltonian that varies in time via adiabatic variation of the crystal Hamiltonian.They show that in the presence of a static and uniform magnetic field B, adiabatic linear response to slow changes in the crystal Hamiltonian induces a spatially uniform macroscopic (minimal-coupling) electronic current density of the form [47]. Unlike above, where the time dependence of J (1) (t) resulted from that of the dynamical electromagnetic field, B is here taken to be static and the t-dependence results from that of the Bloch Hamiltonian itself; t-dependence manifests implicitly in the Bloch energy eigenvectors and eigenvalues.Although their strategy applies more generally, we focus on the special case in which the Bloch Hamiltonian exhibits TRS at every instant t.Then J The explicit expression for Q 3 (k; t), the Chern-Simons 3-form [49] at instant t, is discussed in Sec.IV.It is important to note that Q 3 (k; t) explicitly involves only the Berry connection and curvature that are defined on the vector bundle of occupied Bloch states over the Brillouin zone torus at instant t.That is, at every t, Q 3 (k; t) is a purely geometrical quantity.In Ref. [15] the derivative with respect to time in Eq. ( 10) is taken outside of the Brillouin zone integral to obtain where At each instant t, is determined (modulo 2) by the TRS-induced Z 2 topological classification [10,50,51] of the electronic ground state [52] and the value of α CS (t) is therefore fixed (modulo e 2 /hc) by that topology; if the ground state is Z 2 -even then α CS evaluates to 0 mod e 2 /hc and if it is Z 2 -odd then α CS evaluates to e 2 /2hc mod e 2 /hc [10,11].Importantly, Q 3 (k; t) and indeed α CS (t) are gauge dependent -they are both sensitive to the smooth global frame of the vector bundle of occupied Bloch states over the Brillouin zone torus that is used for their evaluation -which leads to the above described discrete ambiguity.Therefore, even when considering instants t 1 and t 2 between which a band insulating ground state remains in the same topological class, α CS (t 1 ) need not equal α CS (t 2 ) since at each t one has the freedom to choose one of many gauges and with respect to each of those gauge choices a different value of α CS (t) ∈ (e 2 /2hc)Z can result; in this scenario it is only guaranteed that α CS (t 1 ) mod e 2 /hc equals α CS (t 2 ) mod e 2 /hc.To obtain Eq. (10) it must be assumed that the gauge choice used to express Q 3 (k; t) is smooth in k ∈ BZ for each t.And in order to move from Eq. ( 10) to (11), that gauge choice must also be continuous in t [53].The latter condition implies that α CS (t) is continuous in t.If there is TRS at each t then α CS (t) ∈ (e 2 /2hc)Z and thus α CS (t) must be constant in t.Indeed our explicit calculations in Sec.IV demonstrate that once a gauge choice is made at any t, the value of α CS (t) is fixed for all t if that gauge is continuous in t.Then under any TRS preserving adiabatic variation of the Bloch Hamiltonian in time, the right hand side of Eq. (11) evaluates to zero; no bulk currents flow when the parameters of the Hamiltonian change within the space of a given topological phase and the adiabatic approximation is valid.Indeed this result should be expected since Eq. ( 11) is consistent with TRS only if nonzero values for ∂α CS (t)/∂t are activated by some time-reversal symmetry breaking perturbation.
If a Bloch Hamiltonian is adiabatically varied in t without necessarily exhibiting TRS at each t, then α CS (t) need not evaluate to an element of (e 2 /2hc)Z for each t [10] and therefore α CS (t) is not necessarily constant in t.For example, if there is TRS at some initial t i and final t f times, then α CS (t i ), α CS (t f ) ∈ (e 2 /2hc)Z and if the electronic ground state at t i (t f ) is Z 2 -even (Z 2 -odd), then t f ti J (B) (t)dt mod e 2 B/hc equals e 2 B/2hc.Crucially, this ambiguity in t f ti J (B) (t)dt is not a consequence of the gauge dependence of α CS , but rather a manifestation of the fact that the current pumped through the material depends on how the model parameters are changed in time.Indeed if the Bloch Hamiltonian evaluated at t i and at t f equal, then one can define a vector bundle of occupied electronic states over BZ × S 1 and the topology of this structure is characterized by a second Chern invariant which determines [10,54].A similar situation arises in the case of the unperturbed bulk 1D polarization, where a certain first Chern invariant determines the bound charge [55,56].
The linear response result (11) for crystalline insulators that exhibit TRS has been interpreted [15] by identifying the magnetic-field-induced bulk polarization P (B) (t) with α CS (t)B.For a static insulator of interest, its bulk P (B) is then obtained by evaluating α CS (t)B at any time t at which the time-dependent Bloch Hamiltonian describes that insulator.If we accept this identification, then P (B) can differ from the physically realized electric dipole response in a large but finite sample [57].For example, consider a finite crystallite of a non-magnetic Z 2 TI that has TRS everywhere.TRS implies that there is no electric dipole induced by B in linear response.The bulk expression for the polarization is evidently misleading in the absence of magnetic dopants somewhere.Because this identification of P (B) does not imply induced bulk charge or current densities it is still technically acceptable as a bulk quantity, but associating physical significance with it requires additional insight.Indeed the necessity of magnetic surface dopants for the topological magnetoelectric effect to manifest in non-magnetic TIs has been noted [10,11,57,58] previously.A ubiquitous aspect of bulk theories of polarization and orbital magnetization is the importance of interpretation [56].We believe it is useful to formulate an alternative, physically equivalent, perspective to that of Ref. [57].By combining the known topological constraints in bulk with the known symmetry constraints in finite sized systems, we produce below a physically meaningful notion of P (B) in bulk insulators that exhibit TRS, one that adheres to the Peierls principle and implicitly accounts for additional surfacerelated criteria.To that end, the V 2 VI 3 family of TIs are an ideal test bed on which to develop that notion.
We have argued above that even in the strategy (i) of Ref. [15], wherein the Bloch Hamiltonian for a bulk insulator involves crystal parameters that depend on time, α CS (t)B lacks physical consequences if that Hamiltonian exhibits TRS at each t.Thus, in this setting any physically significant identification of P (B) (t) must involve the consideration of surfaces.As a result, any such identification will be ad hoc in that one has to artificially include the role played by the surface in a proposed bulk quantity.This is an aesthetically unappealing, but a seem-ingly unavoidable, feature of any physically meaningful identification of P (B) (t) in bulk insulators with TRS.
One such identification of P (B) (t) could be obtained in a manner similar to that of Essin et al. [15] -by considering the magnetic field dependent current J (B) (t) that flows in a material as its properties are changed -but, rather than considering a time-dependent Bloch Hamiltonian, consider a Hamiltonian that describes a macroscopically large but finite sized material with timedependence only at its surfaces, where magnetic moments order to break TRS and open gaps in the Dirac-like energy dispersion of the surface states.(TRS is maintained in the interior region.)If the effect of interactions that break TRS is local [59], then charge density response to B can only occur near the surfaces.Our interpretation for non-magnetic Z 2 TIs [60], which is supported by later calculations in this paper, imagines replacing the bulk quantity α CS (t) by a new function α expt (t) that compares more directly to experiment [61] by implicitly accounting for the role of surface magnetic dopants.The (physically meaningful) bulk polarization in these materials is where α expt (t) must have the following properties: 1.In the absence of magnetic surface dopants J (B) (t) vanishes in finite samples, which is consistent with the bulk prediction (11).In this case, TRS dictates that α expt (t) ≡ 0 and the identification of P (B) (t) with α CS (t)B for Z 2 TI materials is incorrect.
2. (green regions in Fig. 2) During time intervals over which the magnetic dopant configuration on top and bottom surfaces is opposite and such that the energy dispersion of the surface states remains gapped, the physically realized response is that of the bulk; the prediction (11) that J (B) (t) vanishes and the identification of P (B) (t) with α CS (t)B are correct (the latter modulo e 2 B/hc), i.e. α expt (t) = α CS (t).P (B) (t) is constant in t as long as the surface states remain gapped.
3. (red region in Fig. 2) During time intervals over which the magnetic dopant configuration is such that there is a finite surface density-of-states, a current proportional to B flows through the material.Indeed, the bulk prediction (11) using the adiabatic approximation fails to describe these finite size samples and α expt (t) ̸ = α CS (t).That current, and the pattern of polarization that develops, is sensitive to the magnetic configuration of the surface dopants.For realistic samples with potential and magnetic disorder, the surface gaps will remain closed over a finite range of dopant configurations.Under surface magnetization reversal, currents proportional to B will flow through the material during the time interval over which there is no surface energy gap.α expt (t) now has a contribution that is Fully Gapped TRS TI Fully Gapped TRS TI Ungapped Surfaces 2. Relationship between the bulk magnetoelectric coefficient αexpt(t) (black line) that we propose and αCS(t) (blue line) as obtained from a simple interpretation of the bulk current that is adiabatically induced by slow changes of a Bloch Hamiltonian describing a non-magnetic bulk Z2 TI.The illustrated time interval is imagined to be one over which the magnetization describing a configuration of surface magnetic dopants that initially gaps the energy dispersion of the surface states is changed to its opposite.The red region, where this magnetization vanishes, is given finite width for illustrative purposes.We imagine making a particular gauge choice at some initial time such that αCS(t) exactly agrees with αexpt(t) until a surface gap closes.A purely bulk linear response calculation [15] uses the adiabatic approximation to show that αCS(t) is time-independent and the (magnetic field dependent) bulk current J (B) (t) vanishes (no currents flow) when an instantaneous band insulating state is varied in a manner that maintains TRS.We argue that in finite samples J (B) (t)dt is insensitive to gapped surfaces and is therefore given correctly by the bulk linear response calculation at early and late times within the illustrated time interval.The adiabatic approximation fails to describe finite samples when the surface gap vanishes, in which case the physically realized J (B) (t) is nonzero and the polarization can change.In practice, the sub-interval within which J (B) (t) ̸ = 0 will always be finite because of disorder at the surface.We argue that the interior polarization at a particular point along the surface is a local quantity [10] that changes sign with the local orientation of magnetic dopants, in agreement with TRS.
sensitive to the surface magnetization, which allows the polarization to interpolate between the topologically allowed values of gapped states.
Our interpretation is that during time intervals t ∈ [t i , t f ] over which the surface magnetization is reversed, the value of P (B) (t) changes and these changes are associated with the flow of currents through the bulk of the material.Varying the configuration of magnetic surface dopants in such a way that the total magnetization always vanishes but that surface gaps close then reopen induces bulk currents in much the same way that varying crystal parameters such that at intermediary times TRS is broken in the bulk [62] does; notably, neither mechanism yields a bulk charge density.The change in polarization may be obtained by integrating (∂α expt (t)/∂t)B over time, which, due to TRS in the bulk, will integrate to the difference [63] between two quantized values of α CS (t) that differ in sign.Since in our interpretation we assume that the TRS-induced bulk topology is unaffected by the manipulation of the magnetic surface dopants, α expt (t f ) − α expt (t i ) ∈ (e 2 /hc)Z.If the magnetization of the dopants at t f is opposite to that at t i , we require Under this identification, the magnetoelectric response in non-magnetic Z 2 TIs is nonzero only when magnetic surface dopants are present, solving the TRS conundrum.If there is no surface magnetization, then there will never be a surface energy gap everywhere and no bulk current flows in linear response to a magnetic field.In the following sections of this paper we describe electronic polarization (orbital magnetization) calculations in model non-magnetic Z 2 TIs subject to a uniform dc magnetic (electric) field, which support this interpretation.In particular, in sufficiently thick films the polarization changes and interior currents flow only when the surface-normal magnetization component is varied through zero.

III. MAGNETOELECTRIC RESPONSE IN Z2 TOPOLOGICAL INSULATOR THIN FILMS
In this section we employ a simplified but realistic model of the electronic states near the Fermi energy in thin films of V 2 VI 3 -type non-magnetic Z 2 TIs in which it is possible to account for the presence or absence of magnetic surface dopants.We calculate the magnetoelectric coefficient α me by introducing a magnetic field that is perpendicular to the thin film and accounting for the Landau quantization it produces.We begin in Sec.III A by considering a toy version of this model in which the only electronic states retained are the topologically protected surface states modeled by two coupled 2D Dirac cones, one associated with each surface.The surface magnetic dopants couple to these Dirac-like states separately via an exchange mass m, and for finite thickness films those states are hybridized by an inter-surface tunneling parameter ∆.This model is amenable to analytic analysis, from which it is found that in the limit m → 0 α me → 0, as expected for finite-sized material samples with TRS, and in the limit ∆ → 0 α me → e 2 /2hc, as expected for bulk Z 2 TIs.In Sec.III B we consider a more realistic model, which involves many coupled Dirac cones, and reach the same conclusions.We show that the quantization of α me in the thick film limit can be understood in terms of the topological properties of a Su-Schrieffer-Heeger model that arises in the bulk description.

A. Simplified toy model of the surface states
The analytic calculations in this section illustrate the essential microscopic physics of magnetic-field-induced charge redistribution in finite thickness thin-films of nonmagnetic Z 2 TIs.In materials of this type, the electronic states that are nearest to the Fermi energy occur about the k Γ point (k Γ,x , k Γ,y ) = (0, 0) and are localized at the surface of the sample [64,65].Indeed the bulk topology of a Z 2 TI implies the existence of an odd number of Dirac points at each surface of a finite sample thereof [66,67].A model of a thin film can therefore be constructed using two copies of a 2D k • p Dirac Hamiltonian, one associated with each surface.We account for the (finite) film thickness in the surface-normal direction (taken to be z) by hand by endowing the basis states of each copy of the k • p Hamiltonian with z-dependence to encode its association with a particular surface.
We therefore consider two Dirac Hamiltonians written in a basis of states ū(α,σ),k=kΓ that are orthogonal linear combinations of the k • p basis states and are associated with the bottom (α = 0) and top (α = 1) surfaces of the film.At k Γ , the k • p Hamiltonian will commute with s z , thus the corresponding spin-up and spin-down eigenvalues σ = ±ℏ/2 (≡ ↑, ↓) can be used to identify Hilbert space basis states at k Γ (and therefore label basis states of any k • p Hamiltonian about k Γ ).We include an exchange coupling m α associated with each surface, here describing the effect of magnetic dopants at that surface, that couples to the s z spin components of the states at the surface α; this interaction breaks TRS.We take the exchange mass to be of opposite [68] value at the top and bottom surfaces, m 0 = −m 1 ≡ m [58,69].We also include an inter-surface hybridization parameter ∆.In the basis ( ū(0,↑),kΓ , ū(0,↓),kΓ , ū(1,↑),kΓ , ū(1,↓),kΓ ), the Hamiltonian is specified by [70] where k ± ≡ k x ± ik y .
To account for the presence of a uniform dc magnetic field that is parallel to the surface-normal of the thin film, we first implement the usual prescription [71] to obtain an envelope function approximated (EFA) Hamiltonian that is related to a given k • p Hamiltonian.Following this, we minimally couple the electronic degrees of freedom to the corresponding vector potential.That is, first take ℏk → p(r) ≡ −iℏ∇ followed by p(r) → p mc (r) = −iℏ∇ + e c A(r), where −e is the electronic charge.We take B = (0, 0, −B) and in the Landau gauge A(r) = (0, −Bx, 0).The 2D EFA Hamiltonian that is obtained by employing this prescription to Eq. ( 15) is invariant under translations along y within the film.Thus, we seek energy eigenfunctions of the form Introducing the usual differential ladder operators a ≡ [72], where x ≡ l B q y − x/l B , l 2 B ≡ ℏc/eB, and [a, a † ] = 1, the EFA Hamiltonian that is related to (15) acts on the Φ E,qy (x) via where ω c ≡ √ 2v D /l B , τ a are the Pauli matrices acting on the surface label (i.e.orbital type) degree of freedom, and σ a are the Pauli matrices on spin degree of freedom.
A general eigenfunction of Eq. ( 17) has the form Φ En,qy (x) = Φ n (x), where for integers n > 0 and for n = 0, where χ n,qy (x) = χ n (x) are normalized eigenfunctions of a † a.The corresponding eigenvalues are with the n = 0 eigenvalues non-degenerate and the n > 0 eigenvalues each two-fold degenerate.We adopt the notation Φ ± 0 (x) and Φ ±,1 n (x), Φ ±,2 n (x) for the corresponding eigenfunctions.
We endow the Ψ ± n,qy (x, y) with z-dependence by multiplying the components of (18,19) that are associated with the top (bottom) surface of the thin film by δ(z − z t(b) ), where z t = a/2 and z b = −a/2.Then, taking E F = 0 such that the (Ψ + n,qy (r)) Ψ − n,qy (r) are (un)occupied, and using qy (1) = qy |χ n (x)| 2 dx = eBL x L y /hc, the number of q y per Landau level n, employing Eq. ( 8) yields from which we can identify α me ≡ α zz via P z el = α zz B z .The n = 0 anomalous Landau levels Φ ± 0 (x) are spin polarized, and the action of ( 17) on (19) simplifies such that C ± 0,(b,↑) , C ± 0,(t,↑) is obtained by diagonalizing mτ z + ∆τ x .At half filling, the occupied n = 0 eigenfunction is In the bulk limit of this toy model, the top and bottom surface states decouple and ∆ → 0, in which case Eq. ( 22) yields the quantized magnetoelectric coefficient e 2 /2hc independent of m.Meanwhile in thin films with TRS, which implies m = 0, Eq. ( 22) gives that α n=0 me vanishes.Of course, the occupied n > 0 Landau levels may also contribute to (21).Unfortunately, the eigenvectoreigenvalue equations are not as simple as the n = 0 case.Moreover, due to the two-fold degeneracy of each eigenvalue, the energy eigenvectors are non-unique.One convenient non-orthogonal pair is [73 where N ±,1 n and N ±,2 n are normalization factors.These eigenvectors can be orthogonalized and their contributions to the induced polarization calculated.In the limit ∆ → 0, which is of primary interest, these eigenvectors are indeed orthogonal and are moreover equally weighted combinations of states that are associated with the top and bottom surfaces.Thus, at half filling and in the limit ∆ → 0, it is only the Φ − 0 (x) anomalous Landau level that contributes to the magnetic-field-induced polarization.

B. Coupled Dirac cone model of layered thin films
The results for α me obtained using the simple model described in Sec.III A can be confirmed by considering a more realistic model for finite thickness films of V 2 VI 3type TIs.Materials in this family consist of many stacked 2D layers.They are more accurately modelled by introducing a pair of Dirac cones for each layer, one associated with its top and bottom surface, then coupling the Dirac-like states within the same layer via ∆ S and coupling those associated with the closest surfaces of the nearest-neighbor layers via ∆ D [74,75].This model has the merits that it is readily solved numerically for films that contain many van der Waals layers and that it is readily solved in the presence of a perpendicular external magnetic field, allowing the magnetic field dependence of the electronic polarization to be evaluated explicitly.
The influence of a static and uniform magnetic field that is parallel to the surface-normal of the thin film is accounted for layer-by-layer; a 2D EFA Hamiltonian that is minimally coupled to B is assumed to describe the electronic dynamics within an isolated layer, and is taken as Eq. ( 17) with ∆ → ∆ S and m → M lz for l z ∈ {1, . . ., N } a layer index.We then couple N copies (each indexed by an l z ) of this 2D isolated layer Hamiltonian via ∆ D as described above.We take M lz ̸ = 0 only if l z = 1 or N to account for the presence of magnetic dopants only in the outermost top and bottom layer.As in Sec.III A, we take the orientation of the magnetization describing the Here the exchange splitting m of the Dirac cones associated with the top and bottom layers (resulting from the presence of surface magnetic dopants) is taken to be 0.1∆S and taken to vanish for all other layers.(b) Distribution of the electronic charge density along the stacking direction for the same values of ∆D used in (a), with the film thickness taken to be 20 layers and z = 0 taken to coincide with the middle of the center layer.The dotted (solid) lines identify the contribution of the n = 1 (n = 0) Landau levels to this distribution, which is (anti-)symmetric in z.(c) αme vs. film thickness for various exchange splittings m, with ∆D taken as 1.5∆S;(d) Distribution of the electronic charge density in the stacking direction for the same values of m used in (c) for a 20-layer thin film; (e) αme vs. exchange splitting m due to magnetic surface dopants for films of varying thickness.
configuration of dopants at the top and bottom surface to be opposite (M 1 = −M N ≡ m) [76].
This model can be solved numerically for a varying number N ≥ 2 of layers, which we always take to be even, and a general energy eigenfunction is again of the form [77] Ψ r n,qy (x, y) = e iqyy Φ r n (x)/ L y , where for n = 0, and even (odd) values of j identify components that are associated with the bottom (top) surface of layer number j/2 ((j − 1)/2) enumerated in ascending order along the stacking axis z.We again endow the Ψ r n,qy (x, y) with z-dependence by taking C r n,(j,↑) → C r n,(j,↑) δ(z − z j ); choosing the origin to be at the center of the material, for even values of j we have z j = a(j−N )/2 and z j+1 = z j +a (see Fig. 1 of Lei et al. [75]).The electronic polarization (8) can then be written The results are shown in Fig. 3.In particular, as found in Sec.III A, if m = 0 then α me = 0 any N , whereas if m ̸ = 0 then α me approaches the expected quantized bulk value for N → ∞.That is, The quantization of α me in thick films can be understood by formally examining the model for an infinite number layers.In this case the analog of Eq. ( 17) is where ∆ kz ≡ ∆ S + e ikza ∆ D .To obtain Eq. ( 26) we have taken the basis states ū(α,↑),(kx,ky) → ū(α,↑),(kx,ky,R) for R ∈ aZ as hybrid WFs [17] (in R 3 ) that are spatially localized in z about the surface α in the unit cell at Rz.In the non-magnetic bulk there is one layer per unit cell and therefore α ∈ {0, 1}.Again projecting into the n = 0 subspace, Eq. ( 26) acts on C r 0,(b,↑) , C r 0,(t,↑) via which is a Su-Schrieffer-Heeger (SSH) model with ∆ S and ∆ D playing the role of the hopping parameters.This model has particle-hole symmetry: and therefore at half-filling supports topologically distinct ground states characterized by a Z 2 -valued invariant [78], which happens to be obtained by performing the BZ 1D -integral in the expression for the 1D bulk electronic ground state polarization [30,55].Indeed, the n = 0 contribution to the bulk polarization P z of Eq. ( 26) (which is related to adiabatically induced electronic currents in z [79]) is proportional to that 1D polarization and, in-line with the well-known results of the usual SSH model, at half-filling we find In contrast, the projected Hamiltonian in the n > 0 subspace is symmetric under 1D center-of-inversion (about the layer center) times spin-flip transformation and therefore does not support an electronic polarization [77].

IV. SEMI-ANALYTIC CALCULATION OF α CS IN A 3D TIGHT-BINDING MODEL
In this section our primary aim is to calculate the Chern-Simons coefficient α CS in bulk 3D non-magnetic V 2 VI 3 -type insulators.Effective models for the electronic states near the Fermi energy in these materials have been developed [65,74,75], but since effective models are not generally defined over the entire BZ, they are inadequate for this purpose.The reason for this can be understood by noting that the topology of the vector bundle of occupied Bloch states over the Brillouin zone torus constructed for the true Hamiltonian H(x, p(x)) of a crystalline insulator manifests as an obstruction to the existence of a smooth global gauge thereof [16,80], but smooth local gauges always exist [81].In particular, as explained previously, the linear response expression for α CS as a BZ-integral of the Chern-Simons 3-form is valid only in a smooth global gauge of that bundle.We emphasize this issue in Appendix A, where we explicitly show that in a low-energy effective model an incorrect half quantization of α CS can occur, which is reminiscent of the situation that arises when attempting to calculate the first Chern invariant of a 2D k • p model.Thus, lattice regularization of the previously developed effective models is required, and that is where we begin.

A. Construction of a 3D tight-binding model
In this subsection we construct a lattice regularization of a previously presented [74] effective model for the lowenergy electronic states in bulk 3D V 2 VI 3 -type insulators.Since a formally similar effective model describes anti-ferromagnetic Z 2 TIs including Mn(Sb x Bi 1−x ) 2 X 4 (where X = Se, Te) [75], which will be the focus of future work, we in fact construct a sufficiently general lattice regularization scheme that can be applied to both families of materials.Notably, both families of materials can be viewed as layered compounds, with each layer having 3-fold rotational symmetry about the stacking (z) axis and 2-fold rotational symmetry about an axis perpendicular to the stacking axis [65,82].The bulk crystal structure of the non-magnetic (magnetic) family of materials consists of stacked 5-atom (7-atom) layers, called quintuple (septuple) layers, and has a center-of-inversion symmetry within each such layer [65,82].That layered structure motivated the development of low-energy effective models [74,75] in which discrete 2D k • p continuum models in the planes perpendicular to the stacking axis are coupled to one another with interlayer hopping parameters ∆ D to yield a 1D tight-binding model along the stacking axis for each 2D momentum (see Appendix A); these models can be thought of as the 3D bulk limit of the quasi-2D Hamiltonian employed in Sec.III B. Since the electronic states within each layer are Dirac-like, effective Hamiltonians of this type are called coupled-Dirac cone models.We construct 3D tight-binding models to describe both families of materials by first developing a 2D lattice regularization of the k • p model of an isolated layer (obtained by taking ∆ D = 0), then re-introduce ∆ D as in those past works.

A lattice regularized description of an isolated 2D layer
Two proposed low-energy effective models of 3D V 2 VI 3 -type band insulators include the above described coupled-Dirac cone model [74,75] and a traditional 3D k•p model [65].When restricted to a single layer (i.e. the k x -k y plane), the 2D k • p models that result from these descriptions are unitarily equivalent to first order in k x and k y (when certain parameter values are taken to coincide), as one would expect.Furthermore, it has been observed [83] that the restriction of the 3D k • p model to the k x -k y plane is unitarily equivalent to the Bernevig-Hughes-Zhang (BHZ) model [84] for certain parameter values.Thus, when the coupled-Dirac cone model is employed for an isolated layer of a non-magnetic material (i.e.∆ D = 0 and M lz = 0), it is unitarily equivalent to the BHZ model to first order in k x and k y .In fact, in this limit of the coupled-Dirac cone model, exact unitary equivalence with the BHZ model occurs if we generalize ∆ S → ∆ S − B(k 2 x + k 2 y ) and if ∆ S equals M of BHZ.Rather than generalizing ∆ S → ∆ S − B(k 2 x + k 2 y ), we could instead obtain exact unitary equivalence by restricting the BHZ model to first order in k x and k y by taking B = 0.However, doing so would result in the lattice regularization that we employ below to admit band gap minima at points in the BZ 2D in addition to that at (k x , k y ) = (0, 0), in contrast to the known monolayer band structure [85].Indeed, it will turn out that the parameter regime of immediate interest is B/∆ S < 0.
In their seminal work on 2D Z 2 TIs [84], BHZ present a square lattice regularization of their k • p Hamiltonian.Thus, we seek a square lattice regularization of Eq. (1) of Lei et al. [75] (generalized by taking ∆ S → ∆ S − B(k 2 x + k 2 y )) that, when restricted to the nonmagnetic case and applied to a single layer, is unitarily equivalent to Eq. ( 5) of BHZ [84].In doing so, we generally consider Hamiltonian operators of the form where are tuples of fermionic operators that act in the electronic Fock space.Products in Eq. ( 29) are the usual matrix multiplication.The fermionic operators involved in Eq. ( 29) are of the form ĉI,k , ĉ † I,k , for I a general Wannier type label.We require those operators to be such that the following is satisfied: ψI,k ≡ ĉ † I,k |vac⟩ are Bloch-type vectors [86] that are smooth over BZ d , orthonormal such that ψI,k ψJ,k ′ = δ I,J δ(k −k ′ ), and constitute a basis of the "relevant" electronic Hilbert space [87].The penultimate and final criteria imply that {ĉ I,k , ĉ † J,k } = δ I,J δ(k − k ′ ), which follows from the anti-commutation of the electron field operators themselves.These are the usual criteria that are employed when constructing a tight-binding model and, for a given crystal, there are many sets of vectors ψI,k (and therefore operators ĉI,k , ĉ † I,k ) that satisfy them.An equivalent description can be obtained by noting that a set of Bloch-type vectors ψI,k that satisfy the above criteria bijectively maps to a set of exponentially localized Wannier functions (WFs) [17,88,89], where R ∈ Γ H .We emphasize that in order for the WFs W I,R (r) ≡ ⟨r|W I,R ⟩ to be well-localized, the ψI,k must be smooth over BZ d [17].Typically it is with respect to the fermionic operators ĉI,R , ĉ † I,R defined by |W I,R ⟩ ≡ ĉ † I,R |vac⟩ that a tight-binding Hamiltonian is specified.
In-line with typical tight-binding approaches, we further assume the existence of a Wannier type label I of the form I = (α, σ), where α ∈ {0, 1} is an orbital index that identifies WFs of distinct spatial distribution and σ ∈ {↑, ↓} identifies the s z spin component [90].To coincide with those coupled-Dirac cone effective models [74,75] for which we develop this lattice regularization, we take α = 0 (1) to label states that are associated with the bottom (top) surface of a layer; this identification will not explicitly enter calculations in this paper, but will provide physical motivation for the terms that appear in the 3D Hamiltonian of the following subsection.We further discuss (see preceding footnote) the topological implications related to the assumed existence of a Wannier type label of the form I = (α, σ) in the following subsections, where we also address additional ubiquitous constraints that are imposed on the WFs (30).
A 2D low-energy effective model for the electronic states in a single layer of a MnV 2 VI 4 magnetic (V 2 VI 3 non-magnetic) van der Waals semiconductor is obtained by taking ∆ D = 0 (and M lz = 0) in Eq. (1) of Lei et al. [75].A lattice regularization thereof is motivated by BHZ [84], and is specified by taking d = 2 and the general H (d) (k) of ( 29) to be where A is related to ℏv D and fixed later by fitting the 2D band structure to the k • p dispersion near Γ and where sk x ≡ sin(k x ), ck x ≡ cos(k x ), etc.We have also introduced a layer index l z ∈ Z in anticipation of the generalization to 3D in the following subsection.This is a square lattice regularization since ), thus Γ 2D = span Z ({x, y}); the 2D lattice constant is taken to unity, and k x , k y are taken to be dimensionless.Of course, in the physical materials under consideration Γ 2D is a triangular lattice [65,82].Thus, there likely exists a lattice regularization that better captures the lattice scale physics, but for our purposes a square lattice regularization is sufficient.
In Eq. ( 31) the terms involving M lz act to generate an exchange mass within the layer l z in the magnetic case, and arise from the exchange interactions of the dynamical electronic degrees of freedom deemed relevant with those well below the Fermi energy and thus approximated as static.We assume that this interaction is approximated by a Heisenberg interaction between dynamic spins and the quenched magnetic moments.We assume the direction of the static magnetization is parallel to the stacking axis (the z direction) [75].The terms involving A and ∆ S can be understood as arising from spin-preserving electronic transitions between some relevant WFs of the layer that are associated with its top and bottom surfaces.
The eigenvalues of the single-layer Hamiltonian (31) demonstrate that it is indeed unitarily equivalent to that of BHZ [84] when our M lz = 0, their C = D = 0, and our ∆ S equals their M .If M lz = 0, then we similarly find that all of the energy bands of (31) are degenerate at Γ when ∆ S = 0. Similar degeneracies appear at the other high-symmetry points (k x , k y ) = (π, 0) and (0, π) when ∆ S − 4B = 0, and at (π, π) when ∆ S − 8B = 0. Indeed it was previously found [75] that the materials of immediate interest, for which the band gap at half filling is about Γ, are described by ∆ S > 0 and thus B < 0. At half filling, BHZ find that the zero-temperature ground state of their model is insulating and in a Z 2 -even (odd) phase for M < 0 (0 < M/2B < 2) [84]; varying M/2B across 0 (2) closes and re-opens the band gap(s) at Γ (at (π, 0) and (0, π)), and a topological phase transition occurs.However, this 2D topology is not directly related to that captured by α CS in 3D, which will vanish in the limit of isolated layers independent of the parameters of a single-layer, even though the symmetries on which the classification is based are the same.Our focus in this paper is on the properties of the 3D crystals formed when these 2D layers are stacked and their states hybridized, not on the properties of individual layers.

A non-magnetic 3D multi-layer Hamiltonian
We consider 3D models in which the electronic states associated with adjacent surfaces of nearest-neighbor layers are hybridized by a phenomenological interlayer hopping parameter ∆ D : where we have assumed that all layers are identical apart from a layer-dependent M lz , that each layer is described by the 2D Hamiltonian (31), and that each pair of nearest-neighbor layers are separated by the same distance a in the z direction.In writing Eq. ( 33) we have assumed a non-magnetic or ferromagnetic configuration of the static magnetic moments, in which case ∆ D is always identified with a hopping parameter between WFs associated with different unit cells.For general magnetic configurations there are two issues that need be addressed related to the precise definition of the WFs implicit in defining fermionic operators ĉ(α,σ),R+lzaz , ĉ † (α,σ),R+lzaz that involve the layer label l z .First, depending on {M lz } lz∈Z , the group of translations Γ H under which a single-particle 3D Bloch Hamiltonian H(r, p(r)) -from which the tight-binding model might be obtained -is invariant may not be equal to the group of translations Γ 3D ≡ Γ 2D × aZ under which the multi-layer crystal lattice is invariant (Γ H ⊆ Γ 3D ).Then, although R + l z az ∈ Γ 3D for all l z ∈ Z, generically R + l z az / ∈ Γ H and therefore cannot be used to label WFs or their corresponding operators.In the non-magnetic (M lz = 0) and ferromagnetic (M lz = M l ′ z ) cases, the magnetic and chemical unit cells coincide, Γ H = Γ 3D , and R + l z az ∈ Γ H for all l z ∈ Z.However, in the anti-ferromagnetic case (M lz = −M lz±1 ) this is not so and, strictly speak-ing, Eq. ( 33) requires modification.For completeness, we present the tight-binding Hamiltonian for the antiferromagnetic case in Appendix B and focus on the nonmagnetic case in the remainder of this paper.Second, a set of 3D WFs of a multi-layer crystal (whose corresponding operators appear in Eq. ( 33)) is not generally equal to the set of all l z a 3 -translated 2D WFs of an individual layer (whose corresponding operators appear in Eq. ( 31)) thereof.In principle, to construct an accurate 3D tightbinding model one indeed requires WFs of the 3D Bloch Hamiltonian H(r, p(r)).However, such details will not be relevant in our calculation of α CS as we will not employ actual WFs but rather make simplifying approximations regarding their form.Consequently, this second type of imprecision in writing Eq. ( 33) is inconsequential to our calculation.Meanwhile the line of reasoning leading to Eq. ( 33) yields a physically well-motivated form of Ĥ(3D) reg .Focusing on the non-magnetic case (M lz = 0) such that Eq. ( 33) is valid, we use Eq.(30,31) and for d = 3 and k, k ′ ∈ BZ 3D to recast Eq. ( 33) as a BZ 3Dintegral of the form (29).In this case the general H (d) (k) is taken to be and ∆ k ≡ ∆ S (k x , k y ) + e iakz ∆ D .The eigenvalues of Eq. ( 35) are Fig.4), where Crucially, in this model there is a two-fold degeneracy at each k ∈ BZ 3D , which is a consequence [91] of a centerof-inversion symmetry and a (fermionic) time-reversal symmetry of (35); we explicitly demonstrate these sym-metries in Appendix C. Consequently, eigenvectors of ( 35) are highly non-unique, in a more general sense than the usual k-dependent phase ambiguity associated with Bloch's theorem.However, as described in Appendix A, the coupled-Dirac cone effective model that motivated (35) has a natural set of eigenvectors, and under the substitutions ℏv D k a → Ask a and ∆ S → ∆ S (k x , k y ), those eigenvectors are related to a set of orthogonal eigenvectors of (35), written here in the basis of Bloch-type vectors ψ(α,σ),k that are assumed smooth over BZ 3D (recall the discussion below Eq. ( 29)).The energy eigenvectors in (37) are of Bloch's form, satisfying |ψ n,k ⟩ = |ψ n,k+G ⟩ for any G ∈ Γ * 3D = span Z ({2πx, 2πy, 2π a z}).When ε(k) > 0 for all k ∈ BZ 3D , i.e. when the model is a band insulator at half filling (the case of primary interest), these |ψ n,k ⟩ are smooth over BZ 3D .We can also identify from Eqs. (37) a Γ * 3D -periodic unitary matrix T (k) relating Wannier and energy Bloch-type vectors, Since both ψ(α,σ),k (r) ≡ r ψ(α,σ),k and ψ n,k (r) ≡ ⟨r|ψ n,k ⟩ are of Bloch's form, they have associated Γ 3D -periodic functions ū(α,σ),k (r) ≡ r ū(α,σ),k = (2π) 3/2 e −ik•r ψ(α,σ),k (r) and u n,k (r) ≡ ⟨r|u n,k ⟩ = (2π) 3/2 e −ik•r ψ n,k (r), respectively, which are similarly related to one another by that T (k) and are smooth (modulo exp(−iG • r) for all G ∈ Γ * H ) over BZ 3D .The Bloch energy eigenvectors could alternately be chosen to simultaneously diagonalize the symmetry operators and the Hamiltonian itself.However, topological considerations forbid one from making that choice if the goal is to compute α CS in a gauge defined by those energy eigenvectors, as we now describe.
Many of the surprising single-particle properties of the various types of topological insulators can be understood by studying the structure of a particular (complex) vector bundle [16,80,92], the Hermitian bundle of occupied Bloch states over the Brillouin zone torus denoted V π V − − → BZ d .Band insulators are special in that the Hilbert space of occupied Bloch states at k ∈ BZ d , If there are N fully occupied energy bands then V k ∼ = C N for all k ∈ BZ d .The essential mathematical object is (isomorphic to) k∈BZ d {k} × V k , which is the smooth manifold that is obtained by attaching to each k ∈ BZ d the corresponding V k and equipping this set with a topology and smooth structure such that the natural projection map (k, Although this construction is intuitive, it is convenient to instead consider an isomorphic bundle over BZ d [94] with total space defined by where H is the Hilbert space of Γ-periodic functions, P V (κ) = N n=1 |u nκ ⟩ ⊗ ⟨u nκ | is the projector onto the space of Γ-periodic parts of occupied Bloch functions associated with κ ∈ R d , and the equivalence class It has been proved [80,92] that in any band insulator V π V − − → BZ d is a vector bundle (in particular, a Hermitian bundle) and is therefore locally trivial: for every k ∈ BZ d there exists an open neighborhood U ⊆ BZ d of k over which (π −1 V (U ) In general, V need not be globally trivial; the above is not necessarily satisfied for U = BZ d .However, it has been shown [80,88,92] that time-reversal symmetry implies (V In physics, the topology of V often manifests through the possible frames (also called gauge choices) thereof.A local frame of V over an open subset U ⊆ BZ d can be obtained from a collection of maps ũi in U ⊂ R d (the complete set of representatives of U that is contained in the Wigner-Seitz cell of Γ * ) defined by ũi (κ) ≡ (κ, |ũ i,κ ⟩) that satisfies ∀κ ∈ U : (|ũ i,κ ⟩) i∈{1,...,N } is a basis of RanP V (κ) (i.e. the related Bloch-type functions | ψi,k ⟩ are a basis of V k for k ∈ [κ]); if the |ũ i,κ ⟩ are energy eigenvectors then the gauge is called Hamiltonian.A general result [96] is that a vector bundle E π − → M is locally trivial over an open subset U ⊆ M if and only if there exists a smooth frame of π −1 (U ) π − → U .Then, about every k ∈ BZ d a smooth local frame of V always exists and if there is TRS then a smooth global frame of V exists.For example, if a 2D insulator is characterized by a nonzero first Chern invariant C, then V is not globally trivial; nonzero C acts as an obstruction to the existence of a smooth global gauge [80].In this case, that V is not globally trivial can be understood to manifest in calculations though the fact that the Berry connection cannot be made smooth over BZ 2D and the integral expression for C then returns a nonzero value.The Z 2 invariant relevant in this work is more subtle.TRS implies that V is globally trivial, but it has been shown [16,89,92,97] that a non-trivial Z 2 invariant acts as an obstruction to the existence of a global gauge that is both smooth and timereversal-symmetric.Thus, there need not exist a smooth global Hamiltonian frame of V [98].And since α CS must be calculated with respect to a smooth global gauge of V, in a Z 2 -odd insulator the components of that gauge are topologically forbidden to be energy eigenvectors that simultaneously diagonalize the time-reversal operator.Indeed, the generic nonexistence of a suitable Hamiltonian gauge means that to compute α CS requires [11] the employment of a Wannierization-like process [17] in order to construct Bloch-type functions | ψi,k ⟩ that constitute a smooth global frame of V. Within the model we employ, at half-filling and in the case of a band insulator, timereversal and inversion symmetry together imply that the two occupied Bloch states are globally degenerate.Any frame associated with these bands is therefore Hamiltonian; every smooth global gauge of this V is Hamiltonian.

Band Number Truncation in Tight-Binding Models and Crystal Momentum Dependence of Basis States
Before we can employ this model to compute α CS , there is one more issue to address.Like most tightbinding models that appear in the literature, the Hamiltonian operator identified by Eq. ( 31) and ( 33) is not fully defined since the Bloch-type vectors ψ(α,σ),k that correspond with the operators ĉ(α,σ),k , ĉ † (α,σ),k are not completely specified [99].This can be particularly problematic when calculating, for example, band-diagonal components of the Berry connection, in which case one must explicitly take k-derivatives of the cell-periodic part of Bloch energy eigenvectors.(In k • p models this issue is avoided by construction, since these models are formulated with respect to a constant frame of the space of relevant Bloch states over the subset of BZ d for which the model is accurate [100, 101].)In analytical calculations using tight-binding models it is almost always [102] assumed that the WFs (here the W (α,σ),R ) with respect to which the model is specified are atomic-like, with negligible overlap between WFs that are associated with different unit cells.In this case, the variation of ū(α,σ),k with k is negligible (within the above defined equivalence classification) allowing it to remain unspecified.In this work we assume that the Hamiltonian operator (33) is written in a basis of Bloch-type functions ψ(α,σ),k with corresponding ū(α,σ),k that are constant (modulo a phase) over BZ 3D , ∂ ∂k a ū(α,σ),k ≡ 0 for a = x, y, or z (≡ 1, 2, or 3) [103].While this may seem to be a drastic approximation, and in many cases it is indeed too dras-tic, we provide a physically motivated justification and argue that it is (at least) not mathematically inconsistent with the assumptions already inherent to our model.
A condition that one could potentially require of tightbinding models, in which the full electronic Hilbert space of a crystal is truncated to include only the states associated with a finite number of bands around the Fermi energy, is that the k-dependence of the Bloch states be accurately captured throughout the Brillouin zone.The local k-dependence of the |u n,k ⟩ around some k 0 can be calculated using k • p perturbation theory, which yields where m e is the electron mass.It follows from Eq. ( 39) that a good criteria for the reliability of a tight-binding model is that all momentum matrix elements between the |u m,k0 ⟩ corresponding with included bands and neglected bands are small at all k 0 of interest.When such a set of bands that are dissociated in this sense across the entire Brillouin zone can be identified [104], we say that they satisfy a global isolation condition.When this global isolation condition is satisfied, we can always choose the finite-dimensional basis vectors of the included bands to be independent of k, for example they could be the eigenstates at one particular k 0 .Our expectation is that global isolation holds only when it is implied by chemistry, i.e. the set of energy bands derives from a linear combination of atomic or molecular orbitals.This global isolation condition is probably difficult to satisfy in practice and probably unnecessary for many physically relevant calculations, but it greatly simplifies calculations of Berry connections and related quantities.When we start from a phenomenological tight-binding model, as in this calculation, we know nothing about momentum matrix elements between included and neglected bands, so we have little choice but to assume the global isolation condition.It is probably common in ab initio DFT-derived models that global isolation is not satisfied.In particular, for Z 2 TIs, in which level inversion at one point in the Brillouin zone plays the essential physical role, it will never be possible for the set of occupied bands to satisfy the global isolation condition.
To consider the mathematical implications of taking ∂ ∂k a ū(α,σ),k ≡ 0 over BZ 3D , we return to the bundletheoretic framework.In this paper we always consider the model at half-filling.Thus, at any k ∈ BZ 3D , the ū(α,σ),k ( ψ(α,σ),k ) need not be contained in RanP V (k) (V k ).Therefore, assuming the global isolation condition for (the space spanned by) the ū(α,σ),k does not constrain the topology of V π V − − → BZ 3D .However, this does have implications for the topology of the vector bundle of all relevant Bloch states, which we denote The total space B of that bundle is constructed in a manner similar to V, but now with the fiber at each k ∈ BZ 3D isomorphic to bundle follows from an analogous argument [80,92] to V π V − − → BZ 3D , which can indeed be applied to any set of isolated energy bands.)In constructing the model, we have already explicitly assumed, as one always does in employing tight-binding models, that B π B − − → BZ 3D has a smooth global frame, and is therefore globally trivial.
The assumption that we can label the WFs (30) by I = (α, σ) implies the existence of a smooth global frame of the form ( ū(α,σ),k ) α∈{0,1},σ∈{↑,↓} .In particular, this means that B π B − − → BZ 3D can be decomposed as a product of vector bundles that correspond with spin-up (σ =↑) and spin-down (σ =↓) sectors [105], and that each sector admits a smooth global frame ( ū(α,σ),k ) α∈{0,1} for a given σ.Thus the triple of Chern numbers (since d = 3) that characterizes each sector vanishes and [106,107].Therefore, the topology of B π B − − → BZ 3D does not forbid the existence of a smooth and symmetric global frame.In addition, if we assume the global isolation condition then for any k, k ′ ∈ BZ 3D we have RanP B (k) = RanP B (k ′ ).Then, while we do not make a general existence argument for a frame with components ū(α,σ),k that are constant over BZ 3D , assuming one to exist is (at least) not generically forbidden by the topology of This argument applies to any tight-binding model that is specified with respect to WFs whose type labels are assumed to involve a spin index and satisfies the global isolation condition.

B. Calculation of α CS
As mentioned above, in order to compute α CS as the BZ 3D -integral of the Chern-Simons 3-form, which is valid strictly for band insulators with corresponding V π V − − → BZ 3D that is globally trivial, one is required to work in a smooth global gauge of V. Since Bloch energy eigenvectors are not generically smooth over BZ 3D , particularly in Z 2 TIs, a Hamiltonian gauge choice is not typically viable.In past works [11,15,37,45], a gauge of that type is induced by the choice of WFs, where WFs are there assumed to be constructed from the (un)occupied energy eigenvectors alone.However, the WFs employed to construct tight-binding models are not typically of this type; the gauge ( ū(α,σ),k ) α∈{0,1},σ∈{↑,↓} is not valid for this calculation.
As described in Sec.IV A 2, the model we employ has the positive feature that we are guaranteed a smooth global Hamiltonian gauge due to the combination of TRS and each set of isolated bands being completely degenerate over the entire BZ 3D .Using the |u n,k ⟩ related to (37), we can define a suitable frame u of B π B − − → BZ 3D pointwise by u k ≡ (|u n,k ⟩) n∈{1,2,3,4} ; although u is technically a frame of B, since it is Hamiltonian it projects to a frame of V.This gauge choice corresponds to taking the bulk electronic polarization and orbital magnetization to be defined with respect to the set of WFs |W n,R ⟩ = It has been shown [10,11,15,37,45] that, with respect to any smooth global Hamiltonian gauge u defined by where the sums are over the initially occupied band indices (here v, v ′ , v 1 ∈ {1, 2}) and ξ a nm are the components of the non-Abelian Berry connection induced by u, Although Eq. ( 40) is gauge dependent, transforming from one appropriate gauge of V to another, be it a Hamiltonian gauge or otherwise, can only change its value by an integer multiple of e 2 /hc [10,11]; that is, there is a quantum of indeterminacy associated with α CS .We now calculate α u CS via Eq.( 40).Assuming that ∂ ∂k a ū(α,σ),k ≡ 0 (see Sec. IV A 3), Eq. ( 41) becomes where we have used (38) for α ∈ {1, 2} and σ ∈ {↑, ↓}) and defined the Hermitian matrix T a populated by elements Using Eqs.(37,38) in Eq. ( 43), we obtain All other components are related to those given above.Note in particular that T a is Hermitian, T a nm = (T a mn ) * , and that by explicit calculation it can be shown The relations (45) need not be satisfied for a different choice of energy eigenvectors.Using Eqs.(42,44,45), the integrand of the first term of Eq. ( 40) is here found to be The integrand of the second term of Eq. ( 40) is (in this case) −1/3 of the first.Indeed, applying the general relation ϵ abc ∂ b ξ c nm = iϵ abc r ξ b nr ξ c rm (we now adopt the shorthand ∂ a ≡ ∂/∂k a ) to this case (such that r ∈ {1, 2, 3, 4}) in combination with the relations (45), one can explicitly show that With this, we obtain the final expression We have noted that α u CS is invariant under a change in energy scale and used this to remove the explicit depen- Representative set of values of α u CS obtained from numerical integration of Eq. ( 47) for (arbitrarily chosen) ∆S = 190 meV.From these data the topological phase diagram (Fig. 5) can be deduced.When the values of the crystal parameters are varied such that the band gap does not close, the value to which α u CS evaluates is unchanged.5. Phase diagram plotting the Z2 topological classification of the band insulating electronic ground state of Ĥ(3D) reg (specified in Eq. ( 35) and Sec.IV A 3) at half-filling.The darkest (vertical) dashed lines identify the points in parameter space at which a band gap closing occurs due to one of Eq. (48a,49a) evaluating to zero, the medium dashed lines identify points at which one of Eq. (48b,49b) evaluates to zero, and the lightest dashed lines identify points at which one of Eq. (48c,49c) evaluates to zero.dence on A in Eq. ( 47) through the introduction of scaled parameters B ′ ≡ B/A, ∆ ′ D ≡ ∆ D /A and ∆ ′ S ≡ ∆ S /A, and ε ′ (k) ≡ ε(k)/A.In Table I we list values of α u CS obtained by extrapolating numerical estimates of the integrals over BZ 3D in Eq. ( 47) to convergence for various combinations of model parameters; we were not able to perform an analytical integration.We are careful to avoid sets of parameters for which the band gap (at half-filling) vanishes since Eq. ( 40) applies only to band insulators.
In Sec.III, the parameter B did not appear.This can be understood by recalling that the tight-binding model we employ was developed as a lattice regularization of an effective model that assumes that the low-energy states of the materials of interest are only near (k x , k y ) = (0, 0); this is a property of our lattice model only when B/∆ S < 0 and |∆ S | ≈ |∆ D |.In that case our expectation is that the largest contributions to the integrand of (47) come from near the line (k x , k y ) = (0, 0), since ε(k) is smallest there.Expanding our expression (47) for α u CS to first order in k x and k y about (0, 0), we find that the terms involving B cancel out and that where we have artificially extended the domain of integration for k x , k y from the subset of BZ 2D near the expansion line to R 2 .For B/∆ S < 0 and |∆ D | ≈ |∆ S |, Eq. ( 50) is consistent with Fig. 5.A similar approximate expression for α CS was previously derived by Rosenberg and Franz [108] for models hosting Dirac-like low-energy states.We caution, however, that topological index calculations like this one, which focus only on contributions from regions near certain lines or points in k-space, can fail.An explicit example is provided in Appendix A.

V. ATOMIC-LIKE AND ITINERANT CONTRIBUTIONS TO α CS
In Sec.III we employed a coupled-Dirac cone model of a multi-layer thin film -which had a vanishing magnetic exchange mass (i.e.no magnetic dopants) in the interior layers, but allowed a finite exchange mass in outermost layers -and found (see Fig. 3) that a perpendicular magnetic field can induce a charge density only near the surfaces, that this occurs only in the presence of local magnetic dopants, and that the corresponding magnetoelectric coefficient is quantized in sufficiently thick films only when the dopant configuration in the top and bottom layers is opposite.In Sec.II we argued from physical grounds that a meaningful bulk P (B) , particularly in non-magnetic Z 2 TIs, must (at least) involve an implicit consideration of surfaces to account for these magneticdopant requirements.We now ask whether there is a purely bulk manifestation of those requirements.
To begin, we recall the equivalence of the susceptibilities in Eq. ( 6), which is at least valid in bulk insulators for which V π V − − → BZ 3D is globally trivial.We focus on the macroscopic electronic orbital magnetization M and recall from the modern theory [31,32] (or from other formulations of the microscopic response theory [34]) that where M is called the atomic-like contribution and M the itinerant contribution.(In the modern theory these are typically denoted M LC and M IC , respectively.)In an unperturbed bulk insulator, or when the orbital magnetization is induced by a uniform electric and/or magnetic Maxwell field, M is uniform and static, as are M and M [38].In the modern theory it has been argued [31] that in ferromagnetic (M (0) ̸ = 0) bulk insulators with isolated energy bands, interior and surface currents in finite samples thereof can be used to uniquely determine M (0) and M (0) , respectively [109]; we use the superscript (0) to label spontaneous orbital magnetizations and later use (E) to label magnetization that arises in linear response to the electric field.However, when there are degeneracies within the band structure, the decomposition of M (0) as a sum of M (0) and M (0) is gauge dependent [32].In this case, interior and surface currents do not uniquely determine M (0) and M (0) ; additional data related to the complete set of WFs being employed is required.Still, the partition into contributions associated with interior and surface currents may retain meaning.To investigate this possibility, we explicitly study the expressions for M (E) and M (E) , which, like their unperturbed counterparts, are individually gauge dependent, but their sum M (E) is unique only modulo e 2 E/hc [10,11].
In the approach of the modern theory, the expressions for the atomic-like M (E) and itinerant M (E) contributions to M (E) (Eq.( 6) of Malashevich et al. [45]) are given in terms of states that implicitly involve the electric field E. This is inconvenient for our purposes.As mentioned in Sec.II, there exists a microscopic approach [34], motivated by that of PZW, that identifies atomiclike and itinerant contributions to M by different means, but yields the same expressions in an unperturbed bulk insulator [37].Indeed the microscopic approach reproduces the magnetoelectric susceptibility (9) derived via the modern theory.At linear response to a uniform dc electric field the expressions resulting from those distinct approaches are similar in that M (E) originates from the effect of E on the electronic dynamics directly (called dynamical modifications; see Sec. 3. C. 1 of Mahon and Sipe (MS) [37]), while M (E) originates from a combina-tion of dynamical modifications and changes in the form of the relevant matrix elements due to E (called compositional modifications).A familiar example of this is paramagnetic and diamagnetic response of an atom to an electromagnetic field [110], which can be understood as arising from dynamical and compositional modifications, respectively.Below we employ the expressions derived by MS [37].
Working in the smooth global Hamiltonian gauge u defined by the Bloch energy eigenvectors (37) and using the degeneracy of the energy bands (see Eq. ( 36)), Eq. ( 70) of MS for the atomic-like contribution evaluates to [111] where v, v ′ ∈ {1, 2} and c, c ′ ∈ {3, 4}.Here and elsewhere we often keep the k-dependence of quantities implicit.In this gauge the atomic-like magnetization vanishes as an immediate result of the relations (45) between the band components of ξ a .The itinerant contribution M i(E) has both a dynamical modification M i(E;I) and a compositional modification M i(E;II) (given by Eqs. ( 72) and ( 71) of MS, respectively).In any smooth global Hamiltonian gauge u ′ of our model, M i(E;I) In the particular Hamiltonian gauge u that we employ, Eq. ( 71) of MS evaluates to To reach the final equality in Eq. ( 53) (in our case, m ∈ {1, 2, 3, 4}) we have used (44,45) for ξ a and recognized the result as the analytic expression (47) for α u CS .We find, therefore, that in the gauge u the topological magnetoelectric response is entirely itinerant, i.e. that Now consider working in some other smooth global Hamiltonian gauge u ′ .(Recall the discussion at the end of Sec.IV A 2.) Then at each k ∈ BZ 3D there exists a unitary matrix U (k) relating the components |u ′ n,k ⟩ and |u n,k ⟩ of the gauges u ′ and u [112], where U c,v (k) ≡ 0. The components of the non-Abelian Berry connection induced by u ′ are and are related to those induced by u (denoted in Eq. ( 42) as ξ a nm ≡ (ξ u ) a nm ) by where Notably, U c,v (k) ≡ 0 implies W a cv (k) ≡ 0. Employing Eq. ( 56), we are able to relate the atomic-like and itinerant contributions to the orbital magnetization that are identified with respect to the distinct gauges u and u ′ .In particular, using Eq. ( 52) we find In general, the only restriction on U (k) is that it be smooth and Γ * H -periodic, thus there is no reason to ex-pect ᾱli u ′ to vanish.As an example, consider a gauge transformation U (k) such that U n,m (k) ̸ = 0 only if n = m = 1.Since U (k) is unitary, there exists a smooth function λ(k) such that ∀k ∈ BZ 3D : ∀G ∈ Γ * H : λ(k + G) = λ(k) + 2πj for j ∈ Z, and U 1,1 (k) = e −iλ(k) .Taking for example λ(k) ≡ sin(k x ) cos(ak z ) results in ᾱzx u ′ ̸ = 0 and taking λ(k) ≡ cos(k x ) sin(ak z ) results in ᾱzz u ′ ̸ = 0.Then, since TRS implies that in any smooth global Hamiltonian gauge u ′ we have The only feature that we might think is special about the non-magnetic case considered here is that there exists a smooth global Hamiltonian gauge u in which M i(E) In an related work [77] we investigate this further by considering the lattice regularization of the coupled-Dirac cone model for the MnV 2 VI 4 family of anti-ferromagnetic Z 2 TIs.In that work we show that there does not generically exist a smooth global Hamiltonian gauge in which the topological magnetoelectric response is entirely itinerant.Indeed whether such a gauge exists is related to the geometry of V π V − − → BZ 3D for the insulator under consideration -a smooth global gauge of a vector bundle is equivalent to a global trivialization.
Given that the topological magnetoelectric response in non-magnetic Z 2 TIs is (with important qualifications) itinerant, from these purely bulk considerations alone we might expect it to originate entirely from surface currents.In cases for which this is true, there is good reason to regard the identification of the physically meaningful bulk magnetoelectric response as α CS with suspicion, especially when we know that any finite size material sample will host gapless surface states whose role must also be considered.Surface magnetic dopants can generate gaps in the energy dispersion of the surface states and puts this identification on safer footing.To generate a gap everywhere on the surface, it is also required that the outward normal of the magnetization characterizing the dopant configuration never change sign.For a thin film this leads to the requirement that this magnetization retains a perpendicular component on the side walls and that it have opposite orientations on the top and bottom surfaces.Indeed these are the same magnetic-dopant requirements that were mentioned at the beginning of this section and it is only when they are satisfied that microscopic finite-size calculations agree with the simple interpretation of the bulk magnetoelectric linear response calculation.When these requirement are not satisfied, there is no physically meaningful magnetoelectric response in non-magnetic Z 2 TIs.

VI. SUMMARY AND DISCUSSION
In this paper we study the topological magnetoelectric effect in 3D non-magnetic Z 2 topological insulators.We begin Sec.II with a summary of the considerations that underlie notions of polarization and (orbital) magnetization in material media, and of magnetoelectric response.In doing so, we illustrate the precarious relationship between that response and time-reversal symmetry, in nonmagnetic Z 2 TIs in particular.Focused on materials of that type, we present a slightly modified interpretation of the argument [15] used to derive an expression for the bulk topological magnetoelectric coefficient, one that is manifestly consistent with time-reversal symmetry.We conclude that the polarization adjacent to a particular surface position is activated locally (see below) by timereversal symmetry breaking at that position.This interpretation is supported by explicit calculations in finite sized (Sec.III) and bulk (Sec.V) materials.In the process of formulating this interpretation, we clarify a partially inconsistent conclusion related to the adiabatic current expression from which the magnetoelectric response was previously [15] deduced.The central difference between our interpretation and the present one is that ours accounts for the fact that a quantized magnetoelectric effect occurs in non-magnetic Z 2 TIs only if static magnetic surface dopants are present and are configured such that their out-of-plane magnetization orientations at the top and bottom surfaces are nonzero and opposite.This resolves the tension between the seemingly contradictory roles played by TRS in finite sized and bulk insulators.
The model that we employ in Sec.III is for an infinite quasi-2D thin film with two perfectly flat surfaces on top and bottom.The explicit calculations on which our interpretation is based show that a magnetic field dependent charge density can arise only at those surfaces and depends on the surface normal of the magnetization on that surface.This conclusion can be generalized to arbitrary surfaces by noting that only the sign of the surface normal magnetization component is important, and that the non-locality length along lines where the surface normal changes sign [113], λ = ℏv D /m, is finite.As long as the typical length scales of surface regions with a fixed sign of the surface normal of the magnetization is large compared with λ, the properties we have calculated for uniform surfaces apply to arbitrary surface magnetization profiles, and the term local here refers to averages over finite regions of area λ 2 .One copy of the infinite flat thin film model in Sec.III can be associated with each locally flat spatial region of the sample.Within each such region, our interpretation of the adiabatic current calculation is that ∂P i /∂B i is constrained by bulk topology to have a quantized value.The value that is realized at a given position on a sample boundary is locally determined by the surface magnetization at the same position on the boundary [10].The magnitude of the polarization response to magnetic field tends to be the smallest allowed quantized value, so changes from position to position tend to be changes in sign only.Sign changes in the surface magnetization profile will produce locally ungapped regions on the surface that will support [113] chiral edge states that separate polarization domains.
As we have explained, the topological magnetoelec-tric effect (TME) in our picture is a local surface charge density property.If this is the meaning ascribed to the TME, then the quantum Hall effect that appears in the film as a whole when the magnetizations on top and bottom surfaces have the same orientation can be viewed to be a natural partner of the TME.This magnetic configuration implies opposite dependence of polarization on magnetic field on opposite surfaces of the sample, and via Eq.( 4) to a total charge density in the insulator that varies linearly with magnetic field.A dependence of charge density on magnetic field in a quasi-2D insulator is equivalent to the quantum Hall effect via the Streda formula.(The relationship between the quantum Hall effect in TI thin films observed optically and the TME has been controversial [114][115][116][117][118].) When the magnetizations are anti-parallel on the two surfaces of a thin film, the same argument implies that in our interpretation the field-dependent polarization inside the material is spatially uniform, so that the total charge density is independent of magnetic field.Thermodynamic identities nevertheless imply that the orbital magnetization must vary with the electric field applied across the system, implying that the orbital magnetization is sensitive to an energy shift between top and bottom surfaces in the local density-of-states even though no chiral edge states are present at the Fermi energy, as argued in Ref. [59].
In Sec.V we consider whether these surface magneticdopant qualifications on the magnetoelectric response of non-magnetic Z 2 TIs manifest in the bulk expressions.If they did, one might guess to use the equivalence of the susceptibility tensors (6) and investigate the properties of the atomic-like M and itinerant M contributions to the electric-field-induced magnetic dipole moment.Indeed M was originally associated with surface currents [31].The identification of M and M is generally gauge dependent, even in unperturbed insulators [32], thus associating any physical meaning with that partitioning is suspect.Nevertheless, with respect to a particular smooth global Hamiltonian gauge we analytically demonstrate that the topological response δ il α CS arises entirely through that of M (independent of a particular choice of parameter values, in all regions of parameter space for which there is always a band gap).However, this is not true in a generic Hamiltonian gauge.In a related work [77] we consider anti-ferromagnetic (AFM) Z 2 TIs that exhibit a generalized TRS.In that case we show that no gauge exists for which the response is entirely itinerant (in the above sense).In AFM materials, surface magnetic dopants are unnecessary for a magnetoelectric effect to manifest since the termination of the material already breaks the generalized TRS.Thus, it may be the case that the surface magnetic-dopant qualification in thin films is equivalent to the existence of a smooth global Hamiltonian gauge with respect to which the bulk response is entirely itinerant.Indeed, this may be understood as a type of microscopic bulk-boundary correspondence relating the bulk geometry of V π V − − → BZ 3D to the surface magnetic-dopant qualifications for the TME to manifest in a particular material.Although the nonexistence of certain smooth global gauges has been related [16] to the topology of V π V − − → BZ 3D , we do not suspect this to be the case here; we anticipate that in an AFM TI there could exist some smooth global gauges u for which there exists a finely tuned point in the parameter space where M (E) u equals α u CS E. Thus, whether such a gauge generically exists is taken to be a geometric property of V π V − − → BZ 3D .That magnetic surface dopants activate the topological magnetoelectric effect in thin films of non-magnetic Z 2 TIs can be understood on physical grounds, as we now describe.Assume that the magnetoelectric coefficient in a bulk non-magnetic Z 2 TI is given by α CS = (n + 1/2)e 2 /hc for some n ∈ Z.We seek the corresponding susceptibility tensor α finite (r) in a finite sample thereof, which we assume can be partitioned into contributions from the sample's interior and surface regions, such that α finite (r) = α CS f int (r) + α surf (r)f surf (r); {f int (r), f surf (r), f ext (r)} is a partition of unity over R 3 , where we take f int (r) to have support over the entire interior region of the sample and decay from 1 to 0 on some length scale near its boundary, f surf (r) to be nonvanishing only near the sample boundary, and f ext (r) to identify the region outside of the material.
If we assume that bulk physics always determines the properties of finite-sized samples such that α surf (r) can be neglected, then α finite (r) = α CS f int (r).We now show that this assumption leads to a contradiction.Implementing this assumption in Eq. ( 4), we find [119] ∂r a E b (r) + . . ., ϱ finite (r) = −α CS ∂f int (r) ∂r a B a (r) + . . . .
That is, although in bulk insulators that exhibit TRS the linearly induced charge and current densities are insensitive to α CS , physical consequence of α CS can manifest at the surface of a finite sample thereof.To reach this conclusion we assumed that J (1) finite (r) and ϱ (1) finite (r) are found by restricting the bulk P (1) (r) and M (1) (r) to a finite region of space, and not the bulk J (1) (r) and ϱ (1) (r) themselves.For if it were the latter, then α CS would again not contribute.This is strange; who decides which quantity should be restricted?It has also been pointed out [120] that the above expression for J (1) (r) corresponds to a half-quantized quantum anomalous Hall current at each surface.But how can this be if the sample exhibits TRS? Something has gone wrong.Indeed, the result of Sec.III was that a magnetic-field-induced polarization can only occur in non-magnetic TI thin films when TRS is broken by magnetic dopants at the surface, which is moreover consistent with the usual symmetry analysis.In particular, if m = 0 in Eq. ( 22) prior to taking the bulk limit, then α me = 0.
Thus, the assumption that α surf (r) can be neglected leads to contradictions.This illustrates that, while it may be true that the bulk α CS can manifest at the surface of finite size samples, one must also consider the response of the surface states that exist as a consequence of the bulk topology.This conclusion is consistent with topological field theoretic considerations, where it has been demonstrated that the surface states implied by a nontrivial 3D bulk Z 2 invariant (via the topological magnetoelectric coefficient α CS ) are related to the restoration of the parity anomaly and their response identically cancels that of the bulk [12][13][14].This contrasts the situation for the integer quantum Hall effect related to a non-trivial 2D bulk Chern invariant, in which case the bulk and surface responses add [121][122][123].
Ultimately, a ubiquitous aspect of theories considering polarization and orbital magnetization in crystalline solids has been the importance of interpretation [56].This work illustrates that in order to have a physically meaningful notion of bulk magnetoelectric response, in particular in insulators that exhibit TRS, the implicit involvement of a surface is unavoidable.Otherwise, as has been the case until now, the topological magnetoelectric response generally falls outside of the Peierls paradigm.Our hope is that this work will help to clarify misconceptions related to the topological magnetoelectric effect and thereby assist with its experimental demonstration.We expect that a quantized TME will be observable in non-magnetic TI thin films only when the surface normal of the magnetization does not change direction on either top or bottom surfaces and the magnetic gap induced at the surface is larger than the disorder potential on that surface.Changes in sign of the surface normal and reductions in local gap size on either surface will generically alter the global charge at which the densityof-states reaches a minimum and make it magnetic field dependent, swamping the topological magnetoelectric effect.It follows that the robustness against disorder that is critical to accurate observation of the integer quantum Hall effect, does not apply to the topological magnetoelectric effect.
For illustrative purposes, in this section we perform a too-simplistic calculation of a low-energy approximation of α CS using the bulk (i.e.periodic-in-z) version of the coupled-Dirac cone model introduced in Sec.III.In contrast to the quasi-2D models used in that section, this model is genuinely 3D.Notably, in the smooth local gauge that motivates the eigenvectors (37) of the tightbinding model obtained from a lattice regularization of this model, we will find the same result (50).However, it is generally the case that one does not produce the correct result for α CS using a low-energy effective model alone.Indeed if one somehow knows that the local gauge choice on the BZ 3D subset on which the effective model is valid can be extended to a smooth and periodic gauge over the entire BZ 3D of some lattice regularization thereof, then the correct quantization can result; this was used implicitly in past work [108].Otherwise, incorrect halfquantization of the susceptibility tensor can arise, similar to the situation when too-simplistic attempts are made to calculate the 2D Hall conductivity in a Chern insulator.
The effective low-energy Hamiltonian of 3D bulk nonmagnetic and magnetic multi-layer TIs that was previously studied by Lei et al. [75] is essentially a k • p model about Γ in the b 1 -b 2 (x-y) plane and a lattice model in b 3 (z); by construction, the model is valid only on a subset of the 3D Brillouin zone.The model is written with respect to operators that generate basis vectors v (α,σ),k that are products of 2D k • p basis vectors ū(α,σ),k0 , with the k • p expansion point k 0 taken to be (0, 0), and 1D Bloch-type functions ψ(α,σ),kz (z) that are associated with WFs W (α,σ),na3 (z) that are localized in a 3 ∥ z and taken to coincide with a bottom (top) layer "surface state" for α = even (odd).Then, in this model the operators generate vectors v (α,σ),k = v (α,σ),(0,0,kz) ≡ v (α,σ),kz , and, in fact, the possibility for variation in k z of the associated 1D cell-periodic functions ū(α,σ),kz (z) ∝ e −ikzz ψ(α,σ),kz (z) is neglected.In the non-magnetic case, α ∈ {0, 1} and the effective Hamiltonian is defined by its matrix representation where In this model, the condition of a band insulator at half filling is ∀k z ∈ [−π/a, π/a] : ∆ kz ̸ = 0. Thus, as the ratio |∆ D |/|∆ S | is varied, the band gap can vanish; in particular, if ∆ D = −∆ S then ϵ(0, 0, 0) = 0 and if ∆ D = ∆ S then ϵ(0, 0, π/a) = 0.These are the points in parameter space at which a topological phase transition may occur.
Due to the double degeneracy (A2), at every k for which the model applies, the set of eigenvectors is highly non-unique.However, (A1) has natural solutions in the ak ⊥ → 0, ∞ limits (k ⊥ ≡ k 2 x + k 2 y ), which motivates us to make a particular choice, In the ak ⊥ → 0 limit, the Hamiltonian (A1) involves only inter-layer transitions and the spin index becomes preserved, while in the ak ⊥ → ∞ limit the spin-orbit interaction dominates and the layer index becomes preserved; these limits are evident in the chosen form (A3). Notably, the set of vectors (A3) forms an orthonormal basis of the relevant electronic Hilbert space, which are smooth in k for an insulator (ϵ(k) > 0).From (A3) we can identify a unitary matrix T (k) analogous to that appearing in Eq. ( 38), here with elements Tn,(α,σ) (for a = x, y, or z) analogous to those appearing in Eq. ( 43).Unsurprisingly, we find that these matrices satisfy relations analogous to (45), but now we find Indeed, the matrix elements (A4) are consistent with those (44) that are obtained from the lattice regularized model; taking B = 0 and k x , k y to be small, Eq. ( 44) reduces to Eq. (A4).
In principle, this model, like any k • p description, is insufficient to compute α CS .That is, while low-energy models may accurately describe the physics of the states associated with those crystal momentum k near the k • p expansion point, the topology of V imposes constraints on global properties of Bloch energy eigenvectors and of smooth gauge choices over all of BZ d ; those constraints are generally not relevant on subsets of BZ d .Nevertheless, it is insightful to calculate a low-energy approximation of α CS to see where issues arise.We take the initial electronic state of the insulator to be the zero-temperature ground state and there to be two electrons per unit cell; that is, we take |ϕ 1,k ⟩ and |ϕ 2,k ⟩ to be initially occupied and |ϕ 3,k ⟩ and |ϕ 4,k ⟩ to be initially unoccupied.The components of the Berry connection induced by (A3) can be found via ξa nm (k) = − T a nm (k) (recall the discussion above Eq.(A1)).After some algebra we find To calculate quantities that are generally written as Brillouin zone integrals (such as α CS (40)) via a k • p description, one usually restricts the domain of integration to the subset of BZ 3D on which the model is valid (here to a limited range of k x , k y about (0, 0)) based on the assumption that outside of this region the integrand is negligible.In-line with this assumption, when k x and/or k y are large, the relevant integrand (A5) tends to 0; consequently, the domain of integration in the "k • p dimensions" is usually artificially extended to the entire plane and we do so here.With this, and moving to the cylindrical coordinates (k ⊥ , θ k ⊥ , k z ), we have = denotes an equality under the above described assumptions of the k • p approach.Evaluating the integral, yields the expected quantization of α CS .This is consistent with the result of the tight-binding calculation (see Eq. ( 50) and the surrounding text).
But what happens under a change of gauge?The result is that, in general, a change of gauge shifts the value of α CS by an amount different than n e 2 hc for n ∈ Z, in contradiction to the known behavior for a model defined over the entire Brillouin zone [10,11].This is an artifact of the use of a k • p description, which will not generally yield quantized topological invariants [124].In the case of α CS in particular, the issue is that one does not know whether or not the eigenvectors employed in the gauge choice on the BZ d subset on which the k • p model is valid extend smoothly and periodically over the BZ d of a lattice regularization thereof.For example, in the model studied here, consider a gauge constructed from v ′ n,k = 4 m=1 |v m,kz ⟩ U m,n (k), where This appears harmless within the context of the model (A1) for an insulator (∆ kz ̸ = 0).But, in the lattice regularization presented in the main text, ∆ kz → ∆ k and there may exist k ∈ BZ 3D : |∆ k | = 0 in a band insulator.Then attempting to define Bloch energy eigenvectors over the full Brillouin zone from the v ′ n,k (as was done in the main text (37) for the |v n,k ⟩) is pathological at any k ∈ BZ 3D : |∆ k | = 0; those vectors cannot be used to construct a smooth gauge over BZ 3D .Indeed a calculation of α CS in that gauge would be invalid, since a smooth global gauge was assumed in deriving the BZ 3D integral form of α CS .
There is another, more physical, complication in the case when M lz ̸ = 0 in 3D, which is the possibility of exchange coupling between the dynamical degrees of freedom of a given layer with the static magnetization of all other layers.If we assume that the most relevant exchange interactions for a given WF W (α,σ),R (r) is that with the static magnetization of the layer for which it is associated (via J S ) and that with the static magnetization of the layer with which it is nearest (via J D ) -for α = odd (even) this is the layer below (above) that for which W (α,σ),R (r) is associated -and we focus on the AFM case such that both orbitals within a given layer feel the same M lz±1 = −M lz , we can account for this by taking J S M lz → (J S − J D )M lz ≡ m lz in Eq. (31).Taking m lz = m for l z even and = −m for l z odd, writing (B2) as a BZ 3D integral via Eq.( 30), and using Eq. ( 34 The eigenvalues of (B3) are The double degeneracy of the energy bands at each k ∈ BZ follows from the combination of an inversion and a (fermionic) time-reversal symmetry (see Appendix C).Notice that taking m = 0 in (B4) does not reproduce (A2).This is not surprising because taking this limit implements an incorrect identification of the Bravais lattice for the nonmagnetic Hamiltonian; in deriving (B4) we explicitly use Γ H = Γ 2D ×2aZ rather than the correct Γ H = Γ 3D = Γ 2D ×aZ implemented in deriving (35) in the non-magnetic case.
It follows from Ŝ {ĉ A , ĉ † B } Ŝ −1 = Ŝ {ĉ † A , ĉB } Ŝ −1 = δ A,B that the matrices V S and U S are each others inverse; that is, the matrices U S and V S constructed from the coefficients in Eq. (C1) satisfy U S V S = V S U S = 1, where matrix products are taken as the usual matrix multiplication.If S is (anti-)unitary, then A,B = H A,B (H * A,B ) for S being (anti-)linear.Any anti-linear and anti-unitary symmetry Ŝ of Ĥ might be called a time-reversal symmetry because of its affect on the time-evolution operator Û (t) = e −i Ĥ t/ℏ ; if Ŝ is an anti-unitary symmetry, then Ŝ Û (t) Ŝ −1 = Û (−t).See Footnote 13 of Ryu et al. [125] for a discussion on uniqueness.Nevertheless, there is typically one physically motivated operator named the time-reversal transformation.When S is anti-unitary we have It has been shown [125] that in general U * S U S ∈ {−1, 1}, thus either Ŝ 2 ĉA ( Ŝ −1 ) 2 = ĉA and Ŝ 2 ĉ † A ( Ŝ −1 ) 2 = ĉ † A or Ŝ 2 ĉA ( Ŝ −1 ) 2 = −ĉ A and Ŝ 2 ĉ † A ( Ŝ −1 ) 2 = −ĉ † A are satisfied.This is familiar from the usual analysis on H; i.e. in the case U * S U S = −1, Kramer's theorem results.We now identify some symmetry operators Ŝ (and the corresponding U S ) for which Eq. (C3) is satisfied in the case of the non-magnetic (35) or the anti-ferromagnetic (B3) tight-binding Hamiltonian.

Inversion symmetry
An inversion operator I on the relevant electronic Hilbert space H is taken to be linear and unitary, and satisfy the physically motivated relation I 2 = id H . Any (anti-)linear operator can be specified by its action on a basis of the space in which it acts.For the non-magnetic Hamiltonian we define where the identification of U I follows from comparison with (C1) and using V S B,A = U S * A,B ; equivalently one could consider relations for the ĉ(α,σ),k .To ease notation, we do not include factors to make this unit-less, and take the unitarity to be understood as + + + + + + + + + + + + --------------------

FIG. 3 .
FIG.3.Magnetoelectric coefficient αme vs. film thickness and the distribution of the electronic charge density in the stacking direction (z).(a) αme vs. film thickness for various ∆D.Here the exchange splitting m of the Dirac cones associated with the top and bottom layers (resulting from the presence of surface magnetic dopants) is taken to be 0.1∆S and taken to vanish for all other layers.(b) Distribution of the electronic charge density along the stacking direction for the same values of ∆D used in (a), with the film thickness taken to be 20 layers and z = 0 taken to coincide with the middle of the center layer.The dotted (solid) lines identify the contribution of the n = 1 (n = 0) Landau levels to this distribution, which is (anti-)symmetric in z.(c) αme vs. film thickness for various exchange splittings m, with ∆D taken as 1.5∆S;(d) Distribution of the electronic charge density in the stacking direction for the same values of m used in (c) for a 20-layer thin film; (e) αme vs. exchange splitting m due to magnetic surface dopants for films of varying thickness.
for |∆ D | > |∆ S | (in the Z 2 -odd regime of the model) and α me → 0 for |∆ D | < |∆ S | (in the Z 2 -even regime).Moreover, numerical sums of the n > 0 (n = 0) energy eigenfunction distributions below the Fermi energy at half filling show that the charge density is (anti-)symmetric across the thin film, leading to a (non)vanishing contribution to the ground state dipole moment and thus to α me .This is shown for n = 1 (n = 0) by the dotted (solid) curves in Fig. 3 (b) and (d).