Manipulating synthetic gauge fluxes via multicolor dressing of Rydberg-atom arrays

Arrays of highly excited Rydberg atoms can be used as powerful quantum simulation platforms. Here, we introduce an approach that makes it possible to implement fully controllable effective spin interactions in such systems. We show that optical Rydberg dressing with multicolor laser fields opens up distinct interaction channels that enable complete site-selective control of the induced interactions and favorable scaling with respect to decoherence. We apply this method to generate synthetic gauge fields for Rydberg excitations where the effective magnetic flux can be manipulated at the single-plaquette level by simply varying the phase of the local dressing field. The system can be mapped to a highly anisotropic Heisenberg model, and the resulting spin interaction opens the door for explorations of topological phenomena with nonlocal density interactions. A remarkable consequence of the interaction is the emergence of topologically protected long-range doublons, which exhibit strongly correlated motion in a chiral and robust manner.

Arrays of highly excited Rydberg atoms can be used as powerful quantum simulation platforms. Here, we introduce an approach that makes it possible to implement fully controllable effective spin interactions in such systems. We show that optical Rydberg-dressing with multicolor laser fields opens up distinct interaction channels that enable complete site-selective control of the induced interactions and favorable scaling with respect to decoherence. We apply this method to generate synthetic gauge fields for Rydberg excitations, where the effective magnetic flux can be manipulated at the single-plaquette level by simply varying the phase of the local dressing field. The system can be mapped to a highly anisotropic Heisenberg model, and the resulting spin interaction opens the door for explorations of topological phenomena with nonlocal density interactions. A remarkable consequence of the interaction is the emergence of topologically protected long-range doublons, which exhibit strongly correlated motion in a chiral and robust manner.
Rydberg atoms, held in optical lattices or arrays of optical tweezers, are currently among the most promising and versatile simulation platforms for quantum magnetism [27][28][29][30][31][32][33][34]. Realizing an effective magnetic flux for Rydberg excitations requires a complex-valued exchange interaction between different atomic states. This can be accomplished via a proper tuning of dipolar exchange interactions between different Rydberg states [35]. By applying a real magnetic field, the intrinsic spin-orbit coupling of such interactions [36] can be used to generate an effective spin-lattice model with a nonvanishing Peierls phase that emerges perturbatively to leading order in the strength of the dipole-dipole interaction. The possibility to engineer synthetic gauge fields in interacting spin lattices then holds exciting perspectives for exploring exotic many-body phases and dynamics [37][38][39][40][41].
Here, we describe an approach to implementing synthetic gauge fields via photon-assisted excitation exchange in Rydberg-atom arrays. Dressing the atoms by multicolor laser fields is shown to offer complete optical control of the generated spin-exchange interactions at the level of individual sites. The induced Peierls phase is determined by the relative phase between the applied laser fields, with which the effective gauge flux can be tuned to arbitrary patterns. Together with the nonlocal density interaction between Rydberg excitations, this yields a versatile quantum simulation approach for exploring strongly correlated systems with nontrivial band topologies and finite-range interactions. We illustrate these perspectives by discussing the topologically protected chiral motion of bound pairs [42][43][44] that emerge from the multicolor Rydberg dressing of an atomic square lattice [45].
We consider a two-dimensional array of individually trapped atoms. For a monochromatic Rydberg dressing [46], the ground state |g i of the ith atom (located at r i ) is coupled to a Rydberg state |r i by an off-resonant laser field with Rabi frequency Ω i = |Ω i |e iϕi and detuning ∆ = ω 0 − ω L , where ω 0 and ω L are the atomic transition frequency and laser frequency, respectively. Let us first consider a two-atom model [see Fig. 1 (a)]. In this case, the singly excited pair states |g i r j and |r i g j are degenerate and coupled by two Raman pathways with intermediate states |g i g j and |r i r j [47,48]. Importantly, the contribution from both pathways is not symmetric because the doubly excited state |r i r j is shifted by the van der Waals (vdW) interaction V ij between Rydberg atoms. In the limit of large detunings, |Ω i(j) | ≪ |∆|, |∆ + V ij |, this asymmetry leads to an effective hopping of the Rydberg excitation from site j to site i with a strength Taking into account the phase of each dressing field, this hopping amplitude J ij = |J ij |e iφij becomes complex, with a Peierls phase φ ij = ϕ i − ϕ j given by the sum of the phases ϕ i and -ϕ j of both involved laser fields. However, it is impossible to induce a nonvanishing effective magnetic flux Φ because the same driving field creates and annihilates an excitation on each site, leading to a vanishing net flux as the excitation circulates around a closed loop. This is illustrated in Fig. 1(c) for three atoms whose total flux Φ = (ϕ 2 − ϕ 1 ) + (ϕ 3 − ϕ 2 ) + (ϕ 1 − ϕ 3 ) = 0 vanishes identically.
As we shall show below, a finite magnetic flux can, however, be realized by applying multiple dressing fields with different frequencies (colors). The underlying mechanism exploits the frequency sensitivity of the stimulated excitation exchange, which is only resonant if the two in-volved atoms share a set of dressing fields with identical colors. As illustrated in Fig. 1(d) for a three-atom plaquette, this makes it possible to independently control the phase of the exchange interaction between each atom pair with only three different colors A, B, C of the dressing fields applied, i.e., driving each atom with two different colors opens up a distinct interaction channel for each atom pair. The phase of each laser field can be tuned individually, such that we obtain a finite and tun- Here, we have labeled the dressing fields of different colors by Θ ∈ {A, B, · · · } and denote their phases at a given site i by ϕ iΘ .
We can gain more intuition by considering two atoms that are each dressed by fields with distinct detunings ∆ A and ∆ B [see Fig. 1 (b)]. This implies a finite energy defect δE = ∆ A − ∆ B for the excitation transfer between the states |g i r j and |r i g j . The stimulated emission and reabsorption of photons is, thus, off-resonant by |δE| and will suppress the excitation hopping if |δE| ≫ |J ij |. Under this condition we can introduce distinct interaction channels that only emerge when a pair of atoms is dressed by laser-field with the same color Θ [ij] [49], e.g., Θ [23] = {B, C} ∩ {C, A} = C for the configuration considered in Fig. 1(d). To verify that the crosstalk between different colors can indeed be suppressed, we calculate the dynamics governed by the exact Hamiltonian for two dressed atoms. Here, describes the time-dependent single-particle Hamiltonian of the ith atom, and V ij = C 6 /|r i − r j | 6 denotes the vdW interaction between Rydberg excitations. As shown in Fig. 2(a), the excitation transport is nearly perfect for monochromatic dressing. The inclusion of two additional fields with different colors B and C, does not perturb this transport as long as the corresponding energy offset With feasible experimental parameters [50], this condition can be well fulfilled, while realizing a flux of Φ = π/2 for the three-atom case considered in Fig. 1(d). This is verified by Fig. 3(a), where we find a characteristic chiral motion of the Rydberg excitation around the plaquette based on exact simulation of the dynamics. The simulation is in good agreement with an effective model that neglects any crosstalk between different colors.
We can formulate the corresponding effective Hamiltonian by a perturbative analysis using Floquet theory [17,[51][52][53][54][55][56][57][58]. Introducing bosonic creation operators b † i = |r i g i | with the hard-core constraintb † ib † i = 0, we obtain a single-excitation effective Hamiltonian is the hopping strength of the excitation, and µ i is the on-site chemical potential. All these parameters can be tuned individually by adjusting the dressing fields at each site [49]. Importantly, the effective magnetic flux can be tuned by controlling the phase distribution of the applied dressing fields and does not depend on the laser intensities, lattice configurations, or Rydberg interactions. This makes it possible to generate an arbitrary magnetic flux pattern while implementing the effective Hamiltonian to an arbitrary accuracy by independently suppressing state leakages and crosstalk errors. Furthermore, the U(1) symmetry of the effective Hamiltonian Eq. (3) protects the system against various decoherence processes, e.g., global laser phase noise induced dephasing is eliminated in the resulting decoherence-free subspace [49, 59,60], and the damping caused by decay of the Rydberg states can be mitigated by post-selection as described in the Supplemental Material [49]. In the Supplemental Material [49], we show that the challenging single-site control of the laser intensity and phase distribution can be achieved by a compact optical module, whose complexity does not grow with system size. We also provide explicit schemes for manipulating gauge fluxes in different lattice geometries, including a square lattice, a triangular lattice, and a honeycomb lattice [49]. Altogether, this offers an accurate and scalable approach to synthesizing gauge fields with promising coherence properties. Let us now discuss the behavior of multiple excitations, which feature strong vdW interactions between them. If the initial distance between excitations is sufficiently large such that their interactions V ij are smaller than the detunings [49], the simple Hamiltonian describes the many-body dynamics of the system. Such a Hamiltonian is equivalent to a highly anisotropic Heisenberg model, as the density interaction V ij is two orders of magnitude larger than the hopping strength |J ij |.
The nonlocal density interaction V ij also facilitates the study of many-body dynamics beyond the hard-core constraints. To verify the effective model and its distinct features, we first consider to implement a Harper-Hofstadter ladder. With dressing fields of four different colors introduced in Fig. 3 . We also note that the trajectories for two excitations are no longer symmetric when nonlocal interaction V ij comes into play [61,62]. Such an interesting feature does not appear in systems with only hard-core constraints [Eq.
(3)]. The ability to implement such effective models permits exploration of topological phenomena in the presence of strong and finite-range interactions. As an example [ Fig. 4(a)], we consider two-excitation dynamics in an anisotropic Harper-Hofstadter lattice [63], where the nearest-neighbor hopping take on different strengths J x and J y along x and y directions, respectively. The model can be realized with the similar scheme shown in Fig. 3(b). When the flux Φ is uniform and rational (Φ/2π = p/q), the single-particle spectra split into q gapped topological bands linked by gapless edge modes. Remarkably, the finite-range density interaction can further lead to the emergence of chiral edge-mode bound states with large bond lengths, as we will show below.
Let us consider first the case of two tightly bound excitations, separated by two lattice sites along the y direction as illustrated in Fig. 4(a). When the differences center-of-mass dynamics of the bound pair [ Fig. 4(b)], which moves across the lattice with modified hopping strengths J ′ Most importantly, the effective magnetic flux Φ ′ = 2Φ of the doublon increases to twice the single-particle flux. An interesting case emerges for Φ = 2π/3 such that Φ ′ = 4π/3 (mod 2π) = −Φ. In this case we expect to observe the co-existence of single-excitation chiral edge state and topologically protected edge-mode bound state with opposite chirality.
When the above interaction difference is finite, the situation becomes complicated and an exact diagonalization is required [64,65]. Since the effective Hamiltonian is verified in Fig. 3(d), we will perform the calculation based on Eq. (5). For an infinitely extended lattice along the y direction with a finite number of sites along the x direction, the two-body eigenstate with a y-direction center-of-mass quasimomentum K can be written as where r = y 2 −y 1 and R = (y 1 +y 2 )/2 denote the relative and center-of-mass coordinates along the y direction,b † rν creates a Rydberg excitation at r ν = (x ν , y ν ), and |0 is the vacuum state with all atoms in |g . The associated spectrum of eigenenergies E K is shown in Fig. 4(c), and reveals a number of interesting states. The scattering continuum forms the lowest band with oscillating density of states. Above the scattering continuum, we identify different patterns of bound states. Figures 4(c) and 4(d) show the dispersion and density profile of three types of these states. The type-I bound state corresponds to the one depicted in Fig. 4(a), whose spectrum is split into three energy bands with Chern numbers C = {−1, 2, −1} (from the lowest to the highest) predicted by the centerof-mass motion analysis.
The state, marked as A in Fig. 4(c,I), is located in the upper band and represents a normal bound state in the bulk of the system. The other two states, marked as B and C, lie within the gap between the lowest and the middle bands and are topologically protected bound states that are respectively localized at the right and left edge of the lattice. Similar states can be identified for the type-II bound pair [ Fig. 4(c,II)], which are aligned along the x direction, and form bulk-mode (marked as D) as well as chiral-edge-mode (marked as E) bound states. The finite range of the Rydberg-state interaction also makes it possible to form bound pairs separated by a longer distance, such as type-III state, indicated in Fig. 4(c,III) and shown in panel F of Fig. 4(d).  Fig. 5(b). One clearly finds unidirectional doublon motion that is robust against local defects, and now features opposite chirality compared to the single-excitation transport.  In conclusion, we have described the application of multicolor laser-dressing to generate synthetic gauge fields for Rydberg excitations in neutral atom quantum simulators. Our approach makes it possible to realize arbitrary Peierls phases without compromising the accuracy of the quantum simulation. In particular, the scheme features individual tunability of the magnetic flux, which permits the experimental study of rich physics not yet explored with a non-uniform gauge field, e.g., exotic metal-insulator transition driven by the random flux [66][67][68], and emergence of anyons from a fractal flux pattern [69]. The inherent long-range interacting feature also motives future work on topological phenomena in interacting systems, such as emergent dynamical gauge fields [39,70] and fractional quantum Hall physics [71].
Note added. We recently became aware of two interesting schemes for realizing synthetic gauge fields in Rydberg-atom arrays, respectively based on laserassisted dipole-dipole interactions in a rectangular lattice [72] and Rydberg-dressing induced ground-state interactions in a honeycomb lattice [73].
The Rabi frequencies are chosen according to such that the nearest-neighbor hopping take the same strength. In order to balance the on-site potential µi, the detuning of each field is slightly shifted at each site [49].
The center-of-mass trajectories for the left and right regions in (d) are defined asrL = x≤0 rPr/ x≤0 Pr and rR = x≥0 rPr/ x≥0 Pr, with r = (x, y) and Pr the time-dependent occupation probability at the site r. This Supplemental Material provides additional information for the main text. In Sec. I, detailed derivations of the effective Hamiltonians in different regimes are given by using quantum Floquet theory and perturbation analysis. In Sec. II, we illustrate how to adjust chemical potentials and the non-nearest neighbour hoppings. In Sec. III, we systematically study the decoherence of our scheme. In Sec. IV, a feasible experimental implementation is designed for simulating the Harper-Hofstadter model and we give some considerations about technical challenges in the future.

I. EFFECTIVE HAMILTONIAN
In this section, we provide a Floquet analysis [S1-S4] of our system. By using the generalized van Vleck (GVV) nearly degenerate perturbation theory [S5-S7] to the infinite-dimensional Floquet matrix, time-independent effective Hamiltonians in different regimes are derived.

A. Floquet formulation
Quantum Floquet theory can be applied to quantum systems governed by time-periodic HamiltonianĤ(t) = H(t + T ), that satisfies time-dependent Schördinger equation which possesses a solution in form of ψ(t) = e −iεt/ φ(t), where the so-called Floquet mode φ(t) exhibits the same time periodicity asĤ(t) and ε is the corresponding quasienergy. The eigenvalue problem for quasienergy is given by which can be transformed into an equivalent time-independent problem, with the help of an infinite-dimensional Floquet matrix in an extended Hilbert space [S1, S4]. Following the standard approach, we expand the Hamiltonian H(t) and quasienergy eigenfuction φ(t) into the Fourier components of ω = 2π/T , The extended Hilbert space is then spanned by a complete set of orthonormal basis |α, n = |α ⊗ |n , where |α is a complete set of orthonormal basis states of the original physical system and n denotes the Fourier index ranging from * These authors contributed equally to this work. −∞ to ∞. The matrix element of the time-independent Floquet HamiltonianĤ F can be calculated as where = 1 is assumed unless otherwise noted. The block structure of the Floquet matrixĤ F is shown in Fig. S1, with each block of the same dimension as the physical system, and a physical state belonging to index n can be coupled to another state belonging to index n ′ by matrix elements of blocksĤ ±(n ′ −n) .

B. The generalized van Vleck nearly degenerate perturbation theory
If the Floquet matrix can be divided into an unperturbed partĤ 0 and a perturbation λV, i.e.Ĥ F =Ĥ 0 + λV, the corresponding effective Hamiltonian and its eigenstate solution in the concerned subspace can be expanded in powers of parameter λ,Ĥ Without loss of generality, we set λ = 1 in the following analysis.
For our system,Ĥ 0 =Ĥ eff includes all diagonal terms derived from the transition energy of Rydberg state as well as interactions between Rydberg atoms. Its eigenstates Ψ (0) form the orthonormal basis |α, n described above, which can also be considered as 0th-order eigenstates. The perturbationV represents the remaining nondiagonal coupling terms from driving fields. Thus, employing the GVV method, the 1st-order correction eigenstates for |α, n can be expanded into zeroth-order eigenstates as in the following with 0th-order eigenenergy E |α,n = E α + nω, where E α denotes energy of the physical state |α . Then the lowest few high-order effective Hamiltonian can be computed as follows [S4]:

C. The notations of colors
All colors labeled by Θ in the laser system form a set {Θ|Θ = A, B, · · · }, whose subset {Θ i } ⊆ {Θ} includes colors of dressing fields on specific atom i. Besides, We use Θ [ij] to denote the common colors between two atoms, which also form a subset

D. The full time-dependent Hamiltonian
Under the rotating wave approximation (RWA), the full time-dependent Hamiltonian of our system can be first written asĤ with ω 0 the atomic transition frequency, Ω iΘ and ω Θ the Rabi frequency and laser frequency of dressing field on atom i with color Θ, and V ij the van der Waals interaction between two Rydberg atoms. By using the unitary transformation U = i e −iω0t|ri ri| , the above time-dependent Hamiltonian will reduce toĤ full (t) in the main text.

E. Single-excitation
Though the dimension of the whole Hilbert space L is large, we only need to concern a few relevant subspaces if the excitation number is conserved. For single-excitation, three subspaces should be considered: ground state |G = |g 1 g 2 · · · g N , single-excitation subspace Π 1 = {|Ψ i = |g 1 g 2 · · · r i · · · g N , i = 1, 2, · · · , N }, and double-excitation subspace Π 2 = {|Ψ ij = |g 1 · · · r i · · · r j · · · g N , 1 ≤ i < j ≤ N }. Because the dynamics we focus on occur in the single-excitation subspace, the final effective Hamiltonian can be expressed by states inside Π 1 . The Hamiltonian can be first written aŝ Then, we assume there exist an eligible elementary frequency ω that all laser detunings ∆ Θi in the Hamiltonian can be expreessed as an integer multiple, i.e. ∆ Θi = n Θi ω. Thus, all Fourier componentsĤ n can be found according to Eq. (S3). In the extended Hilbert space, the Floquet matrix can be expressed as diagonal termsĤ 0 and nondiagonal As can be seen, the operatorV −nΘ i (V +nΘ i ) has two effects, one of which is deexciting (exciting) the atom on site i, and the other is to shift the Fourier index by −n Θi (+n Θi ). According to Eq. (S6), 1st-order wavefuction for |Ψ i , N can be calculated as with N an arbitrary integer. Equation. (S7)-(S9) can then be used to calculate the effective Hamiltonian in the extended Hilbert space. It is not difficult to find that the sole nonvanishing term is Ψ (0) |V|Ψ (1) inĤ (2) eff , which contains both diagonal and nondiagonal terms Equation.
(S16) labels the chemical potential correction to an excitation on site i, while hopping dynamics can be understood by further analyzing Eq. (S17). If there exists a channel Θ [ij] between atoms i and j, i.e. n Θi = n Θj ≡ n Θ [ij] , states |Ψ j , N and |Ψ i , N within the same block (labeled by N ) can be resonantly coupled by . (S18) When such a channel is well-defined, we can safely ignore all other nonresonant couplings between different blocks through different colors due to the large energy defect δE = ∆ Θj − ∆ Θi between states |Ψ j , N − n Θi + n Θj and |Ψ i , N . Now, the effective Hamiltonian of single-excitation subspaces within one block can be expressed aŝ which reduces toĤ (1) eff in the main text.

F. Double-excitation
For the case of double-excitation, the subspaces considered include single-excitation subspace Π 1 , double-excitation subspace Π 2 , and triple-excitation subspace Π 3 . The Hamiltonian can be expressed aŝ with all unperturbed diagonal terms in the first line, and the perturbed nondiagonal coupling terms in the second line. The 1st-order correction to the wavefunction |Ψ ij , N can be calculated as Again, the nonvanishing diagonal and nondiagonal terms can be obtained from Ψ (0) |V|Ψ (1) inĤ (2) eff , In fact, states in Π 2 are not quasi-degenerate for all ij due to the existence of vdW interaction V ij . The doubly excited subspace Π 2 can be approximately decomposed into Π 2 = Π ′ 2 ∪ Π ′′ 2 , with Π ′ 2 spanned by the dimer states and Π ′′ 2 the complementary set of Π ′ 2 .
If the initial two excitations i and k are nearest-neighbors (NN), they will be bound together, forming a dimer. The dynamics for this dimer can be described as follows: when one of the excitation i is viewed as fixed, its NN excitation at k can hop to another NN site j (i.e V ij = V ik = V NN ) according to Eq. (S23) as long as there is a channel Θ [jk] connecting them. From Eq. (S23) we can work out this effective exchange interaction as J indicating that excitation k can hop to site j with the help of an auxiliary excitation i. It is noted that this dynamics is sensitive to distances between adjacent atoms because the energy defect of this hopping process is δE ≈ V ij − V ik , thus even a few percent deviations of equal distance can destroy the quasi-degenerate condition.

2
Initially, if two excitations are distant enough such that V ij ≪ ∆, excitations hop to their nearby sites with a strength following the single-excitation condition. However, the presence of a distance-dependent nonvanishing V ij terms in the diagonal prevent the two excitations coming too close, as if there were a repulsive interaction between excitations. In fact, this interaction further distinguishes different subspaces inside Π ′′ 2 . The effective Hamiltonian under this condition becomesĤ (2) eff =Ĥ (1) eff +Ĥ int , withĤ (1) eff the single-excitation effective Hamiltonian, andĤ int = i<j V ijb † ib † jb jbi the nonlocal excitation-excitation interaction.

A. Chemical Potentials
The chemical potential of every site can be adjusted to a more favorable value µ ′ i = µ i −δ i by changing the detunings of laser fields on this site slightly ∆ ′ Θi → (∆ Θi + δ i ), as verified in the case of three atoms [see Fig. S2(a)]. By setting R = 5 µm and keeping other parameters the same as every in the main text, chemical potentials for three sites become significantly different, leading to imperfect chiral motion for a single Rydberg excitation. After adjusting detunings of laser fields on atoms 2 and 3 by δ i = µ i − µ 1 (i = 2, 3) to balance their chemical potentials to the same reference value µ ref = µ 1 , the perfect chiral motion is restored. In practice, it is not as convenient to directly change the frequencies of dressing fields individually. A more feasible method could employ additional far off-resonant addressing laser beams to introduce local energy shift δ i on the specific sites i [S8].

B. The non-nearest neighbouring hoppings
Due to the limited resources available for laser frequency modulation, we should introduce as few colors as possible to assist with experimental implementation. At the same time, this consideration could cause non-nearest neighbor hopping when the number of sites exceeds the number of colors. Fortunately, the strength of non-nearest neighbor hopping relative to the nearest neighbor hopping can be tuned by adjusting distances between atoms while maintaining the same configuration. For example, if we use three colors to connect six atoms [see Fig. S2(b)], four types of hopping arise between atoms: the nearest hopping J 1 (atom 1 ↔ 2, . . . , etc.), the diagonal hopping 2J 2 (atom 1 ↔ 4, . . . , etc.), the next nearest neighbour hopping J 3 (atom 2 ↔ 6, . . . , etc.), and the longer-range hopping J 4 (atom 1 ↔ 6.). Their respective strengths are plotted as a function of lattice spacing d in Fig. S2(b), showing that a convergent result J int = Ω 2 /4∆ when d is small, and J 1 dominates for sufficiently large d.

III. ROBUSTNESS AGAINST DECOHERENCE
In this section, we systematically study the robustness of our scheme against several major decoherence sources, including laser phase noise, thermal Doppler broadening, and decay of the Rydberg state. We show that the effective decoherence rate γ eff is significantly suppressed compared with the bare decoherence rate γ of a single Rydberg qubit, which gives rise to a favorable scaling of the error per site. Monte Carlo simulations, i.e., modeling the noise φ(t) in the exact HamiltonianĤ full during the three-atom chiral dynamics [corresponding to Fig. 3(a) of the main text]. As shown in Fig. S3(a), assuming a large dephasing rate γ = 1MHz, the single-qubit Rabi oscillation rapidly damps over a 10 µs timescale. Nevertheless, the chiral dynamics governed by the effective spin-exchange interaction remain coherent over much longer time scale [see Fig. S3(b)]. The inherent coherence protection relies on the correlation between phase noises from the same laser source, overcoming the otherwise rapid diffusive chiral motion if each atom experiences an independent phase noise [see Fig. S3(c)]. For our multicolor scheme, driving fields of different colors might carry independent noise due to temperature drift or stress of different fibers [see Fig. S6 (b)], fields of the same color, however, always share a common phase noise. Since crosstalk between different colors is already suppressed by the secular dynamics due to the large detuning δ between different colors, the phase noise of the crosstalk is unimportant, and the dynamics remain well protected from the dephasing. We have performed calculations assuming independent phase noise for different color fields, and have confirmed that the coherence is protected [see Fig. S3(d)]. The oscillation damps even slower compared to the case of a global phase noise for all colors [ Fig. S3(b)], which we attribute to the suppression of crosstalk errors (micromotions) due to the uncorrelated phase noises. With the above coherence protection feature, the final decoherence rate is determined by a small bit-flip error rate ∼ (Ω/2∆) 2 γ for each atom, which can be mitigated by choosing a proportionally smaller dressing parameter Ω/∆.

B. Thermal Doppler brodening
The thermal motion of a trapped atom results in an uncertainty ∆ T (Doppler broadening) of the detuning for each experimental realization. For a single qubit, the Ramsey fringes rapidly damp on the time scale 1/∆ T [see Fig. S4(a) with ∆ T = 2π × 44 KHz drawn from Ref. [S10], corresponding to T = 10 µK]. In our scheme, such an uncorrelated dephasing creates a random on-site potential offset on the scale of ∆ T for the effective Hamiltonian. As a result, dephasing can be suppressed if the induced spin-exchange interaction J is larger than ∆ T . To understand this, we can first consider two-atom dynamics, where the eigenenergy of the system is given by ± J 2 + ∆ 2 T ≈ ±(J + ∆ 2 T /J), i.e., the effective dephasing rate is reduced to ∆ T × (∆ T /J). We have confirmed that coherence is maintained for a longer time scale by Monte Carlo simulations of the three-atom chiral dynamics in the presence of Doppler broadenings [see Fig. S4(b)]. For a larger system, the Doppler effect could cause localization of the Rydberg excitation, but the effective decoherence rate is suppressed as well. Assuming a simple Anderson localization, the localization length is on the scale of ξ ∼ J 2 /∆ 2 T . According to Ref. [S11], the dynamics remain ballistic until the localization length is reached. Therefore, the coherence time is approximately given by t c = ξ/J, which corresponds to an effective decoherence rate γ eff = 1/t c ∼ ∆ 2 T /J, in agreement with the two-site analysis. On the technical side, the Doppler broadening can be easily reduced by sideband cooling of the trapped atom, as has already been demonstrated [S12, S13].

C. Rydberg decay
For the Rydberg state |70S 1/2 considered in our paper, its lifetime can be up to τ r = 50 µs in existing setups [S10], which is sufficiently long for the observation of the dynamics induced by the proposed effective spin-spin interactions. The decay of Rydberg excitations can be further reduced by increasing the detuning from the intermediate state in two-photon excitation schemes or using single-photon transitions, to reach τ r = 146 µs at room temperature and τ r = 410 µs at low temperatures by suppressing black-body radiation [S14]. In the following, we show that quantum simulations of physical observables can be purified even when Rydberg state decay remains significant. The decay of Rydberg excitations can be described by the master equation where the coherent spin dynamics are governed by the effective HamiltonianĤ eff that conserves the total Rydberg excitationN r = iσ i rr , while the Rydberg decay causes quantum jumps which decreaseN r by one. Therefore, if we only consider dynamics in the absence of quantum jump, evolution of the system is governed by a non-Hermitian HamiltonianĤ NH = H eff − (iκ/2) iσ i rgσ i gr =Ĥ eff − (iκ/2)N r , whereN r can be replaced by a number N r because the total excitation is conserved [equivalent to the U(1) symmetry] during the no-jump dynamics, i.e., [Ĥ eff ,N r ] = 0. Therefore,Ĥ NH describes the same dynamics as the effective HamiltonianĤ eff after normalization. In practice, the above purification scheme can be implemented by post-selection, which is widely practiced in quantum optics: when evaluating a physical observable in a single experimental run, only measurement outcomes that conserveN r are registered as a successful event. The success probability of the post-selection scales as P = e −Nrκt , while the conditional simulation fidelity is limited by a small non-Hermitian part ∼ (Ω/2∆) 2 κ of the induced interaction J, which can be obtained by replacing ∆ with ∆ − iκ/2 in its expression. The performance of the post-selection is shown in Fig. S5, where a large decay rate κ = 0.  ←→ |nS 1/2 ) as is widely adopted. The 420 nm dressing beam illuminates atoms globally, while the addressing 1013 nm light comes with four distinguishable frequencies assuming different colors. (c) Numerical results of the 1013 nm dressing light for color D on the far-field focal plane. The normalized intensityĨ(x, y) and phase distribution ϕ(x, y) in the region of interest (ROI) are shown respectively.

IV. EXPERIMENTAL IMPLEMENTATION
In this section, we present a scheme for experimentally implementing a flux lattice with multicolor Rydberg dressing as well as some considerations about the techniques.
First, we show how to build up the Harper-Hofstadter lattice system considered in the main text. As illustrated in Fig. S6(b), four dressing colors {A, B, C, D} are created by shifting the frequency of the field from a common laser source with an individual acoustic optical modulator (AOM) and are sufficient to generate individually tunable magnetic flux for arbitrarily large lattice, as one tunes the relative phases between laser fields at different sites for only two colors. This can be achieved by imprinting a computer-generated hologram (CGH) on a spatial light modulator (SLM). To obtain a proper CGH, iterative Fourier transform algorithms can be employed. However, unlike the production of a dipole trap array, for which intensity distribution is the primary concern, we need to also control the phase distribution of the laser field to construct the desired flux lattice [S15, S16]. To achieve such a simultaneous control of the light intensity and its phase distribution in the focal plane, we resort to the conjugate gradient minimization protocol [S16]. As an example, we consider how to realize a uniform magnetic flux in the square lattice shown in Fig. S6(a), for which we choose the phases of the beam arrays in colors B and D to be linearly increasing along the horizontal direction, where Peierls phase corresponding to the nth row is nφ. The generated light pattern for color D is shown in Fig. S6(c), where the desired beam array can be precisely formed, with its phase well defined within the central region of each site (see the black dashed circles) and increases linearly as expected. Taking 87 Rb atom as an example [see Fig. S6(b)], two-photon dressing scheme couples the ground state |5S 1/2 to a Rydberg state |nS 1/2 via the intermediate state |6P 3/2 . Specifically, a 420 nm global laser beam couples the two lower levels (|5S 1/2 ↔ |6P 3/2 ), and four 1013 nm addressing beams couple the two upper levels (|6P 3/2 ↔ |nS 1/2 ). To form different channel s, frequencies of the addressing beams are shifted correspondingly before coupling into fibers. The trapping beams (not shown in figure) can be reused to slightly adjust on-site potential offsets, the working mechanism of which has been explained in Sec. II A.
In the end, we discuss several technical challenges and accuracy of the effective model. First, our scheme requires single-site control of both the laser strength and phase distribution, which needs to be sufficiently homogeneous at each individual site to give a precise hopping strength J and Peierls phase φ. As for laser sources, the number of colors does not grow with system size. For example, to create an arbitrary magnetic flux pattern in a square lattice of an arbitrarily large size, only four different colors are required, with each site addressed by two different colors [see Fig. S7(a)], while for a triangular lattice, the number of colors is reduced to three [see Fig. S7 monochromatic Rydberg dressing schemes. The accuracy of our scheme is limited by two factors: a finite dressing parameter Ω/∆ results in a small bit-flip error ǫ b ∼ (Ω/∆) 2 per site, while a finite detuning δ between different color channels incurs a crosstalk error ǫ c ∼ (J/δ) 2 between neighboring sites, where J ∼ Ω 2 /4∆ denotes the effective hopping strength. For the square lattice shown in Fig. S6(a), assuming color A, B, C, and D respectively have a detuning (∆, ∆+δ, ∆+2δ, ∆+3δ), then a balanced coupling J A = J B = J C = J D = J requires Ω A < Ω B < Ω C < Ω D . For a given strength J, we can obtain the laser intensity I Θ = ( /8παd 2 eff )Ω 2 Θ , (Θ = A, B, C, D) for each color, with α the Fine-Structure Constant and d eff the effective dipole matrix element for two-photon coupling. Then, the total laser power P tot ≈ πw 2 0 N (I A + I B + I C + I D )/2 can be roughly estimated, where w 0 denotes the beam waist of each individual addressing field and N is the size of the system. More specifically, we have i.e., P tot ≈ (w 2 0 N /2αd 2 eff )J 2 (8/ǫ b + 3/ √ ǫ c ). This relation shows the explicit dependency of the total laser power on the errors ǫ b and ǫ c , and can be used to roughly estimate the laser power budget for the experiment. We also note that the design in Fig. S6(b) does not fully utilize the input laser power, and requires a few AOMs and SLMs, the number of which scales linearly with the number of colors. Improvement to reduce the overhead is also possible, e.g., by using different regions of a single SLM to modulate different color fields that are generated from a single AOD and incident obliquely upon the atom array from different directions.