Quantum state preparation of topological chiral spin liquids via Floquet engineering

In condensed matter, Chiral Spin Liquids (CSL) are quantum spin analogs of electronic Fractional Quantum Hall states (in the continuum) or Fractional Chern Insulators (on the lattice). As the latter, CSL are remarquable states of matter, exhibiting topological order and chiral edge modes. Preparing CSL on quantum simulators like cold atom platforms is still an open challenge. Here we propose a simple setup on a finite cluster of spin-1/2 located at the sites of a square lattice. Using a Resonating Valence Bond (RVB) non-chiral spin liquid as initial state on which fast time-modulations of strong nearest-neighbor Heisenberg couplings are applied, following different protocols (out-of-equilibrium quench or semi-adiabatic ramping of the drive), we show the slow emergence of such a CSL phase. An effective Floquet dynamics, obtained from a high-frequency Magnus expansion of the drive Hamiltonian, provides a very accurate and simple framework fully capturing the out-of-equilibrium dynamics. An analysis of the resulting prepared states in term of Projected Entangled Pair states gives further insights on the topological nature of the chiral phase. Finally, we discuss possible applications to quantum computing.

In 1982, Richard Feynman first proposed that one quantum system could be used to simulate the dynamics of other quantum systems of interest [12].Nowadays, quantum simulators based on cold atoms on two-dimensional (2D) optical lattices are being realized experimentally with the realistic perspective of emulating simple models of condensed matter physics [13][14][15][16][17] or studying topological quantum matter [18,19].Rydberg atoms platforms are also used to realize 2D quantum phases of matter [20] like dimer liquids or spin liquids [21,22].Recently, a wavefunction with non-Abelian topological order could be prepared using an adaptive circuit on a trapped-ion quantum processor [23].
Designing atomic platforms of quantum spins with a T-breaking term to emulate CSLs is still an open challenge.Very recently, using Floquet engineering, the emulation of the non-Abelian CSL of the Kitaev honeycomb model in the presence of a chiral term has been proposed theoretically using cold atoms or Rydberg simulators [24,25].Preparing CSL with full spin-rotational symmetry is also very challenging and of great interest.Proposals based on synthetic gauge fields in N-colors SU(N ) Hubbard models on the square lattice with 1 fermion/site were shown to be effective for N ≥ 3 [26].Here we rather focus on the preparation of SU(2)-symmetric spin-1/2 CSLs.
Quench dynamics is an active area of study in recent years [27][28][29][30][31][32], relevant to a wide range of physical systems, including condensed matter systems, ultracold atomic gases, and quantum field theories.While simple quenches provide a way to probe the non-equilibrium dynamics of quantum systems, they cannot serve to prepare quantum phases of matter under specific conditions (low temperature, etc...).In contrast, Floquet engineering consisting of applying a fast and strong periodic drive to the system -described by Floquet theory [33] -is an appealing route towards phase/state preparation [34], and even towards non-equilibrium phases without an equivalent counterpart [35].For example, such a scheme has been proposed to trigger non-trivial energy bands via engineered gauge fields [36].Our aim is here to design some simple protocols based on Floquet theory to prepare genuine spin-1/2 CSLs.For that purpose, we have developed new theoretical tools to compute faithfully the time evolution of our 2D quantum spin systems in order to address Floquet dynamics in various experimental setups, e.g. a sudden quench of the drive or a semi-adiabatic ramping of the drive.Starting from a paradigmatic spin liquid, the Resonating Valence Bond (RVB) state [37], the finite system is slowly driven out of its equilibrium state, and the time evolution of the system becomes very complex, showing rapid oscillations (micromotion).However, using a high-frequency Magnus expansion [38][39][40][41], we analytically construct an effective Floquet Hamiltonian which describes accurately the (stroboscopic) time evolution and makes it much easier to compute numerically.
Under the application of the fast drive, preserving the spin-rotational invariance (singlet character) of the state, we observe a steady increase of the plaquette chirality (more rapid for a quench than for a semi-adiabatic ramping) signaling the emergence of a CSL.Using tensor network techniques in 2D [42], more specifically a chiral Projected Entangled Pair State (PEPS) ansatz [43] (see Appendix A for details), the properties of the final state are uncovered, showing topological features and chiral edge physics characteristic of the Abelian SU(2) 1 CSL [44].
2 Setups and methods

General considerations on time evolution
The considerations below are fairly general and apply to a closed system (i.e. with no external bath) defined by some initial many-body state |Ψ 0 and its unitary evolution |Ψ(t) = U (t, t 0 )|Ψ 0 .In practice, we assume the system is finite.Let us first consider a time-independent Hamiltonian H -typically a Floquet effective Hamiltonian -applied (abruptly) at time t 0 .The unitary time evolution operator U (t 0 + t, t 0 ) is just given by exp (−iHt) (setting ℏ = 1 here) and its action on any state |Φ can be obtained by constructing the Krylov space {H m |Φ } recursively.Then, in principle, one can numerically obtain the exact evolution of the quantum state on a finite system, within machine precision (provided that the computer memory is large enough to accommodate the exponentially large Hilbert space).
In the case of a (more microscopic) time-dependent Hamiltonian H(t), the unitary time evolution operator involves time-ordering, i.e.
In practice, one needs to discretize the time interval by regular time steps τ = t/M .Then, a Trotter-Suzuki (TS) decomposition [45,46] of the unitary time evolution can be performed in term of elementary gates, where t n = t 0 + (n − 1 2 )τ .The time-ordering operator T t ′ can then be seen as the τ → 0 limit of Eq. 1.In practice, τ = t/M is a small time step, significantly smaller than the smallest relevant time scale in the problem, typically the drive period, (τ ≪ T drive ).The TS decomposition involves Trotter errors of order τ 2 which can be controlled by choosing τ small enough.The states |Ψ(τ ) , |Ψ(2τ ) , ... |Ψ(M τ ) are obtained recursively using where using the normalisation condition | Ψ(nτ )|Ψ(nτ ) | = 1 as a test of convergence of the power-expansion up to order m max .The Krylov's basis is computed recursively using

Initial and targeted spin liquid states
We now move more specifically to our system constituted by an assembly of spin-1/2 on a square lattice.The initial prepared state considered in this work is a simple spin-1/2 nearest-neighbor (NN) Resonating Valence Bond (RVB) state [37].It can be viewed as an equal-weight superposition of NN singlet (hard-core) covering of the lattice.Conveniently, this spin liquid has a simple Projected Entangled Pair State (PEPS) representation [47,48] (see Appendix A) which enables to construct it easily on finite systems or on the infinite plane.The NN RVB state has very short-range spin-spin correlations (with a spin correlation length of the order of the lattice spacing) but critical dimer-dimer correlations (on any bipartite lattice).However, any admixture of longer singlets generates a finite dimer correlation length [49,50].In that case, the gapped RVB spin liquid acquires Z 2 topological order, similarly to the Kitaev toric code [51].
Most interestingly, possible (closely related) experimental realizations have been proposed recently using Rydberg platforms [22].Hereafter, we shall consider a 4 × 4 finite torus which conveniently offers lattice translation symmetry.However, we believe our setup proposal and our findings are also relevant to lattice systems with open boundaries, a geometry encountered in experiments.
After applying a (fast) drive over a sufficiently long time interval we wish to prepare the system in an (Abelian) Chiral Spin Liquid (CSL) phase, possibly at low effective temperature (if not in the ground state).It is known that such a phase is realized in a simple chiral spin-1/2 Heisenberg Hamiltonian on the square lattice, which we view as a "target" Hamiltonian, where the last term involves cyclic (let's say anti-clockwise) permutations P ijkl on all plaquettes □.P ijkl performs a shift along the ijkl ring of any spin state This term breaks time-reversal (T) and parity (P).Note that the cyclic permutation on the (oriented) plaquette (i, j, k, l) ≡ (i, i + e x , i + e x + e y , i + e y ) has a simple expression in terms of the spin degrees of freedom, i as can be verified by a direct comparison of the 2 4 × 2 4 matrices in the S z basis on both sides of the above identity.
The other term H 0 is the spin-1/2 quantum antiferromagnet with NN antiferromagnetic coupling J 1 and, optionally, next-NN (NNN) antiferromagnetic coupling J 2 < J 1 , In fact, the phase diagram of this model hosts the desired (gapped) topological CSL phase over a significant region of its parameter space [52][53][54][55], even for a vanishing frustrating coupling J 2 = 0, but for finite λ chiral .Interestingly, the corresponding CSL phase can be represented accurately by PEPS [43,44,53,56], which conveniently give crucial insights on the topological features and on the nature of the chiral edge modes, via the bulk-edge correspondence PEPS theorem [57].The (ideal) monochromatic drive we consider involves dephased sinusoidal modulations of the spin-spin interaction on the NN a, b, c and d bonds around all the sites of one of the two sublattices (let's say A).More precisely, we assume,

Fast chiral periodic drive
where J is the coupling constant and e x , e y are the unit vectors along the crystal axis.
As shown in Fig. 1 the bond pattern is staggered leading to two types of plaquettes adcb and cbad.Interestingly, H drive (t) is chiral, breaking parity P and time-reversal T but not the product PT.Indeed, P (reflection w.r.t. a crystal axis) changes the plaquette bond ordering, e.g.adcb → abcd, i.e. the sign of the phases in Eq. 5. Time-reversal T : t → −t has the same effect.Note that the sign of the chirality in the two types of plaquettes is the same so that the overall chirality is uniform.Note also that the initial state is a spin-singlet and the Hamiltonian in Eq. 5 is invariant under SU(2) spin rotations.The time evolution therefore occurs in the spin-singlet manifold.
It can be noticed that the sequence we are considering has similarities to the driving protocol introduced in Ref. [58] for a tight-binding model of non-interacting particles on the square lattice, in which hopping amplitudes are varied in a spatially homogeneous but chiral, time-periodic way resulting in chiral edge states.
Here the drive is assumed to be very fast compared to the natural time scale ℏ/J, i.e. ℏω ≫ J. Hereafter we set ℏ = 1 (unless useful to keep it) and take ω = 2πf with a frequency f = 10.Hence, the drive time-periodicity T drive is 0.1 (in units of 1/J taken as the unit of time throughout), which is the smallest time scale in the problem.Therefore, to compute the time evolution using a TS decomposition we have to divide the drive period into a sufficiently large number N τ of time steps multiple of τ = T drive /N τ .
Such a dephased monochromatic drive may be implemented on experimental cold atom platforms [19].We shall discuss variations of this setup based on square modulations, also adapted to experiment, later in the conclusion section.Note that the spin-spin interaction may be challenging to realize on analog platforms.Digital (noisy) quantum computers may be an alternative.

Effective Floquet Hamiltonian
In the case of fast drives, different types of high-frequency expansions can be performed in powers of 1/ω.For practical reasons (which will become clear soon), we follow here the Magnus expansion for the stroboscopic Floquet Hamiltonian [41].In that case, the effective (time-independent) Floquet Hamiltonian describes the slow motion at stroboscopic times, T p = pT drive , p ∈ N. The fast micromotion between stroboscopic times is set by the periodic kick operator (not computed) which vanishes at t = T p for all integer p.The drive Hamiltonian has no constant part, (1/T drive ) T drive 0 H drive (t)dt = 0, and just singleharmonics components, where and The leading-order effective Floquet Hamiltonian is therefore [41], where J F = J 2 4ω and P ijkl is the 4-site cyclic permutation which applies to all plaquettes □.The details of the derivation is left in the Appendix B. It is interesting that the next order of the expansion is O(1/ω 3 ) so that we expect H 0 F to capture accurately the stroboscopic motion1 .H 0 F is a chiral Hamiltonian, acquiring its P and T symmetries from the microscopic model H drive : it transforms into its opposite under P or T and is invariant under P T .Interestingly, it is fully symmetric -invariant under lattice translations and under all discrete symmetries of the C 4v point group -in contrast to H drive which bears two types of plaquettes (of same chirality) as seen Fig. 1.The Floquet Hamiltonian is introducing the longest time scale in the problem t F = ℏ/J F = (4ℏω/J) ℏ J .Hence we are in the regime where, with 3 well-separated time (or energy) scales.Next, we shall also consider the case where a constant term κH 0 is added to the drive Hamiltonian leading to a more generic Floquet Hamiltonian H F .For convenience we shall denote |Ψ Floquet (t) the state evolved according to the Floquet Hamiltonian exp (−iH F t)|Ψ RVB , to be distinguished from the state |Ψ RVB following the exact evolution.

Quench and adiabatic evolution
The objective of the setup is to emulate the low-energy physics (possibly the groundstate) of the (static) chiral antiferromagnetic Hamiltonian H target on the square lattice defined in Eq. 3. It is known that the phase diagram of this model bears a wide CSL phase, even for a vanishing frustrating coupling J 2 = 0.While emulating the NN coupling J 1 between the spins can be done experimentally, generating the necessary chiral term is a difficult task.For this aim, we add the Hamiltonian H 0 to the drive Hamiltonian with some (small) amplitude κ to be adjusted, We also consider the possibility of slowly ramping the amplitude J of the drive by some time-dependent factor 0 ≤ λ(t) ≤ 1.Two scenarios will be considered, either (i) a quench setup where λ(t) = 1 from the very beginning of the drive application at t = 0 or (ii) a quasi-adiabatic setup where λ is increased smoothly from 0 at t = 0 to 1 at t = t ramp .
The most interesting regime to consider is the so-called strong (and fast) drive regime κJ 1 ≪ J, i.e. κ ≪ 1 if J and J 1 are both taken as energy reference, J = J 1 = 1.If λ(t) is not too small (i.e.excluding the short initial branching of the drive), the new high-frequency effective Floquet Hamiltonian is simply, which is exactly of the form of the target Hamiltonian (3).Note that in the case of a slow ramp, the Floquet Hamiltonian itself acquires a slow time dependence.When a quench is realized, λ(t) = 1, one can choose κ so that H F = κH target .Using ( 3) and ( 10), one gets Note that, because of the presence of H 0 in Eq. 9, there are corrections to Eq. 10 of order O(1/ω) instead of O(1/ω 3 ) (see Eq. ( 42) of Ref. [41] and Appendix B) .However, these additional terms have amplitudes κJ 1 J/ω i.e, from Eq. 11, of the order of J 1 J 2 /ω 2 and, hence, can be negleted compared to the terms in Eq. 10.

Validity of the Magnus high-frequency expansion
In order to check the relevance of the Floquet Hamiltonian to describe the effective dynamics, we first consider the simplest possible setup, namely the application of the fast chiral drive alone on the NN RVB state, i.e. we set x = 0 and λ(t) = 1 in Eq. 9. Our initial state is invariant under P and T while H F is odd so we expect H F ≡ Ψ Floquet (t)|H F |Ψ Floquet (t) = 0, i.e. the stroboscopic motion occurs in the zero-chirality manifold.Considering a 4 × 4 finite cluster with PBC (torus), we have checked numerically that this property is even exactly realized during the micromotion under H drive (t).Indeed, if one considers any 4-site oriented plaquette (i, j, k, l), at all times, i.e facing triangles in all plaquettes have always opposite chiralities (of small magnitudes).Hence, the overall plaquette chirality Ψ(t)|i(P ijkl − P −1 ijkl )|Ψ(t) is exactly vanishing at all times.To have a more precise comparison between the exact and the effective dynamics, we have computed the infidelity w.r.t. to the initial state, i.e.
, starting from the same initial state |Ψ RVB at t 0 = 0.The results in Fig. 2(a) show an excellent agreement.A zoom at small time in Fig. 2(b) shows a clear oscillation within every time period (micromotion), while the motion at (discrete) stroboscopic times (p/N τ + n)T drive , p fixed, falls on smooth curves when varying n.The result at p = 0, that is at times exact multiples of T drive , compares very well with the calculation obtained using H F , i.e. approximating |Ψ(t) by |Ψ Floquet (t) = exp (−iH F t)|Ψ RVB .
However, a zoom at a longer time t ∼ 30 in Fig. 2(c) shows small deviations between the stroboscopic motion computed with the drive Hamiltonian using two choices of N τ and the effective Floquet motion (computed exactly).These deviations may have two origins: first, the Trotter errors in the time discretization of the drive oscillations and, secondly, the neglected higher order terms in the Floquet Hamiltonian.The truncation of the Floquet Hamiltonian leads to relative errors of order 1/ω 2 , i.e. ∼ 2 × 10 −4 which can be neglected.In contrast, the Trotter errors are still sizeable at N τ = 8.Using N τ = 16 H target given by J 1 = J 2 = 0 and λ chiral = 0.5; (b) κ ̸ = 0 given by Eq. 11, H target given by J 1 = 1, J 2 = 0 and λ chiral = 0.7; (c) Same as (b) but for J 2 = 0.4 and λ chiral = 0.5.Energy levels are represented by black dotted lines and their weight by segments whose 1/2-length is equal to ρ.Note that µ ρ µ = 1.The first singlet excited state, forming together with the GS the pair of relevant CFT states [52], is shown by a green line in (b) and (c).In both cases, note the presence of a low energy singlet state with (π, π) momentum (shown in red).In case (c), green and red states are very close but non-degenerate.
Trotter steps per drive period seems to give already a much better accuracy, typically within 2% of the (exact) effective Floquet motion.Interestingly, we also observe that the amplitude of the micromotion oscillations becomes quite small compared to the average behavior, except at very small time.

Quench of the fast drive
In this section we shall assume the drive Hamiltonian ( 10) is suddenly switched on at t = 0 on the RVB initial state, i.e. we assume λ(t) = 1, independent of t.In that case, the effective motion follows the unitary evolution given by U (t) = exp (−iκH target t), which depends on the control parameters J 1 , J 2 and λ chiral .The case J 1 = J 2 = 0 corresponds to the previously studied case H 0 = 0 i.e. a purely oscillatory drive.If the constant term H 0 is present we set J 1 = 1.

Long-time average and diagonal ensembles
Since the effective Floquet dynamics is expected to be very accurate, one can use it to investigate the physics in the long time limit, after a quench of the drive Hamiltonian, i.e. assuming λ = 1 in Eq. 10.Strictly speaking, the approximate effective Hamiltonian computed from the Magnus expansion does not capture heating effects.The system is expected to heat up to an infinite temperature state in the t → ∞ limit for any finite κ > 0, see e.g., Ref. [35].Note however that on a finite system of N sites, the bandwidth of the many-body energy spectrum scales like W ω ∼ aN J F = (a/4)N J 2 /ω, where a of order 1.Hence, at large enough frequency, the many-body spectrum does not reach the boundary of the Brillouin-Floquet energy zone of extension ω.In our case a ∼ 2.4 (see numerical computations below) and W ω ≃ 0.15J, significantly smaller than ω = 20πJ.We therefore expect an extremely long prethermalization regime [35].
Decomposing the initial (non-chiral) spin liquid state in terms of the Floquet eigenstates |e µ of H F , |Ψ RVB = µ w µ |e µ , the infinite-time average of observables O can be written in terms of the diagonal density matrix ρ diag = µ ρ µ |e µ e µ | [41], where ρ µ = |w µ | 2 .Within a time period, the wavefunction |Ψ(nT + t) is subject to a unitary rotation exp (−iK(t)) where K(t) is a periodic kick operator (see explicit expression in Appendix C).Hence, the infinite-time average in Eq. 12 involves the averaged operator, Note that, for the observables used in this work, the effect of the kick rotations can be neglected so that O µµ ≃ O µµ in Eq. 12.Note also that, since the RVB state is a singlet state invariant under all space group operations, its decomposition in terms of Floquet eigenstates is limited to the most symmetric singlet sector.When κ = 0 (sinusoidal drive alone) the initial RVB state is a highly excited state of H 0 F .Indeed, H 0 F has a symmetric spectrum with ±E µ opposite eigenvalues associated to pairs of complex-conjugate eigenvectors |e 0 µ , − = |e 0 µ , + * .The RVB state being a real (T-invariant) wavefunction, its decomposition is of the form where c 0 µ are real coefficients.The E µ ↔ −E µ symmetry of the weight distribution is indeed clear in Fig. 3(a).One also immediately gets H 0 F ave = 0, in agreement with the fact that the chirality remains zero under time evolution, as noted before.Interestingly, we see that the initial RVB state has a very large overlap with the zero energy states located at the middle of the many-body spectrum of H 0 F .For larger L × L systems of increasing N = L 2 sites, a broader weight distribution is expected, still centered around the spectrum center E = 0, so that the system can be viewed as being at infinite temperature.A more exotic scenario would be realized if the RVB state keeps a finite overlap on the E = 0 state manifold suggesting the existence of low-entangled "scar" states [59,60] (i.e.following the area law) in the center of the spectrum and, hence, a lack of thermalization.
Weight distributions for κ ̸ = 0, i.e. when a constant term is present in the drive Hamiltonian, are shown in Fig. 3(b,c), in the cases corresponding to H target given by J 1 = 1, J 2 = 0 and λ chiral = 0.7, and by J 1 = 1, J 2 = 0.4 and λ chiral = 0.5, respectively.Then, the spectrum of H target = 1 κ H F is no longer symmetric and one gets finite values for the time-averaged plaquette chirality Im P ijkl ave ∼ 0.2496 and 0.2334 in the (b) and (c) cases, respectively.The GS and the 1st singlet excited state (shown as a red line in Figs.3(b,c)) have been identified from a Conformal Field Theory (CFT) construction as the two relevant states forming the CSL doubly-degenerate GS manifold on the infinitesize torus [52].Note that the RVB state has no overlap with the second CFT state which possesses different quantum numbers.However, large weights are found on the GS and on the next lowest fully symmetric singlet eigenstates of H F .

Time evolution at small and intermediate times
We now turn to the actual computation of the time evolution in the case of a quench of the drive.Results obtained on the 4 × 4 torus are shown in Fig. 4(a).The constant H 0 of the quenched Hamiltonian (10) involves either J 1 = 1 and J 2 = 0 or J 1 = 1 and J 2 = 0.4.The amplitude κ of H 0 is set using Eq.11 i.e. κ = J 2 /(4ωλ chiral ), and we choose either λ chiral = 0.7 when J 2 = 0 or λ chiral = 0.5 when J 2 = 0.4, corresponding to the two realistic H target whose spectra are shown in Fig. 3(b,c).We observe a rapid increase of the plaquette chirality, initially vanishing in the RVB state.Note that the amplitude of the micromotion fast oscillations remains quite small, barely visible on the scale of the plot, so that the full time evolution can be compared directly to the effective Floquet dynamics, not only at stroboscopic times.We observe that the effective Floquet approach is able to capture accurately the slow variations of the chirality expectation value up to long times.Even more, we believe that the small deviation at intermediate times with the full time-dependent simulations is primarily due to the Trotter error in the later rather than to the truncation of the Floquet Hamiltonian at order 1/ω.At times t > t F we observe in Fig. 4(b) some strong oscillations around the expected long-time average values computed from Eq. 12.This is clearly a finite size effect and one expects the amplitude of the oscillations to decrease rapidly with increasing system size.In any case, we believe that quenching the drive Hamiltonian is not the correct procedure to achieve the goal of preparing the system in the final CSL state and we now move to a different setup.

Adiabatic ramping of the fast drive
The new setup we shall follow now consists in branching the fast drive in a quasi-adiabatic way.For that we have chosen a simple linear ramp, for t ∈ [0, t ramp ].In a truly infinitely-slow ramping λ(t) from 0 to 1 (e.g.taking t ramp → ∞), the state |Ψ(t) is expected to follow adiabatically the GS of provided that |Ψ(0) is the GS of H 0 .In practice, deviations from this ideal scenario are present.First, the parent Hamiltonian of the RVB state in the thermodynamic limit involves longer range interactions [47,48] than the NN and NNN interactions considered here.However, we have checked that, on a finite system, the RVB state is a good approximation of the exact GS of H 0 for J 2 = 0.4J 1 , with an overlap of order 0.9960 on a 4 × 4 cluster.Secondly, t ramp is finite leading to a deviation from the adiabatic behavior.
It is then interesting to try to quantify the effect of a finite ramp time.Note that, although the branching of the drive is linear in t, the effective coupling to H target in ( 16) is smoother, increasing only as t 2 at small time.Note also that the initial RVB state has a sizeable overlap with the targeted final state, as shown in Fig. 3(c), which should enable an adiabatic evolution in a shorter time.
We have computed the time evolution (using N τ = 16 TS steps per drive period T drive = 0.1) and the emergence of the plaquette chirality for three ramp times of increasing length, t ramp = 80, 160, and 240 as plotted in Fig. 4(a), showing a smooth behavior of the observable, in contrast to the quench setups.As expected, the increase of the chirality is slower and slower when increasing the ramp time but the final chirality reached at t = t ramp becomes larger and larger.varying λ "by hand" from 0 to 1.It should be reached asymptotically by taking the limit t ramp → ∞ keeping the ratio t/t ramp = λ fixed.In fact, the GS chirality of H target is ∼ 0.478 while the achieved chirality at t = t ramp = 320 is only ∼ 0.383.This indicates that longer ramp times are needed to achieve the goal of approaching more closely the GS of H target .However, running the full drive simulation over longer times is expensive in CPUtime and, in fact, not necessary.Indeed, we have checked that for t ramp = 320, the effective evolution with the slowly varying Floquet Hamiltonian H F (t) gives very similar results as seen in Fig. 5(a).In particular, the overlap | Ψ Floquet (t)|Ψ drive (t) | between the two final states at t = t ramp is 0.99908, very close to 1.The effective Floquet dynamics allows to use much longer Trotter steps, typically of order 0.5, compared to T drive /N τ = 6.25 × 10 −3 for the full dynamics, and is therefore a much cheaper alternative for longer times.Hereafter, further results will then be obtained using Floquet dynamics, showing no appreciable difference on the plots and performed easily up to t ramp = 960 (i.e.around 10 000 drive time periods!).For t ramp > 500, while approaching the adiabatic limit, one starts to see oscillations in the time evolution of the chirality which are attributed to the deviation of the initial RVB state from the true GS of H 0 .
To confirm the origin of the spurious oscillations found above, we consider the simple quantum Heisenberg antiferromagnet, setting J 1 = 1, J 2 = 0 in H 0 , and use its GS as initial state, instead of the RVB state.Note that, for this new choice of H 0 , the overlap between the RVB state and the true GS is only 0.9256 on the 4 × 4 torus and, hence, the RVB is no longer an appropriate choice for the initial state.Results are shown in Fig. 5(b), showing a smooth evolution towards the adiabatic limit when increasing t ramp , even at large t ramp , in contrast to Fig. 5(a).

Fidelities of the final CSL state
The results shown in the previous sections strongly suggest that the ramp setup can give efficient implementations of the CSL state.Here we try to analyse more quantitatively the properties of final state.

Comparison with the adiabatic limit
First, we show in Fig. 6, for various values of t ramp , the overlap of the time-evolved state |Ψ tramp (t) with the adiabatic GS |Ψ adiab (λ) , computed at λ = t/t ramp , as a function of the reduced time t/t ramp .As expected, we observe that, for a given ramp time, the fidelity deteriorates as a function of time.However, for a fixed ratio λ = t/t ramp defining a target (chiral) state |Ψ adiab (λ) GS of a target Hamiltonian where λ chiral is replaced by λ 2 λ chiral , the fidelity quickly increases with t ramp .For λ = 1, corresponding to λ chiral = 0.7 in H target , a fidelity of ∼ 0.854 can be obtained using t ramp ≃ 1000.For λ = 0.8, corresponding to λ chiral ≃ 0.45, a fidelity of ∼ 0.953 can be reached with the same value of t ramp , but after a shorter time λ t ramp ≃ 800.

Comparison with chiral PEPS
We expect the final state to be in a CSL phase, in the close neighborhood of the GS of the target Hamiltonian.The later has been shown to be represented faithfully by chiral PEPS ansatze [44,56].It is therefore instructive to compare our final states prepared at t = t ramp with optimized PEPS.We shall consider fully symmetric, i.e.SU(2) and translationally invariant, PEPS with bond virtual spaces ν = 0 ⊕ 1/2 (D = 3), ν = 0 ⊕ 1/2 ⊕ 1/2 (D = 5) and ν = 0 ⊕ 1/2 ⊕ 1 (D = 6).Within these symmetric PEPS manifolds, restricted A 1 + iA 2 chiral sub-manifolds (with same bond dimensions) are designed to represent our Figure 6: Overlap |⟨Ψ tramp (t)|Ψ adiab (λ = t/t ramp )⟩| between the time-evolved state and the adiabatic GS computed at λ = t/t ramp , plotted versus t/t ramp .Here the time-evolved state is obtained using the simplified Floquet dynamics with J 1 = 1, J 2 = 0 and λ chiral = 0.7 and various finite ramp times t ramp up to t ramp = 960 are considered.The GS of H 0 is taken as initial state.The inset is a zoom of the same quantity for t/t ramp ≤ 0.5.CSL (see Appendix A for more details).While generic manifolds possess (n tensors − 1) complex independent parameters, where n tensors is the total number of A 1 and A 2 tensors (see Table 1 for numbers), the chiral PEPS manifolds (with a fixed relative phase between the A 1 and A 2 tensors) have only half degrees of freedom i.e. (n tensors − 1) real parameters to be optimized.
As shown in Fig. 7, as expected the D = 3 PEPS cannot accommodate very well, for increasing t ramp , the increase of entanglement of |Ψ tramp signaled by the rapid decrease of the fidelity w.r.t. the initial state.However, one already sees that the performance of the chiral PEPS (with only two variational parameters) becomes comparable to the generic (symmetric) PEPS for t ramp > 250, suggesting a crossover towards the CSL phase.The extension of the chiral PEPS to D = 5 or D = 6 unfortunately does not improve the overlaps significantly (see Table 1) but a similar trend is found as shown in Fig. 7(b) for D = 6.

Discussions and conclusions
We first discuss possible experimental realisations of our setup, how to realize the fast drive and how to prepare the initial state.

Realizing the drive Hamiltonian
The Heisenberg interactions in Eq. 5 can be realized with cold atoms loaded on an optical lattice [61], for which a microscopic Hubbard-like description applies [62,63].The excellent isolation of the optical lattice from its environment enables to realize ideal unitary dynamics.We assume one atom per site (on average) and a large on-site repulsion U to be in the Mott insulating phase [64], and different hopping terms t v (t) on the four NN bonds v = a, b, c, d of Fig. 1.Our setup is relatively simple as it involves only the modulation in time and space of the NN (effective) Heisenberg couplings.A (selective) modulation in time of the Heisenberg couplings J v (t) of H v can be realized via a (selective) modulation of the hoppings between sites t v (t), i.e. changing the tunneling barriers of the optical lattice appropriately.Note that this procedure is conceptually different from the synthetic gauge field approach [18], in which the chiral term would appear in the Mott insulating phase only in fourth-order of the (complex) hopping t, i.e. λ chiral ∼ t 4 /U 3 .However, in fermionic Hubbard systems the Heisenberg couplings can only be antiferromagnetic (J v = 4(t v (t))2 /U > 0), not ferromagnetic (J v < 0) and, strictly speaking, the singleharmonic modulations of (5) cannot be realized 2 .A more practical alternative would then be to use, on the optical lattice, two-component bosons [65,66] which (i) are more easily brought to low temperatures and (ii) are more versatile than fermions to realize alternating antiferromagnetic and ferromagnetic couplings, tuning the U ↑↑ = U ↓↓ and U ↑↓ on-site Hubbard interactions.
It is interesting to mention that a simplified drive sequence consists in replacing the sinusoidal modulation by a square wave i.e. ], ], This new type of drive sequence introduces higher harmonics 2 nπ (exp (inωt) ± c.c.) in the amplitude modulations of H x and iH y which modify the prefactor J F in Eq. 7 into and the chiral term in Eqs.7 and 10 accordingly.Note also the emergence of chiral terms in the (targeted) Floquet effective Hamiltonian at order O(ω −2 ) instead of O(ω −3 ).Nevertheless, we believe the main ideas/results developed in this work are still valid qualitatively in the case of (17).Finally, we would like to mention that Rydberg atom arrays can also emulate quantum spin systems [67], but keeping SU(2) spin rotational invariance might be challenging.Nevertheless, it would be interesting to investigate whether our setup can be exported to the case of anisotropic XXZ dipolar spin interactions.

Preparation of the initial state
Spin liquids of the RVB type can be realized on cold atom platforms i.e. in Mott insulating ultracold bosons [68].Using a two-component one-dimensional Bose-Hubbard system, a highly entangled Heisenberg antiferromagnet with SU(2) symmetry has been experimentally realized [69].We believe that this promising approach could be extended to two dimensions.Note however that, as shown by our studies, it is important that the initial state in the ramp process is as close as possible to the GS of the Hamiltonian at t = 0 i.e.H(0) = H 0 .Using cold atoms on an optical square lattice, the Mott phase can be realized on fermionic systems of hundred of sites [70,71] with magnetic correlation lengths of several lattice spacings [72].Similarly, systems of two-component bosons could also be brought in Mott phases at large Hubbard repulsions U σσ ′ [66].Owing to rapid experimental progress in decreasing the temperature, we believe the quantum state of such systems may soon approach the (finite size) GS of a plain spin-1/2 NN quantum Heisenberg model (i.e. with J 2 = 0) and could serve as an initial state for the second type of ramping.Tweezer arrays of fermions have also been proposed [73], that could be used to reach the same goal.We point out that the precise nature of the initial state (e.g.spin liquid vs Néel state) is not completely settled on finite size systems which are always quantum-mechanically disordered (i.e. are in a global singlet state) and possess a small gap in their many-body spectrum.In other words, the distinction between different phases is not sharp on finite size systems.

Final remarks
In summary, we have studied different protocols to prepare a (finite size) system of quantum spin-1/2 on a square lattice into a spin-1/2 chiral spin liquid phase.We start from a simple low-entangled (SU(2)-symmetric) initial state like a (non-chiral) spin liquid (e.g. a RVB spin liquid) or the (finite size) GS of a simple NN Heisenberg quantum antiferromagnet which, we believe, could (or will soon) be realized experimentally in platforms of cold atoms on an optical lattice.A simple drive involving dephased periodic modulations of four types of NN bonds is proposed.We focus on the high-frequency/large amplitude (strong drive) regime for which an effective Floquet dynamics is shown to give very accurate results.We argue that our goal cannot be realized by quenching the drive suddenly, but rather by a semi-adiabatic ramping.Quantitative results are given for a periodic 4 × 4 system.In that case, we show that the (finite size) CSL ground state of a relevant target Hamiltonian can be reached with high fidelity typically within less than 10 thousands drive oscillations, still beyond current experimental possibilities.However, high-fidelity counterdiabatic protocols -i.e.far away from the adiabatic limit -respecting experimental limitations (like the limited coherence times) have been proposed [74] and may be adapted to our setup.We would like to emphasize again that it is essential that the system is kept finite and we believe the thermodynamic limit cannot be taken right away for two reasons: first, heating will occur when the many-body Floquet spectrum (expanding with system size as N/ω) will reach the boundary of the Floquet-Brillouin quasi-energy zone of extension ω.This simple argument gives a minimum critical frequency which diverges as √ N for increasing system size; secondly, another issue is the vanishing of the finite size gap in the thermodynamic limit that may lead to a diverging ramp time.
Recent experiments are now performed on systems with hundred sites or even more, and with a box-like confining potential (corresponding theoretically to open boundary conditions).The use of periodic boundary conditions in our system was in fact dictated by convenience, translation invariance making calculations simpler.However, we believe no conceptual difference is to be expected for open boundaries.The case of a harmonic confining potential is more tedious and left for future investigations.In any case, the system size N is an important relevant parameter.Indeed, in the transition region between the non-chiral initial phase and the CSL, the gap should reach a minimum value, typically of order 1/N , or a power of 1/N .Hence, the minimum ramp time to reach a good fidelity should qualitatively scale like N , which could be a limitation for N of a few hundred sites.However we believe that preparation times could be drastically reduced using quantum control algorithms, an interesting issue left for future studies.
Our computation has been performed on the square lattice but similar setups can also be realized in frustrated 2D lattices like kagome or triangular lattices which can be obtained experimentally [75].In that case, dephased (fast) modulations of 3 or 6 groups of equivalent NN bonds can be performed leading to chiral permutations on triangles or hexagons.The resulting topological CSL state and its anyonic excitations may be used for topological quantum computations following Kitaev's proposal [51].
We also note that we have restricted ourselves to the high-frequency regime.However, in Ref. [58] it is argued that there exists novel non-equilibrium chiral phase in the limit of a slow drive.This was then generalized to a Kitaev type model [76].Investigations of our model at lower frequency is of great interest and is currently undertaken looking for a similar phase in an SU(2) symmetric model.
Lastly, we note that the implementation of a dimer spin liquid state (of the toric code universality class) has been proposed on the kagome lattice by using Rydberg atom arrays [67].It would then be interesting to investigate whether a chiral dimer liquid (CDL) could be realized, starting from such an initial state, by Floquet engineering involving dimer degrees of freedom (represented by excited Rydberg atoms in the blockade regime) instead of the original spin-1/2 degrees of freedom considered in this work.

A PEPS representation of spin liquids
PEPS are defined as a tensor network of local onsite tensors of size dD4 where d = 2 is the local physical Hilbert space dimension (spin 1/2) and D is the bond dimension i.e. the number of virtual states on each tensor leg.By contracting the tensors w.r.t. the bond indices, one generates an entangled ansatz (whose entanglement entropy per bond is limited by ln D and, hence, controlled by D).PEPS offer an extremely efficient variational scheme to address local Hamiltonians [77].In addition, PEPS can encode the topological nature of the state and the physics of the edge modes via the PEPS bulk-edge correspondence theorem [78,79].In the infinite-PEPS (iPEPS) algorithm [42] infinite-size systems can be handled.Note that PEPS fail to describe the rapid increase of entanglement entropy (e.g. in the case of a quench, see Ref. [80]) but should still be relevant in the case of an adiabatic (or slow) ramp which leads to a limited growth of entanglement entropy.
Both the initial state and the Floquet Hamiltonian governing the effective dynamics of our problem are SU(2) symmetric and transform according to the trivial representation of the square lattice space group3 .It is therefore relevant to enforce such symmetries at every step of our scheme.A symmetric spin liquid on the square lattice is simply obtained by choosing a C 4v -symmetric site tensor placed uniformly on every lattice sites.As explained in details in Ref. [81], the continuous SU(2) of PEPS can be implemented at the level of local tensors by imposing that the D-dimensional subspace ν of each virtual leg is a direct sum of SU(2) irreducible representations.Working with SU(2)-symmetric PEPS allows to greatly reduce the number of parameters describing the PEPS family and is a major advantage for computations based on optimization.
The NN RVB state [37] which consists of an equal weight superposition of all NN singlet (valence bond) coverings 4 has a simple representation using ν = 0 ⊕ 1/2 virtual bonds (D = 3) [47,48].In fact, the D = 3 symmetric PEPS manifold is represented by two linearly independent fully symmetric (A 1 ) real site tensors (describing the fusion of either one or three virtual spin-1/2 into a physical spin-1/2) and generically describes gapped Z 2 topological RVB spin liquids which include (thanks to "teleportation") longer range singlet bonds [49].The critical NN RVB state is therefore a special point in this two-parameter manifold [50], with exponentially large correlation lengths in its neighborhood [49,50].
The simplest PEPS representation of the (Abelian) spin-1/2 CSL -of the Kalmeyer-Laughlin type [2] -is in fact possible within the ν = 0⊕1/2 virtual space (D = 3) manifold by adding a third site tensor of different A 2 orbital symmetry with a pure-imaginary amplitude, resulting into a A 1 + iA 2 site tensor [43,44].Such an ansatz breaks timereversal T (i → −i) and parity P (A 2 → −A 2 ) while being invariant under their product PT, as expected for a CSL.Also, despite its simplicity, it gives a faithful description of the chiral edge modes -a key feature of the CSL -following closely the prediction of a SU(2) 1 Wess-Zumino-Witten (WZW) conformal field theory (CFT).Moreover, the CSL ground-state [82] of the target Hamiltonian (3) has been successfully studied using similar PEPS [53,56].Increasing the bond dimension from D = 3 to D = 5 (ν = 0 ⊕ 1/2 ⊕ 1/2) or D = 6 (ν = 0 ⊕ 1/2 ⊕ 1) provides a small improvement of the variational energy.We show in Table 1  We start from the expression of the Floquet Hamiltonian at order 1/ω -see Eq. ( 42) in [41] -decomposing H drive as H 1 exp (iωt) + H −1 exp (−iωt), where It is interesting to note that, in this simple case of an harmonic drive alone, corrections only appear at order 1/ω 3 .One can then expand the commutator as [h x , h y ] = i∈A (C i + C ′ i ) where, Using [S α i , S β i ] = iϵ αβγ S γ i we get, After the transformation i + e x → j, C ′ i → C j , we obtain for j ∈ B, C j = (S α j−ex − S α j+ex )[S α j , S β j ](S β j−ey − S β j+ey ) = iS j • (S j+ex × S j+ey ) + iS j • (S j+ey × S j−ex ) + iS j • (S j−ex × S j−ey ) + iS j • (S j−ey × S j+ex ) . ( Because of the small prefactor in the expression of K(t), the unitary remains close to the identity.Most observables are therefore weakly affected during the micromotion w.r.t. to their behaviors under the effective Floquet dynamics.This has been confirmed numerically.

Figure 1 :
Figure 1: Left panel -Sketch of the H drive Hamiltonian given by Eq. 5. Dephased terms H a , H b , H c and H d are represented as blue, yellow, red and green bonds.Note that the bond modulations around the sites A (black dots) and sites B (white dots) of the (bipartite) square lattice are dephased by π.Right panel -Floquet effective Hamiltonian H F given by Eq. 7 acting on square plaquettes as circular anticlockwise (resp.clockwise) permutations with a factor i (resp.−i).

Figure 2 :
Figure 2: Infidelity 1−| Ψ RVB |U (t, 0)|Ψ RVB | w.r.t. the initial state, the NN RVB spin liquid, computed on a 4 × 4 torus in the case of a pure chiral drive H(t) = H drive (t).(a) Comparison of the unitary evolutions with the drive Hamiltonian and with the (time-independent) effective Floquet Hamitonian obtained on a 4×4 cluster.Zooms at small times (b) or around t = 30 (c) are shown.Trotter errors are negligible at short times (b) but become sizeable for N τ = 8 at intermediate times (c).

Fig. 5 (
a) shows similar data as a function of the reduced time t/t ramp .The adiabatic groundstate |Ψ adiab (λ) is obtained by following the GS of the Floquet Hamiltonian (16),

Figure 4 :0
Figure 4: Plaquette chirality vs time for a quench of the drive computed on a 4 × 4 torus.Parameters are J 1 = 1, J 2 = 0 and λ chiral = 0.7 or J 1 = 1, J 2 = 0.4, and λ chiral = 0.5.(a) Short and intermediate times behavior; all computations involving the fast drive are done with N τ = 16 intervals in each drive period but only data at times kT drive and (k+1/2)T drive , k ∈ N, are reported since the amplitude of the oscillations within each time period (micromotion) remains very small.Computations using the effective Floquet dynamics with U (t) = exp (−iκH target t) shown as dashed-dotted curves agree very well.Results involving a ramping of the drive are also shown for comparison.(b,c) Long time behavior obtained using the effective Floquet dynamics; the chirality (blue curve) is shown to oscillate around the average value (black dashed line) computed independently using Eq.(12).The time average 1 t

Figure 5 :
Figure 5: Plaquette chirality vs reduced time t/t ramp for adiabatic ramping of the drive using different values of t ramp , computed on a 4 × 4 torus.(a) The RVB state is used as initial state, approximating the GS of the frustrated Heisenberg quantum antiferromagnet H 0 defined by J 1 = 1 and J 2 = 0.4.κ is given by Eq. 11 with λ chiral = 0.5.All time-dependent computations (shown as continuous curves) are done with N τ = 16 intervals in each drive period.(b) The GS of the Heisenberg quantum antiferromagnet H 0 (J 1 = 1, J 2 = 0) is used as initial state.κ is given by Eq. 11 with λ chiral = 0.7.(a,b) Results using the (weakly timedependent) effective Floquet Hamiltonian of Eq. 10 are shown as dashed-dotted curves.The (ideal) adiabatic limits are shown by green dotted lines.

Figure 7 :
Figure 7: Overlap | Ψ tramp |Ψ PEPS | between the time-evolved state at t = t ramp and the best optimized PEPS states, plotted versus t ramp .(a) J 1 = 1, J 2 = 0.4 and λ chiral = 0.5, and NN RVB state as initial state; (b) J 1 = 1, J 2 = 0 and λ chiral = 0.7, and GS of H 0 as initial state.Different symbols correspond to different PEPS manifold, as indicated in the legend.The overlap | Ψ tramp |Ψ init | with the initial state is shown as a reference.The adiabatic limit expected for t ramp → ∞ is shown on the right of the picture.

Table 1 :
the optimal PEPS overlaps on a 4 × 4 torus with the adiabatic GS considered in this work.Number of A 1 and A 2 tensors and number of variational parameters (fixing one tensor amplitude) used to construct the chiral PEPS, as a function of the bond dimension.The last lines show the corresponding overlaps with the two adiabatic GS considered in this work.B Derivation of the Floquet effective Hamiltonian B.1 H 0 = 0 case