{'fraction_non_alphanumeric': 0.03572258814413018, 'fraction_numerical': 0.011959188944853417, 'mean_word_length': 4.355650850740075, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electromagnetic response of quantum Hall systems in dimensions five and six and beyond

Quantum Hall (QH) states are arguably the most ubiquitous examples of nontrivial topological order, requiring no special symmetry and elegantly characterized by the first Chern number. Their higher dimension generalizations are particularly interesting from both mathematical and phenomenological perspectives, and have attracted recent attention due to a few high profile experimental realizations. In this work, we derive from first principles the electromagnetic response of QH systems in arbitrary number of dimensions, and elaborate on the crucial roles played by their modified phase space density of states under the simultaneous presence of magnetic field and Berry curvature. We provide new mathematical results relating this phase space modification to the non-commutativity of phase space, and show how they are manifested as a Hall conductivity quantized by a higher Chern number. When a Fermi surface is present, additional response currents unrelated to these Chern numbers also appear. This unconventional response can be directly investigated through a few minimal models with specially chosen fluxes. These models, together with more generic 6D QH systems, can be realized in realistic 3D experimental setups like cold atom systems through possibly entangled synthetic dimensions.

While intrinsically 2D phases, QH states can, through special choices of gauge, also be expressed as independent 1D quasiperiodic chains labeled by a synthetic dimensional parameter. This fortuitous dimensional reducibility has allowed for realistic 2D implementations of higher dimensional QH generalizations [51][52][53][54][55]: recent experiments [1,2] have demonstrated signatures of the 4D QH effect with ultracold bosonic atoms in optical superlattices and tunable photonic waveguides. In these setups, a d-dimensional QH system was constructed simply by taking the tensor product of d 2 copies of 2D QH systems [52], with the resultant d-dimensional magnetic flux decomposable into a product of ordinary 2D fluxes. One can alternatively construct nontrivially "entangled" higher-dimensional QH states by introducing additional off-diagonal couplings between the physical and synthetic dimensions [55], as shall be detailed in Sect.IV B. So far, existing proposals and experimental realizations have been limited to 4D QH phases, although QH phases in 5 or more dimensions are just as realizable in principle [56,57].
In all time-reversal broken systems including QH fluids and Chern insulators, small changes in the magnetic field lead to non-perturbative modifications of the lattice band structure due to non-commensurability effects, as epitomized by the Hofstadter butterfly [58]. Yet, their physical electromagnetic responses are of course only perturbatively modified [52,[59][60][61][62][63]. This suggests the strategy of separating a generic time-reversal breaking flux into an "intrinsic" Berry curvature that redefines the band structure, and an "external" magnetic field that only enters the physical response as a perturbing field [52,64]. At the level of semiclassical response dynamics, this simultaneous presence of magnetic field and Berry curvature seemingly leads to a subtle but interesting violation of Liouville's theorem, which can be remedied either by redefining the phase space density of states [65], or by resorting to non-canonical coordinates [66,67]. Both of these resolutions shall give rise to especially profound physical consequences in higher dimensions, as we will detail in this work. Specifically, the interplay of magnetic field and Berry curvature not only produces a quantized third Chern number Hall response in 6 or more dimensions, but also yields new current responses when a Fermi surface is present.
Our paper is structured as follows. In Sect. II, we first derive the semiclassical equations of motion (EOMs) of electronic wavepackets in the presence of both magnetic field and Berry curvature, starting from their Lagrangian. Next we derive new explicit results for the concomitant modification of the phase space density of states in arbitrary number of dimensions. Combining these results, we present the current response terms and re-derive the latter in terms of topological pumping. Following that in Sect. III, we present minimal models for studying interesting non-topological surface current contributions in 5 or 6 dimensions when a Fermi surface is present. We next discuss how these 6D QH systems can be realized with 3D physical systems with 3 synthetic parameters in Sect. IV, before reviewing specifics of their experimental realization in cold atom systems in Sect. V.

II. SEMICLASSICAL TRANSPORT OF ELECTRONIC BANDS WITH NONZERO BERRY CURVATURE AND FLUX
We consider generic d-dimensional electronic systems subject to external perturbing electric and magnetic fields E(r) = E µ (r)e µ = −∇ µ Φ(r) e µ and B µν (r) = ∂ µ A ν (r) − ∂ ν A µ (r), where Φ(r) and A(r) = A µ (r)e µ are the scalar and vector potentials respectively. In dimensions d > 3, the magnetic field cannot be visualized as a vector, but is rather a 2nd-rank tensor field with d(d − 1)/2 independent components. Introducing a lattice regularization if necessary, the energy eigenstates of the system splits into distinct bands on the momentum-space d-torus T d . When intrinsic time-reversal or inversion/particle-hole asymmetry is present, each band |ψ will also acquire a Berry curvature Ω µν (k) = ∂ µ a ν (k) − ∂ ν a µ (k) = i( ∂ kµ ψ|∂ kν ψ − ∂ kν ψ|∂ kµ ψ ), where a µ = −i ψ|∂ kµ ψ is the Berry connection.
In this work, we shall assume non-degenerate bands, such that each band (partially or completed filled) can be treated separately. The current response j = j µ e µ is given by where the electromagnetic fields E, B µν and the Berry curvature Ω µν implicitly affect the quantities in the square parentheses. In the following subsections, we shall detail how exactly they affect the semiclassical wavepacket velocityṙ µ , as well as the phase space density of states D(r, k). The former will be described at the level of classical equations of motion in Sect. II A, and recast in terms of topological pumping in Sect. II D.

A. Semiclassical equations of motion
We first attempt to understand the response properties of an electronic system from the equations of motion of its semi-classical electron wavepackets. We consider wavepackets which are much larger than the lattice spacing (∼ 1Å), such that each of them possesses a single well-defined peak in momentum space. The wavepackets must also be such smaller than spatial variation of the magnetic vector potential A( r), which hence restricts the external magnetic field to be a weak perturbation. The external electric field must also be sufficiently weak for the wave-packet to evolve adiabatically without interband transitions.
Labeling the spatial and momentum-space peaks of a wavepacket as r = r µ e µ and k = k µ e µ , one can write down the Lagrangian (with e and set to unity) through minimal coupling of both the magnetic and Berry gauge fields A µ (r) and a µ (k), with the Hamiltonian H(r, k) simply given by the sum of the Bloch Hamiltonian E(k) and the electrical potential −Φ(r). Solving the Euler-Lagrange equations with Eq. 2, we obtain the equations of motion (EOMs) for the wavepacket [68] (Einstein summation implied from now on) where v µ = ∂E(k) ∂kµ is the group velocity computed from the wavepacket energy of the unperturbed Bloch band. In the above, the momentum and positional peak positions take on symmetrical roles, being linearly accelerated by position-momentum duals v µ and −E µ , and "curved" by analogously dual fields Ω µν and B µν . This position-momentum symmetric formulation pertains to wavepackets in any number of dimensions d, and is one of many examples of position-momentum duality. In related contexts, position-momentum symmetry has been employed to prove results in the free-fermion entanglement spectra of spin textures and thermalized systems [69,70]. Very prevalently, it has also allowed an elegant formulation of 2D QH dynamics in terms of cyclotron and guiding-center coordinates [71,72], which has incidentally been exploited for relating fractional QH phases in lattice and continuum systems [24].
By repeatedly substituting Eq. 4 into Eq. 3, we can obtain the equation of motion for the center-of-mass (CM) velocityṙ up to any desired level of accuracy. For instance, up to the third order in E µ and B µν , Since our eventual purpose is to obtain only the leadingorder terms of each type of transport contribution, we have neglected field-induced shifts in the wavepacket CM, which can be expressed in terms of the quantum geometric metric as in Ref. [73].

B. Modification of density of states
A central foundation of Hamiltonian mechanics is Liouville's theorem, which states that the phase space volume must remain conserved under time evolution, i.e. that d(ln ∆V ) dt = ∇ r ·ṙ + ∇ k ·k should vanish for a suitably defined volume form ∆V . With the EOMs in Eqs. 3 and 4, it is easy to show, at least in 2D [65,66], that the simultaneous existence of Berry curvature and magnetic field lead to non-conservation of the naively defined volume element ∆V naive = ∆r∆k. To derive the properly conserved volume element and remedy this apparent conundrum, we rewrite the EOMs Eqs. 3 and 4 in the language of Hamiltonian dynamics: where {u, v} is the Poisson bracket defined by {u, v} = ω µν ∂ µ u ∂ ν v. ω µν is the inverse of the matrix ω µν , which defines the closed symplectric structure of this Hamiltonian system, and famously gives the volume element that will indeed remain conserved under time evolution [66]. Pf denotes the Pfaffian, the square root of the determinant operator on an even-dimensional antisymmetric matrix. To compute the modification factor Pf(ω µν ), we contract Eq. 6 with ω µν = (ω µν ) −1 and obtain ω µν as the real antisymmetric transformation matrix between (ṙ,k) and (E, v), i.e.
where B and Ω are the magnetic and Berry curvature tensors respectively. With a detailed derivation left to Appendix I, we obtain our main result on the modification to the density of states: where ε µ1µ2...µ d , ε µ1µ2...µ d are generalized Levi-Civita symbols which are equal to the sign of the permutation {µ 1 , ..., µ d } for unique µ j indices, and vanish if there is any repeated µ j indices. As such, for a given number of dimensions d, only 2j-th order terms with j ≤ d 2 exist. Eq. 9 generalizes previous results [52,65] to arbitrary number of dimensions d, with expressions for the first several d given by In the above, there is an implicit sum over all possible subsets of indices whenever the there are more dimensions than indices, i.e. the term containing ε αβγδ will be summed over all 4-element subsets α, β, γ, δ ∈ {1, 2, 3, 4, 5} in the 5D case. This is consistent with the physical intuition that additional synthetic dimensions should only result in new response terms, and not modify existing terms.
C. Leading order current response terms Having discussed the semiclassical EOMs and how they crucially modify the density of states, we are now ready to present results for the electromagnetic response.
In four or fewer dimensions, it is well-established [52] that the response of an insulating electronic system (with completed filled band) is completely quantized, with a Hall response term proportional to the 1st Chern number, and a magneto-electric response term proportional to the 2nd Chern number. Analogous results hold in higher dimensions. Specializing to six dimensions in anticipation of the physically realizations discussed later, a 6D insulating system has a bulk response current: In the above, C µν 1 , C µαβν 2 and the topologically quantized third Chern number C 3 are given by Note that C µν 1 and C µαβν 2 are computed with respect to {µ, ν} and {µ, α, β, ν} hyperplane subsets of the 6D space respectively. Eq. 13 can be put into explicit covariant form if we were to explicitly write it in terms of the generalized Levi-Civita symbols.
In the metallic case (with a Fermi surface ∂Γ), it is known for d ≤ 4 that there exists a distinct surface term producing a non-topological surface current [74] j surf , in addition to corrections to the Chern number contributions from Eq. 13. For five or more dimensions, this surface current takes a much more sophisticated form. Firstly, j µ bulk ceases to be topologically quantized, with its integration region bounded by ∂Γ. Importantly, the total response current j µ = j µ bulk + j µ surf | T 6 →Γ now contains a new contributionj µ surf given by As will be explicated in Appendix II, both the surface current j surf and the topological 3rd Chern number C 3 arise from a careful combination of EOM and density of states modifications (Eqs. 5 and 9), though with the l-th order correction entering only the j surf in 2l + 1 dimensions. Note that Eq. 17 includes only leading (quadratic) order contributions in the magnetic field; we leave the study of its higher order corrections to future work.

D. Semiclassical topological pumping
After examining the wavepacket equations of motion and deriving its semiclassical transport current, which in general does not have to be parallel to the applied fields, it is instructive to re-derive these results for some lower dimensions in terms of topological (Thouless-like) pumping arguments. For the sake of clarity, in this subsection we shall adopt somewhat different notation from the rest of this paper.
Consider the motion of a Bloch wavepacket centered around momentum k and real-space coordinate X k under the influence of a Berry curvature Ω. Semiclassical theory tells us that X k responds to an external electric field E = E s e s according toẊ k = −k × Ω, where k s → k s + E s t from minimal coupling. Hence, we have where the superscript s has been added to emphasize that the time dependence of X is due to an electric field in the directionŝ. This is just the Hall response of wavepackets.

Spectral flow for 2D QH
We next generalize the above Hall response to spectral flow due to topological pumping. For a 2D QH system in the x-s plane, the spectral flow of Wannier function centers is directly related to the Hall pumping of its constituent wavepackets. We write the time-dependent Wannier center coordinate as X(t), which is equal to the Wannier polarization P (k s ) upon the minimal coupling k s → k s + E s t: In the 2nd line, I have used the standard expression for the Wannier polarization [75][76][77]] P X of X(t) due to k s , which for macroscopic displacements can be linearly approximated as in line 3. In the second last line, I have equivalently expressed X(t)−X(0) in terms of the Chern number C xs 1 .
Not surprisingly, the Wannier center X(t) is given, in the last line, by the integral of the centers of all wavepackets in the BZ.

Spectral flow for 4D QH
We next illustrate how to generalize this pumping argument to 4 dimensions. Notably, we shall see that the modification of the density of states must be included to recover the correct response behavior given by the Chern numbers.
With 4 dimensions, there exists nontrivial pumping contributions up to 2nd order in the perturbing fields E and B. Let the 4D QH system span coordinates x, y, s and u. Qualitatively, we know that the combination of electric field E u and magnetic field B ys can lead to the pumping of X(t) through three possible mechanisms: 2. Second order Hall responses mediated by the magnetic field via k s → k s + B sy Y (t) and k y → k y + B ys S(t), where Y (s) and S(t) are themselves pumped by E u , and 3. Second order corrections from phase space modification due to the magnetic field.
Explicitly, we have The above derivation began by considering all the above three momentum-resolved wavepacket contributions to X(t), and arrived at the final results via repeated applications of Eq. 19. The phase space modification factor D d=4 = 1 + B µν Ω µν + higher order ≈ 1 + B ys Ω ys is seen to play a crucial role in producing the term containing Ω ys : without it, the integral giving rise to the 2nd Chern number C 2 will not be complete or even symmetric. Note that higher order corrections must disappear in 4D systems due to antisymmetry. Indeed, we recover the 2nd-order response equation behavior (Eq. 13) derived earlier. While the 1st Chern number term arises simply from the direct Hall response, the 2nd Chern number term arises via a combination of two different second order Hall responses channels, as well as the density of states correction term B ys Ω ys . While the latter contribution cannot have possibly originated from any Hall pumping, its appearance is necessary for constructing the integrand for the 2nd Chern number.

Further interpretations for unentangled higher-dimensional QH
In the case of unentangled (factorizable) 4D states that can be expressed as a tensor product of 2D QH states in x-s and y-u space, only the product Ω xs Ω yu in the integrand of C 2 survives. In this case, the Wannier center flow subject to an electric field E u and magnetic field where k u = E u t by minimal coupling, and P X , P Y are the expressions for Wannier polarization in the X and Y directions. To isolate the 2nd Chern number term from Eq. 21, one can take differences between separate measurements to obtain For the 6D entangled case, one can analogously consider an electric field E v viz. k v = E v t, and magnetic fields B 1 = B sy , B 2 = B uz . The response of the averaged wavepacket due to all 1st, 2nd and 3rd Chern numbers is thus given by Similarly, we can take differences and obtain the final linear approximation holding for large k v (t) = E v t.
E. Non-commutativity and phase-space modification From the above discussions, it has been apparent that electronic transport does not just arise from "bare" electron dynamics, but is also affected by how the phase space distorts around it. This distortion can be interpreted as the "smudging" of the electron itself. Promoting the Poisson brackets in Sect. II B to operator commutators, we see that To leading order, we simply obtain [r µ , , which sets a lower bound for the amount of uncertainty in the phase-space coordinates of an electron, i.e. how much it is "smudged" due to the magnetic and Berry curvature fields. As expected, a magnetic field B causes the electrons to move in space cyclotron orbits, thereby smudging it in momentum-space. Analogously, a Berry flux Ω smudges it in real-space, contributing to a topological obstruction for localized Wannier functions. Nontrivial corrections to these occur precisely when both B and Ω are present. In particular, there will be smudging across the real and momentum space coordinates, with electrons effectively executing phase-space "cyclotron" motion that leads to corrections to the transport.
Due to the antisymmetry of both B and Ω tensors, it can be shown that ω µν and hence all commutators are the quotient of the sum of products of matrix elements of B and Ω with D(r, k) = Pf(ω µν ) = Det(I − BΩ). For instance, where D µ is the density of states factor in the space with dimension µ omitted. In general, the commutators have rather complicated expressions in arbitrary number of dimensions involving the explicit form of (I − ΩB) −1 or (I − BΩ) −1 . Closed form exact expressions can be easily derived for low dimensions however. For d = 2 dimen-sions, [r µ ,r ν ] d=2 = iΩ µν 1 + Ω 12 B 12 (27) [k µ ,k ν ] d=2 = −iB µν 1 + Ω 12 B 12 (28) [r µ ,k ν ] d=2 = iδ µ ν 1 + Ω 12 B 12 (29) as appearing in Refs. [65] and [66]. For d = 4 dimensions, where D d=4 (r, k) is given in Eq. 11. These exact expressions should be contrasted with those in the last line of Eq. 25, which do not contain phase space correction denominators.
FIG. 1. Cartoon picture of electrons moving under the simultaneous influence of magnetic field and Berry curvature. They exhibit "cylotron" orbits in phase space, thereby behaving like "smudged" particles with uncertainties in their coordinates proportional to the commutator between their phase space coordinates.

III. MINIMAL MODELS FOR NOVEL SURFACE CURRENT j µ surf IN 5D AND 6D
In this section, we examine the properties of the new non-topological surface currents through some minimal models with a Fermi surface. These are the simplest models that yield a nonvanishing 2nd order j µ surf , and can in principle be realized with any tight-binding model by placing nonzero flux across the stipulated plaquettes. Note that we will only be focusing on the current contributions of 2nd order in B and Ω, i.e. second and third lines of Eq. 17, which appear only in 5 dimensions and beyond.
A. 5D case Consider a 5D electronic metallic system consisting of 3 spatial dimensions {x, y, z} and 2 synthetic dimensions {s, w}. One possible minimal model consists of just two nonzero components of the magnetic field B µν , and three nonzero components of the Berry curvature Ω µν . Namely, we have uniform magnetic fields and uniform Berry curvatures This configuration gives a non-topological current density j µ surf,5D in the e z direction: We see that j z surf,5D should exist for Γ = T 6 as long as the Berry curvatures are not all uniform, which is generically the case. Note, though that the gradients of Ω, if not already zero, have to be related by the Bianchi identities (see Ref. 25 for a systematic way of making the Berry curvatures uniform). B xs and Ω xz can be implemented as external and intrinsic physical fluxes, while B yw , Ω ys and Ω zw can be implemented as phase factors associated with synthetic parameters y and w. Realistic systems will in general possess other non-vanishing fluxes that change the direction of j z surf,5D , but they will not cause the latter to vanish when the Fermi surface vanishes.

B. 6D case
Analogous to the 5D case, a minimal model for the 6D metallic case with 3 spatial dimensions {x, y, z} and 3 synthetic dimensions {s, w, t} is given by uniform magnetic fields and uniform Berry curvatures Ω zt = −Ω tz (k z , k t ) = 0, which gives rise to a non-topological surface current contribution Both nonzero current components are almost identical, except for the swapping of the roles of x and t components. As before, j surf,6D should not vanish unless the Fermi surface vanish.

IV. REALIZATIONS OF 6D QH STATES IN 3 PHYSICAL DIMENSIONS
Having developed the response formalism for higherdimensional QH systems, we now present a general approach for realizing them in 3-dimensional physical space.
Building going on to 6D QH models that are realized in 3-dimensional space, we first describe how 2D QH models can be realized in one physical dimension plus one synthetic dimension. In Refs. [55,[78][79][80], it was shown that a 2D lattice system under the influence of a magnetic field (Hofstadter model) can, with suitable gauge choice, be expressed as a so-called 1D Aubry-André-Harper (AAH) model with a synthetic parameter. To understand this, consider a square lattice with horizontal (x-direction) hoppings t and vertical hoppings λ. If we choose a gauge such that each vertical hopping (but not horizontal hopping) contains a gauge phase of e 2πibx , the entire lattice will be threaded with a uniform (internal) magnetic field with each square plaquette experiencing a flux of 2πb. Due to translational symmetry in the vertical direction, we can perform a partial Fourier transform such that k represents the lattice momentum in the vertical direction. The hopping Hamiltonian hence looks like which can be interpreted as a collection of decoupled 1D AAH chains given in the penultimate line. Each AAH chain has an onsite energy term with periodicity 1/b and offset k, and a uniform nearest-neighbor hopping t. For a rational flux b = p q with GCD(p, q) = 1, the AAH Hamiltonian splits into an effective q band Hamiltonian. In this sense, a change of magnetic field, no matter how small, changes the band structure non-perturbatively. Indeed, when the flux b becomes irrational, H AAH 1D becomes a quasiperiodic chain with dispersionless Landau Levels. That said, the physical responses remain perfectly smoothly defined, as derived in the previous section where the fields are taken as external perturbations.
A. Unentangled 6D QH system with 3rd Chern number The most direct way of achieving a 6D QH system is to construct it as a tensor product of three 2D QH systems, each represented by a 1D AAH model, i.e.
whose energy spectrum is shown in Fig. 2 for flux b = p/q = 1/4, i.e. with each magnetic unit cell consisting of q = 4 sites. We observe clear localization in the x-surface, in fact at its corners due to C4 symmetry.
As a tensor product of three AAH models, this 6D QH state possess a 3rd Chern number [60] that is the tensor product of three 1st Chern numbers: with the 3 spatial dimensions and their 3 corresponding synthetic dimensions labeled as x, y, z and s, w, t respectively.
B. Entangled 6D QH system with 3rd Chern number Although the 6D system described by Eq. 41 is a bonafide 6D QH system with electronic response given by the earlier sections, its spectral flow properties can be entirely understood in terms of those of its constituent 1st Chern number systems. To generate a 3rd Chern number system that is not such a trivial tensor product, the simplest approach is to mix (or entangle) their directions of spatial modulation [55]: where R is a orthogonal matrix. Since the t-couplings that define the 3D lattice remain unchanged, the R rotation effectively couples the different copies of the AAH models. If the R matrix is replaced by a generic nondegenerate 3-by-3 matrix, we will obtain a 6D QH system that is mixed in a nontrivial though not necessarily symmetrical manner.

V. PHYSICAL REALIZATION AND TOPOLOGICAL RESPONSE MEASUREMENT OF 3RD CHERN NUMBER
In the previous section, we have discussed how periodically modulated AAH systems with synthetic dimensional offsets can possess higher Chern numbers. In the quasiperiodic limit, the Chern bands become flat Landau levels, effectively simulating higher-dimensional QH systems with uniform magnetic field. It is important to note, however, that quasiperiodicity is not a necessary ingredient -Chern bands can emerge even with a relatively small unit cell, as long as the modulations are defined to result in nonvanishing effective flux. For realistic implementation, we consider a 3D optical lattice of ultracold atoms, the 2D version for which has already been experimentally realized [1,81]. Consider an optical potential that is modulated by two sinusoidal profiles in each direction, as illustrated in Fig. 4: where R is the rotation matrix introduced in Sect. IV B, b, b are dissimilar inverse periods and φ i , φ i are modulation potential depths which can be different in different directions i. Due to the interference between the two modulations, we obtain an effective flux that leads to higher Chern bands, as elaborated in Ref.
[1]. When R = I 3×3 , i.e. in Eqs. 43 or 44, the constituent Chern systems are coupled, resulting in an entangled higher Chern (or QH) phase. The synthetic momenta are implemented by the phase offsets k i , which can be controlled by slightly tilting the optical lattice.
For simplicity, 1/b can be taken to be the spacing of the 3D lattice (Fig. 4b), while 1/b shall correspond to an imposed modulation period (Fig. 4c). The potential φ i is given by the blue potential in Fig.4c, which fixes the atoms on lattice sites, while the potential φ i introduces the local modulation as given by the red potential. Taking the s-wave approximation, Eq. 45 describes a spatially modulated tight-binding Hamiltonian similar to Eq. 44, with effective flux of magnitude 2πb eaa in a direction determined by the matrix R, where e is the electronic charge and a, a are the lattice constants in the two directions connected by R.
To detect the higher-dimensional QH transport properties, an external electric field E µ has to be introduced into the cold-atom system. For that, one can implement  Fig. 4(b), with each sphere denoting a cold atom. b) The system is modulated periodically by the sum of two sub-potentials (orange), as given by Eq. 45. The blue potential denotes the on-site modulation we set to the system and the red potential denotes the potential that pins the cold atoms to the lattice sites.
a linear gradient [82] either magnetically [83,84] or optically [85]. The current, whether j surf or the topological contribution, can also be indirectly measured via the center-of-mass velocity of an atomic cloud after accounting for the density of states correction [86,87]. For illustration, we focus on the lowest 6D band, which for appropriate choices of uniform magnetic fluxes is nondegenerate and well-isolated from the higher bands. Suppose all the Berry curvature components are zero except for: Ω zt = −Ω tz (k x , k s ) = 0, We will choose, with no loss of generality, the perturbing (synthetic) electric field to be along the s direction i.e. E = E s e z as in Fig.4 a. However, in 6D, different choices for the perturbing magnetic field will lead to dramatic differences in the observables. If we choose the only nonzero external magnetic field components to be B xw and B yt , the current density j in this experimental setup is explicitly given by: which is schematically shown in Fig.4a. The various Chern numbers C 1 , C 2 and C 3 can then be independently obtained by measuring the various current density components.

VI. CONCLUSION
Starting from the very first principles, we have derived electromagnetic response properties of arbitrarilydimensioned QH systems through the complementary approaches of semiclassical wavepacket dynamics and topological pumping arguments. Together, they illustrate the crucial role of a modified phase space density of states when magnetic field and Berry curvature are simultaneously present. This modification becomes all the more nontrivial in higher dimensions, where it provide the crucially missing pieces in the Chern invariants for the topological response. We derived some new mathematical results for it and interpreted them in terms of noncommuting phase space dynamics.
Notably, we have explored new surface current contributions j surf in dimensions greater than four, which appear in addition to the usual topological terms. We provided minimal models for investigating this enigmatic term, at least up to its leading order. These models, as well as the paradigmatic 6D QH systems, can be realized in realistic 3D experimental setups like cold atom systems through possibly entangled products of AAH chains.
Subsequent to the initial appearance of our manuscript on arXiv, we become aware of a related work (Ref. 88) on higher-dimensional QH systems, which provided a more in-depth description of experimental realization using three-dimensional topological charge pumps in cold atomic systems.

APPENDICES I. METHOD FOR COMPUTING THE DETERMINANT OF A SPECIAL TYPE OF ANTISYMMETRIC MATRIX
In deriving Eq. 9 from the Pfaffian of Eq. 8 in the main text, we need to find the determinant of matrices of the form We will adopt a mathematical induction approach for this computation, and exploit the the antisymmetry of M and the symmetry between the indices of X aibj and Y aibj . First, we give the conditions that the expansion of the determinant should satisfy: • The expansion should be a perfect square.
• X aibj and Y aibj should be symmetric and noninterfering.
• There should not be any odd-ordered term in the expansion.
• Suppose the power of the expansion is m, and the dimensionality of the matrix is n = 2k (k is a nonnegative integer,because obviously n is even), then m = 2(k − 1), k is a nonnegative odd integer; 2k, k is a nonnegative even integer.
The first condition is based on the antisymmetry of the matrix. As we know, the expansion of any antisymmetric determinant must be a square of a polynomial [90] which is called the Pfaffian of the matrix. And it's apparent to get the second condition by observing the characteristics of this matrix, which proclaim that X aibj and Y aibj are symmetric to some extent. The expansion of the determinant stays the same in the case of exchanging all the X and Y with same indexes. Changing any elements of the determinant will not influence the form of the other elements in the expansion. As for the third condition, there has to be one zero at least in the odd-ordered terms, and thus they will not exist in the expansion. The last condition are satisfied when k is even because apparently the power of the expansion equals the dimensionality of the matrix. To discuss the odd-k situation, we need to find out why the highest order term vanishes. Since the highest order term of the determinant will not contain 1, −1 or 0 factors in its determinant expansion, we can deduce that it is decided by these elements: Hence we can find that there are just four ways to form a highest order term. However, considering the definition of determinant: det(M ) = ε i1···in a 1,i1 a 2,i2 · · · a n,in (50) where ε i1···in is a Levi-Civita symbol, we can get these four highest order terms h 1 , h 2 , h 3 and h 4 , which have the same absolute values: and when k is even, these four terms are all nonnegative; When k is odd, two of these four terms are nonnegative which are the opposite number of the other two terms. Thus it's obvious that these four terms cancel each other out, which causes the power of the matrices with an odd k decreases accordingly. Before giving the conjectured solution, we need to introduce the concept of a general Levi-Civita symbol, whose indices don't have to be a permutation consisted of continuous natural number. With α 1 , α 2 , α 3 , · · · α n ∈ N + and p the parity of the permutation in the above sequence, p is odd; 0, if any two indices are equal.
As we have established the four conditions that the expansion should satisfy, we shall now put forward the conjecture solution: ,2) · · · X α ( n where n is the dimensionality of the original matrix M and the set that consists of the indices of X equals the one that consists of the indices of Y in the same term. The reason why we leave an ellipsis behind the whole expansion is that in fact, there exist infinite terms in the expansion but the terms behind all vanish due to the definition of Levi-Civita symbol, which means we only need to calculate this formula to the order that we require, and α ij and β ij satisfy: Moreover,the elements of the matrix X aibj satisfy: X aibj = −X bj ai (55) and f (x) = x is the floor function.
Since we have allowed the dimensionality of the matrix to be any arbitrarily large number, we shall employ mathematical induction to prove that this expansion is exactly the right one.
First, when n = 1, 2, 3, it's apparent that the expansion equals the determinants of the matrices involved. Provided that the proposition is valid when n = k, for n = k + 1 we have: which shows that when: i, j ∈ Z + , i ≤ k + 1, j ≤ k + 1 then: where M n is a n-dimensional antisymmetric matrix derived from zero-padding M when necessary. Apparently, the determinants of M n and M are equal. Hence we only need to consider the extra terms that these new elements in the matrix bring to the determinant. With that, the expansion of det(M 2(k+1) ) must be of the form: where P (X aib k+1 , Y aj b k+1 ) is the polynomial that contains X aib k+1 and Y aj b k+1 (i, j ∈ Z + , i ≤ k + 1, j ≤ k + 1).
Combined with the four rules we mentioned above, we can calculate the k + 1-th order determinant in this way: first we keep the absolute values of all the new Y terms be 1,and whether the specific term is positive or negative depends on its location in the determinant. In our case, the terms above the diagonal terms are positive and the terms below the diagonal terms are negative, which is like: Then we can expand the determinant according to the definition. We know that the result we get is not the general one, which is not easy to get. However considering the symmetry between X terms and Y terms, we can get the general result by adjust the formula artifi-cially, which means that we complement the formula to a symmetric one. For instance, if we have got a term that contains X a1b k+1 , which might as well be named T 1,k+1 , there must exists a term that contains exactly the same elements as T 1,k+1 except changing X a1b k+1 into Y a1b k+1 . Hence we can get: where one of the indices of each term must be the additional (k + 1)th index. Thus we can prove the correctness of our conjectured solution. Back to the problem of the modified phase-space density of states, we just need to consider a matrix of the form which we now know must have a Pfaffian given by Eq. 9.

II. DERIVATION OF 6D CURRENT DENSITY AND TRANSPORT EQUATIONS FROM THIRD-ORDER SEMI-CLASSICS
Here, we fill in the details of the derivation of Eq. 13, specializing to D = 6 dimensions for concreteness. To calculate the response current density of a filled band, we set the integration (phase space) measure to be Ω * i (α,β,γ,δ) =    ε P i1 P i2 P i3 P i4 Ω P i1 P i2 Ω P i3 P i4 , ε P i1 P i2 P i3 P i4 = 0; 0, ε P i1 P i2 P i3 P i4 = 0 P is the matrix of permutations of α, β, γ and δ, i.e.
With the above analysis, one sees that although it has a complicated form, it turns out to be zero under certain circumstances. For instance, if the energy E is an even function of k, each term in the integral will be zero. Also, using the setup in our main text, we find this term turns out to be: where j z surf and j s surf turns out to be zero after integrating, too.
The next set of terms Γ d 6 k (2π) 6 1 64 E ν Ω µν ε αβγδ B αβ B γδ × ε θνλρ Ω θν Ω λρ + 1 2 E δ Ω γδ B νγ Ω µν B βθ Ω βθ + E β Ω αβ B δα Ω γδ B νγ Ω µν can also be simplified using the antisymmetry of the magnetic field strength (B γν = −B νγ ) and the Berry curvature (Ω γν = −Ω νγ ), and consequently for a vanishing Fermi surface reduces to 1 8 where C 3 is the 3rd Chern number as defined in the main text. The final group of terms Γ d 6 k (2π) 6 1 1024 v µ ε ζηθτ κξ B ζη B θτ B κξ × ε σωιφχψ Ω σω Ω ιφ Ω χψ fortuitously all cancel to zero by virtue of the antisymmetry of the magnetic field and of the Berry curvature in 6D, reminiscent of the j surf terms which cancel totally in 4D due to the same symmetries. We conjecture that this final set of terms will contribute to a higher j surf in dimensions greater than 6, although proving that becomes rather tedious. Combining all the above results, we obtain the expressions in Eqs. 13 and 17 in the main text.