# Critical behavior of lattice gauge theory Rydberg simulators from effective
Hamiltonians

Jin Zhang1 S.-W. Tsai2 Y. Meurice3 1 Department of Physics and Chongqing Key
Laboratory for Strongly Coupled Physics, Chongqing University, Chongqing
401331, China  2Department of Physics and Astronomy, University of California,
Riverside, CA 92521, USA  3Department of Physics and Astronomy, University of
Iowa, Iowa City, IA 52242, USA

###### Abstract

We consider multileg ladders of Rydberg atoms which have been proposed as
quantum simulators for the compact Abelian Higgs model (CAHM) in 1+1
dimensions [Y. Meurice, Phys. Rev. D 104, 094513 (2021)] and modified versions
of theses simulators such as triangular prisms. Starting with the physical
Hamiltonian for the analog simulator, we construct translation-invariant
effective Hamiltonians by integrating over the simulator high-energy states
produced by the blockade mechanism when some of the atoms are sufficiently
close to each others. Remarkably, for all the simulators considered, the
effective Hamiltonians have the three types of terms present for the CAHM
(Electric field, matter charge and currents energies) but, in addition, terms
quartic in the electric field. For the two leg ladder, these additional terms
cannot be removed by fine-tuning the adjustable parameters of currently
available devices. For positive detuning, the new terms create highly-
degenerate vacua resulting in a very interesting phase diagram. Using
numerical methods, we demonstrate the close correspondence between the
physical simulator and the effective description for the ground state energy
and real-time evolution. We discuss the phase diagram at fixed geometry with
variable Rabi frequency and detuning and show that a rich variety of phases
can be reached with potential interest in the context of QCD at finite
density. We illustrate how the effective description can be used to design
simulators with desirable properties from the point of view of constructing
hybrid event generators.

##  I Introduction

The idea of using quantum devices to study the real-time evolution of strongly
interacting particles in high-energy and nuclear physics has gained
considerable interest in recent years Dalmonte and Montangero (2016); Bañuls
_et al._ (2020); Kasper  _et al._ (2020); Bauer  _et al._ (2023). In this
context, the possibility of building arrays of Rydberg atoms of significant
size with adjustable geometry and external parameters Bernien  _et al._
(2017); Keesling  _et al._ (2019); Browaeys and Lahaye (2020) offer many new
possibilities for analog simulations of lattice field theory models. In
addition, publicly available interfaces que ; pas  allow users to configure
arrays involving hundreds of Rydberg atoms and run their own experiments. This
sets the path for extensive empirical exploration by scientists who don’t have
direct access to this type of facilities.

This new technology has been used to propose simulations of spin and lattice
gauge theory models Celi  _et al._ (2020); Surace  _et al._ (2020);
Notarnicola  _et al._ (2020); Meurice (2021); Slagle  _et al._ (2021);
Heitritter  _et al._ (2022). One of the simplest model is the Abelian Higgs
model Bazavov  _et al._ (2015); Kuno  _et al._ (2017); González-Cuadra  _et
al._ (2017); Zhang  _et al._ (2018); Meurice (2021); Chanda  _et al._ (2022).
In 1+1 dimensions, and after elimination of the non-compact Brout-Englert-
Higgs mode, the Hamiltonian reads

$\hat{H}_{CAHM}=\frac{U}{2}\sum_{i=1}^{N_{s}}\left(\hat{L}^{z}_{i}\right)^{2}-Y{\sum_{i=1}^{N_{s}-1}}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}-X\sum_{i=1}^{N_{s}}\hat{U}^{x}_{i}\
$ (1)

with $N_{s}$ the number of sites, $\hat{L}^{z}|m\rangle=m|m\rangle$ and $m$ a
discrete electric field quantum number ($-\infty<m<+\infty$), and
$\hat{U}^{x}\equiv\frac{1}{2}(\hat{U}^{+}+\hat{U}^{-})$ with
$\hat{U}^{\pm}|m\rangle=|m\pm 1\rangle$. Note that in the derivation of the
Hamiltonian we obtain a term of the form
$\frac{Y}{2}\sum_{i}(\hat{L}^{z}_{i+1}-\hat{L}^{z}_{i})^{2}$ that accounts for
matter interactions, however for matching purposes which become clear later,
we have reabsorbed the local quadratic couplings in the definition of the
coupling $U$. In practice, we also apply truncations: for a spin-$m_{max}$
truncation we have $\hat{U}^{\pm}|\pm m_{max}\rangle=0$. In the following we
focus on the spin-1 truncation.

It has been argued Meurice (2021) that this spin-1 Hamiltonian can be
approximately simulated using arrays of Rydberg atoms by adapting the optical
lattice construction using ${}^{87}Rb$ atoms separated by controllable (but
not too small) distances, coupled to the excited Rydberg state $|r\rangle$
with a detuning $\Delta$. The ground state is denoted $|g\rangle$ and the two
possible states $|g\rangle$ and $|r\rangle$ can be seen as a qubit with
$n|g\rangle=0,\ n|r\rangle=|r\rangle$. The Hamiltonian for a generic array
reads

$\hat{H}=\frac{\Omega}{2}\sum_{i}(|g_{i}\rangle\langle
r_{i}|+|r_{i}\rangle\langle
g_{i}|)-\Delta\sum_{i}\hat{n}_{i}+\sum_{i<j}V_{ij}\hat{n}_{i}\hat{n}_{j},$ (2)

where the indices label atoms and $V_{ij}=\Omega R_{b}^{6}/r_{ij}^{6}$ for a
distance $r_{ij}$ between the atoms labelled as $i$ and $j$. Note: when
$r=R_{b}$, the Rydberg blockade radius within which states with two or more
Rydberg states have a low probability, we have $V=\Omega$. The simplest
simulator is a two-leg ladder Meurice (2021) with the correspondence:

$|gg\rangle\rightarrow|m=0\rangle,\ |gr\rangle\rightarrow|m=1\rangle,\
|rg\rangle\rightarrow|m=-1\rangle,$ (3)

for the two atoms on a rung. If the distance separating these two atoms is
small enough, the state $|rr\rangle$ is unlikely to appear. In Ref. Meurice
(2021), it was shown that a perfect matching could be obtained for the
individual sites, however, the nearest-neighbor (NN) interactions controlled
by $Y$ could not be matched exactly using the current technology where all the
interactions are repulsive.

In this article, we show that it is possible to construct an effective
Hamiltonian for the simulator provided that the set of atoms used to emulate
the local spin-1 degrees of freedom are close enough, in other words, if the
distance separating them is less than $R_{b}$. We will show that if these sets
form the rungs of a ladder and if the distance between the next-nearest-
neighbor (NNN) rungs is larger than $R_{b}$, the effective Hamiltonian has the
three types of terms found in Eq. (1) with in addition a quartic term of the
form $\sum_{i=1}^{N_{s}-1}(\hat{L}^{z}_{i})^{2}(\hat{L}^{z}_{i+1})^{2}$. The
ladder models considered and their effective Hamiltonians are presented in
Sec. II. Numerical tests are performed in Sec. III. Their phase diagrams are
discussed in Sec. IV. Practical applications involving systems that can be
simulated at facilities such as QuEra are discussed in Sec. V.

##  II Rydberg Simulator

Figure 1: The rectangular multileg ladder of Rydberg atoms, $a_{x(y)}$ is the
lattice spacing in $x(y)$ directions, $i$ labels rungs and $s$ labels legs.
Here it is a checkerboard stripe, where empty circles represent atomic ground
state $|g\rangle$, and the blue solid circles are excited Rydberg states
$|r\rangle$.

The multileg ladder of Rydberg atoms is a promising quantum simulator for
exotic quantum critical phenomena and lattice gauge theories, as we can use
different rung geometries to encode various local degrees of freedom. Figure 1
shows a rectangular multileg ladder of Rydberg atoms with $N_{s}$ rungs and
$N_{l}$ legs. The Hamiltonian of the multileg ladder reads

$\displaystyle\hat{H}_{m\rm{LR}}$ $\displaystyle=$
$\displaystyle\frac{\Omega}{2}\sum_{\begin{subarray}{c}i=1,2,\ldots,N_{s}\\\
s=1,2,\ldots,N_{l}\end{subarray}}(|g_{i,s}\rangle\langle
r_{i,s}|+h.c.)-\Delta\sum_{\begin{subarray}{c}i=1,2,\ldots,N_{s}\\\
s=1,2,\ldots,N_{l}\end{subarray}}\hat{n}_{i,s}$ (4) $\displaystyle+$
$\displaystyle\sum_{(i,s)\neq(i^{\prime},s^{\prime})}V_{i,s;i^{\prime},s^{\prime}}\hat{n}_{i,s}\hat{n}_{i^{\prime},s^{\prime}},$

where $|r_{i,s}\rangle$($|g_{i,s}\rangle$) is the Rydberg excited (ground)
state at the site on the $i$th rung and the $s$th leg,
$\hat{n}_{i,s}=|r_{i,s}\rangle\langle r_{i,s}|$ is the Rydberg number
operator, $\Omega$ and $\Delta$ are the Rabi frequency and detuning,
respectively. The interaction between Rydberg states are long-range van der
Waals repulsive interactions taking the form

$\displaystyle V_{i,s;i^{\prime},s^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{C_{6}}{\left[(i-i^{\prime})^{2}a_{x}^{2}+(s-s^{\prime})^{2}a_{y}^{2}\right]^{3}}$
(5) $\displaystyle=$
$\displaystyle\frac{V_{0}}{\left[(i-i^{\prime})^{2}/\rho^{2}+(s-s^{\prime})^{2}\right]^{3}}$

where $a_{x}$ ($a_{y}$) is the lattice spacing in $x$ ($y$) direction, $C_{6}$
is a constant, $V_{0}=C_{6}/a_{y}^{6}$, and $\rho=a_{y}/a_{x}$ is the inverse
aspect ratio. The Rydberg blockade mechanism means that at most one Rydberg
state is allowed with a significant probability within a sufficiently small
radius, typically the blockade radius $R_{b}$ that is defined by equating the
interaction energy at distance $R_{b}$ to the Rabi frequency,
$C_{6}/R_{b}^{6}=\Omega$, thus $R_{b}=\left(C_{6}/\Omega\right)^{1/6}$. Below
$R_{b}$, the interaction between Rydberg states is so strong that the laser
field cannot excite two Rydberg states simultaneously Šibalić and Adams
(2018). This imposes restrictions on the low-energy excitations for those
building blocks of quantum systems. The Rydberg system is highly programmable.
In our case, both the geometry of rungs and the parameters can be tuned to
encode local spin degrees of freedom, and realize a class of spin chains. In
the following, we label the leg by the spin projection quantum numbers.

###  II.1 Two-leg ladder

Figure 2: The two-leg Rydberg ladder with $N_{s}=2$ rungs. In each rung, the
state with one Rydberg state in the upper leg is labeled by the spin
$|+1\rangle$ state and the one with one Rydberg state in the lower leg is
labeled by the spin $|-1\rangle$ state. The spin $|0\rangle$ state labels the
state with no Rydberg state in the rung. The interactions between Rydberg
states in the same rung, in the same leg, and in different legs are $V_{0}$,
$V_{1}$, and $V_{2}$, respectively.

The one-dimensional Rydberg chain is a special case of the multi-leg ladder
system with $N_{l}=1$. For a two-leg ladder, $N_{l}=2$, there are four degrees
of freedom in each rung: $|g_{i,1}g_{i,2}\rangle$, $|g_{i,1}r_{i,2}\rangle$,
$|r_{i,1}g_{i,2}\rangle$, and $|r_{i,1}r_{i,2}\rangle$. The corresponding on-
rung energy is $0$, $-\Delta$, $-\Delta$, and $V_{0}-2\Delta$, respectively.
In principle, the four states can be used to represent the four projected
states in $z$ direction for spin-$3/2$, then the on-rung interaction can be
expressed by $\sum_{n=1}^{4}A_{n}\left(\hat{L}^{z}_{i}\right)^{n}$, where the
coefficients $A_{n}$ can be found by matching the energy spectrum. However, in
real quantum systems, the most common onsite terms are the linear term that is
coupled to the external field and the quadratic term that is the single-ion
anisotropy. In addition, the Rydberg interaction is strong when the rung size
is smaller than the Rydberg blockade radius such that the
$|r_{i,1}r_{i,2}\rangle$ state is not likely to appear. Here we only consider
a spin-$1$ realization by mapping the first three states to the spin-$1$
projected states $|0\rangle,|+1\rangle,|-1\rangle$, respectively. The onsite
interaction term in spin language is thus
$-\Delta\left(\hat{L}^{z}_{i}\right)^{2}$. The relation of the $z$-component
spin operator to the Rydberg number operator is defined as Meurice (2021)

$\hat{L}^{z}_{i}=\hat{n}_{i,+1}-\hat{n}_{i,-1}.$

If we take the square of this equation, use the property
$\hat{n}^{2}=\hat{n}$, and drop the term $\hat{n}_{i,+1}\hat{n}_{i,-1}$ which
is zero in the low energy sector we obtain effectively

$\left(\hat{L}^{z}_{i}\right)^{2}=\hat{n}_{i,+1}+\hat{n}_{i,-1}$ (6)

Solving for $\hat{n}_{i,m}$ we get

$\displaystyle\hat{n}_{i,+1(-1)}=\left[\left(\hat{L}^{z}_{i}\right)^{2}\pm\hat{L}^{z}_{i}\right]/2.$
(7)

Plugging Eq. 7 into the Rydberg interactions, we can easily write down the
interactions in terms of spin operators. If $\rho\ll 1$, we can just keep the
nearest-neighbor-rung interactions because of fast decaying of the van der
Waals interactions. Figure 2 shows the interactions between nearest-neighbor
(NN) rungs, which contains the Rydberg interactions between atoms in the same
leg $V_{1}=V_{0}\rho^{6}$ and those between atoms in different legs
$V_{2}=V_{0}\rho^{6}/\left(1+\rho^{2}\right)^{3}$. The NN spin interactions
are

$\displaystyle\hat{H}_{\rm
2LR,NN}=\frac{V_{1}-V_{2}}{2}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}+\frac{V_{1}+V_{2}}{2}\left(\hat{L}^{z}_{i}\right)^{2}(\hat{L}^{z}_{i+1})^{2}.$
(8)

For generic values of $\rho$, it is needed to include the long-range
interactions. The interactions between the spin at site $i$ and that at site
$i+k$ takes the same form as Eq. ( 8) by replacing $V_{1},V_{2}$ by
$V_{1}^{(k)}=V_{0}\rho^{6}/k^{6},V_{2}^{(k)}=V_{0}\rho^{6}/\left(k^{2}+\rho^{2}\right)^{3}$.

Finally, noticing that the Rydberg Rabi term can flip the spin projections
between $|0\rangle$ and $|\pm 1\rangle$, but there is no direct flipping
channel between $|+1\rangle$ and $|-1\rangle$, the Rabi term is equivalent to
the spin-1 ladder operator. In summary, if the rung size of the two-leg ladder
is smaller than the Rydberg blockade radius, or $V_{0}\gg\Delta,\Omega$, the
two-leg Rydberg ladder is an effective spin-1 chain

$\displaystyle\hat{H}_{\rm{2LR}}^{\rm{eff}}$ $\displaystyle=$
$\displaystyle-\Delta\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}+\sum_{k}\left(\frac{V_{1}^{(k)}-V_{2}^{(k)}}{2}\sum_{i=1}^{N_{s}-k}\hat{L}^{z}_{i}\hat{L}^{z}_{i+k}\right.$
(9)
$\displaystyle\left.+\frac{V_{1}^{(k)}+V_{2}^{(k)}}{2}\sum_{i=1}^{N_{s}-k}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+k}\right)^{2}\right)$
$\displaystyle+\frac{\Omega}{2}\sum_{i=1}^{N_{s}}\left(\hat{U}^{+}_{i}+\hat{U}^{-}_{i}\right).$

We now consider the case where the long-range interactions have a negligible
effect and keep only the NN interactions. In this situation the effective
Hamiltonian reads

$\displaystyle\hat{H}^{\rm{eff}}_{\rm{2LR}}$ $\displaystyle=$
$\displaystyle-\Delta\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}+\frac{V_{1}-V_{2}}{2}\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
(10)
$\displaystyle+\frac{V_{1}+V_{2}}{2}\sum_{i=1}^{N_{s}-1}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
$\displaystyle+\frac{\Omega}{2}\sum_{i=1}^{N_{s}}\left(\hat{U}^{+}_{i}+\hat{U}^{-}_{i}\right).$

The matching with the target model requires

  * •

$\Delta=-U/2$, note that the sign matters.

  * •

The coefficient for $\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$ is positive for the
simulator (repulsive/antiferromagnetic) but the CAHM has ferromagnetic
interactions. This can be remedied by redefining the observable
$\hat{L}^{z}_{2i+1}\rightarrow-\hat{L}^{z}_{2i+1}$ (staggered interpretation).

  * •

After this redefinition $V_{1}=-V_{2}=Y>0$ but $V_{2}>0$ with current
technology.

  * •

$\Omega=-X$, the sign does not matter because we can redefine the relative
phase between $|g\rangle$ and $|r\rangle$ without physical consequences.

These results agree with two-rung results of Meurice (2021).

###  II.2 Three-leg ladder

Figure 3: The three-leg Rydberg ladder with $N_{s}=2$ rungs. In each rung, the
three states that have only one Rydberg state in the upper, middle, or lower
leg are labeled by the spin $|+1\rangle$, $|0\rangle$, and $|-1\rangle$
states, respectively. In the same rung, the NN interaction is $V_{0}$, and the
NNN interaction is $V^{\prime}_{0}=V_{0}/64$. The inter-rung interactions
between atoms in the same leg, in NN legs, and in NNN legs are $V_{1}$,
$V_{2}$, and $V_{3}$, respectively.

Another scheme to realize spin-one chains is to use a three-leg ladder, where
the three states each with only one Rydberg state in a rung represent the
three spin projection states in $z$ direction Meurice (2021). The
configuration is shown in Fig. 3, where five different types of interactions
$V_{0}$, $V_{0}^{\prime}$, $V_{1}$, $V_{2}$, and $V_{3}$ exist, and we allow a
small offset for the detuning in the middle leg, $\Delta+\Delta_{0}$. There
are eight states in a rung, four of which have more than one Rydberg states,
and one of which has no Rydberg state. We consider two cases here.

In case $1$, the whole rung is within the Rydberg blockade radius
$2a_{y}<R_{b}$, only the four states with less than two Rydberg states are
allowed in the low energy band. We can further tune the detuning $\Delta$ such
that only the three states with one Rydberg state are allowed in the low-
energy band. In case 1, we require
$V_{0},V^{\prime}_{0}\gg\Delta\gg|\Delta_{0}|,|\Omega|$.

In case $2$, Rydberg states in NN legs in the same rung are blockaded by
letting $a_{y}<R_{b}$, but the size of the whole rung is larger than $R_{b}$,
and the state with two Rydberg states in the upper and the lower legs in the
same rung may has a lower energy than those states with one Rydberg state,
that is $-2\Delta+V_{0}^{\prime}<-\Delta$ if $V_{0}^{\prime}=V_{0}/64<\Delta$.
Thus the spin-$1$ sector is not in the lowest energy band. If the energy gap
$\Delta-V_{0}^{\prime}\gg\Omega$, the tunneling from the spin-1 sector to the
ground state is small and the $3$-leg Rydberg system can still simulate the
spin-$1$ dynamics with good accuracy. In case 2, we require
$V_{0},\Delta,|V_{0}-\Delta|\gg V^{\prime}_{0},|\Delta_{0}|,|\Omega|$.

For both cases, the three spin projection states are nearly degenerate and
form an energy band that. They are nearly orthogonal to other Rydberg states.
In both cases, in order to obtain the effective spin-$1$ Hamiltonian, we adapt
the relation between the $z$-component of the spin-$1$ operator and the
Rydberg number operator Meurice (2021), proceed as in Eq. (6), impose the
constraint that there is exactly one Rydberg state for the three spin-1
states, and obtain the equations

$\displaystyle\hat{L}^{z}_{i}$ $\displaystyle=$
$\displaystyle\hat{n}_{i,+1}-\hat{n}_{i,-1}$
$\displaystyle\left(\hat{L}^{z}_{i}\right)^{2}$ $\displaystyle=$
$\displaystyle\hat{n}_{i,+1}+\hat{n}_{i,-1}$ (11)
$\displaystyle\hat{n}_{i,-1}$ $\displaystyle+$
$\displaystyle\hat{n}_{i,0}+\hat{n}_{i,+1}=1$

When the detuning term is a uniform constant $-\Delta$, the constraint among
the occupation numbers in a rung implies that the detuning energy is a
constant. If we allow a small offset $\Delta_{0}$ to the detuning in the
middle leg Meurice (2021), the detuning energy $\hat{H}_{detun.,i}$ at rung
$i$ reads

$\hat{H}_{detun.,i}=\Delta_{0}(\hat{L}^{z}_{i})^{2}-\Delta-\Delta_{0}.$ (12)

The interaction term between NN rungs $\hat{H}_{\rm{3LR,NN}}$ can be rewritten
as

$\displaystyle\hat{H}_{\rm{3LR,NN}}$ $\displaystyle=$
$\displaystyle\left[\left(3V_{1}+V_{3}\right)/2-2V_{2}\right]\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
$\displaystyle+$
$\displaystyle\left[\left(V_{1}-V_{3}\right)/2\right]\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
$\displaystyle+$
$\displaystyle\left(V_{2}-V_{1}\right)\left[\left(\hat{L}^{z}_{i}\right)^{2}+\left(\hat{L}^{z}_{i+1}\right)^{2}\right]+V_{1}.$

The effective Hamiltonian for the Rabi term is different for the two cases
mentioned above. In case $1$, the whole rung is blockaded, the Rabi term can
flip between each two of the three spin states via hopping back and forth
between the empty state and the spin-$1$ sector. Notice that there also
possibility for the atom jumping out from and going back to the same state, so
the Rabi term effectively contributes both a diagonal term and a clock ladder
operator $\hat{C}^{\pm}_{i}$ which provides (anti)cyclic permutations.

In case $2$, the energy of the $|rgr\rangle$ state is about $|\Delta|$ lower
than the states in the spin-$1$ sector ($V_{0}^{\prime}$ is small), while the
energy of the empty state is the same amount higher. The hopping amplitudes of
flipping between $|+1\rangle$ and $|-1\rangle$ via the empty state and that
via the $|rgr\rangle$ state have the nearly the same magnitude but opposite
signs so are canceled. The Rabi term is effectively a spin-$1$ ladder operator
$\hat{U}^{\pm}_{i}$.

We can come to the same conclusion using perturbation theory to calculate the
effective coupling constants. Treating the Rabi term as a perturbation, the
effective Hamiltonian for the Rabi term is

$\displaystyle\langle
e^{\prime}|\hat{H}^{\rm{eff}}_{\rm{3LR},\Omega}|e\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{f\neq e,e^{\prime}}\langle
e^{\prime}|\hat{H}_{\rm{3LR},\Omega}|f\rangle\langle
f|\hat{H}_{\rm{3LR},\Omega}|e\rangle$ (14)
$\displaystyle\left(\frac{1}{E_{e^{\prime}}-E_{f}}-\frac{1}{E_{f}-E_{e}}\right),$

where $|e\rangle$($|f\rangle$) is the state inside(outside) of the spin-1
sector. The effective matrix for the Rabi term in the spin-$1$ sector is

$\displaystyle-\frac{\Omega^{2}}{4}\begin{pmatrix}\mathrm{A}&\Gamma&\Lambda\\\
\Gamma&\mathrm{B}&\Gamma\\\ \Lambda&\Gamma&\mathrm{A}\end{pmatrix}$ (15)

where

$\displaystyle\mathrm{A}$ $\displaystyle=$
$\displaystyle\frac{1}{V_{0}-\Delta-\Delta_{0}}+\frac{1}{V_{0}^{\prime}-\Delta}+\frac{1}{\Delta}$
$\displaystyle\mathrm{B}$ $\displaystyle=$
$\displaystyle\frac{2}{V_{0}-\Delta}+\frac{1}{\Delta+\Delta_{0}}$
$\displaystyle 2\Gamma$ $\displaystyle=$
$\displaystyle\frac{1}{\Delta}+\frac{1}{V_{0}-\Delta}+\frac{1}{\Delta+\Delta_{0}}+\frac{1}{V_{0}-\Delta-\Delta_{0}}$
$\displaystyle\Lambda$ $\displaystyle=$
$\displaystyle\frac{1}{V^{\prime}_{0}-\Delta}+\frac{1}{\Delta}$ (16)

The diagonal part of the effective matrix can be written as
$[(\mathrm{B}-\mathrm{A})\Omega^{2}/4](\hat{L}^{z}_{i})^{2}$ up to a constant.
For case $1$, $\Gamma\approx\Lambda\approx 1/\Delta$, the off-diagonal part is
a three-state clock ladder operator. For case $2$, $\Lambda\approx 0$, the
off-diagonal part is a spin-$1$ ladder operator. Thus effectively

$\displaystyle\hat{H}^{\rm{eff}}_{\rm{3LR},\Omega}=\begin{cases}-\frac{\Omega^{2}}{4\Delta}\sum_{i}\left(\hat{C}_{i}^{+}+\hat{C}_{i}^{-}\right)&\text{for
case 1,}\\\
-\frac{\Omega^{2}\Gamma}{4}\sum_{i}\left(\hat{U}^{+}_{i}+\hat{U}^{-}_{i}\right)&\text{for
case 2.}\end{cases}$ (17)

In summary, one can write down the effective spin-$1$ Hamiltonian for the
three-leg Rydberg ladder

$\displaystyle\hat{H}^{\rm{eff}}_{\rm{3LR}}$ $\displaystyle=$
$\displaystyle\left[\Delta_{0}+\frac{\Omega^{2}\left(\mathrm{B}-\mathrm{A}\right)}{4}\right]\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}$
(18) $\displaystyle+$
$\displaystyle\left(V_{2}-V_{1}\right)\sum_{i=1}^{N_{s}-1}\left[\left(\hat{L}^{z}_{i}\right)^{2}+\left(\hat{L}^{z}_{i+1}\right)^{2}\right]$
$\displaystyle+$
$\displaystyle\frac{V_{1}-V_{3}}{2}\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
$\displaystyle+$
$\displaystyle\left(\frac{3V_{1}+V_{3}}{2}-2V_{2}\right)\sum_{i=1}^{N_{s}-1}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
$\displaystyle+$
$\displaystyle\hat{H}^{\rm{eff}}_{\rm{3LR},\Omega}+\rm{Const.}\hat{I},$

where the constant
$\rm{Const.}=-(\Delta+\Delta_{0})N_{s}-N_{s}\Omega^{2}B/4+V_{1}(N_{s}-1)$.

###  II.3 Triangular prism

Figure 4: Same as Fig. 3, but the middle leg is put out of plane to form a
equilateral triangular prism.

If we allow the middle leg of the three-leg ladder moving out of plane, the
effective Hamiltonian for case 1 does not change as long as the size of the
whole rung is within $R_{b}$. In particular, if the three sites in a rung form
an equilateral triangle (Fig. 4), $V_{2}=V_{3}$, and
$V_{0}=V^{\prime}_{0}\gg\Delta\gg|\Delta_{0}|,|\Omega|$. Thus
$\Gamma\approx\Lambda$, $B\approx A$, the effective Hamiltonian for the
Rydberg simulator is

$\displaystyle\hat{H}^{\rm{eff}}_{\rm{prism}}$ $\displaystyle=$
$\displaystyle\Delta_{0}\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}+\left(V_{2}-V_{1}\right)\sum_{i=1}^{N_{s}-1}\left[\left(\hat{L}^{z}_{i}\right)^{2}+\left(\hat{L}^{z}_{i+1}\right)^{2}\right]$
(19) $\displaystyle+$
$\displaystyle\frac{V_{1}-V_{2}}{2}\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
$\displaystyle+$
$\displaystyle\frac{3(V_{1}-V_{2})}{2}\sum_{i=1}^{N_{s}-1}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
$\displaystyle-$
$\displaystyle\frac{\Omega^{2}V_{0}}{4\Delta(V_{0}-\Delta)}\sum_{i=1}^{N_{s}}\left(\hat{C}^{+}_{i}+\hat{C}^{-}_{i}\right)+\rm{Const.}\hat{I}.$

Here, the coefficient of $\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$ and that of
$(\hat{L}^{z}_{i})^{2}(\hat{L}^{z}_{i+1})^{2}$ have fixed ratio, the physics
originated from the competition between the two terms may be missing in the
Hamiltonian. But the blockade radius required is just half of the one in Fig.
3.

###  II.4 Triangles in the same plane

Figure 5: Same as Fig. 4, but triangles reside in the same plane.

As the three dimensional triangular prism is not easy to realize in
experiment, we can just shift the middle leg in the same plane. As shown in
Fig. 5, we shift the middle leg to left (or right). Now the off-rung
interactions between nearest-neighbor legs depend on its location and take two
values $V_{2},V_{4}$. The effective Hamiltonian reads

$\displaystyle\hat{H}^{\rm{eff}}_{\rm{in-plane}}$ $\displaystyle=$
$\displaystyle\Delta_{0}\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}$
(20) $\displaystyle+$
$\displaystyle\sum_{i=1}^{N_{s}-1}\left[\left(V_{2}-V_{1}\right)\left(\hat{L}^{z}_{i}\right)^{2}+\left(V_{4}-V_{1}\right)\left(\hat{L}^{z}_{i+1}\right)^{2}\right]$
$\displaystyle+$
$\displaystyle\frac{V_{1}-V_{3}}{2}\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
$\displaystyle+$
$\displaystyle\left(\frac{3V_{1}+V_{3}}{2}-V_{2}-V_{4}\right)\sum_{i=1}^{N_{s}-1}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
$\displaystyle-$
$\displaystyle\frac{\Omega^{2}V_{0}}{4\Delta(V_{0}-\Delta)}\sum_{i=1}^{N_{s}}\left(\hat{C}^{+}_{i}+\hat{C}^{-}_{i}\right)+\rm{Const.}\hat{I}.$

By shifting the middle leg, we have more degrees of freedom to tune the
coefficients of the interaction terms: if $V_{0}$, $\Delta$ and $\Omega$ are
fixed, we can tune $\Delta_{0}$, $a_{x}$ and the angle of the triangle to tune
couplings of $(\hat{L}^{z}_{i})^{2}$, $\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$, and
$(\hat{L}^{z}_{i})^{2}(\hat{L}^{z}_{i+1})^{2}$ terms.

##  III Numerical tests

We have derived the low-energy effective Hamiltonian for various Rydberg
ladder systems using perturbation theory. All the effective spin-1
Hamiltonians obtained in Sec. II have the same form

$\displaystyle\hat{H}_{\rm{eff}}$ $\displaystyle=$ $\displaystyle
D\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}+R\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
(21) $\displaystyle+$ $\displaystyle
R^{\prime}\\!\\!\sum_{i=1}^{N_{s}-1}\\!\\!\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}\\!-\\!J\\!\sum_{i=1}^{N_{s}}\left(\hat{U}^{+}_{i}+\hat{U}^{-}_{i}\right),$

where the ladder operator $\hat{U}^{\pm}_{i}$ can be replaced by the clock
raising and lowering operator $\hat{C}^{\pm}_{i}$ for the cases that allow
hopping between $|\pm 1\rangle$ states. In this section, we provide some
numerical calculations to test the validity of our results.

###  III.1 Two-leg ladder

Figure 6: The relative error of the ground-state energy
$E^{\rm{eff}}_{0,\rm{2LR}}$ of the effective Hamiltonian (10) as a function of
$\Omega$. All the van der Waals interactions are kept in the original two-leg
Rydberg ladder to compute the actual ground-state energy $E_{0,\rm{2LR}}$.
Here, $N_{s}=16$, $V_{0}=1000\times 2\pi$ MHz, $\Delta=1\times 2\pi$ MHz, and
two inverse aspect ratios $\rho=0.5,0.4$ are considered. The inset shows the
ground-state energy as a function of $\Omega$. The energy for $\rho=0.4$ is
shifted by $-1$ for better view.

In this case, the effective Hamiltonian is shown in Eq. (10) with $D=-\Delta$,
$R=(V_{2}-V_{1})/2$, $R^{\prime}=(V_{1}+V_{2})/2$, $J=-\Omega/2$. Since the
error of the effective Hamiltonian is of order $\Omega^{2}/(4V_{0})$, we set
$V_{0}$ to a large value. Given that $V_{ij}=C_{6}/r_{ij}^{6}$ with
$C_{6}=858386\times 2\pi$ MHz $\mu$m6, we set $V_{0}=1000\times 2\pi$ MHz by
taking the length of the rung to be $a_{y}=3.083\mu$m. Figure 6 presents the
ground-state energy difference between the two-leg Rydberg ladder
$\hat{H}_{\rm{2LR}}$ and the corresponding effective Hamiltonian (10) as a
function of $\Omega$ for $N_{s}=16$, $\Delta=1\times 2\pi$ MHz. The effective
Hamiltonian only contain NN interactions, while $\hat{H}_{\rm{2LR}}$ contains
all van der Waals interactions. If the inverse aspect ratio $\rho=0.5$, the
relative energy difference is about $18\%$ for small $\Omega=0.2$, and
decreases quickly toward zero as $\Omega$ is increased. For smaller inverse
aspect ratio $\rho=0.4$, the relative error behaves the same way but with an
overall smaller magnitude. Note that the error from perturbation theory is of
order $\Omega^{2}/(4V_{0})$, which should increase as we increase $\Omega$,
but is small compared to the error from the omitted long-range interactions.
For small $\Omega$, the blockade radius $R_{b}/a_{y}=(V_{0}/\Omega)^{1/6}$ is
large and the interactions beyond NN rungs play important roles and are not
negligible. As we increase $\Omega$, the blockade radius $R_{b}$ is decreased
and the effects of long-range interaction fade away. By fixing $a_{y}$ and
increasing $a_{x}$ such that $\rho$ is decreased, $R_{b}$ is becoming
relatively smaller, so are the effects of long-range interaction. For
$\Omega=10$, $R_{b}\approx 2.154a_{y}$, slightly larger than $a_{x}=2a_{y}$
for $\rho=0.5$ and smaller than $a_{x}=2.5a_{y}$ for $\rho=0.4$. We conclude
that it is safe to omit long-range interactions when the blockade radius
$R_{b}$ is small.

###  III.2 Three-leg ladder

We rewrite the Hamiltonian (18) with 00BC as

$\displaystyle\hat{H}_{\rm{eff1}}$ $\displaystyle=$ $\displaystyle
D\sum_{i=1}^{N_{s}}\left(\hat{L}_{i}^{z}\right)^{2}+R\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}$
(22) $\displaystyle+$ $\displaystyle
R^{\prime}\\!\\!\sum_{i=1}^{N_{s}-1}\\!\\!\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}\\!-\\!J\\!\sum_{i=1}^{N_{s}}\left(\hat{U}^{+}_{i}+\hat{U}^{-}_{i}\right),$

where $D=\Delta_{0}+J+2(V_{2}-V_{1})$, $R=(V_{3}-V_{1})/2$,
$R^{\prime}=(V_{3}+3V_{1})/2-2V_{2}$,
$J=\Omega^{2}V_{0}/[4\Delta(V_{0}-\Delta)]$. We can use exact diagonalization
to show that the Rydberg system and the effective Hamiltonian have close
energy density for small system sizes. But because the spin-1 sector is not
the lowest energy band, we cannot check the effective Hamiltonian using normal
DMRG which finds the ground state. We have seen that the effective Hamiltonian
for triangular chain works well, so we believe Hamiltonian (22) well describe
the physics of the spin-1 sector.

We can calculate time evolution of an initial state to check the effective
Hamiltonian. In Fig. 7, we show the time evolution of the initial state
$|0,0,0\rangle$ for $N_{s}=3$, $V_{0}=2\Delta=40\times 2\pi MHz$,
$\Omega=2\times 2\pi MHz$, $\rho=1/3$. Three values of $\Delta_{0}/2\pi MHz$
are considered, $0,0.2,0.4$, which corresponds to $(D,R,R^{\prime},J)=$
$(0.070261,-0.018332,0.011408,0.1),$ $(0.270261,-0.018332,0.011408,0.1),$
$(0.470261,-0.018332,0.011408,0.1)$

Figure 7: Compare the real-time evolution of the initial state $|0,0,0\rangle$
with the Rydberg simulator to that with the effective Hamiltonian (22),
$N_{s}=3$, $\Delta_{0}=(0,0.2,0.4)$, $V_{0}=2\Delta=40$, $\Omega=2$,
$\rho=1/3$. The parameters are in unit of $2\pi MHz$.

###  III.3 Equidistant equilateral triangles in plane

We have shown that the spin-1 effective Hamiltonians of the two-leg and three-
leg Rydberg array have the same form

$\displaystyle\hat{H}_{\rm{S=1}}$ $\displaystyle=$ $\displaystyle
D_{1}\left(L_{1}^{z}\right)^{2}+D\sum_{i=2}^{N_{s}-1}\left(\hat{L}_{i}^{z}\right)^{2}+D_{N_{s}}\left(L_{N_{s}}^{z}\right)^{2}$
(23) $\displaystyle+$ $\displaystyle
R\sum_{i=1}^{N_{s}-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}+R^{\prime}\sum_{i=1}^{N_{s}-1}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
$\displaystyle-$ $\displaystyle
J\sum_{i=1}^{N_{s}}\left(\hat{C}^{+}_{i}+\hat{C}^{-}_{i}\right),$

where $D=\Delta_{0}+V_{2}+V_{4}-2V_{1}$, $R=(V_{3}-V_{1})/2$,
$R^{\prime}=(V_{3}+3V_{1})/2-V_{2}-V_{4}$,
$J=\Omega^{2}V_{0}/[4\Delta(V_{0}-\Delta)]$. We choose
$D_{1}=\Delta_{0}/2+V_{2}-V_{1}$ and $D_{N_{s}}=\Delta_{0}/2+V_{4}-V_{1}$ to
minimize the boundary effects. Because if $D_{1}=D_{N_{s}}=D$, the diagonal
part has an energy
$DN_{s}+R(N_{s}-1)+R^{\prime}(N_{s}-1)=(D+R+R^{\prime})N_{s}-R-R^{\prime}$ for
the ferromagnetic phase, which may has different signs for small and large
system sizes. Calculations for small system may give wrong phase.

We compare the ground-state energy for the original Rydberg system and its
effective Hamiltonian in Fig. 8. Notice that $V_{2}(\rho=0.4)\approx 4$ and
$V_{2}(\rho=1/3)\approx 0.9$, the effective Hamiltonian well describe the low-
energy physics for $V_{0}=2\Delta=100\times 2\pi$ MHz, $\Omega\leq 5\times
2\pi$ MHz. Thus Hamiltonian (22) will also well describe the spin-1 sector for
the three-leg ladder with $\rho\leq 0.5$.

Figure 8: Ground-state energy as a function of $\Omega$ for the Rydberg
simulator and the effective Hamiltonian (23), $N_{s}=8$, $\Delta_{0}=0$,
$V_{0}=2\Delta=100\times 2\pi$ MHz, (top) $\rho=0.4$, (bottom) $\rho=1/3$.

For $\Delta_{0}=0$, $V_{0}=2\Delta=100\times 2\pi$ MHz, and $\rho=0.4$, the
values of ($D$, $R$, $R^{\prime}$) are ($3.26225434954982$,
$-0.07359330845873141$, $-3.335847658008552$). The phase transition can be
characterized by the second-order Rényi entropy.

##  IV Phase Diagrams

All the effective Hamiltonians for ladders have the same form

$\displaystyle\hat{H}_{\rm{eff}}$ $\displaystyle=$ $\displaystyle
D\sum_{i=1}^{L}\left(\hat{L}_{i}^{z}\right)^{2}+R\sum_{i=1}^{L-1}\hat{L}^{z}_{i}\hat{L}^{z}_{i+1}+R^{\prime}\sum_{i=1}^{L-1}\left(\hat{L}^{z}_{i}\right)^{2}\left(\hat{L}^{z}_{i+1}\right)^{2}$
(24)
$\displaystyle-J\sum_{i=1}^{L}\left(\hat{U}^{+}_{i}+\hat{U}^{-}_{i}\right),$

When $D$ and $R^{\prime}$ are large positive, all spins are zero, which is the
large-$D$ phase. When $R$ is large positive, the anti-ferromagnetic (AFM)
phase is favored. When $R$ is large negative, the ferromagnetic (FM) phase is
favored. The competition of the three diagonal terms may also favors the
Rydberg state density wave (RDW) phase: the values of
$\langle(\hat{L}^{z}_{i})^{2}\rangle$ are $1,0,1,0,\ldots$. The order
parameters are $M_{\rm{FM}}=1/L\sum_{i}\hat{L}^{z}_{i}$,
$M_{\rm{AFM}}=1/L\sum_{i}(-1)^{i}\hat{L}^{z}_{i}$, and
$M_{\rm{RDW}}=1/L\sum_{i}(-1)^{i}(\hat{L}^{z}_{i})^{2}$ for FM, AFM, and RDW,
respectively.

###  IV.1 Two-leg ladder: $a_{x}=2a_{y}$

Table 1: The values of order parameters for different phases. The corresponding susceptibilities diverge on the critical lines. |  $M_{\rm{FM}}$ |  $M_{\rm{AFM}}$ |  $M_{\rm{RDW}}$
---|---|---|---
FM |  $\neq 0$ |  $=0$ |  $=0$
AFM |  $=0$ |  $\neq 0$ |  $=0$
FRDW |  $\neq 0$ |  $\neq 0$ |  $\neq 0$
PRDW |  $=0$ |  $=0$ |  $\neq 0$
Disorder |  $=0$ |  $=0$ |  $=0$

$\displaystyle H_{2LR}$ $\displaystyle=$
$\displaystyle\frac{\Omega}{2}\sum_{i;s=1,2}(|g_{is}\rangle\langle
r_{is}|+h.c.)-\Delta\sum_{i;s=1,2}\hat{n}_{i,s}$ (26)
$\displaystyle+\sum_{is\neq
i^{\prime}s^{\prime}}V_{is,i^{\prime}s^{\prime}}\hat{n}_{i,s}\hat{n}_{i^{\prime},s^{\prime}}.$

Labeling the upper leg as spin $|+1\rangle$ state and the lower leg as spin
$|-1\rangle$ state,
$D=-\Delta,R=(V_{1}-V_{2})/2,R^{\prime}=(V_{1}+V_{2})/2,J=-\Omega/2$. This
corresponds to Eq. (10).

Figure 9: Ground-state phase diagram for the effective Hamiltonian (22) of the
two-leg Rydberg ladder. Here $L=512$, $V_{0}=1000\times 2\pi$ MHz,
$\rho=a_{y}/a_{x}=0.5$. The PRDW phase is disordered in even or odd sites, and
the FRDW phase is FM in even or odd sites.

When $\Omega=0$, the diagonal part of the Hamiltonian is in the Large-$D$
phase for $\Delta<0$, the RDW phase for $0<\Delta<2V_{2}$, and the AFM phase
for $\Delta>2V_{2}$. The RDW phase has high degeneracy with states
$|\ldots,\pm 1,0,\pm 1,0,\ldots\rangle$. When $\Omega$ is turned on a little,
the RDW phase splits into the ferromagnetic RDW (FRDW) phase
($|\ldots,1,0,1,0,\ldots\rangle$ or $|\ldots,-1,0,-1,0,\ldots\rangle$) and
paramagnetic RDW (PRDW) phase ($\langle\hat{L}^{z}_{i}\rangle=0$ but
$\langle(\hat{L}^{z}_{i})^{2}\rangle=1$ on even or odd sites). The critical
line between FRDW and PRDW is an Ising critical line. One can perform
degenerate perturbation theory for small $\Omega$ and obtain an effective
Ising model on even or odd sites:

$H_{\rm{Ising}}=\frac{\Omega^{2}}{4}(-J_{eff}\sum_{i}\sigma^{z}_{i}\sigma^{z}_{i+2}-\frac{1}{\Delta}\sum_{i}\sigma^{x}_{i}).$
(27)

with

$J_{eff}=(\frac{1}{\Delta-
V_{1}-V_{2}}-\frac{1}{2\Delta-4V_{2}}-\frac{1}{2\Delta-4V_{1}})$ (28)

One can obtain the critical point between FRDW and PRDW in the
$\Omega\rightarrow 0$ limit by equating the two coefficients. For large
$\Omega$, there maybe quantum phase transitions between ordered phases and the
disordered phase. Table 1 shows how to identify different phases using order
parameters. The spin operators and the Rydberg state operators are related by
Eqs. II.1 and 6.

Figure 10: Finite-size scaling of (top) the peak position of the PRDW
susceptibility and (bottom) the first- and second-order Rényi entropy with
PBC. The results are for the effective Hamiltonian (22) of the two-leg Rydberg
ladder. Here $V_{0}=1000\times 2\pi$ MHz, $\rho=0.5$, $\Delta=3$. The critical
point separates the PRDW phase and the disordered phase, which belongs to the
Ising universality class. Figure 11: Same as Fig. 10, but for
$\Delta=15,\Omega_{c}=1.2531(1)$. The critical point separates FRDW phase and
AFM phase, which is also in Ising universality class.

This suggest to study $\Delta>0$ as an environment effect, staggered
structures (cells with two sites) and to explore outside the region of
validity of $H_{eff.}$

##  V Practical Applications

In this section, we discuss simple examples of real-time evolution for two leg
ladders with 10 atoms (5 sites). This situation can be implemented remotely
using QuEra facilities or simulated with a local SDK que ; aws . In this case,
the system is initialized with the 10 atoms in the ground state. This initial
state is invariant under the transformation
$\hat{L}^{z}_{i}\rightarrow-\hat{L}^{z}_{i}$ or equivalently swapping the two
legs of the ladder and as we turn on the evolution the $<\hat{L}^{z}_{i}>$
remain zero. For this reason, we display the values of $<(L^{z}_{i})^{2}>$ for
the five sites which is also invariant under staggered redefinitions of
$\hat{L}_{i}^{z}$. Note also that the results are left-right symmetric.
Results for the typical QuEra values $\Omega=4\pi MHz$, $\Delta=2\Omega$,
$\rho=0.5$ and the distance between the sites $a_{x}=1R_{b}\simeq 8.7\mu m$
are shown below. The distance between the two atoms on one site is
$a_{y}=0.5R_{b}$ so $|rr\rangle$ at that site are unlikely. The evolution for
a period 0.5 $\mu$s is shown in Fig. 12, with exact diagonalization for the 10
site problem, the QuEra local simulator and the effective Hamiltonian.

Figure 12:  Values of $<(\hat{L}^{z}_{i})^{2}>$ for five sites, $\Omega=4\pi$
MHz, $\Delta=2\Omega$, $\rho=0.5$ $a_{x}=1R_{b}$ and $a_{y}=0.5\times R_{b}$.
The vertical time units are $10^{-8}$ s. Top:exact diagonalization. Middle:
QuEra (local simulator). Bottom: Effective Hamiltonian.

It appears clearly that the three methods give very similar results. A finer
comparison for the first three sites in Fig. 13 confirms that the effective
Hamiltonian results are in excellent agreement.

Figure 13:  Values of $<(\hat{L}^{z}_{i})^{2}>$ for the first three of the
five sites.

We have also explored values of the parameters that are slightly outside the
range of validity of $H_{eff.}$. As a first example we keep $\rho=0.5$, but we
reduce the lattice spacing to $a_{x}=R_{b}/2$. The bloquade mechanism is even
more effective but $a_{x}=R_{b}/2$ implies next to nearest neighbor
interactions not taken into account in the simple $H_{eff.}$ with just NN
interactions as in Eq. (10). Fig. 14 makes clear that the effective
description is a less accurate description of the simulator than in the
previous case, however, there is a qualitative agreement between the exact
diagonalization with 10 atoms and the corresponding $H_{eff.}$.

Figure 14: 5 sites, $\rho=0.5$, $a_{x}=0.5\times R_{b}$, $\Omega=4\pi$MHz,
$\Delta=2\Omega$, vertical time units are $0.01\mu$sec. Top: 10 atom
simulator. Bottom: 5 site $H_{eff.}$.

It should also be noticed that for $a_{x}=R_{b}/2$, $V_{1}$ and $V_{2}$ are 64
times larger than in the previous case and the quartic term becomes very
important. This is signaled by $L_{z}=0$ bands screening the electric field
both in the effective theory and the original simulator. This is illustrated
in Fig. 15.

Figure 15: 5 sites, $\rho=0.5$, $a_{x}=0.5\times R_{b}$, $\Omega=4\pi$MHz,
$\Delta=2\Omega$, vertical time units are $0.01\mu$sec. Top: 10 atom
simulator. Bottom: 5 site $H_{eff.}$.

The screening effect observed in the previous deformation can be partially
remediated by rather increasing $\rho$ in order to have $V_{2}<<V_{1}$. For
this reason we considered the case $\rho=2$, $a_{x}=R_{b}$. The results are
shown in Fig. 16.

Figure 16: 5 sites, $\rho=2$, $a_{x}=1\times R_{b}$, $\Omega=4\pi$MHz,
$\Delta=2\Omega$, vertical time units are $0.01\mu$sec. Top: 10 atom
simulator. Bottom: 5 site $H_{eff.}$.

##  VI Conclusion

In summary we have considered ladder-shaped Rydberg arrays with two or three
atoms per site. Originally these simulators were designed with the idea of
mimicing closely the evolution of the compact Abelian Higgs model. We
constructed an effective Hamiltonian valid when the size of the rungs is small
enough and the distance between the rungs not too small. In all cases, we
found that the effective Hamiltonians have the same three types of terms as
the target model plus an extra quartic term. The effect of the extra term is
significant at positive detuning and is responsible for new quantum phases.
More generally, the ladder models have a very rich phase diagram that is
currently being explored using QuEra que . Matching the effective theory with
the target gauge theory requires $\Delta<0$ (cost for producing electric
field). The new phases appear when $\Delta>0$ ($m=\pm 1$ form the degenerate
ground state). It is possible that a positive detuning could be used as an
environment effect relevant in the context of hybrid hadronization Heitritter
_et al._ (2022). Potential issue with microscopic string breaking generated by
extra terms with a large coefficient suggest to study staggered structures
(cells with two sites) that could be used to describe models different from
the target CAHM.

It should also be noted that here has been interest in inhomogeneous phases
and the Lifshitz regime for QCD at finite temperature and density Pisarski
_et al._ (2019). It has been argued that Kojo  _et al._ (2010) that for
massless quarks in the large $N_{c}$ limit, dimensional reduction occurs and
“chiral spiral” condensates appear Basar and Dunne (2008). In this scenario,
the chiral spiral phase could appear at the end of the crossover line and
separate an ordered phase where chiral symmetry is restored from the hadronic
(confining) phase where chiral symmetry is broken (see Fig. 6 in Pisarski  _et
al._ (2019)). In this context, exploring the possibility of inhomogeneous
phase in simulators is an interesting direction of research.

In the future, it would also be interesting to have a simulator where negative
(attractive) couplings among atoms could be engineered. This would allow us to
cancel the extra terms and to have an effective Hamiltonian identical to the
target Hamiltonian. Similar technological needs to be present for quantum
simulation of the tricritical Ising model (with differrent types of NNN
interactions being negative) Slagle  _et al._ (2021). As we are moving toward
better local control of the individual atoms we should look forward to using
new technologies for lattice gauge theory. It would also be interesting to
compare the manipulation of the three states associated with a rung in our
approach with qutrit simulations Gustafson (2021, 2022) and figure out if the
extra quartic term found here could be understood in the context of Symanzik
improvement Carena  _et al._ (2022a, b)

###### Acknowledgements.

We thank Sergio Cantu, James Corona, Fangli Liu, Kenny Heitritter, Steve
Mrenna, Shengtao Wang and members of QuLAT for helpful discussions. This work
was supported in part by the National Science Foundation (NSF) RAISE-TAQS
under Award Number 1839153 (SWT). Computations were performed using the
computer clusters and data storage resources of the HPCC, which were funded by
grants from NSF (MRI-1429826) and NIH (1S10OD016290-01A1). Y.M. is supported
in part by the Dept. of Energy under Award Number DE-SC0019139. J.Z. is
supported by NSFC under Grants No. 12304172 and No. 12147102, Chongqing
Natural Science Foundation under Grant No. CSTB2023NSCQ-MSX0048, and
Fundamental Research Funds for the Central Universities under Projects No.
2023CDJXY-048 and No. 2020CDJQY-Z003. Y. M. thanks the Amazon Web Services and
S. Hassinger for facilitating remote access to QuEra through the Amazon
Braket.

## References

  * Dalmonte and Montangero (2016) M. Dalmonte and S. Montangero, Contemp. Phys. 57, 388 (2016), arXiv:1602.03776 [cond-mat.quant-gas] .
  * Bañuls  _et al._ (2020) M. C. Bañuls _et al._ , Eur. Phys. J. D 74, 165 (2020), arXiv:1911.00003 [quant-ph] .
  * Kasper  _et al._ (2020) V. Kasper, G. Juzeliunas, M. Lewenstein, F. Jendrzejewski, and E. Zohar, New J. Phys. 22, 103027 (2020), arXiv:2006.01258 [quant-ph] .
  * Bauer  _et al._ (2023) C. W. Bauer _et al._ , PRX Quantum 4, 027001 (2023), arXiv:2204.03381 [quant-ph] .
  * Bernien  _et al._ (2017) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature (London) 551, 579 (2017), arXiv:1707.04344 [quant-ph] .
  * Keesling  _et al._ (2019) A. Keesling _et al._ , Nature 568, 207 (2019), arXiv:1809.05540 [quant-ph] .
  * Browaeys and Lahaye (2020) A. Browaeys and T. Lahaye, Nature Phys. 16, 132 (2020).
  * (8) “Quera,” https://www.quera.com.
  * (9) “Pasqal,” https://www.pasqal.com.
  * Celi  _et al._ (2020) A. Celi, B. Vermersch, O. Viyuela, H. Pichler, M. D. Lukin, and P. Zoller, Phys. Rev. X 10, 021057 (2020), arXiv:1907.03311 [quant-ph] .
  * Surace  _et al._ (2020) F. M. Surace, P. P. Mazza, G. Giudici, A. Lerose, A. Gambassi, and M. Dalmonte, Phys. Rev. X 10, 021041 (2020), arXiv:1902.09551 [cond-mat.quant-gas] .
  * Notarnicola  _et al._ (2020) S. Notarnicola, M. Collura, and S. Montangero, Phys. Rev. Res. 2, 013288 (2020), arXiv:1907.12579 [cond-mat.quant-gas] .
  * Meurice (2021) Y. Meurice, Phys. Rev. D 104, 094513 (2021), arXiv:2107.11366 [quant-ph] .
  * Slagle  _et al._ (2021) K. Slagle, D. Aasen, H. Pichler, R. S. K. Mong, P. Fendley, X. Chen, M. Endres, and J. Alicea, Phys. Rev. B 104, 235109 (2021), arXiv:2108.09309 [cond-mat.str-el] .
  * Heitritter  _et al._ (2022) K. Heitritter, Y. Meurice, and S. Mrenna, (2022), arXiv:2212.02476 [quant-ph] .
  * Bazavov  _et al._ (2015) A. Bazavov, Y. Meurice, S.-W. Tsai, J. Unmuth-Yockey, and J. Zhang, Phys. Rev. D92, 076003 (2015), arXiv:1503.08354 [hep-lat] .
  * Kuno  _et al._ (2017) Y. Kuno, S. Sakane, K. Kasamatsu, I. Ichinose, and T. Matsui, Phys. Rev. D 95, 094507 (2017), arXiv:1605.00333 [cond-mat.quant-gas] .
  * González-Cuadra  _et al._ (2017) D. González-Cuadra, E. Zohar, and J. I. Cirac, New J. Phys. 19, 063038 (2017), arXiv:1702.05492 [quant-ph] .
  * Zhang  _et al._ (2018) J. Zhang, J. Unmuth-Yockey, J. Zeiher, A. Bazavov, S. W. Tsai, and Y. Meurice, Phys. Rev. Lett. 121, 223201 (2018), arXiv:1803.11166 [hep-lat] .
  * Chanda  _et al._ (2022) T. Chanda, M. Lewenstein, J. Zakrzewski, and L. Tagliacozzo, Phys. Rev. Lett. 128, 090601 (2022), arXiv:2107.11656 [hep-th] .
  * Šibalić and Adams (2018) N. Šibalić and C. S. Adams,  _Rydberg Physics_, 2399-2891 (IOP Publishing, 2018).
  * (22) “Amazon braket,” https://aws.amazon.com/braket.
  * Pisarski  _et al._ (2019) R. D. Pisarski, V. V. Skokov, and A. M. Tsvelik, Universe 5, 48 (2019), arXiv:2202.01036 [hep-ph] .
  * Kojo  _et al._ (2010) T. Kojo, Y. Hidaka, L. McLerran, and R. D. Pisarski, Nucl. Phys. A 843, 37 (2010), arXiv:0912.3800 [hep-ph] .
  * Basar and Dunne (2008) G. Basar and G. V. Dunne, Phys. Rev. Lett. 100, 200404 (2008), arXiv:0803.1501 [hep-th] .
  * Gustafson (2021) E. Gustafson, Phys. Rev. D 103, 114505 (2021), arXiv:2104.10136 [quant-ph] .
  * Gustafson (2022) E. Gustafson, (2022), arXiv:2201.04546 [quant-ph] .
  * Carena  _et al._ (2022a) M. Carena, H. Lamm, Y.-Y. Li, and W. Liu, Phys. Rev. Lett. 129, 051601 (2022a), arXiv:2203.02823 [hep-lat] .
  * Carena  _et al._ (2022b) M. Carena, E. J. Gustafson, H. Lamm, Y.-Y. Li, and W. Liu, Phys. Rev. D 106, 114504 (2022b), arXiv:2208.10417 [hep-lat] .
  * Unmuth-Yockey  _et al._ (2018) J. Unmuth-Yockey, J. Zhang, A. Bazavov, Y. Meurice, and S.-W. Tsai, Phys. Rev. D 98, 094511 (2018).
  * Zhang  _et al._ (2021) J. Zhang, Y. Meurice, and S.-W. Tsai, Phys. Rev. B 103, 245137 (2021).

*

##  Appendix A Compact scalar QED

Matching the matrix elements to the target model
$H_{eff}=\hat{H}_{F}+Const.\hat{I}$ ($Const.$ is a constant), we have the
following equations for OBC

$\displaystyle Const.$ $\displaystyle=$
$\displaystyle-\Delta-\Delta_{0}+\frac{\Omega^{2}}{4}\left(\frac{2}{\Delta-
V_{0}}-\frac{1}{\Delta}\right)$ $\displaystyle\frac{U}{2}$ $\displaystyle=$
$\displaystyle\Delta_{0}+\frac{\Omega^{2}V_{0}}{4\Delta(V_{0}-\Delta)}$
$\displaystyle X$ $\displaystyle=$
$\displaystyle\frac{\Omega^{2}V_{0}}{2\Delta(V_{0}-\Delta)}$ (29)

For 00BC,

$\displaystyle Const.$ $\displaystyle=$
$\displaystyle-\Delta-\Delta_{0}+2V_{1}+\frac{\Omega^{2}}{4}\left(\frac{2}{\Delta-
V_{0}}-\frac{1}{\Delta}\right)$ $\displaystyle\frac{U}{2}+Y$ $\displaystyle=$
$\displaystyle\Delta_{0}+2(V_{2}-V_{1})+\frac{\Omega^{2}V_{0}}{4\Delta(V_{0}-\Delta)}$
$\displaystyle X$ $\displaystyle=$
$\displaystyle\frac{\Omega^{2}V_{0}}{2\Delta(V_{0}-\Delta)}$ (30)

We can perform the same procedure for two rungs. The configurations are shown
in Fig. 3. We obtain the following equations for 00BC

$\displaystyle Const.$ $\displaystyle=$
$\displaystyle-(\Delta+\Delta_{0})N_{s}+(N_{s}+1)V_{1}$ (31)
$\displaystyle+\frac{\Omega^{2}}{4}\left(\frac{2}{\Delta-
V_{0}}-\frac{1}{\Delta}\right)N_{s}$ $\displaystyle U$ $\displaystyle=$
$\displaystyle
2\Delta_{0}+2V_{3}-2V_{1}+\frac{\Omega^{2}V_{0}}{2\Delta(V_{0}-\Delta)}$
$\displaystyle Y$ $\displaystyle=$ $\displaystyle 2V_{2}-V_{1}-V_{3}$
$\displaystyle 2(Y+Y^{\prime})$ $\displaystyle=$ $\displaystyle V_{1}-V_{3}$
$\displaystyle X$ $\displaystyle=$
$\displaystyle\frac{\Omega^{2}V_{0}}{2\Delta(V_{0}-\Delta)}$ (32)

Notice that it is impossible to match the boundary terms for OBC, which
requires $V_{1}=V_{3}$. Now we can tune the parameters in Rydberg systems to
obtain the target model Eq. (35). For example, the target model with
$U=0,X=0.02,Y/X=0.35,Y^{\prime}/X=12.627$ corresponds to
$V_{0}/\Delta=2,V_{0}/\Omega=100,\Delta_{0}/\Omega=0.509,\rho=a_{y}/a_{x}=0.431$
in the Rydberg simulator.

We have freedom to change the ratio $V_{0}/\Delta$ to match the same point.
The lowest two energy densities are ($-0.01017156$, $-0.00357107$, $N_{s}=2$),
($-0.00995229$, $-0.00619866$, $N_{s}=3$), ($-0.00983596$, $-0.00730903$,
$N_{s}=4$) for the target model and ($-0.00998068$, $-0.00338451$, $N_{s}=2$),
($-0.00976367$, $-0.00599897$, $N_{s}=3$), ($-0.009649$, $-0.00710675$,
$N_{s}=4$) for the Rydberg simulator. The errors are all of order $10^{-4}$,
consistent with the perturbation theory.

Matching to the target Hamiltonian (35), $U=-2\Delta+2V_{2}$, $Y=-V_{2}$,
$Y^{\prime}=(V_{1}+V_{2})/2$. We can only simulate negative $Y$ using this
two-leg Rydberg ladder. The errors of the effective Hamiltonian are from the
transition from the spin-1 sector to the $|rr\rangle$ state, which if of order
$\Omega^{2}/(4V_{0})$, and from the next-nearest-neighbor interactions $\sim
V_{1}/64$.

The model is the Lattice scalar QED (sQED) on a $N_{s}\times N_{\tau}$
Euclidean spacetime lattice. There are two more columns on two sides, the
$0$th and the ($N_{s}+1$)th columns, controlling the boundary conditions. We
use $00$ boundaray condition ($00$BC), where the plaquette (field) quantum
numbers on the $0$th and the ($N_{s}+1$)th columns are zero thus the total
charge in the system is zero. In the time continuum limit Bazavov  _et al._
(2015); Unmuth-Yockey  _et al._ (2018), we obtain the Hamiltonian in the
charge representation

$\displaystyle\hat{H}_{C}$ $\displaystyle=\frac{U}{2}\sum_{1\leq j,k\leq
N_{s}}c_{jk}S^{z}_{j}S^{z}_{k}+\frac{Y}{2}\sum_{i=1}^{N_{s}+1}(S_{i}^{z})^{2}$
(33)
$\displaystyle-\frac{X}{2}\sum_{i=1}^{N_{s}}(\hat{U}_{i}^{+}\hat{U}_{i+1}^{-}+\hat{U}_{i}^{-}\hat{U}_{i+1}^{+})\
,\sum_{i}S_{i}^{z}=0$

where $c_{jk}=N_{s}+1-\max\\{j,k\\}$. The eigenvalues of the operator $S^{z}$
($S^{z}|n\rangle=n|n\rangle$) are integer charge $n=0,\pm 1,\pm 2,\ldots$
attached on the vertical links, and the eigenstates define basis of the charge
representation. The operator $\hat{U}^{\pm}$ raises (lowers) the charge of a
state by one $\hat{U}^{\pm}|n\rangle=|n\pm 1\rangle$. When the gauge coupling
$U=0$, the spin-$1$ truncation has an infinite-order phase transition from a
gapped phase into a BKT critical line Zhang  _et al._ (2021). For nonzero
gauge coupling, Eq. (33) has unusual long-range interactions, thus it is
difficult to design quantum simulators for it. Using Gauss’s Law, we can go to
the field representation

$\displaystyle\hat{H}_{F}$
$\displaystyle=\frac{U}{2}\sum_{p=1}^{N_{s}}\left(L^{z}_{p}\right)^{2}+\frac{Y}{2}\sum_{p=1}^{N_{s}+1}(L^{z}_{p}-L^{z}_{p-1})^{2}$
(34) $\displaystyle-\frac{X}{2}\sum_{p=1}^{N_{s}}(U^{+}_{p}+U^{-}_{p})\ ,$

where the eigenvalues of $L^{z}$ is the field quantum numbers $m=0,\pm 1,\pm
2,\ldots$ attached to the plaquettes. Thera are $N_{s}$ plaquettes and
$N_{s}+1$ links, $L^{z}_{0}=L^{z}_{N_{s}+1}=0$ under 00BC. For open boundary
conditions (OBCs), $p$ is taken from $2$ to $N_{s}$, the Hamiltonian contains
multiple charge sectors that have total charge $L^{z}_{N_{s}}-L^{z}_{1}$.

The field representation is identical to the charge representation without
applying a truncation. The two Hamiltonians are quite different for small spin
truncations. For one plaquette with two links, the basis in the field
representation is $|-S\rangle$, $|-S+1\rangle$, …, $|S\rangle$. The basis in
the charge representation is $|-S,S\rangle$, $|-S+1,S-1\rangle$, …,
$|S,-S\rangle$. The two Hilbert spaces observe one-to-one mapping. But for two
plaquettes with three links, the dimensions of the two Hilbert spaces are not
the same, e.g. for $S=1$, the basis in the field representation is
$|-1,-1\rangle$, $|-1,0\rangle$, $|-1,1\rangle$, $|0,-1\rangle$,
$|0,0\rangle$, $|0,1\rangle$, $|1,-1\rangle$, $|1,0\rangle$, $|1,1\rangle$,
the corresponding states that satisfy Gauss’s law in the charge representation
are $|-1,0,1\rangle$, $|-1,1,0\rangle$, $|-1,2,-1\rangle$, $|0,-1,1\rangle$,
$|0,0,0\rangle$, $|0,1,-1\rangle$, $|1,-2,1\rangle$, $|1,-1,0\rangle$,
$|1,0,-1\rangle$. The two states $|-1,2,-1\rangle$, $|1,-2,1\rangle$ are
truncated in the charge representation with spin-1 truncation. Generally, for
$N_{s}$ plaquettes with $N_{s}+1$ links, if the states
$|m_{1},m_{2},...,m_{L}\rangle$ with $|m_{p}-m_{p+1}|>S$ in the field
representation have large energy gap to other states, the effective low-energy
basis is a subset of basis in the charge representation with the same spin
truncation.

For spin-1 truncation, we can add a high energy penalty for states
$|\ldots,\pm 1,\mp 1,\ldots\rangle$. We have a new Hamiltonian

$\displaystyle\hat{H}_{F}$
$\displaystyle=\frac{U}{2}\sum_{p=1}^{N_{s}}\left(L^{z}_{p}\right)^{2}+\frac{Y}{2}\sum_{p=1}^{N_{s}+1}(L^{z}_{p}-L^{z}_{p-1})^{2}$
(35)
$\displaystyle-\frac{Y^{\prime}}{2}\sum_{p=1}^{N_{s}+1}L^{z}_{p}L^{z}_{p-1}(L^{z}_{p}-L^{z}_{p-1})^{2}-\frac{X}{2}\sum_{p=1}^{N_{s}}(U^{+}_{p}+U^{-}_{p})$
$\displaystyle=(\frac{U}{2}+Y)\sum_{p=1}^{N_{s}}\left(L^{z}_{p}\right)^{2}-(Y+Y^{\prime})\sum_{p=1}^{N_{s}+1}L^{z}_{p}L^{z}_{p-1}$
(36)
$\displaystyle+Y^{\prime}\sum_{p=1}^{N_{s}+1}(L^{z}_{p})^{2}(L^{z}_{p-1})^{2}-\frac{X}{2}\sum_{p=1}^{N_{s}}(U^{+}_{p}+U^{-}_{p}),$

where the $Y^{\prime}$ term is only nonzero for states $|\ldots,\pm 1,\mp
1,\ldots\rangle$. For OBC, the coefficients of $(\hat{L}^{z}_{p})^{2}$ at
$p=1,N_{s}$ are $U/2+Y/2$ instead of $U/2+Y$. The energy gap in the field
representation will be closer to that in the charge representation by setting
$Y^{\prime}\gg 1$. The states in the charge representation like
$|\ldots,111,\ldots\rangle$ are not in the basis of the field representation,
so the field representation with spin-1 truncation has no infinite-order phase
transitions, either. But it has other interesting critical points and we are
interested in quantum simulation of Eq. (35) using Rydberg ladders. The
raising and lowering operators $U^{\pm}$ can be replaced by clock matrices
$C^{\pm}$ which have additional element connecting $-1$ and $1$. The same idea
can be applied to higher spin truncation to approach the true sQED.

For 00BC,
$\hat{H}_{\rm{eff}}\rightarrow\hat{H}_{\rm{eff}}+(V_{2}-V_{1})[(L^{z}_{1})^{2}+(L^{z}_{N_{s}})^{2}]+2V_{1}\hat{I}$.
Notice that one can tune the first and the last $\Delta_{0}$ locally to
realize 00BC, so we do not really need to add two more Rydberg atoms. This
method can be generalized to arbitrary spin truncation, where $2S+1$ legs are
needed for spin-$S$ truncation. But there is a simpler configuration to
realize spin-1 truncation: the two-leg ladder.

The equations relating the parameters for 00BC are

$\displaystyle Const.$ $\displaystyle=$
$\displaystyle-(\Delta+\Delta_{0})N_{s}+(N_{s}+1)V_{1}+\frac{\Omega^{2}}{4}\left(\frac{2}{\Delta-
V_{0}}-\frac{1}{\Delta}\right)N_{s}$ $\displaystyle U$ $\displaystyle=$
$\displaystyle 2\Delta_{0}+2V_{2}-2V_{1}$ $\displaystyle Y$ $\displaystyle=$
$\displaystyle V_{2}-V_{1}$ $\displaystyle Y^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}Y$ $\displaystyle X$ $\displaystyle=$
$\displaystyle\frac{\Omega^{2}V_{0}}{2\Delta(V_{0}-\Delta)}.$ (37)

In this case, the Rydberg simulator can only simulate the target model with
negative $Y$ and fixed ratio $Y^{\prime}/Y=-3/2$.

$\displaystyle\hat{H}_{Fc}$
$\displaystyle=(\frac{U}{2}+Y)\sum_{p=1}^{N_{s}}\left(L^{z}_{p}\right)^{2}+\frac{Y}{2}\sum_{p=1}^{N_{s}+1}L^{z}_{p}L^{z}_{p-1}$
$\displaystyle-\frac{3Y}{2}\sum_{p=1}^{N_{s}+1}(L^{z}_{p})^{2}(L^{z}_{p-1})^{2}-\frac{X}{2}\sum_{p=1}^{N_{s}}(\hat{C}^{+}_{p}+\hat{C}^{-}_{p}),$
(38)