Floquet chiral hinge modes and their interplay with Weyl physics in a three-dimensional lattice

We demonstrate that a three dimensional time-periodically driven (Floquet) lattice can exhibit chiral hinge states and describe their interplay with Weyl physics. A peculiar type of the hinge states are enforced by the repeated boundary reflections with lateral Goos-H\"{a}nchen like shifts occurring at the second-order boundaries of our system. Such chiral hinge modes coexist in a wide range of parameters regimes with Fermi arc surface states connecting a pair of Weyl points in a two-band model. We find numerically that these modes still preserve their locality along the hinge and their chiral nature in the presence of local defects and other parameter changes. We trace the robustness of such chiral hinge modes to special band structure unique in a Floquet system allowing all the eigenstates to be localized in quasi-one-dimensional regions parallel to each other when open hinge boundaries are introduced. The implementation of a model featuring both the second-order Floquet skin effect and the Weyl physics is straightforward with ultracold atoms in optical superlattices.

In this paper, we show that time-periodic driving of a three dimensional (3D) lattice can give rise to a new type of chiral hinge (i.e.second-order boundary) states.Such hinge modes are associated with a reorganization of the Floquet eigenstates in response to shifting from periodic to open boundary conditions.The phenomena may be intuitive understood starting from a fine-tuned point, where the quasi-energy bands have linear dispersion along only one direction and remain flat along the other two.When open hinge boundaries are introduced, this leads to reorganization of eigenstates of the system into a set states of the quasi-one-dimensional nature orienting parallel to each other and located at various distances from the hinge.Among them, the modes closest to certain lattice hinge correspond to trajectories with uni-directional (chiral) motion within each driving period.Importantly, the counterpart of such a chiral hinge mode with opposite chirality resides on the other side of the lattice.Therefore the chiral transport carried by these hinge modes would be preserved under a small perturbation assuming it does not mix eigenstates spatially separating far away from each other.Such a robustness is verified numerically by observing the particle transport encountering perturbations in the form of parameter change and local defects.
Surface modes exist widely under various circumstances, and it is helpful to mention what may be different in our system.First, unlike the case of (Floquet) higher-order topological phases [25][26][27][28][29][30][31][32][33][34], the hinge states studied here are formed inside the bulk spectrum near zero quasi-energy, as can be seen in Fig. 2 column (2).Their existence and robustness, therefore, requires a new theoretical explanation without invoking a bulk gap.Secondly, the localization of modes at the boundaries resembles the non-Hermitian skin effect [35][36][37][38][39][40][41], with a notable difference being the reduced dimensionality which may be thought of as a higher-order skin effect.Thirdly, different from one-dimensional (1D) periodically driven lattices [42] exhibiting helical modes, in our system the modes with opposite chirality residing at opposite hinges are well spatially separated, so the scattering due to a local perturbation would not reverse the chiral direction of the particle transport.
The aforementioned hinge states coexist with Weyl physics [43][44][45][46][47][48][49][50][51][52] when parameters are tuned away from the limit of uni-directional dispersionless motion.In that case a pair of Weyl points is created at quasienergy π leading to a Fermi-arc surface (i.e.first-order boundary) states.We find numerically that the parameter regime of the hinge states overlaps considerably with that for the Weyl physics, signaling a simultaneous observation of both phenomena in a unified experimental setting.
The model system proposed and studied here consists of a simple, stepwise modulation of tunneling matrix elements involving six steps.It generalizes to three dimensions a twodimensional lattice model introduced by Rudner et.al. [2] for studying anomalous Floquet topolgical insulators (see also Ref. [1]).The latter 2D tunnel modulation has been recently applied to ultracold atoms for the realization of anomalous topological band structures [11].Our 3D model can equally arXiv:2101.08281v3[cond-mat.mes-hall]11 Oct 2021 be implemented with ultracold atoms using such a stepwise tunnel modulation, now in 3D optical superlattices.Furthermore, besides giving rise to a new phenomenon, the unconventional chiral second-order Floquet skin effect, the model proposed here also provides a simple recipe for the implementation of Weyl physics by means of time-periodic driving.This paper is structured as follows.In the next Section II a 3D periodically driven lattice is defined.Subsequently the characteristic features of the bulk and hinge physics are considered in Secs.III and IV.The experimental implementation of our model, using ultracold atoms in modulated superlattices, is discussed in Sec.V, followed by the Concluding Section VI.Additional technical details are presented in four Appendices A -D.

II. MODEL
Consider a bipartite cubic lattice with alternating A-B sublattices in all three Cartesian directions.The lattice is described by a time-periodic Hamiltonian H(t + T ) = H(t), with the driving period T divided into 6 steps.In each step tunneling matrix elements −J between sites r A of sublattice A and neighboring sites r A ± ae µ of sublattice B are switched on, with µ = x, y, z.During the driving period T the tunneling steps appear in a sequence µ± = x+, y+, z+, x−, y−, z−, as illustrated in Fig. 1(a) [53].Within each step the evolution is determined by a coordinate-space Hamiltonian , where −J is the tunneling matrix element, r A specifies the location of sublattices A, and a is the lattice spacing such that r A ± ae µ denotes the locations of sites in sublattice B neighboring to the sublattice A site r A .The tunneling processes occurring in each of the driving steps are characterized by a single dimensionless parameter, the phase The one-cycle evolution operator (or Floquet operator), whose repeated application describes the time-evolution in stroboscopic steps of the driving period T , can be decomposed into terms corresponding to the six driving stages, When dealing with the bulk dynamics we impose periodic boundary conditions in all three spatial directions.The evolution operators for the individual driving steps can then be represented as: where τ 1,2,3 are Pauli matrices associated with the sublattice states A (corresponds to spin-up state) and B (corresponds to spin-down state) and where k µ with µ = x, y, z denotes the Cartesian components of the quasimomentum vector k.Here and in the following, we will use a dimensionless description, where time, energy, length and quasimomentum are given in units of T, /T , a, and /a, respectively.The quasienergies E n,k and the Floquet modes |n, k are defined via the eigenvalue equation for the Floquet operator We first note that the only global symmetry satisfied by the Floquet operator (3)-( 4) is a particle-hole flip Γ = CK, where Ki = −iK is complex conjugation and C = τ 3 the third Pauli matrix.Thus, the system belongs to class D in Altland-Zirnbauer notation [54,55].The Floquet operator satisfies [4,56], and therefore the quasienergies must appear in pairs E 1,k = −E 2,−k .Meanwhile, the system obeys the inversion symmetry PU F (k)P −1 = U F (−k), with P = τ 1 , enforcing that for each band one has E n,k = E n,−k .Together, we see that the Floquet spectrum has pairs of states with quasi-energies E 1,k = −E 2,k .This means possible gaps or nodal points/lines can, modulo 2π, appear only at quasienergy 0 or π.
At k x = k y = k z = ± π 2 (mod π) Eqs. ( 3) and ( 4) yield U F = 1, so that E n,k = 0. Therefore the quasienergy spectrum is always gapless at quasienergy 0 (modulo 2π), regardless of the driving strength φ.Then, the two-band spectrum could only possibly open up a gap at quasienergy equal to π (modulo 2π).To draw a complete phase diagram, we first note that flipping the sign of φ amounts to a particle-hole transformation U F | −φ = CU F | φ C −1 , and from the previous analysis we see that such a flip does not change the spectrum.Furthermore, from e −iφn•τ = τ 0 cos φ − in • τ sin φ, where τ 0 is the identity matrix, the periodicity of the Floquet operator with respect to the parameter φ is clearly seen: In this way, the irreducible parameter range is φ ∈ [0, π/2] as illustrated in Fig. 2.

III. PHASE DIAGRAM FOR PERIODIC BOUNDARY CONDITIONS
Before investigating the system with open boundary conditions and the emergence of chiral hinge modes, let us first consider the case of a translation invariant system with periodic boundary conditions.Here, we would focus on the better understood Weyl physics as an anchor point to draw the phase diagram, in order to use it as a reference frame to discuss the new hinge states in the next Section.We clarify that the phase diagram for Weyl physics overlaps, but does not coincide completely with that for the hinge states.
Let us begin with the topologically trivial high-frequency (weak driving) limit corresponding to φ 1.In that case one can retain only the lowest order terms in φ when expanding the stroboscopic evolution operator U F of Eq. (3), resulting in The spectrum ±2φ µ=x,y,z cos k µ corresponds to that of a static simple cubic lattice artificially described by 2 sublattices, where the bands are folded as the Brillouin zone size is halved.
While the system remains gapless at quasienergy 0 for arbitrary φ, a characteristic feature of the high-frequency (weak driving) regime is a finite energy gap at quasienergy π, resulting from the fact that the band width is proportional to φ and thus is small compared to the dimensionless driving energy ω = 2π.This behavior can be observed in the spectrum for φ = π/8 shown in Fig. 2 at the bottom of column (1).
Increasing the driving strength φ, the band width grows relative to 2π.When φ = π/6 the Floquet band, which is gapless at quasienergy 0, starts to touch itself at quasienergy π and momentum k = 0, as one can see in Fig. 3 (d).Going to the regime φ ∈ (π/6, π/2), the band touching point splits into a pair of Weyl points forming at quasienergy π with topological charges ±1, as shown in Fig. 3 (a).They are located at the quasimomenta k = k 0 d along the diagonal vector d = (1, −1, 1), (7) with (see Appendix A), so that k 0 → ±π/3 as φ → π/2.We observe the emergence of surface Fermi arc states connecting the Weyl points, when comparing the spectrum with full periodic boundary conditions to that with open boundary conditions along z-direction, as illustrated in Fig. 3 (b) and (c) respectively.Note that the reversed process of increasing driving frequency (or equivalently, decreasing φ towards π/6) would correspond to the usual scenario that two Weyl points merge together at k = 0 and the spectrum becomes gapped out.
As the driving strength approaches the fine tuned point, φ = π/2 − ε with ε 1, the Weyl dispersion acquires a highly anisotropic form shown in Fig. 4 (b).The dispersion remains steep along the diagonal coordinate k • d = k x − k y + k z , but becomes increasingly flat in other two directions.Exactly at fine tuning, φ = π/2, the constituent evolution operators (4) reduce to U µ± (k) = −i(cos k µ τ 1 ± sin k µ τ 2 ), and the Floquet stroboscopic operator takes the form U F = −e −i2τ 3 k•d .This provides the quasi-energies where m ∈ Z labels the Floquet bands, and where the upper and lower branches labeled by ± now directly correspond to sublattices A and B, i.e. ± → τ 3 .At this fine-tuned point the Weyl points disappear and the dispersion E k,± is completely flat for the momentum plane normal to d, as one can see in Fig. 4 (a).Hence a particle can only propagate along the diagonal d with a dimensionless velocity v ± = ±2d, depending on whether the particle occupies a site on the sublattice A or B at the beginning of a driving cycle.Note that for fine tuned driving (φ = π/2), the Floquet Hamiltonian is periodic in momentum space only thanks to the periodicity in quasi-energies.Such an effective Hamiltonian is characteristic for periodically driven systems, and does not have a static counterpart with finite range hopping.
The fine-tuned dispersion can be understood by considering the dynamics in real space.For φ = π/2 in each step µ± the particle is fully transferred from a sublattice A site positioned at r A to a neighboring site B situated at r B = r A ± e µ or vice versa.During the six steps composing the driving period, the particle follows the quasi-one-dimensional trajectory shown in Fig. 1 (b).Thus, after completing each period the particle located on a site of sublattice A (B) is transferred by +2d (−2d) to an equivalent site of the same sublattice, giving rise to stroboscopic motion along the diagonal directions ±d at the velocity v ± .
Before leaving this Section, we remark that the disappearing of Weyl points at φ = π/2 is not because they are gapped out (unlike at φ = π/6), but due to the accidental band-gap closure in a whole plane with Weyl points residing on it, as can be seen in Fig. 4 at φ = π/2.Consequently there is no need for the two Weyl points to merge at the same location in ).This can be inferred also from the dot size reflecting the inverse participation ratio with respect to the site basis as a measure for localization, as well as from the color code indicating the mean distance to two opposite hinges and interpolating from blue for one hinge over black near the center of the system to red for the other hinge.The green dots mark the Weyl points.We also plot the real-space densities of various Floquet modes for open boundary conditions along x and y (column 2) or all (column 3) directions.
Brillouin zone when they disappear in this way.More specifically, starting from π/6 < φ < π/2 where Weyl points are well defined, one can enclose a Weyl point with a gapped surface (eigenstates from a certain band) surrounding it, where the Weyl charge is the total Berry flux penetrating out of the surface.At exactly φ = π/2, any such surrounding surface would encounter a divergence of local Berry curvature as the two bands become gapless in a whole plane containing the Weyl point and thus are intersecting at any enclosing surfaces.Weyl points and their charges therefore lose definitions at such a fine-tuned point.When the phase φ exceeds π/2, the two Weyl points reappear with opposite charges compared to the φ < π/2 situation.

IV. CHIRAL HINGE MODES FOR OPEN BOUNDARY CONDITIONS
An intriguing effect shows up when the system is subjected to open boundary conditions with at least two properly chosen boundary planes, where a special type of chiral modes form at certain hinges.In that case eigenstates of the system reorganize into a set states of the quasi-one-dimensional nature orienting parallel to each other and located at various distances from the hinge when open boundaries are introduced.To gain an intuition, it will be useful to start from the finetuned point φ = π/2 where the stroboscopic motion of the particle at the hinge follows a classical trajectory in real space.Subsequently we will demonstrate numerically the robustness For full open boundaries the atomic motion is restricted by planes in all three Cartesian directions.We will also investigate a semi-infinite system restricted by two planes orthogonal to x and y, while keeping periodic boundary conditions along the z direction.Note that the three-fold symmetry x → y, y → z, z → −x would correspond to shifting the time origin by T/6 for the Floquet operator (3), which is merely a gauge transformation and does not change physical observables.Therefore, choosing another combination of open/periodic boundary conditions along different Cartesian directions results in identical phenomena.
Finally, since the boundaries are cut along the orthogonal directions for a cubic lattice, which is not terminating along the directions of Bravais vectors, it is more convenient to label the lattice sites in both sublattices directly by their Cartesian coordinates r = (x, y, z), with x, y, z taking integer values between 1 and L carrying units of the lattice constant.The sites belong to the sublattice A (or B) for even (or odd) values of x + y + z.
A. Fine-tuned driving

Real-space picture
We will commence our discussion of the hinge dynamics at the fine tuned point using the real space picture.Subsequently in Sec.IV A 2 we will present a rigorous analysis of the hinge states by going to the momentum representation along the hinge direction.Let us consider the atomic motion for x ≥ 1 or y ≥ 1 confined by a single boundary plane at x = 1 or y = 1, as illustrated in Fig. 1 (c).Starting from a site deep inside the bulk of, say, sublattice B, the particle will propagate stroboscopically on this sublattice covering a distance −2d in diagonal direction in each driving period, until approaching the boundary x = 1 or y = 1.At the boundary tunneling between the lattice sites B and A cannot occur during one of the six steps of the driving cycle.As a result, the particle changes the sublattice and starts to propagate on the A sublattice in opposite direction.The two colors on each panel of Fig. 1 (c) denote two possible types of such a reflection that are distinguished by the lattice site at which the driving cycle starts leading to the reflection.Specifically a tunneling event is impeded in the first (dark grey) or the fifth (yellow) driving step for a particle situated right at the surface plane x = 1 (y = 1) or next to it x = 2 (y = 2), as marked by small planes in Fig. 1 (c).
Let us next consider a similar motion of the particle in the region x ≥ 1 and y ≥ 1 confined by two intersecting surface planes x = 1 and y = 1.After the particle changes the sublattice B → A at the x = 1 boundary, it travels in reversed direction in steps of 2d until it eventually reaches the y = 1 surface and is again backreflected, this time with the sublattice change A → B. In this way, the particle will move back and forth between the x = 1 and the y = 1 surfaces which share a hinge.
It is now important to note that each reflection is accompanied by a lateral shift along the surface planes, as can be seen in Fig. 1 (c).Such lateral shifts during reflections resemble the Goos-Hänchen effect occurring when light is reflected at the boundary between two media.The lateral shift would have stronger influence when the particle start from a location closer to the hinge x = y = 1 or x = y = L where reflections occur more frequently, as illustrated in Fig. 1 (d).In particular, for particles starting from sublattice B at x = y = 1, the dynamics is completely given by the lateral shift along −z direction, as one can see in the lower part of Fig. 1 (d).The corresponding eigenstate for Floquet operator would be a single chiral mode localized along the hinge propagating along −z.Similarly, at x = y = L, one can verify that another chiral mode associated with sublattice A propagates along +z direction.Note that these two modes are separated by the whole system size, unlike in the 1D case [42] with two modes of opposite chirality residing at the same location and thus are subjected to backscattering on imperfections.
In the next Sec.IV A 2 we will see through rigorous calculation that the whole xy-plane can be divided into two parts by an off-diagonal line joining x = 1, y = L and x = L, y = 1, where eigenstates localized closer to x = y = 1 or x = y = L hinges will carry the chiral motion along −z or +z respectively.Therefore, under perturbations, it will be the eigenstates closest to the center "chirality border line" to lose their chirality along z first, and the hybridization gradually spreads towards the hinges as perturbation strengthens.Importantly, the eigenstates right at the hinges x = y = 1 and x = y = L do not involve any motion in the xy-plane, as can be seen in Fig. 1 (d).They correspond to the genuine chiral hinge modes described by the evolution operator U F over a single period.That means the chirality of such a mode does not depend on the ballistic trajectories finishing several Floquet periods in the xy-plane, and thus the closest hinge mode should be more robust against imperfections.This is verified numerically away from finetuned point in Fig. 2 (2nd column), and also in Figs.7 and  8 in the presence of local defects, as will be discussed more specifically later.
Prior to going to a rigorous analytical calculation, we make a side remark regarding the previous Weyl physics.The general boundary reflections exchanging sublattices A and B break particle-hole symmetry Γ = τ 3 K.But that will not destroy the Weyl physics because it reduces the symmetry class from D to A, which still hosts a Z classification for point defect (i.e.Weyl points) in a 3D Brillouin zone [55].Furthermore, the orthogonal boundary we choose still preserves the inversion symmetry, and therefore Weyl points and their associated Fermi arcs in the open boundary systems would still exhibit an inversion symmetric fashion as found previously.The only thing that might be slightly affected is the symmetry E 1,k = −E 2,k for Fermi arc on the surface, while the Weyl physics phase diagram concerns bulk spectrum and as L → ∞ surface effects would minimize.The hinge states discussed in this Section do not rely on the particle-hole symmetry either.

Floquet states
Before carrying out analytical calculations of the Floquet states at the fine-tuned φ = π/2, in Figs. 5 and 6 we illustrate the full path of a particle projected onto the xy-plane for reflections occurring at x = 1 and y = 1 surfaces.One can see in these figures that the particle returns to the same transverse position (x, y) after making two loops.Such a closed trajectory contains two reflections at each boundary x = 1 and y = 1, corresponding to the two different types of reflections depicted in Fig. 1(b).The back reflected particle propagates in a trajectory situated closer to the hinge or further away from it for the B → A or A → B reflections, respectively.This is similar to changing a track for a train before sending it backwards.Because of such chiral Goos-Hänchen shifts, the particle visits a larger number of B sites than A sites when traveling between the two surface planes.Recalling that the ballistic motion over the B or A sites takes place along the diagonal −d or d, the dynamics projected onto the xy place is accompanied by motion in −z or +z direction depending on the sublattice.Thus one arrives at an overall steady advance in −z direction after the particle completes a closed two loop trajectory in the xy plane containing more sites of the B sublattice.Therefore all eigenstates would occupy localized quasi-1D regions in xy-plane parallel to each other, each carrying certain group velocity along z.
More precisely, as demonstrated in Appendix B, after four reflections by the surfaces, a particle comes back to the original point in the xy-plane but shifts by −2 units in z-direction to an equivalent lattice site.One can similarly verify that a particle comes back to the same (x, y) transverse point after four reflections by the opposite x = L and y = L planes accompanied by a shift by +2 lattice unites along z direction.
In this way, the particle's trajectory starting from any lattice site would roughly covers a quasi-1D region in the xy-plane along e x − e y directions, whose length is proportional to the distance away from the hinge x = y = 1.All complete double loop trajectories in the xy-plane would be associated with a shift by −2 or +2 along z, depending on whether they are closer to the x = y = 1 or x = y = L hinge.Although such a picture applies to a particle situated further away from the hinges, the advance in the −z or +z direction takes place also for trajectories situated very close to the hinge x = y = 1 or x = y = L where the particle is reflected simultaneously from both hinge planes, as illustrated in Fig. 1(d The fact that trajectories represent closed loops in the xyplane allows us to construct analytically the eigenstates of the stroboscopic evolution operator U F with open boundaries in x, y and periodic boundary along z.The details are elaborated in Appendix B while we outline the results and derivations below.We would chiefly describe the hinge related to x = y = 1, while mention certain results for x = y = L directly. Suppose the particle initially occupies a site of the sublattice B with transverse coordinates y = 1 and x = M + 1, so the particle is situated at the hinge surface x = 1 and is M ≥ 0 sites away from the other y = 1 hinge surface, as illustrated in Figs. 5 and 6 for M = 4 and M = 5.In that case, it takes 2M + 1 driving periods for the particle to come back to the initial position (M + 1, 1) in the xy plane, while having shifted by −2 lattice units in z direction, i.e.
After averaging over such a full cycle involving 2M + 1 driving periods, the particle travels with an M-dependent mean A basis vector |s, x, y, k z is characterized by a quasimomenta where k z is defined modulo π, as the lattice periodicity equals to 2 length units in z direction.Equation ( 11) yields The above results are for a particle starting from the hinge surface y = 1, x = M + 1, which traces out a closed loop for its classical trajectory visiting 2M + 1 sites at different stroboscopic moments.In a similar manner, if the particle starts at any intermediate site of the trajectory at stroboscopic moment pT , it will also finish the loop after 2M + 1 driving periods: with p = 0, 1, . . .2M.By superimposing the states |M, k z , p corresponding to all sites involved in the stroboscopic trajectory, one can construct the Floquet hinge states representing eigenstates of the Floquet evolution operator U F : where the index q = 0, 1, . . ., 2M labels these modes.The corresponding quasienergies for the eigenstates ( 16) situated closer to the x = y = 1 reads An analogous dispersion but with an opposite slope (opposite sign for the group velocity along z) can be obtained for the states formed closer to the opposite hinge x = y = L.The dispersion branches for both types of hinge modes reproduce the spectrum shown in the row 1 and column (2) of Fig. 2. The corresponding eigenstates are quasi-1D in xy-plane along e y − e x and are plane-waves along z.
Such a spectrum of the system with open boundary conditions in x and y directions looks completely different from the one for full periodic boundary conditions shown in column (1) of Fig. 2 or Fig. 4 (a), where in the latter case all the bulk modes have the same positive or negative dispersion slope (group velocity) v z = ±2.In contrast, for the beam geometry [column (1) of Fig. 2] the spectrum is now organized into Eq.( 17), representing different quasi-1D eigenstates M sites away from the hinge x = y = 1.Each eigenstate is characterized by a different group velocity decreasing with the distance M from the hinge, where "−" and "+" correspond to the states located relatively closer to x = y = 1 and x = y = L respectively.The red / blue colors in Fig. 2 indicate the mean distance of each mode from the two relevant hinges.The dark red (blue) mode associated with M = 0 is localized directly at the hinge x = y = 1 (x = y = L) and propagates at the largest velocity in negative (positive) z direction.Modes with a smaller slope have larger M and thus are further away from the particular hinge, as indicated by the color.The real-space density of four different hinge modes at φ = π/2 is illustrated in the real-space plot shown in row 1 and column (2) of Fig. 2. In addition, a measure for the degree of localization of a mode |ψ is the inverse participation ratio IPR = j | j|ψ | 4 , with real-space site states | j and the Floquet eigenstates |ψ .It is shown in the spectra of Fig. 2 via the dot size roughly indicating the inverse of the number of sites a mode is spread over.Thus modes localized near the hinges have larger IPR than those that are more spread over in xy-plane locating further away from the hinges, indicating that the latter modes with smaller chiral group velocities are less narrowly localized in xy-plane than the former faster modes located closer to the hinge.In this way, the group velocity v M± along the z direction with open boundaries along x, y is very different from the projected group velocity v z = v • e z of a full periodic system.This is due to the unidirectional bulk group velocity v = ±2d and can be intuitively understood as follows.Each eigenstate of the system with open boundary conditions involves 4 boundary reflections, as can be inferred from Figs. 5 and 6.Such an eigenstate is composed of two pairs of counterpropagating plane waves representing the eigenstates of the original periodic system propagating along the cubic diagonal ±d with the z projection of the group velocity v z = ±2.Consequently the original group velocity in the bulk gets neutralized, and the overall group velocity v M± along z is due to the combined Goos-Hänchen shift over 4 times of boundary reflections.As the Goos-Hänchen process occurs more frequently near the hinges x = y = 1 and x = y = L, the magnitude of group velocity v M,± = ±2/(2M + 1) is maximum at the hinges for M = 0, while close to the central part where M ∝ L, the group velocity v M± becomes rather small, and eventually vanishes in the thermodynamic limit v M± ∝ 1/L → 0. This is quite different from a usual crystal where modes closer to the center of the bulk should be more insensitive to any boundary effects, while in our case one arrives at another situation where v z = ±2 for the periodic boundary conditions while for open boundary conditions one has v M± → 0 for the modes near the bulk center.
The observation that the whole (Floquet) spectrum and the corresponding eigenstates of our system are completely reorganized, when switching from periodic to open boundary conditions, resembles the extensive accumulation of boundary modes featured in the non-Hermitian skin effect [35][36][37][38].Therefore, the formation of an extensive number of reconstructed chiral hinge modes in our Floquet system might be called chiral second-order Floquet skin effect, in analogy to the terminology used for non-Hermitian systems [38].A subtle difference is that the eigenstates of a unitary Floquet evolution operator are orthogonal to each other, unlike eigenstates of non-Hermitian Hamiltonians that that are generally nonorthogonal.It implies the counting rule that on average, there can be at most one Floquet eigenstate per lattice site, so an accumulation of states right at the boundary is not possible in a Floquet system.Note also that the non-Hermitian skin effect is associated with the exceptional points of the non-Hermitian Hamiltonian when the boundaries are introduced [36,37], whereas no exceptional points are formed for periodically driven systems described by the unitary Floquet evolution operators.More details on these issues are available in Appendix C.

B. Beyond fine-tuned driving
Although the above analysis is based on ballistic trajectories at the fine-tuned driving parameter φ = π/2, we expect the chirality of the hinge modes to be robust also against perturba- tions and tuning away from φ = π/2.This applies especially to the hinge states with larger chiral velocities, which are situated closer to the hinge and thus are spatially well separated from counter-propagating modes at the opposite hinge.The chiral hinge modes persist for a rather wide range of φ beyond φ = π/2.This can be observed from the example of φ = π/3 (33.3% detuning, half-way across the Weyl phase transition) displayed in Fig. 2 column (2) around k z = π/2 and E = 0. We can see that the hinge modes at smaller M from the hinge still preserve their chirality along the hinge, as expected previously.In turn, hinge modes with larger M that are closer to the sample center, are gradually mixed with modes of opposite chirality and lose chiral nature together with localization in the xy-plane when φ deviates away from φ = π/2.
We have also considered the exemplary eigenstates for a cube geometry with full open boundary conditions along all Cartesian axes x, y and z presented in column (3) of Fig. 2 showing that for φ = π/2 and φ = π/3 the chiral modes at certain hinges are joined to form a closed loop respecting inversion symmetry of the system [see also Fig. 7 (a)].The six hinges not participating in this closed loop do not carry hinge modes, since their two boundary planes are not connected along the diagonal direction d.Meanwhile, non-hinge modes representing the bulk dynamics all center along the cubic diagonal [column (3) of Fig. 2].
In the previous Sections we have made arguments that at the fine tuning, φ = π/2, the chiral modes locating right at the hinge x = y = 1 (or x = y = L) do not involve the ballistic motion in the xy-plane over many driving periods, and therefore their chiral transporting feature should be preserved even if they become hybridized with nearby modes carrying different or no chirality.We now further test such an expectation by checking the robustness of chiral hinge transports in the presence of local defects in Figs.7 and 8. Here, we simulate the dynamics of a particle in the presence of a defect for a system with open boundary conditions in all three directions representing a generalization of the defect-free situation shown in column (3) of Fig. 2. Figure 7 illustrates the dynamics of a particle initially located at a corner (x, y, z) = (1, 1, 16) of the system, where two orthogonally oriented transporting hinges meet, (a) for the fine tuned situation without a defect and (b) in the presence of a strong potential offset of ∆ = 3π on two neighboring hinge sites (indicated by a green tube) at (x, y, z) = (1, 1, 8), (1,1,9), respectively.The corresponding plots for non-fine-tuned driving with φ = 0.9(π/2) are presented in Fig. 8.We find that despite this strong defect the chiral nature of the hinge modes ensures that no backscattering occurs at the defect and the majority of the wave-packet continues to follow chiral trajectories along the hinges.In Figs. 7 (b) and 8 (b) we also plot representative eigenstates of the system with defects.These eigenstates remain localized at the hinge, with only a small distortion compared to the situation without defect shown in column (3) of Fig. 2. That means the localization properties of the particle carrying the chiral motion along the hinge are not just an ephemeral phenomena.Once the original particle distribution overlaps with such an eigenstate, certain portion of those particles will always be localized along the hinge undergoing constant chiral motion.This is further verified in Fig. 9 for evolution over a long time.It shows the time-evolved state after 100 driving cycles for the non-fine-tuned system (φ = 0.9×π/2) both without defect (a) and with defect (b).Very similar distribution is also found after even longer evolution, e.g. over 1000 driving cycles; the densities are, thus, representative for late-time states in the limit t → ∞.Note that such a process implies the particle have encountered the defects for many (or infinite) times when they loop over the connected hinges repeatedly, without losing their localization the hinge and being scattered away.The densities are representative for late-time states in the limit t → ∞ as similar distribution is found also after 1000 driving cycles.

V. EXPERIMENTAL REALIZATION WITH ULTRACOLD ATOMS IN OPTICAL LATTICES
A. Engineering of the driven lattice Above, we have shown that the proposed modulation of tunneling gives rise to a variety of phenomena, including the robust creation of a pair of Weyl points, unidirectional bulk transport, chiral Goos-Hänchen-like shifts, and the macroscopic accumulation of chiral hinge modes for open boundary conditions corresponding to a chiral second-order Floquet skin effect.The model itself is, nevertheless, rather simple and its implementation with ultracold atoms in optical lattices can be accomplished using standard experimental techniques.All what is needed is a static cubic host lattice potential of equal depth V 0 in each Cartesian direction and a superlattice potential, whose amplitudes along various diagonal lattice directions are modulated in a stepwise fashion in time in order to suppress/allow tunneling along the six different bonds specified by our protocol.This can be achieved using the following optical lattice potential: where only two of the four modulating lasers α ab with a, b = 0, 1 are turned on in each driving step, as shown in Fig. 10.Such a modulation provides the required six-stage driving of the cubic lattice.Note that a similar modulation has recently been implemented in two dimensions [11].

B. Detection of hinge dynamics
To observe the dynamics associated with the hinge modes, one can apply the boxed potential achieved in recent experiments [57][58][59][60][61][62][63].There, thin sheets of laser beams penetrate through the quantum gases creating a steep potential barrier.Three pairs of such beams are imposed in a three-dimensional system, creating the sharp "walls" for the box potential while leaving the central part of the gases homogeneous.
Essentially, such a potential combined with our lattice driving scheme immediately leads to the particle dynamics described in Fig. 7 and Fig. 8.To take into account realistic experimental situations, two modifications are adopted in our following simulations.First, we consider the effect of a relatively "softer" wall for the box potential with where the potential ramps up over a finite distance of roughly 4ξ near the boundaries r (1,2)  µ , see Fig. 11 for instance.The second modification we adopt is that the initial state is not taken to be localized on a single lattice site but described by a gaussian wave packet of finite width, The results of the dynamics are presented in Fig. 11 (c1)-( c3) and (d1)-(d3), for the ideal sharp boundary (as in Fig. 7) and the realistic softer boundaries in experimental setting respectively.We see that the chiral motion snapshots for the ideal/softer boundaries exhibit qualitatively the same characters, signaling that a softer boundary does not cause significant changes.This is in consistent with the previous simulation showing the robustness of hinge states and their resulting chiral dynamics against local defects in Fig. 7 and Fig. 8.The major difference from previous cases, then, derives from the initial state that overlaps with more than one set of eigenstates near the hinge, each with different group velocities as given in Eq. (18).Finally, we mention that some portion of initial particle distribution would reside within the region with significant changes in V box .That portion of the particles could be permanently confined to the initial hinge due to a mechanism similar to Wannier-Stark localization.However, the majority of the particles are still traveling into the connecting hinges, as shown in Fig. 11 (c3) and (d3).
In cold atom experiments, the density profiles are usually detected by taking a certain projection plane, where the integrated (column-averaged) densities are observed.To this end, we point out that the hinge dynamics can be confirmed by observing the density profiles in two perpendicular planes.A schematic plot is given in Fig. 11 (b), corresponding to the dynamics along the hinge x = y ∼ 1.The density profile taken at x − y plane (i.e. the "top" view) would show a localized distribution at the corner, verifying the particles only locate at x = y ∼ 1.Meanwhile, the profile at x − z plane (i.e."side" view) indicates the movement/spreading along z.In a more general situation, i.e. at long time limit with all 6 hinges populated as in Fig. 9, additional image projection planes could be exploited.We also mention that a simultaneous implementation of multiple imaging planes have been applied in experiments [49,50]

C. Detection of Floquet Weyl points
Weyl physics has been explored in recent cold atom experiments and theoretical proposals [46,49,50,64,65], and also extensively in solid state systems [43].Here, we discuss a scheme closely related to a recent experiment [11,66] detecting the spacetime singularities in anomalous Floquet insulators.
First, the band touching at Weyl points can be verified using the Stückelberg interferometry [11,[67][68][69].Such a method measures the smaller gap for the two bands at E ∼ 0 and π, see Ref. [11] for details.Compared with the experiments for insulators, a difference here is that the two bands are, overall, always gapless at E = 0.That means if a global gap is measured, it will prevent us from gaining information about the gaps or band touching at quasienergy E = π.But fortunately, there exists a finite region neighboring to k 0 = π/3(1, −1, 1) where the bands are always gapped at E = for all φ, see Fig. 12 (a1), (a2) for example.Then a local gap closure can be measured near k 0 .The specific measurement for our case can be performed in the following way.An example for k 0 = (π/3 − 0.2)(1, −1, 1) is presented in Fig. 12 (b).Let us denote the quasi-energy of the two Floquet bands at k 0 with E ± (k 0 ) = ±E 0 .Here, we use the branch cut along π in taking the logarithm of Floquet eigenvalues e −iE ± (k 0 ) .They have the same magnitudes and opposite signs due to particle-hole and inversion symmetry as explained in Sec.II.Then, the local gap around E ∼ 0 is ∆ (0) ≡ 2E 0 , while the other gap around the Floquet Brillouin zone boundary E ∼ π is ∆ (π) ≡ 2π−2E 0 .Therefore, ∆ (0) = ∆ (π)  can only occur at 0, π mod 2π.In experiments, one can start from the high frequency limit (φ → 0) where the band width is small compared to 2π and therefore the measured gap always corresponds to ∆ (0) .Slowing down the driving, the gap ∆ (π) shrinks while the other gap ∆ (0) expands.At some point, the two gaps coincide with their magnitudes, as shown in Fig. 12 (b).Since it is always the smaller one of ∆ (0) and ∆ (π) that will show up in experimental measurement, one will observe a cusp shape of the measured gap, i.e. near φ ≈ 0.41 in Fig. 12 (b).One could then imply from the occurrence of the cusp that for φ > 0.41, the experimental data starts to reveal ∆ (π) , whose vanishing at φ ≈ 0.87 shows the existence of the Weyl point around E ∼ π.Similar measurements can be performed for k slightly deviating from k 0 , which will show that at φ = 0.87, ∆ (π) remains finite, proving that the band closure around E ∼ π is a point contact.When the designated φ is slowly approached, one can perform a measurement of the gap at a certain k.A shortcut for our system is that focusing on quasimomenta along the diagonal k 0 = k 0 (1, −1, 1) is sufficient to determine the Weyl point, as discussed in Sec.III.
With the Weyl points determined, one could further apply band tomography [70,71] method for momentum states surrounding a certain Weyl point in order to determine its charge.Note that one does not need the eigenstate information throughout the whole Brillouin zone as the two bands are gapless in certain regions, except for just an arbitrarily small surface wrapping a Weyl point k (Weyl) determined previously.As shown before, near the Weyl points in our model, there exists a finite region where the two bands are fully gapped in both E ∼ 0 and π, which allows for populating eigenstates with bosons at a certain band [11,71].As an example, in Fig. 13 (a) we illustrate the surfaces formed by 6 faces q x,y,z = ±0.5 of a cube, where q = k − k (Weyl) , with k (Weyl) ≈ 0.955 × (1, −1, 1) for φ = π/3, as in Fig. 2. From the full information of the Floquet eigenstates |u n,k given by Eq. ( 5), the Berry curvature penetrating out of a plane normal to the unit vector e µ can be computed as , where ε µνρ is the Levi-Civita symbol, and ± sign denotes that the unit vector penetrating out of the cube is along ±e µ directions.Figure 13 shows momentum resolved Berry curvatures in each wrapping surface and their net fluxes surf d kΩ µ (k) in that plane.

VI. CONCLUSION
In this paper, we have shown that three-dimensional periodically driven lattice systems can experience a complete reconstruction of its eigenstates with drastically different features, when subjected to open boundary conditions.This corresponds to a chiral second-order Floquet skin effect.An intuitive understanding of this effect was given by considering the system at a fine-tuned point of the periodic driving, where the bulk motion can only take place forwards or backwards along a single diagonal direction.As a consequence, for open boundary conditions, the particle is reflected back and forth between hinge-sharing surface planes, indicating a combination of counter-propagating original bulk eigenstates carrying opposite group velocity into new chiral Floquet states associated with the hinge.Thus the original bulk motion becomes neutralized, while the new group velocities are associated with a Goos-Hänchen shift during boundary reflections, which are most significant close to hinges where the particles is reflected more frequently.
This effect offers a alternative mechanism to generate hinge modes in addition to the more well-understood cases such as (Floquet) higher-order topological phases.The robustness of the chiral hinge modes here is more relevant to the special dispersion characteristic of Floquet systems, which allows for all eigenstates to be localized in a quasi-one-dimensional re- gion close to fine-tuned point, while not necessarily require a bulk gap.While we briefly discuss in Appendix D a topological quantity specific to the fine-tuned point φ = π/2, it will be an interesting future direction to explore whether such chiral hinge modes would have a stable and rigorous topological nature.The previous analysis [72] have shown that the requirement of bulk quasi-energy gap becomes less important to characterize topologically protected chiral modes (i.e.Fig. 4 (a) (ii) in [72]).While weak disorder could mix different momentum states and destroy the hinge modes, it will be interesting to see whether a set of more robust hinge modes would emerge in the strongly disordered regime where all modes are localized, and to check the possible crossover or transitions when one increases the disorder strengths.Our current approach and discussions would also be useful for considering observable effects in experiments, similar to the case of Weyl semimetals where surface defects/disorder would partly destroy the Fermi arc but the chiral transport will be partially preserved therein [73].
It is noteworthy that even modes close to bulk center suffer from a complete reconstruction by open boundaries, hosting very different group velocity in the thermodynamic limit.Such an effect resembles the accumulation of boundary modes in systems described by a non-hermitian Hamiltonian.Yet different from this non-Hermitian skin effect, in our system we can have at most one state per lattice site on average, as Floquet states are orthogonal to each other.Thus in our system the accumulation corresponds to the emergence of hingebound modes at increasing distance from the hinge.Another interesting aspect is the competition or interplay between the hinge modes and the emergence of robust Wely points in our system, so the hinge states can co-exist with the Fermi arc surface states.The implementation of the model featuring both the second-order Floquet skin effect and the Weyl physics is straightforward with ultracold atoms in optical superlattices.all 2M + 1 states of the stroboscopic sequence featured in Eqs.

(B3)-(B9)
The origin of such chiral hinge states can be explained as follows.The particle in the sublattice B is reflected to a site of the A sublattice situated closer to the hinge, whereas the particle in the sublattice A is reflected to a site of the sublattice B situated further away from the hinge, as one can see in Figs. 5 and 6, as well as in Eqs.(B3), (B5), (B7), (B9).Consequently the number of B sites visited over all 2M + 1 driving periods (M + 1) exceeds the corresponding number of A sites (M).The four reflections do not yield any total shift of the particle in the z direction.On the other hand, the ballistic motion between sites the same sublattice B (A) is accompanied by a shift by 2 lattice sites in the z (−z) direction for each driving period.This leads to the overall shift of the particle to an equivalent site in the −z direction is due to the difference in the number of the visited B and A sites after 2M + 1 driving periods.
Appendix C: Non-Hermitian Hamiltonian corresponding to stroboscopic operator Recently it was suggested [74] to associate a non-Hermitian Hamiltonian H NH (k) to the momentum space stroboscopic evolution operator iU F (k).Let us consider such a non-Hermitian Hamiltonian for our 3D periodically driven lattice For the fine tuned driving (φ = π/2) the bulk stroboscopic evolution operator corresponds to a non-Hermitian Hamiltonian describing a unidirectional transfer between the lattice sites along the diagonal d = (1, −1, 1) and in the opposite direction −d for the sublattices A and B, respectively: (C3) The bulk non-Hermitian Hamiltonian (C2) supplied with the open boundary conditions (C3) describes a unidirectional coupling between unconnected linear chain of the A or B sites terminating at the hinge planes.The eigenstates of each linear chain represent non-Hermitian skin modes which are localized at one end of the chain depending on the direction of asymmetric hopping [36,37].In the present situation such skin modes would be trivially localised on different planes of the hinge for the chains comprising different sublattice sites A or B, and no chiral motion is obtained along the hinge.Yet the open boundary conditions (C3) are not sufficient to properly represent the boundary behavior of a particle in our periodically driven lattice.In fact, bulk non-Hermitian Hamiltonian (C2) supplied with the boundary conditions (C3) is no longer a unitary operator.Thus one can not associate such an non-Hermitian operator with the evolution operator, in contradiction with Eq. (C1).The unitarity is restored by adding to H bulk NH extra terms [in addition to the conditions (C3)] to include effects of the chiral backward reflection at the hinge planes to a neighboring linear chain composed of the sites of another sublattice.These terms are described by Eqs.(B3), (B5), (B7) and (B9), and correspond to the dashed lines in Figs.5-6.The chiral reflection appears due to the micromotion during the driving period involving six steps (see Fig. 1(b)), and is a characteristic feature of the periodically driven 3D lattice.As discussed in the previous Appendix B, the backreflection leads to the chiral motion along the hinge.

Appendix D: Tentative discussions on topological origins
An interesting question is whether the robust chiral hinge states are a consequence of topological properties of the driven system.However, as we see from column (2) of Fig. 2, the quasi-energies of hinge states are fully mixed with the bulk spectrum, and therefore no traditional topological band theories for gapped or semimetallic systems apply.Here, we apply an approach similar to that used in refs.[1] and [42] for computing topological invariant at a fine-tuned parameter point, where eigenstates are analytically obtainable, and a topological invariant based on eigenstate projections can be calculated.
In section IV A we have obtained equation ( 15) describing the evolution over 2M + 1 driving periods of the Mth hinge state |M, k z , p at the fine tuned point.Using this equation one can define the the quasienergy winding number [1] for the Mth hinge state via the Floquet evolution operator over the 2M + 1 driving periods: where with p = 0, 1, . . .2M.In Eq. (D1) the integration over k z extends over one Brillouin zone of width π, as the distance between two non-equivalent lattice sites equals to 2 in z direction.A similar procedure can be applied to the opposite hinge at (x, y) = (L, L), where the hinge modes shown in blue in column (2) of Fig. 2 are characterized by the opposite group velocity and thus the winding number is opposite W M = −1.
The rigorous quantization of the topological invariant W M is associated with the fine tuning, φ = π/2.However, the spatial separation between hinge modes of opposite chirality allows to preserve their chiral character also away from the fine tuning point φ = π/2, as one can see in column (2) of Fig. 2. Thus, the formation of bulk states via the mixing of hinge modes of opposite chirality happens mostly in the center of the system, where hinge modes corresponding to large M and small chiral velocities lie close by to their counter propagating partners associated with the opposite hinge.In turn, the states with the largest chiral velocity, which are situated close to the hinge and far away from counter-propagating modes of the opposite hinge, are much less affected by a small detuning.
In this way, tuning away from φ = π/2 we find a crossover (rather than a topological transition) in which the chiral hinge modes are gradually destroyed, as can be observed in the realspace plots in columns (2) and (3) of Fig. 2. Note that already a small deviation from the fine tuned point φ = π/2 destroys the chiral hinge states in a narrow region near k z = 0 and E = π, as can be seen from the spectrum shown for φ = π/3 in column (2) of Fig. 2. In this spectral surface Fermi arc states are formed, which equally provide definite chiral transport, yet around the surfaces rather than the hinges.Thus a fraction of the chiral hinge states is transformed into Fermi arc surface states in the vicinity of the Weyl points.The latter states extend to a larger and larger spectral area as the detuning increases.
FIG. 1.(a) Bonds connected during the driving steps 1 to 6.(b) Bulk trajectories within a Floquet cycle at the fine tuned point φ = π/2.Depending on the starting sublattice, particles will travel in opposite directions along the cubic diagonal d = (1, −1, 1) within a period.(c) Trajectories at the same φ = π/2 but with a surface termination (open boundary).Particles starting from sublattice A (or B) near the y = 1 (or x = 1) surface would have their dynamics impeded by the open boundary at a certain driving step.That results in a switch of sublattice after a Floquet cycle, and therefore the direction of the trajectory is reversed after the reflection from the surface.(d) The hinge formed by two intersecting terminating surfaces renders unidirectional modes at φ = π/2.The figure shows two of the modes closest to the hinge starting at a site of B (lower plot) or A (upper plot) sublattices directly at the hinge.Each color for the arrow denotes one Floquet cycle.The lower-panel of (d) shows the mode situated directly at the hinge which undergoes a uni-directional motion along z within a period.

φ ( 1 )FIG. 2 .
FIG. 2. Table showing energy spectra and Floquet modes.The phases of the system are indicated in the vertical axis on the left.The table's rows refer to the driving parameters φ = π/8, π/3, π/2, corresponding to the metallic phase, the Weyl semimetal/hinge phase, and the fine-tuned point, respectively.The columns represent periodic/open boundary conditions (PBC/OBC) along the specified directions.The system sizes are L x = L y = 40 in column (1), 20 in (2), and 16 in (3).Columns (1) and (2) contain energy spectra projected onto k z .Although the spectra for fully periodic and open boundary conditions in x and z direction are almost identical in the metallic phase (φ = π/8), they are very different in the WSM/hinge phase (φ = π/3) and at fine tuning (φ = π/2).The difference is a result of the formation of chiral hinge modes for open boundary conditions in x and z directions (column 2).This can be inferred also from the dot size reflecting the inverse participation ratio with respect to the site basis as a measure for localization, as well as from the color code indicating the mean distance to two opposite hinges and interpolating from blue for one hinge over black near the center of the system to red for the other hinge.The green dots mark the Weyl points.We also plot the real-space densities of various Floquet modes for open boundary conditions along x and y (column 2) or all (column 3) directions.
FIG. 3. (a) The different phases can be distinguished by the presence or absence of Weyl points at quasienergy π.For φ = π/8, in the metallic phase (φ < π/6), no Weyl points are present.At the transition, φ = π/6, the band touching point appears at k = 0 (green circle).For π/6 < φ < π/2 the band touching splits into two Weyl points of opposite charge (green circles with ± sign) that separate along the diagonal k x = −k y = k z , as shown for φ = π/3.At fine tuning φ = π/2 the dispersion becomes flat in two directions and the Weyl points disappear to reappear for φ > π/2 with reversed charges (signs in green dots correspond to the limit φ = π/2 − 0).(b-d) Quasienergy spectra versus k x , k y for periodic boundary conditions in x and y and either periodic (b,d) or open (c) boundary conditions in z direction.A surface Fermi arc can be observed for φ = π/3 by comparing the spectra with periodic (b) and open (c) boundary conditions the z direction.The orange surface denotes the contour formed by quasi-energies closest to E = π at each (k x , k y ).For φ = π/6 (d) a pair of Weyl points reduces to a single band touching point at the quasi-energy E = π.
FIG. 4. (a) Anisotropic 1D-like dispersion along the coordinate k x − k y + k z at the fine-tuned point φ = π/2.(b) For a small detuning, φ = π/2 − 0.1, the dispersion is no longer completely flat in other two directions, and a pair of non-equivalent Weyl points is formed along the diagonal at k = ±k 0 (1, −1, 1).Yet the dispersion remains highly anisotropic.

FIG. 5 .
FIG.5.An example of a fine-tuned stroboscopic motion of a particle at the lower-left hinge.The picture shows the projection of the particle's trajectory in the xy plane.The sites of the B and A sublattices are marked in blue and red, respectively.The particle is initially at the site of the sublattice B characterized by the coordinates x = M+1 and y = 1, with even M = 4.This corresponds to the lower row (y = 1) and the fifth column (x = 5).Subsequently the particle undergoes the stroboscopic evolution described by Eqs.(B3)-(B9) in Appendix B. Dashed lines with arrows show stroboscopic reflections from the planes x = 1 or y = 1.Bulk ballistic trajectories over one/two driving periods are indicated by thin/thick solid arrows.The particle returns to its initial site after 2M + 1 = 9 periods.

FIG. 6 .
FIG.6.Like in Fig.5the particle is initially at the site of the sublattice B with x = M + 1 and y = 1, but now with odd M = 5.The particle returns to the initial site after 2M + 1 = 11 periods.
velocity of v M− = −2/ (2M + 1) along the hinge axis z (see Appendix B for more details).Here |s, x, y, z describes a particle located at site (x, y, z) belonging to sublattice s = B = −1 or s = A = 1 with s = (−1) x+y+z .Let us now consider periodic boundary conditions in z direction with z = 1, 2, . . .2N z , i.e. |s, x, y, z + 2N z = |s, x, y, z , while keeping open boundary conditions in x and y.It is then convenient to introduce hybrid position-momentum basis states |s, x, y, k z = 1 √ N z z e ik z z |s, x, y, z .

20 FIG. 7 . 20 FIG. 8 .
FIG. 7. The dynamics of a particle initially localized at the site (x, y, z) = (1, 1, 16) for a system of 16 × 16 × 16 sites with full open boundary conditions and fine tuned φ = π/2.The squared wave function at different times is reflected in the opacity of the plotted dots.Different colors indicate time.(a) Without defect.(b) In the presence of a potential defect of energy ∆ = 3π at the two sites marked by the green tube.Additionally, we plot the squared wave function of one hinge mode.(c) Snapshots of the time evolution in the presence of the defect at different times.
FIG.9.Density distribution of a particle initially localized at the corner site (x, y, z) = (1, 1, 16) after an evolution over 100 driving cycles for the non-fine-tuned parameter φ = 0.9 × π/2, without defect (a) and with a defect (b) [corresponding to the parameters of Fig.8 (b)].The densities are representative for late-time states in the limit t → ∞ as similar distribution is found also after 1000 driving cycles.

FIG. 10 .
FIG.10.The stepwise modulation of the dimensionless superlattice amplitudes α ab according to the protocol given in the table, gives rise to different dimerizations of the cubic lattice in each driving step, that enables tunneling along the desired bonds. .
FIG. 11.The dynamics of particles in a box potential with different softness of the boundary.The initial density distribution takes a Gaussian profile spreading over several lattice sites.The parameters are φ = π/2, V b T/6 = 7.5π.The initial Gaussian profile has the center (x 0 , y 0 , z 0 ) = (4, 4, 12) and width s 0 = 0.75.
FIG.12.Simulation of gap measurements using Stückelberg interferometry.(a1) and (a2) Contours for quasi-energy gap at E = 0, for φ = π/8 and φ = π/3 respectively.(b) The gaps at E ∼ 0 and E ∼ π, and the measured gap which takes the smaller one of the two.
FIG. 13.The Berry curvatures for the surfaces wrapping a Weyl point.Adding the net Berry curvatures up we have 2π.
r A |A, r A + 2d A, r A | − i r B |B, r B − 2d, B, r B | .(C2)The open boundary conditions for the hinge corresponding to x ≥ 1 and y ≥ 1 are obtained by imposing a constraint on the state-vectors entering the real space non-Hermitian Hamiltonian (C2) |s, r s = 0 for r s • e x,y ≤ 0, with s = A, B .