Axion coupling in the hybrid Wannier representation

Many magnetic point-group symmetries induce a topological classification on crystalline insulators, dividing them into those that have a nonzero quantized Chern-Simons magnetoelectric coupling (“axion-odd” or “topological”) and those that do not (“axion-even” or “trivial”). For time-reversal or inversion symmetries, the resulting topological state is usually denoted as a “strong topological insulator” or an “axion insulator,” respectively, but many other symmetries can also protect this “axion Z2” index. Topological states are often insightfully characterized by choosing one crystallographic direction of interest and inspecting the hybrid Wannier (or equivalently, the non-Abelian Wilson-loop) band structure, considered as a function of the two-dimensional Brillouin zone in the orthogonal directions. Here, we systematically classify the axion-quantizing symmetries and explore the implications of such symmetries on the Wannier band structure. Conversely, we clarify the conditions under which the axion Z2 index can be deduced from the Wannier band structure. In particular, we identify cases in which a counting of Dirac touchings between Wannier bands, or a calculation of the Chern number of certain Wannier bands, or both, allows for a direct determination of the axion Z2 index. We also discuss when such symmetries impose a “flow” on the Wannier bands, such that they are always glued to higher and lower bands by degeneracies somewhere in the projected Brillouin zone, and the related question of when the corresponding surfaces can remain gapped, thus exhibiting a half-quantized surface anomalous Hall conductivity. Our formal arguments are confirmed and illustrated in the context of tight-binding models for several paradigmatic axion-odd symmetries, including time reversal, inversion, simple mirror, and glide mirror symmetries.


I. INTRODUCTION
Recent progress in the classification of topological insulators (TIs) has been dramatic. Building on early work on the quantum anomalous Hall (QAH) effect [1], the first phase in these developments focused on systems with timereversal (TR) symmetry. In this context, two-dimensional (2D) systems are known as spin-Hall insulators, and threedimensional (3D) insulators are classified as weak or strong TIs [2]. While QAH insulators are characterized by an integer topological invariant, the TR-protected systems are described by Z 2 invariants {0, 1} (or {+, −}), corresponding to trivial and topological states, respectively. The additional presence of inversion symmetry was shown by Fu and Kane [3] to provide a simple means of determining the strong Z 2 index, based on counting the numbers of odd-parity states at highsymmetry points in the Brillouin zone (BZ).
The basic notion of the topological classification is the assignment of two insulators to the same class if and only if they can be connected by a continuous deformation of the Hamiltonian without gap closure, while preserving a given set of symmetries. We say that the symmetries "protect" the topological classification. At a formal level, the list of protecting symmetries can include chiral and particle-hole symmetries [4,5]. However, while relevant in the context of exotic phases such as superconductivity, these are rarely of interest for the description of ordinary insulators at the singleparticle level.
In 2011, Fu introduced the concept of topological crystalline insulators [6], where the set of symmetries used for the classification includes point symmetries, such as rotation or mirror operations, in addition to TR. The natural development of this trend is to consider the space group-or, for TR-broken systems, the magnetic space group-in setting up the topological classification. Here there has been dramatic progress in recent years. On the one hand, we have "theoretically complete classifications" based on K-theory that consider the effect of crystalline symmetries [7][8][9][10][11][12][13][14][15]. On the other hand, approaches based on symmetry indicators of elementary band representations [16][17][18][19][20][21] have emerged that are particularly well suited for high-throughput searches of topological materials [22][23][24][25].
Another very useful tool for characterizing topological states is the hybrid Wannier (HW) representation; see Ref. [26] for a recent overview. Here one selects one crystallographic direction, sayẑ, along which to transform from k space to real space via the standard construction of onedimensional (1D) maximally localized Wannier functions and their centers. This is done independently at each k y in 2D, or at each (k x , k y ) in 3D. Thus, the HW functions are localized in real space along z but remain as extended Bloch-like functions in the orthogonal direction(s). The Wannier-center positions can be obtained from a parallel transport analysis performed independently for each k-point string encircling the BZ in the k z direction. This kind of analysis also goes under the name of the "Wilson loop" [27][28][29], the generalized non-Abelian Berry phases corresponding to the HW centers are frequently referred to as "Wilson-loop eigenvalues." For a d-dimensional insulator, the flow of these HW centers as a function of wave vector in the orthogonal (d − 1)dimensional BZ, which we shall refer to as the Wannier band structure [30] often proves to be a very useful tool for determining its topology [31][32][33][34][35][36]. For example, the integer Chern number of a 2D QAH system, the Z 2 index of a 2D TRinvariant insulator, and the strong and weak indices of a 3D TR-invariant insulator, are easily diagnosed via an inspection of the flow of the HW centers. They can also provide strong hints as to the location in the BZ of the band inversion responsible for the topological state, and to the flow of surface energy bands for surfaces orthogonal to the wannierization direction.
One interesting aspect of the 3D TR-invariant class of insulators is that the Z 2 index that distinguishes between strong-TI and trivial (or weak-TI) systems also serves as an "axion Z 2 index" [37][38][39]. The axion coupling is expressed in terms of a phase angle θ , the "axion angle," that constrains the possible values of the anomalous Hall conductivity (AHC) on any insulating surface, as measured relative to an outward unit normal, to take values (e 2 /h)(N − θ/2π ) for integer N. Since θ is mapped into −θ under TR, and is only well defined modulo 2π , the presence of TR allows only two possible values of θ , 0 and π , which correspond to the trivial and topological Z 2 phases, respectively. This means that any insulating surface of a strong TI must exhibit a half-integer surface AHC, as measured in units of e 2 /h. If all facets of a crystallite are insulating and have the same value of surface AHC (i.e., the same integer N), then the crystallite as a whole behaves like a magnetoelectric material with an isotropic linear magnetoelectric coupling of θ = (e 2 /h)(θ/2π − N ). In this sense, θ is understood as describing a formal isotropic magnetoelectric coupling.
In fact, this formal coupling is π for ordinary strong TIs such Bi 2 Se 3 and its relatives, even though they normally show no magnetoelectric effect. This remains consistent with the previous considerations because the surfaces are necessarily metallic unless TR is broken at the surface. The topologically protected Dirac cones on the surface provide a canceling contribution, and the surface AHC vanishes (as it must if TR is present).
However, it was shown that inversion symmetry also protects the axion Z 2 index [38,40,41]. When TR is absent and the system is axion-odd, such systems are generally known as "axion insulators" [42]. A simple criterion was given by Turner et al. for determining the axion Z 2 index in the presence of inversion symmetry [41], generalizing the parity-counting analysis of Fu and Kane [3] to the TR-broken case. Because the symmetry protecting the bulk topology (i.e., constraining the values of θ ) is inversion instead of TR, and because inversion is never a good symmetry at any surface, it is much easier for the surface of an axion insulator to remain insulating. In some cases other symmetries may be present in addition and may force some facets to be metallic, and there is always the possibility that nontopological surface states will be present at the Fermi energy. Nevertheless, if the goal is to find insulators that naturally display a half-integer surface AHC response, then it appears that axion insulators are much more promising. Unfortunately, physically realized examples of 3D bulk crystals that behave as axion insulators have been very difficult to find. To date the closest to a physical realization seems to the antiferromagnetic topological insulator [43] MnBi 2 Te 4 , which has been the subject of much recent interest [44][45][46][47][48][49].
Other symmetry operations also reverse the sign of θ and thus support an axion Z 2 classification. These include simple mirrors and glide mirrors, rotoinversions, and time-reversed rotations and screws. In general, an insulating material whose magnetic space group contains such axion-odd operators has an axion Z 2 index, and if that index is nontrivial, then the material is guaranteed to have a half-integer QAH response on any insulating surface.
In this paper, we investigate the symmetry constraints imposed by the presence of axion-odd symmetry operations on the Wannier band structure, and more specifically, the additional constraints associated with the axion-odd topological state. Conversely, we show how the axion Z 2 index can often be deduced from an inspection of the Wannier band structure. In many cases, this involves only a visual inspection of the Wannier band structure, possibly including a counting of Dirac nodes between certain bands. In other cases it may require a calculation of the Chern number of a subset of Wannier bands, or a more complicated computation of Berry fluxes on truncated Wannier bands and Berry phases on the truncation boundaries.
The paper is organized as follows. In Sec. II we review the formal definition of the axion coupling as an integral of the Chern-Simons (CS) three-form over the 3D Brillouin zone (3DBZ). We also introduce the HW representation and the definitions of Berry potential and curvature on the Wannier bands. Then in Sec. III we review the expression for the axion coupling in the HW representation, which was originally derived only in the case of isolated Wannier bands. We generalize these expressions to the case that the bands touch to form isolated groups, and to the case that the entire Wannier band structure is connected, without yet focusing on symmetries that can quantize the axion angle. This is done in Sec. IV, where we classify the axion-odd symmetry operations according to whether or not they reverse theẑ axis, and whether or not they involve nonsymmorphic fractional translations alongẑ. For each of these classes, we discuss the simplifications that occur because of cancellations in the expression for the axion θ , and clarify when and how the axion index can be determined from an inspection of the Wannier band structure. Section V provides a discussion of which symmetries allow insulating versus metallic surfaces in axionodd insulators, or gapped versus gapless Wannier bands, and a consideration of the connections with higher-order topology. In Sec. VI we discuss several common cases, such as the presence of TR or inversion by itself, or the existence of a mirror or glide mirror symmetry. We illustrate each of these cases with explicit calculations for a tight-binding model embodying the symmetry in question. We summarize and conclude in Sec. VII. Finally, we also provide an Appendix in which the case of inversion symmetry is treated by an alternative approach involving counting of parity eigenvalues of the Wannier bands at high-symmetry k points, leading to conclusions in agreement with the analysis in Sec. IV.

II. PRELIMINARIES
In this section we review the expression for the axion coupling in the usual Bloch representation, and then we introduce the hybrid Wannier representation. This will set the stage and define the notation for the main developments that follow.

A. Axion coupling in the Bloch representation
We start from the usual expression for the axion coupling derived from the second Chern number, a topological invariant defined in 4D, through a dimensional reduction procedure. One obtains the expression for the axion coupling of a 3D insulator in the Bloch representation [37,38]. Here k = (k x , k y , k z ) runs over the 3DBZ and A α is a shorthand notation for the non-Abelian Berry connection A α,nm (k) = u nk |i∂ k α |u mk obtained from the cellperiodic Bloch functions |u nk . The trace in Eq. (1) is over the occupied valence bands indexed by m or n. The integrand of Eq. (1) is known as the CS three-form. Being a (pseudo)scalar quantity characterizing the ground state of a 3D insulator, one might expect the axion coupling θ to be gauge invariant, but this is not the case. In fact, a gauge transformation, that is, a k-dependent unitary transformation U mn (k) that mixes occupied bands, can cause θ → θ + 2π N for some integer N. This implies that the only well-defined part of the bulk axion coupling lives in an interval of 2π .
Physically, the axion coupling describes an isotropic contribution to the formal magnetoelectric (ME) tensor [2,50] In other words, an orbital magnetic field will induce a parallel polarization, or equivalently, an electric field will induce a parallel orbital magnetization, with a constant of proportionality given by We described this as a "formal" ME response because it manifests itself through the appearance of a surface AHC on any insulating surface facet that is given by In this exact relation, the ambiguity modulo 2π in θ is consistent with a freedom to prepare insulating surfaces with values of σ surf AHC differing by the quantum e 2 /h, e.g., by changing the Chern number of some surface bands, or of adding or deleting a surface layer with a nonzero Chern number [37,[50][51][52][53]. If all surfaces adopt the same branch choice-i.e., the same value of σ surf AHC -then the sample as a whole exhibits a true magnetoelectric response of −α CS , where the quantized part of the response has been absorbed into the branch choice for α CS . This phenomenon is a higher-dimensional analog of the modern theory of electric polarization [54], where the 2π ambiguity of the Berry phase reflects the inability to define the bulk polarization, or to predict the bound charge density of an insulating surface, except modulo a quantum.

B. The hybrid Wannier representation
For a given crystalline 3D insulator, we choose a reciprocal-lattice direction along which the HW representation will be constructed, and we orient the Cartesian axes such that the corresponding real-space direction is z. Given a gauge for the Bloch states |ψ nkk z , the corresponding HW states are expressed as Here k = (k x , k y ) is the wave vector in the perpendicular plane, i.e., in the projected two-dimensional Brillouin zone (2DBZ); c = 2π/b is the lattice constant along z, with b the magnitude of the shortest reciprocal lattice vector along z; l is an index that runs over unit cells along z; and n runs over the J occupied bands in the insulator, which is the same as the number of HW functions at each k in one vertical unit cell. Henceforth, the multiband gauge of the Bloch wave functions is always taken such that the |h lnk are maximally localized along z [55]. The HW wave functions h lnk (r) = r|h lnk are thus Bloch-like and cell-periodic [56] in the in-plane directions, while at each k they are the maximally localized Wannier functions of the effective 1D Hamiltonian H k at that k. Their centers z ln (k) = h lnk |z|h lnk (5) form the Wannier "bands" (or "sheets"). These are periodic in real space along z, z ln (k) = z 0n (k) + lc, as well as periodic in the in-plane reciprocal space, z ln (k) = z ln (k + G), where G is an in-plane reciprocal lattice vector. We shall frequently drop the explicit k dependence in the following. As noted earlier, this is essentially the same construction as that of the non-Abelian Wilson loop. The Berry connection is defined in the HW representation as where ∂ α = ∂/∂ k α with α = {x, y}. As a reminder, l and l run over unit cells in the z direction, and n and n label Wannier bands within the unit cell. Since there are J occupied bands in our insulator, n and n run from 1 to J. We shall also need to make use of the Wannier-band-diagonal Berry curvature Periodicity implies that A α ln,l n = A α 0n,(l −l )n (8) and therefore ln,ln = 0n,0n . Next, we discuss the gauge dependence of the quantities introduced above, as first presented in Ref. [57]. The Wannier bands are predetermined by the maximal localization procedure, so except in the case of degeneracy between Wannier bands [58], the most general gauge transformation that we need to consider is a Wannier-band-dependent phase twist This leads to new Berry connections The Berry curvature, however, is gauge invariant, which follows from the first line of Eq. (7) using We shall also have reason to define the (k-dependent) quantity which is fully gauge invariant, i.e., ln,l n = ln,l n .
To see this, note that the (z l n − z ln ) prefactor in Eq. (11) insures that only off-diagonal elements of A α contribute, and that the product of A x ln,l n and A y l n ,ln leads to a cancellation of the phase factors in Eq. (10).
Finally, we discuss the connectivity (sometimes called the "flow") of the Wannier bands. We restrict ourselves to the case that all bulk Chern indices are zero. We refer to the vanishing of the Chern number in the x-y plane as the "in-plane Chern constraint." The vanishing of the other two Chern indices guarantees that each Wannier band returns to itself (not to higher or lower partners) as one traverses the 2DBZ by a reciprocal lattice vector. Of course, it also returns to itself on any closed loop that does not wind by a reciprocal lattice vector, since all such loops are contractible [59]. Importantly, these considerations allow us to label the Wannier bands globally by integers that we take as increasing alongẑ.
A gap is said to exist between a pair of adjacent Wannier bands if these are not connected by degeneracies anywhere in the 2DBZ. A Wannier band is said to be isolated if a gap exists above and below it. A connected group of Wannier bands is a set of adjacent bands that are connected by degeneracies, but that are separated by gaps above and below the group. The HW sheet structure as a whole is said to be connected if there are no gaps; otherwise it is disconnected and is composed of M internally connected groups separated by M gaps per unit cell along z. These features are illustrated in Fig. 1.
While all the quantities in Eqs. (15) and (17) are gauge invariant, so that there is no ambiguity on this account, there is an important ambiguity of a different kind in θ z . It comes from the freedom to choose the set of Wannier bands assigned to the home unit cell; this affects θ z because of the appearance of z n in Eq. (17). By contrast, only differences of z values appear in Eq. (15), so this term is fully unambiguous.
To understand the ambiguity in θ z , recall that the labels n = {1, 2, . . . , J} simply count the bands in ascending order within the home unit cell. (As a reminder, we are considering here the case that all Wannier bands are isolated.) But how do we choose the unit cell? Which band shall we label as n = 1? One way to think about this is that the choice of home unit cell corresponds to the choice of one of the J gaps as the "primary" one, such that the counting starts with the Wannier band just above this gap. A different choice of primary gap just has the effect of shifting some subset of the Wannier bands along z by distance c. For each shifted band n, θ z in Eq. (17) gets shifted by where C n is the Chern number computed over this Wannier band. Since C n is necessarily an integer, this means that θ z , and also the total θ as given by Eq. (13), is only well-defined modulo 2π . Actually, such an ambiguity is expected, since we know on other grounds that θ is only well-defined modulo 2π , so this is not problematic. However, it is still something that we have to anticipate and deal with in the analysis that follows.
To clarify why θ xy is unaffected by the choice of unit cell, it is instructive to rewrite Eq. (15) as where we have introduced a condensed index notation μ = (ln), the sums over μ and μ both run over all Wannier bands in a large system of N unit cells, and μμ = i(z μ − z μ )A x μ A y μ . In other words, θ xy is just the sum of all elements taken per unit cell alongẑ. With this perspective, it is obvious that the ln labeling of the Wannier bands is irrelevant for θ xy . By contrast, it is impossible to write an expression similar to Eq. (19) for θ z , since the appearance of z μ itself, rather than the difference z μ − z μ , would render the average over N cells ill-defined.

Composite groups of bands
We now consider the case that at least some of the bands form internally connected composite groups, but the Wannier band structure as a whole remains disconnected. An example of a system of this type is shown in Fig. 1(b), where there are M = 2 connected groups, each consisting of a pair of bands joined by a nodal point. We do not expect such nodal points to appear generically, since accidental degeneracies have codimension three, and thus require fine tuning. However, such nodal degeneracies may sometimes be induced by symmetry, often occurring at the time-reversal invariant momenta (TRIM), or at other high-symmetry points or lines, in the 2DBZ. When such degeneracies are present, the evaluation of Eq. (13) becomes problematic because nn can diverge in the vicinity of the degeneracies between Wannier bands. In the case of a Dirac node, defined as an isolated nodal touching with linear dispersion of z n (k) close to the node, nn has a δ-function singularity at the node. To solve this problem, we can go over to a formulation in terms of a gauge-covariant treatment of the Berry curvature within each group.
To see this, let the ath group (a = 1, . . . , M) be composed of M a Wannier bands, and combine terms such that the contribution of group a is taken to be As a reminder, z n , nn , and nn refer to contributions coming from the home cell l = 0 (and l = 0) only. Then the total Chern-Simons coupling is where is the sum of Eq. (20) over groups a, and is identical to θ xy in Eq. (15) except that the prime on the sum indicates the omission of all terms with Wannier bands l n belonging to the same group as 0n. The sum over groups in Eq. (22) counts isolated bands by treating them as groups with M a = 1, and Eq. (13) with Eqs. (15) and (17) are recovered if all bands are isolated. But now the quantity inside the parentheses in Eq. (20) can be simplified by writing it as where is the diagonal element of a gauge-covariant Berry curvature matrix. Substituting Eq. (25) into Eq. (24) easily demonstrates the equivalence of these expressions. The Chern number of a connected group of bands can then be expressed in terms of the gauge-covariant Berry curvature as The advantage of this formulation is that nn remains a smooth and divergence-free function of k in the vicinity of degeneracies between bands in the group. This is a wellknown feature of the gauge-covariant Berry curvature and can be seen by expressing it as where∂ α |h n = Q a ∂ α |h n is the gauge-covariant derivative of the HW state and Q a = 1 − n∈a |h n h n | is the projector onto all states other than the Wannier bands in group a in the home unit cell. (The philosophy is similar to that used in defining the gauge-covariant derivative∂ α |u nk = Q nk ∂ α |u nk of the Bloch energy eigenstate |u nk , where Q nk = 1 − m =n |u mk u mk |, except that here we work with the spectrum of z, not that of H.) The projector Q a eliminates from∂ α |h n the divergences that would be present in ∂ α |h n arising from mixing between Wannier bands inside the same group.
In summary, in this part we have argued that Eqs. (21)-(24) form a robust set of equations that can be used to evaluate the CS coupling even in the presence of composite groups of internally-degenerate Wannier bands.
In case there is a doubt about the correctness of this expression, we can consider the application of a small symmetrylowering perturbation λV that gaps out the degeneracies between bands within the group. In this case we know Eq. (13) is correct, and we can argue that it is equivalent to Eq. (21), and then take the limit λ → 0. Since Eq. (21) is insensitive to degeneracies within the group, we can conclude that its evaluation at λ = 0 provides the correct CS coupling.

B. Introduction of a cutting surface
When it comes time to consider the case of a fully connected Wannier band structure, it will not be possible to assign Wannier bands to the home unit cell without cutting through the Wannier bands in some way. To prepare for this, we begin by returning to the case of isolated Wannier bands as in Sec. III A 1, and we define a "cutting surface" z cut (k) that is smooth and periodic in k, such that all Wannier bands with z cut (k) < z ln (k) < z cut (k) + c are assigned to the home cell l = 0. As a reference, if the cutting surface lies entirely inside the primary gap, as illustrated in Fig. 2(a), then the definition of the home unit cell is the same as it was in Sec. III A.
Instead, let z cut (k) be increased such that it cuts through one or more of the Wannier bands, as illustrated for a single Wannier band in Fig. 2(b). This has the effect that z n is shifted upwards by c for all bands lying below the new cut. Following an earlier argument, the contribution of Wannier bands lying entirely below the cut is changed by 2π times a Chern integer, making no change to θ modulo 2π . As for band n pierced by the cut, we define S n to be the region of the 2D plane for which z n lies below the cut, and C n is its boundary, i.e., the intersection with the cut, as shown in Fig. 2(c). The reassignment of the HW states inside S n shifts z n upwards by c, so that the term θ z in Eq. (17) changes by an amount where is the Berry phase evaluated on C n . In general, the boundary C n could be multiply connected, in which case the sum over loop Berry phases is implied in Eq. (28). The total change in θ z is then where φ cut is the total Berry phase from all Wannier bands that intersect z cut (k). Of course, this quantity is only well defined modulo 2π , but this is not an issue since the axion coupling has the same indeterminacy. However, the expression for θ xy in Eq. (15) is unchanged, since no matter the labeling, all pairs of Wannier bands at the same k eventually enter the sum in Eq. (15) in the same way as before. The overall change in Eq. (13) is then just θ = −φ cut . To correct for this, we just have to add back a piece to cancel Eq. (30), and we arrive at As a reminder, θ z is still evaluated as in Eq. (17), but now with the band label n running from 1 to J independently at each k beginning with the first band above the cutting surface, and the Berry-phase term φ cut accounts for the contributions of the loops of intersection of the cutting surface with the bands. Equation (31) is one of the principal results of the present work.

C. Connected Wannier band structure
The result in Eq. (31) was derived for the case that all bands are isolated, but we now wish to consider a fully connected Wannier band structure, in which no bands are isolated and no gaps occur. FIG. 3. Sketch of degeneracy region (shaded) associated with two Wannier bands z n (k) near a point node, and its projection onto the disk-shaped region R α in the 2DBZ (bottom).

Degeneracy regions
As a first step, we revise Eq. (31) for the case of composite groups considered in Sec. III A 2. In particular, we consider the case that the cutting surface z cut (k) cuts through one or more of the connected bands within a group (but avoiding degeneracies). Unfortunately, we cannot simply cut through a connected group that is being treated using the gaugecovariant formulation of Sec. III A 2, because φ cut is defined as the Berry phase on an individual Wannier band; this is equal to the integral of the nn over the enclosed area, but not that of nn .
To remedy this problem, we can adopt a more restrictive treatment of degeneracies, as follows. Wherever there is a degeneracy between Wannier bands, we identify a degeneracy region (DR) surrounding the degeneracy. The αth DR consists of a set of M α adjacent bands (M α 2) that are involved in the degeneracy, and a small region R α surrounding the degeneracy in the 2DBZ, as illustrated in Fig. 3. If a cutting surface z cut (k) is present, then we shall insist that it be chosen to avoid all the DRs, so that each DR lies entirely inside the home unit cell. We then collect together the terms in Eqs. (13)-(15) that only involve bands inside the DR to get a contribution Here we have followed the same approach leading to Eq. (20), but now restricting the k integral only to the region R α , and the sum to run only over the M α bands involved in the degeneracy. We then write the total Chern-Simons axion coupling in Eq. (31) as where and The double prime on the first sum in Eq. (34) indicates that all terms coming from within a DR are to be omitted, and in Eq. (35) it excludes terms where both bands lie in the same DR in the same cell; these omissions are compensated by the second term in Eq. (34). Equation (33) is the desired formula, which remains robust in the presence of degeneracies, including in the connected case in which no gaps are present. It provides the formal solution to the problem of expressing the axion coupling θ in the HW representation, even in the case of fully connected Wannier bands, and is one of our principal results.
While Eq. (33) might be somewhat awkward to implement in practice, involving as it does the choice of some DRs that need to be treated differently while integrating over the 2DBZ and summing over bands, we shall mainly be interested below in cases in which some symmetry is present that quantizes θ to 0 or π . In such cases, we shall insist on choosing the cutting surface and the DRs in such a way as to respect those symmetries, so that the same kinds of symmetry arguments used for the HW sheet contributions can also be used for the DR contributions. Thus, in cases where the θ z contribution would vanish in the absence of DRs, θ z vanishes in their presence as well. In other cases, we shall argue shortly that one can take a limit in which the size of the DRs goes to zero. Either way, an explicit calculation of the contribution of a DR can typically be avoided. In fact, it often happens that only the last term φ cut survives in Eq. (33), so that the axion coupling is given just by the Berry phase on the cutting loop (or the total Berry phase in the case of multiple loops).
Of course, whenever the Wannier band structure is actually disconnected, it is simpler to return to Eq. (21), where no cutting surface is needed and the bands can be indexed by counting from above some chosen gap. Again, symmetry will often allow us to decide the value of θ based on rather general features of the Wannier band structure in this case as well.

Shrinking the degeneracy regions
Up to this point, we have avoided specifying the nature of the degeneracies between Wannier bands, which in general could occur at point nodes, or along lines, or even over a plane spanning the 2DBZ, depending on the type of symmetries present. Henceforth, we will focus on point touchings of two or more Wannier bands, commenting only occasionally on the case of higher-dimensional degeneracies. Then, from a formal point of view, we can shrink the size of the αth DR surrounding a point node to a disk of some small radius in the 2DBZ. This may be problematic computationally, since the immediate vicinity of the DR may become difficult to treat without the gauge-covariant formulation, but as a formal manipulation it is permissible. Furthermore, we can argue that the contribution θ DR α vanishes in the limit that → 0, since nn remains finite as the degeneracy is approached (since it only "sees" Wannier bands outside the degenerate group), and the area of the disk goes to zero.
Formally speaking, then, we can simply neglect the contributions from the DRs in the small-DR limit. This will prove useful in analyzing the contributions to θ in the presence of certain symmetries, as we shall see.

IV. SYMMETRY CONSIDERATIONS
For the remainder of the manuscript, we restrict ourselves to insulating systems with axion-odd symmetries, which are classified by the axion Z 2 index. In this section we apply the formalism developed above to such cases. We consider three different classes of axion-odd symmetries, explaining when and how the axion Z 2 index can be inferred from an inspection of the Wannier band structure. When it cannot, we indicate what is the additional information that is needed to determine the axion Z 2 index in the HW representation.

A. Axion-odd insulators
For the remainder of the manuscript we restrict ourselves to insulating systems with symmetries that protect the quantization of θ to 0 or π . This symmetry could be time reversal (TR), in which case a spinful system with θ = π is usually denoted as a strong topological insulator (TI). It also could be inversion I, in which case the system is usually called an axion insulator. However, many other symmetries can quantize θ , such that θ/π = 0 or 1 defines an "axion Z 2 index." We use the term "axion-odd insulators" to refer to systems in which the nontrivial Z 2 index is protected by one of these symmetries, with TR-protected strong TIs and inversion-protected axion insulators as special cases.
In this section we survey the symmetries that can quantize the axion coupling, and describe their consequences for the Wannier band structure. Moreover, we show that in many cases it is possible to determine the axion Z 2 index without recourse to the expressions given in Sec. II. In particular, it is often enough just to have a knowledge of some elementary features of the Wannier bands, such as the type and number of touchings between bands, or the total Chern numbers of certain bands or band groups.
To decide whether θ is quantized to 0 or π by a set of crystal symmetries-i.e., whether the axion Z 2 index is protected-we can just look at whether there are any elements in the magnetic point group that reverse the sign of θ . We shall call these the "axion-odd" symmetries, and they are composed of the proper rotations composed with TR and the improper rotations not composed with TR. If one or more of these symmetries is present in the magnetic point group, then θ is quantized to be 0 or π , i.e., the Z 2 index exists.
In such a case, we can argue as follows that the θ xy term in Eq. (15) must vanish. Like the θ z term, θ xy has its sign reversed by any axion-odd symmetry operation. However, unlike θ z , θ xy has no quantum of ambiguity, as we saw at the end of Sec. III A 1. As a result, θ xy is immediately forced to vanish in the presence of such a symmetry, while θ z is not. The restricted sums θ xy and θ xy in Eqs. (21) and (33), respectively, will vanish as well, since the assignment of bands to composite groups automatically respects the symmetry, and we chose the shapes and locations of the degeneracy regions to respect the symmetry as well. For this reason, we ignore the TABLE I. Axion-odd point-group symmetries and their decomposition as g = g g ⊥ .
θ xy terms in the considerations that follow, and the remaining question is whether the θ z term, taken together with φ cut if a cutting surface is present, yields 0 or π . In many cases the magnetic point group may contain several axion-odd symmetry operations, but any one of them is enough to quantize θ . Therefore, we shall just consider each kind of axion-odd point symmetry in turn, and study its consequences, while keeping in mind that other symmetries may exist as well.
Let g be the axion-odd magnetic point symmetry in question. For this g, we choose our Cartesian frame such that the directionẑ is either invariant or reversed by g, and then we carry out the wannierization along theẑ direction. In some cases this choice can be made in more than one way; for example, for a mirror,ẑ could be chosen to lie in, or normal to, the mirror plane. In such cases more than one avenue of investigation may be available, with one possibly being more advantageous, but for now we just assume that some such choice has been made.
In this chosen frame, the possibilities for the spatial part of the point-group operator can be enumerated as follows. Ifẑ is not reversed, then it is the identity E ; a proper rotation C n by 2π/n about theẑ axis; or a reflection M d across a plane containing theẑ axis. Ifẑ is reversed, then it is the inversion I; a mirror M z across the x-y plane; an improper rotation S n = M z C n for n = {3, 4, 6}; or a twofold rotationC 2 about an axis lying in the x-y plane. We attach a prime to denote composition with TR, so that E represents TR itself. This establishes our notations for the symmetry operations. We also decompose g = g g ⊥ , where g is E or M z , and g ⊥ is the in-plane part of the operator, including the TR component if present. For orientation, we list in Table I the axion-odd pointgroup symmetry operations g along with their decomposition into g and g ⊥ .
Keeping in mind that we are restricting ourselves to the case that g is axion-odd, we can distinguish three general cases: (1) Operation g reversesẑ. We assume a choice of z-origin such that the corresponding space-group operation has no fractional translation alongẑ [61].
(2) Operation g preservesẑ, and the corresponding spacegroup operation has no fractional translation alongẑ.
(3) Operation g preservesẑ, and the corresponding spacegroup operation involves a fractional translation cẑ/m. The first case corresponds to g = M z , while the second and third pertain to g = E . In all cases there may also be a fractional (nonsymmorphic) in-plane translation τ ⊥ in the space-group operation associated with g; these play a role in determining when and where point nodes of degeneracy may occur between neighboring Wannier bands but will not concern us here. If point nodes do occur, then these will typically be located at z = 0 or z = c/2 for g = M z , and at general locations for g = E .
We denote the axion-odd space-group operations as being either z-reversing, z-preserving, or z-nonsymmorphic according to the cases enumerated above. In the following subsections we shall consider each of these cases in turn.
Before proceeding, we list the transformation rules for the Wannier bands and their Berry curvatures. Under an operation {g|cẑ/m + τ ⊥ } of the space group, a Wannier band transforms as [31] where the minus sign applies when g ⊥ includes TR. The Berry curvature is a pseudovector pointing alongẑ, and as a result it transforms as and the same rules apply to the gauge-covariant Berry curvature nn .

Symmetry operation reversesẑ
For the operations that reverseẑ, listed in the first four rows of Table I, it follows from Eq. (36) that each HW center at (z, k) in band n has a partner at (−z, ±g ⊥ k) in some band n (again, the minus sign applies when g ⊥ includes TR). Moreover, according to Eq. (37) the Berry curvatures of the Wannier bands are identical at these locations. At least in a simple situation such as that in Sec. III A 1, it therefore appears at first sight that all contributions to the θ z term given by Eq. (17) will cancel in pairs when summing over all Wannier bands and integrating over the 2DBZ, because of the sign reversal of z. However, it is not always possible to choose the home unit cell in such a way that the cancellation in Eq. (17) is complete.
a. Disconnected case. To see this, first consider the case of a disconnected Wannier band structure. As discussed earlier, we arrange the Wannier bands into a set of internally connected groups that are isolated from one another by gaps, and we work with Eqs. (22) and (24), which were formulated for this case. Clearly the arrangement of HW groups has to respect the z ↔ −z symmetry, so we can proceed as follows.
We organize the internally connected band groups into three collections as follows. First, we see whether there is a connected group centered at z = 0, and if so, we call it the "origin-centered" group. In general, this can be done by listing all connected groups that pass through z = 0. If the number of such groups is odd, then the central one is the origin-centered group, and if not, then there is no origincentered group. A similar construction at z = c/2 allows us to identify the "boundary-centered" group, if there is one. We FIG. 4. Sketches of hybrid Wannier band structures along a path in k-space that captures all the degeneracies of a system with an axion-odd symmetry that reversesẑ. (a) Disconnected Wannier band structure, decomposed into an origin-centered group (blue shading), a boundary-centered group (red shading), and an uncentered collection (yellow shading). The axion index depends on whether the Wannier bands in the red region have an odd Chern index, which in turn is determined by a counting of Dirac nodes; here θ = 0 because the number of nodes in that region is even. (b) Connected Wannier band structure, decomposed into a thin slice centered at z = c/2 (red shading) and the remainder (blue shading). Here θ = π since there is an odd number of Dirac cones in the red region. then see whether there are any remaining HW groups lying between the origin-centered and boundary-centered groups (or between the corresponding gaps if these groups are absent). If so, then we collect these, together with their z → −z symmetry partners lying below the origin-centered group, into an "uncentered" collection of band groups. This leads to a unique procedure for defining the home unit cell as composed of the union of the origin-centered group, the boundarycentered group, and the uncentered collection. Figure 4(a) shows an example in which there are origin-centered and boundary-centered groups composed of two Wannier bands each, and an uncentered collection accounting for two more Wannier bands, for a system with six occupied valence bands.
We can also obtain the contributions of these three collections to various physical properties. Let C OC , C BC , and C UC denote the total in-plane Chern numbers of the origincentered, boundary-centered, and uncentered collections, respectively, as computed from Eq. (26). We can also obtain the contribution of each collection to θ z in Eq. (22). We argued earlier that θ xy = 0 in the presence of any axion-odd symmetry, so we denote the θ z contributions as simply θ OC , θ BC , and θ UC , respectively, where we have also dropped the prime for brevity. Then the total axion coupling is just Now the contribution θ OC from the origin-centered group takes the form given in Eq. (24), and this clearly vanishes in view of the cancellations between contributions at z and −z inside this Wannier band group (recall that nn has the same sign in the canceling pair). Similar arguments imply that θ UC also vanishes, since the groups below and above z = 0 cancel in pairs. However, there is no such cancellation for the Wannier bands in the boundary-centered group, because their symmetry-mapped partners are centered about z =−c/2 and are outside the chosen home unit cell (see Fig. 4(a), where the home cell consists of the union of the four shaded regions). Instead, each Wannier band at (z, k) in the boundary-centered group has a partner in the same group at (c − z, ±g ⊥ k). Again 155130-9 nn is identical at both locations, so when inserting z n into Eq. (24), we make no mistake if we treat both as located at the average location z = c/2. But then Eq. (24) becomes where Eq. (26) has been used. Since this is the only contribution, we have that θ = π if and only if the Chern index of the boundary-centered group is an odd integer. Now recall that we assumed that all bulk Chern numbers vanish in our material, so that C BC + C OC + C UC = 0. Clearly C UC is even, since the uncentered collection is always composed of groups that contribute in pairs. Thus, C BC is odd if and only if C OC is odd. That is, the existence of an odd-Chern boundary-centered group also implies the existence of an odd-Chern origin-centered group. We are thus guaranteed to get the same result using the boundary-centered or origincentered group in Eq. (39). The same conclusion follows from the freedom to shift the origin by c/2 along z.
b. Counting of Dirac nodes. In some cases it may be possible to deduce the Chern number of a boundary-or origincentered group by inspecting the nodal touching of Wannier bands. For example, suppose there are only two Wannier bands in the boundary-centered group, and that the locus of degeneracies between these bands consists only of some number K of Dirac point nodes [62]. In view of Eq. (36), these are typically found at z = c/2 in the Wannier direction, and at locations obeying ±g ⊥ k = k in the 2DBZ. For example, they may occur at some of the four projected TRIM (PTRIM) for g= E , or at other high-symmetry points or even at generic positions in the 2DBZ for other symmetries.
Because we are considering symmetries that preserve the Berry curvature while interchanging z and c − z, and thus interchanging the two bands, we know that each Wannier band carries the same Berry flux 0 . Now, if the symmetry is weakly broken in such a way as to gap the Dirac nodes, then an additional Berry flux of ±π is transferred to each Wannier band at each of the nodes. However, this has to result in an integer Chern number, so we conclude that 0 must be 0 or π (mod 2π ) if the number K of Dirac nodes is even or odd, respectively. Recalling Eq. (39), this means that the axion Z 2 index is odd only if K is odd. The same analysis can be applied to the origin-centered group.
This argument easily generalizes to the case that the number M BC of Wannier bands in the group is even, still assuming that the only touchings are Dirac-like. That is, the Chern number is odd if and only if the total number K tot of Dirac nodes is odd. In this counting procedure, a node involving an L-fold degeneracy is counted as L/2 Dirac nodes.
If the number M BC of Wannier bands in the group is odd, then we have to treat the central band separately. We let j = (M BC + 1)/2 be the index of this band, and assume that it is regular enough to have a well defined Chern number C j . Each nodal touching involving this band will also involve L/2 pairs of neighbors from among the remaining (M BC − 1) bands. These can be treated as before, providing L/2 Dirac nodes at one touching location, and a total number K tot when summed over the 2DBZ. The conclusion is that the total Chern index of the group, and thus the axion Z 2 index, is odd if and only if C j + K tot is odd. This case is not quite as convenient, because it requires the evaluation of the Chern index of the central band, in addition to a simple counting of Dirac nodes.
In this analysis, we have assumed that the locus of contact between Wannier bands consists only of point Dirac nodes. If cases arise involving high-order (e.g., quadratic) point nodes, or nodal lines or regions, then the analysis given above would need to be reconsidered.
c. Connected case. Essentially identical results involving the counting of Dirac nodes, and possibly the calculation of the Chern index of a central band, can be derived for the case of a connected Wannier band structure. To do so, we first establish a "nominal cell boundary" z nom (k) located near z =−c/2 as follows. Let N be the number of Wannier bands passing through the point (−c/2,¯ ). If N is odd, then let z nom (k) be identified with the central one of these bands; if N is nonzero and even, let it be the average of the central two bands; and if N = 0, let it be the average of the next lowest and next highest bands about z =−c/2. Then z nom (k) is a surface centered at −c/2 that respects the symmetries of the system and passes through the point (−c/2,¯ ).
If N is even, then z nom (k) lies between two Wannier bands. We assumed a connected band structure, so even if N = 0 these bands must touch at one or more nodes, which lie on the surface z nom (k) by construction. If N is odd, then z nom (k) is identified with a Wannier band, which again touches with the next higher band by assumption. In either case, let the cutting surface be where δ > 0 is a small vertical shift. With this convention, the unit cell consists of all portions of the Wannier bands lying between z cut (k) and z cut (k) + c.
We first discuss the case of even N. An example is sketched if Fig. 4(b); there the symmetry is such that z nom (k) is perfectly flat at z =−c/2, and the unit cell consists of the union of the two shaded regions shown there. Regarding the degeneracies themselves, we let of Sec. III C 2 go to zero faster than δ → 0, so that the cutting surface avoids the DRs. Then θ z reduces to the first term of Eq. (34). This in turn can be broken into two contributions: one from the blue region z n (k) ∈ [z nom (k) + δ, z nom (k) + c − δ], and another from the red region Fig. 4(b). The first region is centered about z = 0, so all contributions cancel in pairs, while the latter region vanishes in the limit δ → 0, yielding a vanishing contribution to θ z [63]. Thus, θ z = 0. We argued earlier that the θ xy -type terms always vanish in the presence of axion-odd symmetries, so in fact the first two terms in Eq. (33) both vanish, and we are left with That is, the axion coupling is given by the sum of Berry phases of the vanishingly small loops of intersection of the Wannier bands with the cutting surface of Eq. (40) around the nodal points as δ → 0. If these are simple Dirac nodes lying on z nom (k), then we know that each such Berry phase is π , and it follows that the system is axion-odd if and only if the number of such Dirac nodes is odd, as in Fig. 4(b).

155130-10
If the number N of degenerate bands at (−c/2,¯ ) is odd, then arguments like those above yield where C j is the Chern number of the central band. Again, if the nodes are simple Dirac nodes with their nodal points lying on z j (k), then the system is axion-odd if and only if the sum of this central Chern index and the number of Dirac cones is odd. Note that we could equally well have chosen to work with a nominal cell boundary centered around z = 0 instead of z =−c/2; as in the disconnected case, this follows from the freedom to shift the origin by c/2 along z.

Symmorphic operation preservingẑ
Here we consider the case that g preserves the sense ofẑ, so g ⊥ = g has to be either TR itself, a time-reversed rotation about the z axis, or a simple reflection about a plane containing this axis, as summarized in Table I. For the moment we assume there is no associated fractional translation alongẑ; that case will be considered in the next subsection. Now the Berry curvature contribution from an area element at z on one Wannier band gets mapped onto one of opposite sign on the same sheet at the same z, giving canceling contributions to θ z .
In the disconnected case, this implies that θ = 0 and the system is always in the axion-even phase, since any gap can be chosen as the primary one defining the unit cell. Conversely, if the system is axion-odd, then the Wannier bands must be connected.
For the connected case, any cutting surface z cut (k) that respects the in-plane symmetries will automatically give θ z = 0, since the same kind of cancellations occur. In this case, the axion Z 2 index is just given by φ cut in Eq. (33). We can arbitrarily choose one Wannier band j, make the cut at z j (k) + δ for vanishingly small δ, and count the total Berry phase of the small loops of intersection of the Wannier bands with z cut (k). In the case that these are simple Dirac nodes, the system is axion-odd if and only if their number is odd.

Nonsymmorphic operation preservingẑ
Finally, we again consider theẑ-preserving case, but now assuming that the corresponding space-group operation involves a fractional translation c/m alongẑ for some integer m.  , k) always has a partner at (z + c/2, ±g ⊥ k), with the minus sign applying when g ⊥ involves TR. Moreover, all of the operations M d , E , C 2 , C 3 , C 4 , and C 6 reverse the sign of the Berry curvature , which is thus equal in magnitude but opposite in sign between these two partners.
If the band structure is disconnected, then there must be an even number of gaps. We choose a primary gap and divide the Wannier bands into two disconnected groups related to each other by Eq. (36): a "lower group" consisting of the first J/2 bands above the primary gap, and an "upper group" consisting of the rest. Consider a band n in the lower group; its contribution z n (k) nn (k) to the integral in Eq. (24) can be paired with a contribution at k = ±g ⊥ k in the upper group. The sum of these two contributions, when combined with the −1/c prefactor of Eq. (24), yields nn (k)/2. Since we know that the θ xy term does not contribute, Eqs. (21) and (22) then yield where LG is the total flux of Berry curvature in the lower group, with C LG being the corresponding Chern number (an integer). The system is thus axion-odd if and only if this Chern index is odd.
To prepare for the connected case, let us return to the disconnected case for a moment. We start with a cutting surface z (0) cut (k) that lies entirely in the primary gap, so that the lower group has the total Berry flux LG of Eq. (44), which we now relabel as (0) LG . We then raise the cutting surface to some chosen new location z cut (k) that cuts through some of the bands at the bottom of the lower group. (We also insist that z cut (k) avoids any DRs.) The new lower group now extends from z cut (k) to z cut (k) + c/2. As a result, the total Berry flux LG in this group changes for two reasons: because of the omission of terms below z cut (k), which is compensated by the addition of a term φ cut ; and by the addition of new contributions near z cut (k) + c/2, which is compensated by the addition of the same φ cut because of the sign reversal of nn . Thus, we find that (0) The quantity on the right side of this equation is 2π times an integer, even though the individual terms are not, and the system is axion-odd if and only if this integer is odd. Finally, we argue that this same formula applies to the connected case, since following in the spirit of the discussion in Sec. III B, we can imagine that we introduce a perturbation that opens a gap, place the cutting surface there, raise the cutting surface, and then close the gap. A concise representation of the final result is This formula should be correct for any choice of cutting surface, as long as it avoids the DRs.
In the case of Dirac cones connecting a pair of bands, we can choose a cutting surface just above the average of these bands, so that φ cut could again be obtained as π times the number of such Dirac cones. However, unlike in Secs. IV B 1 and IV B 2, here we cannot avoid doing an explicit calculation of the total Berry flux in half the unit cell. In this sense, Eq. (46)  omit left-handed screws, since they behave in the same way as right-handed ones. ) We start by assuming a disconnected band structure, and consider the example of {C 4 |c/4}. The number of gaps must be a multiple of four; we choose one such gap to define the unit cell in the usual way. Now {C 2 |c/2}, which is the square of {C 4 |c/4}, is also in the symmetry group. It simply translates the Wannier bands by a half lattice vector alongẑ with an extra C 2 rotation without affecting the value of the Berry curvature on a given patch of the Wannier band, which is just carried along to its new location. From the point of view of Eqs. (17) or (24) for the contribution to θ , the rotation component is irrelevant, since it does not affect the z n or nn value. Thus, for our purposes we can think of {C 2 |c/2} as defining a smaller "unit subcell" of height c = c/2, and we can compute θ by focusing on just one subcell containing J = J/2 Wannier bands.
For the case of isolated bands, for example, Eq. (17) becomes (the fact that we count half as many bands is compensated by the factor of two from the replacement of c by c in the prefactor), and a similar modification would apply to Eqs. (20) and (22) in the case of composite groups. Now the action of {C 4 |c/4} in the subcell is entirely analogous to the action of {C 2 |c/2} in the full cell; this is one of the cases that we studied above, and the same conclusions apply. That is, we conclude that θ is given by Eq. (44), where C LG is now interpreted as the total Chern number coming from the bottom half of the subcell, i.e., from the first J/4 bands. The case of a connected band structure under {C 4 |c/4} is handled using the same analogy to the {C 2 |c/2} case, and a formula like Eq. (46) applies.
The same strategy can be applied to reduce the remaining operations to previously studied ones in a similar way. All of these remaining operations, namely, {C 3 |c/3}, {C 6 |c/3}, {C 3 |c/6}, and {C 6 |c/6}, have the property that their fourth power is (modulo full translations) just a simple threefold screw, which divides the unit cell into three subcells. Thus, we can restrict our attention to a subcell of height c = c/3 containing J = J/3 bands.
For the first two operations, {C 3 |c/3} and {C 6 |c/3}, there are no remaining fractional translations within the subcell. Their third powers correspond to E and C 2 , respectively, so in these cases the analysis of the subcell is analogous to the symmorphic case discussed in Sec. IV B 2, and the conclusions found there apply. Specifically, θ is trivial in the disconnected case, and a counting of Dirac nodes can provide the value of θ in the connected case.
For the last two operations, {C 3 |c/6} and {C 6 |c/6}, the subcell of height c/3 is subdivided by a further half translation of c/6. The situation is similar to the {C 4 |c/4} case discussed above, except that now the subcell is of size c/3 instead of c/2. This subcell has lower and upper portions containing J/6 Wannier bands each, related to each other by an operation containing TR. Thus, the analysis of the subcell in these cases is analogous to the that of the full cell in the nonsymmorphic {E |c/2} and {C 2 |c/2} cases, respectively. In the disconnected case, Eq. (44) applies with C LG interpreted as the Chern number of the first J/6 bands, and in the connected case a formula like Eq. (46) once again applies [65].

Summary
In summary, for the case of a disconnected band structure, the rules for determining the axion index are quite simple. For a z-reversing symmetry operator, we test for the presence of a boundary-centered group with an odd Chern number, or equivalently, an origin-centered one; if present, the system is axion-odd. For symmetry operators that preserve the sign ofẑ, we can consider three cases. If there is no associated fractional translation alongẑ, then the system is forced to be axion-trivial. For the case of a half translation, the system is axion-odd if and only if the total Chern number in a half unit cell is odd. For smaller fractional translations, we can define a subcell that tiles the full cell under simple screw rotations, and apply the earlier analysis to the subcell.
The rules for the connected case are more complicated, but in the simplest case that the connection is via Dirac nodes, the results can be expressed in terms of a counting of Dirac nodes, and may also require the calculation of the total Berry flux of a band or subgroup of Wannier bands in the unit cell.
Note that we may have some choice regarding the crystallographic direction along which to perform the wannierization. In the case of I or E , any convenient primitive reciprocal direction will do. In the case of a mirror, we may choose either an axis normal to or lying in the mirror plane, applying the analysis of Secs. IV B 1 or IV B 2, respectively. (If θ = π , then the Wannier bands are necessarily connected for the choice of wannierization in the mirror plane.) The case of C 2 rotations provides a similar choice of options, i.e., the HW construction can be done either normal to the C 2 axis (Sec. IV B 1) or along it (Sec. IV B 2). Of course, there may be more than one axionodd symmetry in the magnetic point group, in which case one also has a freedom to choose which operation to select. Thus, it may often be possible to simplify the determination of the axion index by an appropriate choice of symmetry and setting.

V. BULK-BOUNDARY CORRESPONDENCE
Even though the axion coupling θ is a bulk quantity, it has important implications for surfaces. A given surface facet of an axion-odd insulator will either be metallic with an odd number of Dirac cones and a vanishing surface AHC, or it will be gapped with a half-quantized surface AHC of (N + 1/2)e 2 /h. Furthermore, if all surfaces are gapped, then θ = π implies a higher-order bulk-boundary correspondence. That is, both bulk and surfaces are gapped, but one-dimensional chiral channels are present.
Previous work has shown that there is often a close relationship between the Wannier band structure in a given crystallographic direction, and the surface band structure for a surface facet normal to that direction. For example, the flow of the Wannier bands for a weak or strong TI in a 3D TRinvariant insulator is reflected in the presence and character of the surface energy bands on the corresponding surface facet. Thus, it is of interest to explore the relationship between the Wannier band structure and the corresponding surface band structure in axion-odd insulators.
In Secs. V A and V B we address two closely related questions concerning spectral flow in an axion-odd insulator. First, we identify the symmetry conditions under which topologically protected metallic states must exist on a given surface. Conversely, when are insulating surfaces allowed to appear? We then discuss the conditions under which the flow of the bulk Wannier bands in a given direction is topologically protected. Conversely, when can an axion-odd insulator have a gapped Wannier band structure? In Sec. V C we discuss the identification of axion-odd insulators as second-order topological insulators, discuss the conditions on the occurrence and the number of chiral channels on a given hinge, comment on their robustness, and make connections with previous literature on second-order topological phases.

A. Protected flow of the surface energy bands
To answer the first question posed above, we ask whether the surface magnetic point group contains any axion-odd symmetry operations. If so, then the surface AHC would be forced to vanish (see below). However, this is inconsistent with the half-integer surface AHC that, according to Eq. (3), must occur for any insulating surface of our assumed axionodd bulk insulator. Thus, under these conditions the surface is necessarily metallic. In this case, Eq. (3) becomes [53] where the extra phase angle φ is the Berry phase summed over all the Fermi loops in the surface BZ. Since each Dirac cone contributes a Berry phase of π , a topological metallic surface must have an odd number of them to cancel the θ = π contribution from the bulk. If no axion-odd symmetries are present at the surface, then the surface is allowed to be insulating. Of course, whether it is really insulating or not will depend on details of the surface electronic structure for the particular surface termination, but the bulk topology no longer requires a metallic state.
One immediate consequence is that if TR itself is a symmetry of the bulk, and also of the surface of interest, then this surface must be metallic. This is the well-known statement that strong TIs have metallic surface states on all surfaces, unless TR is somehow broken on the surface. By contrast, if inversion is the only axion-odd symmetry protecting the axionic topology of the bulk, then all surfaces are allowed to be insulating, because no surface preserves inversion symmetry.
More complicated situations can be analyzed as follows. For a given facet on the surface of an axion-odd crystal, we first reduce from the 3D bulk magnetic space group to the 2D magnetic space group that describes the highest symmetry that the surface could possibly have. This will consist of all operations of the 3D group that leave the surface normal (chosen asẑ) invariant, and that involve no translation (either by a full or a fractional lattice vector) alongẑ; the axion-odd operations among these, if any, are of the z-preserving type. We then construct the surface magnetic point group in the usual way, by listing the point operations appearing in the space group. Finally, we omit any operations that are not axion-odd. The remaining operations that may be present are the z-preserving operations in Table I: TR itself, the C n timereversed rotations aboutẑ (n = {2, 3, 4, 6}), and reflections M d about a plane whose unit normal lies in the surface. All these symmetries flip the sign of the surface AHC (a pseudovector pointing alongẑ), forcing it to vanish. Hence, if any of these symmetries are present at the surface, then the surface must be metallic, with φ = π in Eq. (48).
Of course, any given surface may have lower symmetry than the one allowed by the above considerations. For example, a bulk-allowed M d mirror symmetry may be spontaneously broken by the formation of a symmetry-lowering surface reconstruction, or TR may be broken at the surface by the spontaneous appearance of magnetic order. In general, if any of the axion-odd symmetries survive at the surface, then it must be metallic; otherwise, it can be insulating.
We mentioned earlier that if inversion is the only axion-odd bulk symmetry, then all surfaces are allowed to be insulating. One may ask whether this is also true for any other symmetries. The answer is yes. Consider, for example, the case that an improper rotation S 3 , S 4 , or S 6 about the z axis is a symmetry; then for any facet normaln, the component of n on the (x, y) plane is rotated by the symmetry, and the z component is reversed, so there are no surfaces that obey this symmetry. The time-reversed screw rotations also have this property. This time they do preserven for surfaces normal to the rotation axis, but the fractional translation that comes with the screw operation is always inconsistent with such a surface.

B. Protected flow of the bulk Wannier bands
As we have seen, a surface of unit normalẑ on an axionodd insulator is required to be metallic if any bulk axion-odd symmetries of the z-preserving type are preserved at the surface. Assume that such a high-symmetry metallic surface has been prepared. Following Refs. [66,67], the surface spectrum can be continuously deformed into the bulk Wannier spectrum obtained by wannierizing alongẑ. Briefly, the actual surface band structure is connected to a model surface obtained by energetic flattening and spatial truncation, and then the abrupt truncation is replaced by a crossover region whose width is allowed to diverge. In this limit, these authors show that the Wannier band structure is recovered. Since the surface spectrum flows by hypothesis, the Wannier spectrum must flow as well.
This argument, based on the bulk-boundary correspondence, reproduces the conclusions of Sec. IV B 2, where purely bulk considerations led to the conclusion that the Wannier bands must be connected if a z-preserving axion-odd symmetry is present. These considerations also imply that the flow should occur via an odd number of Dirac nodes between adjacent Wannier bands.
We emphasize, however, that the bulk-boundary argument does not necessarily work in reverse. That is, if the surface band structure does not flow, this does not necessarily imply the same for the bulk Wannier bands. Sometimes a z-reversing or z-nonsymmorphic symmetry can lead to flow of the Wannier bands in a "fragile" sense. By this we mean that the flow can be destroyed, without breaking any axion-odd symmetries, by adding some weakly-coupled trivial bands to the valence manifold; this behavior was demonstrated in Ref. [68] for the case of inversion symmetry. This is again consistent with the analysis of Secs. IV B 1 and IV B 3, where both connected and disconnected Wannier band structures were found to be compatible with θ = π . This will be illustrated for the case of a pyrochlore model with inversion symmetry in Sec. VI B 2, and for the case of a model with glide mirror symmetry in Sec. VI D 2. Since z-reversing and z-nonsymmorphic operations are never symmetries of the surface, there is never any such fragile protection of the flow of the surface energy bands in these cases.

C. Higher-order topology
Here we revisit the bulk-boundary correspondence for an axion-odd insulator, now relaxing the requirement that the surface facet in question is itself invariant under the symmetry, as was assumed in Sec. V A. We consider a macroscopic crystallite and only require that the boundary as a whole is left invariant under an axion-odd symmetry.
Under these conditions, facets that are mapped to themselves under the symmetry are necessarily metallic, but those that are not can be insulating. When that happens, an insulating facet gets mapped into a sequence of n−1 other facets for a symmetry operator of order n; we take these to comprise an "orbit." All facets in the orbit have the same magnitude of surface AHC, but the sign alternates from facet to facet. There could be more than one orbit in general. Whenever insulating facets with different values of surface AHC meet, the difference must be an integer multiple of e 2 /h, with that integer corresponding to the number of chiral modes propagating along the connecting hinge [69]. Since the sign of the AHC cannot possibly be uniform on all facets, some hinge modes are necessarily present.
Such axion-odd crystals thus constitute examples of second-order TIs [70][71][72][73][74][75][76][77][78][79], i.e., crystals whose bulk topology in d = 3 dimensions gives rise to protected boundary states of dimension d −2 = 1. For example, consider an axion-odd insulator with C 4 symmetry [73], and take a crystallite with a tetragonal shape. The top and bottom surfaces are required to be metallic because they are mapped onto themselves by C 4 symmetry. Instead, the four side surfaces can be insulating, and we assume they are. In that case they form an orbit in which the AHC of adjacent side facets alternates between ±(N + 1/2)e 2 /h for some fixed integer N. This results in 2N + 1 chiral modes on each of the four vertical hinges, whose chiralities change sign from one hinge to the next.
A similar story applies to the case of a tetragonally shaped crystallite invariant under not C 4 but S 4 symmetry [78]. In this case no facet is left invariant under the axion-odd symmetry S 4 , so that all six facets are allowed to be insulating. The surface AHC alternates between ±(N 1 + 1/2)e 2 /h on the side-surface orbit, and between ±(N 2 + 1/2)e 2 /h on the topand-bottom orbit, where the integers N 1 and N 2 depend on the surface preparation. The pattern of hinge modes will depend on the integers N 1 and N 2 , but no choice of these integers can prevent the appearance of hinge modes.
An intrinsic [12,80,81] second-order TI is one whose bulk topology is such that metallic hinge modes are required if the entire crystallite, including the crystal termination, is invariant under the topology-protecting symmetries. It follows that the axion-odd insulators discussed here belong to this class. By contrast, the classification of extrinsic [12,80,81] secondorder TIs defines topological equivalence based on continuous deformations that preserve bulk and surface gaps [12]. In this context, we note that the chiral hinge modes on axionodd insulators are robust against small deformations of the Hamiltonian that preserve bulk and surface gaps, even if they break the axion-odd protecting symmetries [69]. We also note that in some circumstances, chiral edge channels can appear not only on hinges, but also on steps and domain-walls within a single facet [34,43]. However, sufficiently strong surface deformations can always remove the chiral channels, leading to the topological magnetoelectric effect. This can be arranged, for example, by appropriately "decorating" the surfaces of the crystallite with quantum anomalous Hall layers in such a way that the boundary has a uniform half-quantized AHC [34]. In any case, it is clear that the Wannier band structure cannot predict when and where extrinsic hinge states will appear, since it is a purely bulk construction.

VI. SPECIAL CASES AND NUMERICAL TESTS
In this section we consider four axion-odd symmetry operations that deserve special emphasis: TR (i.e., E ); inversion I; mirror (either M z or M d , depending on the wannierization direction); and glide mirror {M d |c/2}. In each case we shall introduce an elementary tight-binding model that can be used to illustrate some of the behaviors expected from theory.

General considerations
Here we consider the case that TR itself, E , is a symmetry of the system. This is a z-preserving symmetry as discussed in Sec. IV B 2, and the conclusions given there apply directly.
We first consider the case of spinor electrons, so that (E ) 2 = −1. As is well known, all energy bands of a 2D system appear as Kramers pairs at all four TRIM. In a similar way, the Wannier bands of a 3D system all appear as Kramers degeneracies at the four PTRIM, as discussed in Ref. [31]. Therefore, any given Wannier band j is either Kramers degenerate with band j − 1, or with band j + 1, at each of the four PTRIM. Let us call these "down-touchings" and "up-touchings," respectively; their numbers are N d ( j) and N u ( j) = 4 − N d ( j), respectively. The next band clearly has , there are three cases to consider: (1) N d alternates between 0 and 4 as bands are counted. The Wannier bands are thus "glued together" in pairs, and are generically gapped from the next higher or next lower pair [82]. From this it follows that the axion Z 2 index is trivial. As discussed in Ref. [31], it also follows that the Z 2 indices are trivial on the four TR-invariant planes k x = 0, k x = π , k y = 0, and k y = π in the 3DBZ. This system is either topologically trivial, or it is a weak topological insulator with indices (0;001), i.e., equivalent to a stack of quantum spin Hall (QSH) layers along z. (2) N d = 2 for all bands. The Wannier band structure is connected; there are no gaps. These systems correspond to weak topological insulator configurations, this time corresponding to a stacking of QSH layers along x or y or in the diagonal xy direction. According to the discussion in Sec. IV B 2, we can choose a cutting surface just above any chosen band; the two Dirac cones cut by this surface contribute an even integer times π , and the system is again axion-trivial.
(3) N d alternates between 1 and 3 as bands are counted. The Wannier band structure is again connected. A cutting surface is again chosen just above one Wannier band; since there will be 1 or 3 Dirac cones cut by this surface, each carrying a Berry phase of π , the system is a strong topological insulator and it is axion-odd.
These conclusions are consistent with the well-known properties of TR-protected topological insulators, as discussed, for example, in Ref. [31].
For the case of scalar particles, (E ) 2 = +1, the Kramers degeneracies are not enforced, and the Wannier band structure will generically be disconnected, in which case θ = 0. In case a connected Wannier band structure is enforced by degeneracies associated with additional symmetries, the Dirac node counting procedure described at the end of Sec. IV B 2 should still apply.

Fu-Kane-Mele model
To illustrate how the trivial, weak-TI, and strong-TI phases manifest themselves in the HW representation, we use the Fu-Kane-Mele (FKM) model [83]. In this model, spinor s-like orbitals interact with spin-dependent and spin-independent hoppings on a diamond lattice of conventional lattice constant a. The Hamiltonian contains two terms, The first term describes spin-independent first-neighbor hoppings, denoted as i j ; the hopping amplitudes are t i j = t 0 for all bonds except those aligned along [111], for which t i j = t 0 + t. The second term, which models the effect of spinorbit coupling, consists of spin-dependent second-neighbor hoppings denoted as i j ; σ = (σ x , σ y , σ z ) are the Pauli matrices, and Trivial insulator, to put the system in a trivial-insulator, weak-TI, or strong-TI phase at half filling [31]. For our numerical tests we pick a trivial insulator with = 2.5, a weak TI with t = −1.0, and a strong TI with t = 1.0. Their Wannier band structures are shown in Fig. 5; in the top panels the wannierization is along (111), and in the bottom panels it is along (001).
As explained earlier, in a trivial insulator the number N d of Kramers down-touchings at the PTRIM is expected to alternate between 0 and 4 irrespective of the wannierization direction, and this is indeed what is observed in the left panels of Fig. 5.
As for the middle panels of Fig. 5, we note that a weak TI can always be thought of as a stack of QSH layers in some direction determined by the weak indices [83]. In this example the weak-TI phase corresponds to a stack of QSH layers along (111), which should not come as a surprise since we have weakened the hoppings along [111]. Hence, when wannierizing along (111) we expect a gapped spectrum with N d alternating between 0 and 4, as seen in Fig. 5(b). However, for other wannierization directions we expect a connected spectrum with N d = 2 for all bands, as seen in Fig. 5(e).
Finally, the behavior of the right panels of Fig. 5 is as predicted for a strong TI: the Wannier band structure is necessarily connected, with N d alternating between one and three.

General considerations
Next we consider the case that a simple inversion I is a symmetry of the system, so that it is denoted as an "axion insulator" if it is axion-odd. This is an example of a zreversing axion-odd symmetry, and the conclusions derived in Sec. IV B 1 apply. To summarize, we found that in the disconnected case, θ = π if and only if there is a boundarycentered group and its Chern number C BC is odd. (The same applies to the origin-centered group.) Unlike TR, the presence of inversion symmetry does not require an axion-odd insulator to have a connected Wannier band structure [34,68], as discussed in Sec. V B. Regardless of whether it does or not, we can frequently deduce the Z 2 index from a node-counting argument, e.g., by focusing on the set of Wannier bands that touch z =−c/2 at¯ . If the number of bands in this group is even, then we place a cutting surface infinitesimally above the average of the central two bands, and we conclude that the Z 2 index is odd if and only if the number of Dirac nodes sliced in this way is odd. If the number of Wannier bands is odd, then we need information about the Chern number of the central band in addition. The same considerations apply to an analysis centered around z = 0. We note in passing that the work of Alexandradinata et al. [32] provides a complementary view of the Wannier (Wilson-loop) band structure for 2D insulators.
The axion Z 2 index can also be obtained from a very different approach based on parity counting. In a pathbreaking work, Turner et al. [41] derived a set of rules for deducing many topological properties of a centrosymmetric crystalline material based on a counting of the odd-parity eigenvalues of the occupied Bloch states at the eight TRIM in the 3DBZ. They showed that the system can be insulating only if the total number of odd-parity states is an even integer, and moreover it is axion-even or axion-odd depending whether this number is of the form 4n or 4n + 2, where n is an integer.
In the Appendix, we carry this analysis over to the Wannier band structure. Here, the original eight TRIM of the 3DBZ project onto four PTRIM in the 2DBZ. At each of these k, the Wannier bands can be classified as being of even or odd parity at z = 0; even or odd parity at z = c/2; or appearing in pairs at ±z. We develop counting rules inherited from the 3D parity counting of the Bloch states. Consistent with the analysis above, we find that the axion Z 2 index can often be determined directly from an inspection of the Wannier band structure, even without evaluating the parities of the HW states. In some cases, however, some parity eigenvalues do need to be evaluated, allowing a determination of whether the boundary-centered group has an odd Chern number as discussed above. The details of this analysis are deferred to the Appendix. TABLE II. Number of odd-parity Bloch states at the eight TRIM for three representative cases: "Trivial" phase with (t, λ, ) = (1.0, 0.1, 3.0); "Axion I" phase with (t, λ, ) = (1.0, 0.1, 0.6); and "Axion II" phase with (t, λ, ) = (−1.0, −0.1, 0.6). The TRIM are labeled as in Ref. [42]. Trivial  0  2  1  3  12  Axion I  0  2  0  4  10  Axion II  0  2  2  2  14

Pyrochlore model
For the case of inversion symmetry we borrow an illustrative class of centrosymmetric tight-binding models discussed in Ref. [34]. These models consist of one pair of spinor basis orbitals per site of the pyrochlore lattice, as might be used to describe the J eff = 1/2 manifold of the pyrochlore iridates. The Hamiltonian reads where the spin indices are suppressed. Just as for the FKM model, the first two terms describe spin-independent and spindependent first-neighbor hoppings, with amplitudes t and λ, respectively. In the latter,ĝ i j is a unit vector normal to both the hopping and the potential gradient directions; see Ref. [34] for details. The third term is a Zeeman term responsible for breaking TR symmetry, with the vectorsn i (i = 1, 2, 3, 4) describing the directions of the magnetic moments on the pyrochlore lattice sites. We choose the all-in-all-out magnetic configuration [34] and consider the model at half filling.
Whether the system is axion-even or axion-odd can easily be determined using the parity criteria for inversionsymmetric insulators [41]. That is, the system is even or odd if the total number of odd-parity Bloch states at the TRIM is of the form 4n or 4n + 2, respectively, where n is an integer. Table II shows how the odd-parity states are distributed among the eight TRIM for three representative choices of parameters (see the caption) corresponding to different regions of the phase diagram. Their total number is 12, 10, and 14, indicating axion-even, axion-odd, and axion-odd topology, respectively. We refer to them below as the trivial, Axion I, and Axion II cases, respectively.
We present the corresponding Wannier bands in Fig. 6. Along (111) the pyrochlore lattice consists of alternating triangular and kagome layers containing one and three atoms, respectively. This is reflected in the Wannier band structure of the trivial phase in Fig. 6(a), since in that phase electrons are localized on the atomic sites. Similarly, along (001) we have a stack of tetragonal layers rotated by 45 • with respect to each other, and in Fig. 6(d) we see that the OC and BC groups have two Dirac nodes each, implying θ = 0.
In the axion-insulator phase, the Wannier band structure can be connected or disconnected depending on the choice of parameters, as illustrated in the middle and right panels of Fig. 6. The two axion phases differ by four odd-parity states (see Table II Table II). II, is another illustration of the kind of fragile protection that was discussed in Sec. V B. For example, it is evident that the addition of a pair of trivial flat Wannier bands near z =±0.25c would convert the hybrid Wannier band structure in Fig. 6(b) into a disconnected one after hybridization with the crossing bands would occur.
In the connected Axion I case, the number of Dirac nodes at the nominal cell boundary is an odd integer, as per Eq. (41). Indeed, in Fig. 6(b) there are five nodes at z = c/2 and three at z = 0, while in Fig. 6(e) there are three nodes at both z = c/2 and z = 0.
For the disconnected band structures in the left and right panels of Fig. 6, we showed in Eq. (39) that θ = 0 or π depending on whether the Chern index of the boundary-centered group is even or odd. While it is possible in principle to compute Berry curvatures and Chern numbers from the bulk Wannier bands [57,60], here for convenience we carry out the analysis based on a slab configuration chosen thick enough that the central Wannier bands are bulk-like. In particular, we construct slabs with a thickness of ten unit cells alonĝ z (111). Slab HW states maximally-localized alongẑ are obtained as eigenstates of the projected position operator P k z P k [55,60,84]. For a sufficiently thick slab, the ones located far from the surfaces are indistinguishable from the bulk HW states. We then organize the corresponding eigenvaluesthe Wannier bands of the slab-into origin-centered and boundary-centered groups, i.e., centered at integer and halfinteger z /c , respectively. Finally, we calculate explicitly the Chern index of each group of Wannier bands using Eq. (39). In Fig. 7 we confirm that in the disconnected Axion II phase the Chern index of the boundary-centered groups is an odd integer (C = −1), while in the trivial phase it is even (C = 0).
For the Axion II phase, we could have reached the same conclusion from either Fig. 6(c) or 6(f) by counting Dirac nodes, since the origin centered and boundary centered groups have two Wannier bands each. For the trivial phase this strategy cannot be used for the case of 3+1 bands in Fig. 6(a), but it does does apply to the (001) wannierization with 2+2 bands in Fig. 6(d), confirming the trivial topology.

General considerations
In Sec. IV A we considered two types of mirrors relative to a given wannierization directionẑ. Whenẑ lies on the mirror plane, the mirror is of type M d , and is a symmorphic z-preserving operation as discussed in Sec. IV B 2. Whenẑ is perpendicular to the mirror plane, we have an M z mirror reflection, a z-reversing axion-odd symmetry of the kind presented in Sec. IV B 1.
The case of M d is similar to that of TR, while the case of M z is similar to that of inversion, but with an important difference. In the case of inversion, the Bloch states only have well-defined parities at the eight TRIM, which are the reciprocal-space points that are carried onto themselves by inversion. By the same token, the parity analysis of the HW FIG. 7. Chern index for connected groups of slab Wannier bands of the pyrochlore model, plotted as a function of their average positions z . Red and blue markers correspond to groups with integer and half-integer z /c , respectively. We consider slabs with a thickness of 10 unit cells along (111) in (a) the trivial phase of Fig. 6(a), and (b) the axion II phase of Fig. 6(c); in the axion phase, the top layer was "exfoliated" to remove unwanted nontopological surface states. In the trivial phase, both the groups with integer z /c -one band-and those with half-integer z /c -three bands-have zero Chern index. In the axion II phase both groups comprise two bands each, and far from the surfaces their Chern indices are −1 and +1, respectively. states discussed above is only applicable at the four PTRIM. As a result, Dirac touchings between Wannier bands may occur at z = 0 or z = c/2 at the four PTRIM but are generically absent elsewhere unless enforced by additional symmetries.
The M z operation, however, acts as the identity on k x and k y , so that all k in the 2DBZ follow "mirror parity" rules similar to those discussed above for the PTRIM under inversion. That is, at any k, the Wannier bands can be decomposed into those centered at z = 0 with even or odd parity labels; those centered at z = c/2 with even or odd parity labels about their center; and mirror-symmetric pairs of Wannier bands at ±z. A consequence is that Wannier bands can be pinned exactly at z = 0 or z = c/2 over the entire 2DBZ.
Moreover, in the presence of mirror symmetry, a new kind of topological invariant known as the "mirror Chern number" comes into play [85][86][87]. The mirror Chern number on the k z = 0 symmetry plane is normally defined in reciprocal space in terms of the difference of Chern indices of the even-and odd-parity mirror subspaces on this plane, whose points are left invariant under M z . Except for certain centered lattices [88], a second mirror Chern number can be defined in the same way on the k z = π/c plane as well.
In a separate collaboration involving some of us, we have investigated the Wannier bands in the presence of M z symmetry in some detail, clarifying the generic behaviors that are expected, and discussing the rules for deducing not only the axion Z 2 index, but also the mirror Chern number, from the Wannier band structure. The work will be published elsewhere [89].

Alternating Haldane model
To illustrate the case of mirror symmetry we use a spinless tight-binding model introduced in Ref. [60], which we refer to as the "alternating Haldane model." This model consists of layers of the Haldane model [1] stacked along the (001) direction, with two layers per cell. The Hamiltonian can be expressed as where the first term describes isolated layers, the second term couples adjacent layers, and "H.c." stands for "Hermitian conjugate." The Hamiltonian is constructed in such a way that in the decoupled limit adjacent layers have opposite Chern numbers. Specifically, factor in the third term is responsible for reversing the Chern number between adjacent layers. The first term also contains a (−1) p factor that reverses the on-site energies; as a result, adjacent layers are not simply time-reversed as in the case of "antiferromagnetic TIs" with {E |c/2} symmetry [43].
The above model was used in Ref. [60] to pump a quantum of axion coupling during a slow cyclic evolution of the Hamiltonian along the path t 1 = 4.0, parameterized by the angle φ. Here we are interested in the configurations φ = 0, π where t 3,p becomes independent of p, and as a result the model acquires mirror symmetry about the planes of the layers. At φ = 0 the half-filled system is axioneven, and at φ = π it is axion-odd [60]. We first wannierize the axion-even and odd systems along the direction (001) normal to the mirror plane, as done in Appendix A of Ref. [60]. Within this setting the mirror operation becomes M z , and the resulting Wannier bands are shown in the top panels of Fig. 8. The earlier discussion indicates that a consequence of M z symmetry is that entire bands can be pinned exactly at z = 0 or z = c/2 over the full 2DBZ. As is clear from the figure, this is the case here. In fact, with just two occupied bands, our model is so simple that this is all we have, with one Wannier band pinned at z = 0 and the other at z = c/2. Just as for the inversion-symmetric case, the axion angle θ is 0 or π depending on whether the Chern index of the boundary-centered group is even or odd. This is confirmed by FIG. 9. Chern index of each isolated slab Wannier band of the alternating Haldane model, plotted as a function of its average position z. We consider slabs with a thickness of 10 unit cells along (001) in (a) the trivial phase of Fig. 8(a), and (b) the axion phase of Fig. 8(b). In the trivial phase each Wannier band carries a zero Chern index, while in the axion-odd phase the Chern index alternates between −1 and +1.
inspection of Fig. 9, where we show the Chern numbers of the Wannier bands in a slab geometry.
Finally, we consider the case of wannierization along an axis lying in the mirror plane, specifically, along (010), which we denote as z . Now the mirror operation is of type M d , and the calculated bulk Wannier bands are displayed in the bottom panels of Fig. 8. Figure 8(c) shows a disconnected band structure, and consulting Sec. IV B 2 we conclude that the system is axion-trivial. Conversely, in the axion-odd phase the Wannier bands are now required to flow, each with an odd number N d of down-touchings. This is indeed what happens in Fig. 8(d), where N d = 1 for all bands.

General considerations
A Z 2 classification of 3D insulators with glide mirror symmetry was introduced in Refs. [90] and [91], and was later argued to be equivalent to the axion Z 2 classification [88]. Interestingly, glide mirror symmetry realizes all three types of axion-odd space-group symmetries, depending on the choice of wannierization direction. Ifẑ is chosen as the normal direction to the reflection plane, then the glide operation becomes {M z |τ ⊥ }; recalling that our classification discards in-plane fractional translations τ ⊥ , here the quantizing symmetry is simply M z , a z-reversing operation, as in Sec. IV B 1. Instead, we can chooseẑ to lie in the reflection plane and be perpendicular to the half-lattice translation. In this frame the glide operation is {M d |τ ⊥ }; again discarding τ ⊥ , the axionodd operation becomes M d , a symmorphic z-preserving operation, per Sec. IV B 2. Finally, ifẑ is chosen along the halflattice translation, then the glide operation becomes {M d |c/2}, a nonsymmorphic z-preserving operation, corresponding to Sec. IV B 3.
In the numerical tests below we shall consider the latter nonsymmorphic scenario, which we have not encountered in our previous examples, as well as the z-reversing case. We will omit the symmorphic z-preserving configuration, since we have already encountered this kind of M d symmetry in our study of the alternating Haldane model.

Model
As an example, we take the four-band tight-binding model at half filling described in Ref. [90] (an equivalent model was also used in Ref. [88]). The simple orthorhombic a×b×c unit cell contains two pairs of orbitals with reduced coordinates (0,0,0) and (0.5,0,0), and the Hamiltonian is expressed in k space as The Pauli matrices τ can be considered to represent either spin or orbital degrees of freedom, and the matrices ρ are associated with sublattice degrees of freedom. The above model is invariant under the glide operation {M z |a/2} consisting of a reflection about the z = 0 plane followed by a half-lattice translation alongx. Because it is so simple, the model has two additional axion-odd symmetries, {M y |a/2} and M x , whose presence obscures the role of {M z |a/2} in quantizing θ . For example, the mirror M x forces the (100)-oriented Wannier bands to be completely flat. To break these extra symmetries, we displace the orbitals to (0,0,0.1) and (0.5, 0, −0.1), and include an additional term in the Hamiltonian, chosen weak enough so as not to cause any band inversion. At half filling, the two lowest bands are filled. The model has two adjustable parameters; we set φ = 0.4, and vary m to access different phases. Figure 10 shows the Wannier band structures calculated with m = 3.5 (trivial phase) and m = 2.0 (axion phase), for two different wannierization directions. In the upper panels the wannierization is along the half-lattice translation (along  x), and in the lower panels it is normal to the reflection plane (alongẑ). (For this model, instead of reorienting the axes so that the wannierization direction is always along z as in the previous examples, we have chosen to keep the axes fixed.) In the upper panels of Fig. 10, the glide operation is {M d |a/2} in the notation of Table I. This is a nonsymmorphic operation that preserves the wannierization direction, of the type discussed in Sec. IV B 3. In axion-odd phases protected by such symmetries, the Wannier spectrum is not required to flow. (This is in contrast to axion-odd phases protected by a simple mirror M d , which must exhibit flow as illustrated in Fig. 8(d) for the alternating Haldane model.) Indeed, the spectrum of the present model is gapped not only in the trivial phase of Fig. 10(a), but also in the axion-odd phase of Fig. 10(b). We can therefore use Eq. (44) to determine the axion index, taking for C LG the Chern number of either band. The Chern numbers of the two bands are (0,0) in Fig. 10(a) and (−1, +1) in Fig. 10(b), leading to θ = 0 and θ = π, respectively.
Consider now the lower panels of Fig. 10, where the glide operation returns to {M z |a/2} in the notation of Table I. This is a z-reversing operation, like spatial inversion or a simple mirror M z , so we refer to Sec. IV B 1. Now the Wannier spectrum is gapped in the trivial phase of Fig. 10(c), but connected in the axion phase of Fig. 10(d). (This is different from the alternating Haldane model with M z symmetry, which exhibited gapped Wannier spectra in both phases.) In the trivial phase of Fig. 10(c), the two isolated bands have vanishing Chern numbers, so that Eq. (39) gives θ = 0. For the axion phase of Fig. 10(d), which exhibits "fragile" flow, we resort to Eq. (41), and since there are single Dirac cones at integer and half-integer values of z n /c, we conclude that θ = π . At variance with our previous examples of connected band structures, here the Dirac nodes are located at generic low-symmetry points in the 2DBZ, not at high-symmetry points or lines.

VII. DISCUSSION AND OPEN QUESTIONS
In this work we have surveyed the symmetry operations that can protect the axion Z 2 index in 3D insulators, inducing a topological classification into those that are trivial (axioneven, θ = 0) and those that are topological (axion-odd, θ = π ). When considered in the context of the hybrid Wannier construction performed along a given crystallographic direction, the relevant axion-odd symmetries are classified into those that involve reversal along the wannierization direction and those that preserve it, with the latter case subdivided into those that do not involve an additional fractional translation along this direction and those that do.
We then systematically explored the consequences of such symmetries for the Wannier bands, and conversely, the means by which the axion Z 2 index can often be deduced from characteristics of the Wannier band structure. In particular, we described the conditions under which a counting of Dirac touchings between Wannier bands, and/or the calculation of the Chern number of some central Wannier bands, allows for a direct determination of the axion index. We then confirmed our conclusions and illustrated their consequences for several paradigmatic cases, namely time reversal, inversion, mirror, and glide-mirror symmetries.
For some of the cases considered here (including time reversal, time reversal composed with either a half-lattice translation [43] or a twofold rotation [90], inversion [40,41], and glide-mirror symmetry [88,90,91]), the Z 2 classification has been previously discussed. However, this was often done on a case-by-case basis by defining the Z 2 index in a way that was tailored to each specific symmetry, and it was not easily generalizable to other axion-odd symmetries. Our work systematizes these considerations, placing them in the context of a complete classification of axion-odd symmetries and their influences on Wannier band structures.
Our analysis shows that the Wannier bands are sometimes topologically required to display "flow," i.e., to have touchings between all neighboring bands. Even when not required, such touchings are often present, usually resulting from Kramers or symmetry-induced degeneracies at high-symmetry points in the projected 2DBZ. To carry forward our analysis, we typically assumed that the only type of touchings between Wannier bands were point Dirac nodes, i.e., displaying a linear splitting with k away from the nodal point. One interesting avenue for future investigation is to consider whether there are some symmetries, or combinations of symmetries, that can give rise to higher-order point nodes or nodal lines or loops. In such cases, our prescriptions for deducing the axion Z 2 index from the nodal structure of the Wannier bands would need to be reconsidered. (An extreme example is the case of a simple M z mirror, encountered in Sec. VI C, which can sometimes produce two or more entirely flat degenerate Wannier bands.) Magnetic axion insulators offer great promise in spintronic applications [92], ranging from the high-temperature quantum anomalous Hall effect to chiral Majorana fermions [93]. Unfortunately, and in spite of recent progress, axion-odd insulators protected by symmetries other than TR still elude experimental observation. As mentioned in the Introduction, it appears that MnBi 2 Te 4 should be a realization of an antiferromagnetic topological insulator [43] with an axion-odd state protected by both inversion and time-reversed fractional translational symmetries, but the most recent experiments do not observe the expected gapping of the surface Dirac cone [44][45][46][47][48][49]. Nevertheless, the interest of the research community has been aroused, and searches are currently under way for magnetic axion-odd insulators in the broader family (MnBi 2 Te 4 )(Bi 2 Te 3 ) m [94][95][96][97][98], as well as in new candidates such as EuIn 2 As 2 [99][100][101] and NpBi [25].
A mainstream search strategy is to focus on topological transitions driven by band inversion at high-symmetry points in the Brillouin zone, with strong spin-orbit coupling playing a role in driving the band inversion and/or gapping the system after the band inversion has occurred. More radical scenarios that do not involve spin-orbit coupling, or that depend on the spontaneous formation of highly correlated many-body states, might also prove viable. Recent progress in theoretical screening of databases of magnetic materials to identify those with potential topological properties is also encouraging [22][23][24][25].
Our analysis serves as an alternative tool for the identification and classification of axion-odd topological insulators, providing an approach that is complementary to methods based on K-theory classifications, or on symmetry indicators of elementary band representations. The calculation of Wannier bands is easily automated, and is therefore well suited for high-throughput computational searches for new axion-odd insulators. By providing a pictorial perspective, our approach may help provide intuitive clarity about the consequences of the symmetries that are present, especially in connection with the bulk-boundary correspondence. We are hopeful that these hybrid-Wannier based tools will facilitate the discovery of material realizations of novel topological states exhibiting a quantized CS axion coupling and half-integer surface quantum anomalous Hall conductivity.

APPENDIX A: INVERSION SYMMETRY AND PARITY COUNTING
In an important paper by Turner et al. [41] derived a set of rules for extracting information about the topological state of a 3D centrosymmetric crystal from a counting of the oddparity occupied eigenstates at the TRIM. Specifically, they showed that the total number of odd-parity states summed over the occupied bands at the eight TRIM must be even in an insulator, and that the axion Z 2 index of the system is even or odd if the total number of such odd-parity states is of the form 4n or 4n + 2, respectively, where n is an integer. We also call attention to Refs. [40] and [32] for alternative viewpoints on the topology of centrosymmetric insulators.
Here we develop a connection between the numbers of odd-parity Bloch states at the eight TRIM and the corresponding numbers of odd-parity HW states at the four PTRIM, and use this to develop criteria for deciding the axion Z 2 index based on an inspection of the Wannier band structure. The results are consistent with those of Sec. VI B, thus providing an alternative derivation of the results presented there.
As a preliminary exercise, we begin by providing a pedagogical development of some important rules that underlie the theory of Ref. [41], beginning with a 1D centrosymmetric insulator and working our way to two and three dimensions.

Parity counting for 1D insulators
Consider a 1D centrosymmetric insulator of lattice constant c, with its two inversion centers located at the origin z = 0, denoted as location A, and at z = c/2, denoted as B. In reciprocal space, there are J occupied bands at each of the two TRIM at (k = 0) and X (k = π ; reduced wave vectors are used here). Let N o μ be the number of odd-parity Bloch eigenstates at μ = or μ = X , with N o = μ N o μ summed over and X , where the parities of the Bloch states are defined relative to the chosen origin at A. Now consider the maximally localized Wannier centers for this 1D crystal and the corresponding Wannier functions. In view of the inversion symmetry, the general description is that we find M o A and M e A Wannier states centered at A with odd and even parity, respectively, about A; M o B and M e B states centered at B with odd and even parity, respectively, about B; and M P pairs of Wannier states centered at ±z, distinct from both A and B. The total number of occupied states is thus We can also write down the total number N o of odd-parity Bloch states summed over the two TRIM as This follows because an odd-parity Wannier state at A contributes odd-parity Bloch states at both and X ; an evenparity one at A contributes none at or X ; an odd-parity one at B contributes an odd-parity state only at ; an even-parity one at B contributes an odd-parity state only at X ; and the pair contributes one odd-parity state at each of and X . From Eqs. (A1) and (A2) it follows trivially that These considerations provide an elementary derivation of a 1D polarization-parity relation. A 1D centrosymmetric insulator has a polarization Z 2 index, which we shall refer to as Z pol 2 , defined as 0 or 1 if the total Berry phase of the occupied Bloch bands is 0 or π , or the electronic contribution to the polarization is 0 or e/2, respectively. Each Wannier center at location B contributes one unit to Z pol 2 , while those at A do not contribute, and neither do the pairs. Thus, That is, the Bloch parity sum determines whether the 1D electronic polarization is trivial (Wannier charge is centered at A) or not (centered at B). Put another way, the polarization Z 2 index and the parity Z 2 index are the same in 1D.

Parity counting for 2D insulators
The authors of Ref. [41] also derived a Chern-parity relation for a 2D insulator. Let N o 2D be the number of odd-parity occupied Bloch eigenstates obtained by summing over the four TRIM of the 2DBZ. They showed that where C is the total Chern number of the 2D insulator (again summed over all occupied bands). This follows from the considerations of the previous subsection. Let the system extend in the y-z plane, and we wannierize along z as a function of k y . The polarization P z (k y ) is well defined and varies smoothly for intermediate k y , only taking quantized values {0, e/2} at 0 and π . As is well known, the Chern number C describes the pumping of polarization P z as k y is taken adiabatically from 0 to 2π . A change of P z by −e, or equivalently a change of the Berry phase by 2π , corresponds to a unit Chern number. In view of the inversion symmetry, half of the pumping of P z occurs as k y traverses from 0 to π , and the other half between π and 2π . Thus, the value of Z pol 2 changes from k y = 0 to k y = π if and only if the Chern number is odd.

155130-21
However, the earlier 1D analysis shows that the value of Z pol 2 at k y = 0 (k y = π ) corresponds to the number of odd-parity occupied eigenstates summed over the two TRIM projecting onto k y = 0 (k y = π ). Thus, it follows that the value of Z pol 2 changes from k y = 0 to k y = π , and hence the Chern number is odd, if and only if N o 2D , the sum of odd Bloch parities over all four TRIM, is odd. This proves Eq. (A5).
Our considerations will be limited to the case of insulators with C = 0, in which case Eq. (A5) implies that N o 2D must be even. This fact will be useful later.

Parity counting for 3D insulators
We now apply the rules derived above to the consideration of a 3D centrosymmetric insulating crystal. Let C z be the Chern index on the plane k z = 0 or k z = π ; these are equal since for an insulator C z must be the same at at all k z . Now N o 3D , the total number of odd-parity Bloch states at the eight TRIM, is just the sum of those at k z = 0 and k z = π , and using Eq. (A5), we have that This provides a rederivation of one of the important results given in Ref. [41]. As those authors pointed out, this implies that any system with N o 3D odd cannot be insulating; it must have nodes of degeneracy between the nominal valence and conduction bands, as in a Weyl semimetal. (In this work we assume vanishing Chern numbers, but the above result is general.) The authors of Ref. [41] derived another result that is central here, namely, that the axion Z 2 index, which we now denote as Z axi 2 in analogy with Z pol 2 in 1D, is given by In other words, if N o 3D is of the form 4n for integer n, then the system is axion-even, while if N o 3D = 4n + 2, then it is axionodd. This is a rather deep result which we shall not attempt to rederive here. Our goal now is to use this result to determine the axion Z 2 index from an inspection of the Wannier band structure.
We choose to wannierize along z and retain the (k x , k y ) Bloch labels. Let μ label the four PTRIM in the 2DBZ in (k x , k y ) space; each of these corresponds to a pair of the 3D TRIM separated by π along k z . For a given PTRIM μ, the string of k points extending from k z = 0 to 2π describes a 1D centrosymmetric insulator. We can then apply the analysis of Appendix Sec. A 1, rewriting Eq. (A3) as where N o μ is the total number of odd-parity Bloch states at the two associated TRIM. This motivates us to define as the excess number of odd over even-parity HW states pinned at location A, summed over the four PTRIM. Then, noting that J is an integer, Eq. (A8) implies that Because of Eq. (A6), we know that a 3D insulator must have an even value of N o 3D , so m A /2 is an integer. Combining Eq. (A10) with Eq. (A7), we arrive at the first major result of this Appendix, namely, that the axion Z 2 index is simply given by That is, the system is axion-odd if and only if m A /2 is an odd integer.
Again we point out that the choice of origin at A was arbitrary; placing the origin at B leads to the conclusion that the axion Z 2 index is also determined by m B = 4 μ=1 (M o μB − M e μB ), which must equal m A mod 4. At one level, Eq. (A11) provides an answer to the question of how to extract the axion Z 2 index within the HW representation. One simply looks for Wannier bands pinned at z = 0 (or z = c/2) at the PTRIM, computes the parities of the corresponding HW states, and sums the parities over the four PTRIM. Half this number gives Z axi 2 . Unfortunately, this approach does require an extraction of the parities of these HW states, so the needed information does not appear to be contained in the Wannier band structure itself. However, we show next that this is not always true; a simple count of the number of such pinned states, without a knowledge of their parities, is sometimes sufficient.

Node counting for 3D insulators
Again we restrict ourselves to 3D centrosymmetric insulators with all zero Chern indices. It is useful to repeat the definition of Eq. (A9) for each PTRIM individually, and also to define the corresponding count p μA of the total number of Wannier bands pinned at A, via With this notation, Eq. (A8) becomes Now consider two neighboring PTRIM μ and μ ; these define a 2D centrosymmetric system corresponding to the k-space plane in 3D passing through the four TRIM projecting onto these two PTRIM. Following the discussion at the end of Sec. A 2 of this Appendix, the zero-Chern assumption requires that the total number of odd-parity Bloch states on these four TRIM, i.e., N o μ + N o μ , must be even. Summing Eq. (A14) over μ and μ and noting that J is an integer, it follows that m μA + m μ A is even. But Eqs. (A12) and (A13) differ by an even integer, so p μA + p μ A is also even. That is, p μA and p μ A are either both even or both be odd. Continuing to apply this principle to other pairs of PTRIM, we come to an important constraint, namely that the number p μA of Wannier bands pinned at A is either even at all four PTRIM, or odd at all four PTRIM. We consider these two cases in turn.
We first assume that p μA is even at all four PTRIM. To make further progress, we assume henceforth that at any given PTRIM, the HW states pinned at A are all of the same parity. We call this the "uniform parity" assumption. This is a rather natural assumption, since states of even and odd parity about A will generically have a nonzero z matrix element between them. The maximal localization procedure corresponds to diagonalizing the z operator in the occupied band space [55,84], and this will generically result in hybridization and splitting to form a pair of HW states at ±z, so that these states will no longer be pinned at A. If this process repeats until the only states left at A are all of the same parity, then our uniform parity assumption holds [102]. This procedure may break down if symmetries other than inversion are present (e.g., the even and odd parity states at one of the TRIM fail to split because they belong to different irreps under rotation about the z axis), or in case the system has been fine-tuned to force the splitting to vanish. Therefore, we treat it simply as an assumption in the following. Now with this assumption, |m μA | = p μA and both are even, so m μA = p μA mod 4, and defining p A = 4 μ=1 p μA , Eq. (A11) can be replaced by Thus, the axion Z 2 index is obtained just by counting the number of HW states at the PTRIM that are pinned at A! Again, the choice of origin at A was arbitrary, so the same applies to the counting of HW states pinned at B. Assuming these PTRIM degeneracies are point nodes of contact between Wannier bands, p A /2 is just the number of Dirac crossings at the PTRIM at A, and the system is axion-odd if and only if this number of Dirac crossings is odd. The same applies to the Dirac nodes at B. This is in perfect agreement with the conclusions of Sec. IV B 1 as expressed in Eq. (41) and the discussion following it.
In case the number p μA of Wannier bands pinned at A is odd at each of the PTRIM, and again assuming uniform parity of the HW states at each of these four points individually, we can use a similar analysis to evaluate the contribution from all except the central Wannier band; the contribution of these is again given by Eq. (A15). To obtain the contribution of the central Wannier band, we need to know the sum of the parities over the four PTRIM for this band. If we have access to the parities, then we can do this directly. Alternatively, a computation of the Chern number of this band will also suffice, since Eq. (A5) tells us that the Chern index is odd if and only if the parity sum is odd. In this way we again arrive at Eq. (42), which involves the computation one Chern number in addition to a counting of the Dirac cones passing through z = 0.
Once again we note that the same analysis can equally well be applied at location B, and the choice of inversion center to use for the analysis can be made on the basis of convenience. In any case, we have succeeded in arriving at the same conclusions as were presented in the main text at the end of Sec. IV B 1, but from a very different point of view.