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

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, ${\hat{H}}_{\mathrm{H}} = {\hat{\mathbf{p}}^2}/{(2m_r)} - {e^2}/(4\pi{\epsilon}_0|{\hat{\mathbf{r}}}|)$, where we ignore electronic spin degrees of freedom and m_r is the reduced electron mass, −e the electron charge, ₀ the vacuum permittivity and $\hat{\mathbf{p}}$ and ${\hat{\mathbf{r}}}$ are the momentum and position operators, respectively.

A convenient approach to solve for the bound states of ${\hat{H}}_{\mathrm{H}}$ is to exploit the conservation of orbital angular momentum ${\hat{\mathbf{L}}}$ and Runge–Lenz (RL) vector ${\hat{{\mathbf{A}}}}$ (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

$$\begin{align} {\hat{\mathbf{J}}}_a = \frac{1}{2}\left( {\hat{\mathbf{L}}} -{\hat{\mathbf{A}}}\right)\,\,\mathrm{and}\,\,{\hat{\mathbf{J}}}_b = \frac{1}{2}\left( {\hat{\mathbf{L}}} +{\hat{\mathbf{A}}}\right)\,, \end{align} \tag{ 1 }$$

and therefore obey $[{\hat{J}}_{a,i}, {\hat{J}}_{a,j}] = i{\epsilon}_{ijk} {\hat{J}}_{a,k}$ and $[{\hat{J}}_{b,i},{\hat{J}}_{b,j}] = i{\epsilon}_{ijk} {\hat{J}}_{b,k}$, furthermore, the angular momentum raising and lowering operators are defined as ${\hat{J}}_{a(b),\pm} = {\hat{J}}_{a(b),x} \pm i {\hat{J}}_{a(b),y}$. Note, in this paper we use units where $\hbar = 1$. Due to the orthogonality between the orbital angular momentum and the RL vector (${\hat{\mathbf{L}}} \cdot {\hat{{\mathbf{A}}}} = 0$), the lengths of the two angular momenta are constrained, ${{\hat{\mathbf{J}}}_{a}^{2}} = {{\hat{\mathbf{J}}}_{b}^{2}}$.

Within this formalism, the Hamiltonian can be recast in the simple form ${\hat{H}}_{\mathrm{H}} = -{\mathrm{Ry}}/[2({\hat{\mathbf{J}}}_a^2 + {\hat{\mathbf{J}}}_b^2)+1]$, see appendix A, where ${\mathrm{Ry}} = e^2/(8\pi{\epsilon}_0a_0)$ is the Rydberg energy and $a_0 = 4\pi{\epsilon}_0/(e^2m_r)$ the Bohr radius. A natural basis for the hydrogen atom eigenstates is therefore given by states $\vert J,m_a,m_b\rangle$, where

$$\begin{align} {\hat{\mathbf{J}}}_{a\,(b)}^2\vert J,m_a,m_b\rangle = J(J+1)\vert J,m_a,m_b\rangle \quad \mathrm{and} \quad {\hat{J}}_{a\,(b),z} \vert J,m_a,m_b\rangle = m_{a\,(b)}\vert J,m_a,m_b\rangle, \end{align} \tag{ 2 }$$

with $J = 0,1/2,1,\dots$ and $m_{a\,(b)}\in\{-J,\,-J+1,\,\dots,\,J\}$. The corresponding energies of the Hamiltonian take the form $E_n = -{\mathrm{Ry}}/(2J+1)^2 = -{\mathrm{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² degenerate states $\vert J,m_a,m_b\rangle$, 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 $\hat{\boldsymbol{\mu}} = - e{\hat{\mathbf{r}}}$ whose transition matrix elements, within a single n manifold, can be identified as

$$\begin{align} \hat{\boldsymbol{\mu}} = -3n\, e a_0\, {\hat{{\mathbf{A}}}}/2 = 3n\,e a_0\,\big({\hat{\mathbf{J}}}_a - {\hat{\mathbf{J}}}_b\big)/2\,, \end{align} \tag{ 3 }$$

see [1, 71–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 ${\hat{\mathbf{J}}}_a$ and ${\hat{\mathbf{J}}}_b$.

The coupling of the electron to weak static external electric ${{\mathbf{F}}}$ and magnetic B fields lifts the degeneracy of states with fixed n, and is described by the Hamiltonian

$$\begin{align} {\hat{H}}_\mathrm{\mathbf{B},{\mathbf{F}}} = \mu_B {\hat{\mathbf{L}}}\cdot\mathbf{B}-\hat{\boldsymbol{\mu}}\cdot{\mathbf{F}}\,, \end{align} \tag{ 4 }$$

where $\mu_B = e/(2m_r)$ is the Bohr magneton. For a manifold of hydrogen-like states with a specific n, the ${\hat{\mathbf{L}}}$ and $\hat{\boldsymbol{\mu}}$ can be replaced by angular momentum operators according to equations (1) and (3), respectively. The replacement of $\hat{\boldsymbol{\mu}}$ is valid below the Ingris-Teller limit $|{\boldsymbol{F}}| \lt {F_{\mathrm{IT}}} = 2{\mathrm{Ry}}/(3\,ea_0\,n^5)$, i.e. when couplings to adjacent $n^{^{\prime}} = n\pm 1$ manifolds are weak with respect to the gap $E_{n}-{E_{n\pm 1}}$ [3, 74, 75]. Analogously, diamagnetic couplings are negligible in the limit $|\mathbf{B}|\lt 2{\mathrm{Ry}}/(\mu_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, ${{\mathbf{F}}} = F {\mathbf{e}}_z$ and ${\mathbf{B}} = B {\mathbf{e}}_z$, see figure 1(a), such that the projection $m_l = m_a + m_b$ of the orbital angular momentum $\hat L_z$ remains a good quantum number. Within the limits discussed in the previous paragraph the Hamiltonian from equation (4) becomes

$$\begin{align} {\hat{H}}_n = -{\omega}_a\,{\hat{J}}_{a,z} + {\omega}_b\,{\hat{J}}_{b,z}\,, \end{align} \tag{ 5 }$$

with ${\omega}_a = {\omega}_S - {\omega}_Z$ and ${\omega}_b = {\omega}_S + {\omega}_Z$, where ${\omega}_{S} = 3n\,e a_0\,F/2$ is the Stark splitting and ${\omega}_Z = \mu_BB$ 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 ⁸⁷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⁻¹, and for magnetic fields of hundreds of Gauss, and principal quantum numbers in the range $30 \lt n \lt 70$, we have ${\omega}_Z\lt{\omega}_S$ such that ${\omega}_{a(b)}\gt0$. Within this parameter regime the absolute values of ω_a and ω_b can be of the order of $2\pi\times(10^2 - 10^3)\, \mathrm{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 ${\mathbf{F}}_a = F_a{\mathbf{e}}_-{\textrm{exp}}(-i{\omega}_{a}^{\scriptscriptstyle{\mathrm{MW}}}t) + \mathrm{c.c.}$, and ${\mathbf{F}}_b = F_b{\mathbf{e}}_+{\textrm{exp}}(-i{\omega}_{b}^{\scriptscriptstyle{\mathrm{MW}}}t) + \mathrm{c.c.}$, with frequencies ${\omega}_{a(b)}^{\scriptscriptstyle{\mathrm{MW}}}$, electric field amplitudes $F_{a(b)}$ and polarizations $\mathbf{e}_{\pm} = \mp(\mathbf{e}_{x} \pm i \mathbf{e}_{y})/\sqrt{2}$. The electric field of the MWs couple to the dipole operator $\hat{\boldsymbol{\mu}}$, 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 operators

$$\begin{align} {\hat{H}}_{\mathrm{MW}}=\sum_{\sigma=a,b}\left[-\Delta_{\sigma}\,{\hat{J}}_{\sigma,z} + \left({\Omega}_\sigma{\hat{J}}_{\sigma,+} + {\Omega}^*_\sigma{\hat{J}}_{\sigma,-} \right)\right], \end{align} \tag{ 6 }$$

where we dropped fast oscillating terms and the rotating frame transformation is defined by $\hat R = {\textrm{exp}}[-i({\omega}_a^{\scriptscriptstyle{\mathrm{MW}}}{\hat{J}}_{a,z}-{\omega}_{b}^{\scriptscriptstyle{\mathrm{MW}}}{\hat{J}}_{b,z})t]$. Here, $\Delta_{a} = {\omega}_a-{\omega}_{a}^{\scriptscriptstyle{\mathrm{MW}}}$ ($\Delta_{b} = {\omega}_{b}^{\scriptscriptstyle{\mathrm{MW}}}-{\omega}_{b}$) and ${\Omega}_{a(b)} = 3n\,ea_0\, F_{a(b)}/(2\sqrt{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 ${\omega}_a-{\omega}_b$ via the magnetic field.

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 ${\mathbf{R}}_{ij} = R_{ij}\, \mathbf{e}_{ij}$, here R_ij is the inter-atomic separation and $\mathbf{e}_{ij} = ({\textrm{cos}}{\phi}_{ij}\,{\textrm{sin}}{\theta}_{ij}, {\textrm{sin}}{\phi}_{ij} {\textrm{sin}}{\theta}_{ij},$ ${\textrm{cos}}{\theta}_{ij})^\mathrm{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]

$$\begin{align*} {\hat{H}}^{ij}_{\mathrm{dd}} = \frac{1}{4\pi{\epsilon}_0 R_{ij}^3} \left[\hat{\boldsymbol{\mu}}^{(i)}\cdot \hat{\boldsymbol{\mu}}^{(j)} - 3 (\hat{\boldsymbol{\mu}}^{(i)}\cdot \mathbf{e}_{ij}) (\hat{\boldsymbol{\mu}}^{(j)}\cdot \mathbf{e}_{ij}) \right], \end{align*}$$

where $\hat{\boldsymbol{\mu}}^{(i)}$ and $\hat{\boldsymbol{\mu}}^{(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^{^{\prime}} = n\pm 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 as

$$\begin{align} {\hat{H}}^{ij}_{\mathrm{dd}}= V_{ij} \Big[&{\hat{\mathbf{J}}}^{(i)}_a\cdot {\hat{\mathbf{J}}}^{(j)}_a - 3 ({\hat{\mathbf{J}}}^{(i)}_a\cdot \mathbf{e}_{ij})({\hat{\mathbf{J}}}^{(j)}_a\cdot \mathbf{e}_{ij})\nonumber\\+ &{\hat{\mathbf{J}}}^{(i)}_b\cdot {\hat{\mathbf{J}}}^{(j)}_b - 3 ({\hat{\mathbf{J}}}^{(i)}_b\cdot \mathbf{e}_{ij}) ({\hat{\mathbf{J}}}^{(j)}_b\cdot \mathbf{e}_{ij})\nonumber\\-&{\hat{\mathbf{J}}}^{(i)}_a\cdot {\hat{\mathbf{J}}}^{(j)}_b + 3 ({\hat{\mathbf{J}}}^{(i)}_a\cdot \mathbf{e}_{ij})({\hat{\mathbf{J}}}^{(j)}_b\cdot \mathbf{e}_{ij})\nonumber\\-& {\hat{\mathbf{J}}}^{(i)}_b\cdot {\hat{\mathbf{J}}}^{(j)}_a + 3 ({\hat{\mathbf{J}}}^{(i)}_b\cdot \mathbf{e}_{ij}) ({\hat{\mathbf{J}}}^{(j)}_a\cdot \mathbf{e}_{ij})\Big], \end{align} \tag{ 7 }$$

with $V_{ij} = (3n\,e a_0)^2/(16\pi{\epsilon}_0R_{ij}^3)$. Corrections to the above Hamiltonian, like van der Waals interactions mediated by states with principal quantum number $n^{^{\prime}}\neq 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–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.

The Rydberg manifolds of states with $n\gg 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} = -{\mathrm{Ry}}/(n-\delta_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]

$$\begin{align} {\hat{H}}_n^\mathrm{\,QD} = {\hat{H}}_n +\sum_{\substack{|m_l|\leqslant l \\~~ l\lt l^*}}\vert n,l,m_l\rangle \langle n,l,m_l\vert (E_{n,l} - E_n), \end{align} \tag{ 8 }$$

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 ${\hat{H}}_n^\mathrm{QD}$ do not form a perfectly regular rhombus structure anymore because the inner region $|{m_{l}}| \lt {l^{*}}$ presents an irregular spacing [3, 94]. However, for $n\gg 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, $\vert J, m_a, m_b\rangle\otimes\vert s,m_s\rangle$ with the electron spin state $\vert s,m_s\rangle = \vert 1/2, \pm1/2\rangle$, 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 $\vert s,m_s\rangle = \vert 1/2, 1/2\rangle$.

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–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 \gt 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 ${\omega}_{a(b)}$, and therefore enter as dephasing. Experimentally demonstrated electric field stabilizations on the level of tens of µV cm⁻¹ [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–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.

The first and largest manifold that we consider is the triangular manifold defined for a fixed n by

$$\begin{align} \mathcal H^t = \big\{\vert \psi^t_{m_a,m_b}\rangle\big\},\, \mathrm{with}\,\, |\psi^t_{m_a,m_b}\rangle \equiv |J, m_a, m_b \rangle\,, \end{align} \tag{ 9 }$$

and $m_l\geqslant 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.

The second manifold that we consider is a subset of the triangular manifold $\mathcal H^t$ and is defined by the states with maximally-polarized angular momentum $\hat{\mathbf J}_b$, namely

$$\begin{align} \mathcal H^a = \big\{\vert \psi^{a}_{m_a}\rangle \big\},\, \mathrm{with}\,\, |\psi^{a}_{m_a}\rangle \equiv |J, m_a, m_b = J \rangle \,, \end{align} \tag{ 10 }$$

and $m_a+m_b\geqslant l^*$, which is highlighted in red in figure 1(c). This manifold is closed under dipole–dipole interactions governed by $H_{\mathrm{dd}}^{(ij)}$ as all processes leaving $\mathcal{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^{^{\prime}} \gt n$ manifolds. This additional ingredient will provide the kinetic term required for the simulation of QFTs.

The last manifold considered in our discussion is the vertical manifold defined by

$$\begin{align} \mathcal H^v = \big\{\vert \psi^{v}_{m_a}\rangle \big\}\,, \mathrm{with}\,\, |\psi^{v}_{m_a}\rangle \equiv |J, m_a, m_b = l^*-m_a \rangle , \end{align} \tag{ 11 }$$

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 $\mathcal H^a$ manifold.

The triangular manifold $\mathcal 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.

3.4.1. 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, ${\theta}_{ij} = \pi/2$. Furthermore, we focus on the regime $|{\omega}_{a(b)}| \gg J^2 V_{ij}$, where the interaction term ${\hat{J}}^{(i)}_{a,+}{\hat{J}}^{(j)}_{a,+}$ and ${\hat{J}}^{(i)}_{b,+}{\hat{J}}^{(j)}_{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

$$\begin{align} {\hat{H}}^t=\sum_i\,&\left({\hat{H}}^{(i)}_{\mathrm{MW}} + {\hat{H}}^{(i)}_\mathrm{P}\right)+\frac{1}{2}\sum_{i\neq j} V_{ij}\bigg[\left({\hat{J}}_{a,z}^{(i)}- {\hat{J}}_{b,z}^{(i)}\right)\left({\hat{J}}_{a,z}^{(j)}-{\hat{J}}_{b,z}^{(j)}\right)\nonumber\\&-\frac{1}{4} \sum_\sigma\, \left({\hat{J}}^{(i)}_{\sigma,+}{\hat{J}}^{(j)}_{\sigma,-} + \mathrm{H.c.} \right) + \frac{3}{2} \left(e^{i2{\phi}_{ij}} {\hat{J}}^{(i)}_{a,+} {\hat{J}}^{(j)}_{b,+} + \mathrm{H.c.}\right)\bigg], \end{align} \tag{ 12 }$$

where we transformed to the rotating frame defined below equation (6) under the resonant condition $\Delta_a + \Delta_b = {\omega}_a - {\omega}_b$. The term ${{\hat{H}}^{(i)}_{\mathrm{P}}}$ is the ponderomotive coupling for the ith atom defined as

$$\begin{align} {\hat{H}}_\mathrm{P} =\sum_{\kappa_a,\kappa_b} \lambda_{\kappa_a, \kappa_b} ({\hat{J}}_{a,+})^{\kappa_a} ({\hat{J}}_{b,+})^{\kappa_b} +\mathrm{H.c.}\,, \end{align} \tag{ 13 }$$

where $\lambda_{\kappa_a, \kappa_b}$ denotes complex coupling strengths. The ponderomotive coupling in equation (13) corresponds to a transfer of orbital angular momentum by $\kappa = \kappa_a + \kappa_b$, where $\kappa_{a(b)}$ can take on negative and positive values [thus with the identification $({\hat{J}}_{a(b),+})^{-1} \rightarrow {\hat{J}}_{a(b),-}$]. For the special case $\kappa_b = 0$ implying $\kappa = \kappa_a$, the ponderomotive Hamiltonian ${\hat{H}}_{\mathrm{P}}$ reduces to ${\lambda}_{\kappa}({\hat{J}}_{a,+})^{\kappa} + {\lambda}_{\kappa}^*({\hat{J}}_{a,-})^{\kappa}$, where we use the simplified notation ${\lambda}_{\kappa} \equiv \lambda_{\kappa,0}$. In this case, this coupling Hamiltonian can be used, for instance, to generate a two-axis twisting term [102], i.e. $i({\hat{J}}_{a,+})^2 -i({\hat{J}}_{a,-})^2$ for κ = 2.

3.4.2. Ponderomotive manipulation

In order to realize the nonlinear coupling in equation (13), we employ the ponderomotive manipulation techniques [62, 103–105]. The core idea is to interfere two co-propagating (along the z-direction) hollow Laguerre-Gauss laser beams with different frequencies ω₁ and ω₂ and non-zero orbital angular momentum. The corresponding interference pattern forms a ponderomotive potential $U_\kappa(\mathbf{r})e^{-i\kappa {\phi}}e^{-i\delta{\omega} t}+\mathrm{H.c.}\,$, rotating at the beating frequency $\delta{\omega} = {\omega}_1-{\omega}_2$, which can be arranged to be in the MW frequency domain to match the resonance condition $\delta{\omega} = -\kappa_a\, {\omega}_a + \kappa_b\, {\omega}_b$. When centered at the atomic position, this spatially-rotating ponderomotive potential transfers $\kappa = \kappa_a + \kappa_b$ quanta of orbital angular momentum, see figure 2(a).

Zoom In Zoom Out Reset image size

For simplicity, we focus here on the case $\kappa_b = 0$ and therefore $\kappa = \kappa_a$. The transition matrix elements of the ponderomotive coupling $U_{\kappa,m_a,m_b} \equiv \langle {\psi}^{t}_{m_a+\kappa,m_b}\vert U_\kappa(\mathbf{r})e^{-i\kappa {\phi}}\vert \psi^t_{m_a,m_b}\rangle$ are shown in figure 2(b) for $\kappa = 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 $U_{\kappa,m_a,m_b} = {\lambda}_{\kappa} \langle \psi^t_{m_a+\kappa,m_b}\vert({\hat{J}}_{a,+})^{\kappa} \vert \psi^t_{m_a,m_b}\rangle$. 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_{\kappa,m_a}$ scale as $n^{2\kappa}$ for $|m_a|\ll J$, and are in the range $2\pi\times(10^2-10^3)\,\mathrm{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 $|\kappa_a| + |\kappa_b|$ is increased. For $|\kappa_a| + |\kappa_b| \lesssim 6$ and a combined laser power of hundred mW per site give rise to coupling strength of hundreds of $2\pi\times$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.

3.4.3. 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 $\vert \psi^t_{J,J}\rangle$. We thus perform a Holstein–Primakoff transformation ${\hat{J}}_{-,\sigma}^{(i)} = \sqrt{2J}\sqrt{1-\hat n_{i\sigma}/(2J)}\,\hat c_{i\sigma}$ and ${\hat{J}}_{\sigma,z}^{(i)} = 2(J-\hat n_{i\sigma})$, where $\hat c^{}_{i\sigma}$ ($\hat c^\dagger_{i\sigma}$) are standard bosonic annihilation (creation) operators satisfying the commutation relations $[\hat c^{}_{i\sigma},\hat c^\dagger_{j\sigma^{^{\prime}}}] = \delta_{\sigma\sigma^{^{\prime}}}\delta_{ij}$ and $\hat n_{i\sigma} = \hat c^\dagger_{i\sigma} \hat c^{}_{i\sigma}$, with $\sigma = a,b$ playing the role of pseudo-spin.

In this framework, the circular level serves as the vacuum for both bosonic modes, $\vert 0,0\rangle\equiv\vert \psi^t_{J,J}\rangle$, and the Hamiltonian (without MW and ponderomotive couplings) takes the form of a generalized spinful Bose–Hubbard model

$$\begin{align} {\hat{H}}^t_{\mathrm{HP}} = \sum_{\sigma=a,b}\Bigg[-&\sum_i\,\Delta_\sigma\,\hat n_{i\sigma} - \frac{1}{2}\sum_{i\neq j} h_{ij}\, \bigg(\hat c_{i\sigma}^\dagger \hat c_{j\sigma}^{}+\mathrm{H.c.}\bigg)\Bigg]\\ \quad + \frac{1}{2}&\sum_{i\neq j} \left[\left(w_{ij}\hat c^\dagger_{ia}\hat c^\dagger_{jb} +\mathrm{H.c.}\right)+U_{ij}(\hat n_{ia}-\hat n_{ib})(\hat n_{ja}-\hat n_{jb})\right],\nonumber \end{align} \tag{ 14 }$$

which is valid up to corrections $O(1/J)$. Here $\Delta_\sigma$ are onsite energies, $h_{ij} = JV_{ij}/2$ describes intra-species hopping, $w_{ij} = 3JV_{ij}{\textrm{exp}}[i2{\phi}_{ij}]$ pair gain/loss processes, and $U_{ij} = 2V_{ij}$ 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_{ij}|\lesssim-\Delta_\sigma$ the low energy dynamics (with respect to $-J\Delta_\sigma$) 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 $\vert \psi^t_{J,J}\rangle$, which can be initialized through well established circularization methods, i.e. via the adiabatic rapid passage method [2, 107–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].

The edge manifold $\mathcal H^a$ is defined in equation (10) by having ${\hat{\mathbf{J}}}_b$ maximally polarized, as highlighted in figure 1(c). Notice that we could define an analogous manifold $\mathcal H^b$ by considering the states with ${\hat{\mathbf{J}}}_a$ maximally polarized. Choosing the static external electric and magnetic fields such that $|{\omega}_{a(b)}|,$ $|{\omega}_a \pm{\omega}_b| \gg J^2 V_{ij}$, this manifold is closed under dipole–dipole interaction induced processes governed by ${\hat{H}}^{ij}_\mathrm{dd}$ [112]. Since $\mathcal H^a$ is a subset of $\mathcal H^t$, it inherits many of the properties discussed before, such as the ponderomotive control. Additionally, we can also construct a nonlinear squeezing term

$$\begin{equation} {\hat{H}}_{\mathrm{SQ}} = \chi ({\hat{J}}_{a,z})^2\,, \end{equation} \tag{ 15 }$$

where χ is the squeezing strength. This can be generated by off-resonantly coupling the $\mathcal H^a$ manifold with principal quantum number n to another $\mathcal H^a$ manifold with $n^{^{\prime}}\gt 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 $\Delta_{m_a} = \Delta_0 + \delta{\omega}_a m_a$ of each two-level system varies with m_a , where Δ₀ is an overall detuning and $\delta{\omega}_a = {\omega}_a(n^{^{\prime}}) - {\omega}_a(n)$ is the differential level shift. Moreover, the coupling strength ${\Omega}_{m_a}$ depends smoothly on m_a through the dipole transition matrix elements. In the regime $\Delta_{m_a} \gg {\Omega}_{m_a}$ the coupling leads to an AC Stark shift ${\Omega}_{m_a}^2/(4\Delta_{m_a})$ that can be expanded in a power series of m_a when $\Delta_0\gg \delta{\omega}_a m_a$, namely ${\Omega}_{m_a}^2/(4\Delta_{m_a}) = \chi m_a^2+ O(m_a^4)$, where the constant and linear terms are absorbed in Δ₀ and $\delta{\omega}_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 Δ₀ and ${\Omega}_{m_a}$, as discussed in appendix E.

For moderate values of the MW intensity, χ can be on the order of $2\pi\times 100\,\mathrm{kHz}$, thus giving rise to level shifts as large as $2\pi\times 10\,\mathrm{MHz}$ for $|m_a|\sim 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 ${\omega}_a \gg {\Omega}_{m_a}^2/\Delta_{m_a}$.

By taking into account all the terms discussed so far, the Hamiltonian for the $\mathcal{H}^a$ manifold takes the form

$$\begin{align} {\hat{H}}^a= &\sum_i\, \left\{-\Delta_a {\hat{J}}_{a,z}^{(i)} + \chi \big({\hat{J}}_{a,z}^{(i)}\big)^2 + \sum_\kappa \left[\lambda_{\kappa} \big({\hat{J}}_{a,+}^{(i)}\big)^{\kappa} + \mathrm{H.c.}\right] \right\} \nonumber\\ \quad&+\frac{1}{2}\sum_{i \neq j}V_{ij}\left[ {\hat{J}}^{(i)}_{a,z} {\hat{J}}^{(j)}_{a,z}-\frac{1}{4}\left({\hat{J}}^{(i)}_{a,+} {\hat{J}}^{(j)}_{a,-}+\mathrm{H.c.}\right)\right]\,. \end{align} \tag{ 16 }$$

The first line describes linear and nonlinear single-particle terms. As previously mentioned, the κ = 1 term can also be obtained by MW engineering ($\lambda_1 = {\Omega}_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 (${\theta}_{ij} = \pi/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 $\mathcal H^a$ manifold require specific schemes involving the states in $\mathcal H^v$, which we discuss in detail in the next subsection.

The vertical manifold $\mathcal 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 $\mathcal H^v$. Second, a mapping scheme that we demonstrate below allows to transfer arbitrary states from $\mathcal H^v$ to $\mathcal H^a$ and vice versa with high fidelity. As a consequence, the manifold $\mathcal H^v$ provides the capabilities to perform state preparation and readout on $\mathcal H^a$, which make the Hamiltonian in equation (16) a particularly accessible many-body model.

The presence of low orbital angular momentum components in $\vert \psi^v_{m_a}\rangle$, see figure 3(a), offers a unique opportunity for laser coupling the atomic ground state to the Rydberg manifold $\mathcal{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 $\mathcal{H}^v$ manifold of erbium in turn contains an admixture of Rydberg g-states, which allows direct laser coupling from the ground state $|\mathrm{gs}\rangle$ to specific states in $\mathcal{H}^v$, selected by the laser frequency, see figure 3(b). For multiple laser frequencies, this coupling is described in the rotating frame by the Hamiltonian

$$\begin{align} {\hat{H}}^{\,vg} = \sum_{m_a} \left( \frac{{\Omega}_{m_a, \mathrm{gs}}}{2}\, |\psi_{m_a}^v\rangle \langle \mathrm{gs}| + \mathrm{H.c.} \right) \,, \end{align} \tag{ 17 }$$

where ${\Omega}_{m_a, \mathrm{gs}}$ are the individual Rabi frequencies that satisfy $|{\Omega}_{m_a, \mathrm{gs}}|\ll |{\omega}_S|$. This laser adressability can also be used to perform projective readout within the $\mathcal{H}^v$ manifold employing quantum gas microscope techniques [113, 114]. Note that the connectivity offered by the individual ${\Omega}_{m_a, \mathrm{gs}}$ in equation (17) provides a convenient setting to perform holonomic quantum computing operations within the $\mathcal{H}^v$ manifold [115].

Zoom In Zoom Out Reset image size

We now proceed to discuss a scheme realizing the state transfer

$$\begin{align} \hat U\,\,^{\,va}:\,\, |\psi_{m_a}^v\rangle \longleftrightarrow |\psi_{m_a}^a \rangle \quad \forall m_a \,, \end{align} \tag{ 18 }$$

via a adiabatic rapid passage that generalizes the one originally developed to prepare circular Rydberg levels [2]. We start by noting that a state $\vert \psi^v_{m_a}\rangle$ is connected to the final state $\vert \psi^a_{m_a}\rangle$ 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 ${\hat{H}}_{\mathrm{MW}}(t)$ (equation (6)) with ${\Omega}_a = 0$, such that only couplings among the highlighted states in figure 3(b) are present. During the driving sweep, the detuning $\Delta_b(t)$ has to change sign and the Rabi coupling ${\Omega}_b(t)$ has to be turned on and off again in order to adiabatically ($|{\Omega}_b(t)|\ll |{\omega}_{a(b)}|$) connect the different states $\vert \psi^v_{m_a}\rangle$ and $\vert \psi^a_{m_a}\rangle$. Note that in general the individual states acquire a state-dependent dynamical phase that can be canceled by an appropriate transformation in $\mathcal H^v$.

In the following, we numerically demonstrate the adiabatic rapid passage $\hat U\,\,^{va}$ for a selected state ($m_a = 0$) and a particular driving protocol $\Delta_b(t) = \Delta_0{\textrm{cos}}(\pi\,t/T)$ and ${\Omega}_b(t) = {\Omega}_0{\textrm{sin}}(\pi\,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 $\vert \psi^v_{m_a}\rangle$ changes into the target final state $\vert \psi^a_{m_a}\rangle$, 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 \lt -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|\lt 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 $\mathcal{H}^v$ to their respective counterpart in $\mathcal{H}^a$. For Δ₀ and Ω₀ on the order of $2\pi\times 10\,\mathrm{MHz}$ and T on the order of a few $\mu\mathrm{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.

Standard quantum-mechanical position $\hat{\varphi}$ and momentum $\hat{\pi}$ operators obey the canonical commutation relation $\left[\hat{\varphi}, \hat{\pi}\right] = i$. For quantum mechanics on a unit circle, it will be more convenient to consider phase operators $e^{i\hat{\varphi}}$ satisfying

$$\begin{align} \left[\hat{\pi},e^{i\hat{\varphi}} \right] = e^{i\hat{\varphi}} \;, \end{align} \tag{ 19 }$$

where the momentum operator $\hat \pi$ takes discrete eigenvalues $m_\pi \in \mathbb{Z}$ on its eigenstates $|m_\pi\rangle$.

The commutation relation (19) shows that $e^{i\hat{\varphi}}$, similar to a spin operator ${\hat{J}}_+$, acts as a raising operator on $|m_\pi\rangle$. Taking into account the normalization of the operators, this suggests to identify

$$\begin{align} \hat{\pi} \leftrightarrow {\hat{J}}_z \;, \qquad e^{\pm i\hat{\varphi}} \leftrightarrow \frac{1}{\sqrt{J(J+1)}}{\hat{J}}_\pm \;. \end{align} \tag{ 20 }$$

The $(2J+1)$-dimensional Hilbert space spanned by eigenstates $|m\rangle$ of ${\hat{J}}_z$ with $m = -J, \dots J$ thus 'regularizes' the infinite-dimensional Hilbert space of the continuous variables by cutting off the spectrum of $\hat{\pi}$ at a maximum value $|m_\pi| \leqslant J$ [120]. The identification becomes exact for states around m = 0 in the limit $J\, \gg 1$.

We now numerically demonstrate how the above identification converges to a continuous variable representation in the large spin limit $J\gg 1$. Here, we focus on a single Rydberg atom in the subset of states $\mathcal{H}^a$ described by the Hamiltonian in equation (16). Choosing $\Delta_a = 0$ and defining ${\Omega}^{^{\prime}} = -2{\Omega}_a \sqrt{J(J+1)}$ and $\lambda^{^{\prime}} = -2 {\lambda}_{\kappa} e^{i {\theta}} [J(J+1)]^{\kappa/2}$ for a fixed κ > 1, $\lambda^{^{\prime}}\gt0$ and a phase θ, the single Rydberg atom Hamiltonian in the continuous variable limit becomes

$$\begin{equation} {\hat{H}}_{J\rightarrow\infty} = \chi \, \hat \pi^2 - {\Omega}^{^{\prime}} {\textrm{cos}} {\hat{\varphi}} - \lambda^{^{\prime}} {\textrm{cos}}(\kappa{\hat{\varphi}}+{\theta}) \;. \end{equation} \tag{ 21 }$$

The Hamiltonian in equation (21) describes a particle of mass $(2\chi)^{-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 ${\hat{\varphi}}$ and solving the eigenvalue problem in Fourier space, thus yielding the result shown in figure 4(b) for $\lambda^{^{\prime}} = 0$. Depending on the ratio ${\Omega}^{^{\prime}}/\chi$, the low-energy physics of this model ranges from the free particle regime (${\Omega}^{^{\prime}} \ll \chi$) to the harmonic oscillator regime (${\Omega}^{^{\prime}} \gg \chi$).

Zoom In Zoom Out Reset image size

We now consider an intermediate regime, where the nonlinear character of ${\textrm{cos}}{\hat{\varphi}}$ 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 \gtrsim 30$. In figure 4(d), we further test the effect of $\lambda^{^{\prime}}$ 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 $\lambda^{^{\prime}}$.

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.

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]

$$\begin{align} {\hat{H}}_\mathrm{SG} = \int \mathrm{d} x \, \bigg\{\frac{1}{2}[\hat\pi(x)]^2 + \frac{1}{2}[\partial_x{\hat{\varphi}}(x)]^2 -\frac{M_0^2}{\beta^2}{\textrm{cos}}[\beta {\hat{\varphi}}(x)]\bigg\} \;, \end{align} \tag{ 22 }$$

where the fields obey $[{\hat{\varphi}}(x),\,\hat\pi(x^{^{\prime}})] = i\delta(x-x^{^{\prime}})$, M₀ 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 $\beta^2 = 8\pi$, 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 \lt \beta^2 \lt 8\pi$), the low-energy degrees of freedoms are fermions interacting attractively (repulsively) for $\beta^2 \lt 4\pi$ ($\beta^2 \gt 4\pi$), which can be identified with quantized solitons of mass $M \overset{\beta\rightarrow 0}{\longrightarrow} 8 M_0/\beta$. In the attractive regime, the solitons form bound states, so-called breathers. The lightest one of them has a mass $m_1\overset{\beta\rightarrow 0}{\longrightarrow} 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 $\mathcal 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

$$\begin{align} {\hat{H}}_{\text{SG}}^{(\text{lat})} = &\sum_i \left[\chi \hat{\pi}_i^2 - \lambda^{^{\prime}} {\textrm{cos}}(\kappa\hat{\varphi}_i) - \sum_{j> i}\frac{V_{ij}^{^{\prime}}}{2}{\textrm{cos}} \left(\hat{\varphi}_{j}-{\hat{\varphi}}_i\right) \right] \,, \end{align} \tag{ 23 }$$

where we rescaled the interaction $V_{ij}^{^{\prime}} = V_{ij} J(J+1)$, while keeping $\Delta_a = 0$, and $\lambda^{^{\prime}} = -2{\lambda}_{\kappa} [J(J+1)]^{\kappa/2}$ for a single value of κ, which corresponds to either a MW (κ = 1) or ponderomotive ($\kappa \geqslant 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^{^{\prime}}_{i,i+1} \equiv V^{^{\prime}}_\text{nn}$ and rescale the continuous variables and Hamiltonian appropriately (see appendix G). The correspondence is achieved by the identification of parameters as follows

$$\begin{align} \frac{V^{^{\prime}}_\text{nn}}{\chi} = \frac{4\kappa^4}{\beta^4}\;, \qquad \frac{\lambda^{^{\prime}}}{\chi} = \frac{2\kappa^2(\ell M_0)^2}{\beta^4} \;, \end{align} \tag{ 24 }$$

where we introduced a short-distance scale $\ell$ that sets a UV-cutoff $\Lambda \propto 1/\ell$ for the lattice regularization, $\ell M_0 = \kappa \sqrt{2\lambda^{^{\prime}}/V^{^{\prime}}_\text{nn}} \rightarrow 0$. With these identifications, all the regimes of the theory, $0 \lt \beta^2 = 2\kappa^2 \sqrt{\chi/V^{^{\prime}}_\text{nn}} \lt \infty$, are in principle accessible. In particular, in appendix G we present realistic experimental parameter estimates, which indicate that values $\beta^2 \gtrsim 8\pi$ 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.

4.2.2. 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^{^{\prime}}_\text{nn} = 10 \chi$ and κ = 1, i.e. $\beta^2 = \sqrt{0.4} \ll 8\pi$, and vary $\lambda^{^{\prime}}/\chi$, thereby effectively scanning different values of $\ell M_0$. In this parameter regime, we expect the model to be well approximated by a quadratic (free) theory obtained by replacing ${\textrm{cos}} {\hat{\varphi}}_i \approx 1 - {\hat{\varphi}}_i^2/2$, and similarly for ${\textrm{cos}} \left(\hat{\varphi}_{i+1} - {\hat{\varphi}}_i\right)$. After diagonalizing the quadratic theory (see, e.g. [128] and appendix H), we obtain the single-particle dispersion relation ${\omega}_q = \sqrt{2\chi(\lambda^{^{\prime}}+V^{^{\prime}}_{\text{nn}}-V^{^{\prime}}_{\text{nn}}{\textrm{cos}} q)}$ with $q \in [0,2\pi]$, where we 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 $\langle {\hat{\varphi}}_i \rangle = 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 ${\omega}_q \approx \sqrt{2\chi[\lambda^{^{\prime}}- V_\text{dd}^{^{\prime}}(c_2^{} q^2 + c_2^{^{\prime}}\, q^2 {\mathrm{log}} q)]}$, with $c_2 \approx -0.739$ and $c_2^{^{\prime}}\approx 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 $\Delta E = {\omega}_0 = \sqrt{2\chi\lambda^{^{\prime}}}$ is thus unaffected by the long-range tails.

Zoom In Zoom Out Reset image size

We now study the gap $\Delta E$ of the model in more detail. As previously stated, we expect the Ising terms' contribution to vanish in the $J\gg 1$ limit. This is confirmed in figure 5(b), where we show the dependence of the energy gap $\Delta E$ in the finite-spin model as a function of J with $J = 5, \dots, 12$, both with and without Ising terms. Comparing against the predicted result ${\omega}_0 = \sqrt{2\chi\lambda^{^{\prime}}}$ 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 ${\textrm{cos}}{\hat{\varphi}}_i \rightarrow {\textrm{cos}}({\hat{\varphi}}_i+{\alpha})$, which can be obtained in the spin model by taking, for example, a MW field in equation (6) with complex Rabi frequency ${\Omega}_a \rightarrow {\Omega}_a e^{i{\alpha}}$. For small α, this corresponds to a small uniform displacement of the spin expectation value away from $\langle {\hat{J}}_y^{(i)} \rangle = 0$. Quenching to α = 0, the energy gap can be extracted by measuring the oscillation frequency of $\langle {\hat{J}}_y^{(i)} \rangle (t)$.

The leading contribution of the induced dynamics within linear-response, i.e. for $|{\alpha}| \ll \pi$, is $\langle {\hat{J}}_y^{(i)} \rangle \propto \langle {\hat{\varphi}}_i \rangle + O\left(\langle {\hat{\varphi}}_i^3 \rangle \right)$, 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 $\lambda^{^{\prime}}/\chi$ values, as shown in figure 5(c). The initial state preparation can be performed using the scheme outlined in section 3.6 for large $\lambda^{^{\prime}}$ followed by an adiabatic ramp to the final value of $\lambda^{^{\prime}}/V^{^{\prime}}_\mathrm{nn}$ and $\lambda^{^{\prime}}/\chi$, which is protected by the many-body gap. The same scheme in section 3.6 also allows to measure the value of $\langle {\hat{J}}_y^{(i)} \rangle$ 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 $\langle {\hat{J}}_+^{(i)} \rangle \propto \langle e^{i{\hat{\varphi}}_i}\rangle = \langle {\textrm{cos}}{\hat{\varphi}}_i\rangle$ and compare it against the free theory prediction $\langle {\textrm{cos}}{\hat{\varphi}}_i\rangle = {\textrm{exp}}\left[ -\frac{\chi}{2N}\sum_q \frac{1}{{\omega}_q} \right]$. The plot displays a systematic difference between the numerical ED result and the free theory prediction. As shown in the inset for $\lambda^{^{\prime}} = 50\chi$, 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 $\langle {\hat{J}}_+^{(i)} {\hat{J}}_-^{(j)} \rangle$, 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–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].

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

$$\begin{align} {\hat{H}}_{\text{mSG}}= \int \mathrm{d}x \, \bigg\{\frac{1}{2} [\hat \pi(x)]^2 + \frac{1}{2} [\partial_x {\hat{\varphi}}(x)]^2 + \frac{M^2}{2} [{\hat{\varphi}}(x)]^2 - u {\textrm{cos}} [\beta {\hat{\varphi}}(x) + {\theta}]\bigg\} \;, \end{align} \tag{ 25 }$$

provides a dual low-energy effective description of the massive Schwinger model [119, 121, 136], described by

$$\begin{align} {\hat{H}}_\text{QED} = \int \mathrm{d}x \, \bigg\{\frac{1}{2} &\left[\hat E(x) + \frac{e{\theta}}{\pi}\right]^2 + \hat {{\psi}}^\dagger(x) \gamma^0 \left[-i \gamma^1 \hat D_x + m \right] \hat \psi(x) \bigg\} \;. \end{align} \tag{ 26 }$$

Here, $\gamma^{0(1)}$ are the gamma matrices in two space-time dimensions, we abbreviated $\hat D_x = \partial_x - ie\hat A(x)$, the fermionic spinors obey $\left\{\hat\psi_j(x), \hat \psi^\dagger_{j^{^{\prime}}} (x^{^{\prime}})\right\} =$ $\delta_{jj^{^{\prime}}} \delta(x-x^{^{\prime}})$ and the gauge fields fulfill $\left[\hat A(x), \hat E(x^{^{\prime}})\right] = i \delta(x-x^{^{\prime}})$. The duality requires to fix the mass M, the coupling u and the parameter β as

$$\begin{align} \beta = 2\sqrt{\pi}, \qquad M^2 = \frac{e^2}{\pi}, \qquad u = \frac{{\textrm{exp}}({\gamma})}{2\pi} \Lambda m, \end{align} \tag{ 27 }$$

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)) $\Lambda \gg e, m$. In this limit, one obtains a QFT described by two dimensionless parameters, the coupling $0 \leqslant e/m \leqslant \infty$ and the topological angle ${\theta} \in \left[-\pi,\pi\right]$.

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\gt0$, the gauge fields mediate a density–density interaction $\propto \hat \rho (x) \hat \rho (x^{^{\prime}})$ among the charge densities $\hat \rho = e \hat \psi^\dagger \hat \psi$ with a Coulomb potential growing with distance as $\propto |x-x^{^{\prime}}|$. 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 ${\hat{\varphi}}$, which can be identified with the electric field $\hat E / e$. In the limit $e/m = \infty$, 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 $\pm e {\theta}/ \pi$ via the Gauss law constraint $\partial_x \hat E = \hat \rho$. The physics is periodic in θ with period 2π (which is manifest in ${\hat{H}}_{\mathrm{mSG}}$) because any background charge $Q = \pm ne$ with integer n is screened by the production of particle pairs induced by the corresponding background electric field, $\langle \hat E \rangle = Q$. For the special case of ${\theta} = \pm \pi$ (without additional external charges), the ground state of the model undergoes a second-order quantum phase transition at $(m/e)_c \approx 0.3335(2)$ [137], with order parameter $\langle \hat E \rangle$ (or $\langle {\hat{\varphi}} \rangle$), which falls in the Ising universality class and corresponds to the spontaneous breaking of the $\hat E \rightarrow - \hat E$ (or ${\hat{\varphi}} \rightarrow - {\hat{\varphi}}$) $\mathbb Z_2$ symmetry. For fixed $m/e \gt (m/e)_c$, the system undergoes a first-order phase transition at ${\theta} = \pm\pi$, see figure 6(a) for a sketch of the phase diagram. Further details can be found in [119, 136].

Zoom In Zoom Out Reset image size

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 reads

$$\begin{align} {\hat{H}}_{\text{mSG}}^{(\mathrm{lat})} = \sum_i \left[\chi \hat{\pi}_i^2 - \lambda^{^{\prime}}_\kappa {\textrm{cos}}(\kappa\hat{\varphi}_i+{\theta}) -{\Omega}^{^{\prime}}{\textrm{cos}}(\hat{\varphi}_i) -\sum_{j> i}\frac{V_{ij}^{^{\prime}}}{2}{\textrm{cos}} \left(\hat{\varphi}_{j}-{\hat{\varphi}}_i\right) \right] \;, \end{align} \tag{ 28 }$$

where we use the identification $V_{ij}^{^{\prime}} = V_{ij} J(J+1)$, ${\Omega}^{^{\prime}} = -{\Omega}_a [2\sqrt{J(J+1)}]$ and $\lambda^{^{\prime}}_\kappa = -2{\lambda}_{\kappa} [J(J+1)]^{\kappa/2}$, and neglect the Ising interactions. Our target model is again approached in the continuum limit with the identifications

$$\begin{align} \frac{V^{^{\prime}}_\mathrm{nn}}{\chi} = \frac{\kappa^4}{(2\pi)^2} \,, \qquad \frac{\lambda^{^{\prime}}_\kappa}{\chi} = \frac{\kappa^2\ell^2{\textrm{exp}}(\gamma)\Lambda m}{(2\pi)^2}\,, \qquad \frac{{\Omega}^{^{\prime}}}{\chi} = \frac{\kappa^4e^2 \ell^2}{(2\pi)^3}\,, \end{align} \tag{ 29 }$$

and keeping only nearest-neighbor terms $V^{^{\prime}}_\mathrm{nn}\equiv V^{^{\prime}}_{i,i+1}$. As discussed in appendix I, the massive SG model in equation (25) is recovered under the conditions

$$\begin{align} \chi,\, {\lambda}_{\kappa}^{^{\prime}} \ll {\Omega}^{^{\prime}} \ll V_\mathrm{nn}^{^{\prime}} \;, \qquad \lambda^{^{\prime}} \ll \sqrt{\chi V_\mathrm{nn}^{^{\prime}}}\,. \end{align} \tag{ 30 }$$

We emphasize that a suitable choice of parameters in principle allows us to probe the whole non-trivial range of the theory with $0 \lt e/m \lt \infty$. A specific realistic example of parameters is provided in appendix I, where we employ a spin length of J = 30, yielding $e/m \approx 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 $\mathcal H^a$ manifold.

4.3.2. Numerical illustration: capacitor discharge

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 $\vert i\rangle$ with energy ω_i to a final state $\vert f\rangle$ with energy ω_f is determined by the spontaneous emission rates $\gamma_{i, f} = {\omega}_{i,f}^3|\langle i\vert\hat{\boldsymbol{\mu}}\vert f\rangle|^2/(3\pi{\epsilon}_0 c^3)$ [153], where ${\omega}_{i,f} = {\omega}_i- {\omega}_f$ is the transition frequency and $\langle i\vert\hat{\boldsymbol{\mu}}\vert f\rangle$ the dipole transition matrix element. The total lifetime τ_i of $\vert i\rangle$ is then given by $\tau^{-1}_{i} = \sum_{f\lt i} \gamma_{i, f}$, where the $\sum_{f\lt i}$ is limited to final states with ${\omega}_f \lt {\omega}_i$.

In figure 8(a) we present the numerically calculated lifetimes for the manifold $\mathcal{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⁵ [61] and it gradually changes to the numerically calculated n⁴ scaling for states with low orbital angular momentum components. The n⁴ 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³ [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.

Zoom In Zoom Out Reset image size

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–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–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_{ij}J(J+1)$, as considered in section 4 of the main text, gives rise to an additional potential $U( R_i) = 3R_iV_{ij}J(J+1)/|R_{ij}|$ 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 = \sqrt{1/(M_a{\omega}_t)}$ (M_a is the mass of the atom) is small compared to ω_t , i.e. $U(l_t)\ll {\omega}_t$, then at short times (${\omega}_t t/(2\pi)\ll 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_{ij}J(J+1) = 2\pi\times 300\,\mathrm{kHz}$ for n = 41 and an inter-atomic separation of $R_{ij} = 17\,\mu\mathrm{m}$. Assuming a typical trapping frequency of ${\omega}_t = 2\pi\times40\,\mathrm{kHz}$ and strontium atoms ($M_a \approx 87 \mathrm{u}$) gives rise to $l_t = 54\,\mathrm{nm}$ for which we get $U(l_t) \approx 2\pi\times 3 \, \mathrm{kHz}\ll{\omega}_t$. For these parameter values, $n(10t_c) \approx 0.05$ excitations are created after ten interaction cycles of duration $t_c = 2\pi/[V_{ij}J(J+1)]$.

The replacement of the dipole operator by the angular momentum operators ${\hat{\mathbf{J}}}_{a}$ and ${\hat{\mathbf{J}}}_{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^{^{\prime}} \neq n$ manifolds in the presence of a static external electric field. For electric fields below the Ingris Teller limit ($F\lt F_\mathrm{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 ${\hat{J}}_{a,+}^{(i)}{\hat{J}}_{a,-}^{(j)}+{\hat{J}}_{a,-}^{(i)}{\hat{J}}_{a,+}^{(j)}$, relevant for the realization of the Sine-Gordon and massive Schwinger model, the relative mismatch due to manifold mixing is 10⁻³ for n = 41.

Zoom In Zoom Out Reset image size

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^{^{\prime}}, n^{^{\prime\prime}}\neq 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_{ij}\gtrsim 10\,\mu\mathrm 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 $\mathcal{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).

The vector potential of a linearly polarized LG laser beam with orbital angular momentum $\delta m$ and mode number p, propagating along the z-axis, is in the Lorenz gauge given by [157, 158]

$$\begin{align} {\mathbf{A}}_{p,\delta m}({\mathbf{r}},\, t) =&\, \hat{\mathbf{e}}_x\, A_{p,\delta m} \, u_{p,\delta m}({\mathbf{r}})e^{i{\alpha}} \,e^{-i\delta m{\phi}}\,e^{-i(kz-{\omega} t)}/2 + \mathrm{c.c.}, \end{align} \tag{ D1 }$$

where we used cylindrical coordinates ($\rho = \sqrt{x^2+y^2}$ and $\tan {\phi} = y/x$). Here ω is the light frequency, $k = 2\pi/\lambda$ is the magnitude of the corresponding wave vector, λ is the wavelength and α is a phase. Moreover, $A_{p,\delta m}$ denotes the field amplitude and $u_{p,\delta m}({\mathbf{r}})$ is a dimensionless field distribution function satisfying the Helmholtz equation. In the paraxial limit $u_{p,\delta m}({\mathbf{r}})$ can be expressed in terms of associated Laguerre polynomials $L_p^{(|\delta m|)}(x)$ as

$$\begin{align} u_{p,\delta m}(\rho, z) & = C_{p,\delta m}\,\frac{w_0}{w(z)}\,\bigg(\frac{\sqrt{2}\rho}{w(z)} \bigg)^{|\delta m|} \,L_p^{(|\delta m|)} \bigg(\frac{2\rho^2}{w(z)^2} \bigg) {\textrm{exp}}\bigg(\frac{-\rho^2}{w(z)^2} \bigg)\nonumber\\ & \quad \times {\textrm{exp}} \left(i\left[\frac{k\rho^2z}{2(z^2+z_R^2)}-(2p+|\delta m|+1)\mathrm{atan}\frac{z}{z_R} \right] \right). \end{align} \tag{ D2 }$$

The beam waist is defined as $w(z) = w_0\sqrt{1+z^2/z^2_R}$, where w₀ denotes the waist in the focal plane and $z_R = \pi w_0^2/\lambda$ is the Rayleigh length. The normalization constant $C_{p,\delta m} = \sqrt{p!\,2/(p+|\delta m|)!}$ is chosen such that $\int\,\mathrm{d\,}A |u(\rho,z)|^2 = w_0^2 \pi$. Due to this normalization the time averaged electric beam power $P = c{\epsilon}_0{\omega}^2|A_{p,\delta m}|^2w_0^2\pi/4$ only depends on the field amplitude and not on the LG mode numbers.

The ponderomotive potential that a highly excited Rydberg electron experiences from a fast oscillating laser field is given by [159] $U_\mathrm{PM}({\mathbf{R}}+\mathbf{r},t) = e^2|{\mathbf{A}}(\mathbf{R}+\mathbf{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₀, R can be replaced by its mean value. Furthermore, if the laser beam is focused onto the center of the trap, i.e. $\langle \mathbf{R}\rangle = 0$, which can be achieved by re-using the trapping light optics, the ponderomotive potential becomes approximately independent of the core position $U_\mathrm{PM}({\mathbf{R}}+\mathbf{r},t)\approx U_\mathrm{PM}(\mathbf{r},t)$.

As already mentioned in the main text, the core idea to couple states with multiple orbital angular momentum $\kappa\geqslant 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–105]. Using equation (D1), the vector potential of the two co-propagating LG beams is given by

$$\begin{align} {\mathbf{A}}_\mathrm{LG} ({\mathbf{r}},\, t) =\, \hat{\mathbf{e}}_x\, A_\mathrm{LG} \big[u_{\delta m_1}(\rho,z)e^{-i\delta m_1 {\phi}}e^{-i({\omega}_1 t - k_1 z-{\alpha}_1)} + u_{\delta m_2}(\rho,z)e^{i\delta m_2 {\phi}}e^{-i({\omega}_2 t-k_2 z-{\alpha}_2)} + \mathrm{c.c.}\big]/2, \end{align} \tag{ D3 }$$

where the subscripts 1 and 2 are labels for the first and second LG beam. The field distribution functions $u_{\delta m_{1(2)}}$ are in the most general case a superposition of multiple LG modes $u_{p,\delta m_{1(2)}}$ from equation (D2), with fixed $\delta 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

$$\begin{align} U_\mathrm{PM}(\mathbf{r},t) \approx \big[ U_c({\mathbf{r}}) + (U_\kappa({\mathbf{r}})e^{-i\kappa{\phi}}e^{-i\delta{\omega} t} + \mathrm{c.c.})\big], \end{align} \tag{ D4 }$$

where $\delta{\omega} = {\omega}_1 - {\omega}_2$ is the oscillation frequency of the interference term and $\kappa = \delta m_1 + \delta m_2$ is the transferred orbital AM. The spatial shapes of the constant and oscillatory terms are given by

$$\begin{align} U_c({\mathbf{r}}) &= U_0\,\big[|u_{\delta m_1}(\rho,z)|^2 + |u_{\delta m_2}(\rho,z)|^2\big]/2\nonumber\\ U_\kappa({\mathbf{r}}) &= U_0e^{i({\alpha}_1-{\alpha}_2)}\,u_{\delta m_1}(\rho,z)\,u^*_{\delta m_2}(\rho,z)/4, \end{align} \tag{ D5 }$$

with $U_0 = e^2P_\mathrm{tot}/(16\pi^3m_e{\epsilon}_0c^3)\,(\lambda^2/w_0^2)$, where $P_\mathrm{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

$$\begin{align} U^{\kappa_a,\kappa_b}_{m_a ,m_b}=\langle \psi^t_{m_a+\kappa_a, m_b+\kappa_b}\vert U_\kappa(\rho,z) e^{-i\kappa{\phi}} \vert \psi^t_{m_a,m_b}\rangle\;, \end{align} \tag{ D6 }$$

that satisfy orbital angular momentum conservation $\kappa_a+\kappa_b = \kappa$. Note, $\kappa_{a(b)}$ can take on negative and positive values. For the determination of the transition matrix element the functional behavior of $U_\kappa(\rho,z)$ is of crucial importance. In the limit where the beam waist w₀ 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_\kappa(\rho,z = 0)\approx C_{\kappa+2\eta}\,\rho^{|\kappa|+2\eta}$, where the leading order $|\kappa| + 2\eta$, with $\eta\geqslant 0$, is determined by the LG beam modes.

Zoom In Zoom Out Reset image size

We now claim that, within the manifold $\mathcal{H}^t$, the matrix elements of the ponderomotive potential near the phase vortex $c_{\kappa+2\eta}\,\rho^{|\kappa|+2\eta}e^{-i\kappa{\phi}}$ are exactly proportional to those of $({\hat{J}}_{a,+})^{\kappa_a} ({\hat{J}}_{b,+})^{\kappa_b}$, with $|\kappa_a| + |\kappa_b| = \kappa+ 2\eta$ and where we use the convention $(J_{a(b),+})^{-\kappa_{a(b)}} = (J_{a(b),-})^{\kappa_{a(b)}}$. More specifically,

$$\begin{align} W^{\kappa_a,\kappa_b}_{m_a,m_b} &\propto n^{|\kappa|+2\eta} M^{\kappa_a,\kappa_b}_{m_a,m_b}, \end{align} \tag{ D7 }$$

with

$$\begin{align*} W^{\kappa_a,\kappa_b}_{m_a,m_b} &= \langle \psi^t_{m_a+\kappa_a,m_b+\kappa_b}\vert(\rho/a_0)^{|\kappa|+2\eta} e^{-i\kappa{\phi}} \vert \psi^t_{m_a,m_b}\rangle, \end{align*}$$

and

$$\begin{align*} M^{\kappa_a,\kappa_b}_{m_a,m_b} &= \langle \psi^t_{m_a+\kappa_a,m_b+\kappa_b}\vert({\hat{J}}_{a,+})^{\kappa_a} ({\hat{J}}_{b,+})^{\kappa_b} \vert \psi^t_{m_a,m_b}\rangle. \end{align*}$$

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 $\kappa_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².

To analyze the ponderomotive transition matrix elements $U^{\kappa_a,\kappa_b}_{m_a ,m_b}$ from equation (D6) it is important to take into account the full shape of $U_\kappa(\rho,z)$. As the functional behavior of $U_\kappa(\rho,z)$ is determined by the LG beam modes, it is sufficient to analyze $u_{\delta m_1}(\rho,0)\,u^*_{\delta m_2}(\rho,0)$. A series expansion of the field distribution functions at z = 0 and around ρ = 0 yields in the most general case $u_{\delta m_1}(\rho,0)\,u^*_{\delta m_2}(\rho,0) =$ $c_{\kappa+2\eta}\rho^{|\kappa|+2\eta} + c_{\kappa+2\eta}^{(2)}\rho^{|\kappa|+2\eta + 2}+\dots\;$. As only the leading term $c_{\kappa+2\eta}\rho^{|\kappa|+2\eta}$ 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_{\delta m_{1(2)}} = \sum^{N-1}_{p = 0} a_{p}^{1(2)} u_{p,\delta m_{1(2)}}$, with $\sum (a_p^{1(2)})^2 = 1$, which allows to eliminate $2N-1$ higher orders $u_{\delta m_1}(\rho,0)\,u^*_{\delta m_2}(\rho,0) = c_{\kappa+2\eta}\rho^{|\kappa|+2\eta}(1 + O(\rho^{2N}))$. Hence, the scheme presented here becomes applicable for larger n without increasing the beam waist w₀. However, with increasing N the maximum of the ponderomotive potential moves further away from the vortex and, therefore, $U^{\kappa_a,\kappa_b}_{m_a ,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_{\kappa_a,m_a,m_b} \equiv U^{\kappa_a,\kappa_b = 0}_{m_a,m_b}$, the LG beams are chosen as $u_{\delta m_1}(\rho, z) = u_{\delta m_2}(\rho, z)$, with $\delta m_1 = \delta m_2 = \kappa_a/2$ and η = 0. The number of modes used for $\kappa_a = 2$ and 4 is N = 2 and for $\kappa_a = 6$ it is N = 3.

For the same LG laser beams the transition matrix elements $U^{\kappa_a,\kappa_b = 0}_{m_a = 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\kappa_a}$ scaling as the large Rydberg wavefunction probes the whole ponderomotive potential. Note, the scaling $n^{2\kappa_a}$ is obtained from equation (D7) and the fact that the angular momentum transition matrix element $U^{\kappa_a,\kappa_b = 0}_{m_a = 0,m_b = J}$ scales for $n\rightarrow \infty$ as $n^{\kappa_a}$. We furthermore point out that by tuning the relative laser phase ${\alpha}_1-{\alpha}_2$ the coupling matrix elements $U^{\kappa_a,\kappa_b}_{m_a = 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 $\delta {\omega} = -\kappa_a{\omega}_a + \kappa_b {\omega}_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 $\Gamma_T = I\sigma_T/{\omega}$, where $I = 2P/(\pi{\omega}_0^2)$ is the mean beam intensity and $\sigma_T = (8\pi/3)e^4/(4\pi{\epsilon}_0m_ec^2)^2$ is the Thomson scattering cross-section [161]. For $\lambda = 1300\,\mathrm{nm}$ and $w_0 = \lambda/2$ we obtain $\Gamma_T/P = 2\pi\times 0.7\,\mathrm{Hz}/\mathrm{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.