High-dimensional SO(4)-symmetric Rydberg manifolds for quantum simulation

We develop a toolbox for manipulating arrays of Rydberg atoms prepared in high-dimensional hydrogen-like manifolds in the regime of linear Stark and Zeeman effect. We exploit the SO(4) symmetry to characterize the action of static electric and magnetic fields as well as microwave and optical fields on the well-structured manifolds of states with principal quantum number $n$. This enables us to construct generalized large-spin Heisenberg models for which we develop state-preparation and readout schemes. Due to the available large internal Hilbert space, these models provide a natural framework for the quantum simulation of Quantum Field Theories, which we illustrate for the case of the sine-Gordon and massive Schwinger models. Moreover, these high-dimensional manifolds also offer the opportunity to perform quantum information processing operations for qudit-based quantum computing, which we exemplify with an entangling gate and a state-transfer protocol for the states in the neighborhood of the circular Rydberg level.


Introduction
In the early days of quantum mechanics, Wolfgang Pauli proposed an elegant algebraic solution of the hydrogen problem based on the SO(4) symmetry of the Hamiltonian for an electron moving in a Coulomb potential [1]. This provided not only a derivation of quantized energy levels for Hydrogen described by the Rydberg formula E n = −Ry/n 2 , with n = 1, 2, . . . the principal quantum number, but also explained the n 2 -fold (orbital) degeneracy of the n-manifolds. In (weak) static external electric and magnetic fields, these degeneracies are lifted resulting in a linear Stark and Zeeman effect [2,3], leading to a large and regularly structured internal state space associated with each value of n.
Here we look at such high-dimensional and well-structured internal atomic state spaces, as provided by hydrogen-like Rydberg states of multi-electron atoms, as a novel opportunity for both storing and manipulating quantum information for quantum-simulation and computing. This involves the control of a given n manifold of single Rydberg atoms with external electromagnetic fields, as well as entangling operations via strong dipolar interactions. Availability and control of large internal Hilbert spaces is of interest in quantum simulation of general scalar field theories [4], gauge theories of high-energy physics [5], condensed matter models with large spins [6][7][8][9][10][11], or the emulation of synthetic spatial dimensions [12][13][14][15][16]. In view of the recent experimental advances with neutral atoms laser-excited to Rydberg states [17,18], we will explore some of these opportunities from a theory perspective and with Pauli's algebraic formalism providing a framework to treat the single-and many-body problem.
Rydberg atoms are already known as one of the most promising neutral-atom platforms for quantuminformation and simulation purposes in various scenarios [17][18][19][20]. These include the storage of qubits (or qudits) in ground states where fast and high fidelity entangling gates can be performed via interactions between atoms that are laser-coupled to Rydberg states [21][22][23][24][25][26][27][28][29][30][31][32]. Moreover, Rydberg excited atoms provide a direct realization of interacting spin-1/2 models by encoding a spin between a ground state and Rydberg  87 Rb, where energies corresponding to n = 21 are colored in blue. The large amount of discrete energy levels emerging from highly degenerate F = 0 manifold makes the spectrum appear continuous. Notice that for non-hydrogenic atoms low orbital angular momentum states are shifted away from the Stark manifold (d-states are shown). (c) Energy eigenvalues for the n = 21 manifold sorted with respect to the quantum number m l for a typical electric field value (indicated by the dashed line in (c) for F = F IT /2) F IT > F > 0 and parallel magnetic field B > 0, with ω Z = 0.1ωS, forming a rhombus shaped level structure. The manifold of states is described by two large angular momenta, where the application of raising operatorsĴ a,+ andĴ b,+ connects neighboring states along the diagonal direction indicated by the corresponding arrow. For non-hydrogenic atoms the highly regular spacing of the rhombus for |m l | < l * (blue shaded area) is disturbed as low angular momentum states are shifted to lower energies due to quantum defects. For 87 Rb the quantum defects are significantly nonzero for l ≲ 3 [63,64], and for sufficiently large electric fields as considered here, only s, p, d-states are effectively missing and therefore we have l * = 3. The bordered and shaded areas indicate different manifolds of states, namely H t (orange), H a (red) and H v (green), which are used in this paper for quantum simulation purposes. The inset is a zoom of the level structure where the eigenstates and energy differences ωa and ω b are explicitly indicated.
In contrast, we now investigate opportunities offered by the larger Rydberg n-manifold of states of alkali, alkaline earth and lanthanide atoms trapped in optical tweezers in the regime of a linear Stark and Zeeman effect, as illustrated in figure 1. We use the SO(4) Lie algebra represented by two large angular momentum operators, which we denote byĴ a andĴ b , to characterize the hydrogen-like high orbital angular momentum states and the action of the external electromagnetic fields. Within this framework, we show how to encode and manipulate quantum states in these high-dimensional Rydberg manifolds. Moreover, this algebraic representation in terms of two angular momenta allows us to translate the dipole-dipole interactions between Rydberg atoms into generalized (large angular momentum) Heisenberg models, for which we discuss state preparation and readout schemes.
By exploiting the mapping of large angular momentum operators into conjugated (continuous) variables, we show that such Heisenberg models can be used to simulate quantum field theories (QFTs). Our approach enables the simulation of far-from-equilibrium dynamics, which is notoriously difficult to simulate classically and largely inaccessible to implementations based on low-energy approximations. We illustrate our ideas with the sine-Gordon model in near equilibrium and its massive extension, which is dual to 1+1D QED, namely the Schwinger model in out-of-equilibrium situations. Furthermore, we also discuss opportunities for quantum information processing offered by the Rydberg n-manifolds. We exemplify this for states near the circular level by providing a state-transfer protocol that can be employed to realize entangling gates for qudit-based quantum computation.
A timely opportunity to implement the results of our work is provided by recent Rydberg experiments with open inner shell atoms such as Er [65]. First, the hydrogen-like high orbital angular momentum states can be directly accessed by laser light from low-lying atomic states due to the admixing of the core configurations with large angular momentum components. Second, alkaline earth and lanthanides offer the possibility to optically trap the atoms in Rydberg states via the valence core electrons and the corresponding ion core polarizability [57,[66][67][68].
This article is structured as follows. In section 2, we introduce the SO(4) formalism, the coupling to external fields, the form of dipole-dipole interactions and the role of quantum defects. In section 3, we discuss how to encode and manipulate quantum information in the n-manifolds of Rydberg atoms, and we show the corresponding generalized Heisenberg models for Rydberg atoms in tweezer arrays. In section 4, we demonstrate how our approach can be used to simulate continuous variable systems, and in particular we discuss the paradigmatic sine-Gordon model and its massive extension. In section 5, we discuss a protocol to perform state transfer and entangling gates between two atoms within the manifold of states near the circular level. We conclude our work in section 6 with an outlook on future perspectives.

SO(4)-symmetric Rydberg manifolds
In this work, we consider atoms trapped in tweezer arrays (e.g. alkali, earth-alkaline and lanthanide atoms) that are excited to hydrogen-like Rydberg states with principal quantum number n in the regime of linear Stark and Zeeman effect. The corresponding manifold of n 2 levels displays a clean and regular structure originating from an underlying SO(4) symmetry described by two large angular momentaĴ a andĴ b , as illustrated in figure 1, which provide an efficient description to treat the single-and many-body problem. In this section, we introduce the angular momentum structure of the Hydrogen problem, the resulting Rydberg state manipulation with electromagnetic fields and interactions, and discuss how these results translate to experimentally relevant multi-electron atoms.

Algebraic solution of the hydrogen problem
The SO(4) symmetry of the hydrogen atom [1,69,70] allows for an elegant algebraic description of the Rydberg manifold of states with a fixed principal quantum number n. In particular we consider the non-relativistic hydrogen Hamiltonian,Ĥ H =p 2 /(2m r ) − e 2 /(4πϵ 0 |r|), where we ignore electronic spin degrees of freedom and m r is the reduced electron mass, −e the electron charge, ε 0 the vacuum permittivity andp andr are the momentum and position operators, respectively.
A convenient approach to solve for the bound states ofĤ H is to exploit the conservation of orbital angular momentumL and Runge-Lenz (RL) vectorÂ (see appendix A for details). These two sets of conserved quantities generate the SO(4) symmetry of the hydrogen Hamiltonian. As the SO(4) group is doubly covered by SU(2)×SU(2), we can construct two commuting angular momentum operators that read and therefore obey [Ĵ a,i ,Ĵ a, j ] = iϵ i jkĴa,k and [Ĵ b,i ,Ĵ b, j ] = iϵ i jkĴb,k , furthermore, the angular momentum raising and lowering operators are defined asĴ a(b),± =Ĵ a(b),x ± iĴ a(b),y . Note, in this paper we use units where ℏ = 1. Due to the orthogonality between the orbital angular momentum and the RL vector (L ·Â = 0), the lengths of the two angular momenta are constrained,Ĵ 2 a =Ĵ 2 b . Within this formalism, the Hamiltonian can be recast in the simple formĤ H = −Ry/[2(Ĵ 2 a +Ĵ 2 b ) + 1], see appendix A, where Ry = e 2 /(8πϵ 0 a 0 ) is the Rydberg energy and a 0 = 4πϵ 0 /(e 2 m r ) the Bohr radius. A natural basis for the hydrogen atom eigenstates is therefore given by states |J, m a , m b ⟩, wherê with J = 0, 1/2, 1, . . . and m a (b) ∈ {−J, −J + 1, . . . , J}. The corresponding energies of the Hamiltonian take the form E n = −Ry/(2J + 1) 2 = −Ry/n 2 from which we can read off the relation between the principal quantum number n and angular momentum length J = (n − 1)/2. Therefore, a manifold of states with fixed n hosts n 2 degenerate states |J, m a , m b ⟩, which can be labeled by the quantum numbers of two angular momenta. A final result of the angular momentum formalism discussed in this subsection is the form of the electric dipole operatorμ = −er whose transition matrix elements, within a single n manifold, can be identified aŝ see [1,[71][72][73] and appendix A. In the following, we will exploit this relation to analyze the angular momentum level structure in the presence of external electric fields and to express dipole-dipole interactions in terms ofĴ a andĴ b .

Coupling to external fields
The coupling of the electron to weak static external electric F and magnetic B fields lifts the degeneracy of states with fixed n, and is described by the Hamiltonian where µ B = e/(2m r ) is the Bohr magneton. For a manifold of hydrogen-like states with a specific n, theL and µ can be replaced by angular momentum operators according to equations (1) and (3), respectively. The replacement ofμ is valid below the Ingris-Teller limit |F| < F IT = 2Ry/(3 ea 0 n 5 ), i.e. when couplings to adjacent n ′ = n ± 1 manifolds are weak with respect to the gap E n − E n±1 [3,74,75]. Analogously, diamagnetic couplings are negligible in the limit |B| < 2Ry/(µ B n 4 ) [76].
In the following analysis, we focus on a single manifold of states with fixed n and take the two fields to be parallel and oriented along the z-axis, F = Fe z and B = Be z , see figure 1(a), such that the projection m l = m a + m b of the orbital angular momentumL z remains a good quantum number. Within the limits discussed in the previous paragraph the Hamiltonian from equation (4) becomeŝ with ω a = ω S − ω Z and ω b = ω S + ω Z , where ω S = 3n ea 0 F/2 is the Stark splitting and ω Z = µ B B is the Zeeman splitting (Larmor frequency). The linear Stark and Zeeman effect lifts the energy levels degeneracy, as shown in figure 1(b) as a function of the electric field strength for 87 Rb. The resulting spectrum, at fixed electric field, takes a regular spacing that forms a rhombic-shaped level structure [74,77,78], as shown in figure 1(c). Deviations around m l = 0 of the regular spacing occur due to core corrections for non-hydrogenic atoms, which we discuss below in section 2.4. For positive values of the electric field of a few V cm −1 , and for magnetic fields of hundreds of Gauss, and principal quantum numbers in the range 30 < n < 70, we have ω Z < ω S such that ω a(b) > 0. Within this parameter regime the absolute values of ω a and ω b can be of the order of 2π × (10 2 − 10 3 ) MHz, which allows for direct coupling between the angular momentum states with microwave (MW) radiation.
Let us consider two circularly polarized MW fields whose electric field amplitudes are given by , electric field amplitudes F a (b) and polarizations e ± = ∓(e x ± i e y )/ √ 2. The electric field of the MWs couple to the dipole operatorμ, which for hydrogen-like states is represented by equation (3). By combining the effect of the static and time-dependent fields and after going to the rotating frame, equations (4) and (5) give rise to the most general Hamiltonian linear in angular momentum operatorŝ where we dropped fast oscillating terms and the rotating frame transformation is defined bŷ 2) are independently tunable detunings and complex Rabi frequencies, respectively. We note that such MW control has already been proposed [2] and demonstrated in recent experiments [79]. Deviations from perfectly circularly polarized MW fields can give rise to unwanted transitions, which can be suppressed by tuning the level spacing imbalance ω a − ω b via the magnetic field.

Dipole-dipole interactions
We now turn our attention to the role and form of interactions. Let us consider a pair of atoms i and j, with a relative position R i j = R i j e i j , here R ij is the inter-atomic separation and e i j = (cosϕ i j sinθ i j , sinϕ i j sinθ i j , cosθ i j ) T is the corresponding unit vector, where θ ij and ϕ ij are the azimuthal and polar angle, respectively, see figure 1(a). For distances of a few to tens of micrometers, i.e. relevant for optical tweezer experiments, the interaction between two atoms is dominated by the electric dipole-dipole coupling [80,81] whereμ (i) andμ ( j) are the dipole operators of atom i and j, respectively.
We consider the situation where the two atoms i and j are initially prepared in a Rydberg n-manifold, and dipole-dipole interactions are weak as compared to the energy gap to adjacent n ′ = n ± 1-manifolds, such that no transitions to other n ′ -manifolds occur. Within these criteria, the dipole operator can be expressed in terms of angular momentum operators through the application of equation (3), and the interaction Hamiltonian can be compactly written aŝ with V i j = (3n ea 0 ) 2 /(16πϵ 0 R 3 i j ). Corrections to the above Hamiltonian, like van der Waals interactions mediated by states with principal quantum number n ′ ̸ = n, are discussed in appendix C.
Similarly to other dipolar systems [82], the spatial anisotropy (θ ij and ϕ ij ) of interactions and the geometrical arrangement of the atoms can be exploited to tune the dipole-dipole interactions. Furthermore, time-dependent methods [83][84][85] can also be applied to control and design the interaction Hamiltonian, as recently experimentally demonstrated for magnetic dipolar atoms [86] and for Rydberg atoms [87,88]. This will therefore allow us to realize and tune a whole class of interaction Hamiltonians.

Nonhydrogenic atoms and the electron spin
The Rydberg manifolds of states with n ≫ 1 of experimentally relevant alkali, alkaline earth and lanthanide atoms follow, up to core corrections, the SO(4) symmetric hydrogenic results presented in section 2.1. In particular, only a few low orbital angular momentum states, which are typically l = s, p, d(, f ) states, penetrate the core region and are therefore energetically shifted away from the n manifold under consideration. Their energy shift is described by quantum defect theory [89,90] and is given by E n,l = −Ry/(n − δ l ) 2 , where δ l is the quantum defect [91].
Including the quantum defects, the Rydberg Hamiltonian from equation (5) for a given n manifold becomes [92] where the threshold value l * is dictated by the magnitude of the quantum defects, which is different for each atomic species [93]. For example, for the case of Rb, and sufficiently strong electric fields, l * = 3 while for Er typically l * = 4. In contrast to the hydrogenic case, the energies ofĤ QD n do not form a perfectly regular rhombus structure anymore because the inner region |m l | < l * presents an irregular spacing [3,94]. However, for n ≫ 1 the large majority of states follow the hydrogenic description, as shown in figure 1(c) for the n = 21 manifold.
The treatment of non-relativistic Hydrogen has thus far ignored the effect of spin-orbit coupling, which in fact should be taken into account for realistic atoms. Nevertheless, spin-orbit coupling primarily only affects the low l states which are shifted away from a Rydberg manifold n due to quantum defects. The spin orbit coupling of high l states instead, is negligible for sufficiently large magnetic fields, due to a strong-field spin Zeeman (or Paschen-Back) effect [95]. Hence, the rhombus manifold, which consists of only states with high l, decouples from the electron spin an thus comes in two copies, |J, m a , m b ⟩ ⊗ |s, m s ⟩ with the electron spin state |s, m s ⟩ = |1/2, ±1/2⟩, energetically resolvable through the spin Zeeman splitting. This provides in principle a further degree of freedom that doubles the available state space, but in the remainder of this paper we will focus on a single spin state |s, m s ⟩ = |1/2, 1/2⟩. Let us finally comment on possible decoherence mechanisms of the hydrogen-like Rydberg states of a given n-manifold. First, radiative decay due to spontaneous emission and black body radiation is a well known source of decoherence in Rydberg experiments [96][97][98]. However, similar to experiments with circular Rydberg states [61,99], black body radiation can be suppressed in a cryogenic environment, while at the same time spontaneous emission is strongly reduced for all states with m l > l * since the primary decay channels to energetically low-lying states are dipole-forbidden (for further details see appendix B.1). Radiative decay rates can therefore safely be assumed to be well below all other relevant energy scales in this work.
Further sources of decoherence can for instance originate from fluctuations in electric and magnetic fields, which directly affect the energy levels ω a(b) , and therefore enter as dephasing. Experimentally demonstrated electric field stabilizations on the level of tens of µV cm −1 [100] and magnetic field stabilization of tens of µG [101] should however be sufficient to render these decoherence effects negligible. Another source of decoherence for tweezer trapped Rydberg atoms is motional dephasing, due to the trapping potential being dependent on the internal state. Trapping of the core-electron, as is possible for alkaline-earth and lanthanide atoms, would remove this problem by providing state insensitive traps [66][67][68], see appendix B.2. Furthermore, atoms internally prepared in the Rydberg state experience (state dependent) dipole-dipole forces, which leads to heating of the initially laser-cooled [50,51] motional state. These heating rates are negligible for the inter-atom distances of few tens of µm and for the time scales considered in this manuscript, see appendix B.2.
In summary, Rydberg states with a fixed principal quantum number n of experimentally relevant atomic species provide a natural setting for quantum simulation and quantum information processing with a large internal state space. Except for very few states that are affected by quantum defects, the Rydberg n-manifold inherits the physics of the hydrogen states. This includes the single-particle couplings to external fields, see equation (6), and two-body dipole-dipole interactions, see equation (7), giving rise to generalized Heisenberg models for arrays of many atoms. In the next section, we consider different subsets of the n-manifold and illustrate several concrete examples of many-body Hamiltonians that can be naturally implemented with this platform.

Engineering of quantum many-body models
In this section, we consider different subsets of the n-manifold in order to engineer specific examples of many-body Hamiltonians for the hydrogen-like states, thus avoiding the quantum defect modified region. More specifically, we discuss in the following two manifolds that we employ for quantum simulation and quantum information processing applications, namely the triangular and edge manifold. In addition, we discuss a manifold that can be employed for state preparation and measurement, the vertical manifold.

Triangular manifold
The first and largest manifold that we consider is the triangular manifold defined for a fixed n by and m l ⩾ l * . This defines the set of all states unaffected by quantum defects on the right side of the rhombus and highlighted by the orange triangle in figure 1(c). In addition to dipole-dipole interactions and single-particle control, we will demonstrate how to design nonlinear terms in the angular momentum operators with ponderomotive manipulation techniques that can be used for different kinds of state manipulation and engineering schemes. Later in this work, these nonlinear terms will be instrumental for the analog simulation of QFTs. Furthermore, near the circular level we derive an effective description by using the Holstein-Primakoff transformation, which enables us to perform a state-transfer operation between pairs of atoms as a building block for qudit-based quantum computing. The quantum defect modified region at low m l can be avoided by choosing the system parameters such that the dynamics is restricted to the tip of the diamond, as discussed in section 3.4.

Edge manifold
The second manifold that we consider is a subset of the triangular manifold H t and is defined by the states with maximally-polarized angular momentumĴ b , namely and m a + m b ⩾ l * , which is highlighted in red in figure 1(c). This manifold is closed under dipole-dipole interactions governed by H dd as all processes leaving H a are energetically suppressed by the Stark-and Zeeman-splitting, as discussed in section 3.5. In addition to the control offered by the triangular manifold, here we can also engineer additional nonlinear squeezing terms by off-resonantly MW coupling between different n ′ > n manifolds. This additional ingredient will provide the kinetic term required for the simulation of QFTs.

Vertical manifold
The last manifold considered in our discussion is the vertical manifold defined by and fixed m l = l * , which defines a set of states immediately next to the quantum defect modified region highlighted in green in figure 1(c). Although this set of states is not closed under dipole-dipole interactions, its unique properties can be exploited for state preparation and readout for the H a manifold.

Triangular manifold H t
The triangular manifold H t is defined in equation (9) and is highlighted in orange in figure 1(c). Here below we discuss a specific example of a many-body Hamiltonian, the ponderomotive manipulation technique to engineer a class of nonlinear terms and we provide an effective theory for the states near the circular level.

Model Hamiltonian
To be specific, we consider the case where all atoms are placed in a plane perpendicular to the static electric and magnetic fields, b,+ in equation (7) are energetically suppressed by the Stark-and Zeeman-splitting.
By considering the single-particle and many-body terms introduced in the previous section in equations (6) and (7), with the addition of the ponderomotive drive that we detail in the next subsection, the many-body Hamiltonian takes the form where we transformed to the rotating frame defined below equation (6) under the resonant condition P is the ponderomotive coupling for the ith atom defined aŝ where λ κa,κ b denotes complex coupling strengths. The ponderomotive coupling in equation (13) corresponds to a transfer of orbital angular momentum by κ = κ a + κ b , where κ a(b) can take on negative and positive values [thus with the identification (Ĵ a(b), For the special case κ b = 0 implying κ = κ a , the ponderomotive HamiltonianĤ P reduces to λ κ (Ĵ a,+ ) κ + λ * κ (Ĵ a,− ) κ , where we use the simplified notation λ κ ≡ λ κ,0 . In this case, this coupling Hamiltonian can be used, for instance, to generate a two-axis twisting term [102], i.e. i(Ĵ a,+ ) 2 − i(Ĵ a,− ) 2 for κ = 2.

Ponderomotive manipulation
In order to realize the nonlinear coupling in equation (13), we employ the ponderomotive manipulation techniques [62,[103][104][105]. The core idea is to interfere two co-propagating (along the z-direction) hollow Laguerre-Gauss laser beams with different frequencies ω 1 and ω 2 and non-zero orbital angular momentum. The corresponding interference pattern forms a ponderomotive potential U κ (r)e −iκϕ e −iδωt + H.c. , rotating at the beating frequency δω = ω 1 − ω 2 , which can be arranged to be in the MW frequency domain to match the resonance condition δω = −κ a ω a + κ b ω b . When centered at the atomic position, this spatially-rotating ponderomotive potential transfers κ = κ a + κ b quanta of orbital angular momentum, see figure 2(a).
For simplicity, we focus here on the case κ b = 0 and therefore κ = κ a . The transition matrix elements of the ponderomotive coupling U κ,ma,m b ≡ ⟨ψ t ma+κ,m b |U κ (r)e −iκϕ |ψ t ma,m b ⟩ are shown in figure 2(b) for κ = 2, 4, 6, m b = J and different values of n. As discussed in appendix D, when the beam waist is large compared to the Rydberg orbital radius, the leading order expression of these matrix elements is For simplicity, we extract the coefficient λ κ by a least square fit and show the accuracy of this procedure in figure 2(b). The degree of accuracy can be controlled by increasing the laser beam waist, which requires more laser power to maintain the same strength for the transition amplitudes. The parameter λ κ can be complex and is controlled by the relative phase of the Laguerre-Gauss beams. Furthermore, the coupling matrix elements U κ,ma scale as n 2κ for |m a | ≪ J, and are in the range 2π × (10 2 − 10 3 ) kHz for hundreds of mW of laser power. Further details of these calculations and precise beam parameters are discussed in appendix D.
We point out that for fixed laser power and n achievable coupling strength decrease as |κ a | + |κ b | is increased. For |κ a | + |κ b | ≲ 6 and a combined laser power of hundred mW per site give rise to coupling strength of hundreds of 2π×kHz, see figure 2(b). Moreover, notice that the expression (13) also contains the linear expression (6) obtained via a MW driving as a special case. However, as the ponderomotive coupling is driven by laser light, it directly provides a way to locally engineer these couplings, i.e. at the level of a single Rydberg atom, which can be exploited for digital quantum simulation purposes.

Effective theory near the circular state
A model Hamiltonian similar to equation (12) can also be realized with coupled atomic ensembles [106], where a single large angular momentum of length J is constructed out of 2 J two-level atoms. In our case two large angular momenta, both of length J, are encoded in a single atom with principal quantum number n = 2J + 1, which can be manipulated with external fields enabling, for example, the preparation of squeezed states, see equation (13). While atomic ensembles are usually coupled via atom-light interactions, here pairs of angular momenta are naturally coupled by dipole-dipole forces.
Analogous to atomic ensembles, where states close to the maximally stretched state are described by continuous variables [106], we now focus on the Hamiltonian (12) on a subspace of states close to the maximally stretched (circular) level |ψ t J,J ⟩. We thus perform a Holstein-Primakoff transformation J In this framework, the circular level serves as the vacuum for both bosonic modes, |0, 0⟩ ≡ |ψ t J,J ⟩, and the Hamiltonian (without MW and ponderomotive couplings) takes the form of a generalized spinful Bose-Hubbard model pair gain/loss processes, and U i j = 2V i j density-density interactions. In contrast to standard Bose-Hubbard models the pair gain/loss term w ij violates total particle number conservation and can be used to generate parametric squeezing. Further below in section 5 we outline another application of equation (14), where we employ the hopping term h ij for quantum state transfer and an entangling gate between a pair of atoms. We note that this Hamiltonian is valid inside the triangular manifold and not in the quantum defect modified region, see figure 1(c). In the limit |w i j | ≲ −∆ σ the low energy dynamics (with respect to −J∆ σ ) is constrained to the tip of the triangle and this condition is met. Furthermore, interaction terms causing transitions into the quantum defect modified region are typically energetically suppressed due to the irregular level spacing in this region. The physics of equation (14) can be explored starting from, for example, the vacuum state |ψ t J,J ⟩, which can be initialized through well established circularization methods, i.e. via the adiabatic rapid passage method [2,[107][108][109], the crossed electric and mangnetic fields method [3,110], multi MW photon transfer [79,111] or direct ponderomotive coupling [105]. Another possible scheme to prepare more general product states is outlined in section 3.6, where we describe a generalization of the adiabatic rapid passage protocol used to prepare circular levels. Readout can be performed by applying energy-selective MW pulses that transfer the population of individual levels to different n ′ -manifolds, which are subsequently resolved by an electric field ionization measurement [79].

Edge manifold H a
The edge manifold H a is defined in equation (10) by havingĴ b maximally polarized, as highlighted in figure 1(c). Notice that we could define an analogous manifold H b by considering the states witĥ J a maximally polarized. Choosing the static external electric and magnetic fields such that [112]. Since H a is a subset of H t , it inherits many of the properties discussed before, such as the ponderomotive control. Additionally, we can also construct a nonlinear squeezing term where χ is the squeezing strength. This can be generated by off-resonantly coupling the H a manifold with principal quantum number n to another H a manifold with n ′ > n, through MW radiation.
In the blue detuned regime, see figure 2(c), and z-polarized MW fields, the system can be treated as a collection of independent two-level systems, one for each m a . As shown in figure 2(c), the detuning ∆ ma = ∆ 0 + δω a m a of each two-level system varies with m a , where ∆ 0 is an overall detuning and δω a = ω a (n ′ ) − ω a (n) is the differential level shift. Moreover, the coupling strength Ω ma depends smoothly on m a through the dipole transition matrix elements. In the regime ∆ ma ≫ Ω ma the coupling leads to an AC Stark shift Ω 2 ma /(4∆ ma ) that can be expanded in a power series of m a when ∆ 0 ≫ δω a m a , namely where the constant and linear terms are absorbed in ∆ 0 and δω a , respectively. The (quadratic) leading term of this expansion therefore gives rise to the anticipated squeezing term, as shown in the upper panel of figure 2(d). Cubic terms can be canceled exactly by carefully picking ∆ 0 and Ω ma , as discussed in appendix E.
For moderate values of the MW intensity, χ can be on the order of 2π × 100 kHz, thus giving rise to level shifts as large as 2π × 10 MHz for |m a | ∼ J. For such values of the level shifts, a small admixture of states from the n ′ manifold is present, see the lower panel of figure 2(d). Moreover, imperfect MW polarization can lead to unwanted two-photon Raman processes, which can be made off-resonant by tuning ω a ≫ Ω 2 ma /∆ ma . By taking into account all the terms discussed so far, the Hamiltonian for the H a manifold takes the form The first line describes linear and nonlinear single-particle terms. As previously mentioned, the κ = 1 term can also be obtained by MW engineering (λ 1 = Ω a ), see equation (6). The second line describes long-range dipolar interactions under the assumption that the atoms are placed in a plane perpendicular to the static electric and magnetic fields (θ i j = π/2), thus realizing a 'large-spin' XXZ model. In section 4, we will discuss two examples of specific many-body models obtained from the Hamiltonian (16) that provide a direct application of the platform developed in this work. State preparation and readout for the H a manifold require specific schemes involving the states in H v , which we discuss in detail in the next subsection.

Vertical manifold H v
The vertical manifold H v defined in equation (11) and highlighted in green in figure 1(c) plays a fundamental role in our system, for multiple reasons. First, this manifold can be accessed from atomic ground states using a single or multiple laser transitions (depending on the atomic species), thus enabling to coherently (de)populate H v . Second, a mapping scheme that we demonstrate below allows to transfer arbitrary states from H v to H a and vice versa with high fidelity. As a consequence, the manifold H v provides the capabilities to perform state preparation and readout on H a , which make the Hamiltonian in equation (16) a particularly accessible many-body model.
The presence of low orbital angular momentum components in |ψ v ma ⟩, see figure 3(a), offers a unique opportunity for laser coupling the atomic ground state to the Rydberg manifold H v . In particular, for lanthanides like erbium the valence electrons have f -orbital character in their electronic ground state due to a submerged shell structure [65]. The H v manifold of erbium in turn contains an admixture of Rydberg g-states, which allows direct laser coupling from the ground state |gs⟩ to specific states in H v , selected by the The inset shows the instantaneous eigenenergies of the adiabatic rapid passage path in the rotating frame defined above equation (6) for a single ma = 0 and for parameters n = 31, laser frequency, see figure 3(b). For multiple laser frequencies, this coupling is described in the rotating frame by the HamiltonianĤ where Ω ma,gs are the individual Rabi frequencies that satisfy |Ω ma,gs | ≪ |ω S |. This laser adressability can also be used to perform projective readout within the H v manifold employing quantum gas microscope techniques [113,114]. Note that the connectivity offered by the individual Ω ma,gs in equation (17) provides a convenient setting to perform holonomic quantum computing operations within the H v manifold [115]. We now proceed to discuss a scheme realizing the state transfer via a adiabatic rapid passage that generalizes the one originally developed to prepare circular Rydberg levels [2]. We start by noting that a state |ψ v ma ⟩ is connected to the final state |ψ a ma ⟩ for any m a by a multi-photon MW transition (see the highlighted red arrow and levels in figure 3(b)), while transitions into the defect region are energetically suppressed. The required time-dependent drive is based on the MW HamiltonianĤ MW (t) (equation (6)) with Ω a = 0, such that only couplings among the highlighted states in figure 3(b) are present. During the driving sweep, the detuning ∆ b (t) has to change sign and the Rabi coupling Ω b (t) has to be turned on and off again in order to adiabatically (|Ω b (t)| ≪ |ω a(b) |) connect the different states |ψ v ma ⟩ and |ψ a ma ⟩. Note that in general the individual states acquire a state-dependent dynamical phase that can be canceled by an appropriate transformation in H v .
In the following, we numerically demonstrate the adiabatic rapid passageÛ va for a selected state (m a = 0) and a particular driving protocol ∆ b (t) = ∆ 0 cos(π t/T) and Ω b (t) = Ω 0 sin(π t/T). In the inset of figure 3(b), we show the instantaneous energy levels in the rotating frame defined above equation (6) as a function of time. For t = 0, several energy levels are present. The ones with positive energy (E > 0) are equally spaced and correspond to the levels highlighted in red in figure 3(b). Levels with different m a are not shown, as they are not coupled by the MW field. As the adiabatic sweep progresses, the initial state |ψ v ma ⟩ changes into the target final state |ψ a ma ⟩, as indicated by the red line in the inset. However, one can notice that several level crossings take place during the time evolution. These occur with states that have negative energy E < 0 at t = 0, corresponding to states with m l < −l * and the same value of m a , i.e. located on the left side of the rhombus level structure. Additionally, level crossings with states from the defect modified region (|m l | < l * ) can also be present. All these crossings are in fact avoided crossings with a very narrow gap that is determined by multiple off-resonant transitions through the defect region. Thus, they can be diabatically crossed via Landau-Zener tunneling without hindering the sweep protocol.
Due to the highly regular Rydberg level spacing, we remark that the passage connects all the states in H v to their respective counterpart in H a . For ∆ 0 and Ω 0 on the order of 2π × 10 MHz and T on the order of a few µs, the passage reaches numerically calculated fidelities above 99% for all the states.
In summary, the hydrogen-like Rydberg manifold discussed in this paper can be accessed by laser exciting atoms to Rydberg states with small orbital angular momentum components and a subsequent transfer to states with large orbital angular momentum components. This allows to explore the physics offered by the engineered many-body models introduced in sections 3.4 and 3.5. In the following Sections we show how these models can be used to simulate QFTs and how controlled entanglement between pairs of atoms can be generated.

Application: quantum simulation
The central feature of the simulator discussed in this work is the high-dimensional and regularly structured Rydberg manifold leading to a many-body problem described by Heisenberg models with large spins (or angular momenta). A promising application of these models is the simulation of QFTs [4], which we discuss in this section. For instance, a scalar quantum field is represented by pairs of conjugate variables on each spatial point, which span an infinite-dimensional space. These can be conveniently represented by the large spins of the Rydberg n-manifold, a procedure that naturally requires a truncation set by the finite (large) spin length.
Our construction provides a bottom-up approach to simulating QFTs, differently from other systems, e.g. ultracold atomic gases [116], where QFTs appear as a low-energy description. The Rydberg platform offers the opportunity to avoid unwanted thermal effects inherent to atomic gases, thus allowing us to access the zero temperature limit of these theories and to freely choose different geometries or dimensionalities. Moreover, the programmability and control of the platform allows for the preparation of initial states of interest, e.g. product states or other states obtained through adiabatic sweeps, the readout of the corresponding out-of-equilibrium quantum dynamics via quench protocols and the exploration of phase diagrams from the weakly-to the strongly-interacting regime. In particular, our approach paves the way for simulating non-equilibrium dynamics of scalar quantum fields in real-time, something which is notoriously hard for classical simulators and in general also not accessible with implementations based on low-energy approximations.
In this section, we illustrate these possibilities with two case studies, namely the massless and massive sine-Gordon (SG) model. The first case is a paradigmatic and ubiquitous integrable scalar QFT that plays an important role in many areas, ranging from high-energy physics [117] to condensed matter [118]. The second case is a scalar QFT that is dual to 1+1D quantum electrodynamics (QED), namely the massive Schwinger [119]. The quantum simulation of such a theory requires a continuous variable representation of the scalar fields, which are naturally embedded in the Rydberg manifold, as we show below.
In section 4.1, we first discuss how the platform's unique properties, namely the availability of spins with large spin length J ≫ 1, can be used to faithfully approximate the physics of continuous variables. In the following section 4.2, we show that the naturally occurring interactions enable us to realize a lattice regularization of the SG model. For this example, we focus on close-to equilibrium properties in order to provide a clean benchmark of our setup. Finally in section 4.3, we turn our attention to gauge field theories, more precisely QED in one spatial dimension, the so-called massive Schwinger model. We discuss how to realize a dual lattice formulation of this model using the Rydberg architecture. Here, we demonstrate how to probe nontrivial gauge theory dynamics, relevant for investigating processes like pair production and string breaking. These demonstrations serve as a qualitative illustration for potential applications in non-equilibrium situations that are inaccessible by other means.

Conjugated continuous variables
Standard quantum-mechanical positionφ and momentumπ operators obey the canonical commutation relation [φ,π] = i. For quantum mechanics on a unit circle, it will be more convenient to consider phase operators e iφ satisfying π, e iφ = e iφ , where the momentum operatorπ takes discrete eigenvalues m π ∈ Z on its eigenstates |m π ⟩. The commutation relation (19) shows that e iφ , similar to a spin operatorĴ + , acts as a raising operator on |m π ⟩. Taking into account the normalization of the operators, this suggests to identify The (2J + 1)-dimensional Hilbert space spanned by eigenstates |m⟩ ofĴ z with m = −J, . . . J thus 'regularizes' the infinite-dimensional Hilbert space of the continuous variables by cutting off the spectrum ofπ at a maximum value |m π | ⩽ J [120]. The identification becomes exact for states around m = 0 in the limit J ≫ 1. We now numerically demonstrate how the above identification converges to a continuous variable representation in the large spin limit J ≫ 1. Here, we focus on a single Rydberg atom in the subset of states H a described by the Hamiltonian in equation (16). Choosing ∆ a = 0 and defining Ω ′ = −2Ω a J(J + 1) and λ ′ = −2λ κ e i θ [J(J + 1)] κ/2 for a fixed κ > 1, λ ′ > 0 and a phase θ, the single Rydberg atom Hamiltonian in the continuous variable limit becomeŝ The Hamiltonian in equation (21) describes a particle of mass (2χ) −1 moving on a ring and subject to a periodic potential, as depicted in figure 4(a). The spectrum of the continuum theory (21) can be numerically computed by exploiting the periodicity ofφ and solving the eigenvalue problem in Fourier space, thus yielding the result shown in figure 4(b) for λ ′ = 0. Depending on the ratio Ω ′ /χ, the low-energy physics of this model ranges from the free particle regime (Ω ′ ≪ χ) to the harmonic oscillator regime (Ω ′ ≫ χ).
We now consider an intermediate regime, where the nonlinear character of cosφ is evident, and calculate the spectrum of the spin model for increasing values of the angular momentum J, as shown in figure 4(c). For the chosen parameters and the selected energy range, the truncation accurately reproduces both bound and a few unbound states of the model with continuous variables, equation (21), when J ≳ 30. In figure 4(d), we further test the effect of λ ′ on the spectrum for a fixed value of the spin magnitude J = 50, and find that the spin model also captures the target model including the higher-harmonic potential, controlled by λ ′ .
We have thus verified that for experimentally relevant spin lengths J (or principal quantum numbers n = 2J + 1) the large spin captures the physics of a continuous variables theory in a single Rydberg atom. This capability constitutes the basis to engineer a QFT simulator using an array of interacting Rydberg atoms, which we illustrate in the following subsection for the SG model.

Sine-Gordon model 4.2.1. Overview of the model and Rydberg implementation
The SG model is a paradigmatic QFT that is described by the Hamiltonian [117,121] where the fields obey [φ(x),π(x ′ )] = iδ(x − x ′ ), M 0 has units of energy, distances are measured in units of inverse energy and β is a dimensionless parameter. Let us briefly summarize some key features of the SG model (see, e.g. [122]) and references therein). At β 2 = 8π, the theory exhibits a Berezinskii-Kosterlitz-Thouless (BKT) transition, separating a gapless phase at large β from a gapped phase at small β. In the gapped phase (0 < β 2 < 8π), the low-energy degrees of freedoms are fermions interacting attractively (repulsively) for β 2 < 4π (β 2 > 4π), which can be identified with quantized solitons of mass M β→0 −→ 8M 0 /β. In the attractive regime, the solitons form bound states, so-called breathers. The lightest one of them has a mass m 1 β→0 −→ M 0 (see appendix F for general expressions of the masses). For the implementation of the SG model we consider a one-dimensional array of tweezer-trapped Rydberg atoms excited to the manifold H a . The system is therefore described by the many-body Hamiltonian (16), which in the large J limit allows for a continuous variable description through equation (20) that readsĤ where we rescaled the interaction V ′ i j = V i j J(J + 1), while keeping ∆ a = 0, and λ ′ = −2λ κ [J(J + 1)] κ/2 for a single value of κ, which corresponds to either a MW (κ = 1) or ponderomotive (κ ⩾ 1) coupling, respectively. Notice that in the large J limit, the Ising terms have been neglected, which we benchmark in the next subsection.
In order to recover the continuum QFT in equation (22), we take only nearest-neighbor terms V ′ i,i+1 ≡ V ′ nn and rescale the continuous variables and Hamiltonian appropriately (see appendix G). The correspondence is achieved by the identification of parameters as follows where we introduced a short-distance scale ℓ that sets a UV-cutoff Λ ∝ 1/ℓ for the lattice regularization, ℓM 0 = κ 2λ ′ /V ′ nn → 0. With these identifications, all the regimes of the theory, 0 < β 2 = 2κ 2 χ/V ′ nn < ∞, are in principle accessible. In particular, in appendix G we present realistic experimental parameter estimates, which indicate that values β 2 ≳ 8π can be reached with a spin length of J = 20. This suggests that this platform can be employed to investigate the critical properties of the model across the quantum phase transition.
Before turning to a benchmark analysis of our SG model implementation, let us briefly point out some existing proposals and experimental realizations. The SG model can be realized with ultracold atoms, for example, by using binary mixtures [123,124] or tunnel-coupled superfluids [125], and the theory has been tested in the semiclassical regime by measuring high-order correlation functions [116]. Recently, it has been proposed to realize a lattice regularization of the quantum SG model using superconducting circuits [126], which is most closely related to our approach with large spins discussed here.

Numerical benchmark in the massive phase
We now illustrate how the lattice SG physics emerges from the large spin approximation. In particular, we test the influence of finite spin length, the presence of Ising terms and discuss the effect of long-range tails. For this, we numerically analyze an array of N = 5 Rydberg atoms assuming periodic boundary conditions [127] and compare with analytical results for the low-energy excitation spectrum and ground state expectation values. Further below, we also briefly discuss how to measure the many-body gap through a realistic quench protocol.
We perform our analysis deep in the gapped phase, setting V ′ nn = 10χ and κ = 1, i.e. β 2 = √ 0.4 ≪ 8π, and vary λ ′ /χ, thereby effectively scanning different values of ℓM 0 . In this parameter regime, we expect the model to be well approximated by a quadratic (free) theory obtained by replacing cosφ i ≈ 1 −φ 2 i /2, and similarly for cos (φ i+1 −φ i ). After diagonalizing the quadratic theory (see, e.g. [128] and appendix H), we obtain the single-particle dispersion relation (c) Energy gap for a wide range of λ ′ /χ obtained via the quench protocol with α = π/50 described in the text (cross markers) and compared with the ED gap (empty circles). The plot shows also the free theory prediction (dashed line) and the exact lowest breather mass mB = m1 (dotted line). In the inset, the oscillation dynamics of ⟨Ĵy⟩ from which we extract the gap for λ ′ = 10χ. (d) Vertex operator expectation value ⟨e iφ ⟩ ↔ ⟨Ĵx⟩/ J(J + 1) obtained with ED (full circles) and within the free theory (dashed line). The discrepancy is mainly originating from the small spin length, as shown in the inset for λ ′ = 50χ. ED in panels (b)-(d) is performed for J = 10 and the site indices are dropped as we considered a translation invariant system. truncated the interactions to nearest neighbors. These low-energy massive excitations with dispersion ω q correspond to coherent displacements of the phase variable from the equilibrium position ⟨φ i ⟩ = 0, which can be interpreted as gapped spin waves of the original spin model.
As illustrated in figure 5(a) for finite J, the numerically calculated dispersion obtained by neglecting Ising terms and long-range contributions is in perfect agreement with the analytical formula for ω q . We find that including the Ising terms has the effect of reducing the overall bandwidth of the dispersion relation. In figure 5(a), we also show the effect of the long-range dipolar tail, which is calculated semi-analytically in appendix H by using Ewald's resummation [129]. For q ≈ 0, we find ω q ≈ 2χ[λ ′ − V ′ dd (c 2 q 2 + c ′ 2 q 2 logq)], with c 2 ≈ −0.739 and c ′ 2 ≈ 1/2. The long-range dipolar interactions have therefore the effect of increasing the bandwidth of the excitation spectrum, which gains a non-analytical sub-leading contribution at small momenta. Note that the energy gap ∆E = ω 0 = √ 2χλ ′ is thus unaffected by the long-range tails. We now study the gap ∆E of the model in more detail. As previously stated, we expect the Ising terms' contribution to vanish in the J ≫ 1 limit. This is confirmed in figure 5(b), where we show the dependence of the energy gap ∆E in the finite-spin model as a function of J with J = 5, . . . , 12, both with and without Ising terms. Comparing against the predicted result ω 0 = √ 2χλ ′ of the free theory, we find that both finite-spin results overestimate the value of the gap. These systematic errors vanish smoothly as J is increased, demonstrating that the finite spin model reproduces the expected gap in the large J limit. Given the regular J-dependence, we note that a finite spin-length scaling analysis can also be performed in an experimental realization by repeating appropriate measurements in different n manifolds. In principle, this allows to extrapolate the asymptotic value of the gap, thereby improving the prediction of the Rydberg quantum simulator.
In a quantum simulation experiment based on the Rydberg platform discussed in this work, the gap can be accessed by the following quench protocol. First, we prepare the ground state of equation (23) with cosφ i → cos(φ i + α), which can be obtained in the spin model by taking, for example, a MW field in equation (6) with complex Rabi frequency Ω a → Ω a e iα . For small α, this corresponds to a small uniform displacement of the spin expectation value away from ⟨Ĵ (i) y ⟩ = 0. Quenching to α = 0, the energy gap can be extracted by measuring the oscillation frequency of ⟨Ĵ The leading contribution of the induced dynamics within linear-response, i.e. for |α| ≪ π, is ⟨Ĵ 3 i ⟩ , which therefore gives direct access to the coherent displacement of the phase variable. The determination of the many-body gap can be accurately determined for a wide range of λ ′ /χ values, as shown in figure 5(c). The initial state preparation can be performed using the scheme outlined in section 3.6 for large λ ′ followed by an adiabatic ramp to the final value of λ ′ /V ′ nn and λ ′ /χ, which is protected by the many-body gap. The same scheme in section 3.6 also allows to measure the value of ⟨Ĵ (i) y ⟩ after the time evolution.
Several other types of correlation functions are also measurable via the scheme detailed in section 3.6. We give another example in figure 5(d), where we compute the ground-state expectation value of the vertex operator via ⟨Ĵ (i) + ⟩ ∝ ⟨e iφ i ⟩ = ⟨cosφ i ⟩ and compare it against the free theory prediction ωq . The plot displays a systematic difference between the numerical ED result and the free theory prediction. As shown in the inset for λ ′ = 50χ, this is mainly originating from the finite spin length (J = 10) used in the ED simulations. Note that the observables discussed so far can be measured by only using global MW control. However, more general correlation functions as ⟨Ĵ (i) +Ĵ ( j) − ⟩, with i ̸ = j, require single site control, which can be achieved by using electric field spatial gradients for MW fields [130] or local laser addressing, see section 3.4.2.
The main discussion of this section has focused on benchmarking low-energy SG properties in the massive phase. However, the Rydberg simulator offers the unique opportunity to prepare more general high-fidelity pure states and study their non-equilibrium quantum dynamics, e.g. after quenching the model parameters. An intriguing property of the SG model is the presence of topological excitations, or solitons, which are relevant to characterize the BKT transition. A possible application is the simulation of such high-energy metastable states [131][132][133]. This will provide a testbed to study questions related to false-vacuum decay [134], which plays a central role in cosmological models as well as in first-order phase transitions [135].

Massive Schwinger model 4.3.1. Overview of the model and Rydberg implementation
We now come to our final quantum simulation application, namely QED in one spatial dimension, also known as the massive Schwinger model. Here we make use of the fact that a non-integrable modification of the SG model, the so-called massive SG model, described by the Hamiltonian provides a dual low-energy effective description of the massive Schwinger model [119,121,136], described byĤ QED =ˆdx Here, γ 0(1) are the gamma matrices in two space-time dimensions, we abbreviatedD The duality requires to fix the mass M, the coupling u and the parameter β as where e denotes the electric charge, m is the free fermion mass, and γ ≈ 0.577 is the Euler-Mascheroni constant. The above duality is valid in the continuum limit for a UV cutoff (imposed on the Hamiltonian (25)) Λ ≫ e, m. In this limit, one obtains a QFT described by two dimensionless parameters, the coupling 0 ⩽ e/m ⩽ ∞ and the topological angle θ ∈ [−π, π].
The Schwinger model is a well-studied toy model that shares several qualitative properties with more complicated gauge theories such as quantum chromodynamics, for example confinement. Let us briefly discuss some key features of the model [119], starting with the case of θ = 0. For e/m = 0, the gauge fields trivially decouple and we obtain a theory of non-interacting massive fermions. For finite e/m > 0, the gauge fields mediate a density-density interaction ∝ρ(x)ρ(x ′ ) among the charge densitiesρ = eψ †ψ with a Coulomb potential growing with distance as ∝ |x − x ′ |. Fermions of opposite charge thus interact attractively and are confined into bosonic bound states. In the dual theory these are directly described by the fieldφ, which can be identified with the electric fieldÊ/e. In the limit e/m = ∞, namely u/M 2 = 0, these bosons propagate freely with a mass given by M.
A finite value of θ corresponds to a background electric field, which is identical to the field created by two static background charges ±eθ/π via the Gauss law constraint ∂ xÊ =ρ. The physics is periodic in θ with period 2π (which is manifest inĤ mSG ) because any background charge Q = ±ne with integer n is screened by the production of particle pairs induced by the corresponding background electric field, ⟨Ê⟩ = Q. For the special case of θ = ±π (without additional external charges), the ground state of the model undergoes a second-order quantum phase transition at (m/e) c ≈ 0.3335(2) [137], with order parameter ⟨Ê⟩ (or ⟨φ⟩), which falls in the Ising universality class and corresponds to the spontaneous breaking of theÊ → −Ê (or φ → −φ) Z 2 symmetry. For fixed m/e > (m/e) c , the system undergoes a first-order phase transition at θ = ±π, see figure 6(a) for a sketch of the phase diagram. Further details can be found in [119,136].
We obtain a lattice version of equation (25) by repeating a similar discretization as described in the previous section for the SG model. This requires to include the MW field and the ponderomotive couplings simultaneously. In particular, starting from the many-body Hamiltonian in equation (16), in the large J limit the continuous variable description readŝ where we use the identification , and neglect the Ising interactions. Our target model is again approached in the continuum limit with the identifications and keeping only nearest-neighbor terms V ′ nn ≡ V ′ i,i+1 . As discussed in appendix I, the massive SG model in equation (25) is recovered under the conditions We emphasize that a suitable choice of parameters in principle allows us to probe the whole non-trivial range of the theory with 0 < e/m < ∞. A specific realistic example of parameters is provided in appendix I, where we employ a spin length of J = 30, yielding e/m ≈ 4.5.
In recent years, there has been a growing interest in quantum simulating gauge theories. Experimentally, real-time dynamics of the massive Schwinger model has been realized digitally on trapped-ion [138] and super-conducting qubit [139] quantum computers. Analog simulations of directly related models have been performed with Rydberg arrays [140], cold mixtures [141] and Bose-Hubbard systems [142]. Our proposal differs from all of these existing realizations since we directly target the dual scalar field theory, which can be naturally implemented with Rydberg atoms in the H a manifold.

Numerical illustration: capacitor discharge
In the remainder of this section, we illustrate how our platform allows us to probe paradigmatic gauge-theory phenomena. Let us consider a system of size D with a given background electric field determined by θ ̸ = 0. This effectively forms a charged capacitor with plates separated by distance D. In general, the ground state of this system for sufficiently small e/m exhibits a nonvanishing electric field ⟨Ê(x)⟩ ̸ = 0, which is inhomogeneous due to the finite system size. Turning off the bias field, θ → 0, which corresponds to a quantum quench as indicated in figure 6(a), induces a discharge of the capacitor as dynamical charges are produced and accelerated across the capacitor [143], thus reducing the background electric field, as sketched in figure 6(b) (see [144] for another motivation to study quenches of the θ angle). The dynamical processes are reminiscent of Schwinger pair production and string breaking in strong fields, which leads to plasma oscillations due to the coupling between the gauge field and the fermionic matter [145]. In the dual picture, these processes can be probed by monitoring the damped oscillations of ⟨Ê⟩ = −e⟨φ⟩/(2π) together with the charge density ⟨ρ⟩ = −e⟨∂ xφ ⟩/(2π). In general, studying these dynamical processes is computationally hard, which motivates us to address this problem with an analog quantum simulator.
While an experimental realization of our proposed setup will be able to faithfully simulate the full dynamics in the massive SG model, we necessarily restrict ourselves to a relatively small lattice system with small spin lengths to illustrate the qualitative behavior in a numerical simulation. To be specific, we consider a lattice with N = 7 atoms, spin length J = 4 and κ = 4. The result of a quench θ = π/2 → 0 is shown in figures 6(c) and (d), where the chosen parameters are indicated in the caption. ⟩, which correspond to the lattice version of the electric field ⟨Ê i ⟩ and the charge density ⟨ρ i ⟩, respectively, up to appropriate proportionality constants. Even for the small system size and finite spin length considered here, we observe the expected damped oscillations of the electric field and the corresponding charge dynamics. These are driven by the interplay of pair creation and the subsequent interaction among the fermions with the gauge fields. At later times, we find partial revivals that we attribute to finite-size effects, thus prohibiting us from studying the long-time limit in our simulations. Nevertheless, as shown in figure 6(d), the plasma oscillations show an inhomogeneous damping, consistent with the expectation that the electric field decays in the middle of the capacitor [143], while a finite charge density accumulates at the system's boundaries, leading to an effective screening effect, reminiscent of string breaking. In contrast to our benchmark simulations, an experimental realization can be scaled up to large spin lengths J ≫ 1 and large lattice sizes, namely N ≫ 1 Rydberg atoms, with the perspective to reach the QFT limit. Such a realization will enable a quantitative analysis of the intricate behavior towards later times, which can help to improve our general understanding of the thermalization process in gauge theories [146,147].

Quantum information prospects
The large and regularly structured Rydberg manifold together with the high level of external control and manipulation available for such a system offer an opportunity to exploit this platform for qudit-based quantum information processing. As an illustration we present an entangling gateÛ E i j between a pair of atoms i and j, where the quantum information of atom i (j) is stored in its H a (H b defined in the beginning of section 3.5) manifold as depicted in figure 7(a). A state transfer operationÛ ST i j moves the quantum information from the b-degree of freedom of atom j to the b-degree of freedom of atom i, where a general single particle operationÛ SP i acting on atom i subsequently entangles the two subspaces, before transferring the b degree of freedom back to atom j.
We now discuss the state transfer gateÛ ST i j that is inspired by the one developed for coupled harmonic oscillators, see for example [148]. The state transfer protocol can be realized near the circular level |0, 0⟩ by exploiting the dipole-dipole interactions, see section 3.4.3, as we show at the end of this section. To this end, consider the initial two-atom state |ψ 0 ⟩ = |n i , 0⟩ i ⊗ |0, n j ⟩ j , where the atom i encodes information in the H a manifold and the atom j encodes information in the H b manifold. In order to transfer the information from atom j to atom i, we consider the , see equation (14), valid for n i , n j ≪ J, which achieves the desired state transfer, see appendix J for further details. Generalization to arbitrary initial superposition states is straightforward. Entanglement between the two atoms can then be generated by constructing a more general entangling gateÛ E i j defined by the sequenceÛ see figure 7(a), where the single-particle unitaryÛ SP i required for the entanglement generation between the two atoms can be realized, for example, by the ponderomotive manipulation techniques discussed in section 3.4.2.
The implementation ofÛ ST i j relies on the dipole-dipole interaction Hamiltonian, which we approximated withĤ ST i j . From inspecting the effective dipole-dipole interaction Hamiltonian, equation (14), we see that undesired pair-creation terms can be eliminated by tuning |ω a − ω b | ≫ JV i j . Furthermore, we can isolate the b hopping (or flip-flop) terms, by using electric and magnetic gradient fields that shift the a hopping processes out of resonance by making ω a spatially dependent. Finally, density-density (or Ising) interactions can be neglected in the limit n i , n j ≪ J.
In figure 7(b), we quantify the accuracy of the state transfer protocol by plotting the gate fidelity F n i ,n j = ⟨ψ T |exp(−iĤ t T)|ψ 0 ⟩. Here we take several possible initial states and we include only Ising and b flip-flop terms inĤ t . We find that a high fidelity F n i ,n j ≳ 0.99 can be obtained for n i , n j ≲ 3 and J = 50. The main source of error is originating from the finite spin length, which affects both the Ising and the flip-flop terms via the Holstein-Primakoff transformation. The fidelity may be further improved by using optimal control techniques [149] or by adapting the scheme demonstrated in [87,88] to our system that could suppress the effect of the Ising terms.

Conclusions
In this work, we have presented an hardware-efficient quantum simulation Rydberg toolbox for atoms in tweezer arrays based on the high-dimensional manifold of states with principal quantum number n in the regime of linear Stark and Zeeman effect. We have exploited the SO(4) symmetry of the hydrogen-like states to characterize these energy levels in terms of two angular momenta of large length, which have been conveniently used to represent the action of external fields (static electric and magnetic fields as well as optical and MW fields).
Within this formalism, we have shown that the many-body problem is represented by generalized 'large-spin' Heisenberg models whose large local Hilbert space can be used to encode conjugate (continuous) variables to simulate QFTs. In particular, we have illustrated how to realize the 1D SG model in the massless and massive case. Besides the quantum simulation applications, these high-dimensional manifolds and their control via external fields can also be exploited for qudit-based quantum information processing, which we have exemplified with an entangling gate and a state-transfer protocol involving states in the vicinity of the circular level. For such applications it is crucial to quantitatively verify the Rydberg simulator, which defines one of the major challenges in quantum simulation [150]. One possible approach to analyze the experimentally realized processor is via Hamiltonian learning [151].
Our results offer the opportunity to simulate more general QFTs in higher dimensions, gauge degrees of freedom with large occupation numbers or large-spin Heisenberg models in arbitrary lattice geometries. As our analysis has focused on encoding a local Hilbert space on a single n-manifold, a possible extension is to study models where multiple manifolds are considered or where the electron spin is also included. Furthermore, other possible future directions include to employ these Rydberg states to encode synthetic dimensions [15] or to achieve quantum-enhanced sensing [102,152].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Algebraic solution of the Hydrogen problem
Here, we summarize further details of the algebraic solution of the hydrogen atom outlined in the main text section 2.1.
As already pointed out in the main text, the non-relativistic hydrogen Hamiltonian readŝ H H =p 2 /2 − e/4πϵ 0 |r|, and a convenient approach to solve for the eigenvalues is to use constants of motion. In particular, the orbital angular momentumL =r ×p and the Runge-Lenz (RL) vector are commuting withĤ H . At the classical level, the RL vector is a conserved quantity of the Kepler problem, it is aligned with the major axis of a closed orbit and its magnitude is proportional to its eccentricity. Therefore, the RL vector is orthogonal to the orbital angular momentum vector, which is also inherited by the quantum problem, namelyL ·Â 0 = 0, and will be used below.
As we are interested in bound state solutions (E < 0), it is useful to redefine the RL vector such that we obtain commutation relations of the form = iϵ i jkÂk that form the SO(4) Lie algebra. It is convenient to change the basis of the algebra to which doubly covers SO(4). The commutation relations of the two angular momentaĴ a andĴ b are From the nontrivial relationĤ H (L 2 +Â 2 + 1) = −Ry, it is possible to rewrite the hydrogen atom Hamiltonian asĤ H = −Ry/[2(Ĵ 2 a +Ĵ 2 b ) + 1]. Furthermore, the orthogonality between the orbital angular momentum and the RL vector (L ·Â = 0) also fixes the length of the two angular momenta, J 2 a =Ĵ 2 b = (L 2 +Â 2 )/4, such that an obvious basis for the hydrogen atom eigenstates is |J, m a , m b ⟩, wherê it is straightforward to notice that the matrix elements of the left-hand side of this equation within a specific J manifold vanish because they are all eigenstates ofĤ H with the same energy E n . As a result, we obtain that the electric dipole operatorμ = −er projected on a specific n manifold takes the simple form see [71][72][73] for further details of the derivation.

Appendix B. Decoherence
In this appendix we discuss decoherence channels of the Rydberg levels discussed in the main text section 2.4. In particular, we present spontaneous emission rates and discuss motional decoherence.

B.1. Spontaneous emission
In a cryogenic environment (i.e. for temperatures below a few K [61]) the lifetime of the Rydberg levels is fundamentally limited by spontaneous emission to lower lying levels. The decay rate from an initial state |i⟩ with energy ω i to a final state | f ⟩ with energy ω f is determined by the spontaneous emission rates where ω i, f = ω i − ω f is the transition frequency and ⟨i|μ| f ⟩ the dipole transition matrix element. The total lifetime τ i of |i⟩ is then given by τ −1 In figure 8(a) we present the numerically calculated lifetimes for the manifold H a (introduced in the main text section 3) for different principal quantum numbers n. The relatively large lifetime of the circular level (m a = (n − 1)/2) is to a large extent preserved for states with lower orbital angular momentum components. In figure 8(b) we monitor the scaling of the lifetime with n. The lifetime scaling of the circular level is n 5 [61] and it gradually changes to the numerically calculated n 4 scaling for states with low orbital angular momentum components. The n 4 lifetime scaling of states with low angular momentum components can be understood qualitatively as follows: the lifetime of low orbital angular momentum states typically scales as n 3 [96]; if these states are not affected by quantum defects, in the presence of an external electric field they are distributed among O(n) Stark states [69] and, therefore, their lifetime is enhanced by a factor n.

B.2. Motional decoherence
Here we summarize decoherence mechanisms associated with motional effects. First we discuss magic (state independent) trapping of the Rydberg states considered in this paper and then estimate motional excitation rates due to dipole-dipole forces.
It has been proposed [62,103] and experimentally demonstrated [154][155][156] that Rydberg states can be spatially confined using ponderomotive trapping techniques. For hollow bottle beam tweezers the size of the Rydberg wavefunctions is typically comparable to the radial extent of the tweezer trap and, hence, the ponderomotive potential typically depends on the Rydberg state. Therefore, finding magical trapping conditions that avoid internal state dephasing is challenging with ponderomotive traps.
An alternatively route available for alkaline earth and lanthanide atoms is to use the polarizability of the optically active core to trap the atom with a red-detuned Gaussian tweezer beam [66][67][68]. In particular, when the trapping light is almost resonant with a core transition, the core polarizability becomes much greater than the ponderomotive potential (from the trapping light) experienced by the Rydberg electron. Therefore the trapping potential becomes to a large extent independent of the Rydberg state [62,68], which enables magic trapping for many different Rydberg states.
Next, we estimate the creation of vibrational excitations due to dipole-dipole forces acting on the nucleus inside a harmonic trap with frequency ω t . An interaction with effective strength V i j J(J + 1), as considered in section 4 of the main text, gives rise to an additional potential U(R i ) = 3R i V i j J(J + 1)/|R i j | acting on the ith atom, where R i denotes the position of the ith nucleus measured from the center of the trap. If the strength of U(R i ) at the characteristic trapping length scale l h = 1/(M a ω t ) (M a is the mass of the atom) is small compared to ω t , i.e. U(l t ) ≪ ω t , then at short times (ω t t/(2π) ≪ 1) excitations are created as n(t) = U(l t ) 2 t 2 /8, where we assumed that the nucleus is initialized in the ground state of the oscillator, which can be achieved by sideband cooling [50,51].
As an example we consider the parameters used to realize the SG model from section 4 of the main text, for which we have an interaction strength of V i j J(J + 1) = 2π × 300 kHz for n = 41 and an inter-atomic separation of R i j = 17 µm. Assuming a typical trapping frequency of ω t = 2π × 40 kHz and strontium atoms (M a ≈ 87u) gives rise to l t = 54 nm for which we get U(l t ) ≈ 2π × 3 kHz ≪ ω t . For these parameter values, n(10t c ) ≈ 0.05 excitations are created after ten interaction cycles of duration t c = 2π/[V i j J(J + 1)].

Appendix C. Corrections to dipole-dipole interactions
In the following appendix we discuss corrections to the dipole-dipole interaction HamiltonianĤ i j dd introduced in the main text section 2.3.

C.1. Dipole-dipole interactions
The replacement of the dipole operator by the angular momentum operatorsĴ a andĴ a (see main text equation (3)) underlies the assumption that the eigenstates of the hydrogen Hamiltonian with principal quantum number n are not mixed with eigenstates from different n ′ ̸ = n manifolds in the presence of a static external electric field. For electric fields below the Ingris Teller limit (F < F IT ) this mixing of different n-manifolds is small and, therefore, the dipole operator and also the dipole-dipole interactions are well described by the formalism introduced in the main text section 2.1. In figures 9(a) and (b) the dipole-dipole interaction matrix elements including and excluding manifold mixing are compared to each other for different principal quantum numbers. In particular, for the termĴ a,+ , relevant for the realization of the Sine-Gordon and massive Schwinger model, the relative mismatch due to manifold mixing is 10 −3 for n = 41.

C.2. Van der Waals interactions
Non-zero dipole transition matrix elements between different n-manifolds give rise to dipole-dipole couplings from pair states with principal quantum number n to pair states with different principal quantum numbers n ′ , n ′ ′ ̸ = n. Since the single particle energy mismatch between the pair states is typically considerably larger than the interaction matrix elements, these couplings give rise to second order corrections, namely Van der Waals (VdW) interactions. In figures 9(a) and (b) we display the strength of the VdW interaction matrix elements as a function of interatomic separation R ij for different n, calculated with degenerate perturbation theory as discussed in the appendix of [61] or in the main text of [62]. For R i j ≳ 10 µm and n = 41, as considered for the analog simulation of QFTs (see main text section 4), VdW interactions are considerably smaller than direct dipole-dipole interactions, and can therefore safely be neglected.
For the edge manifold H a (see main text equation (10)), additional VdW terms, mediated by pair states contained within the n manifold under consideration, have to be taken into account. Since the single particle energy imbalance is controlled by the static electric and magnetic fields these couplings can be rendered negligible for field strengths and interatomic separations relevant for this manuscript, see figures 9(a) and (b).

Appendix D. Ponderomotive coupling
In the following we discuss the ponderomotive manipulation techniques introduced in section 3.4 of the main text in more detail. The discussion is structured as follows: first we introduce Laguerre-Gauss (LG) laser beams and then we discuss the ponderomotive potential generated by interfering two LG beams.

D.1. Laguerre-Gauss laser beam
The vector potential of a linearly polarized LG laser beam with orbital angular momentum δm and mode number p, propagating along the z-axis, is in the Lorenz gauge given by [157,158] A p,δm (r, t) =ê x A p,δm u p,δm (r)e iα e −iδmϕ e −i(kz−ωt) /2 + c.c., where we used cylindrical coordinates (ρ = x 2 + y 2 and tan ϕ = y/x).
Here ω is the light frequency, k = 2π/λ is the magnitude of the corresponding wave vector, λ is the wavelength and α is a phase. Moreover, A p,δm denotes the field amplitude and u p,δm (r) is a dimensionless field distribution function satisfying the Helmholtz equation. In the paraxial limit u p,δm (r) can be expressed in terms of associated Laguerre polynomials L (|δm|) p (x) as The beam waist is defined as w(z) = w 0 1 + z 2 /z 2 R , where w 0 denotes the waist in the focal plane and z R = π w 2 0 /λ is the Rayleigh length. The normalization constant C p,δm = p! 2/( p + |δm|)! is chosen such that´d A|u(ρ, z)| 2 = w 2 0 π. Due to this normalization the time averaged electric beam power P = cϵ 0 ω 2 |A p,δm | 2 w 2 0 π/4 only depends on the field amplitude and not on the LG mode numbers.

D.2. Ponderomotive potential
The ponderomotive potential that a highly excited Rydberg electron experiences from a fast oscillating laser field is given by [159] U PM (R + r, t) = e 2 |A(R + r, t)| 2 /(2m e ), where e and m e are the electron charge and mass, respectively. The vector potential of the laser field is evaluated at the position of the electron in the laboratory frame, i.e. R is the position of the core and r is the position of the electron relative to the core. In the tight trapping limit, when the size of the vibrational core wave-packet is much smaller than the LG waist w 0 , R can be replaced by its mean value. Furthermore, if the laser beam is focused onto the center of the trap, i.e. ⟨R⟩ = 0, which can be achieved by re-using the trapping light optics, the ponderomotive potential becomes approximately independent of the core position U PM (R + r, t) ≈ U PM (r, t).
As already mentioned in the main text, the core idea to couple states with multiple orbital angular momentum κ ⩾ 1 difference is to use the ponderomotive potential of two co-propagating LG laser beams with different oscillation frequencies and non-zero orbital angular momentum [62,[103][104][105]. Using equation (D1), the vector potential of the two co-propagating LG beams is given by where the subscripts 1 and 2 are labels for the first and second LG beam. The field distribution functions u δm 1(2) are in the most general case a superposition of multiple LG modes u p,δm 1(2) from equation (D2), with fixed δm 1(2) . For simplicity, the field amplitudes are assumed to be the same for both beams. The corresponding ponderomotive potential, after dropping fast oscillating terms, becomes where δω = ω 1 − ω 2 is the oscillation frequency of the interference term and κ = δm 1 + δm 2 is the transferred orbital AM. The spatial shapes of the constant and oscillatory terms are given by with U 0 = e 2 P tot /(16π 3 m e ϵ 0 c 3 ) (λ 2 /w 2 0 ), where P tot is the combined power of both beams. The contributions of the ponderomotive potential are twofold: the time independent term gives rise to a level dependent energy shift and the time dependent term can couple states that differ by κ orbital angular momenta.
Let us start by analyzing the time dependent contribution. In particular, we are interested in the transition matrix elements of the form that satisfy orbital angular momentum conservation κ a + κ b = κ. Note, κ a(b) can take on negative and positive values. For the determination of the transition matrix element the functional behavior of U κ (ρ, z) is of crucial importance. In the limit where the beam waist w 0 is greater than the typical extension of the Rydberg wave-function, see figure 10(a), we can expand the ponderomotive potential in the focal plane z = 0 around the phase vortex at ρ = 0, and obtain U κ (ρ, z = 0) ≈ C κ+2η ρ |κ|+2η , where the leading order |κ| + 2η, with η ⩾ 0, is determined by the LG beam modes. We now claim that, within the manifold H t , the matrix elements of the ponderomotive potential near the phase vortex c κ+2η ρ |κ|+2η e −iκϕ are exactly proportional to those of (Ĵ a,+ ) κa (Ĵ b,+ ) κ b , with |κ a | + |κ b | = κ + 2η and where we use the convention (J a(b),+ ) −κ a(b) = (J a(b),− ) κ a(b) . More specifically, with As it is quite tedious to find an exact algebraic relation, we instead confirmed the relation (D7) numerically and found perfect agreement. In figure 10(b) we show such a numerical comparison for the special case of κ b = 0. Figure 10(c) showcases the scaling with n. Note, the relation (D7) was inspired by the findings of [72], who derived similar results for the matrix elements of z 2 .
To analyze the ponderomotive transition matrix elements U κa,κ b ma,m b from equation (D6) it is important to take into account the full shape of U κ (ρ, z). As the functional behavior of U κ (ρ, z) is determined by the LG beam modes, it is sufficient to analyze u δm1 (ρ, 0) u * δm2 (ρ, 0). A series expansion of the field distribution functions at z = 0 and around ρ = 0 yields in the most general case u δm1 (ρ, 0) u * δm2 (ρ, 0) = c κ+2η ρ |κ|+2η + c (2) κ+2η ρ |κ|+2η+2 + . . . . As only the leading term c κ+2η ρ |κ|+2η generates the desired nonlinearity, higher orders have to be suppressed, which becomes essential for larger principal quantum numbers, as the increased size of the Rydberg-orbit starts to explore the full beam shape, see figure 10(a).
To suppress the higher orders we consider further beam shaping beyond a single LG beam mode, which can be achieved for instance by a spatial light modulator [160]. Alternatively one could use super-positions of N different LG modes, i.e. u δm 1(2) = N−1 p=0 a 1(2) p u p,δm 1(2) , with (a 1(2) p ) 2 = 1, which allows to eliminate 2N − 1 higher orders u δm1 (ρ, 0) u * δm2 (ρ, 0) = c κ+2η ρ |κ|+2η (1 + O(ρ 2N )). Hence, the scheme presented here becomes applicable for larger n without increasing the beam waist w 0 . However, with increasing N the maximum of the ponderomotive potential moves further away from the vortex and, therefore, U κa,κ b ma,m b for fixed beam power is reduced.
For the transition matrix elements presented in figure 2(b) in the main text, where we used the simplified notation U κa,ma,m b ≡ U κa,κ b =0 ma,m b , the LG beams are chosen as u δm1 (ρ, z) = u δm2 (ρ, z), with δm 1 = δm 2 = κ a /2 and η = 0. The number of modes used for κ a = 2 and 4 is N = 2 and for κ a = 6 it is N = 3.
For the same LG laser beams the transition matrix elements U κa,κ b =0 ma=0,m b =J as a function of the n are presented in figure 10(d). For large principal quantum numbers the transition matrix elements deviate from the predicted n 2κa scaling as the large Rydberg wavefunction probes the whole ponderomotive potential. Note, the scaling n 2κa is obtained from equation (D7) and the fact that the angular momentum transition matrix element U κa,κ b =0 ma=0,m b =J scales for n → ∞ as n κa . We furthermore point out that by tuning the relative laser phase α 1 − α 2 the coupling matrix elements U κa,κ b ma=0,m b can be complex. Moreover, a pair of co-propagating LG beams generates a single nonlinear term and in order to drive the engineered transitions resonantly the laser frequencies have to satisfy δω = −κ a ω a + κ b ω b .
Let us briefly return to the first term of equation (D5). In contrast to the time dependent term, the constant part of the ponderomotive potential generates state dependent energy shifts. These shifts can in principle be used to engineer nonlinearities like higher powers of J a,z and J b,z . However, the strength of the shifts turns out to be relatively weak for reasonable laser powers of hundredths of mW per site, when compared to the nonlinearities generated by off-resonant MW coupling.
Finally, we remark on decoherence effects induced by the LG beams due to Thomson scattering of the nearly free Rydberg electron [62]. The scattering rate is given by Γ T = Iσ T /ω, where I = 2P/(πω 2 0 ) is the mean beam intensity and σ T = (8π/3)e 4 /(4πϵ 0 m e c 2 ) 2 is the Thomson scattering cross-section [161]. For λ = 1300 nm and w 0 = λ/2 we obtain Γ T /P = 2π × 0.7 Hz/mW, which is typically much less than the achieved coupling rates, see figure 2(b). Furthermore, we note that for the range of parameters considered here the Rydberg electron does not probe the high intensity regions (see figure 10(a)) and, therefore, the calculated scattering rates are overestimated.

Appendix E. Offresonant MW coupling
This appendix provides further details on the engineered spin-squeezing termĤ SQ = χ J 2 z,a introduced in the main text equation (15). To engineer this nonlinearity we consider z-polarized MW light that couples the manifold H a with principal quantum number n to the levels of the manifold n ′ = n + δn, with δn > 0. If the MW field is blue detuned with respect to the eigen-energies of the n ′ manifold, only couplings to the manifold H a of n ′ , i.e. the states |ψ ′ a ma ⟩, are important, see main text figure 2(c). The underlying Hamiltonian is that of multiple independent two-level systems, one for each m a , and is in a frame rotating at the MW frequency given bŷ with m ′ a = m a − δn/2. Here, ∆ ma = ∆ 0 − m a δω a is the level dependent detuning given in terms of the bare MW detuning ∆ 0 and the differential Stark shift δω a = ω a (n ′ ) − ω a (n) = 3 ea 0 δnF/2, see figure (2)c. Moreover, Ω ma = 2F MW ⟨ψ ′ a ma |μ z |ψ a ma ⟩ is the Rabi frequency, where F MW is the electric field amplitude of the MW.
As already mentioned in the main text, in the far off-resonant limit (∆ ma ≫ Ω ma ), the dressed eigenstates | ψ a ma ⟩, which are adiabatically connected to |ψ a ma ⟩ for Ω ma → 0, experience a state dependent AC-Stark shift Ω 2 ma / 4∆ ma . Furthermore, in the limit ∆ 0 ≫ m a δω a and due to the smooth dependence of Ω ma on m a the AC-Stark shift is expandable in powers of m a , where prefactors χ and ξ are in the following determined numerically. Note, the constant and linear part are absorbed in ∆ MW and δω a , respectively. The strength χ of the nonlinearity is limited by the electric field ionization threshold, for which a lower bound is the Ingris Teller limit (|F MW | < F IT ). This limitation originates from the fact that the MW also couples to the permanent dipole moments ⟨ψ a ma |μ z |ψ a ma ⟩ and thus induces a modulation of ω S , if |F MW | exceeds F IT the modulation leads to MW assisted ionization [162]. For |F MW | < F IT , the modulation averages out because the MW frequency is much larger than the modulation amplitude 2F MW ⟨ψ a ma |μ z |ψ a ma ⟩ for the range of parameters considered in this paper.
Within this limitation achievable prefactors χ can be on the order of hundreds of 2π× kHz, see figure 11(a) upper panel. We point out that for every value of χ there exists a point in ∆ 0 , Ω 0 parameter space where the cubic contribution ξ of equation (E2) vanishes, see figure 11(a) middle panel. Note, second order Stark corrections also give rise to nonlinearities (see figure 11(a) upper panel Ω 0 = 0), but these rapidly decrease with increasing n as 1/n −6 [72,163]. We also note that a misalignment of the MW polarization induces Raman transitions between neighboring levels, which are however far off-resonant in the limit |ω a | ≫ Ω ma . As a result, χ is largely insensitive to polarization errors (below 10% of F MW ).
Due to the MW coupling the dressed eigen-states | ψ a ma ⟩ acquire character of the n ′ manifold. The corresponding purity is quantified by |c ma | 2 = |⟨ψ a ma | ψ a ma ⟩| 2 and decreases as m a → −J, since the effective MW detuning is reduced, see lower panel of figure 11(d). Furthermore, in the lowest panel of figure 11(a) we monitor |c ma | 2 for m a = 0 and different values of the detuning and driving strength. The admixture of Rydberg levels from n ′ to | ψ a ma ⟩ modifies the dipole transition matrix elements, which thus deviate from the angular momentum formalism introduced in equation (3). As an example of this deviation we present ⟨ ψ a ma+1 |μ + | ψ a ma ⟩, whereμ + =μ x + iμ y for different MW parameters in figure 11(b). The modifications of the dipole matrix elements are expected to be in practice negligible for the quantum simulation application presented in the main text section 4.

Appendix F. Exact results for the SG model
For completeness, we provide exact expressions for the spectrum of the SG model in this appendix, see e.g. [126] and references therein. In the massive phase (β 2 < 8π), there exist solitons with mass where we abbreviated ξ = β 2 /(8π − β 2 ). For β 2 < 4π, the theory has additional bound states, called breathers. The number of different breathers depends the value of ξ and they can be labeled by the integers n = 1, 2, . . . , ⌊1/ξ ⌋. The mass of the nth breather is given by

Appendix G. Implementing the sine-Gordon model with large spins
In the following we outline how to simulate the SG model (22) with the Rydberg simulator. As pointed out in the main text, we start with a one-dimensional array of N equally spaced Rydberg atoms. Isolating the subset of states H a (see section 3), the system is described by the many-body Hamiltonian in equation (16). Within the large spin limit (J → ∞), we can employ the identification shown in equation (20) for every lattice site i. In order to obtain the desired limit J → ∞, we rescale the interaction V ′ i j = V i j J(J + 1), while keeping ∆ a = 0, and λ ′ = −2λ κ [J(J + 1)] κ/2 . For the special case κ = 1, this can be achieved by using a MW field as discussed before. The ponderomotive technique instead allows to achieve couplings with κ > 1, which gives access to a larger parameter regime of the SG model, as we discuss below. The resulting continuous variable lattice model readsĤ where the Ising interactions are neglected in the limit J → ∞. The continuum SG model arises by keeping only the nearest-neighbor terms V ′ i,i+1 ≡ V ′ nn . To see this, we rescale the fields as κφ → βφ,π/κ → (ℓ/β)π and the Hamiltonian asĤ → (2χκ 2 ℓ/β 2 )Ĥ. Here, we introduced a short-distance scale ℓ that sets a UV-cutoff Λ ∝ 1/ℓ [136]. As a consequence, we find that the lattice regularizationĤ (lat) SG approximatesĤ SG upon identifying the atomic couplings with SG parameters according to To summarize, our implementation of the SG model is formally valid in the continuum limit (ℓM 0 = κ 2λ ′ /V ′ nn → 0), together with sufficiently large spin length (J → ∞) and a large number of Rydberg atoms (N → ∞). This requires tuning the atomic parameters according to equations (G2), which in principle allows to probe all regimes of the theory with 0 < β 2 = 2κ 2 χ/V ′ nn < ∞. In practice, the range of parameters is mainly limited by the maximum value of χ ≲ 2π × 200 kHz. For example, κ = 1, J = 20 and choosing χ = 2π × 200 kHz, V ′ nn = 2π × 100 kHz gives β 2 ≈ π. Moreover, for λ ′ = 2π × 5 kHz, we then find ℓM 0 ≈ 0.3. From these estimates, we thus expect that several tens of atoms are sufficient to simulate the SG QFT with reasonable accuracy up to β 2 ≲ π, that is deep in the massive phase. Larger values of β 2 can be achieved by using the ponderomotive coupling with κ > 1. For example, taking κ = 4, χ = 2π × 200 kHz, V ′ nn = 2π × 300 kHz and λ ′ = 2π × 3 kHz gives β 2 ≳ 8π and ℓM 0 ≈ 0.6. With these parameters one can then investigate the critical properties of the model across the quantum phase transition.
where the quantum state of the oscillator b is completely transferred to a. The generalization of this protocol to the Rydberg manifold is straightforward.