Engineering and probing non-Abelian chiral spin liquids using periodically driven ultracold atoms

We propose a scheme to implement Kitaev's honeycomb model with cold atoms, based on a periodic (Floquet) drive, in view of realizing and probing non-Abelian chiral spin liquids using quantum simulators. We derive the effective Hamiltonian to leading order in the inverse-frequency expansion, and show that the drive opens up a topological gap in the spectrum without mixing the effective Majorana and vortex degrees of freedom. We address the challenge of probing the physics of Majorana fermions, while having only access to the original composite spin degrees of freedom. Specifically, we propose to detect the properties of the chiral spin liquid phase using gap spectroscopy and edge quenches in the presence of the Floquet drive. The resulting chiral edge signal, which relates to the thermal Hall effect associated with neutral Majorana currents, is found to be robust for realistically-prepared states. By combining strong interactions with Floquet engineering, our work paves the way for future studies of non-Abelian excitations and quantized thermal transport using quantum simulators.


I. INTRODUCTION
Quantum simulation is emerging as one of the most promising pillars of quantum technology, as it enables the observation of phenomena predicted to occur in models that are difficult to encounter in nature. A central ingredient for emulating the behavior of quantum systems is the ability to engineer the underlying Hamiltonian and probe its physics. Over the last decade, a toolbox was developed based on periodic (Floquet) drives, with the aim of imprinting novel physical properties by dressing the states of static systems [1][2][3].
In this work, we propose a protocol to engineer the timereversal-broken Kitaev honeycomb model using periodicallydriven ultracold atoms. We further address the significant challenge of probing the topological properties associated with neutral, particle-nonconserving Majorana modes in quantum simulators, with only limited access to the spin degrees of freedom. To this end, we propose nonequilibrium protocols combining slow ramps with abrupt quenches in the presence of the Floquet evolution, to reveal the chiral spin liquid phase, and measure the edge heat currents as its hallmark signature. We close by discussing state preparation, opportunities for a physical implementation, potential challenges and outlooks.

II. FLOQUET ENGINEERING THE KITAEV HONEYCOMB MODEL
Consider a periodic drive that implements the following Floquet unitary (i.e., the time-evolution operator over one period T of the driving sequence): where T =2π/ω D defines the relation between the drive period and frequency; see Fig. 1(d). Physical implementations of this driving sequence, using quantum simulators, are discussed in Sec. IV. A straightforward application of the van Vleck inversefrequency expansion [1][2][3] shows that, in the high-frequency limit, the dynamics is stroboscopically generated by the effective Hamiltonian H eff and kick operator K eff , according to with [see Appendix A] Here J α =−J α /3, and [ijk] αβγ = i, j α j, k γ denotes three neighboring sites where all associated spin operators are different: α =β =γ; see Fig. 6. We point out that the kick operator is crucial for capturing the correct stroboscopic dynamics [1,2,53]. While the infinite-frequency contribution to H eff is readily recognizable as the Kitaev model [first term in Eq. (4)], the (1/ω D )-correction opens up a chiral topological gap ∆ between the two Majorana bands (at the Dirac point); see Fig. 1(b) and Fig. 2(a). Here, the gap magnitude is ∆=3 √ 3g/4, with g=πJ 2 /(3ω D ) for isotropic systems (J α =J). Note that time-reversal symmetry is broken by the periodic drive; however, unlike applying an external magnetic field [29], one can verify that vortices remain dispersionless under this Floquet drive to any order in the drive frequency [ Fig. 1(b)], since the plaquette operators [29] commute with the generators of the three unitaries in Eq. (2). Hence, U F preserves the vortex structure of the initial state, see Appendix B.
Related Floquet-Kitaev physics has recently been explored in Josephson junction arrays [54], and in solid state systems [55,56], setting the focus on anomalous topological phases of matter [57][58][59], and the characterization of edge modes and their properties [60].
In the next Sec. III, we place the emphasis on the design of practical probes tailored for cold-atom simulators. Specifically, we first demonstrate the existence of an energy gap in the spectrum of H eff , using a practical spectroscopic approach. Then, we present a nonequilibrium protocol to probe the Majorana chiral edge transport, based on local spin-spincorrelation measurements. Finally, we discuss a state preparation sequence that allows us to demonstrate the topological properties of the system starting from a realistically-prepared state. We show numerical simulations of these protocols in the presence of the Floquet driving sequence that generates the model.

III. PROBING FLOQUET-KITAEV PHYSICS
Let us start by listing a number of outstanding issues, which are related to present-day limitations in measuring and prob-ing chiral spin liquids using cold-atom simulators. (i) The Kitaev honeycomb model appears deceptively similar to the emblematic Haldane model of Chern insulators [61][62][63][64], when written in terms of Majorana fermions [29]; however, it is substantially more challenging to probe and analyze, since Majorana excitations are fractionalized excitations which do not conserve the number of particles. In particular, it is not possible to access the topological properties of this model through mass/particle transport measurements and traditional spectroscopic probes. (ii) Majorana particles emerge as an effective degree of freedom: in contrast, the fundamental degree of freedom, accessible on the quantum simulator, is the spin. Thus, the question arises as to which quantity to measure in order to observe the topological properties associated with the effective (Floquet-Kitaev) Hamiltonian. (iii) A different, yet equally important, issue concerns the state preparation: can one have access to topological signatures from a realisticallyprepared state? Last, (iv), we emphasize the critical role of the Floquet drive itself, which ensures that the vortex and Majorana degrees of freedom remain decoupled throughout the dynamics; this important feature of our scheme is instrumental for accessing the physics of chiral spin liquids in the quantumsimulation framework.
To exclude extraneous features associated with a non-ideal initial state, we begin our analysis in Secs. III A and III B by using the rigorously determined ground state of the system as the initial state for our proposed Floquet protocol. Subsequently, in Sec. III C, we examine the preparation of the initial state and its influence on the system's response.

A. Gap spectroscopy
We first discuss the detection of the bulk gap in the spectrum of the driven spin system [Eq. (2)], which opens as a consequence of the (chiral) periodic drive; see Fig. 2(b) and Section II. Since the Majorana fermions are emergent degrees of freedom, probing their dispersion directly poses a significant challenge for present-day quantum simulators. The difficulty arises from the presence of vortices to which an arbitrary probe would also generically couple. Here, we address this question by showing how a practical spectroscopic probe, applied directly to the spin degrees of freedom, can accurately detect the gap opening.
As we show below, a spectroscopic probe can be simply realized through a sinusoidal modulation of the (uniform) coupling between the z-bonds, where ω p and A p are the probe's frequency and amplitude, respectively. In our simulations, we find that this protocol works properly in the time-scale separated regime ω −1 D ∝∆∼ω p ω D . Figure 2c (inset) shows the number of excited particles as a function of time, N exc (t), starting from the ground state |G and applying simultaneously the Floquet drive and the spectroscopic probe [Eq. (5)]. This quantity is proportional to the spectroscopic signal that one would detect in an experimental context. The main panel shows the dependence on the probe frequency ω p at the measurement time t * , chosen such that the resonant signal is detectable; in practice, this amounts to Jt * ∼10 2 . One observes clearly discernible resonant peaks corresponding to transition frequencies ω res between the ground and excited many-body states, with a peak positioned at 2∆. The physical origin of this factor 2 is intimately related to the fractional character of the Majorana quasiparticles: it can be traced back to the restriction that excitations can only be created in pairs in the Kitaev model [29]. A rigorous proof for the appearance of this factor 2 can be found in Appendix C.
In Fig. 2(d), we display the resonant frequency for different values of ω −1 D . Floquet theory [solid black line] predicts a gap closing in the infinite-frequency limit; see Eq. (4). The extracted gap clearly follows this law for large enough systems, with deviations due to higher-order frequency correc-tions discernible at large values of ω −1 D . The numerical data shows a reasonable agreement for different system sizes, both for cylindrical (circles) and planar (crosses) geometries. In very small systems with boundaries (planar geometry), we note that finite-size effects can cause ambiguity in identifying the peak corresponding to the bulk gap, see Appendix E. In Fig. 2(d), we explicitly compare different system sizes in order to help identify a suitable regime for a potential experimental realization.

B. Probing chiral Majorana edge transport
The detection of the bulk gap opening is insufficient to demonstrate the topological nature of the ground state associated with the effective time-reversal broken Kitaev Hamiltonian H eff . In chiral spin liquids, a natural signature of topological (non-Abelian) excitations is provided by the quantized thermal Hall effect [30,65]. This phenomenon, which arises from the chiral transport of heat along the edge of the system, was recently measured in "Kitaev" materials [32,34,38]. Motivated by such energy transport measurements [30,65], we propose a non-equilibrium probe to detect the chiral transport of heat along the edge [ Fig. 3a], and hence to reveal the chiral spin liquid nature of our Floquet-engineered system.
Let us consider a system with boundaries prepared in the ground state of the isotropic effective Hamiltonian H eff (i.e. J α =J). We now perform a quench by adding a local perturbation on a single z-bond, i, j z z , located on the edge of the system: H z →H z +J q S z i S z jz ; see Fig. 3a and Fig. 10. The modified Floquet unitary, which includes this local perturbation, is denoted byŨ F . We let the system evolve under the quenched Floquet drive for a fixed number of cycles ∈{1, 2, . . . , N T }, and we measure the local excess of energy in every unit cell located on an edge: where the unit cell m is located on the edge, and i, j x x , i, j y y , i, j z z denote the sets of links connected to m; see Fig. 10. Measuring E m ( ) stroboscopically at each drive cycle, ∈{1, 2, . . . , N T }, gives rise to a timetrace of data. The chiral nature of the edge dynamics can be revealed by post-processing this time-trace using a discrete Fourier transform [66], where L y is the length of the edge considered. Indeed, we verified that this Fourier spectrum clearly captures the chiral dispersion relation of the Majorana edge modes; see Appendix F.
A closer investigation reveals that the chiral signal is already detectable in the Fourier spectrum of the local zzcorrelations on the edge: where the edge z-bond i, j z is shown in Fig. 10. This quantity is a constituent part of the local energy E m in Eq. (6);  however, it requires fewer measurements. Therefore, for the sake of experimental simplicity, we focus our study on the behavior of C m . We note that this detection protocol requires unitary evolution (set by the modified Floquet sequence) and a single final measurement, which is particularly suitable for analog quantum simulators. In particular, it does not require the application of additional spin operators to create excitations on the edge of the sample, and whose effects have to be subsequently measured. Figure 3(b),(d) shows the Fourier spectrum |A(k, ω)| extracted from C m ( ), as computed from a numerical simulation of the nonequilibrium (Floquet) quench, for two system sizes. In these simulations, we consider a cylindrical geometry, and we apply the quench protocol on opposite edges of the cylinder. We observe a clear signal at low frequencies, which corresponds to the excitation of the chiral edge mode.
Importantly, the slope of the chiral signal in (ω, k) space corresponds to twice the group velocity v edge of the Majorana edge modes [67]; see Appendix C. We show the dependence of the extracted velocity on the Floquet drive frequency ω D in Fig. 3(c); see Appendix I for details on the fitting proce-dure. In particular, we verify that the chiral signal disappears (i.e., the edge mode's velocity tends to zero) in the infinitefrequency limit, where the topological gap closes; see Eq. (4). Moreover, we have also verified that the chiral edge transport detected by our probe vanishes upon transitioning from the chiral spin liquid to the non-chiral phase of the phase diagram [ Fig. 1(c)]. This confirms that the chiral edge signal is a genuine signature of the chiral spin liquid phase, and not a mere consequence of the circular "polarization" of the drive, cf Fig. 12.
Interestingly, one can use this measurement as an alternative to the spectroscopic measurement discussed in Sec. III A: indeed, one can estimate the size of the topological bulk gap from the extracted velocity, using the relation v edge =3∆/(2π).
Finally, to demonstrate a direct correspondence between the chirality of the edge mode and that of the Floquet drive, we reverse the "polarization", , of the periodic driving sequence, and observe a reversal in the propagation direction of the excitations in the Fourier spectrum; compare the top and bottom panels in Fig. 3(b) and (d). A close examination shows that |A (k, ω)| =|A (k, −ω)|; this implies that the Floquet quench excites a small fraction of non-chiral modes. Since these undesired excitations grow when increasing the quench strength J q , it is preferable to work in the weak-quench regime J q /J 1.
While the proposed measurement protocol is intimately related to the thermal Hall effect associated with non-Abelian excitations in chiral spin liquids [30,65], detecting the quantization of the thermal Hall conductivity remains an outstanding challenge in the context of quantum simulation; see Refs [68][69][70][71] regarding possible heat-transport measurements in quantum gases.

C. State preparation sequence
So far, we have shown how to detect hallmark features of the chiral spin liquid by assuming that the system can be initialized in the vortex-free ground state. However, preparing topological many-body states (which are entangled over long distances) from product states is a highly nontrivial task per se [72][73][74][75][76]. In particular, adiabatic quantum state preparation would require going through a topological phase transition, which is associated with a gap closing point in the thermodynamic limit. In quantum simulators, one can nevertheless exploit the finite size of the system (and hence, the possible existence of a finite-size gap at the topological phase transition) to perform adiabatic quantum state preparation [7,77,78].
To address this relevant aspect in detail, we build on existing ideas for preparing the ground state of the toric code [79][80][81]. At present, this is feasible using digital quantum simulators and requires the appropriate sequential application of high-fidelity Hardamard and CNOT gates [13], or entangling gates combined with mid-circuit measurements [82]; to the best of our knowledge, the optimal state preparation of the toric code ground state is currently an open problem in analog quantum simulators [83].
The toric code ground state describes the low-energy physics of our system in the close vicinity of the point J x,y =0 [29]; see Fig. 1(c). Thus, starting from this (Abelian topological) gapped state, we can gradually cross the critical point [77] and end up in the desired (non-Abelian) gapped state, in the presence of the Floquet drive. This amounts to slowly ramping-up the couplings: over a time t ramp T , while keeping J z =J fixed. We note that the ramp terminates at the isotropic point J α =J, and that the bulk gap closes at the critical point [ Fig. 4a], which inevitably causes excitations. Nevertheless, we find that the latter are not detrimental to the chiral edge signal introduced in Sec. III B for realistic finite system sizes, as we now discuss. We apply the ramp (8) in the presence of the evolution generated by U F (J α (t)) for a duration t ramp , and then perform the quench U F →Ũ F following the protocol described in Sec. III B. Figure 4(b) shows the corresponding Fourier spectrum A(k, ω), which displays a noticeable chiral signal. Besides, one clearly observes additional (undesired) bulk excitations, which are created due to the finite ramp duration.
We point out that, even in the ideal "adiabatic" limit where t ramp →∞, a long observation time t obs =N T T is required for accurately extracting the chiral signal from the Fourier spectrum. This is illustrated in Fig. 4(c), which shows the convergence of the extracted edge mode velocity as one increases the observation time.
Moreover, we emphasize that the finite duration of coldatom experiments limits the duration of the ramp t ramp used to prepare the state of interest. Figure 4(d) shows that a ramp duration of the order of Jt ramp ≈150 suffices to clearly detect the chiral signal (v edge =0), as we illustrate by comparing the results obtained for Floquet sequences of opposite chirality [circle vs. diamond markers]. We verify that the edge mode velocity decreases when increasing the drive frequency, v edge ∝ω −1 D , in agreement with theory; see the grey datapoints in Fig. 4(d) and Fig. 16.
We stress that the ramp time required for adiabatic quantum state preparation necessarily diverges in the thermodynamic limit. Nevertheless, our simulations demonstrate that clear signatures of the chiral spin liquid should be detectable in realistic system sizes and within reasonable time scales.

IV. PHYSICAL IMPLEMENTATION
Having discussed the detection of topological signatures associated with the effective Hamiltonian H eff , let us now turn to the details of its physical implementation. We argue that it is possible to either use ultracold fermions (with Hubbard interactions), dipolar interactions between molecules or Rydberg-atom-based spin systems for this purpose.
A first cold-atom realization of the Kitaev model (not relying on high-frequency periodic drives) was proposed in Ref. [84]. However, this proposal remains challenging to implement in practice, since it requires a large number of laser beams and a two-photon Raman coupling. Moreover, it remains unclear whether it can be generalized to open up the chiral topological gap, while preserving the vortex configuration.

A. Implementation of the Floquet drive
Consider the spinful Fermi-Hubbard model on a honeycomb lattice: where J σ denotes the strength of nearest-neighbor speciesdependent hopping, and U is the onsite interaction strength.
The fermionic operators obey the anti-commutation relation {c iσ , c † jσ }=δ ij δ σσ . Recall that, in the strongly-interacting atomic limit, J σ U , double occupancies are suppressed, and the spectral function exhibits gaps of order U , separating the so-called Hubbard bands. In this regime, the system is described by the Heisenberg model [84][85][86][87]: where J zz =−2(J 2 ↑ +J 2 ↓ )/U is the Ising interaction, and J ⊥ =−4J ↑ J ↓ /U is the exchange interaction. The spin operators are defined in terms of the fermion operators as Freezing one of the spin species (i.e., suppressing maximally its hopping matrix element, e.g., J ↓ =0) inhibits the exchange interactions, J ⊥ =0, and leaves only the global nearest-neighbor Ising z-interaction J zz . Alternatively, an Ising-only interaction can be achieved by applying a magnetic field gradient to suppress the exchange interactions, and then turning on a weak resonant periodic drive to enable the z-interaction term [88]; we remark that the latter is a resonant nonequilibrium scheme, which may produce unwanted heating, when combined with the primary Floquet drive from Eq. (2). Recently, it was also shown that a controllable anisotropy in the Heisenberg couplings can be realized using a system of Rydberg arrays [89], or by means of Floquet engineering in superconducting qubits [90] and ultracold molecules [91].
With this idea in mind, the basic steps behind the fermionic implementation of the Kitaev model can be summarized as follows: First, we eliminate the spin-exchange term, which ensures that each bond interacts along a single spin direction (Ising-z). Second, we apply a periodic three-step modulation of the hopping matrix element J ↑ (t) to isolate the different spatial bond types (x, y, and z) in time. Finally, strong singleparticle π/2-pulses in the spin channel, applied between the steps of the periodic drive, rotate the z-interactions into the desired type (x or y). Applying the periodic step drive at a moderately high drive frequency results in the Hamiltonian from Eq. (4).
Let us now explain the procedure in more detail. In the first step, a time-periodic spatial modulation of the non-zero hopping matrix element, J ↑ (t+T )=J ↑ (t), can be used to single out one type of bonds per Floquet step [see Fig. 1 This can be achieved using a drive similar to Ref. [92]. In this scheme the strength of the tunnel coupling along the three different directions is modulated as a function of time by changing the intensities of the individual laser beams that form the underlying honeycomb lattice. The main challenge for experimental realizations will be to implement such a scheme  12). In the high-frequency limit, J ↑ U ωD, applying the inversefrequency expansion (i)→(ii) leads to a modified static fermion Hamiltonian with nearest-neighbor hopping terms (also in the spin channel) of strength J ↑ /3. A Schrieffer-Wolff transformation (iii) then gives the effective spin-1/2 Hamiltonian H eff from Eq. (4) with interaction strength J=2J 2 ↑ /(9U ), which governs the physics of the lower Hubbard band stroboscopically. Alternatively, (i)→(iv), in the stronglyinteracting limit J ↑ ωD U , we first apply the Schrieffer-Wolff transformation to arrive at the spin-1/2 Hamiltonian (13), with periodic in time, spatially oscillating z-interaction Jzz,ij(t)=2J 2 ↑,ij (t)/U . A subsequent application of the inverse-frequency expansion then gives rise to H eff with interaction strength J=2J 2 ↑ /(3U ).
Step (iv) can also be independently taken as the starting point for the implementation in a spin-1/2 quantum simulator (red shaded oval). with state-dependent optical potentials, that enable a statedependent manipulation of the tunnel couplings. In principle, such a drive can also be implemented with quantum gas microscopes [93] using tightly-focused optical tweezer beams, whose position and strength can be controlled dynamically, as recently demonstrated [94,95]. Again, the main challenge lies in the realization of the required state dependence. The interaction strength U is kept constant on all bonds throughout the drive. This periodic modulation gives rise to a piece-wise constant periodic Hubbard Hamiltonian H Hubbard (t)=(H z Hubbard , H y Hubbard , H x Hubbard ); the superscript α=x, y, z here refers to the spatial bonds in the honeycomb lattice [see Fig. 5b (i)]. In the atomic limit, J ↑ U , all terms H α Hubbard give rise to a nearest-neighbor Ising zinteraction on the x, y, z-links, respectively.
At this stage, in the limit J ↑ U , the system is effectively described by a spin-1/2 Ising model, whose Ising interaction is periodically modulated to switch between the different bonds on the honeycomb lattice [see Fig. 5b (iv)]: (13) Incidentally, Eq. (13) can also be taken as the starting point for a spin-1/2 implementation [ Fig. 5(b), red oval], e.g., using Rydberg simulators [96] operated in the nearest-neighbor interaction regime [97].
To generate the Floquet unitary in Eq. (2), it remains to apply global on-site π/2-rotations in the spin channel; if they occur at a time scale much faster than the inverse interaction strength J zz , this will consecutively rotate the zinto yand x-interactions; see Fig. 5(a). Due to the Euler angles theorem, it suffices to implement π/2 pulses along two out of the three spin axes.
Notice that the discussion so far tacitly assumed the limit J ↑ ω D U , which places the drive frequency in between the two lowest Hubbard bands. Hence, the above analysis is equivalent to first applying a Shrieffer-Wolff transformation (SWT) [to eliminate the largest energy scale U , see Fig. 5b (i)→(iv)], followed by the inverse-frequency expansion (IFE) [second-largest scale ω D ]; see Fig. 5b (iv)→(v). The effective Ising interaction is then given by J z =J zz /3=−2J 2 ↑ /(3U ). By contrast, in the high-frequency limit, J ↑ U ω D , we have to first apply the IFE [ Fig. 5b (i)→(ii)] and then the SWT, which leads to J z =J zz /9=−2J 2 ↑ /(9U ) [Fig. 5b (ii)→(iii)]. Note that the Ising interaction strength J z is larger [98] -by a factor of three -in the regime J ↑ ω D U , due to a dressing by the Floquet drive, as compared to the high-frequency regime. We note that these effective couplings are related to those from Eq. (2) via J z =−J z /3, etc.
Setting ζ=ω D /U , the two limits can be nicely unified by the more general expression, valid away from resonances U =lω D (l∈N): which is obtained by reconciling the SWT and the IFE into a single framework [99], see Appendix A. Figure 5(c) shows that the stroboscopic magnetization dynamics of the Fermi-Hubbard model on a single honeycomb cell, following the above Floquet protocol [solid lines], agrees well with the dynamics of H eff [dotted lines], for a wide range of non-resonant frequencies [100]. Notice the proper, ζ-renormalized scaling of the x-axis, which is essential for capturing the correct effective z-interaction across two orders of magnitude in the drive frequency (top to bottom panels).
Finally, we note that the strongly-interacting regime is particularly appealing for Floquet engineering, since the condition J ↑ U readily ensures the high-frequency limit w.r.t. the Majorana bandwidth of the Kitaev model: The anisotropy in the spin couplings can be controlled by setting the durations of the corresponding periodic pulses, respectively. Thus, in platforms where modulating J z periodically in time is infeasible, the spectroscopic drive from Sec. III A can be effectively implemented by slowly changing in time the duration of only the z-bonds pulselength in the Floquet unitary (2). The same idea can be used to implement the slow ramp in Sec. III C; this will also enable the study the phase transition between phases A and B in the time-reversal broken Kitaev model. Finally, note that to realize the local quench on selected z-bonds from Sec. III B, one may use a quantum gas microscope [93,101,102] to imprint the perturbation; alternatively, if the atoms are captured in tweezers with controllable positions [95], one could change the bond distance locally, which will affect the corresponding interaction strength.

C. Potential challenges
The most prominent challenge for a potential experimental realization arises from the need to stabilize the global onsite spin rotations out to a large number of driving cycles. However, this could be potentially solved by adding spinecho pulses to cancel dephasing, as recently demonstrated with molecules [91]. A second issue may occur due to the necessity to keep the temperature of the state at the order of, or below, the topological gap, which is only a fraction of the superexchange energy scale. Ultracold fermionic systems have recently reached low enough temperatures to observe antiferromagnetism [103][104][105] in the half-filled Fermi-Hubbard model. However, to reach low enough temperatures in order to observe ordered states, even lower temperatures are needed. The expected energy scales are of similar order as compared to the ones that appear in the current work [106]. There are several promising techniques that have been proposed in order to reach these temperature scales [107,108]. and we anticipate that despite this significant challenge cooling techniques will continue to improve. This will be beneficial for a number of fermionic quantum simulation experiments. Moreover, as we have shown above, detecting the chiral energy currents using nonequilibrium probes can still be robust to a certain amount of excitations. Finally, the gap size is controlled by the drive frequency; we observed numerically that ω D ≈8J is about the smallest drive frequency for which resonances within the Majorana bands are still suppressed.

V. CONCLUDING REMARKS
An important issue inherent to the Floquet engineering of strongly-interacting systems concerns detrimental heating processes [81]. While a comprehensive study of energy absorption and entropy creation in this nonequilibrium setting is beyond the scope of this work, we expect thermalization to infinite-temperature in the driven spin system to be absent for frequencies above the single-particle Majorana bandwidth: this follows from the conservation of the plaquette quantum numbers under the Floquet drive, which renders the Floquet Hamiltonian effectively single-particle, and makes the present model particularly appealing and promising for quantum simulation. That said, secondary heating effects can still be caused by the presence of resonant processes with higher bands, or due to the creation of doubly occupied sites in the fermionic realization.
Looking forward, two exciting directions concern (i) the investigation of local operations that would allow for the braiding of vortices and Majorana degrees of freedom (both in the bulk and on the edge) [25][26][27]109]; having dispersionless (i.e., immobile) vortices at our disposal is expected to prove instrumental in this context; (ii) the design of probes in view of demonstrating the half-quantization of the thermal Hall conductivity associated with chiral edge transport. Going beyond the present effective Kitaev Hamiltonian H eff , we note that it is well within the scope of our Floquet engineering protocol to incorporate long-range interactions, or even add an additional Heisenberg term to the effective Hamiltonian; this will enable the quantum simulation of a larger class of Hamiltonians describing Kitaev materials. Thus, by analyzing genuinely nonequilibrium protocols for probing and engineering topological Hamiltonians, our work paves the way for investigating Floquet topological phases of matter without static counterparts [58,59] using quantum simulators, and demonstrates the possibility to reconcile Floquet engineering protocols with strong inter-actions to study strongly-correlated quantum phases of matter.
Note added: While finalizing this manuscript, we became aware of a similar proposal based on Rydberg digital quantum simulation [96].
Acknowledgments: We acknowledge discussions with C. Hickey where H eff is the time-independent Floquet Hamiltonian and K eff (t) = K eff (t + T ) is the time-periodic generator of micromotion, also known as the kick operator. In the van Vleck IFE, K eff satisfies the boundary condition In this section, we will use the primed interaction strength J α , which is related to the one used in the main text by J α = −J α /3. The Floquet unitary is generated by the time-periodic Hamiltonian with the piece-wise constant, time-periodic three-step function We can expand the effective Hamiltonian and kick operator to leading-order in the inverse-frequency ω D , using the weighted time-ordered integrals: where we fix the convention (x mod 2)∈[−1, 1). The righthand side equations make use of the Fourier expansion of the Hamiltonian, given by H(t) = ∈Z e +i ωt H . For the spin-1/2 implementation from Eq. (2), we readily obtain (with J = −J /3): Fig. 6] (A7)    The data show clearly that all corrections to the desired order are captured by Eq. (4). Note that the transformation generated by the kick operator is crucial for the validity of the test. The system size is a single honeycomb cell, and all interactions are isotropic J α = J .
with the piece-wise linear T -periodic function The three-body terms in H (1) eff encompassed by the notation [ijk] αβγ are shown in Fig. 6. In turn, by investigating the scaling with ω D , in Fig. 7 we show that the above expressions for the effective Hamiltonian and kick operator are complete.
In the fermion implementation, additional terms can arise due to higher-order terms of the Schrieffer-Wolff transformation, of order O(U −2 ), O(U −1 ω −1 D ), O(ω −2 D ), etc. We leave their investigation to a future study.
Instead, here we focus on the derivation of the effective interaction-renormalized Ising coupling, cf. Eq. (14). To this end, we apply the generalized Schrieffer-Wolff transformation (SWT): using similar arguments to the derivation of Eq. (31) in the Supplemental Material to Ref. [99], we pretend that the drive is resonant with the interaction, U = lω D (l ∈ N), and do analytic continuation away from the resonance in the end. In particular, the effective Ising coupling in the generalized SWT is given by where c are the Fourier coefficients of the Floquet drive c(t), and we used the relation l = U/ω D a few times. Notice that the last equality is only valid sufficiently far away from resonance, i.e., for U = lω D . Using the Fourier coefficients for the periodic three-step drive, gives the final expression (A11) In the infinite-frequency limit, U ω D , only the = 0 harmonic survives, which leads to J z = −2J 2 ↑ /(9U ). By contrast, in the strongly-interacting limit, ω D U , all harmonics contribute equally; using the relation ∞ =−∞ sin( π/3) π 2 =1/3, then leads to J z = −2J 2 ↑ /(3U ). Since the derivation in Eq. (A9) is, to a certain degree, ambiguous, and given the heuristic character of the analytic continuation, we performed numerical tests, shown in Fig. 5c, In this Appendix, we provide more details on the representation of the Kitaev honeycomb model in terms of Majorana fermions, which plays a central role in this work. We focus on the isotropic case J α = J, but the generalization to the anisotropic model is straightforward.
For the time-reversal broken Kitaev model, the effective Hamiltonian used in the main text is [29] By using the Jordan-Wigner transformation [110], one can transform the spins to Majorana fermions c andc, which obey {c i ,c j } = 0 and {c i , c j } = 2δ ij = {c i ,c j }. Then, the Hamiltonian becomes (see Refs. [30,110] and Fig. 8) Here, η b = ic jcjz , takes the value ±1 on each z bond b = jj z z , p denotes a single plaquette of the honeycomb lattice, b 1 and b 2 are the two z bond on the plaquette, and j denote the sites belonging to one sublattice ("filled circles"); see For the Kitaev model, it is well known that, for each plaquette p [see Fig. 8], there exists a local conserved quantity, the plaquette operator [29,110]: Here, σ α j = 2S α j are the Pauli matrices. The eigenvalues of W p are ±1 and under the Jordan-Wigner transformation, the plaquette operator can be written as [110] By convention, states with W p = 1 on all plaquettes are defined to have a vortex-(or flux-)free configuration [29].
Since [H, W p ] = 0, any state of the original many-body spin system can be written as a product state over the vortex (W ) and Majorana (M ) sectors. In particular, for the ground state, we have In the following, we discuss the two sectors separately.

Vortex sector
Let us start with the vortex sector. According to Lieb's theorem [29,111], in the thermodynamic limit the η-dependent many-body ground state of the original spin system, is given by the vortex-free field configuration |G W =+1 . For the planar geometry used in this work, η b eigenvalues defined on the same row should have the same sign, which can be defined as η r with r indicating the rth row (see Fig. 8). Then, since there are L y − 1 rows of η for a planar geometry with L y unit cells along the y-direction, the vortex-free configuration is actually highly degenerate, with the degeneracy being 2 Ly−1 . Nevertheless, these different η b configurations for a system with a vortex-free planar geometry are actually equal to each other because they can be transformed into the all η b = 1 configuration by the gauge transformation with j r denoting the row index for site j; see Fig. 8. Similarly, for the vortex-free cylinder system, there are 2 Ly different configurations of η r that keep the system vortex-free. Nevertheless, they can be divided into two different configurations, defined by the even or odd number of rows with η r = −1. Moreover, in the case of all η b = 1, if, instead of a periodic boundary condition on the cylinder, we use an anti-periodic boundary condition, we will obtain a configuration with only one row of η b = −1. This means that the two different η b configurations have the same topological properties in the thermodynamic limit.
Considering that our probe does not alter the η r configuration and the states with different η r configurations are orthogonal, the detected signals from different η r configurations do not interfere with each other. Therefore, we can fix the gauge, and focus on the case where η b = 1.

Majorana sector
Replacing η b by their eigenvalues, a system with 2N sites gives rise to a 2N × 2N skew-symmetric matrix A jj = −A j j , which defines a quadratic Hamiltonian describing the Majorana sector: The many-body Majorana ground state |G M can be constructed from the single-particle modes. For a 2N × 2N single-particle Majorana Hamiltonian H M , the creation operator for the i-th single-particle eigenmode can be written as Here, f i j is the jth expansion coefficient of the ith eigenmode. Then, the ground state in the Majorana sector is given by [30] with the eigenvalues corresponding to the eigenmodes γ † i (i = 1, . . . , N ) being negative.
Assuming translation invariance (all η b = 1), we can also derive the Majorana spectrum for a system with periodic boundary conditions in momentum space. To specify the unit cell notation in the Majorana representation, we represent the site index j of the Majorana operator c j , as the composite index (m, l), with l running over the two sites in the unit cell m. Following Ref. [30] .

(B13)
Here, the distance between two nearby sites sets the unit length, and q x and q y are the two component of the quasimomentum vector q. The dispersion relation reads as From this equation, the Dirac points can be shown to be located at q = ± 4π 3 √ 3 , 0 , with the energy gap ∆ = 3 √ 3g/4.
Hence, experiments on finite-size systems, aiming to detect the topological properties of the Majorana bands, should select the underlying lattice geometry carefully. For instance, to hit the Dirac point (where the Berry curvature of the Majorana bands is strongest) in small systems, the number of the unit cells in each direction should be a multiple of 3. In this section, we recall the single-particle formalism we use to treat the Majorana sector of the model. We also explain the origin of the additional conversion factor of 2 between the dynamically extracted data for (∆, v edge ) in the main text, and their actual (theory) values.
To time-evolve the system, consider again a generic quadratic Hamiltonian H M in the Majorana representation, Using the normalization factor 1/4, this operator fulfills the relation [29] [ where the commutator on the left-hand-side acts on Fock space, while the commutator on the right-hand-side -on the 2N × 2N dimensional space associated with the skewsymmetric matrix A. Then, the time evolution of Majorana modes in the Heisenberg picture can be calculated as with Q = e −tA .
We are now in a position to explain the origin of the conversion factor of 2 between the dynamically extracted gap and edge velocity, and their actual theoretical values. For our Majorana system, it follows from the relation Q = e −At that the dynamics is described by the eigenvalues of the matrix iA [instead of iA/4, as one may naively think due to Eq. (B7)]. This factor is important since it scales time in Eq. (C3). Let us now contrast this with the static properties of the Majorana system. For instance, to get the many-body gap, we use Eq. (B12), which states that the gap is determined by the eigenvalues of the matrix iA/2 [30]. Comparing the dynamics generated by iA, to the spectral properties obtained from iA/2, we find the additional factor of 2, required for the conversion between the two. The physical origin of this factor is intimately related to the fractional character of the Majorana quasiparticles: it can be traced back to the restriction that excitations can only be created in pairs in the Kitaev model [29]. Finally, let us mention that we also confirmed the existence of this factor in the full spin system, Eq. (4), by extracting the many-body gap from spectroscopy.

Appendix D: Floquet drive in the Majorana sector
In the presence of the Floquet drive, the time evolution operator over one Floquet period reads as where H α = − j,j α J α S α j S α j . Note that the plaquette operators W p [Eq. (B4)] commute with all H α , which means they also commute with U F . Hence, the Floquet drive does not mix the vortex and Majorana sectors, and we can numerically simulate the dynamics of the spin system in the Majorana sector.
Hence, we can restrict the analysis to the dispersive Majorana sector: which is obtained from the original spin Hamiltonian H α = − j,j α J α S α j S α j using the Jordan-Wigner transformation, cf. Sec. B.
Following Eq. (C3), the time evolution of a Majorana operator becomes Note that the order of unitaries in the definition of A eff is reversed compared to U F . With this definition of A eff , the Floquet Hamiltonian is given by There is no handy closed-form expression for the exact matrix A eff . However, for sufficiently high drive frequencies, an approximation can be constructed using the inverse-frequency expansion.
Appendix E: Details on the protocol for gap spectroscopy Here we discuss in detail the spectroscopy protocol in the presence of the Floquet drive. We then go on and show supplementary results to the main text. We prepare the system in the Majorana ground state of the isotropic effective Hamiltonian H eff , whose eigenmodes are annihilated by the operators γ i . We then evolve the system in the presence of the Floquet drive and simultaneously apply the spectroscopy protocol; this gives rise to the modified Floquet unitary where the z-bond strength changes at a frequency ω p ω D , according to: with A p the (dimensionless) probe strength, and -the Floquet cycle. Treating the spectroscopic probe as a constant during the cycle is justified by the scaling of the topological gap we want to detect: ω p ∼ ∆ ∝ ω −1 D . Thus, at times t n = nT (n ∈ N), the time-evolution operator is given by The spectroscopic signal is proportional to the number of excitation particles N exc , defined as where N is the number of occupied Majorana modes in the initial state (out of 2N modes in total), and γ i is defined in terms of the real-space Majorana operators in Eq. (B8). A straightforward application of Eq. (C3) leads to the expression where the 2N × 2N correlation matrix C has the matrix elements C ij = G M |c i c j |G M , and is an N -dimensional vector of time-evolved expansion coefficients, cf. Eq. (B8). Using Eq. (E5), we calculate the dependence of N exc on the probe frequency ω p for different geometries, cf. Fig. 9. For the planar geometry, we see a pronounced finite-size effect in small systems: for the 6 x × 6 y system shown in Fig. 9(a), the lowest resonant peak in the spectroscopy signal falls inside the gap due to the presence of boundary state excitations. Yet, for a large planar system, e.g., 18 x × 18 y , the influence of the boundary states becomes weaker, and the lowest resonance edge Figure 10. Schematic representation of the notation for the quench protocol: m labels a unit cell at the edge of the system; the site i belongs to the unit cell m, while i, jx x, i, jy y , and i, jz z indicate the bonds connected to i that are considered in the expression for the local unit cell energy Hm, cf. Eq. (6). peak shown in Fig. 9(b) gets close to 2∆ [see Fig. 2(d)]. On the other hand, the influence of the boundary states is less pronounced in the cylinder geometry. For the 12 x × 6 y and 18 x × 6 y systems shown in Fig. 9(c) and (d), both of their first resonant peaks are close to 2∆. Due to finite-size effect, the resonant peak in the larger system is closer to the theorypredicted gap value, when compared to the smaller system, especially for a system with a small gap [see also Fig. 2(d)]. From this equation, one finds that the expectation is just the number of particles (i.e., the occupation) in the m-th unit cell times the corresponding energy, which is a good description of the energy in the m-th unit cell.
To get the chiral edge-state transport, we take a system with a boundary, and prepare it in the ground state of H eff , cf. Eq. (B9). Then we apply a perturbation H p as a quench to the effective Hamiltonian, acting in the middle of the edge: with J q being the perturbation strength, and ij z denotes the z-bond in the middle of the edge. The corresponding Floquet unitary after the quench thus becomes Finally, we evolve the system withŨ F , and measure the local unit cell energy H m at the m-th unit cell at the edge: which can be evaluated in practice with the help of Eq. (C3). Similarly, in Sec. III B of the main text we defined the dynamics of the local zz-correlator, C m ( ), after the quench.
The Fourier spectrum of E m (t) can be obtained from where, L y is the number of unit cells along the y direction, N T is the number of Floquet cycles used to do the Fourier transformation, k j y = 2πj/L y , and ω = 2π /(T N T ). In the main text, to simplify notation, we used the quantities k and ω. We note that since E m ( ) is real-valued, we have the reflection symmetry |A(k, ω)| = |A(−k, −ω)|.
In our calculation of H m , considering that the strength g of the three-spin correction term is much smaller than J, we drop all terms proportional to g in Eq. (B2) to simplify the calculation. We show the obtained Fourier spectrum for the local unit cell energy in Fig. 11 (a) and (b), corresponding to drive frequency ω D /J = 8, 20, respectively. The data show a strong signal due to the excitation of the chiral edge modes, which agrees well with twice the theoretical group velocity v th edge = 3∆/(2π), as shown by the thin black dashed line in the figures. Note that in Sec. III B, we only showed results for the zz-correlation [Eq. (7)], instead of the local unit cell energy for the purpose of simplifying the experimental realization. Here, for completeness, we also show data for the zzcorrelation in Fig. 11 (c) and (d) [with ω D /J = 8, 50, respectively]. Comparing Fig. 11(a) and (c), one finds that the spectrum for the zz-correlation may not show all excitations, but it captures the major chiral signal of the local unit cell energy results; moreover, the signal is cleaner for the zz-correlations, which is advantageous in determining the value of v edge using a linear fit [see Sec. I]. Especially when the topological gap is small, e.g., in Fig. 11(d) where ω D /J = 50, the agreement of the chiral signal with the theoretical prediction for v th edge is particularly clear. However, the corresponding signal for the local unit cell energy is surrounded by additional excitations as shown by the Fig. 11(b).
To further demonstrate that the chiral signal arises from the properties of the chiral spin liquid phase B, instead of being a mere consequence of the circular "polarization" of the Floquet drive, we display the Fourier spectrum of the boundary zzcorrelation at the phase transition point, i.e., J x = 0.5J = J y , and J z = J in Fig. 12(a), as well as for one point in the bulk of the non-chiral toric code phase A z (J x = 0.25J = J y , and J z = 1), see Fig. 12(b). As expected, no chiral signal is observable outside the chiral spin liquid phase.
Finally, we point out that we use a relatively small value J q /J = 0.008 in the numerical simulations in order to avoid unwanted additional excitations that may interfere with the process of extracting the edge-state velocity. With respect to the feasibility of a potential experimental realization, it is important to note that this value can be increased. Let us briefly examine the role of the quench strength J q /J. The obtained results, shown in Fig. 13, reveal that increasing J q /J introduces additional "noise" to the Fourier spectrum, but the dominant chiral signal and the extracted edge-state velocity v edge remain qualitatively unchanged even for values up to J q /J 0.6. where the notation is the same as in Fig. 10, and ρ(t) denotes the density matrix at time t. At t = 0, the latter can be obtained from the ground state Eq. (B9). Using the Fourier transform the zz-correlation becomes Here, N y is the number of the sites in the periodic direction, and δ z = r m+1 − r m is the distance between two nearby zbonds at the edge.
In the interaction representation, ρ I (t) = e iH0t ρ(0)e −iH0t , we have where, H 0 = H eff is given in Eq. (D4), and H is the perturbation Hamiltonian with interaction picture representation H I (t) = e iH0t H e −iH0t . Then, defining i(c † q,2 c k+q,1 e i(k+q)δz − c † q,1 c k+q,2 e −iqδz ) = c † q c k+q , (H5) and using the linear response approximation, i.e., replacing ρ I (τ ) by ρ I (0), one obtains In doing so we neglected the term Tr[e iH0t c † q,2 c k+q,1 e −iH0t ρ 0 ] since it vanishes for q = 0. We note that for a cylinder geometry with an even number 2N x of sites in one unit cell, the Hamiltonian has exactly N x positive eigenvalues and N x negative eigenvalues for each momentum q. To make this explicit in the notation of the eigenstate γ q,j , we define the index j to range from −N x to −1 and from 1 to N x . Thus, the Hamiltonian can be written as with qj the corresponding eigenvalues. Moreover, due to the property of Majorana fermions, γ † q,j = γ −q,−j , we have the relation q,j = − −q,−j . Therefore, {γ qj , γ † q j } = δ qq δ jj and {γ qj , γ q j } = δ −qq δ −jj .
With these definitions, the perturbation term can be written as and one can show that Using Eq. (H7), the time evolution of H I (τ ) can be obtained as Moreover, for our perturbation H = J p (ic m,2 c m+1,1 − ic m+1,1 c m,2 )/8 with m at the center of the edge, V j,j k+q,k can be computed as Here,f l * qj is the eigenmode component of the Majorana fermion in momentum space, i.e., c † ql = jf l * qj γ † qj .
Then, using the identity Tr(γ † kj γ k j ρ 0 ) = δ jj δ kk for j, j < 0, we obtain and we dropped the time-independent terms along the way. From Eq. (H14), one can calculate the Fourier spectrum by using Eq. (F7) with = t/T . We can then compare the resulting expression to our exact numerical simulation, shown in Fig. 3(b). The results are shown in Fig. 15, showing an almost perfect agreement. This allows us to analytically analyze the dominant signal, from which we extract the edge velocity. To this end, we employ the delta-function identity ∞ 0 e iωt dt = πδ(ω) to calculate the Fourier spectrum, in the limit t → ∞. This gives C(k, ω) = δ(2 q,j − 2 k+q,j + ω) −2 q,j + 2 k+q,j −F j j q,k+q δ(2 q,j − 2 k+q,j + ω) 2 k+q,j − 2 q,j V j,j k+q,q .
We recall from Fig. 2(a) in the main text, that when q y = π, the eigenstates are highly degenerate and have energy ±J/4 (crossing point of all curves). There are N x − 1 degenerate energy levels in each of the positive and negative dispersion branches. Hence, C(k, ω) will be dominated by those values of ω that match the energy difference in the argument of the delta functions, so that all degenerate states contribute. To account for this degenerate contribution, one can take one of the two states in the delta function to be at the degeneracy point (i.e., q=π,±j = ±J/4, with j > 1), and the other at the edge state (i.e., k+q=π+k,±1 = ±v edge k). Then, the two delta functions in the equation above become δ(2 q,j − 2 k+q,j + ω) = δ 1 2 J + 2v edge k − ω , δ(2 q,j − 2 k+q,j + ω) = δ 1 2 J + 2v edge k + ω . (H16) We verified that this result matches our numerical simulation; see Fig. 15. In particular, we conclude that the correlation between the highly degenerate state and the edge state dominates the signal |A(k, ω)|. Importantly, the result in Eq. (H16) shows that the dominant spectroscopic signal can be used to probe the dispersion relation ω = (1/2)J + 2v edge k, and hence, to extract the velocity of the edge mode v edge , as we now describe in the following section.
Appendix I: Extracting the edge mode velocity Finally, we turn our attention to the fitting procedure that we use to extract (twice) the chiral edge mode velocity v edge from the Fourier spectrum defined in Eq. (F7).
Consider the Fourier spectrum as a dataset of triples D = { k j y , ω , |A(k j y , ω )| }, where to each Fourier momentum and frequency point, we associate the corresponding value of |A(k j y , ω )|. Our goal is to extract the slope of the most dominant part of the Fourier signal, which corresponds to the discernible chiral pattern. Thus, in constructing the dataset D, in order to focus the signal on the dominant chiral region, we only consider points (k j y , ω ) for which the strength of the Fourier spectrum is larger than 40% of the maximum available signal |A(k j y , ω )|. Comparing this problem to ordinary linear regression, here we have to take into account the strength of the Fourier spectrum. Thus, we apply a weighted linear regression, where each (k j y , ω ) point is additionally multiplied by its strength |A(k y , ω)|. The corresponding cost function reads as L(v edge ) = This cost function L can be interpreted as the variance or error of the fit. Thus, we can use the value L(v edge ) to define the error bars on the extracted value for v edge . Note that, in principle, one can also leave the ω-axis intercept as a fitting parameter; instead, we use the value 1/2 which is obtained from the linear response analysis (cf. Eq. (H16)); we also verified this independently from the fit in the clean theoretical regime of large N T and t ramp . By adopting this value, we reduced the fitting process to only one parameter, namely, the slope. Thus, a limited number of critical data points are sufficient for the fitting, and these can be obtained by selecting appropriate hyperparameters; the procedure maintains a high degree of robustness as long as the cutoff is sufficiently high and the introduced noisy data points are negligible. Figure 16 shows the dependence of the extracted edge velocity on the drive frequency ω D with different cutoffs; we use the scaling N T = 50ω D /J +1 to keep the total physical observation time t obs = T N T approximately the same throughout the drive frequency axis, and Jt ramp =250. By using these parameters and ignoring the k = 0 point, the fitting procedure shows a decent resolution when the chiral gap and the cutoff are sufficiently large (i.e., for J/ω D > 0.06 and cutoff value not less than 0.4). However, when the gap or the cutoff value is small, the effect of bulk states excited during the ramp cannot be disregarded these excitation points then lead to large error bars as exhibited by the J/ω D < 0.06 region or the cutoff value of 0.2 (red data points in the figure). Thus, potential experiments should target the parameter regime with the largest accessible chiral gap and choose a relatively large cutoff in the fitting.