Neutral atom arrays driven into Rydberg states constitute a promising approach for realizing programmable quantum systems. Enabled by strong interactions associated with Rydberg blockade, they allow for simulation of complex spin models and quantum dynamics. We introduce a new Floquet engineering technique for systems in the blockade regime that provides control over novel forms of interactions and entanglement dynamics in such systems. Our approach is based on timedependent control of Rydberg laser detuning and leverages perturbations around periodic manybody trajectories as resources for operator spreading. These time-evolved operators are utilized as a basis for engineering interactions in the effective Hamiltonian describing the stroboscopic evolution.

As an example, we show how our method can be used to engineer strong spin exchange, consistent with the blockade, in a one-dimensional chain, enabling the exploration of gapless Luttinger liquid phases. In addition, we demonstrate that combining gapless excitations with Rydberg blockade can lead to dynamic generation of large-scale multi-partite entanglement. Experimental feasibility and possible generalizations are discussed.

Introduction.-Programmable quantum simulators provide unique insights into complex many-body systems. They can be used for explorations of strongly correlated quantum phases of matter [1][2][3][4][5], non-equilibrium quantum dynamics [6][7][8][9][10][11][12][13][14][15], many-body entanglement [8,[16][17][18][19][20], and quantum metrology [13,14,[21][22][23][24]. Neutral atom arrays are a promising approach to realizing programmable quantum simulators [10,[25][26][27], where tunable atom trapping geometry along with strong interactions resulting in Rydberg blockade allow one to generate a variety of strongly correlated spin models. Coherent laser excitation into Rydberg states generates dynamics within the accessible Hilbert space, similar to that provided by a global transverse field in the Ising model. These constrained dynamics result in new physical phenomena such as quantum many-body scars [10,12,28], which evade thermalization starting from certain product initial states. At the same time, extending the toolbox of Rydberg quantum simulation in the blockade regime to dynamical generators beyond simple transversal fields is an open challenge. This is important, for instance, for realizing spin liquids [3,[29][30][31] and lattice gauge theories [32][33][34][35][36], for steering the dynamics of many-body states via counterdiabatic terms [37][38][39], realizing new types of quantum optimization algorithms [40][41][42] and for recent efforts to generate metrologically useful entanglement in systems with quantum many-body scars [43,44].

Motivated by these considerations, in this Letter we introduce a technique for Floquet engineering [45][46][47][48][49][50][51][52][53][54] that employs time-dependent control to realize effective Rydberg-blockaded models with versatile interactions. While Floquet engineering is widely utilized for interacting spin systems [55][56][57][58][59][60], conventional techniques rely on local Pauli frame transformations that generally violate the blockade constraint. Our approach, illustrated in Fig. 1 (b), leverages driven, periodic many-body trajectories originally discovered in the context of stabiliz-ing quantum many-body scars [12,61,62]. The complex micromotion of these trajectories serves as a resource for programmable Hamiltonian engineering, since perturbations applied during the periodic drive act at stroboscopic times via an effective Hamiltonian generated by time-evolved operators [61]. The resulting class of timeevolved operators forms a basis for the realization of novel interactions with tunable coefficients, which are not accessible in the static native Hamiltonian. Using this approach, we show how blockade-consistent spin exchange interactions can be engineered, enabling the investigation of gapless phases with emergent particle number conservation [63]. Moreover, in this new regime, we demonstrate the dynamic generation of structured multi-partite entanglement from Néel (Z 2 ) product initial states and explain this effect in terms of the dynamics of domain walls.

Hamiltonian engineering.-The key idea of this work can be understood by considering driven PXP model illustrated in Fig. 1 (a), which describes a one-dimensional atom chain with periodic boundary conditions, driven into the Rydberg state with fixed Rabi frequency Ω under idealized nearest-neighbor Rydberg blockade, with time-dependent global detuning ∆(t):

Here the operators

Our objective is to realize off-diagonal number conserving processes beyond single spin flips, utilizing the intrinsic controls Ω and ∆(t). In a conventional (static) approach (∆(t) = ∆) to this prob- We consider the driven PXP model as an approximate description of Rydberg atoms in optical tweezers, with Rabi frequency Ω and τ -periodic time-dependent global detuning ∆(t). Our protocol consists of π-pulses that realize a many-body echo (green) and deliberately placed perturbations (red). The echo realizes evolution under HPXP for effective time -τ 4 ≤ t ′ ≤ τ 4 , returning to t ′ = 0 at stroboscopic times. (b) The resulting dynamics perturbs around the periodic trajectories of the many-body echo, illustrated as micromotion in the full blockade-constrained Hilbert space. At stroboscopic times nτ , an effective Floquet Hamiltonian HF generates evolution within a constant energy submanifold. In our main application, the stroboscopic dynamics evolves a Néel (Z2) initial state towards a highly entangled GHZ state. Inset: Micromotion (black) and stroboscopic dynamics (red) of the Rydberg density ⟨n⟩ on a L = 16 periodic chain. Despite large oscillations within Floquet cycles, the stroboscopic evolution approximately conserves ⟨n⟩. (c) Our approach realizes effective models HF with tunable control over the dynamics of domain wall excitations on top of the Néel order. This includes a chemical potential J, blockade-consistent spin exchange interactions h that act as two-site hopping for domain walls, as well as creation/annihilation terms g of domain wall pairs. The regime of small g provides a mechanism for generating long-range multi-partite entanglement via growing superpositions of alternate Z2 orders.

lem, one typically relies on large detuning ∆ ≫ Ω, where multi-body spin flips emerge perturbatively in powers of (Ω/∆) n . The relative strengths of such processes is generally weak and cannot be tuned independently.

Our dynamic protocol circumvents this restriction by modulating laser detuning ∆(t) and leveraging a many-body spin echo. Specifically, since H PXP anticommutes with the operator i σ z i = e iπN , π-pulses of the global detuning effectively reverse its sign [28,61]: e iπN H PXP e iπN = -H PXP . This property enables a dynamical decoupling in the strongly interacting PXP model by means of a simple pulse sequence (n ∈ N),

Within each Floquet period τ , the system evolves forward under H PXP for τ /4, backward for τ /2, and forward again for τ /4. At a given time t ≤ τ , the system has thus undergone an effective evolution by H PXP for a time t ′ = t ′ (t) = ||t-τ /4|-τ /2|-τ /4, where t ′ ∈ [-τ 4 , τ 4 ], as shown in Fig. 1 (a) (see Supplemental Material [64], Sec. I A). As t ′ (τ ) = 0, the system exhibits periodic revivals at stroboscopic times nτ .

To generate a nontrivial effective evolution, we utilize complex micromotion along the periodic trajectory, as illustrated in Fig. 1 (b). Specifically, we introduce global detuning perturbations consisting of τ -periodic discrete pulses around the echo protocol of Eq. (3), ∆(t) = ∆ 0 (t) + ∆(t) with ∆(t) = j,n ∆j δ (tt jnτ ). At stroboscopic times, these perturbations translate into evolution under a static, local effective Floquet Hamiltonian H F , which holds up to an exponentially long prethermal timescale T p ≳ (τ /| ∆|)e cp/| ∆| , where | ∆| = j | ∆j | and c p > 0 [65,66]. In the interaction picture with respect to the perfect echo evolution for ∆ = 0, we obtain the leading contributions to H F through a Floquet-Magnus expansion [61,67] (see [64], Secs. I B-I C), resulting in

We note that H

F are constructed from Rydberg number operators conjugated by evolution under H PXP for the effective times t ′ j = t ′ (t j ) of the echo protocol, Ñ (t ′ ) ≡ e it ′ Ω 2 HPXP N e -it ′ Ω 2 HPXP . Consequently, operator spreading under the micromotion generated by H PXP induces interactions in the form of n-nested commutators (t ′ ) n [H PXP , ..., [H PXP , N ]...] in Ñ (t ′ ) and thus H F . Coefficients of these terms in H F are controlled by the locations t j and weights ∆j of the pulses. Pulses with t ′ j = 0 couple to the bare Rydberg number operator N . Further, symmetric weights of pulses at ±t ′ j ensures that all terms containing an odd number of commutators vanish in H (0) F , which thus conserves Rydberg number parity. Hence, for small t ′ j the leading non-trivial contribution to H (0) F appears at order (t ′ j ) 2 , which contains nearestneighbor pair flips. Single spin flips appear only in H (1)

Based on this intuition, we introduce the following detuning perturbations, parameterized in terms of the (dimensionless) variables γ, θ, ϵ, see Fig. 1 (a),

Here, the first two pulses both occur at effective time t ′ = 0, while the latter occur at t ′ = ± τ 4 . Inserting into Eq. ( 4) and expanding Ñ (τ /4) perturbatively for short periods Ωτ 4 ≪ 1, we obtain an approximate closed-form expression for the effective Hamiltonian,

Here, we have kept all terms to quadratic order in Ωτ, ϵ, γ, θ (see details in [64], Sec. I C);

Our Floquet protocol provides flexible relative tunability of the coefficients J, h, g by controlling the period Ωτ and parameters ϵ, γ, θ. Thus, our construction extends the capabilities of the Rydberg quantum simulator to blockade models with independent control over quasiparticle number conservation, motion, and creation/annihilation processes. In particular, this enables access to exchange-dominated regimes h ≳ g, which we explore in the following.

Controlled multi-partite entanglement.-We consider dynamics starting from a Néel state, |Ψ(t = 0)⟩ = |Z 2 ⟩. Excitations on top of this state can be viewed as domain walls, created in pairs by H PXP , see Fig. 1 (c); J sets a chemical potential and H PXYP generates (two-site) hopping of domain walls. We fix the parameters of the Floquet drive to Ωτ /2π = 0.77 (close to the scar period of the bare PXP model [10]), and ϵ = -0.45, γ = 1.0, θ = 0.15. This translates to an effective model of Eq. ( 6) with J ≈ 0.225, h ≈ 0.068, g ≈ -0.017, where the rate of domain wall hopping is stronger than of creation/annihilation. As shown in Fig. 1 (b), the average Rydberg density ⟨n(t)⟩ varies rapidly during micromotion, but evolves only slowly at stroboscopic times. Even though the parameters perturbing the echo are sizeable, the effective Hamiltonian Eq. ( 6) nonetheless provides a good description of the stroboscopic evolution as demonstrated in Fig. 2 (a). Interestingly, while the density of Rydberg excitations remains high, the staggered magnetization j (-1) j ⟨σ z j (t)⟩ washes out. As this happens, the system develops growing connected correlations

. Within the correlated region, a superposition of alternate Z 2 orders emerges, akin to patches of GHZ states. Once these large scale fluctuations reach the system size (here: L = 16), the state develops large overlap with an antiferromagnetic GHZ state, defined as

2 ⟩ , which we quantify via the fidelity max ϕ |⟨GHZ|Ψ(t)⟩| 2 [18,68], see Fig. 2 (c). As system size increases, the time of maximum GHZ fidelity changes roughly linearly in L, suggesting that the relevant processes are not exponentially suppressed. At the same time, the corresponding peak height decreases with L, indicating that correlations are less likely to reach the full system size.

These results can be understood based on the effective model Eq. ( 6), as sketched in Fig. 1 (c): gH PXP slowly creates a superposition of the initial Z 2 state with states containing pairs of domain walls, which then move rapidly via the strong hH PXYP interaction. The propagating pair carries a growing string of the alternate Z ′ 2 order, thus producing antiferromagnetic GHZ-like correlations. On a periodic chain, the pair may re-annihilate at the antipodal point to form the state |GHZ⟩. Based on this picture, we expect that the coherent spread of correlations persists over a timescale t * ∼ 1/g, beyond which it is interrupted by the emergence of additional domain walls. The size of the GHZ-like patch l * is determined by the rate of hopping v * ∼ h within this timescale, l * = v * t * ∼ h/g. For g ∼ 1/L, this size is l * ∼ L, and a full chain GHZ state may form in a time linear in system size.

We confirm these predictions using numerical analysis of the Quantum Fisher Information (QFI) density, which quantifies the metrological potential of the pure state |Ψ(t)⟩, evolving under H F , with respect to a stag-

. In particular, a QFI density F Q /L > m implies at least (m+1)body entanglement [69], such that F Q /L > 1 indicates non-classical correlations and F Q /L = L corresponds to a maximally entangled GHZ state, which saturates the Heisenberg limit. In Fig. 3 (b), iTEBD simulations of an infinite chain demonstrate that the QFI obeys the predicted scalings for the formation time t * and maximum 6) for g/h ≲ 1 weakly couples the Néel states to pairs of domain walls with momentum ±k and quasiparticle dispersion ε k = -2h cos k. Depending on the energy offset between Néel states and the dispersive band, this coupling can be on-or off-resonant. (b) Dynamics under HF for small g generates a maximum QFI density (see text) of order ∼ h/g within a time t * ∼ 1/g. This is seen through a scaling collapse for the dynamics of an infinite size chain simulated via iTEBD for different values of h, g at fixed J = 2h. (c) Phase diagram of HF at the integrable point g = 0, with gapped paramagnetic and Z2 phases, as well as a Luttinger liquid that is stable for K < 1/2 (in the presence of Z2 number parity), and unstable for K > 1/2 (towards a gapped topological state [63,70]). The Luttinger parameter K (black line) is obtained from the Bethe ansatz solution. (d) iTEBD simulation of the QFI density from a Z2 initial state for weakly broken integrability at small g. The QFI functions as a dynamical probe of the transition between Z2 phase and Luttinger liquid in c), and grows at a rate tied to the domain wall dispersion.

size l * of the multi-partite entangled regions.

Luttinger liquid dynamics-Due to the small value of g, the entanglement features observed in the previous section may be understood as a dynamical probe of the low energy properties of the constrained model Eq. ( 6) with U (1) symmetry at g = 0. In particular, H F | g=0 is known to be integrable [71], and single domain walls form a band of quasiparticles with dispersion ε k = -2h cos(k). The Néel Z 2 -states are offset from the center of this band by an energy J -2h due to chemical potential and H ZIZ term, see Fig. 3 (a). The ground state phase diagram of this model, which we calculate from Bethe ansatz integral equations in [64], Sec. II A, is shown in Fig. 3 (c). For J/h ∈ (-3, 6), the system is in a Luttinger liquid phase with gapless domain wall excitations. Outside this regime, the ground state transitions into a gapped paramagnetic (J/h < -3) or Z 2 (J/h > 6) phase.

From this phase diagram, we see that the protocol of the previous section corresponds to a quantum quench of the Z 2 state into the Luttinger liquid phase. For small g, the initial Z 2 state couples resonantly to pairs of domain walls with momenta |k, -k⟩ such that J -2h + 2ε k = O(g), which mediate the growing entanglement as described above. In particular, this coupling is proportional to the group velocity of the cosine dispersion, ⟨k, -k|gH P XP |Z 2 ⟩ ∼ 2g sin(k) (see derivation in [64], Sec. II B). For J/h outside the Luttinger liquid phase, domain walls can only be created virtually (i.e. off-resonantly), suppressing the growth of the QFI. As a consequence, dynamics of multipartite entanglement [72], quantified by the QFI density, directly probes the transition from the gapped Z 2 symmetry-breaking phase to the gapless Luttinger liquid, with its early time growth reflecting the quasiparticle group velocity. We confirm this prediction numerically by computing the dynamics of the QFI under H F of Eq. ( 6) using iTEBD for different values of J/h, see Fig. 3 (d). The growth of the QFI accurately captures the phase boundary at J/h = 6 and shows an enhanced rate towards the center of the Luttinger liquid phase. We have also verified these features in a direct simulation of the Floquet protocol for finite systems, see [64], Sec. III.

Robustness and Implementation in Rydberg arrays.-Successfully realizing and controlling H F in Eq. ( 6) relies on the expansions in small perturbations ϵ, γ, θ and drive periods Ωτ , which also tune the rate of dynamics. Thus, an optimal choice of drive parameters must consider the trade-off between a high-accuracy target Hamiltonian and the realistic constraint of observing significant dynamics within finite coherence time of an experimental device, an interplay we explore in [64], Sec. IV. In addition, in experiment, Rydberg atoms interact via van der Waals interactions

blockade radius R b . To demonstrate that our approach applies qualitatively, we construct a pure spin exchange model (-2ϵ = γ = 2θ ⇒ J = g = 0, h ̸ = 0) and numerically evaluate the quantum walk of a single Rydberg excitation resulting from Floquet evolution including V vdW at R b = 1.5. We further add a small constant detuning to mitigate the long range tail of V vdW [12] and use Gaussian pulse profiles with finite width, suitable for realistic control hardware with limited rate of detuning modulation; Fig. 4 (a). Comparison with the corresponding PXP result in Fig. 4 (b) shows very good qualitative agreement over an extended duration. Quantitatively, the effective hopping is even stronger in the experimentally relevant scenario, likely due to residual contributions from finite pulse width and long range tails. A detailed analysis of these effects is left for future work.

Discussion & Outlook.-We have introduced a Floquet protocol for systems of Rydberg atoms that exploits periodic trajectories of quantum states, enabling versatile Hamiltonian engineering. As an application, we realized models with emergent particle number conservation and dominant blockade-consistent exchange interactions in one dimension, exploring previously inaccessible gapless Luttinger liquid phases. In particular, we found that combining Rydberg blockade and gapless domain wall excitations leads to the generation of long-range, multi-

0.0 0.2 0.4 0.6 0.8 time Ωt/(2π) -8 -4 0 detuning ∆(t) (a) 0 20 40 0 4 8 12 site i n i 0 10 20 30 time Ωt/(2π) 0 4 8 12 site i 0.00 0.25 0.50 0.75 1.00 Rydberg PXP (b) partite entanglement upon evolving Z 2 product states.

Although our discussion focuses on one-dimensional systems, generalization to other geometries is natural. For instance, models akin to Eq. ( 6) are relevant to Rydberg spin liquids in two dimensions [30,31], as well as to achieving quantum speedup in combinatorial optimization by enabling delocalization in the adiabatic algorithm [41]. Moreover, going beyond short evolution times of the operator Ñ (t) provides access to even higher-body spin interactions, an approach we employ in Ref. [73] to study the dynamics of two-dimensional lattice gauge theories. We further emphasize that our Floquet scheme uses only simple global controls, but may be extended by incorporating site-resolved detuning fields [15], which could allow exploration of chiral interactions [74]. It would also be interesting to study the connection of our protocol with time-dependent methods for state preparation, such as counter-diabatic driving [75][76][77] and trajectory optimization [39]. Finally, we note that our Floquet protocol can be extended to other quantum simulation platforms, contingent on time-dependent control to engineer non-trivial periodic trajectories, as is available in dipolar interacting systems [55][56][57][58]60], trapped ions [78,79], neutral atoms [80,81], or superconducting devices [82][83][84].

ensures that all terms containing an odd number of commutators in Eq. (S11) vanish, which preserves Rydberg number parity. Moreover, for θ = γ/2, this symmetry persists at all orders of the Floquet-Magnus expansion, in a weakly rotated basis. Specifically at this point, ∆(t) becomes periodic in τ /2, and H(t) possesses period-τ /2 "twisted time-translation symmetry" with respect to the operator X τ /2 : H(t + τ /2) = X † τ /2 H(t)X τ /2 . Using the formalism of Ref. [86], the resulting Floquet unitary can be approximated as U (τ /2) ≈ VX τ /2 e -iDτ /2 V † , and

2 ≈ Ve -iDτ V † for some unitary frame transformation V perturbatively close to identity, and effective Hamiltonian D that commutes with X τ /2 = e iπN , i.e. has emergent Rydberg number parity symmetry. Indeed, at the θ = γ/2 point, we obtain

where the symmetrization Ñ τ 4 + Ñτ 4 ensures that D has parity symmetry at both zeroth and first orders (as well as higher orders not computed here). V implements corrections to this effective Hamiltonian beyond the zeroth order, thereby recovering the leading order contributions to H F in the original frame, computed in Eqs. (S18,S19).

In a regime of small Floquet periods Ωτ /4 ≪ 1, we may perform a perturbative expansion of the time-evolved Ñ (t) operators,

where H PYP ≡ i P i-1 σ y i P i+1 enacts blockaded local spin-flips with phase, H PXYP ≡

is the blockaded nearest-neighbor spin-exchange interaction, and H PZP ≡ i P i-1 σ z i P i+1 is a diagonal term that can be decomposed into the total Rydberg number operator and next-nearest neighbor ground-ground and Rydberg-Rydberg repulsion:

Inserting into Eqs. (S18,S19), we thus obtain a closed-form expression of the effective Floquet Hamiltonian:

by keeping terms up to quadratic order in Ωτ, ϵ, γ, θ.

The effective Hamiltonian H F described in Eq. (S22) features total Rydberg number conservation at g = 0, and is known to be integrable at this point [71]. Furthermore, Ref. [63] explored this U (1) symmetric model in the absence of the H ZIZ -term: As a function of the chemical potential J, the model exhibits gapped paramagnetic and Néel-ordered phases, as well as an intermediate gapless Luttinger liquid phase. The Luttinger liquid is robust to U (1)-breaking single spin fluctuations when the Luttinger parameter K < 1/8, and is stabilized by Z 2 parity symmetry for 1/8 < K < 1/2. For K > 1/2, the Luttinger liquid is unstable towards U (1)-breaking perturbations even in the presence of Z 2 symmetry and flows to a gapped topological phase.

Here, we adapt a similar analysis for H F | g=0 which additionally includes the H ZIZ term. Importantly, this term maintains integrability: Our system corresponds to a class of constrained XXZ models that Ref. [71] provides exact Bethe ansatz solutions for (concretely, t = 1 and ∆ = -0.5 in [71]), which we review here for completeness. As the total Rydberg number N ∈ {0, . . . , L/2} is a conserved quantity, each N -particle sector of the model can be diagonalized separately, by eigenstates labeled by quasi-momenta {k} = k 1 , k 2 , . . . , k N satisfying the following Bethe ansatz equations [71,87]:

1 + e ik l + e i(kj +k l ) , (S23)

and with the following total energy and momentum:

Low-energy properties of this model can be studied in the thermodynamic limit, through integral equations for quasi-momenta distributions specified by the Bethe equations. In particular, the Luttinger parameter K as a function of Rydberg density n 0 = N/L in the ground state reads:

where η(U ) and U 0 are determined by integral equations The Luttinger parameter K (Eq. (S26)), the ground state energy E|J=0 (Eq. (S30)) in the absence of chemical potential, and the chemical potential J (Eq. (S38)) as a function of ground state Rydberg density n0. These quantities are obtained from the numerical solution of the integral equations in Eq. (S22).

in Ref. [71]:

which are subject to the constraint

Moreover, the ground state Rydberg density n 0 at a given chemical potential J is determined by minimizing the ground state energy as a function of Rydberg density. The ground state energy E for our model is related to that of the unconstrained XXZ chain, Ẽ, provided in Ref. [88]:

where the factor of (1n 0 ) accounts for the reduced effective length of the chain due to the blockade constraint. This energy is minimized when

Employing Leibniz integral rule on Eq. (S29) gives

Furthermore, Ref. [88] derives

Combining these, we arrive at the ground state condition

where we define

and thus establish the relation between chemical potential J and ground state Rydberg density n 0 , as desired:

By numerically integrating these equations using the method of quadratures, as shown in Fig. S1, we obtain the Luttinger parameter K as a function of chemical potential J in Fig. 3 (c) of the main text. Along with the ground state Rydberg density n 0 , this enables characterizing the ground state phases of H F | g=0 , sweeping J/h: The model exhibits a gapped Z 2 phase with n 0 = 1/2 for J/h > 6 and a gapped paramagnet with n 0 = 0 for J/h < -3. Between -3 < J/h < 6, where 0 < n 0 < 1/2, a gapless Luttinger liquid emerges. Following the arguments of Ref. [63], the Luttinger liquid is stable to general U (1) symmetry breaking perturbations for K < 1/8 and stable to U (1) symmetry breaking perturbations that preserve Z 2 parity symmetry for 1/8 < K < 1/2. For K > 1/2, the Luttinger liquid is unstable towards such perturbations and flows to the gapped topological phase of the Kitaev chain [89].

In the main text, we considered quantum quenches of the Z 2 state into the Luttinger liquid, enabled by turning on a small U (1) symmetry breaking field g that couples the Z 2 state to other particle number sectors. This small field generates a low density of domain walls on top of the Z 2 state, and majority of dynamics can be studied in the 0-domain wall (Néel) and 2-domain wall sectors at sufficiently early times.

In order to study dynamics in this regime, we switch to a new description of H F | g=0 in Eq. (S22) in terms of bond variables, with domain walls on top of a Z 2 vacuum defined as the particle degrees of freedom. One subtlety of working with domain walls is to track the global Z 2 gauge degree of freedom: There are two distinct vacuum states |Z 2 ⟩ and |Z In the single-particle (single domain wall) sector, the Hamiltonian acts as

and can be exactly diagonalized using a plane wave ansatz

Accordingly, the single particle dispersion Ẽk is given by

To solve the two-particle sector, we use the Bethe ansatz

where we make use of the fact that domain walls come in even-odd pairs and their relative ordering is fixed. On an infinite chain, the eigenenergies Ẽk,k ′ with respect to H F | g=0,J=0 can be determined from the stationary Schrödinger equation for |uu ′ | ≫ 1, when the two domain walls are far apart. This results in

Once the chemical potential J is incorporated into H F | g=0 , the vacuum and two-domain wall state energies read:

With Ẽk,k ′ given by Eq. (S44), the scattering phase S(k, k ′ ) can be determined by projecting the stationary Schrödinger equation onto the state |u, o; u, e⟩, where two domain walls are located in the same unit cell u. Specifically,

Equipped with the two-particle eigenstates |k; k ′ ⟩, we introduce a weak U (1) symmetry breaking perturbation gH PXP and investigate the resulting dynamics starting from the |Z 2 ⟩ product initial state. In particular, the PXP perturbation connects |Z 2 ⟩ to the two-particle sector via

Moreover, as both the Hamiltonian and initial state are two-site translationally invariant, the only non-vanishing couplings are to states |k; k ′ ⟩ with k ′ = -k, which exhibit the same translational invariance. In particular, we consider the relevant matrix elements between |Z 2 ⟩ and the normalized two-domain-wall eigenstates N bc L |k; -k⟩, where N bc is an O(1) normalization constant that depends on the boundary conditions; N obc = 2, N pbc = √ 2. They are given by

and we note that λ(k) is proportional to the group velocity 2h sin(k) of the single domain wall dispersion. We thus see that the dynamics from |Z 2 ⟩ probes the lowenergy spectrum of H F and couples to the two-domainwall band at strength λ(k) and energy offset δ(k) ≡ E k,-k -E ∅ = J -2h -4h cos(k). When δ(k) ≲ λ(k), domain walls are created slowly but resonantly, generating coherence between the two different Z 2 orders and leading to the growing Quantum Fisher Information observed in Fig. 3 (d) of the main text.

As shown in the previous section and verified in Fig. 3 of the main text, the creation of multipartite entangle- 45, θ = 0.15 while varying γ, which effectively sweeps J for fixed h, g in the effective Hamiltonian Eq. (S22). Within a suitable range of γ values, dynamics starting from a Néel state generates domain walls that wash out the initial Z2 order and lead to a build-up of large Quantum Fisher information density. (d) Effective Hamiltonian coefficients J, h, g evaluated for these drive parameters as a function of γ. The Z2 state is resonant with a segment of the two-domain wall band within |J -2h| < |4h|, which indeed aligns with the range of γ values exhibiting growth of large multipartite entanglement.

ment in the dynamics of the effective Hamiltonian H F is due to a weak but resonant coupling between the initial Z 2 state and a low-energy two-domain-wall band. As such, the entanglement dynamics acts as a probe for the transition between a gapped Z 2 state and a gapless Luttinger liquid of domain walls in H F .

Here, we verify numerically that these features are indeed also present in the corresponding stroboscopic dynamics of the full Floquet time evolution, even for large perturbations around the many-body echo. Specifically, we consider the Floquet protocol with fixed parameters and J = g = 0, see Eq. (S22). The decay curves, which incorporate a phenomenological decay t * = 100Ω -1 , are presented as a function of (a) number of Floquet cycles, (b) physical time in units of the Rabi frequency, and (c) effective time in units of the hopping strength h of the effective Hamiltonian. The results demonstrate the tradeoff between robust drives with small period τ and perturbation ϵ, and efficient drives with larger perturbations that generate significant many-body dynamics within shorter amounts of physical time. Our definition for the effective coherence time incorporates both aspects. (d) Coherence time tc (see Eq. (S50)) of the Floquet protocol as a function of drive period τ and perturbation strength |ϵ| for the drive parameterization considered in (a-c). We estimate an optimal h × tc ≈ 10 in units of h. Red curve depicts a contour of constant h.

tion |Ψ(nτ )⟩ and target effective Hamiltonian evolution |Ψ F (nτ )⟩ for the quantum walk of a single Rydberg excitation. This is achieved by choosing the drive parameters as -2ϵ = γ = -θ, which generates a pure spin exchange model with h = -ϵΩ 2 τ 32 and J = g = 0 in (S22). Using state fidelity as a metric, we extract an effective coherence time t c from the decay |⟨Ψ F (nτ )|Ψ(nτ )⟩| 2 e -L(nτ /t * ) 2 ∼ e -L(nτ /tc) 2 , (S50) which incorporates a phenomenological decay constant t * = 15(2π/Ω) to model finite coherence time of an experimental device. (The choice of an overall Gaussian decay is motivated by the numerically evaluated fidelities displayed in Fig. S4 (a-c).) Then, we vary the pulse parameters τ and |ϵ| to maximize h × t c , which measures how much hopping occurs before the system decoheres (see Fig. S4 (d)).

Here, we highlight the tradeoff between robustness and efficiency inherent to our Floquet protocol: On one hand, S7 detuning profiles with small drive period τ and small perturbation |ϵ| result in a high fidelity between the stroboscopic dynamics and the evolution under H F . On the other hand, the corresponding hopping h ∼ ϵτ in H F is small, thus leading to slow dynamics that is eventually limited by physical coherence times. To demonstrate this tradeoff, we show the decay of the many-body fidelity of Eq. (S50) for varying drive parameters and different choices of units of time. As expected, smaller perturbations ϵ and Floquet periods τ result in a higher fidelity after a given number of Floquet cycles n, see Fig. S4 (a). However, when plotted against physical time in units of the Rabi frequency Ωτ n, the relative fidelities for drives with different periods τ are altered significantly, see Fig. S4 (b). Finally, identifying the time h×nτ in units of the hopping strength h of the effective Hamiltonian H F as the most relevant scale for the quantum simulation of H F , we see in Fig. S4 (c) that Floquet protocols with sizeable parameters ϵ, τ can outperform protocols with very weak perturbations around the manybody echo. Consequently, the Floquet coherence time h × t c extracted from Eq. (S50) exhibits a non-monotonic dependence on pulse parameters |ϵ| and τ , as shown in Fig. S4 (d), indicating an optimal choice of drive profile for the relevant experimental timescales.