{'fraction_non_alphanumeric': 0.036345865641493454, 'fraction_numerical': 0.008770443734931908, 'mean_word_length': 4.428025748037066, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Competing Adiabatic Thouless Pumps in Enlarged Parameter Spaces

The transfer of conserved charges through insulating matter via smooth deformations of the Hamiltonian is known as quantum adiabatic, or Thouless, pumping. Central to this phenomenon are Hamiltonians whose insulating gap is controlled by a multi-dimensional (usually two-dimensional) parameter space in which paths can be defined for adiabatic changes in the Hamiltonian, i.e., without closing the gap. Here, we extend the concept of Thouless pumps of band insulators by considering a larger, three-dimensional parameter space. We show that the connectivity of this parameter space is crucial for defining quantum pumps, demonstrating that, as opposed to the conventional two-dimensional case, pumped quantities depend not only on the initial and final points of Hamiltonian evolution but also on the class of the chosen path and preserved symmetries. As such, we distinguish the scenarios of closed/open paths of Hamiltonian evolution, finding that different closed cycles can lead to the pumping of different quantum numbers, and that different open paths may point to distinct scenarios for surface physics. As explicit examples, we consider models similar to simple models used to describe topological insulators, but with doubled degrees of freedom compared to a minimal topological insulator model. The extra fermionic flavors from doubling allow for extra gapping terms/adiabatic parameters - besides the usual topological mass which preserves the topology-protecting discrete symmetries - generating an enlarged adiabatic parameter-space. We consider cases in one and three \emph{spatial} dimensions, and our results in three dimensions may be realized in the context of crystalline topological insulators, as we briefly discuss.


I. INTRODUCTION
Insulating matter is often thought of as a transportinert environment. It was thus quite surprising that charge currents could be established in an insulator at zero temperature and the concept of 'adiabatic pumping', as described by Thouless, now pervades the literature 1-3 . Such pumps characterize whether a periodic variation of a set of Hamiltonian parameters of a gapped system leads to the transport of (quantized) charge, that is conserved due to the presence of a symmetry.
To illustrate, the conventional example considers a 1D insulator (with lattice constant a) under a slow, adiabatic translation -parametrized by a phase shift δφ(t) -of the underlying periodic potential: V (x) → V (x + δφ(t)a/2π) 1 . During one cycle, δφ(t) : 0 → 2π, an integer amount of electron charge is pumped through the 1D system. This integer is equal to the 2D Chern number of the 1D Bloch Hamiltonian parameterized by the 1D momentum k and δφ. Other examples consist of spin pumps, defined by the spin Chern number when some component of the electron spin is conserved, and discrete Z 2 pumps (e.g., Kramers' pair, fermion parity) defined by a bulk Z 2 invariant 4-6 . Higher dimensional generalizations of the charge pump are also possible, the first of which being the magneto-electric polarizability pump, fixed by the second Chern number [7][8][9] .
The ubiquitous relation between quantum pumping and bulk Chern numbers connects the notion of adiabatic transport in gapped systems to the concept of topological band insulators. Such topological phases constitute a class of gapped systems which generically display gap-less surface states, when in the presence of certain symmetries such as time-reversal. These surface states are inherently robust to disorder, persisting even in the presence of interactions as long as the relevant symmetries are not broken (explicitly or spontaneously) and no surface topological order is developed. This class of systems displays a short-ranged profile of spatial entanglement, and constitutes what are now known as 'symmetry protected topological phases' (SPTs) 10,11 .
In addition to adiabatic pumping during closed cycles, topological responses (electromagnetic, thermal, etc.) of SPTs may be derived by considering adiabatic transformations of the gapped Hamiltonian on special open paths in parameter space. These paths consist of adiabatic (gapped) interpolations between trivial and topological symmetry-preserving reference Hamiltonians 7,12 . In order to adiabatically connect the two end points, somewhere (or everywhere) in-between, the protecting symmetry must be broken. For example, for a 1D insulator with inversion symmetry, a trivial phase has a vanishing charge polarization, modulo an integer charge. Deforming the corresponding Hamiltonian parameters into the topological regime changes the polarization to its topological value of e/2 (again modulo integer charge). Comparing the initial and final states, the difference is a halfinteger polarization. A symmetry preserving interpolation between the Hamiltonians of an inversion symmetric system in different topological phases generically implies the closure of the system gap somewhere during the interpolation. This gapless region is singular in the sense that it renders the notion of 'adiabatic transformation' ill defined. An adiabatic/gapped interpolation between the two reference points, therefore, demands the introduction of parameters which break the protecting symmetry of the Hamiltonian (inversion in the present example). If the interpolation is continued onward to cycle back to the initial trivial state this forms a closed path and the polarization in the final state can only differ from the polarization of the initial state by an integer number of electrons; hence, indicating a quantized number of electrons pumped.
From the previous discussion, the adjective 'adiabatic' should be understood in the present context as implying that for all allowed values of the parameters defining the Hamiltonian, the gap should remain open. The description of the interesting pumping path discussed above demands a minimal set of two parameters, the original symmetry preserving gap term that switches between the topological phases, and a symmetry breaking gap term. Interpolations define one-dimensional trajectories in this non-simply connected two-dimensional space (the point with both gap terms vanishing must be removed as the system is gapless there).
In this article we focus on generalizing the notion of adiabatic quantum transport and responses to enlarged nodal parameter spaces. We consider, concretely, parameter spaces in one higher dimension (i.e., three parameters instead of two.) that are non-simply connected, and allow for different classes of trajectories which cannot be continuously deformed into one another. The specific case we consider will be a three-dimensional parameter space with several nodal/gapless lines. Interestingly, even when paths connect the same two reference Hamiltonians, the different possible classes can lead to different response properties, and ultimately the pumping of different conserved charges. This should be contrasted with the aforementioned cases -the difference in the electromagnetic response (e.g., difference in polarization) defined for an open path with two symmetry-preserving end-points depends only on the end-points themselves, but not on the choice of path (modulo integer numbers of electrons). Accordingly, all adiabatic cycles of a given Hamiltonian would pump the same quantum numbers. This result depends on the mathematical form of the topological invariant and the topology of the pumping parameter space itself. Here we will show that changing the topology of the parameter space can lead to interesting new features.
While our results will apply in general, for concreteness, we focus on time-reversal invariant gapped systems in 1D and 3D (spatial dimensions), with U (1) charge conservation symmetry. The models we consider are doubled versions of the usual topological insulator minimal Dirac models in those dimensions. Because of the doubled degrees of freedom, we find that the natural pumping parameter spaces in these models can display gapless singular lines due to a competition between incompatible gapping parameters (mass terms). The resulting pumping transport and responses are then found to depend on which singularities are encircled during the adiabatic cycles, i.e., they depend on the path and not just the end points. We approach our analysis by several different methods, all with matching results. Of note, for one of our methods we employ the Mañes-Bardeen form of the Wess-Zumino action [13][14][15] to compute the action change upon adiabatic transformations of the Hamiltonian, which is a powerful technique that has not been applied in this context.
The models on which we focus are essentially doubled versions of 1D and 3D Z 2 topological insulators. As such they are trivial insulators according to the 10-fold classification table 16,17 . Indeed, the doubling allows extra possible mass terms that can be chosen to gap any surface states without breaking the protecting discrete symmetries. Thus, while the pumping processes we consider are stable, the topological phenomena and connection to the related electromagnetic responses are not reliable since the parent topological insulator phases are rendered trivial. The response physics of the models we describe, however, may have relevance in the context of crystalline topological insulators 18,19 , in which case the instabilities will be removed by requiring the preservation of, e.g., mirror symmetries.
The paper is organized as follows. To illustrate the basic ideas, we start in Section II with a review of adiabatic quantum pumping in 1D, considering the example of the spinless Su-Schrieffer-Heeger model 1,20 . In Section III, we proceed to explain the important concepts for Thouless pumping processes in our higher dimensional parameter spaces. Section IV is devoted to our first example in 1D. We consider a spinfull Su-Schrieffer-Heeger chain and demonstrate the details of the generalized adiabatic pumps through a four-fold analysis. First by a perturbative continuum field theory computation, second in a microscopic lattice picture and, third, in a nonperturbative approach using the Mañes-Bardeen formula for the anomaly. Finally, we employ a bosonization approach to further solidify the results.
Next, in Section V, we consider a second example of a 3D time-reversal invariant insulator. In this case we again consider the perturbative field theory computation first. We then proceed to understand the problem from the point of view of the surface state properties in the presence of a magnetic field. Finally, although the bosonization picture is not available, we also consider a non-perturbative approach from the Mañes-Bardeen form of the chiral anomaly. We then conclude with a discussion of future directions. In the appendices we add some details of the perturbative calculations.

II. ADIABATIC QUANTUM PUMPING IN 1D
We begin our work with an in-depth discussion and review of the 1D insulator example described in the introduction. We take a 1D insulator coupled to external, adiabatic perturbations represented by a set of n parameters {θ i } which enter the Bloch Hamiltonian H(k, {θ i }), parametrized by the crystal momentum k for translationally invariant systems. H(k) must include a minimum of two bands if it is to represent an insulator; for all of the models in this article there will be an even number of bands.
The adiabatic deformation condition implies that for all values of the parameters in the set {θ i }, the energy gap between these bands remains open. In this framework, if a closed path is traversed in the {θ i } parameter space, then a quantized amount of electric charge will be pumped from one side of the 1D wire to the other 1 . While this statement applies in quite general settings, it is important to note in passing that it fails if electric charge is not conserved during the adiabatic process. For now we assume strict charge conservation, and as such, we would discover that during a cyclic adiabatic process, an integer number of electric charges (possibly zero) would be transferred through the sample.
For most insulators, and most choices of adiabatic perturbation cycles, the amount of transferred charge vanishes. To find cases when the charge transferred is nonzero we need to define a few quantities: (i) the periodic part of the Bloch functions |u α (k, where α labels the band indices and ∂ 0 = ∂ k , ∂ i = ∂ θi ; and, finally, (iii) the first Chern number where the k integral is over the Brillouin zone (BZ). The θ i integral is over the 1-dimensional curve C traversed during the adiabatic process in the n-dimensional {θ i } space, and F αβ µν is the Berry curvature of the occupied bands. Only in the special case when the Chern number (which is an integer by definition) is non-zero is there a finite charge pumping equal to ∆Q = eC 1 .
To explicitly illustrate charge pumping we will use the canonical model, i.e., the spinless Su-Schrieffer-Heeger model 20 . This is a non-interacting dimerized chain with two atoms per unit cell labeled by A/B. The Hamiltonian for the electrons in such a lattice can be written and can be supplemented by an onsite energy term for each atom The Bloch Hamiltonian of H 1 + H ∆ (i.e., the Rice-Mele model 21 ) is in a basis (c Ak , c Bk ) T . Notice that the onsite energy Hamiltonian H ∆ breaks both inversion (P=τ x ) and charge conjugation (C=τ z K, with K the complex conjugation operator) symmetries, and that both of these symmetries quantize the polarization and must be broken in order for one to continuously pump charge. Now, suppose that we have some external control over the parameters of this model which we parameterize with an angular phase θ as a curve in a two-dimensional space (t, ∆): t ≡ m cos θ and ∆ ≡ m sin θ. We will adjust θ such that the perturbations are always adiabatic. In principle, this angle can vary as a function of position and/or time, and the system will respond accordingly. Linear response theory (the Kubo formula) then dictates that the corresponding current density is equal to the time derivative of the charge polarization of the 1D system, Meanwhile, by the continuity equation, the charge density becomes The charge polarization for our system is given by the solid angle subtended by the curve d (k) = t cos k + m cos θ, t sin k, m sin θ which, for m t , gives simply P (θ) eθ/2π. In fact, in this limit it is easy to see that as θ → θ + 2πn then n charges are pumped.
Let us now review the connection to the electromagnetic response of 1D inversion 22,23 symmetric insulators 7 . We have seen from the explicit calculations above that charge density is bound to spatial variations of θ, and charge currents flow in response to a time-dependent θ. This topological response is identical to the Goldstone-Wilczek response 24 and is captured by a θ-term effective action where F µν is the electromagnetic field-strength tensor, and E x is the x-component of the electric field. To see how the response is encoded, one takes the functional derivative We recognize these equations from 1D electromagnetism, and we can identify eθ(x, t)/2π as the electric charge polarization of the 1D insulator. We would like to compare the difference in the electromagnetic response of a trivial insulator and a topological insulator. A trivial insulator has an integer charge polarization and hence, for a homogeneous system, θ = 2πn for some integer n. On the other hand, a topological insulator will have a contribution from a half-integer polarization: θ = 2π(q + 1/2) for some integer q. If we take the end points of our adiabatic path to be these two symmetrypreserving systems, then the topological response in this case turns out to be the difference 7 The integer contribution 2π(q − n) to the bulk polarization can be removed by stacking the system up with other trivial insulators, and hence the topological response is due to the θ = π contribution. Indeed, for systems with inversion symmetry, θ is quantized to be a multiple of π, and there is a bulk Z 2 topological invariant that distinguishes even (trivial) and odd (non-trivial) multiples of π 7,22,23 .
We also see one more interesting feature. Since the response equations depend on spatial gradients of θ, there is some consequence of gapped spatial evolution (as opposed to adiabatic time evolution), i.e., a gapped interface or boundary where θ varies with position. The response equations dictate the quantum numbers bound to regions of the system where θ is changing in space and the system remains gapped. For the present example, the gapped interpolation between the interior and exterior of a material where θ changes determines the amount of bound electric charge at the interface. This is equivalent to the usual boundary charge theorem for polarized insulators and the result of Su, Schrieffer and Heeger, (and Jackiw and Rebbi earlier 25 ) that solitons in their (spinless) model bind localized half-charges.
To summarize, we have seen that closed adiabatic paths can lead to a pumping process, while open, symmetry-breaking paths with symmetric end points can serve to determine topological (electromagnetic) responses. Finally, these considerations also determine bound charges/states in regions where the adiabatic parameters vary in space.

III. THOULESS PUMPS WITH LINE SINGULARITIES
We now consider an enlarged parameter space. We will be as general as possible to illustrate how the reasoning works in systems as arbitrary as possible. The subsequent sections will discuss specific examples.
So far the effective pumping parameter space has been a two-dimensional plane parametrized by polar coordinates m and θ with the origin m = 0 removed. Every other point besides the origin represents a gapped Hamiltonian. Additionally, in this parameter space, only the x-axis, with θ = 0 or π, represents Hamiltonians obeying charge conjugation and inversion symmetries. We now enrich the parameter space by doubling the number of fermionic degrees of freedom. In 1D insulators, a simple model for such a system consists of a spinfull Su-Schrieffer-Heeger (SSH) chain, built from two copies (four internal degrees of freedom) of the previously discussed model. In 3D, a simple model consists of a time-reversal invariant 3D gapped Dirac model with an additional doubling (eight internal degrees of freedom). Hence, it represents two copies of a usual 3D topological insulator. Both doubled systems preserve the symmetries of the originalprior to doubling-ones.
Since these systems can be described by Dirac Hamiltonians, let us consider then for concreteness a general Dirac-type insulator described by the Bloch Hamiltonian It depends on momenta k and three parameters which we call m, m 5 and ∆; this exhausts the scenarios that we consider in this article. The adiabatic parameters control the gap in the system and dictate the possible classes of adiabatic paths.
To be explicit, we take m to represent the usual "topological insulator gap", or TI mass. It conventionally preserves any requisite symmetries, and will act diagonally in the two fermion subspaces resulting from the doubling. In the presence of the protecting discrete symmetries, its value (or sign in a continuum model) effectively defines the topological regime of the model. We can use the TI mass to define the reference points, namely (m, m 5 , ∆) = (m, 0, 0) and (−m, 0, 0) , where by convention m > 0. For our regularization convention, the former will be a symmetry-preserving trivial insulator. The latter would usually represent a symmetry preserving topological insulator if the system were not doubled; the doubling trivializes any Z 2 topological invariants. We note that a continuous evolution between the reference points in a straight line in parameter space, flipping the sign of m, has a vanishing gap at m = 0.
The mass term m 5 is chosen to break all of the symmetries that would protect a topologically non-trivial (short-range entangled) phase. Hence, one can always find a gapped interpolation between the two reference points even when constrained to the 2D parameter space (m, m 5 ). To be concrete, we choose m 5 to act identically within each copy of the doubled system as the usual axial mass. The symmetry it breaks is dimension dependent, e.g., in 1D it breaks charge-conjugation and inversion, and in 3D it breaks time-reversal and inversion. This is the usual symmetry-breaking mass term we considered before, when defining the charge-pumping operation in the 1D spinless SSH model.
Additionally, since we have two copies of a Z 2 TIthus a topologically trivial system-the remaining mass ∆ can be chosen such that it may or may not preserve the protecting symmetries; either way one can still adiabatically connect the two reference points, now in the generic 3D parameter space (m, m 5 , ∆). We will choose ∆ such that it couples the two copies. Doubling a Dirac model always introduces three distinct possibilities for such a mass term. For our analysis it will be always chosen as to break the protecting symmetry, and possibly other discrete anti-unitary symmetries not required for topological stability. For example, in our 1D model it is chosen to break both charge-conjugation, and it breaks time-reversal symmetry, the latter not required for topological stability.
As constructed, both of these masses are compatible with m, i.e., their presence enhances the insulating gap (the corresponding mass matrices anti-commute with the TI mass). They, however, compete among themselves. Physically this means that if both are non-vanishing, there exist certain regions of the parameter space such that their contributions to the gap cancel each other. Mathematically, it means that the mass matrices commute with each other.
For a non-interacting, translation invariant Dirac system, this can be seen explicitly by looking at the spectrum at k = 0, with both masses included: Clearly, when m 5 = ∆, the gap reduces to the TI mass contribution alone, which always vanishes at some point along the interpolation path. In practice, the competing mass terms introduce new singular lines into the 3D parameter space in the m = 0 plane when m 5 = ±∆. These lines will play an important role below, but let us briefly mention what the result would be if they were not there. If both masses were compatible with each other, instead of competing, then the spectrum would be This is only singular at the origin (0, 0, 0). Since any closed path in this 3D space is contractible, there cannot be any non-trivial, quantized pumping processes (assuming that there is no enforced symmetry that forces the path to lie fixed in a plane). Since each path with end points at our two reference points can be continuously deformed to all others with the same end points, we would also expect that every path would determine the same difference in the response between the two reference states. Thus, if all three mass terms were compatible, the properties determined by closed paths would be trivial and the response would not depend on the choice of an open path.
Let us return to the case of interest with competing mass terms. As mentioned in the previous section, one is able to determine the difference in the topological electromagnetic response between the two phases of the single 1D SSH chain by comparing the effective actions at the two reference points. The situation is more interesting with the higher-dimensional parameter space. In general, to adiabatically extract the electromagnetic response difference between the reference  states, we must take paths between them which penetrate into the discrete-symmetry-breaking regions of parameter space. In order to unambiguously define charge transfer these paths must also preserve the corresponding continuous symmetries responsible for conserving the relevant charges along the pumping path.
In Fig. 1 we illustrate some relevant scenarios in our higher dimensional space. They consist of the horizontal and vertical trajectories I and II, in the (m, m 5 ) and (m, ∆) planes, respectively. Trajectories III and IV correspond to tilting the previous trajectories I and II towards the singular line m 5 = ∆. The standard procedure to obtain the response in 1D and 3D is a generalized Goldstone-Wilczek calculation 24 . One introduces an interpolation parameter θ, such that (m, m 5 , ∆) (θ) ≡ (m (θ) , m 5 (θ) , ∆ (θ)). Correspondingly, this gives rise to a Hamiltonian computed at an arbitrary point of a given trajectory in the parameter space One's goal then becomes to compute the change in the effective action for the gauge fields coupled to the conserved charges of the original Hamiltonian, induced by changes in the Hamiltonian parametrized by θ. From such effective actions, one can distil the distinct pumping phenomena, as described in Sec.II. Regarding such effective actions, it is important to consider parametrized Hamiltonians which does not pertain to a single point in the parameter space uniformly through space and time. In other words, it may be that the parametrization value may be dependent on the space-time position, θ = θ(x, t). If one defines the reference states as (m, m 5 , ∆) (0) = (m, 0, 0) and (m, m 5 , ∆) (π) = (−m, 0, 0), the Hamiltonian in distinct positions of space-time may pertain to distinct reference states. Such possibility is crucial for us, in order to be able to describe domain walls between the reference Hamiltonians with opposite TI mass. Such domain walls, as discussed, can bind localized modes carrying specific quantum numbers, which one may extract also from the effective gauge theory.
The possibility of space-time fluctuations in θ introduces difficulties in our derivations. Domain walls are represented by strongly space-time varying functions θ. In order to control calculations of the pumping and charge-response effects in our effective theories, even in the presence of such defects 5 , we consider building up the changes in our Hamiltonian adiabatically in infinitesimal steps, all the way up to the profile θ(x, t) desired. To realize this, we introduce a path parameter α. Considering a small variation α → α + δα, with α(x, t) varying slowly in spacetime, we can find the change in the effective action of the conserved gauge fields to first order in δα. We then combine each infinitesimal build up to form a complete functional dependence, θ(x, t). A more detailed account of this procedure is presented in the Appendix.
We have shown examples of different classes of adiabatic paths and the consequences of traversing the various classes of paths will be analyzed in detail for the chosen examples. Let us, however, make a few more general comments first. As discussed, the two reference points respect the imposed discrete symmetries and, naively, one would expect that the difference in the response behavior between them should only depend on the initial and final Hamiltonians. Yet, we will see that there is a layer of subtlety as to how the response is calculated when there are competing mass terms. Although, due to doubling, both reference points correspond to trivial phases, with respect to the 10-fold classification table, the information obtained by comparing states of opposite TI mass signs aids in the determination of the quantum numbers pumped in closed paths, as well as the degrees of freedom trapped in solitonic defects of the parametrization variable. These results are then relevant for both trivial and topological phases.
Finally, the interpolation between reference Hamiltonians can be extended to adiabatic cycles (relevant both in trivial and topological systems since pumping processes do not depend on the starting Hamiltonian). Some generically different classes of adiabatic cycles arise, as described in Fig. 2. In the doubled 1D case, for example, the open paths I and II of Fig. 1 constitute transfers of spin-1/2 and a full unit (due to doubling) of electronic charge respectively. To constitute a true pumping (with matching initial and final ground state subspaces), one may stitch these paths together, resulting in the loop II a of Fig. 2, in which the quantum numbers of a "full electron" are pumped. The closed trajectories I a and I b then correspond to pumping of pairs of unit-charges and spin-1/2 moments, respectively. As discussed, trajectories II a,b combine the separate fractional pumping, with respect to the microscopic contents of a unit cell, of a single charge and a unit of spin-1/2, thus constituting the "full electron" pump when traversed.
In what follows we will discuss in detail the above results, using specific examples in (1+1)D and in (3+1)D.

A. Generalities
To illustrate a generalized Thouless pump in (1+1)D we will use two copies of the SSH model by adding a second copy with equal Hamiltonian. For definiteness, we will consider the two copies to represent fermions with up and down flavors of spin. The Hamiltonian in the continuum limit can be written as follows (14) in the basis Ψ T = ψ A↑ ψ B↑ ψ A↓ ψ B↓ , with ψ Aσ , ψ Bσ fermionic annihilation operators of spin σ. In Eq. (14), A x is the electromagnetic vector potential, m is the mass generated by a staggered hopping between the lattice sites ("TI mass"), and m 5 is the mass generated by a staggered onsite energy. The artificial doubling allows the introduction of three extra mass terms ∆ = ∆ 1 , ∆ 2 , ∆ 3 (i.e., there exist three new 4 × 4 matrices which anti-commute with the kinetic part of the Hamiltonian) which couple the different spin subspaces and split their degeneracy. Table I. Symmetry properties of the masses for the (1+1)D doubled SSH model under particle-hole (C) and time-reversal (T ) transformations. The corresponding symmetry is broken/preserved for ×/ marks, respectively We choose the Dirac matrices: where σ i and τ i are Pauli matrices in the spin and sublattice (or orbital) Hilbert spaces, and Kronecker products are implicit. The matrices σ 0 and τ 0 are 2 × 2 identity matrices. In the absence of Λ 5 and ∆ masses, the model displays time-reversal and particle-hole symmetries, with the operators T = iσ y τ 0 K and C = τ z K, where K is the complex conjugation operator. The remaining masses break these discrete symmetries in a pattern described in Table I. The commutation relations, with a = 0, 1; i, j = x, y, z, imply that, although all masses are compatible with the TI mass m, the mass m 5 and the set ∆ actually compete. In the absence of m, this can lead to gapless spectra when |m 5 | = |∆|.
Let us consider the effective electromagnetic action and quantum pumping processes in this system. First, let us look at the situation with ∆ = 0, m 5 = 0. This corresponds to just a pair of disconnected linear chains, and each chain is known to have topologically distinct phases for m > 0 and m < 0. To find the topological response from each independent chain, one interpolates between these phases using H (α), with α an angle varying slowly in space-time as discussed in Sec. III (the complete, doubled chain, system is trivial overall and will give a trivial total response). In order to achieve an adiabatic interpolation, the gap must not close as α is varied. Hence, we must avoid the point m = m 5 = ∆ = 0 and the lines m = 0, m 5 = ±∆. From this parameterization we can see that this model contains all the features discussed in Sec. III.
As an example of a class of adiabatic paths, consider the following interpolation with α changing from 0 to an arbitrary angle θ. This corresponds to the "pure" axial-mass path, labeled path I in Fig. 1. The fermionic path integration for the change in effective action due to an infinitesimal α, described in detail in Appendix A, is nothing but a Goldstone-Wilczek-like calculation and gives rise to an axion term. Building up the contributions of infinitesimal changes in α from 0 to an arbitrary finite angle θ, one finds where Θ 1D is the standard 1D "θ-term" This action dictates both the adiabatic pumping process (by considering full evolutions of θ by 2π), and quantum numbers bound to solitons, like domain-walls in the TI mass. In particular, the latter case can be described by is the Heaviside step function. In practice, this implies that each wire contributes charge ±e/2 at the ends of the system, yielding in total an integer charge bound to the edge (and likely unobservable in realistic systems because of the charge polarization ambiguity in lattice models 26 ). One would naively expect the response difference between reference Hamiltonians at θ(x, t) = 0 and θ(x, t) = π to be a property associated only with the initial and final Hamiltonians H (0) and H (π). The path taken in the parameter space, however, can be easily seen to be relevant in our model. For example, considering the path II in Fig. 1 (for ∆ = ∆ zẑ , say), one lets ∆ z = m and defines, Following the same Goldstone-Wilczek type of calculation as before, one finds for the effective action which is clearly different from Eq. 18. If the initial and final Hamiltonians H (0) and H (π) are the same, we would expect that any continuous deformation of the path in parameter space should result in the same response difference for the two cases. The discrepancy between (18) and (21) can only arise from the fact that the masses m 5 and ∆ z compete. This is so as one cannot deform path I to II without passing through a gapless point/line. Further support for this idea can be found by calculating the response difference on the tilted paths III and IV in Fig. 1. One can do this by considering, for example, +Γ 0 m cos α + Λ 5 m sin α + Λ z ∆ z sin α and computing corrections in the effective action perturbatively in ∆ z /m. As long as this dimensionless parameter is small, we do not cross the gapless line. The procedure is described in detail in Appendix A. An analogous calculation can be done for path IV, perturbatively in m 5 /m. What one finds is that the effective actions do not change from their respective original behaviors, i.e., Table 3 summarizes the above results. They suggest that, indeed, the change in the behavior between paths I and II is due to the discontinuous jumps at singular lines in parameter space, as opposed to a smooth continuous evolution. There is, however, a subtlety in this discussion of the electromagnetic response. As far as the topological response properties of lattice fermionic systems go, S I ef f is essentially equivalent to S II ef f since they both represent the response properties of trivial insulators. To put this another way, the integer charge polarization from Eq. 18 can be removed by a gauge transformation on the Bloch wavefunctions. Closed adiabatic evolutions of the Hamiltonian, however, still lead to pumping of quasi-particles. Closing a circular path along trajectory I pumps twice an electric charge in the present scenario.
Remarkably, these considerations show that a closed path evolution of the Hamiltonian around the gapless line, i.e., going forward on path I and reverse on path II will not bring the system back to its original state, which indicates a ground state degeneracy, and a type of Berry phase holonomy that is being gathered in this process. As a further consequence, we expect some other 'hidden' charge pumping during this process. In fact, the spinfull SSH Hamiltonian allows for a set of SU(2) conserved charges, i.e., just the spin components themselves. The different ∆ masses break this SU(2) spin invariance, but an arbitrary path in parameter space through ∆ still preserves a spin-U(1) symmetry for rotations along a given axis. Fixing a single conserved component at, say, S z and gauging this symmetry by introducing a gauge field A Sz ν , path II now induces a finite response. One finds, following the regular Goldstone-Wilczek calculation, where θ I,II corresponds to adiabatic evolution of the families of Hamiltonians along the I or II trajectories. Also, although we chose spin conservation along the zdirection, an arbitrary direction could have been chosen, with an appropriately parametrized corresponding rotation. These results suggest that the competition of mass terms is related to the known physics of spin-charge separation in spinfull SSH chains. To further explore these phenomena, we now proceed to consider a microscopic picture of this problem, which is convenient to study the physics of edge bound modes. Figure 3. Summary of the electromagnetic θ-term responses calculated in each path of Fig. 1 for the spin-doubled 1D system.

B. Microscopic picture
We can shed some light on the above findings by considering an exactly solvable lattice model of the continuum model above, still parameterized by m, m 5 and ∆. The lattice (doubled) SSH model is where σ = ± for up/down spins, and η is the deviation from the (normalized to 1) initial hopping amplitude. We have only included the ∆ z spin-dependent mass term, though we will mention the other spin mixing terms later.
The TI mass parameter m is fixed by η.
It is enough to study the limit η = −1, beginning with m 5 = ∆ z = 0. At this point in parameter space, the correlation length vanishes, and the eigenstates are formed by which leave the exact zero-modes a 1σ and b N σ at the ends of the open chain. The presence of the spin degrees of freedom from the doubling enforces additional degeneracies in the ground state of the open chain. In this flat band limit, we can write down the boundary theory with m 5 and ∆ z as perturbations, which simply reads in the basis |1 a ↑ |1 a ↓ |N b ↑ |N b ↓ , where 1 a , N b represent the sites on the left/right end.  ∆z, one sees that the system is polarized (b) m5 ∆z, one now sees the spin polarization; Blue and green circles correspond to filled states in the bulk and edges of the chains respectively. Effects of finite values of m5, ∆z on the site densities are represented by the relative size of the circles. In the thermodynamic limit with m5, ∆ → 0, the staggered density pattern disappear but the filled ground states remain the same.
From (28), the scenarios of bound states implied by the axion action (25) can be studied easily, and are described pictorially in Fig. 4. At half filling, i.e., when the bulk periodic chain is insulating, two out of the four boundary modes will be occupied. When |m 5 | > |∆ z |, the lowest energy states reside on a single edge. This leads to a ground-state occupation of the states on the given edge, and hence a non-vanishing charge polarization of the system, which then breaks C and inversion symmetries. On the other hand, for |m 5 | < |∆ z | the ground-state occupation has one state occupied on each edge, and for our choice of a ∆ z mass term these states have opposite spin projection in z direction. This system has no charge polarization, and hence no boundary charge, but it does have overall dangling spins at the edges.
The finite polarization state may be understood as a spontaneously symmetry broken state. Concretely, one considers finite m 5 and ∆ and fills the two lowest energy boundary modes. In the thermodynamic limit with an infinite chain, the states on the boundaries are isolated and cannot tunnel between each end. If we take limits ∆ z → 0 first, followed by m 5 → 0, we are led to a spontaneously symmetry broken state with nonvanishing charge polarization due to the gapped bulk. Inverting the order of the limits, however, leads to a state with dangling spins at the edges. This state presents no spontaneously broken symmetries; the only possible broken symmetry would be time reversal, however, as the two opposite spin states are degenerate at each edge, any superposition is possible for the ground state, and timereversal symmetry is not broken. Under a spin U(1) symmetry conservation (in the z-direction, as we considered before), all that can be said is that the domain wall in the TI mass carries a single unit of spin 1/2.
Regarding the other spin mixing terms, if we consider non-zero values of ∆ x , ∆ y , as well as ∆ z , we may find a different direction for the end state spin polarization. Nevertheless, if m 5 → 0 and then the ∆ i → 0 limit is taken, each end electron forms a spin-1/2 degeneracy with full SU (2) invariance. This also matches the result of the effective action with gauged SU(2) symmetry.
This picture is the hallmark of spin-charge separation for SSH models with spin. [27][28][29] Domain wall excitations in this system can be chargons or holons, with charge (±e), or spinons, with spin-1/2. These results are not limited to the flat band limit, and one may verify them analytically by considering Eq. (26) with η = 1. One can simulate the edge by a domain wall with a kink in the TI mass. The other edge is simulated by an anti-kink. Introduction and projection of m 5 and ∆ z into these states recovers the discussed results. Even though the charged kinks carry integer charge they are still fractionalized objects compared to the local degrees of freedom in the unit cell since the fundamental excitation is an electron with both charge and spin. For example, adding a 'full' electron to a positively charged kink takes away its charge but converts it to a spin-1/2 kink. Thus, these fundamental kink defects will contain spin or charge, but not both. This type of spin-charge separation is also seen in flux defects in the quantum spin Hall state 30,31 .
The microscopic picture of spin-charge separation is linked to the two planar (non-tilted) adiabatic paths by our action Eq. (25) (in the limit of spin-z conservation, for simplicity). Suppose now we take a closed adiabatic path which encloses a singular line. We start in a symmetry preserving state with m 5 = ∆ z = 0 and then perform a π-rotation in trajectory I (m − m 5 -plane) and then a −π-rotation in trajectory II (m − ∆ z plane) that brings us back to our starting point. We can ask if this system is different than, say, taking a full 2π-rotation along trajectory I (m − m 5 -plane.) If we consider these adiabatic paths as spatial evolutions/interfaces, then from the analysis of the bound quantum numbers, in the first case we would have both a bound chargon and a spinon, i.e., an additional full electron, and in the latter case we would have either two chargons or two holons. Hence the different classes of adiabatic paths give rise to physically distinct pumping processes, which are captured by Eq. (25) when considering domain-walls in the corresponding θ I,II angles.

C. Adiabatic Pumping and the Mañes-Bardeen formula
We can also consider the adiabatic evolutions of the gapped 1D Hamiltonians in a more systematic way. Effective actions carrying topological content can be deduced with minimal labor by deploying an analysis pioneered by Witten 13 and further developed in Refs. 14 and 15. Given a system of fully gapped fermions, their currents are defined by variations with respect to external gauge fields via the effective action for these fields. Such effective actions may carry topological terms, which are defined by corresponding anomalies; the results of Refs. 13-15 indicate how to construct these topological terms systematically, for any type of fermionic mass fields.
To apply this method, let us re-analyze the doubled SSH model Hamiltonian (14). Noticing that iΓ 0 Λ 5 =Γ defines a chiral symmetry, one may rewrite the fermion operators as Ψ → χ R,L = [(1 ±Γ)/2]Ψ. We will start from the most symmetric system and introduce the analysis of the pumping paths by a symmetry breaking procedure. Thus, the starting point is the action of the doubled massless relativistic fermions in 1+1D, Here,χ i = χ † i γ 0 and the Lagrangian Dirac matrices are Also, the chiral matrix reads, explicitly,Γ = σ 0 τ y . In general we can introduce a pair of non-Abelian gauge fields R, L, which gauge all the symmetries of this model. Let us analyze this problem in depth. The fields display a vector U(1) symmetry χ L = e iθ χ L and χ R = e iθ χ R and an axial U(1) symmetry, χ L = e iΓθ χ L = e −iθ χ L and χ R = e iΓθ χ R = e iθ χ R . Furthermore, they display an SU (2) + × SU (2) − symmetry. To see this, consider the following matrices In our particular representation, Notice that λ a =Γt a and that the t a = S a matrices are nothing but the spin operators from (24). The commutation relations obeyed by these matrices are [t a ,Γ] = [λ a ,Γ] = 0, as well as The SU (2) + ×SU (2) − symmetry is manifestly generated as Now we will begin introducing mass terms, and hence start breaking the symmetries. First, we consider adding the TI mass, which breaks the axial symmetries locally and globally. In fact, this mass term also breaks the SU (2) + × SU (2) − symmetry to its diagonal SU (2) subgroup, g + = g − = g and g + g − = e itaφ a . We thus must constrain the non-Abelian gauge fields R = L ≡ A and the action reads The terms on the second line represent the TI mass. The remaining symmetry of the problem is now SU (2)×U (1), and we gauge it as where A S and A e are the spin and electromagnetic gauge fields. It is our desire to consider transformations of the original Hamiltonian which send m → −m. For this to be accomplished, we will have to consider transformation paths which break the discrete symmetries of the Hamiltonian, while conserving charges for transport. This can be achieved, in the present context, by considering where generically U = e i(θ I +θ II ·t) .
We are now set up to consider the pumping processes. The two limiting cases of interest are: (i) path I rotations parameterized as θ I : 0 → π, |θ II | = 0, and (ii) path II rotations parameterized as θ I : 0, θ II = (0 → π)ẑ. Notice that path II breaks the SU(2) symmetry, but not fully. There is always a projection of the charges along the path which is conserved, which we fix to be t 3 = S z for concreteness. We can also consider a general interpolation between these two limiting cases, as whereα is a unit vector pointing in the SU(2) t direction corresponding to a choice of a conserved U(1) subgroup. Hence, the interpolation between positive and negative TI masses m is fixed by α : 0 → π, while the interpolation between paths I and II is fixed by φ : 0 → π/2. The gauge fields along these paths are reduced to just These are deduced from the conserved U(1) charges, spin and electromagnetic, and are the fields which couple to the conserved charges whose pumping we probe.
At this point, we are ready to introduce the gauged Wess-Zumino action of Mañes-Bardeen. This action captures the topological part of the response of gapped fermions to variations in its mass-generating fields (as well as external gauge fields). Such a topological contribution is fully determined by gauge invariance and anomaly inflows, which are related to the computation of triangle graphs as those used in our previous diagrammatic analysis (for a recent in-depth analysis, see Ref. 32 and references therein). This connects Mañes-Bardeen's Wess-Zumino action directly to our problem at hand; we will soon see that the method has computational advantages over the diagrammatic approach when considering general interpolation paths between the reference points in the presence of all conserved charges gauge fields. In 1+1D, the gauged Wess-Zumino action of Mañes-Bardeen reads, in differential form notation, where ω 3 is the Chern-Simons form where we have restored the gauge fields L and R for a moment, M 2 is the 1+1D manifold, M 3 is the extended manifold for the Wess-Zumino-Witten term, and the normalization constant C can be fixed, for example, by considering the usual (anomalous) chiral rotation, whose result is known, and then comparing the pre-factors. When considering our gapping terms we must fix L = R = A, and this action simplifies to Since the only space-time dependent function in U is α, see Eq. (43), one can write dU = U π/2 dα, where U π/2 ≡ U (φ, π/2 + α). Also, the result simplifies since the remaining gauge field components commute among themselves, and with U . Finally, the usual chiral anomaly in 1+1D for a doubled Dirac system implies Ctr [1] = 1 2π , where tr[1] = 4 in the present matrix dimensionality. Building up infinitesimal changes in the Wess-Zumino action slowly by choosing α → α + δα with dα ∼ 0 at each step, and taking contributions to first order in δα only, we allow α to evolve to a full profile θ(x, t) and the full Wess-Zumino action is given by where the charge, θ e , and spin, θ s , 'axion' fields read To show that we obtain the expected results we plot in Fig. 5 the values corresponding to full evolutions m → −m by using a fully uniform θ(x, t) = π. We see that as a function of φ, the pre-factors jump exactly by π when φ crosses the singular line value at π/4. For path I, one sees the build up of electric charge; for path II, spin is pumped, as already discussed.
The arguments of this section generalize for the 3+1D case below, but, first, we will consider a final analysis of the 1+1D system via bosonization.

D. Adiabatic Pumping Using Bosonization
Since we are considering 1+1D systems, this allows for a non-perturbative verification of the previous results using a bosonization picture of the continuum model above 11 . Considering (14), we may proceed with the standard bosonization recipe, πφ Lσ for bosonic fields φ R,L . Generically, Klein factors would appear when analyzing the mass terms Λ x or Λ y , but to simplify our discussion we will only keep Λ z in this section, as to avoid discussing the non-Abelian bosonization of the full SU(2) symmetric system. We can define the charge and spin boson fields: With a normalization of the boson fields by √ 2π, the Hamiltonian (14), with ∆ x,y = 0, is mapped into the Lagrangian where A µ = A 0 , A 1 is the electromagnetic vector potential, A Sz µ is the corresponding spin-gauge field, and a is a short-distance cut-off scale.
One of the advantages of the bosonized picture is that we can see explicitly how the electromagnetic response of the system is associated only with a charge-like bosonic field φ c , and that there is no coupling between φ s and A µ . Hence, we would not expect that shifts of φ s would generate any electromagnetic response, something that is not as obvious in the perturbative field theory calculations of the Dirac fermions above. A similar analysis holds or the spin-gauge field response.
Since the ground state configurations in the bosonized picture are simple, this allows us to study the properties of the full parameter space, even where perturbative calculations of the Dirac fermions would break down close to the gapless singular lines in the adiabatic parameter space. At low-energies, a semi-classical study of the problem is enough for our purposes. First, the kinetic terms guarantee that the ground states should be homogeneous constant fields, φ c = const. and φ s = const. . Now, let us first consider m 5 = ∆ z = 0. One reads from the potential that there are two degenerate minima for each sign of m, namely if m < 0. (54) As one might expect, m = 0 leaves an ill-defined (gapless) ground state. We note that to adiabatically connect the m > 0 pair to the m < 0 pair we need to consider a path with finite m 5 or ∆ z in the Hamiltonian parameter space.
A possible question which may arise at this point is: why does the bosonized version of the closed/periodic (1+1)D problem have two ground states (for a given mass sign) while the original free fermion problem only appears to have one? The answer is that in our original fermionic problem we did not discuss the possibilities of different types of periodic boundary conditions for the fermions (i.e., different spin structures). Indeed, distinct pinned bosonic ground states can be accessed by piercing the closed fermionic chain by π−fluxes of the charge/spin gauge fields. Such fluxes flip the boundary conditions from periodic to anti-periodic, hence flipping the parity of the fermionic ground states. For our system there are two parity operators, giving four possibilities for the ground states. These are the states manifest in the present bosonic language. Now let us consider adiabatic evolution. We can nicely parameterize our paths of interest by the potential With α evolving from 0 to π, we may consider the two different sets of π-rotations from Section IV A. If ∆ z = 0, this evolution corresponds to the usual chiral rotation which, in the bosonization formalism, is translated to Then, from (51), and we see that this transformation connects the ground states (0, π) → (π, π) or the ground states (π, 0) → (0, 0). In either case, the change in the Lagrangian is simply |∆L| = e µν ∂ ν A µ which, as expected, corresponds to twice the 1D θ-term, i.e., the same result found in Eq. 18. Hence, rotating φ c by 2π recovers the adiabatic charge pumping mechanism that was discussed above in a fermion language. Analogously, if m 5 = 0 and one follows the evolution in the m − ∆ z -plane, one finds The ground states are now connected as (0, 0) → (π, π) or (π, 0) → (0, π), and there is no change in φ c . In fact, the additionally imposed U (1) spin gauge field coupled to the conserved spin component generates a θ-term for the spin gauge field, and shifting φ s by 2π acts as an adiabatic spin pump. Now, if one considers finite values of both m 5 and ∆ z , as we evolve α, the minima in (φ c , φ s )-space will start to trace more complex trajectories, straying away from the horizontal and vertical lines of the two special cases above. However, as long as |m 5 | = |∆ z | then the starting and ending points match the results for max(|m 5 |, |∆ z |). Also, despite the degeneracy of the ground state for a fixed sign of m, mixing between the ground states never happens. On the other hand, if |m 5 | = |∆ z | then during the path, when α = π/2, we find a valley of degenerate minima connecting the original two minima of the  Figure 7. Evolution of the ground states for m = 1. Symmetry allows us to consider only values of 0 ≤ φc, φs ≤ π. The dark blue circles correspond to the minima of V (φc, φs, α) for m > 0 while the white circles correspond to the minima for m < 0. The red dashed lines show the trajectories of the minima, red dots, as α is swept from 0 to π (a) m5 = 1, ∆z = 0.25, during the evolution the minima move mostly in the horizontals but connects the two pairs of ground states uniquely. (b) ∆z = 1, m5 = 0.25, In this case, the minima move mostly in vertical but again the ground states are always well defined. system, thus allowing for a non-adiabatic change in the ground state, and making the evolution ill defined, as expected. These considerations are summarized in Figs. 6 and 7.
Interestingly, we can also consider the hybrid path discussed earlier where we first rotate in the m − m 5 -plane by π and then in the m − ∆ z -plane by −π to return to the starting point. In this case we see that during this process we have the transformation Hence, for a fixed sign of m, this process keeps us in the same degenerate ground state subspace but switches among the two ground states. Thus, we can see that the singular line in our parameter space acts as a kind of non-Abelian Berry flux that acts to switch the ground state. Encircling the singular line twice will return us to the original ground state.

V. 3+1D DOUBLED TI
So far we have studied generalized adiabatic pumping processes in doubled 1D systems. For this, we computed electromagnetic-and spin-gauge effective responses for 1D insulating systems via adiabatic transformations of the Hamiltonian. From the doubled model we were able to develop inequivalent classes of paths in parameter space. Depending on the path chosen, the effective action implied that the adiabatic transformation led to a charge-pump, a spin-pump, or a "full-electron" pump.
The general results and procedures in 1D and 3D are very similar. From the point of view of electromagnetic responses, adiabatic pumping in 3D is signaled by the Chern-Simons axion coupling where A is the electromagnetic vector potential, and θ parametrizes the adiabatic deformation of the Hamiltonian. The corresponding currents derived from this effective action are To generate a response one needs both a space or time variation of the adiabatic parameter, and an electric or magnetic field, respectively. The physics of this axion coupling is exactly the topological magneto-electric effect that exists in 3D time-reversal invariant topological insulators 7 . In this context, the corresponding pumping process obtained by the evolution of the θ parameter has been dubbed a "Chern-Simons axion adiabatic pumping" 33 . An simple way to understand this response can be developed in the presence of a uniform magnetic field. In this case, the response becomes analogous to copies of the response for a 1D Z 2 topological insulator, i.e., the usual Goldstone-Wilczek response. The number of copies is simply equal to the number of flux quanta. We will adopt this idea and use it to explore the different quantum numbers pumped in this context.

A. Generalities
Consider doubling a minimal model for a 3D TI Hamiltonian and including additional mass terms. For example, we study the model now with the basis Table II. Symmetry properties of the masses for the (3+1)D doubled TI model under charge-conjugation (C) and timereversal (T ) transformations. The corresponding symmetry is broken/preserved for ×/ marks, respectively. Our choice of symmetry structure leads to an exactly mirrored behavior between the C and T columns with respect to the 1D case.
In this basis, σ i represents spin, τ i are orbital degrees of freedom, and ρ i are the doubling-related degrees of freedom (which we do not specify explicitly, but may be of orbital content also, and will henceforth be referred to as "doubling-orbital" degrees of freedom). This model has discrete time-reversal and particle-hole symmetries with operators T = iσ y K and C = ρ y τ y σ y K. Table II summarizes the symmetry properties of the extra mass terms with respect to these operations. These additional mass terms will define the pumping parameter space as discussed below. As an aside, if one were to identify H 3D with a BdG doubled Hamiltonian, such that C were a strict particle-hole redundancy, then this model is just a 3D topological insulator with a chiral mass m 5 , and ∆ y and ∆ x are the real and imaginary parts of a proximity-coupled s-wave superconducting order parameter. As the electromagnetic charge U(1) symmetry is broken in this case, the concept of pumping is not as well defined and we leave further comments on this to the conclusion. Another interpretation for this model is that of a topological crystalline insulator. This system has a trivial strong invariant, but will have non-vanishing mirror Chern numbers for the mirror operators given bŷ m x = iΓ x Λ 5 ,m y = iΓ y Λ 5 , andm z = iΓ z Λ 5 . We can see this by restricting the model to the k x = 0, k y = 0 or k z = 0 planes respectively and calculating the mirror Chern number which we find to be C x = 2, C y = 2, or C z = 2 respectively, if we keep our previous choice of regularization where a negative mass represents the topological phase. This implies that on surfaces withm x ,m y orm z symmetry there will be two surface Dirac cones stabilized by mirror symmetry that lie along a mirror line; recall that only an odd number of cones is stabilized if we just have time-reversal symmetry.
When exploring the response and pumping processes, many of the results from Sec. IV A hold in the 3D case with few modifications. The extra mass terms are again compatible with the usual TI mass, but compete with each other. We find that, of example, the evolution through path I computed from the standard triangle anomaly using the m 5 mass gives twice the standard θterm, i.e., Figure 8. Summary of the responses calculated in each path of Fig. 1 for the particle-hole doubled (3+1)D case. where All of the responses can be calculated similarly to our method in 1D, and we explicitly show this in Appendix A. Table 8 summarizes the results obtained for the doubled 3D case using these perturbative methods. In what follows, we give a microscopic interpretation for the (3+1)D case.

B. 3D TI surfaces in a magnetic field
In this section we aim to connect the results found in (3+1)D in Sec. V A with those in Sec. IV for the (1+1)D system, so that we can easily interpret the types of quantum numbers are pumped on each class of paths in our adiabatic parameter space. We will consider this in the context of the doubled 3D TI response (65) in the case of a uniform magnetic field, say, in the z direction, If we defineĒ z =´dxdyE z /A xy to be the average of the electric field in the xy plane with area A xy , we find withθ = eΦ 2π θ = N Φ θ, Φ = B z A xy is the magnetic flux through the plane, and N Φ is the number of flux quanta in units of Φ 0 = 2π/e(= h/e). Because of this form of the effective action we expect that the surfaces of the doubled 3D TI in the presence of a uniform magnetic field must behave analogously to (N Φ copies of) the ends of the doubled SSH model. In particular, the spin and charge pumps considered in Sec. IV should have an analogy in (3+1)D.
To see this we can solve the surface Dirac Hamiltonian explicitly when in the presence of the uniform perpendicular magnetic field. The surface physics of the 3D TI is known for being described by massless Dirac fermions, and the surface modes may be found by considering a TI mass with a domain wall as a simple model for the surface and looking for states localized around z = 0. An ansatz for these localized states takes the form where N normalizes the exponential factor, and φ is a normalized spinor satisfying For the doubled system, the spinors are of the form φ = (φ 1 , φ 2 ) T , where φ 1,2 are four-component spinor wavefunctions of only the x, y coordinates. This reduces the matrix dimension of the problem by two, whose solution allows for the definition of a basis for the localized states and the basis Now we may project the original Hamiltonian (61) in the new basis, to obtain the surface Hamiltoniañ where ⊥ fixes only components in the xy plane. The projected matrices areΓ whereσ and ρ are 2×2 matrices in the basis (72). Notice that ρ still represents the doubling degrees of freedom of the original Hilbert space, while, from (71), the states χ 1,2 have well defined spin (up/down respectively). The gauge field, being z-independent, is unaffected by zaveraging Now we may study the Landau level problem of the surface, our aim being to find the low-energy degrees of freedom bound to the surfaces in a magnetic field, and then to consider the effects of masses m 5 , and ∆ perturbatively as we did for the end states of the (1+1)D system. For the surface, the blocks of the doubled degrees of freedom decouple and one has to solve the eigenvalue problem for with i = 1, 2 fixing the two surface copies, and where we have replacedÃ with A since they are equivalent here. The solution to this problem is known and is analogous to studying the Landau level problem in graphene 34 . Diagonalization may be achieved by considering A in the symmetric gauge A = − r×B 2 = B 2 (−y, x, 0) and writing complex coordinates ξ = x + iy. With these one defines the ladder operators acting on each copy of the surface, such that withσ ± ≡ 1 2 (σ x ± iσ y ) and the cyclotron frequency in the present units reads ω c = √ 2B. Squaring the Hamiltonian shows that the number operators O † i O i are good quantum numbers and one can solve the Hamiltonian to find the energies to be E ±N = ±ω c √ N . Of importance for our analysis are just the, lowest energy, zero modes where ξ, ξ * |0 i is the zero-mode wave function in the ith copy. Notice that in our surface of choice, theχ i,0 spinor is totally polarized in the upper χ 1 space, which we see from the previous analysis, carries spin ↑. Hence, both copies of the surface state zero-modes are totally up-spin polarized; they will be down polarized on the other surface. Now we can treat the mass terms perturbatively in the subspace of degenerate zero modes, including each copy. The projection into the zero-mode subspace implies Assuming that the surface we considered was the upper one, we write the low-energy surface Hamiltonian where u stands for the upper surface. To closely compare with our (1+1)D lattice model results we need to also consider the lower surface. The projection for the lower surface can be found in practice by simply inverting the sign of the external magnetic field, such that spin will then point anti-parallel to the surface. We may finally write in the basis |u ↑ , 1 |u ↑ , 2 |l ↓ , 1 |l ↓ , 2 , where u ↑ , l ↓ label the up-spin upper surface and down-spin lower surface.
With this form of the perturbative Hamiltonian we could repeat a discussion analogous to the one in Sec. IV B for the 1D case. In fact, comparing H surf derived in here with H edge from Sec. IV B shows the exact same result, the only difference being the quantum numbers carried by the end modes, and the fact that there is a Landau level degeneracy in 3D equal to the number of magnetic flux quanta penetrating the surface. In the present case, the interplay between the strength of m 5 and ∆ z determines whether we will have spin polarization on the surface, or a polarization in the doubling-orbital degrees of freedom.
Going back to our interest in pumping processes, these surface effects imply that evolutions through path Ib, from Fig.2 should correspond to pumps of the (surface) electromagnetic Chern-Simons coefficient, i.e., an electromagnetic surface Hall conductivity, i.e., one complete cycle pumps a Hall conductivity of |2e 2 /h|. As we will see in the next section, this effect is also accompanied by a pumping in the surface Chern Simons term for the gauge field coupling to the conserved U (1) charge of the doubled-orbital degree of freedom. Interestingly, the corresponding effect for path Ia is that the electromagnetic field may pump a mixed Chern Simons term between the electromagnetic and doubled-orbital gauge fields. In this case, electric (magnetic) fields of one gauge field will drive a Hall current (trap charge) of the other gauge field. This suggests that the bulk action in this case should contain a mixed axion term between gauge fields corresponding to the electromagnetic field and the conserved corresponding charge from the doubled-orbital since this translates to a mixed Chern-Simons term between these gauge fields on the surface. We will see in the next Section that this is realized precisely. Finally, evolutions through paths IIa and IIb should correspond to a combination of pumps between the regular Hall (electromagnetic and doubledorbital) coefficients and the mixed Hall ones, with opposite signs between IIa and IIb.

C. Mañes-Bardeen Anomaly in 3D
In three spatial dimensions, the bosonization analysis is not easily available. However, it is possible to repeat the anomaly analysis derived from the Bardeen form of the gauged Wess-Zumino action 14,15,35 . This will allow for a non-perturbative check on our perturbative response and pumping results, e.g., to confirm there is no axion-term build up in ∆-mass type evolutions of the Hamiltonian.
The analysis of the gauge symmetries, the action, and paths follows exactly as in Sec.IV C, as that is dimension independent. In our representation of the 3D system the SU(2) matrices now read while the axial U(1) one reads Now one may consider interpolation paths which breaks the doubling-orbital SU(2) symmetry generated by t i into a doubling-orbital U(1) Considering a general interpolation between path I and path II by a path which preserves t 3 , one writes where A o3 is the gauge field corresponding to the conserved SU(2) component related to the doubled-orbital degree of freedom.
In this basis, we may perform an analysis, similar to what we carried out in the 1+1D case, and introduce the Bardeen form of the Wess-Zumino action in the presence of the U matrix, which captures the system's response due to anomalies. In 3+1 dimensions the result is much more involved, and reads, for vector gauge fields 35 where M 3 and M 4 represent the (3+1)D TI manifold and extended (4+1)D manifold for Wess-Zumino-Witten term, respectively. Plugging our specific form for the Umatrix one achieves great simplification. The constant C can be determined from the usual chiral anomaly for a doubled TI system via CT r [1] = 1 120π 2 . We finally find S W Z = 2 8π 2ˆθ e ( x, t) (dA e dA e + dA o3 dA o3 ) Notice that the θ e and θ s functions are the same as the ones computed in the 1D scenario in IV C, except now the full evolution profile θ( x, t) may now fluctuate in (3+1)D space, not only in (1+1)D (a note on notation, do not confuse the functions θ e and θ s with the profile θ of evolution of the parameter α. Compare with equations (49) and (50)). As a consequence, path I binds axion terms for both electromagnetic and doubling-orbital U(1) gauge fields, while path II binds mixed axion terms between the electromagnetic and doubling-orbital gauge fields.. As a check, the integrand of θ e in the limiting cases of φ = 0, π/2, and π/2 − (for infinitesimal ) can be analytically shown, to first order in , to reduce exactly the results found from the perturbative diagrammatic approach of the calculation of the electromagnetic part of the action. More generally, for the evolution to a nonuniform θ( x, t) one may use these expressions to study surface effects. For example, if we evolve the system to a profile θ( x, t) equal to π for z > 0, however locked to 0 in for z < 0, then from Fig. 5, one may write where Θ is the Heaviside step-function. Plugging into the axion effective actions and integrating by parts will turn the spatial-dependent step functions into Dirac deltas which act to bind remaining part of the effective action to the domain-wall/surface (i.e., z = 0). The effective description on the surface for evolution on path I of Fig.1 is electromagnetic and doubled-orbital Chern-Simons terms, while for path II is mixed Chern-Simons between them, as expected from our previous discussions.

VI. CONCLUSIONS
Many decades ago, insulating systems were thought to be uninteresting from the point of view of transport physics. However, Thouless pumps are one example that demonstrate exotic charge transport phenomena in insulators. By introducing multiple copies of canonical minimal models for Thouless pumps in (1+1)D and (3+1)D, we extended the parameter space for adiabatic Hamiltonian modifications into a larger, multiply-connected space. This is possible because the extra degrees of freedom allow for the addition of different gapping terms to the Hamiltonian which can compete with one another depending on their symmetry properties. Regions of parameter space where the competing terms compensate each other completely define gapless singular regions that must be avoided during pumping processes.
We found that, remarkably, the effective electromagnetic action built in an adiabatic evolution in an open path in such a parameter space depends on the class of path chosen and its corresponding broken symmetry, instead of simply the starting and ending reference Hamiltonians. Full cycle evolutions back to the starting point in parameter space may then define quantum pumps of different quantum numbers and depend on the class of closed loop. In the (1+1)D context with an added spin U(1) symmetry, we found three different pumping processes (i) a charge pump, (ii) a spin pump, and (iii) for a path encircling a singular gapless line, the pumping of a full electron (charge plus spin). From a microscopic point of view, the adiabatic evolutions of the Hamiltonian lead to different, degenerate, ground states introduced by the degrees of freedom from doubling.
The (3+1)D case led to similar results as in (1+1)D, with some reinterpretation of the meaning of the degrees of freedom of the doubled system. In particular, the possibility of realizing the pumps in a mirror-symmetric crystalline topological insulator context implies that the SU(2) (valley) gauge field discussed above may be generated by strain fields. The considerations in this 3D scenario are reminiscent of, and enforce, the concept of Chern-Simons axion pumping described in Ref. 33. We also note that, when viewed from the open path/effective response theory point of view, our Hamiltonian transformations present a realization of the T and S transformations for Sp(2N, Z) actions proposed in Ref. 36 for general (3+1)D U (1) N Abelian gauge theories (with N = 2 here). These generalize electric-magnetic duality on the space of conformally invariant boundary conditions for a free U (1) N -flavor Abelian (3+1)D gauge theory. The Sp(2N, Z) has three generators, of "T-,S-, and GL-"types which should act at the (3+1)D system boundary. While rotations of type of path-I correspond to T-type transformations, which simply adds to the conserved currents the Hodge dual of the field strength without changing the underlying theory (i.e., it transforms the action by adding the corresponding surface Chern-Simons terms), S-type transformations can be realized by rotations of type-II in our formalism. This adds to the surface theory, as we have seen, mixed Chern-Simons terms, which in fact make the background gauge fields actually dynamical. Searching then for actual physical realizations of such systems as presented here, with surfaces gapped spontaneously or explicitly with order parameters in the doubling degrees of freedom may be an interesting route to novel topological phenomena.
which we expand to first order in δα considering α ∼ const.. This is arranged so that we can use small transformations of the Hamiltonian and develop them, stepby-step, into a full profile θ(x, t), as discussed in the Introduction. Importantly, θ(x, t) may have a strongly fluctuating space-time dependence (like a step function). For the full profile θ, we then have Path I (14) with ∆ = 0, m 5 = m, and introduce a path connecting the different insulating phases as in (17), which we repeat here for convenience for the reader,

Consider Hamiltonian
with iγ 5 = Γ 2 Λ 5 = iΓ 1 . Now let us introduce a small change in α → α + δα and unwind the mass phase by rotating the spinors as Equality of the partition functions before and after the fermion rotation leads to the "naive" chiral Ward-Takahashi identity implying the conservation of the chiral current for vanishing TI mass. We know that this conservation is actually spoiled due to the chiral anomaly, captured by a non-invariance of the Jacobian. In the present approach, α is considered, overall, a constant and builds no Jacobian. On then just needs to compute the effective action in gradient expansion in δα (Goldstone-Wilczek calculation.) Up to the neglected gradients of α, the Hamiltonian thus reads: The Lagrangian becomes The resolvent reads and the potential is Terms in the effective gauge action at first order in δα give the (topological) contributions to the average of the chiral current. After integrating out the fermions with (A6) in the path integral, the effective action to second order in the potential V is given by Keeping one contribution from the dynamical δα, and one from the gauge field A, reduces the effective action to The trace over the spin matrix σ 0 is trivial and gives a factor of 2. Also, the only non-vanishing trace is the one with τ x τ y τ z , with tr [γ µ γ ν γ 5 ] = 4 µν . We find tr [G 0 (k) γ µ G 0 (q + k) γ 5 ] = −tr γ σ k σ + m k 2 − m 2 γ µ γ ν (q + k) ν + m (q + k) 2 − m 2 γ 5 To the lowest order in q (gradient expansion) we have δS topo = 4m 2 µνˆq A µ (−q) δα (q) q νˆk 1 in agreement with (18).
In computing the approximate Green function we use that, given a matrix A and a unitary transformation U that diagonalizes it as U AU † = D, we can write and since we can use Taylor expansion to every eigenvalue of D (or 1 + D), The Green function is, to first order in ∆, Now we notice that Γ 5 and γ 5 commute and may be diagonalized simultaneously. So we may use (A18), with G 00 = γ µ k µ + m k 2 − m 2 and Also, the potential may be written This allows us to separate the first order in ∆ (over m) correction to the action We have shown that δS 0 topo gives (twice) the expected theta-term. Our goal is to see if there is any correction due to δS 1 topo . Realizing that V φ1 = σ z −i∆γ 5 e −iγ5φ cos φδφ G 01 = σ z (Φ 1 γ 5 − Φ 2 (1 + 2mG 00 )) it is easy to see that, since the traces over the σ and the τ matrices decouple, δS 1 topo vanishes (since the σ z traces vanish).