Non-Abelian Bloch oscillations in higher-order topological insulators

Bloch oscillations (BOs) are a fundamental phenomenon by which a wave packet undergoes a periodic motion in a lattice when subjected to an external force. Observed in a wide range of synthetic lattice systems, BOs are intrinsically related to the geometric and topological properties of the underlying band structure. This has established BOs as a prominent tool for the detection of Berry phase effects, including those described by non-Abelian gauge fields. In this work, we unveil a unique topological effect that manifests in the BOs of higher-order topological insulators through the interplay of non-Abelian Berry curvature and quantized Wilson loops. It is characterized by an oscillating Hall drift that is synchronized with a topologically-protected inter-band beating and a multiplied Bloch period. We identify the origin of this synchronization mechanism through a quantum dance of Wannier centers. Our work paves the way to the experimental detection of non-Abelian topological properties in synthetic matter through the measurement of Berry phases and center-of-mass displacements.

T he quest for topological quantization laws has been a central theme in the exploration of topological quantum matter 1,2 , which originated from the discovery of the quantum Hall effect 3 . In the last decade, the development of topological materials has led to the observation of fascinating quantized effects, including the half-integer quantum Hall effect 4 and the quantization of Faraday and Kerr rotations 5 in topological insulators 6,7 , as well as half-integer thermal Hall conductance in spin liquids 8 and quantum Hall states 9 . In parallel, the engineering of synthetic topological systems has allowed for the realization of quantized pumps 10,11 , and revealed quantized Hall drifts [12][13][14] , circular dichroism 15 , and linking numbers 16 .
In this context, Bloch oscillations (BOs) [17][18][19][20][21][22] have emerged as a powerful tool for the detection of geometric and topological properties in synthetic lattice systems [23][24][25][26] , hence providing access to quantized observables. Indeed, transporting a wavepacket across the Brillouin zone (BZ) can be used to explore various geometric features of Bloch bands, including the local Berry curvature 23 and the Wilson loop of non-Abelian connections 25 . This strategy has been exploited to extract the Berry phase 27,28 , the Berry curvature 29,30 , the Chern number [12][13][14] , and quantized Wilson loops 31 in ultracold matter and photonics.
The Wilson loop measurement of ref. 31 highlighted a fundamental relation between two intriguing properties of multi-band systems: the quantization of Wilson loops, a topological property related to the Wilczek-Zee connection 32 , and the existence of "multiple Bloch oscillations", which are characterized by a multiplied Bloch period 31,[33][34][35][36][37][38][39] . The effect investigated in ref. 31 was eventually identified as an instance of "topological Bloch oscillations", whose general framework was proposed in ref. 40 based on the space groups of crystals and its implications on the quantization of geometric quantities (Zak phases differences). The Bloch period multiplier appears as a topological invariant, protected by crystalline symmetries, thus making multiple BOs genuinely topological. Furthermore, when a Wannier representation of the bands is possible, Zak phases correspond to the positions of charges within the unit cell, namely the Wannier centers. As a consequence, a Zak-Wannier duality allows to connect BOs to the relative phases acquired by charges within a classical point-charge picture. More recently, BOs displaying topologically protected sub-oscillations have also been found in periodically driven systems in the context of quantum walks 41 .
In this work, we identify a distinct topological effect that manifests in the BOs of higher-order topological insulators (HOTIs). These newly discovered systems belong to the family of topological crystalline insulators [42][43][44][45][46] , i.e., gapped quantum systems characterized by crystal symmetries; they are characterized by quantized multipole moments in the bulk and unusual topologically protected states (e.g., corner or hinge modes) on their boundaries; see refs. [47][48][49][50][51][52][53][54][55][56][57][58][59][60] . Considering the prototypical Benalcazar-Bernevig-Hughes (BBH) model 47,48 , we unveil a phenomenon by which multiple BOs take the form of an oscillating Hall drift, accompanied with a synchronized inter-band beating, for special directions of the applied force, as summarized in Fig. 1a. Although the Hall motion is attributed to the finite non-Abelian Berry curvature of the degenerate band structure, the inter-band beating captured by the Wilson loop is shown to be topologically protected by winding numbers. The synchronization of real-space motion and inter-band dynamics is elucidated through a quantum Rabi oscillation of Wannier centers. Finally, we observe that detached helical edge states are present on specific boundaries, compatible with the special symmetry axes associated with the topological BOs. A topological transition signaled by the sign change of the identified winding numbers and by the corresponding appearence/disapperance of these states is identified.
Overall, our results demonstrate the rich interplay of non-Abelian gauge structures and winding numbers in the topological BOs of HOTI's, but also establish BOs as a powerful probe for non-Abelian topological properties in quantum matter.

Results
Model and symmetries. We consider the BBH model, as introduced in refs. 47,48 . It consists of a square lattice with alternating hopping amplitudes J 1 and J 2 in the two spatial directions and a π flux per plaquette, as depicted in Fig. 1b. We have introduced the flux by Peierls phases on the vertical links but other conventions can be used without affecting the results of this work. The model is represented by a chiral-symmetric Hamiltonian of the form where the 4 × 4 Dirac matrices are written in the chiral basis Γ i = −σ 2 ⊗ σ i for i = 1, …, 3 and Γ 4 ¼ σ 1 I. This model has twofold . The eigenfunctions of the lowest two bands read u 1 The full expressions of the d i (k)'s for the isotropic BBH model is d 1 ðkÞ ¼ ðJ 1 À J 2 Þ sinðk y =2Þ, d 2 ðkÞ ¼ ÀðJ 1 þ J 2 Þ cosðk y =2Þ, d 3 ðkÞ ¼ ðJ 1 À J 2 Þ sinðk x =2Þ and d 4 ðkÞ ¼ ÀðJ 1 þ J 2 Þ cosðk x =2Þ and the corresponding energy dispersion is displayed in Fig. 1c. Here, we take the periodicity d = 2a = 1, where a is the lattice spacing. The ordering of the unit cell sites chosen to represent the model in Eq. (1) is indicated in Fig. 1b. Notice that the chosen basis takes into account the geometric shape of the unit cell. As a consequence, the Hamiltonian is not Bloch invariant, namely H (k + G) ≠ H(k), with G a reciprocal lattice vector.
The presence of time-reversal symmetryT, withT 2 ¼ 1 and chiral symmetryŜ, represented by Γ 0 ¼ σ 3 I, sets the model into the BDI class 1 . Moreover, several crystalline symmetries are h i of a gaussian wavepacket experiences a sign-changing Hall drift under the applied force F, while displaying a synchronized beating within two occupied bands. This synchronized effect, represented by the metronome, is topologically protected by the winding number w. b BBH model: a square lattice with π flux and staggered hopping amplitudes J 1 and J 2 . Vertical bonds with arrows correspond to a Peierls phase π in the hopping amplitudes. c Band structure of the model. Each band is twofold degenerate. d Brillouin zones and paths C and C exhibiting topological BOs.
The two mirror symmetries guarantee that inversion (C 2 ) is also a symmetry of the model, withĈ 2 ¼M xMy . Moreover, the presence of mirror and rotation symmetries allows us to define a pair of mirror symmetries with respect to the diagonal axes of the lattice,M xy 1 M yĈ4 andM x y ¼ ÀM xĈ4 . This is one of the central ingredients allowing for topologically protected BOs, as we will show below.
Owing to the two non-commuting mirror symmetriesM x and M y , the BBH model is a quadrupole insulator that has a quantized quadrupole moment in the bulk, vanishing bulk polarization and corner charges 47 . The non-commutation of the mirror symmetries also provides a non-vanishing non-Abelian Berry curvature Ω xy ðkÞ ¼ ∂ k x A y À ∂ k y A x À i½A x ; A y of the twofold degenerate lowest (or highest) bands 61,62 , where A is the non-Abelian Berry connection 32 . However, the total Chern number of the degenerate bands remains zero due to time-reversal. It is then possible to define Wannier functions jν α x;k y i and jν α y;k x i, with α = 1, 2 numbering the bands below the energy gap, which are eigenstates of the position operatorsPxP andPŷP projected onto the lowest two bands, respectively 47,48 . The non-commutation of the mirror symmetries (and therefore of the projected position operators) forces the use of hybrid Wannier functions, namely Wannier states that can only be maximally localized in one direction 63,64 . Furthermore, it provides a necessary condition to have gapped Wannier bands, namely Wannier centers that are displaced from each other at every momentum k x or k y . In ref. 47,48 , the Wannier gap has been exploited to define a winding of the Wannier states (nested Wilson loop) as a condition to have a quantized quadrupole moment in the bulk, which can be revealed from the Wannier-Stark spectrum 65 .
Winding numbers. We now prove that a nontrivial topological structure captured by novel winding numbers characterizes the BBH model along the diagonal paths of the BZ, C and C, which are shown in Fig. 1d. In order to emphasize the generality of these results, we hereby consider a generic Dirac-like model [Eq. (1)] without specifying the components of the d(k) vector. We assume that all previously discussed symmetries are satisfied with the additional constraint that each d i function only depends on one component of the momentum k, namely we assume that . Such constraint is satisfied by the BBH model. Mirror symmetries impose that d 1 and d 3 are odd functions whereas d 2 and d 4 are even. We then find that along C (i.e., for k = k x = k y ), the diagonal mirror symmetry represented by the operatorM xy requires d 3 (k) = d 1 (k) and d 4 (k) = d 2 (k), while along C (i.e., for k = k x = −k y ), the symmetry operatorM x y requires d 3 (k) = −d 1 (k) and d 4 (k) = d 2 (k). We then conclude that only two components of d are independent and we therefore define the vector dðkÞ ðd 1 ðkÞ; d 2 ðkÞÞ. After writing the Hamiltonian in its chiral where ε 12 = −ε 21 = 1, and where we used thed vector of the BBH model in the last step. We therefore conclude that the quantities w C and w C count how many times the vectord winds over the closed paths C and C, respectively. The quantized windings w C and w C are here protected by the crystalline symmetriesM x ,M y , andĈ 4 , as shown in Supplementary Note 1. These symmetries also imply that the two topological invariants are not independent. We point out that similar winding numbers have been introduced in chiralsymmetric one-dimensional topological superconductors 66 .
Finally, the sign change of the winding numbers at the gap closing point J 1 = J 2 signals a phase transition. We will show below that the transition corresponds to the appearance of detached helical edge states. Let us now focus on the BOs of the BBH model and the role played by the quantized winding numbers discussed above.
Topological BOs: band-population dynamics. We consider a wavepacket obtained as a superposition of the lowest two bands and centered at k, which we write as u k ðtÞ j i¼ η 1 ðtÞ u 1 Owing to the degeneracy of the states u 1;2 k , other parametrizations can be chosen. As shown in Supplementary Note 2, this gauge ambiguity can be removed by weakly breaking time-reversal symmetry, thus splitting the two states in energy. Under an applied homogeneous and constant force F, which makes the crystal momentum change linearly in time, _ k ¼ F, the bands occupation evolves according to 67 where the matrix elements of the Berry connection are defined as Here, the force is assumed to be weak enough so that transitions to upper bands are neglected.
We can formally solve Eq. (3) as ηðtÞ ¼ expðÀi where we have denoted as k i and k f the initial and final momenta of the BO, respectively. For a closed path C 0 with k f = k i + G, where G is a reciprocal lattice vector, the bands population dynamics is determined by the Wilson loop matrix A Á dkÞ. Importantly, the winding numbers w Cð CÞ that we have previously introduced appear in the Wilson loops defined along the diagonal paths C and C, as W Cð CÞ ¼ expðið2π=4Þw Cð CÞ σ 1ð3Þ Þ, with w Cð CÞ ¼ ±1. From this, we obtain that BOs require four loops in momentum space in order to map the wavefunction back to itself, namely ½W Cð CÞ 4 ¼ I. Notice that the degeneracy of the bands brings a trivial dynamical phase that does not influence the internal band-population dynamics.
According to the classification of topological BOs discussed in ref. 40 , rotational symmetriesĈ n can quantize BOs with a force applied orthogonal to the rotational symmetry axis. This is (partially) the case here, withĈ 4 providing period-four BOs. However,Ĉ 4 symmetry alone is not sufficient to quantize the BOs. Additional symmetries, namelyM x andM y , are required in order to have a protected winding number along the paths C and C, see Supplementary Note 1. In the general framework presented in ref. 40 , a "Wannier-Zak" relation is demonstrated when mirror symmetries commute. As a consequence, the Zak phase winding that appears in the Wilson loop has a one-to-one correspondence with the position of the Wannier centers. This is well described by independently evolving point charges within a classical picture. In our case, such a direct correspondence is not possible owing to the non-vanishing Berry curvature and we will see below that the physical consequences of this feature appear on the real-space motion of the wavepacket.
Topological BOs: real-space dynamics. Let us now consider the real-space motion of the wavepacket's center-of mass, which satisfies the following semiclassical equations 67 Here, Ω xy denotes the SU(2) Berry curvature, whose components are shown in Fig. 2a for the BBH model; they satisfy the following conditions: Ω 11 xy ¼ ÀΩ 22 xy , Re Ω 12 xy ¼ Re Ω 21 xy and Im Ω 12 xy ¼ ÀIm Ω 21 xy . One anticipates from the accumulation of Berry curvature near the M point of the BZ that the paths C and C may display nontrivial features also in the real-space dynamics and not only in the band population beating discussed above. As we shall explain in detail below, the wavepacket experiences a transverse Hall drift that changes sign after each BO, thus bringing the center-of mass position back to its initial point after two BOs. This behavior is synchronized and tightly connected with the band-population dynamics captured by the Wilson loop.
The twofold degeneracy of the bands allows us to parametrize the evolving state u k ðtÞ j ion the Bloch sphere as η 1 ðtÞ ¼ cos θðtÞ and η 2 ðtÞ ¼ sin θðtÞe iϕðtÞ . We can therefore rewrite the anomalous velocity as η y Ω xy η ¼ ðjη 1 j 2 À jη 2 j 2 ÞΩ 11 xy þ sin 2θ cos ϕ Re Ω 12 xy À sin 2θ sin ϕ Im Ω 12 xy : On the C path, the angle ϕ is a constant of motion, namely _ ϕ ¼ 0. This means that the Bloch vector is confined to a meridian of the Bloch sphere. Moreover, since Re  4)). Once the path has crossed the M point, the occupations are flipped but so is also the sign of the Berry curvature, thus the Hall displacement continues with the same sign until the wavepacket reaches the Γ point. As the bands occupations have exchanged, a second BO will experience an opposite Hall drift and bring back the wavepacket to its initial position, as shown in Fig. 2c, d. From this, we obtain that the real-space motion is a witness of the non-Abelian band dynamics. For θ(0) ≠ 0, the offdiagonal component of the Berry curvature also contributes, and it fully suppresses the Hall displacement when θ(0) = π/4 since the anomalous Hall velocity vanishes identically: the two bands are equally populated, no band exchange takes place and therefore the positive and negative deflections compensate each other.
On the C path, the Berry curvature has only off-diagonal components and the Hall dynamics is determined by the relative phase ϕ, whereas θ is a constant of motion. In this case, the populations of the two bands do not exchange over time but the relative phase does by an angle π. The transverse displacement as a function of θ(0) and ϕ(0) for the two paths is shown in Fig. 2e.
We have thus found that the center-of mass of the wavepacket displays a period-two BO instead of a period-four one. This difference with respect to the Wilson loop analysis occurs because after two BOs, the wavefunction has picked up an overall phase 2 × (π/2), which does not appear in observablesÔ , such as for the center-of mass position.
Atomic limit: single-plaquette dynamics. In order to elucidate the role of Wannier functions in the BOs analyzed here and the synchronization between real-space and band-population dynamics, we consider the instructive atomic limit with J 2 = 0, where we can study the time dynamics of a single plaquette. The lowest energy eigenstates read u 1 j i ¼ ð1=2; 1=2; 0; 1= ffiffi ffi 2 p Þ T and We construct the position operatorr ¼ P i ðr i À r 0 Þ r i j i r i h j by setting the spatial origin at the plaquette center. We obtain the matricesx=a ¼ diagð1=2; À1=2; À1=2; 1=2Þ andŷ=a ¼ diagð1=2; À1=2; 1=2; À1=2Þ. Let us now callP the projector operator on the states u 1 j i and u 2 j i, from which we can construct the projected position operatorsPxP x P ¼ P α;β¼1;2 u α j i u α jxju β u β andPŷP ŷ P ¼ P α;β¼1;2 u α j i u α jŷju β u β .
It follows that ½x P ;ŷ P ≠ 0, whereas fx P ;ŷ P g ¼ 0. The eigenfunctions of the projected position operators are the Wannier functions jν x;y i and the corresponding eigenvalues are the Wannier centers 48,61,68 , which read here ν x = ν y = ±a/4. In the presence of an external tilt (or electric field), the perturbative Hamiltonian governing the dynamics for small values of the force F readŝ The projected Hamiltonian reveals how the external force induces a quantum dynamics between the eigenstates of non-commuting position operators, in the form of a Rabi oscillation. This result is in sharp contrast with the classical point-charge picture introduced in ref. 40 , which is valid for commuting position operators. We can now diagonalizeĤ F and we find the spectrum . The Rabi period can be easily obtained as T R ¼ 2π=ðFa= ffiffi ffi 2 p Þ. This solution is general and it does not depend on the direction of the force. Besides, we can always rotate the coordinate system in order to have one axis parallel to the force and one axis orthogonal to it, r ∥ and r ⊥ , and reduce the Hamiltonian toĤ F ¼ F kr k P . Then, the corresponding time dynamics can be represented by the eigenstates ofr ? P , namely the Wannier functions obtained by diagonalizingr ? P . As a consequence, we observe a transverse dynamics compared with the direction of the applied force F, as shown in Fig. 3a, b. However, the Wannier centers dynamics is not directly connected to the BOs and its period does not have to be the same as the Rabi period of the Wannier centers. For example, let us consider a BO with F y = 0. The periodicity of the BO occurs at the discrete times T B = 2πn/dF x , for n 2 Z þ where d = 2a. There is no solution that satisfies T B = T R . However, if we take F x =±F y we find that n = 2 provides T B = T R . Therefore, a force oriented along the diagonal axes allows to synchronize the Wannier centers dynamics with the BOs, whereas the other directions yield out-of-sync oscillations that does not bring the wavepacket back to its initial position at integer multiples of the fundamental Bloch period.
Away from the atomic limit, we can still use Wannier functions as a complete basis to express the wavepacket. A direct calculation (see Fig. 3c) shows that along the paths C and C, the Wannier centers remain gapped and their spectrum flat, namely they are equispaced along the entire path. We interpret this fact as a witness that the Wannier centers can be thought as oscillators with the same oscillation frequency (i.e., displacement), as in the atomic limit represented by Eq. (6), thus keeping the same oscillatory motion while changing the momentum k C or k C . In conclusion, Wannier centers perform a quantum periodic Edge states. The quantized winding numbers w C and w C , which we have previously identified along the paths C and C, indicate that a topological transition takes place when J 1 = J 2 . Here, we show that an open system with edges along the diagonals x ± y displays detached helical edge states. From Fig. 4a, we notice that near the atomic limit J 2 → 0, the bulk has gapped states at energies E b $ ± ffiffi ffi 2 p J 1 . The edge displays disconnected single sites at energy E s~0 and trimers, with energies E t1~0 and E t2 ¼ ± ffiffi ffi 2 p J 1 . Thus, a pair of zero energy (E s and E t1 ) modes exists at the edge.
To compute the dispersion relation of the edge modes, it is convenient to consider a cylindrical geometry. In this case, we find that the edge modes become helical, see Fig. 4b. An example of such states is shown in Fig. 4c.

Discussion
In this work, we have shown a new type of multiple BOs that is connected to the quantum beating of Wannier centers and we have identified HOTIs as a model where this effect can be observed. By studying the BBH model, we have shown that the Wilson loop imposes period-four oscillations and the centerof-mass motion displays an anomalous Hall displacement over one period of oscillation. We have connected these features to the crystalline symmetries of the model and we have identified quantized winding numbers that protect the topological BOs.
Moreover, we have shown that detached helical edge states emerge in an open system with the required symmetries.
Our results can be observed with cold atoms 70,71 , where flux engineering can be achieved through time-dependent protocols 72,73 and where the staggered hopping amplitudes requires a bipartite lattice 27,74 . Interferometric and tomographic methods can be exploited to measure the Wilson loop winding 25,31,75 and real-space cloud imaging makes possible to measure the center-of-mass displacement 76 . A fundamental question concerns the preparation of the initial state, owing to the degenerate nature of the bands. As shown in Supplementary Note 2 the bands can be split by slightly breaking time-reversal symmetry. In this case, it is possible to prepare a non-degenerate Bose-Einstein condensate (BEC) at the Γ point. When projected onto the eigenstates of the BBH model, this state is peaked at specific values of θ and ϕ. One can then obtain the desired superposition of the two zero-momentum modes (the BEC and the gapped mode) by a coherent coupling through an external driving. The subsequent BOs require that the applied force has a magnitude that is larger than the band separation to effectively recover the band degeneracy during the BOs.
In the context of photonics, our results can be investigated by using optical waveguides 77 , where it has been recently possible to realize synthetic π flux 78,79 . In this platform, the input laser profile can be inprinted in order to map the degenerate manifold of states at the Γ point that are parametrized by the angles θ and ϕ. It is then possible to reconstruct the Wilson loop dynamics by measuring the output field phase profile, whereas the Hall displacement is obtained from the spatial profile of the field intensity.
As a perspective of our work, it would be interesting to generalize our results to other two-and three-dimensional topological crystalline insulators and consider corrections to the semiclassical equations, e.g., involving the quantum metric once an inhomogeneous electric field or a harmonic trap potential are introduced 80,81 . Finally, given the role played by the initial state in the observation of the anomalous Hall displacement, BOs can be thought as a tool to witness the phenomenology of symmetry-broken condensates where the ground state degeneracy has been removed by interactions (Di Liberto et al., in prepration).

Methods
The SU(2) Berry curvature, defined as Ω xy ðkÞ ¼ ∂ k x A y À ∂ k y A x À i½A x ; A y , reads Real-space wavepacket dynamics. To validate the semiclassical real-space dynamics, we have numerically simulated the evolution of a real-space wavepacket using a finite size L x × L y lattice. We start by constructing a Gaussian wavepacket centered at k 0 = Γ = (0, 0) of the form ψ r;μ ðt ¼ 0Þ where μ indicates the sublattice degree of freedom in the unit cell and r μ is its spatial position. In the simulations we take a grid of k pts × k pts points in k space within an interval k ∈ [−3σ, 3σ] × [−3σ, 3σ]. We evolve the state with the realspace HamiltonianĤ and calculate the observable r num ðtÞ ¼ ψðtÞjrjψðtÞ h i , witĥ r P r ðr À r 0 Þ r j i r h j. We also study the exact evolution of the Bloch wave vector, by considering the velocity operatorvðkÞ ¼ ∂ĤðkÞ=∂k This method corresponds to solving the Schrödinger equation in k space. We find the displacement by integration r exact ðtÞ ¼ R t 0 dt vðtÞ þ r 0 .
Edge states. Here we derive the effective theory at the edge by considering periodic boundary conditions along the y 0 direction (see Fig. 5) for J 1 ≈ J 2 . In this stripe geometry, we have to double the unit cell to correctly represent the lattice periodicity, which reads d 0 ¼ 2 ffiffi ffi 2 p a. In chiral form, the Hamiltonian reads where Q s ðkÞ ¼ ÀJ 1 ÀJ 1 ÀJ 2 ÀJ 2 e Àik x J 1 ÀJ 1 J 2 e Àik y ÀJ 2 e Àiðk x þk y Þ ÀJ 2 e iðk x þk y Þ ÀJ 2 e ik y ÀJ 1 ÀJ 1 Here, we are taking units d 0 ¼ 1 and we are also using the convention that the lattice points within the unit cell are all sitting in the center of the unit cell. The Hamiltonian H s (k) is therefore in Bloch form. In order to build an effective theory near zero energy 69 , let us take J 2 = (1 + m)J 1 , with |m| ≪ 1. We can then split the Hamiltonian intoĤ s ðkÞ ¼Ĥ s ðk ¼ 0Þ þV s ðkÞ, whereV s ðkÞ is expanded to lowest order in k. The zeroth order termĤ s ðk ¼ 0Þ can be diagonalized and we find four eigenvectors v i 0 with energy E ¼ ± ffiffi ffi 2 p mJ 1 , which we use as a basis for the effective theory, and four high-energy states v i e that we neglect. We can then construct the projection operatorP s ¼ where we have rearranged the order of the components to have the Hamiltonian in block-diagonal form and we have defined in units where J 1 = 1.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The code that supports the plots within this paper are available from the corresponding author upon reasonable request.