Floquet Second-Order Topological Phases in Momentum Space

Higher-order topological phases (HOTPs) are characterized by symmetry-protected bound states at the corners or hinges of the system. In this work, we reveal a momentum-space counterpart of HOTPs in time-periodic driven systems, which are demonstrated in a two-dimensional extension of the quantum double-kicked rotor. The found Floquet HOTPs are protected by chiral symmetry and characterized by a pair of topological invariants, which could take arbitrarily large integer values with the increase of kicking strengths. These topological numbers are shown to be measurable from the chiral dynamics of wave packets. Under open boundary conditions, multiple quartets Floquet corner modes with zero and π quasienergies emerge in the system and coexist with delocalized bulk states at the same quasienergies, forming second-order Floquet topological bound states in the continuum. The number of these corner modes is further counted by the bulk topological invariants according to the relation of bulk-corner correspondence. Our findings thus extend the study of HOTPs to momentum-space lattices and further uncover the richness of HOTPs and corner-localized bound states in continuum in Floquet systems.

In this manuscript, we investigate a periodically kicked rotor in two spatial dimensions, whose momentum space could form a two-dimensional (2D) discrete lattice holding rich Floquet HOTPs. In Section 2, we introduce the Hamiltonian of the system and obtain its Floquet operator under the quantum resonance condition. Based on the symmetry analysis of the model, we construct a pair of integer topological invariants (w 0 , w π ) in Section 3, which could fully characterize the Floquet HOTPs that are protected by the chiral symmetry of the system. These Floquet HOTPs are further shown to be able to possess arbitrary large topological numbers with the increase of kicking strengths. In Section 4, we show that these topological invariants could be dynamically probed by measuring the time-averaged mean chiral displacements of wave packets in two-dimension. Under the open boundary conditions, we find many quartets of Floquet corner modes at zero and π quasienergies in Section 5. The numbers of these corner modes are predicted by the bulk topological invariants (w 0 , w π ), yielding the bulk-corner correspondence of Floquet HOTPs in momentum space. Moreover, the zero and π Floquet corner modes are found to be embedded in the continuous bulk bands of delocalized states, forming corner-localized Floquet bound states in the continuum that are originated from higher-order Floquet topology. We summarize our results and discuss potential future directions in Section 6.

Model
In this section, we introduce a representative driven lattice model, which could possess rich Floquet HOTPs in momentum space. Our system can be viewed as a twodimensional extension of the double kicked rotor (or lattice) model [109][110][111][112], which describes a quantum particle kicked twice by a periodic potential at different times within each driving period. The time-dependent Hamiltonian of the system takes the form Here, (x,ŷ) and (p x ,p y ) are the position and momentum operators of the particle along x and y directions. T is the driving period. τ ∈ (0, T) controls the time delay between the two kicks inside a driving period. κ 1,3 and κ 2,4 are kicking strengths of the potentials along x and y directions. φ x,y ∈ [0, 2π) describe the phase differences between the kicking potentials applied at t = T and t = T + τ in the 's driving period. The quantities in Equations (2)-(4) are all set in dimensionless units. In experiments, the model Hamiltonian H may be realized in cold atom systems, where the kicking potentials could be implemented by optical-lattice potentials with relative phase shifts [111,112]. Within a given driving period (e.g., from t = 0 − to t = T + 0 − ), the dynamics of the system is therefore governed by a kickV applied at t = 0, followed by the free evolutionĤ 0 from t = 0 → τ, a second kickŴ at t = τ, and finally the free evolutionĤ 0 over a time duration T − τ. The Floquet operator, which describes the evolution of the system over such a complete driving period, is then given byÛ Due to the periodicity of kicking potentialsV andŴ inx andŷ, the eigenvalues of momentum operatorsp x,y take the formsh(n x,y + k x,y ), where n x,y ∈ Z, and k x,y ∈ [0, 1) are the quasimomenta. A Floquet system with periodicity in momentum space is achieved by setting k x,y = 0, which maybe realized experimentally by a Bose-Einstein condensate with large coherence width [113][114][115][116]. Under this condition, we can identify the momentum operatorsp x,y in Equation (2) ashn x,y with integer eigenvalues n x,y . The Floquet operator in Equation (5) then takes the explicit form where we introduce K j = κ j T/h for j = 1, 2, 3, 4 as rescaled dimensionless kicking strengths. Furthermore, under the quantum resonance conditionhT = 4π that has been considered experimentally [113][114][115][116][117][118][119], we obtain the two-dimensional extension of on-resonance double kicked rotor (or lattice) model, whose Floquet operator readŝ It is then clear that, oncehτ = 2π p/q, with p and q being coprime integers, the Floquet operatorÛ will have translational symmetries in bothn x andn y with the common period q, i.e., a periodic crystal structure in the momentum space of the two-dimensional on-resonance double-kicked lattice. In one-dimensional (1D) descendant models of Equation (7), rich first-order Floquet topological phases have been discovered, which are characterized by large Chern (winding numbers), multiple chiral (dispersionless) edge modes and topologically quantized acceleration in momentum space [107,[120][121][122][123]. These discoveries further motivate us to explore HOTPs in the 2D on-resonance double-kicked lattice model. To obtain a minimal version of our model with nontrivial higher-order topology, we choose the time delay between the two kicks to be τ = T/4, which implies thathτ = π. Moreover, fixing the phase differences at φ x = φ y = π/2, the Floquet operator of the 2D on-resonance double-kicked lattice reduces tô U = e i π 2 (n 2 x +n 2 y ) e −i[K 2 cos(x)+K 4 cos(ŷ)] e −i π 2 (n 2 x +n 2 y ) e i[K 1 sin(x)+K 3 sin(ŷ)] .
By solving the Floquet eigenvalue equationsÛ x |ψ x = e −iE x |ψ x andÛ y |ψ y = e −iE y |ψ y , we could obtain the eigenstates |ψ = |ψ x ⊗ |ψ y ofÛ with eigenphases (quasienergies) E = (E x + E y ) mod 2π. This observation immediately allows us to deduce the possible origin of higher-order topology in the 2D on-resonance double-kicked lattice. That is, if |ψ x and |ψ y are edge modes ofÛ x andÛ y with eigenphases (E x , E y ) = (0, 0) or (E x , E y ) = (π, π), they will be coupled to form a corner mode ofÛ with eigenphase E = E x + E y = 0, i.e., a Floquet corner zero mode. Similarly, if |ψ x and |ψ y are edge modes ofÛ x andÛ y with eigenphases (E x , E y ) = (0, π) or (E x , E y ) = (π, 0), they will be coupled to form a corner mode ofÛ with eigenphase E = E x + E y = π, i.e., a Floquet corner π mode. These are the two types of topological corner modes that could appear in Floquet HOTPs of our model, and their numbers are determined by the numbers of edge modes in the subsystems described byÛ x andÛ y , which are further determined by the winding numbers ofÛ x andÛ y according to the principle of bulk-edge correspondence [124,125].
In the following section, we construct the bulk topological invariants for Floquet HOTPs in the 2D on-resonance double kicked lattice based on these analysis and establish the bulk topological phase diagram of the system.

Topological Invariants and Phase Diagram
Since the Floquet operatorÛ in Equation (8) possesses a tensor product structure, its spectrum and eigenstates are known once the eigenvalue equationsÛ j |ψ j = e −iE j |ψ j for j = x, y are solved. Inserting the identity operators I j = ∑ j |n j n j | in the momentum space, and performing Fourier transforms from the momentum to quasiposition (the conserved quantity due to the translational symmetry n j → n j + 2 in the momentum lattice) representation, the Floquet operatorÛ j can be expressed in the form of where {|θ j } is the eigenbasis of quasiposition with θ j ∈ [0, 2π) and j = x, y. Explicitly, the Floquet matrices U x (θ x ) and U y (θ y ) are given by with shorthand notations Here, σ x,y,z and τ x,y,z are Pauli matrices acting on two sublattice degrees of freedom along the x and y directions in the momentum lattice (see [121,122] for derivation details of U j (θ j ) for 1D descendant models of the 2D on-resonance double-kicked lattice). The standard characterization of 1D Floquet topological phases is achieved by introducing a pair of symmetric time frames upon similarity transformations [124,125]. For our model, there are two symmetric time frames for both U x (θ x ) and U y (θ y ). Putting together, there are in total four such time frames for the Floquet matrix U(θ x , θ y ) = U x (θ x ) ⊗ U y (θ y ) of the 2D system. In these time frames, U(θ x , θ y ) takes the form where α = 1, 2, β = 3, 4, and The auxiliary matrices F, G, F and G are explicitly given by (see [121,122] for derivation details of these matrices for 1D descendant models of our system) With these considerations, it is straightforward to verify that U αβ (θ x , θ y ) possesses the chiral symmetry Γ = σ z ⊗ τ z for all α = 1, 2 and β = 3, 4, in the sense that Γ 2 = 1 and This symmetry then allows us to characterize the Floquet HOTPs of our system by integer topological invariants [121,124,125].
To construct these topological numbers for our system, we take the Taylor expansion for each term of the Floquet matrices in Equations (17) and (18), and recombine the relevant terms. The resulting Floquet matrices take the forms where α = 1, 2 and β = 3, 4. The eigenphase dispersions E x (θ x ) and E y (θ y ) along the two different dimensions are given by The explicit expressions of unit vectors [n αx (θ x ), n αy (θ x )] and [n βx (θ y ), n βy (θ y )] are summarized in Appendix B. In the quasiposition representation, the 2D on-resonance doublekicked lattice then possesses four bulk eigenphase (quasienergy) bands, whose dispersion relations are given by where s, s = ±. The system could undergo topological phase transitions when these bands touch and separate at the quasienergies zero and π.
In previous studies [121,122], it has been demonstrated that a 1D system described by the Floquet operator U α (θ x ) or U β (θ y ) possesses a topological winding number where ν = {α, β}, j = x, y and the winding angle ϕ ν (θ j ) in the ν's time frame is defined as Using these winding numbers, we can further construct two pairs of invariants (w 0x , w πx ) and (w 0y , w πy ) for the Floquet subsystems described byÛ x andÛ y , respectively [124,125]. They are related to the values of w ν through the relations These invariants always take integer quantized values. They provide a complete characterization of the topological phases in 1D Floquet systems with chiral (sublattice) symmetry [124,125]. Moreover, under the open boundary condition, the invariants w 0x (w 0y ) and w πx (w πy ) could predict the numbers of topological edge modes with quasienergies zero and π in the subsystem described byÛ x (Û y ) [121,122], and therefore also capturing the bulk-edge correspondence of these Floquet subsystems. For our 2D system, the Floquet HOTPs can be characterized by appropriate combinations of these 1D topological numbers. Specially, referring to our analysis on how the zero and π Floquet edge modes can be coupled to form corner modes in the last section, we introduce a pair of topological invariants (w 0 , w π ) for the 2D on-resonance double-kicked lattice, which are defined as w π ≡ |w 0x w πy | + |w πx w 0y |. (34) It is clear that (w 0 , w π ) ∈ Z × Z due to the quantized nature of w 0j and w πj (j = x, y). Furthermore, as demonstrated in Section 5, the values of w 0 and w π could correctly count the numbers of Floquet corner modes with quasienergies zero and π in the momentum space of our system. Therefore, the invariants (w 0 , w π ) could provide us with a complete characterization of Floquet HOTPs in the 2D on-resonance double-kicked lattice and other chiral symmetric 2D lattice models whose Floquet operators can be expressed in the form ofÛ =Û x ⊗Û y . We also notice that (w 0 , w π ) = (0, 0) only when the subsystems described byÛ x andÛ y are both topologically nontrivial. The Floquet HOTPs of the 2D on-resonance double-kicked lattice are thus originated from the nontrivial cooperation of topological natures of the two subsystems in lower dimensions.
In the remaining part of this section, based on the evaluation of invariants w 0 and w π in Equations (33) and (34), we present topological phase diagrams of the 2D on-resonance double-kicked lattice for two typical situations. In the first case, we show the phase diagram of the system with respect to the kicking strengths (K 2 , K 4 ) in Figure 1. We observe that with the increase of these kicking strengths, a series of topological phase transitions can be induced, which each of them being accompanied by the quantized jump of w 0 or w π . At large values of (K 2 , K 4 ), we further obtain rich Floquet HOTPs characterized by large values of (w 0 , w π ). It is not hard to verify that there is no upper bound for the values of these invariants if the kicking strengths keep increasing. Therefore, the 2D on-resonance double-kicked lattice serves as a good candidate to generate Floquet secondorder topological phases in momentum space with arbitrarily large topological invariants. Besides, according to the tensor product structure of Floquet operator in Equation (16), the boundaries separating different Floquet HOTPs are determined by the gapless conditions along one dimension of the lattice. More explicitly, these phase boundaries are determined by the conditions cos[E x (θ x )] = ±1 and cos[E y (θ y )] = ±1, or equivalently with m i ∈ Z and |m i π/K i | ≤ 1 for i = 1, 2, 3, 4. In Figure 1, our choice of system parameters yield m 1 = m 3 = 0, and the phase boundaries are reduced to straight lines according to Equation (35), which is also consistent with the numerical results.
In the second case, we present the phase diagram of the 2D on-resonance doublekicked lattice versus kicking strengths (K 3 , K 4 ) in Figure 2. With the increase of these kicking strengths, we also observe rich Floquet HOTPs featured by large invariants (w 0 , w π ), together with multiple transitions between these phases followed by quantized changes of (w 0 , w π ). Numerically, we have also checked the phase diagrams of the system versus any one pair of kicking parameters (K i , K i ), with the other pair being fixed for i, i = 1, 2, 3, 4, and obtain similar patterns for the topological phases and phase transitions. Therefore, we conclude that the 2D on-resonance double-kicked lattice indeed possesses rich Floquet HOTPs, which are characterized by a pair integer topological invariants (w 0 , w π ). These invariants could not only predict the numbers of Floquet zero and π corner modes in the system under open boundary conditions, but also be detectable experimentally from the dynamics of wavepackets, as presented in the following sections.

Mean Chiral Displacements
In this section, we sketch a dynamical approach that can be used to probe the invariants (w 0 , w π ) of Floquet HOTPs in our system. This approach is based on the detection of a 2D extension of the time-averaged mean chiral displacement, which is introduced in [82,95]. In a 2D lattice (either in position or in momentum space), we define the chiral displacement operator of the dynamical evolution in the time frame (α, β) aŝ Here, t counts the number of evolution periods. The Floquet operatorÛ αβ =Û α ⊗Û β , as defined in Equation (16) for α = 1, 2 and β = 3, 4.N x andN y are unit-cell position operators along x and y directions of the lattice. Γ x and Γ y describe chiral symmetries of the descendant systemsÛ α andÛ β along x and y directions, respectively. For our 2D on-resonance double-kicked lattice, we have Γ x = σ z and Γ y = τ z , whose tensor product gives the chiral symmetry operator Γ of the system. To extract the topological winding numbers ofÛ αβ , we initialize the system in a fully polarized state at the central unit cell of the lattice. The initial state vector then takes the form where |0 x (|0 y ) is the eigenstate ofN x (N y ) with eigenvalue N x = 0 (N y = 0). | ↑ x and | ↑ y are the eigenvectors of Γ x and Γ y with eigenvalues +1. After the evolution over a number of t's driving periods in the time frame (α, β), the mean chiral displacement of initial state |ψ 0 is given by the expectation value Since |ψ 0 is not an eigenstate ofÛ αβ , C αβ (t) is expected to be an oscillating function of time. To extract the topological information from C αβ (t), we take the average of C αβ (t) over many driving periods, which in the long-time limit yields the stroboscopic time-averaged mean chiral displacement Following the derivation steps as detailed in [82], it can be shown that for α = 1, 2 and β = 3, 4, Here, w α and w β are the winding numbers defined in Equation (29). With the help of Equations (31) and (32), we can recombine the time-averaged mean chiral displacements to obtain the products of winding numbers w 0x w 0y , w 0x w πy , w πx w 0y and w πx w πy , which finally yield the invariants (w 0 , w π ) of the Floquet HOTPs. From now on, we denote the recombined time-averaged mean chiral displacements that are related to w 0 and w π as C 0 and C π , respectively (see Appendix C for their explicit expressions).
To verify the relations between the time-averaged mean chiral displacements and the topological invariants of the 2D on-resonance double-kicked lattice, we compute and compare (w 0 , w π ) and (C 0 , C π ) for a typical case in the remaining part of this section. In Figure 3, we present the topological invariants and mean chiral displacements with respect to the kicking strength K 4 for evolutions over t = 50 driving periods. We observe that the time-averaged mean chiral displacements (C 0 , C π ) indeed take nearly quantized values in each Floquet HOTPs of the 2D on-resonance double-kicked lattice, which demonstrate the relation between them and the topological invariants of the system, i.e., (C 0 , C π ) = (w 0 , w π ) (see Appendix C for derivation details of this relation). Besides, the values of (C 0 , C π ) also possess quantized jumps around all topological phase transition points (K 4 = π, 2π, 3π, 4π in Figure 3). Therefore, the mean chiral displacements could also serve as a dynamical prob to the phase transitions between different Floquet HOTPs. The deviations of (C 0 , C π ) from quantization are due to finite-time effects, which will gradually go to zero with the increase of the number of driving periods t. Numerically, we checked that the time-averaged mean chiral displacements (C 0 , C π ) remain close to quantization for t = 20 driving periods, which should be within reach in current or near-term experiments in photonic [126,127] and cold atom [128,129] systems. . Topological invariants (w 0 , w π ) and time-averaged mean chiral displacements (C 0 , C π ) of the 2D on-resonance double-kicked lattice versus the kicking strength K 4 . The other system parameters are set as K 1 = 0.5π, K 2 = 3.5π and K 3 = 0.5π. The mean chiral displacements are averaged over t = 50 driving periods to generate the results for (C 0 , C π ).
In the following section, we demonstrate another topological signature of Floquet HOTPs, i.e., the symmetry protected bound states localized around the corners of the momentum-space lattice, and relate their numbers to the topological invariants (w 0 , w π ) of the system.

Floquet Topological Corner Bound States in Continuum
A defining feature of HOTPs in D spatial dimensions is symmetry-protected states localized along its (D − d)-dimensional boundaries, where d > 1. For a 2D lattice studied in this work, such bound states could appear at the geometric corners of the system. Moreover, since the Floquet bands of a chiral symmetric system come in positive and negative pairs ±E, they could touch and separate at the quasienergies zero and ±π. Therefore, in principle, there could be two distinct types of Floquet corner modes appearing at these quasienergies, which are called Floquet zero and π corner modes. In the 2D on-resonance double-kicked lattice described byÛ =Û x ⊗Û y , since these corner modes are formed by the coupling between doubly degenerate edge modes of 1D descendant systemsÛ x andÛ y , their numbers (N 0 , N π ) are always integer multiples of four. Furthermore, the numbers (N 0 , N π ) of Floquet zero and π corner modes tend to be connected with the topological invariants (w 0 , w π ) of the bulk through the relation (N 0 , N π ) = 4(w 0 , w π ). (41) Equation (41) (9) and (10) in momentum lattice representations as (see Appendix A for details) 2i ∑ ny (|n y n y +1|−h.c.) .
The quasienergies E and Floquet eigenstates |ψ of the system are then obtained by solving the eigenvalue equationÛ|ψ = e −iE |ψ , where the diagonalization ofÛ can be performed separately forÛ x andÛ y due to the tensor product structure of the Floquet operator U =Û x ⊗Û y .
In Figure 4, we present the quasienergy spectrum E and number of Floquet corner modes N 0 and N π at zero and π quasienergies in the momentum space of our model. In Figure 4a, we observe an almost continuous Floquet spectrum with no gaps around the quasienergies E = 0 and E = π. This is different from the situations usually observed in 2D static and Floquet HOTPs, where corner modes are separated from bulk states by spectral gaps. However, by evaluating the inverse participation ratio IPR ≡ ∑ n x ,n y |ψ(n x , n y )| 4 (44) for all the Floquet eigenstates ofÛ, we find different numbers of corner modes N 0 and N π at the quasienergies zero and π in distinct Floquet HOTPs of the system. Their numbers are presented together with the topological invariants w 0 and w π in Figure 4b. Note that the inverse participation ratio of corner modes differ from bulk and 1D edge states of the system in their order of magnitudes, and can thus be numerically distinguished from one another. Furthermore, we observe that the relation (N 0 , N π ) = 4(w 0 , w π ) holds throughout the considered parameter regime, validating the bulk-corner correspondence of our model as established in Equation (41). Besides, with the increase of K 4 , the system undergoes a series of topological phase transitions, yielding Floquet HOTPs with more and more zero and π corner modes. The 2D on-resonance double-kicked lattice thus provides us with a nice platform to investigate Floquet HOTPs with multiple corner states and strong topological signatures in momentum space. For completeness, we also checked the quasienergy spectrum and corner modes in other parameter regions of the system and obtained results that are consistent with the above descriptions. In Figure 5, we present the quasienergy spectrum, inverse participation ratio and corner modes of the system for a typical situation. In Figure 5a, we observe that the Floquet spectrum E spread throughout the first quasienergy Brillouin zone, and no gaps can be observed around E = 0, ±π. However, in Figure 5b, we find three clear peaks in the inverse participation ratio versus the quasienergy E around E = 0 and E = ±π, which indicates the existence of localized bound states in the system at these quasienergies. By investigating the data, we obtain eight (four) such bound states at E = 0 (E = ±π), with their probability distributions shown explicitly in Figure 5d,e [ Figure 5c]. We see that all of these bound states are indeed localized around the corners of the system, and their numbers are determined precisely by the bulk-corner correspondence relation in Equation (41) for the given set of system parameters. Therefore, these corner states originate from the higher-order topology of our system. They represent topologically degenerate corner modes in momentum space of the system, which are protected by the chiral symmetry Γ. Besides, these corner modes coexist with extended bulk states at the same quasienergies. We therefore refer to them as corner-localized Floquet topological bound states in continuum, in the sense that they do not hybridize with the surrounding bulk states even without a bulk band gap. This observation also extends the scope of bulk-corner correspondence of Floquet HOTPs to more general situations, in which the symmetry-protected corner modes can not only be found in spectral gaps, but also appear within topological Floquet bands. Note in passing that the corner-localized bound states in continuum with zero energy have been discovered before in static HOTPs [130][131][132]. By contrast, the corner-localized Floquet bound states in continuum subject to different topological characterizations, and the corner bound states in continuum with quasienergy E = π are unique to Floquet HOTPs found in this work. Experimentally, the implementation of open boundary conditions might be achieved in Bose-Einstein condensates by a setup similar to the one employed in the realization of topological quantum walk in momentum space [129]. Meanwhile, by applying the mapping introduced in [121] from momentum to position space lattices, the 2D on-resonance double-kicked lattice may also be realized in real space, and the corner modes may then be probed in setups such as photonic waveguide arrays [132]. Figure 5. Quasienergy spectrum, inverse participation ratio and Floquet zero/π corner modes of the 2D on−resonance double−kicked lattice. System parameters are set as K 1 = 0.5π, K 2 = 3.5π, K 3 = 0.5π and K 4 = 1.5π. The size of the lattice is L x = L y = 300. (a) n denotes the index of the state. The peaks around E = 0 (E = ±π) in (b) correspond to the inverse participation ratio of zero (π) Floquet corner modes, whose probability distributions are shown in panels (c-e).
To demonstrate the robustness of Floquet second-order topological phases found in our system, we add disorder to the nearest neighbor hopping amplitudes in the momentum lattice, i.e., by letting the kicking strengths K 1,2 → K 1,2 + r 1,2 (n x ) and K 3,4 → K 3,4 + r 3,4 (n y ) to be dependent on the lattice site indices n x and n y . Here, r 1,2 (n x ) and r 3,4 (n y ) are uniformly distributed random numbers in the range [−W/2, W/2] at different site indices (n x , n y ). In numerical calculations, we set W = 0.1 as the strength of disorder and present a representative example of the Floquet spectrum and corner states with such hopping disorder in Figure 6. We observe that with disorder, the Floquet corner modes are still preserved, as shown in Figure 6c-e, and their quasienergies are also pinned at E = 0 and E = ±π. They are thus well-defined Floquet zero and π corner modes and their presence demonstrate the robustness of Floquet second-order topological phases of our system with disorder. Furthermore, the total number of corner modes with quasienergy zero (π) is eight (four), as shown in Figure 6c-e, which is precisely predicted by the bulk-corner correspondence relation in Equation (41). Therefore, the bulk topological invariants w 0 and w π introduced in Equations (33) and (34) could also describe the Floquet second-order topological phases of our system in the presence of weak hopping disorder. Figure 6. Quasienergy spectrum, inverse participation ratio and Floquet zero/π corner modes of the 2D on−resonance double−kicked lattice with hopping disorder. System parameters are set as K 1 = 0.5π + r 1 , K 2 = 3.5π + r 2 , K 3 = 0.5π + r 3 and K 4 = 1.5π + r 4 , where r 1,2 (r 3,4 ) take random values in the range [−W/2, W/2] with W = 0.1 at different nearest neighbor bonds along n x (n y ) direction. The size of the lattice is L x = L y = 300. (a) n denotes the index of the state. The peaks around E = 0 (E = −π) in (b) correspond to the inverse participation ratio of zero (π) Floquet corner modes in the presence of disorder, whose probability distributions are shown in (c-e).

Conclusions
In this work, we find rich Floquet second-order topological phases in a two-dimensional extension of the on-resonance double kicked rotor. Each of the Floquet phases is characterized by a pair of integer topological invariants (w 0 , w π ), which can be extracted from the dynamics of the system in four distinct symmetric time frames. The values of invariants w 0 and w π take quantized jumps whenever the system undergoes a transition between different phases. Furthermore, Floquet second-order topological phases characterized by arbitrarily large (w 0 , w π ) can be found in principle with the increase of kicking strengths. Experimentally, the invariants (w 0 , w π ) could also be obtained by measuring the timeaveraged mean chiral displacements of initially localized wavepackets in different time frames. Under open boundary conditions, corner-localized modes with quasienergies zero and π are found to be coexisting with extended bulk states at the same quasienergies, realizing second-order Floquet topological bound states in continuum. The numbers of these corner modes are further determined by the bulk topological invariants (w 0 , w π ), leading to the bulk-corner correspondence of Floquet second-order topological phases. Experimentally, the proposed model may also be engineered in two-dimensional photonic or cold atom systems, where the Floquet corner modes could be imaged either in realspace [132] or in momentum space [129]. Putting together, we uncover a unique set of Floquet second-order topological phases in momentum space, which are featured by large topological invariants, rich phase diagrams and multiple Floquet topological corner bound states in continuum. In future studies, it would be interesting to extend our results to Floquet HOTPs in higher dimensions, different symmetry classes and superconducting systems. The possibility of realizing topological time crystals by superposing the zero and π Floquet corner modes and the potential applications of these corner bound states in Floquet quantum computing tasks also deserve more systematic explorations. Acknowledgments: L.Z. acknowledges Jiangbin Gong for helpful comments.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Floquet Operator in Momentum Representation
In this appendix, we give the expressions of Floquet operatorsÛ x andÛ y in Equations (9) and (10) in the momentum-lattice representation and further derive their forms in the symmetric time frames.
In the momentum basis {|n x }, we can express each term in Equation (9) of the main text as U x± ≡ e ±(π/2)n 2 where we insert the identityÎ x = ∑ n x ∈Z |n x n x | in order to arrive at each expression, as well as write the basis with even and odd lattice indices separately. Since Equation (9) is invariant under the translation over two sites inn x , we could decompose the momentum space lattice into two chains containing only odd and even sites, denoted by sublattice indices A (for odd sites) and B (for even sites). A unit cell of the momentum space lattice now contains two sublattice sites. We can then introduce Pauli matrices in the sublattice representation as The raising and lowering operators in sublattice basis can also be written as In such a bipartite lattice representation, Equations (A1)-(A3) can be expressed as where x is the unit cell index alongn x -direction in momentum space, and we make the identifications for n x ∈ Z. In momentum space, we can then write Equation (9) aŝ Following the same procedure, we can also express Equation (10) in momentum representation asÛ where U y± = e ±iπ/4 e ±i(π/4) ∑ y ∈Z | y y |τ z , U y3 = e i(K 3 /2i) ∑ y ∈Z (| y y |τ + +| y y +1|τ − −H.c.) , (A13) U y4 = e −i(K 4 /2) ∑ y ∈Z (| y y |τ + +| y y +1|τ − +H.c.) , with τ ± = (τ x ± iτ y )/2 and τ x,y,z being Pauli matrices acting in the subspace of the two sublattices alongn y -direction. The Floquet operators in quasiposition representations, as given by Equations (12) and (13) of the main text, can further be obtained by performing the Fourier transform | x,y = (1/ N x,y ) ∑ θ x,y e −iθ x,y x,y |θ x,y , with N x,y being the number of unit cells alongn x -andn y -directions.
To obtain the Floquet operatorsÛ α andÛ β in symmetric time frames α = 1, 2, β = 3, 4 and in the momentum-lattice representation, which are used in the calculation of the chiral displacement in Equation (36), one can simply perform similarity transformations to Equations (A10) and (A11), yieldinĝ

Appendix C. Relations between Mean Chiral Displacements and Topological Invariants
In this appendix, we present the explicit relationship between the time-averaged mean chiral displacements and the topological invariants of our model. Following Equations (31), (32) and (40) in the main text, it is straightforward to verify the following equalities between the time-averaged mean chiral displacements and winding numbers in different combinations of symmetric frames C 13 + C 14 + C 23 + C 24 = w 0x w 0y , (A27) Combining these equations with the definitions of w 0 and w π in Equations (33) and (34) of the main text, we find Therefore, by measuring the time-averaged mean chiral displacements C αβ for α = 1, 2 and β = 3, 4, we could obtain the invariants (w 0 , w π ) of Floquet HOTPs from the wavepacket dynamics in different time frames.