$q$th-root non-Hermitian Floquet topological insulators

Floquet phases of matter have attracted great attention due to their dynamical and topological nature that are unique to nonequilibrium settings. In this work, we introduce a generic way of taking any integer $q$th-root of the evolution operator $U$ that describes Floquet topological matter. We further apply our $q$th-rooting procedure to obtain $2^n$th- and $3^n$th-root first- and second-order non-Hermitian Floquet topological insulators (FTIs). There, we explicitly demonstrate the presence of multiple edge and corner modes at fractional quasienergies $\pm(0,1,...2^{n})\pi/2^{n}$ and $\pm(0,1,...,3^{n})\pi/3^{n}$, whose numbers are highly controllable and capturable by the topological invariants of their parent systems. Notably, we observe non-Hermiticity induced fractional-quasienergy corner modes and the coexistence of non-Hermitian skin effect with fractional-quasienergy edge states. Our findings thus establish a framework of constructing an intriguing class of topological matter in Floquet open systems.

In a periodically driven system, the central object for the description of topological properties is the Floquet operator, which is the evolution operator of the system over a complete driving period T . Taking the square-root of such a propagator for the purpose of generating its topological descendant is, however, a highly nontrivial task. This can be seen by writing the Floquet operator However, such a trial of generating square-root Floquet topological phases tends out to be problematic and useless. First, the exact form of H eff can be rather complicated (usually including driving-induced long-range coupling terms), not physically obtainable, or even insufficient to describe Floquet phases with no static counterparts such as those possessing anomalous Floquet edge modes [42][43][44]. Second, there are no transparent ways to find H eff from H(t), i.e., the parameters in H eff are usually nonlinear combinations of physical parameters in H(t), such that simply reducing the parameters of H(t) by half could not yield H eff /2. Even obtained, the H eff and H eff /2 describe essentially the same physical system up to a global constant, and no new physics are expected to emerge following such a halving process. Therefore, the straightforward operation, U = e −i T ħ h H eff 2 , does not generate a desired square-root of the parent system U.
To resolve this puzzle, a nontrivial route of taking the square-root for U is introduced [45], which closely follows the original idea of Dirac by adding internal degrees of freedom for electrons before taking the square-root of their relativistic wave equation. However, the general applicability of this idea to the construction of Floquet models beyond taking 2 n th-root has not been revealed. Motivated by this gap of knowledge, we propose a generic procedure to yield a variety of qthroot Floquet phases, where q is any arbitrary integer, not necessarily in the form of 2 n . This is achieved by utilizing a Z q generalization of Pauli matrices as ancillary degrees of freedom. While our construction is applicable to any periodically driven systems, we focus on two timely examples of non-Hermitian Floquet matter as case studies.
This paper is structured as follows. In Sec. 2, we recap the strategy of Ref. [45], generalize it to the construction of any qth-root Floquet system, and elaborate the application of this general construction for the case of q = 3. In Sec. 3, we introduce two typical models of first-and secondorder non-Hermitian Floquet topological insulators, whose square-and cubic-root descendants are studied in detail in Sec. 4 as an application of our method. In Sec. 5, we sum up our results and discuss potential future directions.

Theory
We first review the approach to take the nontrivial square-root of a Floquet system. We set the Planck constant ħ h = 1 and driving period T = 1 throughout. Following Ref. [45], we write the one-period evolution (Floquet) operator of any time-periodic system as where H(t) is the system's Hamiltonian. The procedure of Ref. [45] is to first enlarge Hilbert space of H(t) by introducing a pseudospin-1/2 degree of freedom with the corresponding Pauli matrices τ x, y,z . A two-step Hamiltonian is next defined in the enlarged Hilbert space as where ∈ . τ 0 is the identity in the pseudospin-1/2 subspace. 0 is the identity in the Hilbert space of H(t). The Floquet operator of the evolution in the enlarged Hilbert space reads Performing the Taylor expansion and introducing τ ± = (τ x ± iτ y )/2, we find and Since U 1 U 2 = U and U 2 U 1 = U 2 U U −1 2 are related by a similarity transformation, they describe the same parent Floquet system up to a half-period shift of the initial evolution time. The system described by U 2 1/2 can thus be viewed as two equivalent copies of U up to a global phase shift π. Therefore, U 2 1/2 and U are expected to share the same topological features concerning the stroboscopic dynamics. We could view U 1/2 as a nontrivial square-root of U in the spirit of Dirac's taking square-root to reach his equation for electrons [28]. The dynamical and topological properties of U can further be carried over to U 1/2 , which are confirmed by explicit studies of Floquet topological superconductors and time crystals [45].
Iterating the same procedure, we can construct the 2 n th-root of U, i.e., U 1/2 n for any n ∈ + . For example, we could generate U 1/4 by letting 2 ), H 2 = πτ y ⊗ 0 , and defining in a further enlarged Hilbert space where τ y,z and τ 0 are Pauli matrices and identity matrix acting in the subspace of an extra pseudospin-1/2. 0 denotes the identity in the Hilbert space of H 1,2 . The resulting Floquet operator, then defines the nontrivial 4th-root of U. It is straightforward to verify that whose diagonal blocks contain four equivalent copies of U up to a unitary transformation and a global phase π. The extension of the above approach to find any qth-root of U can be achieved by introducing higher-dimensional pseudospin degrees of freedom, i.e., the generalized q × q Pauli matrices where ω = e i2π/q . These operators satisfy where η 0 is the identity matrix acting in the pseudospin subspace. Our qth-rooting procedure can then be executed in two steps. First, given any time-periodic Hamiltonian H(t), we divide the Floquet operator into q time-steps, i.e., Next, we define a two-step Hamiltonian of the form Figure 1: Schematic of a cubic-root system obtained from the procedure of Eq. (12). In the parent system, a particle (marked by black star) evolves under H 3 , H 2 , and H 1 over the course of a period. In the corresponding cubic-root system, a particle living in a given subsystem effectively evolves under the same three Hamiltonians only when viewed over the course of three periods.
where M It then follows that the associated Floquet operator is such that U q 1/q is block diagonal and consists of all q permutations of q =1 U . This shows that the resulting system indeed represents the qth-rooted version of U.
Intuitively, the above construction can be understood as follows. First, note that Eq. (12) is defined on a system consisting of n subsystems. Consider a particle initially living in the jth subsystem. During the first half of the period, evolution under Eq. (12) amounts to transporting the particle towards subsystem j − 1 mod q. By noting that P (q) j represents a projection onto the jth subsystem, it then follows that Eq. (12) further evolves the particle under H j−1 (t) during the second half of the period. In the next Floquet cycle, the particle continues moving to subsystem j − 2 mod q, followed by the half-period evolution under H j−2 (t). As the process continues, at the end of q periods, the particle returns to the subsystem j mod q, while accumulating U j · · · U q U 1 · · · U j−2 U j−1 , which is unitarily equivalent to U.
Having demonstrated the generality of our construction, we will focus on square-root and cubic-root systems for brevity in the remainder of this paper. To this end, we will now present an explicit application of the above construction to obtain a nontrivial cubic-root of a system relevant to the case studies below. Specifically, such a parent system follows a three-step periodically quenched drive, whose time-dependent Hamiltonian takes the form with the corresponding Floquet operator of Note that the Floquet operator associated with a system following a two-step periodically quenched drive, whose Hamiltonian switches between h 1 and h 2 after every half period, can also be cast into the form of Eq. (15) by shifting the initial evolution time from t = 0 to 3/4 and identifying H 1 , H 3 = 3h 1 /4 and H 2 = 3h 2 /2. In both cases, the cubic-root of U in Eq. (15) can be obtained ac- We further identifyH j = with j = 0, 1, 2, as well as η x,z in the explicit matrix forms Indeed, it can be directly verified that Eq. (16) satisfies the algebra of Eq. (10). The corresponding Floquet operator of the cubic-root model then reads with That is, the three diagonal blocks of U 3 1/3 only differ from one another by the starting time of evolution and describe equivalent Floquet systems in stroboscopic dynamics concerning the spectral and topological properties. This implies that U 1/3 is indeed a nontrivial cubic-root of the parent system U in Eq. (15). The presented cubic root procedure is schematically depicted in Fig. 1.
We now discuss how the rooted Floquet system could inherit the symmetry protected edge states of the parent model while altering their quasienergies to rational fractions of 2π. A key symmetry that is relevant to the topological characterization of the parent systems considered in this work is the chiral symmetry (CS). If a general Floquet operator U possesses the CS, there is a unitary operator Γ such that Therefore, Γ |ψ〉 is an eigenstate of U with quasienergy −E. Now if there is an eigenstate |ψ〉 with E = 0 or π, the CS enforces the presence of another eigenstate Γ |ψ〉 also at E = 0 or π (E = ±π are identified as the same quasienergy since E is defined mod 2π), yielding eigenstate degeneracy at the center or boundary of the quasienergy Brillouin zone E ∈ [−π, π]. If such eigenmodes appear at the edge or corner of the system, we obtain CS-protected degenerate edge or corner modes of U.
For the square-root system U 1/2 , we already see that U 2 1/2 is block diagonal and its two diagonal blocks share the same spectral and topological properties with the parent model U. If |ψ 〉 is an eigenstate of U 1/2 with quasienergy E , i.e., U 1/2 |ψ 〉 = e −iE |ψ 〉, it is straightforward to see that |ψ 〉. U 1/2 and U 2 1/2 thus share the same eigenbasis. When the parent model U possesses the CS Γ , the diagonal blocks of U 2 1/2 possess the CS, such that U 2 1/2 is chiral symmetric with respect to Γ = τ z ⊗ Γ . Degenerate topological edge/corner modes of U 2 1/2 can thus only appear at E = 0, ±π = 2E mod 2π. This implies that in the square-root system U 1/2 , we could only find topological edge/corner modes at the quasienergies E = 0, ±π/2, ±π, which are indeed protected by the CS Γ of the parent model. Interestingly, the degenerate eigenmodes at E = ±π/2 are present only in the U 1/2 and are thus unique to the square-root Floquet system.
While the topological protection of the edge/corner modes in U 1/2 can be understood from the presence of chiral symmetry in its corresponding parent system U, it would also be insightful to discuss the protecting symmetries that arise at the level of U 1/2 directly. To this end, we first note that Γ = τ z ⊗ Γ is also a chiral symmetry with respect toŨ 1/2 = e −i π 4 U 1/2 e i π 4 , i.e., U 1/2 under the shift in the initial time from t = 0 to t = 1/4. Similar to its parent counterpart, such a chiral symmetry is responsible for protecting E = 0, ±π quasienergy edge states in the square-root system. Next, we identify an additional symmetry Γ 1/2 = τ z ⊗ I (I being the identity operator), which acts only within the enlarged degree of freedom and is thus referred to as the "subchiral" symmetry [79]. Such a symmetry operates as In this case, a quasienergy which satisfies E±π = −E, i.e., E = ±π/2, is necessarily twofold degenerate due to the product Γ Γ 1/2 . The associated quasienergy eigenstates could further be chosen to be simultaneous ±1 eigenstates of Γ Γ 1/2 . This is automatically the case for the quasienergy ±π/2 edge/corner states. In particular, since ±1 eigenstates of Γ Γ 1/2 correspond to states localized at two opposite edges/corners, the discreteness of Γ Γ 1/2 eigenstates pins such edge/corner states at quasienergy ±π/2 in the presence of symmetry-preserving perturbations. This completes the symmetry protection analysis of quasienergy ±π/2 edge/corner states in the square-root system. The above argument can be easily extended to conclude that, for any qth-root version of the system, U q 1/q also possesses the CS with respect to Γ = η z ⊗Γ . A generalized "subchiral" symmetry can further be identified as Γ 1/q = η z ⊗I, which operates as Γ 1/q U 1/q Γ † 1/q = ω † U 1/q and thus forces the quasienergies of U 1/q to form a cluster of E −2π j/q with j = 0, 1, · · · , q−1. In the case of cubicroot Floquet systems, which are explicitly studied below, both symmetries lead to the protection of degenerate edge/corner modes at E = 0, ±π/3, ±2π/3, ±π. Following the same routine, we can deduce that if the parent model U possesses the CS Γ , the existence of its edge/corner modes at the quasienergies E = 0, π guarantees the presence of degenerate edge/corner states at the quasienergies (0, 1, ...2 n )π/2 n and (0, 1, ..., 3 n )π/3 n of the systems described by U 1/2 n and U 1/3 n , respectively. Notably, the boundary modes appearing at the fractional quasienergies pπ/q with p < q and (p, q) being co-prime integers are, to the best of our knowledge, not identified by previous studies on the symmetry classification and bulk-boundary correspondence of Floquet systems. They are thus a unique product of our qth-root procedure operated on Floquet operators. In Sec. 4, we will apply our theory to explicitly construct square/cubic-root first-and second-order non-Hermitian FTIs based on the parent models defined in the following section.

Parent models
In this section, we introduce two non-Hermitian Floquet topological insulator (FTI) models that will be taking square-and cubic-roots. Detailed investigations of these parent models can be found in Refs. [62] and [66]. All system parameters below are assumed to be properly scaled and set in dimensionless units.
The first model of our consideration describes a non-Hermitian FTI with rich topological phase diagrams and arbitrarily many degenerate edge modes in the presence of Floquet NHSE [66]. Its time-dependent Hamiltonian is where t denotes time, ∈ , and Here n ∈ is the unit cell index. σ x, y,z are Pauli matrices acting on the two sublattices in each unit cell. J 1,2 and iλ describe symmetric and asymmetric parts of intercell hopping amplitudes. µ is the intracell coupling strength. The Floquet operator U = e −i 1 2 H 1 e −i 1 2 H 2 that governs the evolution of the system over a driving period (e.g., from t = 0 to 1) is nonunitary once λ = 0. This yields a model that could possess non-Hermitian FTI phases, which are characterized by integer or half-integer quantized topological invariants under the periodic boundary conditions (PBC) [66]. Under the open boundary conditions (OBC), the CS of the model Γ = N ⊗ σ z (N is the number of unit cells and N is an N × N identity) allows multiple edge modes to appear in pairs at the quasienergies zero and π, whose numbers can be determined by the OBC bulk winding numbers ν 0 and ν π (see Sec. A for their definitions). These edge modes are further found to coexist with sufficient amounts of bulk states localized around both edges of the system due to the NHSE [66].
The second model that we will employ describes a non-Hermitian Floquet second-order topological insulator (FSOTI), which could possess multiple quartets of corner-localized states at real quasienergies zero and π [62]. The Hamiltonian of the model takes the form of The x and y are identity matrices for the basis along x and y directions of the lattice. σ x, y,z are Pauli matrices and σ − = (σ x − iσ y )/2. m, n ∈ are unit cell indices along the two spatial dimensions. ∆ and J 1,2 describe hopping amplitudes between nearest neighbor cells along the x and y directions. µ characterizes the strength of an onsite potential bias. Gain and loss are introduced to make the system non-Hermitian by setting µ = u + i v, with u, v ∈ and v = 0. The Floquet operator of the system takes the form U = e −i 1 2 H 1 e −i 1 2 H 2 , whose spectrum under the OBC features fourfold degenerate topological corner modes at zero and π quasienergies. The numbers of these corner modes n 0 and n π are related to a pair of bulk topological winding numbers ν 0 and ν π of U (see Sec. B for their definitions) through a bulk-corner correspondence relation (n 0 , n π ) = 4(|ν 0 |, |ν π |) [62]. The fourfold degeneracy of Floquet corner modes at the quasienergies E = 0, π is protected by the CS Γ = σ z ⊗ σ y of the two-dimensional system described by U under the PBC [62].
Applying the procedure of Sec. 2, we will obtain the square and cubic roots of the two Floquet models introduced here, and unveil the intriguing topological features of the resulting systems in the following section. As will be demonstrated, our qth-root procedure endows the non-Hermitian Floquet phases in the above two parent models with even richer topological properties.

Results
In Sec. 4.1, we present square-and cubic-root non-Hermitian FTIs generated by the first model in Sec. 3, which will be shown to possess multiple and tunable numbers of degenerate edge modes with the quasienergies π/2, π/3 and 2π/3 that could survive under the NHSE. In Sec. 4.2, we discuss square-and cubic-root non-Hermitian FSOTIs yielded by the second model in Sec. 3, which hold non-Hermiticity induced quartets of topological corner modes at the π/2, π/3 and 2π/3 quasienergies.

Square/Cubic-root non-Hermitian FTIs
We now apply the procedure in Sec. 2 to find the square-and cubic-roots of the first model in Sec. 3. In the lattice representation, the square-root Floquet system is obtained by identifying U 1 = e −iH 1 /2 and U 2 = e −iH 2 /2 in Eq. (4), where the H 1 and H 2 are given by Eqs. (19) and (20), respectively. The Floquet operator U 1/2 is then derived following Eq. (1), i.e., To obtain the cubic-root model, we may identifyH 1 =H 3 = 3H 1 /4 andH 2 = 3H 2 /2 in Eq. (17), where H 1 and H 2 are defined by Eqs. (19) and (20), respectively. This then leads to the Floquet operator Solving the eigenvalue equations U 1/2(1/3) |ψ〉 = e −iE |ψ〉 under the OBC, with E being the quasienergy, provides us with all bulk and edge states of the square-(cubic-) root Floquet system.
As an important note, if there is a pair of degenerate edge modes with zero-quasienergy in the parent model U = e −iH 1 /2 e −iH 2 /2 , their quasienergies will be shifted to π in U 2 1/2 according  to Eq. (5), yielding edge modes at quasienergies E = ±π/2 in the system described by U 1/2 . On the other hand, if a pair of degenerate edge modes appears at E = π in the parent model, their quasienergies will become 0 (mod 2π) in U 2 1/2 , leading to edge states at E = 0, ±π in the squareroot model. Following the same routine, we deduce that the edge modes at zero (π) quasienergy in U = e −iH 2 /4 e −iH 1 /2 e −iH 2 /4 generate edge states with E = 0, ±2π/3 (E = ±π/3, ±π) in the system described by U 1/3 after taking the cubic root. Now if we could relate the numbers of zero and π edge modes in the parent system U to its topological invariants, these invariants should also predict the numbers of zero, π/2, π/3, 2π/3 and π modes in the square-and cubic-root systems if the symmetry that protects their quantization is not broken during the process of taking roots.
To showcase the fractional-quasienergy edge modes in the spectrum in a more transparent manner, we introduce the gap function F with respect to a quasienergy , which is defined as Note that the E in Eq. (27) is the collection of all quasienergies obtained by diagonalizing the Floquet operator of the system under consideration. It is clear that once there is an edge state with real quasienergy that resides in a gap on the complex plane, we would have F = 0 for that state and F > 0 for all other bulk states. To locate the expected edge states of U 1/2 and U 1/3 , we choose = 0, π/2, π and = 0, π/3, 2π/3, π for them, respectively, in the following numerical calculations.
In Fig. 2, we present the quasienergy (Floquet) spectrum and gap functions of the first model in Sec. 3 and its square-root descendant under both the PBC and OBC. The quasienergies and gap functions of the parent model in Figs. 2(a) and 2(c) are reproduced from Ref. [66]. A clear distinction between the spectrum under PBC (gray dots in the background) and OBC (blue dots) can be observed especially around the phase transition points, implying the presence of NHSE in the system. To retrieve the bulk-edge correspondence, a pair of open-boundary winding numbers (ν 0 , ν π ) is introduced in Ref. [66] and reviewed in Sec. A, which correctly counts the number of zero-and π-quasienergy edge modes n 0 and n π in the parent model through the relation (n 0 , n π ) = 2(|ν 0 |, |ν π |). Here n E denotes the number of edge states at the quasienergy E. According to our square-root procedure, the edge modes at the quasienergies E = 0, ±π (E = ±π/2) are generated by taking the square-root of the π (zero) Floquet edge modes. Therefore, we arrive at the following bulk-edge correspondence for the square-root FTIs described by U 1/2 , i.e., n π/2 = 2|ν 0 |, n 0 = n π = 2|ν π |, where n π/2 means the number of degenerate edge states at E = ±π/2. These relations are readily confirmed by comparing the spectrum and gap functions presented in Figs. 2(b,d) and Figs. 2(a,c). Notably, with the increase of hopping amplitude J 1 , we observe a series of gap closing and topological phase transitions in the square-root model. After each transition, the number of edge modes n 0 , n π or n π/2 is found to be increased by 2 even in the presence of NHSE. Specifically, we find (n π/2 , n 0 , n π ) = (0, 0, 0), (2, 0, 0), (2, 2, 2), (2, 4, 4), (4, 4, 4), (6, 4, 4), (6, 6, 6) with the increase of J 1 in Fig. 2(d), meanwhile the winding numbers are (ν 0 , ν π ) = (0, 0), (1, 0), (1, −1), (1, −2), (2, −2), (3, −2), (3, −3) according to the calculation reported in Ref. [66]. This process could continue with the further increase of J 1 . We can thus in principle obtain arbitrarily many topological edge modes at fractional quasienergies E = ±π/2 in our square-root non-Hermitian Floquet system. This highlights the universal advantage of Floquet engineering in generating unique nonequilibrium states with strong topological signatures.  In Figs. 3(a) and 3(b), we further show the Floquet spectrum and gap function of the cubicroot non-Hermitian FTI. In addition to edge states at the quasienergies E = 0, ±π, we also observe degenerate edge modes at fractional quasienergies E = ±π/3, ±2π/3. Recall that Eq. (26) cubes to a block diagonal matrix consisting of multiple copies of Floquet operator of the parent model U. Therefore, the edge states at quasienergies (0, ±2π/3) [(±π/3, ±π)] are indeed descendants of the zero [π] edge modes in the parent model, whose numbers are counted by ν 0 [ν π ] [66]. We then obtain the bulk-edge correspondence for cubic-root non-Hermitian FTIs as With the increase of J 1 , the cubic-root system could also undergo a series of topological phase transitions, with each of them being accompanied by the increase of either n 0 and n 2π/3 or n π/3 and n π by two. We can thus obtain arbitrarily many π/3 and 2π/3 degenerate edge modes by tuning the single driving parameter J 1 even in the presence of NHSE. Since it has been demonstrated that Floquet edge states could be utilized to construct boundary discrete time crystals (DTCs) [22,80], we expect the emergence of unique non-Hermitian Floquet boundary DTCs through the superposition of (±π/3, ±2π/3) edge modes and other edge states in the cubic-root FTIs. For completeness, we present the spectrum and gap function in Figs. 3(c) and 3(d) for the fourth-root non-Hermitian FTI, which is constructed by applying the procedure in Eqs. (6) and (7) to the first model in Sec. 3. The resulting system holds Floquet edge states at E = ± π/4 with = 0, ..., 4. Similar to our analysis of U 1/2 and U 1/3 , the number of these edge modes are related to the bulk topological invariants (ν 0 , ν π ) of the parent Floquet system via n π/4 = n 3π/4 = 2|ν 0 |, n 0 = n π/2 = n π = 2|ν π |.

Square/Cubic-root non-Hermitian Floquet second order topological insulators
We next generate square-and cubic-root non-Hermitian Floquet second-order topological insulators (FSOTIs) by applying our theory in Sec. 2 to the second model in Sec. 3. To find the square-root model, we identify U 1 = e −iH 1 /2 and U 2 = e −iH 2 /2 in Eq. (4), where the H 1 and H 2 are given by Eq. (21) in Sec. 3. The resulting square-root Floquet operator reads  Similarly, to obtain the cubic-root model, we identifyH 1 =H 3 = 3H 1 /4 andH 2 = 3H 2 /2, where the H 1 and H 2 are given by Eq. (21). This leads to the cubic-root Floquet operator Recall that the parent Floquet system U = e −i 1 2 H 1 e −i 1 2 H 2 possesses multiple and non-Hermiticity induced fourfold degenerate corner modes at zero and π quasienergies, which are obtained by solving the eigenvalue equation U|ψ〉 = e −iE |ψ〉. The numbers of these corner modes are related to a pair of topological invariants introduced in Ref. [62] (see also Sec. B). Since the process of taking the square (cubic) root of U does not break the protecting CS of these corner modes, we expect the topological invariants of U to be able to predict the numbers of corner modes at the (0, π/2, π) [(0, π/3, 2π/3, π)] quasienergies of The spectra of U 1/2 and U 1/3 are presented in Fig. 5. In Figs. 5(a) and 5(c), we observe states at the fractional-quasienergies E = ±π/2 and E = ±π/3, ±2π/3 for the square-root and cubic-root systems respectively, whose spatial profiles are found to be localized around the four corners of the lattice. The numbers of these corner modes can further be controlled by changing the real part of onsite potential u, as can be seen clearly in Figs. 5(b) and 5(d). The critical values (u 1 , u 2 ) where the number of corner modes change correspond to topological phase transition points, which are determined by the gapless condition of the parent model [62], i.e., Moreover, since the π/2 (zero and π) modes of U 1/2 are inherited from the zero (π) modes of the parent model U, their numbers are determined by the winding numbers (ν 0 , ν π ) of the second model in Sec. 3 (see also Sec. B) through the bulk-corner correspondence relations n π/2 = 4|ν 0 |, n 0 = n π = 4|ν π |.
Finally, in Fig. 7, we present the gap functions and spatial profiles of Floquet corner modes at the quasienergies π/2, π/3 and 2π/3. Their numbers are found to precisely coincide with the bulk-corner correspondence relations in Eqs. (34) and (35). The non-Hermiticity enriched higher-order topology in rooted Floquet systems may also assist us to engineer unique DTCs and quantum computing schemes with the multiple quartets of corner modes at different fractionalquasienergies that are robust to the perturbation of environment.

Conclusion
In summary, we proposed a systematic approach to construct the qth-root of any periodically driven system and presented 2 n th-and 3 n th-root Floquet topological insulators as explicit examples. The latter are shown to exhibit degenerate edge/corner modes at fractional quasienergies π 2 n (0, 1, ..., 2 n ) and π 3 n (0, 1, ..., 3 n ), whose topological nature is inherited from their 2 n th-and 3 n thpower parent systems. Square-and cubic-root non-Hermitian Floquet topological insulators with multiple and tunable topological edge/corner states at quasienergies π/2, π/3 and 2π/3 were investigated in details. Further connections were made between the number of these edge/corner modes and the bulk topological invariants of parent systems, yielding the bulk-edge/corner correspondence in two classes of rooted Floquet topological insulators. Intriguingly, non-Hermitian effects are found to induce more corner modes with fractional-quasienergies and generate multiple edge states coexisting with the non-Hermitian skin effect in rooted systems. Our discoveries thus uncover a unique class of topological phases that originates from the cooperation among driving, non-Hermiticity and the process of taking the nontrivial roots of Floquet systems.
From the experimental perspective, the obtained systems from our qth-rooting procedure are expected to be implementable with the same setups employed for realizing their parent models. To this end, the additional degrees of freedom required in our qth-rooting procedure can be principally implemented by coupling multiple copies of the parent system. In the context of the non-Hermitian Floquet first and second order topological insulators explored in this paper, their corresponding square-and cubic-root systems can thus be realized in principle via setups like photonic quantum walks [78,[81][82][83][84]. For example, the anisotropic hopping amplitude and non-Hermitian lattice potential can be implemented by introducing controlled optical losses with acousto-optical modulators in coupled optical fibre loops [78]. Moreover, the winding numbers used in characterizing their topology can in principle be experimentally probed via measuring mean chiral displacements [59,61,85] or time-averaged spin textures [60,86], which can also be conducted in similar photonic setups [81][82][83][84].
In future work, it would be interesting to apply our scheme to realize qth-root chiral Floquet topological insulators and gapless topological phases in higher spatial dimensions. The application of our approach to systems with many-body interactions is also expected to be fruitful. In particular, it was recently shown that the interplay between interaction and periodic driving may promote 2π/2 n modes into Z 2 n parafermions [87]. Other families of 2π/q modes obtained in this work thus open avenues for exploring different types of Floquet parafermions not covered in Ref. [87]. In particular, Z 3 parafermions, which are expected to arise in systems with 2π/3 modes when subjected to appropriate interactions, form a main ingredient for constructing the powerful Fibonacci anyons [88] that enable topologically protected universal quantum computation.

A Topological invariants of the non-Hermitian FTI
Here we briefly recap the open-boundary winding numbers (OBWNs) introduced in Ref. [66]. They will be used to establish the bulk-edge correspondence for the first parent model in Sec. 3 and its qth-root descendants in Sec. 4.1. We first consider the dynamics of the model in two symmetric time frames, where the Floquet operator U = e −iH 1 /2 e −iH 2 /2 is transformed to U a = e −iH 2 /4 e −iH 1 /2 e −iH 2 /4 and U b = e −iH 1 /4 e −iH 2 /2 e −iH 1 /4 . Next we define the Q-matrix in the time frame α (= a, b) as Q α = j (|ψ + α j 〉〈ψ + α j |−|ψ − α j 〉〈ψ − α j |). The right and left Floquet eigenvectors |ψ ± α j 〉 and 〈ψ ± α j | satisfy the eigenvalue equations U α |ψ ± α j 〉 = e −i(±E j ) |ψ ± α j 〉 and 〈ψ ± α j |U α = 〈ψ ± α j |e −i(±E j ) . An OBWN for U α is then defined as ν α = Tr B (Γ Q α [Q α , X ])/L B . Here Γ is the chiral symmetry (CS) operator. X is the unit cell position operator. For a system with L lattice sites, we decompose it into a bulk region and two edge regions at the left and right. The trace Tr B is taken over the bulk region, which contains L B lattice sites. The length of each edge region is L E = (L− L B )/2, which should be chosen properly in order to avoid the obstruction of non-Hermitian skin effect. Finally, we define two OBWNs for a 1D non-Hermitian FTI with CS as ν 0 = (ν a +ν b )/2 and ν π = (ν a −ν b )/2. These winding numbers can only take integer values. They are further related to the numbers of Floquet edge modes at zero and π quasienergies n 0 and n π through the relations (n 0 , n π ) = 2(|ν 0 |, |ν π |). Following our analysis in the main text, (ν 0 , ν π ) could also count the numbers of fractional-quasienergy edge modes in the qth-root descendants of the parent model U.

B Topological invariants of the non-Hermitian FSOTI
Here we summarize the definition of bulk winding numbers for the second parent model in Sec. 3 of main text. Following Ref. [62], we first transform the Floquet operator U into two symmetric time frames a and b by shifting the initial time of evolution from t = 0 to t = 1/4 and t = 3/4 respectively. The Floquet operators in these time frames take the forms U a = e −iH 2 /4 e −iH 1 /2 e −iH 2 /4 and U b = e −iH 1 /4 e −iH 2 /2 e −iH 1 /4 . Performing Fourier transforms from position to momentum representations, we obtain U α = k x ,k y |k x , k y 〉U α (k x , k y )〈k x , k y | with α = a, b. In the tensor product The H x (k x ), H y1 (k y ) and H y2 (k y ) are Fourier transforms of the Eqs. (22)- (24) in the main text. U α (k x , k y ) has the CS in the sense that Γ U α (k x , k y )Γ = U −1 α (k x , k y ) for α = a, b, where Γ = σ z ⊗σ y . In our model, U 0 simply describes the evolution operator of a Su-Schrieffer-Heeger model in its topological flat band limit, which possesses a winding number w = 1. Taking the Taylor expansion of U α (k y ) yields U α (k y ) = cos(E) − i(d αx σ x + d αz σ z ), for which another winding number can be defined as . Put together, we obtain a pair of winding numbers (ν a , ν b ) = w × (w a , w b ) for the Floquet operators (U a , U b ). Their combination results in the integer topological invariants ν 0 = (ν a + ν b )/2 and ν π = (ν a − ν b )/2 of the twodimensional parent system U, which are related to the numbers of Floquet corner modes at zero and π quasienergies n 0 and n π through the relations (n 0 , n π ) = 4(|ν 0 |, |ν π |) [62]. Following the analysis in the main text, (ν 0 , ν π ) could also determine the numbers of fractional-quasienergy corner modes in the qth-root descendants of the parent model U so long as the chiral symmetry Γ is preserved.

C Stability to disorder
In this section, we demonstrate the stability of qth-root Floquet topological phases to disorder through numerical calculations. For the non-Hermitian FTI, we add random intracell coupling terms H 1d = n W n |n〉〈n| ⊗ σ y to H 1 and H 2d = n W n |n〉〈n| ⊗ σ x to H 2 in Eqs. (19) and (20), respectively. Here W n , W n take different values for different unit cells n and vary randomly in the range of [−W, W ]. The form of disorder terms H 1d and H 2d are chosen to be general enough and also to ensure that the chiral symmetries of the parent and the qth-root systems are preserved. In Fig. 8, we show the quasienergy spectrum and gap functions of the qth-root non-Hermitian FTI for q = 2, 3, 4 with the disorder amplitude W = 0.2 (comparable with the minimal energy scale of the clean system). The results show that the degenerate edge modes at different fractional quasienergies in the bulk spectrum gaps are well preserved under the impact of disorder. In Fig. 9, we further show the E = π/2, π/3, 2π/3, π/4, 2π/4, 3π/4 edge modes and their spatial profiles for q = 2, 3, 4 with the same disorder amplitude W = 0.2. It is clear that these fractional quasienergy edge modes indeed survive in the disordered system and are well localized around the sample boundary. The results presented in Figs. 8 and 9 are obtained for one disorder realization.   Figure 9: The real part of Floquet spectrum and fractional-quasienergy edge modes of the qth-root non-Hermitian FTIs for q = 2, 3, 4 with disorder. The notations, length of lattice and system parameters are the same as those used in Fig. 4.
We also checked a number of other disorder realizations in the calculation and found no observable difference.
For the 2D model, we introduce disorder by adding H 1d = x ⊗ n W n |n〉〈n|⊗σ z and H 2d = x ⊗ n W n |n〉〈n|⊗σ x to H y1 and H y2 in Eqs. (23) and (24), respectively. The W n and W n take different values for different cell indices n and vary randomly in the range of [−W, W]. We also choose the disorder terms in H 1d and H 2d to be general enough and to make sure that the chiral symmetries of the parent and the qth-root models are retained. In Figs. 10 and 11, we present the gap functions and the spatial profiles of fractional-quasienergy corner modes of the qth-root non-Hermitian FSOTIs for q = 2, 3 and the disorder amplitude W = 0.2 (comparable with the minimal energy scale of the clean system). The numerical results clearly suggest that the qth-root second order topological phases and their accompanying corner states in our system are robust to perturbations induced by symmetry-preserving disorder. The results presented in Figs. 10 and 11 are obtained for one disorder realization. We also checked a number of other disorder realizations in the calculation and found no observable difference.