# Extending generalized Debye analysis to long timescale magnetic relaxation

Jeremy D. Hilgar
Department of Chemistry and Biochemistry
University of California, San Diego
La Jolla, California
&Aaron K. Butts
Department of Chemistry and Biochemistry
University of California, San Diego
La Jolla, California
&Jeffrey D. Rinehart*
Department of Chemistry and Biochemistry
University of California, San Diego
La Jolla, California
<EMAIL_ADDRESS>

###### Abstract

As the ability to generate magnetic anisotropy in molecular materials
continues to hit new milestones, concerted effort has shifted towards
understanding, and potentially controlling, the mechanisms of magnetic
relaxation across a large time and temperature space. Slow magnetic relaxation
in molecules is highly temperature-, field-, and environment-dependent with
the relevant timescale easily traversing ten orders of magnitude for current
single-molecule magnets (SMM). The prospect of synthetic control over the
nature of (and transition probabilities between) magnetic states make
unraveling the underlying mechanisms an important yet daunting challenge.
Currently, instrumental considerations dictate that the characteristic
relaxation time, $\tau$, is determined by separate methods depending on the
timescale of interest. Static and dynamic probe fields are used for long and
short timescales, respectively. Each method captures a distinct, non-
overlapping time range, and experimental differences lead to the possibility
of fundamentally different meanings for $\tau$ being plotted and fitted
globally as a function of temperature. Herein, we present a method to generate
long-timescale waveforms with standard vibrating sample magnetometry (VSM)
instrumentation, allowing extension of alternating current (AC) magnetic
impedance measurements to SMMs and other superparamagnets with arbitrarily
long relaxation time.

_K_ eywords single-molecule magnetism $\cdot$ superparamagnetic relaxation
$\cdot$ generalized Debye analysis $\cdot$ Cole-Cole Plot $\cdot$ lanthanide
magnetism

## Introduction

Single-molecule magnets (SMMs) are a class of zero-dimensional materials that
utilize the spin-orbit interaction to generate an axially anisotropic angular
momentum. Due to a preferred orientation along the anisotropy axis, the
magnetization vector exhibits slow relaxation, as if it were a superparamagnet
with a classic double well potential replaced by discrete quantum
states.sessoliHighspinMoleculesMn12O121993;
sessoliMagneticBistabilityMetalion1993 As a form of superparamagnetism,
single-molecule magnetism cannot result in magnetic ordering. In recent years,
however, the timescale of magnetic relaxation has extended drastically such
that single magnetic ions are now shown to retain magnetization information
for timescales exceeding 1 s in the liquid nitrogen temperature regime. At
temperatures well below the superparamagnetic blocking temperature, TB, these
highly anisotropic systems can exhibit "frozen" dynamics due to their
extremely long relaxation times, allowing the performance of typical permanent
magnet characterization experiments, such as magnetization vs. field
measurements that display hysteretic behavior despite the lack of a permanent
magnetic ground state. Magnetization vs. field scans contain information about
the field dependence of quantum tunneling as well as the maximal magnetization
and coercivity of the material. They do not, however, contain quantitative
information about the number and time constants of the various relaxation
processes contributing to the magnetic relaxation. This time-domain data is
readily obtained as phase-dependent AC magnetic susceptibilities between
approximately 0.001 and 1 s (AC relaxation), yet its collection becomes
problematic at longer timescales. The standard method for determining long
timescale relaxation time constants (DC relaxation) is to apply a large static
magnetic field to magnetize a macroscopic amount of microcrystalline SMM
material, remove the field as quickly as possible, and track the decay towards
zero magnetization as a function of time. These data are a valuable source of
information, yet have several drawbacks that can obscure the nature of the
isolated molecular zero-field relaxation. First, there is a significant gap in
measurement timescale between that captured by AC and DC relaxation (generally
1–10 min., depending on instrumentation capability). Second, commercial
instrumentation is generally only capable of setting zero field within $\pm$
30 Oe, leading to an inherent decay offset. In very sensitive systems, this
small field offset can drastically alter the time constant or even completely
change the dominant magnetic relaxation mechanism. Third, the scan from high
field to zero field can instigate the transfer of energy between the
(ostensibly isolated) molecules leading to propagating, multi-molecular
relaxation processes. This last point especially contributes to multi- and
stretched exponential dynamics that can obscure the molecular zero-field
dynamics of interest for comparison to AC susceptibility relaxation data.
Given the influx of new long-timescale SMMs and the need for analysis tools
that leverage existing instrumentation, we sought to devise a method to extend
the dynamic magnetization measurement methodology of AC susceptibility to
longer timescales. Importantly, this design required compatibility with
commonly-employed research magnetometers, similarity to the input signal
characteristics of commonly used AC waveforms for probing magnetic relaxation,
and extendability of the measurement to arbitrarily-long relaxation regimes.
Herein we present a technique fitting these criteria whereby phase-dependent
magnetic susceptibilities are extracted from standard vibrating sample
magnetometer (VSM) measurements. To provide an initial demonstration of the
method’s efficacy, we examine the relaxation behavior of a variant of the SMM
"erbocene" sandwich motif, [K(18-c-6)][Er(hdcCOT)2] (1, hdcCOT =
hexahydrodicyclopentacyclooctatetraenide dianion, 18-c-6 = 18-crown-6) that
has not been previously magnetically characterized.

## Results and Discussion

### Synthesis and solid-state structure

Preparation of compound 1 was carried out under air-free conditions using a
procedure similar to that employed for its previously reported lithio analog,
[Li(THF)(DME)][Er(hdcCOT)2].hillerParamagneticNMRAnalysis2016a Neutral hdcCOT
was synthesized from 1,6-heptadiyne via a Ni0-catalyzed dimerization
reactionwenderNickelCatalyzedCycloadditions2007b and was subsequently reduced
using potassium graphite (KC8). The dipotassium salt of the alkyl-substituted
COT dianion (K2hdcCOT) reacted readily with a suspension of erbium trichloride
in THF to give K[Er(hdcCOT)2]; $\delta$H(500 MHz; THF-d8) –19.3 (br), –34.4
(br), –43.5 (br), –53.5 (br), –61.9 (br), –213.3 (br), –223.9 (br), –269.3
(br), and –343.3 (br) (Figures S1–S3). Removal of KCl and addition of 18-c-6
(to facilitate crystallization) yielded a solution which, when layered with
pentane through vapor diffusion, deposited bright-yellow crystalline rods of
1. X-ray analysis of these crystals showed that the solid-state structure of 1
includes an Er3+ ion in a homoleptic coordination environment with two
hdcCOT2- ligands coordinated in a 90°staggered conformation (Figure 1). Steric
bulk is presumably responsible for the observed conformation about the metal
center, however this bulk imparts no appreciable elongation of the average Er-
hdcCOTcentroid distance relative to an unsubstituted analog
(CCDCgroom_cambridge_2016 Identifier YIWTUVmeihausMagneticBlocking102013;
ungurFinetuningLocalSymmetry2014). The potassium ion in 1 supports a
$\kappa{}^{6}$ interaction to 18-c-6 as well as an $\eta{}^{4}$ interaction to
one face of a hdcCOT2- ring and the latter interaction lowers the
hdcCOT–Er–hdcCOT sandwich angle from linearity to 178°.

Figure 1: Solid-state structure of 1 with spheres representing Er (pink), K
(purple), O (red), and C (gray). Hydrogen atoms and lattice THF molecules have
been omitted for clarity.

### Ab initio electronic structure

Previous studies on [Er(COT)2]– meihausMagneticBlocking102013;
ungurFinetuningLocalSymmetry2014; leroySandwichComplexAxial2014;
leroyCouplingStrategiesEnhance2014 and
[Er(COT)]+jiangOrganometallicSingleIonMagnet2011;
mengBoratabenzeneCyclooctatetraenylLanthanide2016;
hilgarFerromagneticCouplingChloridebridged2017; hilgarMetalLigandPair2018;
chenSoftPhosphorusAtom2018; heEnantiopureBenzamidinateCyclooctatetraene2018;
hilgarMillionfoldRelaxationTime2019 motifs have demonstrated that a remarkably
consistent single-ion anisotropy is generated at Er3+ when suitably equatorial
coordination environments are selected. The efficacy by which the planar COT
dianion in particular can stabilize prolate, high-moment $m_{J}$ states on
Er3+ has been well-rationalized by varying levels of theory and herein we
modelled magnetic properties of 1 using complete active space, self-consistent
field (CASSCF) methods (Figure S4). Decomposition of the calculated ground
state spin-orbit wavefunctions into an $m_{J}$ basis reveals that the lowered
symmetry around the Er3+ center in 1 relative to unsubstituted [Er(COT)2]–
species plays a negligible role toward wavefunction mixing and consequently
the ground Kramers doublet is highly axial with 97.7% $m_{J}=\pm{}15/2$
character. The first two excited doublets (94.7% $m_{J}=\pm{}13/2$ and 93.1%
$m_{J}=\pm{}1/2$) are predicted to lie 169.2 and 211.3 cm–1 above the ground
doublet. Ground to first and second excited Kramers doublet splittings have
previously been calculated for
[K(18-c-6)][Er(COT)2]ungurFinetuningLocalSymmetry2014; interestingly, we
observe no appreciable difference in calculated splittings between 1 and its
unsubstituted analog and we anticipate alkyl-substitutions on COT play a
marginal role in affecting low-energy spin-orbit states on erbium. Transverse
magnetic moment matrix elements connecting states within the $J=15/2$ manifold
have also been calculated. As discussed
elsewhere,ungurInterplayStronglyAnisotropic2013 the square of these matrix
elements are roughly proportional to transition rates between the states that
they connect and the effective barrier of magnetic reversal ($U_{eff}$) in the
high-temperature limit can be approximated as the shortest path between ground
states with non-negligible matrix elements. Calculated matrix elements between
the ground doublet and first excited doublet are small (3.3$\times$10–6 and
5.1$\times$10–4 $\mu_{B}$, respectively) which indicate that QTM between these
states should be suppressed. Matrix elements connecting the second excited
doublets ($m_{J}\sim\pm{}1/2$) and cross-terms between the first and second
excited doublets (e.g. $m_{J}\sim+13/2$ and $m_{J}\sim-1/2$) are larger (3.2
and 4.1$\times$10–2 $\mu_{B}$, respectively) and likely relaxation pathways in
the high temperature limit will involve thermally assisted quantum tunneling
of the magnetization (QTM) or Orbach relaxation modes through the second
excited Kramers doublet.

### Magnetic characterization

#### Static magnetism

Zero-field cooled static magnetic susceptibilities for 1 were measured over
2–300 K with a 1000 Oe bias field and the data were plotted as the molar
susceptibility times temperature product ($\chi{}T$) vs. temperature (Figure
S5). At 300 K, $\chi{}T=11.69$ emu mol-1 K, which is close to the expected
value for an Er3+ ensemble with equal populations across the $J=15/2$ spin-
orbit manifold ($\chi{}T=11.49$ emu mol-1 K). As $T$ is lowered from 300 K,
$\chi{}T$ shows a small decline as this manifold becomes thermally
depopulated. At 12 K, $\chi{}T=11.24$ emu mol-1 K; below this temperature the
data display a precipitous drop to 1.60 emu mol-1 K at 2 K. This drop is
indicative of magnetic blocking on the timescale of the DC scan and thus
isothermal magnetization measurements ($H_{max}=|7|$ T, $~{}10$ Oe sec–1 ramp
rate) were collected to further probe for superparamagnetic behavior. At 2 K,
the magnetization saturates at 4.98 $\mu_{B}$ mol–1 when the field is swept
above 3.5 T. Consistent with the ZFC results, hysteretic behavior is observed
when the field is swept back from 7 T. Near 0 T, the magnetization drops
abruptly to 1.25 $\mu_{B}$ mol–1 and the application of negative fields yields
a coercivity $H_{C}$ = 1.1 T on this timescale. Zero-field magnetization loss
and waist-restricted hysteresis are commonly observed in SMMs when transverse
crystal field (CF) components enable fast relaxation via QTM between ground
states.liuSymmetryStrategiesHigh2018 Symmetry optimization strategies can be
employed to mitigate this effect and in the case of [Er(COT)2]– species the
psuedo-C$\infty$ axis minimizes transverse CF components to yield a ground
doublet within which QTM is largely quenched. The striking zero-field
magnetization loss in 1 mimics well the behavior observed in
[K(18-c-6)][Er(COT)2] which, through a magnetic dilution study with yttrium,
was shown to be facilitated by a bulk magnetic avalanche
effect.meihausMagneticBlocking102013

#### Dynamic magnetism

Consistent with SMM behavior, the AC susceptibility phase-shift of 1 shows a
clear dependence on drive-field frequency below 24 K and short time-scale
relaxation times were extracted from these data by simultaneously fitting the
in-phase ($\chi{}\prime$) and out-of-phase ($\chi{}\prime{}\prime{}$) signals
($H_{DC}$ = 0 Oe, $f$ = 0.1–1000 Hz, $T$ = 12–24 K) to a generalized Debye
relaxation model (Figure 2, Equations 1 and 2).
gatteschiMolecularNanomagnets2006

$\chi{}\prime{}(\omega{})=\chi{}_{S}+(\chi{}_{T}-\chi{}_{S})\frac{1+(\omega{}\tau{})^{1-\alpha{}}sin(\pi{}\alpha{}/2)}{1+2(\omega{}\tau{})^{1-\alpha{}}sin(\pi{}\alpha{}/2)+(\omega{}\tau{})^{2-2\alpha{}}}$
(1)
$\chi{}\prime{}\prime{}(\omega{})=(\chi{}_{T}-\chi{}_{S})\frac{(\omega{}\tau{})^{1-\alpha{}}cos(\pi{}\alpha{}/2)}{1+2(\omega{}\tau{})^{1-\alpha{}}sin(\pi{}\alpha{}/2)+(\omega{}\tau{})^{2-2\alpha{}}}$
(2)

Figure 2: Cole-Cole plot of 1 between 2–24 K. Colored points are
susceptibilities measured via standard AC measurements (circles) and extracted
from Fourier analysis of VSM data (triangles). Black lines represent fits a
generalized Debye model (Equations 1 & 2).

Fitted eccentricities are low ($\alpha{}_{max}=0.25(1)$) indicating a single
relaxation time with a narrow distribution is associated with each temperature
across the measured frequency range. A plot of $ln(\tau{})$ vs. $1/T$ (Figure
5) shows that an over-barrier Orbach relaxation mechanism is operant in this
temperature range and least-squares fitting to an Arrhenius Law yielded an
effective barrier height $U_{eff}$ = 147.7(7) cm–1 ($\tau{}_{0}$ =
3.1(1)$\times$10–8 s). A previous analysis on [K(18-c-6)][Er(COT)2] yielded an
effective barrier $U_{eff}$ = 147 cm–1; this striking similarity further
corroborates our ab initio results which suggested that the alkyl
substitutions on the COT ring play a negligible role toward modifying the
crystal field environment around Er3+.

Preliminary studies on relaxation dynamics longer than the MPMS3 AC timescale
($<0.1$ Hz) were conducted between 2–10 K by measuring the zero-field DC
relaxation profile of a sample previously in equilibrium with a 7 T field.
Consistent with the presence of a small remnant field in the MPMS3
superconducting magnet, final moments are negative and their magnitudes range
from 0.37–1.8% relative to the moment at $t$ = 0 s. This residual field is
typically on the order of 20–30 Oe and two approaches used to minimize it
include oscillating the field to zero and resetting the magnet (see e.g. "MPMS
Application Note 1014-208," Quantum Design, 2001). We note that neither
approach is feasible during the DC relaxation experiment. Given that multiple
distinct processes or a distribution of processes are typically responsible
for SMM relaxation at a given temperature, several forms of exponential decay
models are routinely used to extract relaxation times in the long time-scale
regime. Herein we fit DC relaxation data of 1 (Figure 3, Table 1) using a
stretched exponential (Equation 3, 0 < $\beta{}$ < 1) and a bi-exponential
model (Equation 4).

$M(t)=M_{f}+(M_{0}-M_{f})exp\left(\frac{-t}{\tau_{str}}\right)^{\beta}$ (3)
$M(t)=M_{f}+(M_{0}-M_{f})\left[A_{1,2}exp\left(\frac{-t}{\tau_{1}}\right)+(1-A_{1,2})exp\left(\frac{-t}{\tau_{2}}\right)\right]$
(4)

Figure 3: DC relaxation data of 1 between 2–10 K. Colored points are moments
measured from near-zero-field DC scans. Black lines are fits to (a) stretched
and (b) bi-exponential models.

Corroborating with the observation of multiple potential relaxation mechanisms
in the magnetization data, we found that a simple exponential equation
(Equation 3, $\beta{}$ = 1) does not adequately model the relaxation behavior
of 1 below 10 K. Disagreement between data and fits become more prominent as
temperature is lowered and it is evident that the extracted relaxation time is
biased toward the short time-scale. Agreement between the data and fit can be
improved with the addition of a stretching parameter, $\beta{}$ to the
exponential. Similar to the $\alpha{}$ parameter of the generalized Debye
equation, $\beta{}$ accounts for a distribution about an average relaxation
process. Incorporation of $\beta{}$ results in the extraction of relaxation
times with wide distributions at low T ($\beta{}$2 K = 0.505(9)). A similar
quality fit was obtained when two exponential relaxation processes were used
to model the data. Interestingly, this model reveals that relaxation of 1 can
be rationalized by two processes operating on significantly different
timescales (i.e. separated by an order of magnitude at 2 K) with the long
relaxation time $\tau_{1}$ and short relaxation time $\tau_{2}$ bracketing the
relaxation times obtained previously with a stretched exponential model. The
presence of a small residual field (vide supra) during the DC relaxation
measurement provides a potential explanation for the observed long timescale
process. This field can lift Kramers degeneracy between bistable orientations
of the ground state and concomitantly hinder QTM. The addition of further
constants and weighting factors would undoubtedly improve the fit, yet the
ambiguity of the existing parameters and their importance to the intrinsic
zero-field molecular relaxation make further parameterization of dubious
value.

Table 1: Long time-scale relaxation times of 1 extracted from various models.

| Equations 1 & 2 | Equation 3 | Equation 4
---|---|---|---
T (K) |  $\tau_{VSM}$ (s) | $\alpha$ |  $\tau_{str}$ (s) | $\beta$ |  $\tau_{1}$ (s) |  $\tau_{2}$ (s) | $A_{1,2}$
$2$ | $54(4)$ | $0.18(4)$ | $145(3)$ | $0.505(9)$ | $1093(29)$ | $87(1)$ | $0.222(4)$
$4$ | $43(3)$ | $0.18(4)$ | $119(3)$ | $0.58(1)$ | $621(18)$ | $69(1)$ | $0.252(6)$
$6$ | $37(3)$ | $0.17(4)$ | $106(2)$ | $0.62(1)$ | $442(9)$ | $61.8(8)$ | $0.271(5)$
$8$ | $33(2)$ | $0.15(5)$ | $81(1)$ | $0.74(1)$ | $241(5)$ | $52.7(6)$ | $0.271(7)$
$10$ | $13.4(8)$ | $0.10(3)$ | $36.3(4)$ | $0.88(1)$ | $122(8)$ | $30.6(4)$ | $0.13(1)$

Parentheses following fitted values represent 95% confidence intervals.

#### VSM waveform analysis

Prompted by the desire for more quantitative information about long-timescale
relaxation, we sought a method that could more faithfully model the zero-field
relaxation dynamics of 1. Ideally the method would be well-suited for
extracting $\tau{}$ information from non-ensemble magnetic processes and would
furthermore be capable of directly probing for the existence of multiple such
relaxation processes. Naturally, Debye model fitting of complex AC magnetic
susceptibilities fit these criteria well and is the standard when relaxation
times fall within common magnet modulation coil limits ($\nu$ = 0.1–1500 Hz).
An interesting question that follows from the preceding is whether complex
susceptibilities can be reliably extracted from data collected from waveforms
constructed in either DC or VSM scan modes. Although magnetic field controls
have a limited linear sweep rate ($R_{H_{max}}$ = 1.6 Oe sec–1 for the Quantum
Design MPMS3 magnetometer), drive fields of arbitrarily low frequency can be
generated if an appropriate waveform is selected. To test whether phase
information could be extracted from 1 outside the typical AC scan range, we
developed a sequence which measures VSM moments generated by a square-wave
magnetic drive field oscillating about 0 Oe ($H_{max}$ = 8 Oe, $R_{H}$ = 30 Oe
sec–1). Our initial studies on this approach yielded profiles of average
moment significantly different from 0 emu mol–1, indicating a persistent
offset from zero bias field. A single magnet reset operation before data
collection was sufficient to remove this residual field, allowing measurement
of periodic moments with no appreciable offset (Figure 4a–c). To optimize
sequence run times, magnetic field cycle counts were scaled with frequency
(e.g. 20 cycles were measured at 0.0264 Hz and 1 cycle was measured at 0.00013
Hz). An example truncated MPMS3 sequence file can be found in Listing 1 in the
SI. The discrete Fourier transform was applied to $H$ vs. $t$ and $M$ vs. $t$
data which yielded complex-valued frequency space representations. Plots of
the absolute value of the transform vs. frequency revealed spectra with a
primary peak at the fundamental drive field frequency followed by overtone
peaks of diminishing intensity. We set
$\chi=|FT_{D}(M_{max})|/|FT_{D}(F_{max})$, where $|FT_{D}(M_{max})|$ and
$|FT_{D}(F_{max})|$ are absolute values at the fundamental frequency of the
complex moment and field spectra, respectively. The in-phase ($\chi{}\prime$)
and out-of-phase ($\chi{}\prime\prime$) components of the magnetic
susceptibility were then calculated as $\chi{}\prime=\chi{}cos(\phi{})$ and
$\chi{}\prime\prime=\chi{}sin(\phi{})$, with $\phi{}$ being the phase angle
between the field and moment spectra at the fundamental frequency. Using this
method, complex susceptibilities were extracted from VSM moment waveforms of 1
collected at temperatures between 2–10 K using square-wave drive fields with
periods between 38–7690 s.

Figure 4: Representative raw VSM waveform data (top row) and discrete Fourier
transform of VSM waveform data (bottom row) with period lengths 7688 s (a and
a*), 968 s (b and b*), and 128 s (c and c*). Black lines are field values and
purple lines are measured moments.

Susceptibility parameters extracted from Fourier analysis of the VSM waveform
data (Figure 2) show a clear phase dependence on drive field frequency. In an
approach analogous to that taken with standard AC susceptibilities, we fit in-
phase and out-of-phase components of this data to a generalized Debye model
(Equations 1 & 2, Table 1) and obtained $\alpha{}$ values close to 0
indicating that only a single characteristic relaxation time exists at each
temperature. An Arrhenius plot of fitted relaxation times (Figure 5) shows
that $\tau{}$ is weakly temperature dependent between 2–8 K and at 10 K
$\tau{}$ begins to drop considerably as it approaches the high-temperature
Orbach relaxation regime. We note that the relaxation time measured at 10 K
via this method nearly matches the relaxation time predicted by an Arrhenius
fit to the AC susceptibility data. This provides a strong indication that
these methods are complementary probes of molecular magnetic relaxation for 1.
A comparison between relaxation times extracted from VSM waveforms and DC
relaxation measurements reveal that each fitting method yields relaxation time
profiles that deviate sublinearly from Arrhenius behavior. However, $\tau{}$
values obtained from exponential fits of DC data across the measured
temperature range are consistently higher than the corresponding VSM waveform
measurements. If we assume the Debye analysis to provide the most accurate
$\tau{}$, DC relaxation times extracted with Equation 3 ($\tau_{str}$) and
Equation 4 ($\tau_{2}$) represent average errors of 270% and 175%,
respectively. The long component of Equation 4, $\tau{}_{1}$ represents an
even more extreme error of 1261%.

Figure 5: Standard AC relaxation times (circles), VSM waveform Debye fit
(triangles), stretched exponential fit (squares), and bi-exponential fit
(diamonds) of magnetic relaxation data for 1 with fit to relaxation model S1
(black line, parameters in Table S3). (Inset) Plot of $\tau{}$ vs. 1/T showing
region of crossover between fitting methods. Black line again shows the
relaxation model fit.

## Conclusions

In this work we have utilized the slow magnetic relaxation behavior of 1 to
test the efficacy of a new technique for extracting characteristic relaxation
times that are beyond the scope of standard AC susceptibility experiments.
This proof-of-concept work has shown that the VSM waveform method is viable
and yields results that are a consistent continuation of the fast timescale
dynamics. We hope that given the wide availability of the required
instrumentation and the current interest in understanding and controlling the
rate of magnetic relaxation mechanisms access to another characterization tool
will prove useful. In addition to the expanded timescale, the methods
described herein may prove useful for targeting specific relaxation mechanisms
that occur at non-zero field. For example, with the field tuned to a level-
crossing, quantum-tunneling transitions could be studied in isolation from
competing relaxation mechanisms.

## Acknowledgements

This research was funded through the Office of Naval Research Young
Investigator Award N00014-16-1-2917. The authors thank Drs. Milan Gembicky and
Anthony L. Spek (crystallography) and Dr. Anthony Mrse (NMR) for their expert
assistance.

## References

Supporting Information

###### Contents

  1. 1 Preparative Details
    1. 1.1 General Considerations
    2. 1.2 [K(18-c-6)][Er(hdcCOT)2] (1)
  2. 2 Sample Characterization
    1. 2.1 1H NMR of K[Er(hdcCOT)2]
    2. 2.2 Crystallographic Methods
    3. 2.3 Computational Details
    4. 2.4 Magnetic Data Collection

##  1 Preparative Details

###  1.1 General Considerations

Manipulations involving the synthesis of 1 were carried out in a nitrogen-
atmosphere glovebox. Tetrahydrofuran (THF) and pentane were dried on activated
alumina columns and stored over a 1:1 mixture of 3 and 4 Åamolecular sieves.
1,6-heptadiyne (purchased from Alfa Aesar) was degassed using freeze-pump-thaw
cycles prior to use. Nickel(II) bromide DME complex (Combi-Blocks), 18-crown-6
(18-c-6, Aldrich), anhydrous erbium trichloride (Aldrich), and potassium
graphite (KC8, Strem) were used as received. THF-d8 was purchased from Sigma-
Aldrich and was degassed using freeze-pump-thaw cycles, dried over stage 0 NaK
on silica, and filtered before use. Neutral proligand hdcCOT was synthesized
according to a literature procedure.wender_nickel0-catalyzed_2007-1 1H NMR
spectra were collected at –20 °C on a Jeol ECA 500 spectrometer. CHN elemental
analysis was conducted by Midwest Microlab, Indianapolis, IN.

###  1.2 [K(18-c-6)][Er(hdcCOT)2] (1)

To a $-50$ °C stirring suspension of KC8 (0.368 g, 2.72 mmol, 4 mL THF) was
added a solution of hdcCOT (0.239 g, 1.30 mmol, 6 mL THF). The reaction
mixture was allowed to warm to room temperature and was stirred for a total of
16 hours. Graphite was removed via centrifugation and subsequent filtration
through a glass-fiber filter yielded a dark brown homogeneous solution.
Concentration in vacuo and cooling to $-50$ °C yielded large yellow plates of
K2hdcCOT over the course of 24 hours. The mother liquor was removed from these
plates and they were subsequently dried in vacuo; mass obtained: 0.261 g
(Yield: 76.5%). K2hdcCOT (0.182 g, 0.693 mmol) and ErCl3 (0.095 g, 0.347 mmol)
were combined in 20 mL THF and stirred for 16 hours. The resulting yellow
suspension was centrifuged and the supernatant was separated and dried in
vacuo to yield 0.149 g (74.7%) of K[Er(hdcCOT)2] as a bright yellow powder.
18-crown-6 (0.0635 g, 0.240 mmol) was combined with this powder (0.1382 g,
0.240 mmol) in 20 mL THF and the resulting solution was stirred for 16 hours.
Vapor diffusion of pentane into this solution yielded yellow rods of 1 (200
mg, Yield: 85.0%). CHN analysis (calculated, found) for C40H56ErK: C (57.25,
54.26); H (6.73, 6.51); N (0, 0).

##  2 Sample Characterization

###  2.1 1H NMR of K[Er(hdcCOT)2]

Figure S1: 1H NMR spectrum of K[Er(hdcCOT)2] in THF-d8 at –20 °C. Labelled
signals correspond to THF residual peaks. Figure S2: Downfield 1H NMR spectrum
of K[Er(hdcCOT)2] in THF-d8 at –20 °C. Labelled signals correspond to THF
residual peaks. Unlabelled signals are taken to be decomposition products.
Figure S3: Upfield 1H NMR spectrum of K[Er(hdcCOT)2] in THF-d8 at –20 °C.
Labelled signals correspond hdcCOT protons.

###  2.2 Crystallographic Methods

Single crystal diffraction data for 1 was collected at 100 K on a Bruker
$\kappa$ Diffractometer using a Ga(Kα) METALJET source and a PHOTON II Area
Detector. Data integration was carried out using SAINT and output intensities
were corrected for Lorentz and air absorption effects. Additional absorption
corrections were applied using SADABS. The structure was solved in space group
No. 14 (P21/n) using direct methods with the SHELXTsheldrick_shelxt_2015
program and anisotropic atom positions were refined against F2 data using the
SHELXLsheldrick_crystal_2015 program. Olex2 was used during the refinement
stage as a graphical front-end.dolomanov_olex2_2009 Real and imaginary
anomalous dispersion coefficients for Ga(Kα) radiation were taken from the
Brennan and CowanbrennanSuiteProgramsCalculating1992 model. The position of
all hydrogen atoms were determined using a riding model. Supplementary
crystallographic data can be accessed from the Cambridge Crystallographic Data
Center, CCDC 1939762.

###  2.3 Computational Details

Ab initio electronic structure modelling was carried out at the CASSCF level
using the MOLCAS 8.2 software suite. Input atom coordinates were taken from
crystallographic data and used without further geometry optimization. K,
18-c-6, and solvent THF were excluded from the input geometry. Basis functions
of the ANO-RCC type were generated with the SEWARD module and the quality of a
specific atomic basis function was determined as a function of the atom’s
distance from the Er3+ ion (Er: ANO-RCC-VTZP; atoms bound to Er: ANO-RCC-VDZP;
all other atoms: ANO-RCC-VDZ). To save disk space and reduce calculation cost
two-electron integrals were Cholesky decomposed ($10^{-6}$ cutoff). A
7-orbital, 11-electron activate space (CAS(11,7)) was selected for the CASSCF
calculation which was carried out with the RASSCF module. In this space all 35
configuration-interaction (CI) roots of spin multiplicity 4 and all 112 CI
roots of spin multiplicity 2 were included. Spin-orbit matrix elements between
CAS output wavefunctions were calculated with the RASSI module. SINGLE_ANISO
was used to calculate relevant magnetic properties based on these
multiconfigurational SCF results.

Figure S4: Calculated spectrum of the four lowest energy Kramers states for 1. States (black lines) are labelled by their largest $\pm{}m_{J}$ component (left) and the percentage of that component (right). Transverse magnetic moment matrix elements (dotted lines) are colored according to their respective values (colorbar). Table S1: $J=15/2$ manifold energy spectrum of 1. KD (i) | Mz ($\mu{}_{B}$) |  $E$ (cm–1)
---|---|---
G | 8.976951952 | 0.000
1 | 7.769140677 | 169.229
2 | 0.613462182 | 211.322
3 | 1.834519491 | 298.637
4 | 6.547320548 | 411.560
5 | 3.031744486 | 433.319
6 | 4.722208148 | 542.787
7 | 5.730749227 | 550.231
Table S2: Selected average magnetic moment matrix elements between the $J=15/2$ multiplets of 1. KD (i) | $\langle{}i_{\uparrow}|\hat{M}_{avg}|i_{\downarrow}\rangle{}$ | $\langle{}i_{\uparrow}|\hat{M}_{avg}|i_{\uparrow}+1\rangle{}$ | $\langle{}i_{\uparrow}|\hat{M}_{avg}|i_{\downarrow}+1\rangle{}$
---|---|---|---
G | 0.327451096777E-05 | 0.155011707285E+01 | 0.131694298009E-04
1 | 0.506016013299E-03 | 0.928060135350E-01 | 0.405995456082E-01
2 | 0.318682266522E+01 | 0.316661251988E+01 | 0.720938920092E-01
3 | 0.600396549586E-01 | 0.134946515237E+00 | 0.941722228670E-02
4 | 0.142404299502E-02 | 0.242852634761E+00 | 0.973779023853E-02
5 | 0.485270066726E-01 | 0.285000005858E+01 | 0.300470628704E-01
6 | 0.549287282880E-01 | 0.266558195629E+01 | 0.429984200984E-01
7 | 0.508798096471E-01 |  |

###  2.4 Magnetic Data Collection

Magnetic data were collected under DC scan and VSM scan modes using a Quantum
Design MPMS 3 SQUID Magnetometer with equipped AC susceptibility attachment.
Samples were loaded in custom quartz tubes (D&G Glassblowing Inc.) which were
subsequently flame-sealed under static vacuum. To all samples was added a
portion of melted eicosane wax to abate sample torquing and facilitate thermal
conductivity. Diamagnetic contributions from the sample and eicosane were
subtracted from all static moment data using Pascal’s
constants.bain_diamagnetic_2008 Magnetic relaxation data were collected in DC
scan mode after first equilibrating the sample at a given temperature to a 7 T
field then ramping the field (700 Oe sec–1) to 0 T. Details related to the
Fourier analysis of long-timescale magnetic data are discussed in the main
text. MPMS3 data parsing, fitting, and plotting was performed with MATLAB; the
object-oriented code package and documentation used for all processes is
available at https://www.github.com/RinehartGroup/qdsquid-dataplot.

$\tau^{-1}=\tau_{0}^{-1}exp\left(\frac{-U_{eff}}{k_{B}T}\right)+\tau_{qtm}^{-1}+CT^{2}$
(S1)

Relaxation mechanism equation

Listing 1: Truncated MPMS3 MultiVu Sequence

⬇

1Set Temperature 50K at 50K/min, Fast Settle

2Wait For Temperature, Delay 0 sec, No Action

3

4Magnet Reset

5Wait For Delay 600 secs (10.0 mins), No Action

6

7Set Temperature 2K at 50K/min. Fast Settle

8Wait For Temperature, Delay 120 secs (2.0 mins), No Action

9

10! REMARK - Frequency ~ 0.0264 Hz

11MPMS3 Measure for 0.5 sec at 1 mm every 0 sec Auto-Tracking

12Wait For Delay 2 secs, No Action

13Scan Time 0.0 secs in 20 steps

14 Set Magnetic Field 8.0Oe at 30.00Oe/sec, Linear, Stable

15 Wait For Field, Delay 15 secs, No Action

16 Set Magnetic Field -8.0Oe at 30.00Oe/sec, Linear, Stable

17 Wait For Field, Delay 15 secs, No Action

18End Scan

19Set Magnetic Field 0.0Oe at 30.00Oe/sec, Linear, Stable

20Wait For Field, Delay 0 secs, No Action

21Stop Measurements

22

23Wait For Delay 60 secs (1.0 mins), No Action

24

25! REMARK - Frequency ~ 0.0147 Hz

26MPMS3 Measure for 0.5 sec at 1 mm every 0 sec Auto-Tracking

27Wait For Delay 2 secs, No Action

28Scan Time 0.0 secs in 16 steps

29 Set Magnetic Field 8.0Oe at 30.00Oe/sec, Linear, Stable

30 Wait For Field, Delay 22 secs, No Action

31 Set Magnetic Field -8.0Oe at 30.00Oe/sec, Linear, Stable

32 Wait For Field, Delay 22 secs, No Action

33End Scan

34Set Magnetic Field 0.0Oe at 30.00Oe/sec, Linear, Stable

35Wait For Field, Delay 0 secs, No Action

36Stop Measurements

37

38! REMARK - Frequencies truncated

39

40Wait For Delay 240 secs (4.0 mins), No Action

41

42! REMARK - Frequency ~ 0.00013 Hz

43MPMS3 Measure for 0.5 sec at 1 mm every 0 sec Auto-Tracking

44Wait For Delay 2 secs, No Action

45Set Magnetic Field 8.0Oe at 30.00Oe/sec, Linear, Stable

46Wait For Field, Delay 3840 secs (1.1 hours), No Action

47Set Magnetic Field -8.0Oe at 30.00Oe/sec, Linear, Stable

48Wait For Field, Delay 3840 secs (1.1 hours), No Action

49Set Magnetic Field 0.0Oe at 30.00Oe/sec, Linear, Stable

50Wait For Field, Delay 0 secs, No Action

51Stop Measurements

52

53Set Temperature 300K at 50K/min. Fast Settle

54Wait For Temperature, Delay 1 secs, No Action

Table S3: Fitting values for relaxation mechanism parameters. Note these parameters are meant to highlight differences in the form of the relaxation data determined by different methods and are purely phenomenological. |  $U_{eff}$ (cm–1) |  $\tau_{0}$ (s) |  $\tau_{qtm}$ (s) | C (K–2s–1)
---|---|---|---|---
$\tau_{AC}$ | 147.7(7) | 3.1(1)$\times$10–8 |  |
$\tau_{AC}+\tau_{VSM}$ | 142(3) | 4.7(1.0)$\times$10–8 | 54(7) | 2.2(8)$\times$10-4
$\tau_{AC}+\tau_{str}$ | 146(5) | 3.5(1.1)$\times$10–8 | 141(31) | 6.2(6.2)$\times$10–5
$\tau_{AC}+\tau_{1}$ | 154(13) | 1.8(1.0)$\times$10–8 | 1209(621) | 4.0(4.2)$\times$10–5
$\tau_{AC}+\tau_{2}$ | 145(6) | 3.6(1.2)$\times$10–8 | 80(17) | 6.9(9.7)$\times$10–5
Figure S5: Magnetic susceptibility of 1 between 2 to 300 K. Colored dots are
data measured under a 1000 Oe field and the dotted line is the theoretical
$\chi{}T$ value for a free Er3+ ion. Figure S6: Isothermal magnetization of 1
from –7 to 7 T (Rate = 10 Oe sec–1). Figure S7: In-phase magnetic
susceptibility of 1 bewteen 2 (fuschia) and 24 (yellow) K. Colored points are
susceptibilities measured via standard AC measurements (circles) and extracted
from Fourier analysis of VSM data (triangles). Black lines represent fits to a
generalized Debye model (Equations 1 & 2). Figure S8: Out-of-phase magnetic
susceptibility of 1 bewteen 2 (fuschia) and 24 (yellow) K. Colored points are
susceptibilities measured via standard AC measurements (circles) and extracted
from Fourier analysis of VSM data (triangles). Black lines represent fits a to
generalized Debye model (Equations 1 & 2).

## References


# Lee–Yang zeros in the Rydberg atoms

Chengshu Li<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua
University, Beijing 100084, China  Fan Yang  Institute for Advanced Study,
Tsinghua University, Beijing 100084, China

###### Abstract

Lee–Yang (LY) zeros play a fundamental role in the formulation of statistical
physics in terms of (grand) partition functions, and assume theoretical
significance for the phenomenon of phase transitions. In this paper, motivated
by recent progress in cold Rydberg atom experiments, we explore the LY zeros
in classical Rydberg blockade models. We find that the distribution of zeros
of partition functions for these models in one dimension (1d) can be obtained
analytically. We prove that all the LY zeros are real and negative for such
models with arbitrary blockade radii. Therefore, no phase transitions happen
in 1d classical Rydberg chains. We investigate how the zeros redistribute as
one interpolates between different blockade radii. We also discuss possible
experimental measurements of these zeros.

Keywords: Lee–Yang zeros, Rydberg atom, statistical mechanics

##  I Introduction

Much of modern statistical physics builds on the concept of (grand) partition
functions. For finite systems in general, they are positive, analytic
functions in terms of physical quantities, and, in particular, are smooth
functions with regard to temperature $T$ when the latter is finite.
Historically, this fact caused some confusion as to how discontinuity, as
occurs in phase transitions, can emerge from such a smooth function, before it
became clear that phase transitions can only happen in the thermodynamic
limit. However, it was not until the work of Lee and Yang [1, 2] when a more
detailed description of the relation between the analyticity of partition
functions and phase transitions was put forward.

Exactly 70 years ago in 1952, in two seminal papers [1, 2], Lee and Yang
proposed to study the complex zeros of the (grand) partition function as a
polynomial of some sort of fugacity, and showed how the distribution of such
zeros in the thermodynamic limit relates to the existence and properties of
phase transitions. They proved that if there exists a region containing part
of the positive real axis with no zeros inside, in the corresponding regime of
physical parameters no phase transition can happen. Conversely, if the zeros
accumulate to a point on the positive real axis in the thermodynamic limit,
phase transitions will generally ensue. As a concrete example, they considered
the Ising model with arbitrary ferromagnetic couplings, and proved that all
the zeros are located on the unit circle with the fugacity defined as
$\exp(-\beta h)$, where $h$ is the magnetic field. A direct corollary is that
phase transitions are not possible unless $h=0$ in the Ising model.

The original construction of Lee and Yang focused on spin-1/2 Ising
ferromagnets, or equivalently attractive lattice gases. Generalizations of
this program to higher spins and Heisenberg models soon followed [3, 4, 5, 6,
7, 8, 9, 10]. In the next decades, results on other geometries and
interactions, in particular antiferromagnetic Ising models, became available
[11, 12, 13, 14, 15, 16, 17, 18]. There are also more recent theoretical
progresses as well as generalizations to dynamical and/or quantum systems [19,
20, 21, 22, 23]. The endeavor of characterizing complex zeros of real
polynomials with positive coefficients has attracted interest from the
mathematics community alike [24, 25, 26, 27].

While originally a purely theoretical discussion with no direct experimental
relevance envisioned, Lee–Yang (LY) zeros turn out experimentally measurable
using nuclear magnetic resonance (NMR) techniques [28, 29]. The idea is to
carefully couple a collection of spins to a probe spin, and extract the zeros
from its quantum evolution. Recently, cold Rydberg atoms have attracted
extensive research efforts for their versatile control and measurements
together with strong dipole interactions. Two major frontiers are unfolding —
First, various novel phases and phase transitions emerge from geometries
enabled by single-atom manipulations, including one dimensional (1d) chains
[30, 31] and two dimensional (2d) square [32, 33, 34, 35] and kagome lattices
[36, 37, 38, 39, 40]. Even in the seemingly innocuous 1d case, there are some
debates on the nature of the phase transitions [41, 42, 43, 44, 45, 46], a
caricature of the rich and subtle physics in such systems. Second, the
constrained Hilbert space drastically changes the dynamical behavior, leading
to violations of the eigenstate thermolization hypothesis (ETH) [30], which is
in turn closely related to the quantum many body scar states [47, 48].

Motivated by these exciting progresses, we study how the LY zeros distribute
in classical Rydberg atom systems. Here by “classical” we mean that all the
terms in the Hamiltonian commute with each other, yet the Hilbert space is
still constrained by the Rydberg blockade. One of the important features of
these models, as we discover in our analysis, is that the distribution of
their LY zeros can be obtained rigorously. This is in contrast to general
cases, where the distribution of LY zeros does not have a clear pattern. We
provide a rigorous proof that in 1d for any blockade radius the zeros are real
and negative. Therefore, no phase transitions occur at finite temperature.
This is in agreement with the general belief that classical phase transitions
do not happen in 1d systems at finite temperature. Our analysis adds the
classical 1d Rydberg chain to a small number of models whose distribution of
zeros of the partition function can be obtained analytically, making the study
of the mathematical structure of such models interesting in its own right.
Beyond rigorous analytical discussions, we also numerically investigate how
the zeros redistribute as the Rydberg blockade interpolates between different
radii. Finally, inspired by previous experimental efforts, we propose a way to
measure the LY zeros of the Rydberg chain, a seemingly abstract quantity, in
cold atom experiments, bringing our predictions to physical reality.

##  II Rydberg blockade Hamiltonian, partition function, and LY zeros

We begin with the classical Rydberg blockade Hamiltonian,

$H=-\Delta\sum_{i}n_{i},$ (1)

where $n_{i}=0,1$ is the number of Rydberg excitation on each site $i$ and
$\Delta$ is the detuning. Generally there is a Rabi oscillation term that
couples the ground state and the Rydberg state, which we assume to be small
and ignore here. The effect of the Rydberg blockade is taken into account by
requiring that no two atoms can be excited within a radius of $r$, or that
there can be at most one Rydberg atom in each consecutive $(r+1)$ sites. For a
chain with $n$ sites and open boundary conditions, the number of different
configurations with $m$ Rydberg atoms, $F_{n}^{m}$, is

$F_{n}^{m}=\begin{pmatrix}n-r(m-1)\\\ m\end{pmatrix}$ (2)

from elementary combinatorics. These coefficients satisfy a recursion relation

$F_{n+1}^{m}=F_{n}^{m}+F_{n-r}^{m-1}.$ (3)

To see this, notice that either the $(n+1)$-th site is unoccupied, in which
case one is free to put $m$ Rydberg atoms in the first $n$ sites, giving the
first term, or the $(n+1)$-th site is occupied, blocking sites $n-r+1$ to $n$,
giving the second term. To proceed we introduce the fugacity
$y=\exp(\beta\Delta)$, and the partition function is

$Z_{n}(y)=\sum_{m=0}^{\lfloor(n+r)/(r+1)\rfloor}F_{n}^{m}y^{m},$ (4)

where $\lfloor\cdot\rfloor$ denotes the floor function. The partition function
in turn satisfies

$Z_{n+1}(y)=Z_{n}(y)+yZ_{n-r}(y),$ (5)

together with initial conditions

$Z_{n}(y)=1+ny,\ 1\leq n\leq r+1.$ (6)

From now on we deem $y$ as a complex variable and explore the complex roots
(zeros) of the partition function as a polynomial in $y$. The first main
result we obtain is

######  Proposition 1.

For any $r,n$, all the zeros of the partition function as a polynomial of $y$
are real and negative.

The proof is elementary yet slightly involved. Readers who are not interested
in the proof can skip to the next section.

Figure 1: An illustration of the proof of Lemma 1 for the case $r=2$. Blue
dots denote the zeros and orange triangles (green squares) indicate that the
polynomial is positive (negative) at that point. As a result of the recursion
relation (5), the zeros in each group of three polynomials (e.g.,
$y_{1,1},y_{1,2},y_{1,3}$) determine the sign of the polynomials in the next
group, which in turn guarantees the existence of zeros in between. A few
examples are in order. We label the values of the polynomials with Greek
letters as shown in the figure. First, $\delta<0$ and $\epsilon=0$ give
$\gamma<0$, which together with $\alpha>0$ fixes $y_{2,4}$. Then, $\iota=0$
and $\kappa>0$ give $\theta<0$, which together with $\lambda>0$ fixes
$y_{1,4}$. As a final example, $\eta<0$ and $\iota=0$ give $\zeta<0$, which
together with $\beta>0$ fixes $y_{2,6}$. Note that we only mark the zeros and
the points which are _directly_ used to fix the zeros, e.g., $\alpha$,
$\beta$, $\gamma$, $\epsilon$, $\zeta$, $\theta$, $\iota$, and $\lambda$.

The main idea of the proof is as follows, which is similar to the one given in
[12]: for any fixed $r$, we already know that the partition function is a
polynomial of degree $\lfloor(n+r)/(r+1)\rfloor$. That is, as $n$ increases,
the first $(r+1)$ polynomials have degree one, the next $(r+1)$ ones degree
two, and so on. The claim is proved by finding $\lfloor(n+r)/(r+1)\rfloor$
real roots of the polynomials in each group. This follows from the following
lemma.

######  Lemma 1.

For each $i\geq 0$ and $j=1,\cdots,r+1$, the $i+1$ roots of
$Z_{n=(r+1)i+j}(y)$ are real and different. Furthermore, when they are ordered
by $0>y_{1,n}>y_{2,n}>\cdots>y_{i+1,n}$, they satisfy

$\begin{split}&y_{1,(r+1)i+r+1}>y_{1,(r+1)i+r}>\cdots>y_{1,(r+1)i+1}\\\
>{}&y_{2,(r+1)i+r+1}>y_{2,(r+1)i+r}>\cdots>y_{2,(r+1)i+1}\\\ >{}&\cdots\\\
>{}&y_{i+1,(r+1)i+r+1}>y_{i+1,(r+1)i+r}>\cdots>y_{i+1,(r+1)i+1}.\end{split}$
(7)

###### Proof.

We first prove the case $r=2$ for simplicity, see Fig. 1, but the same proof
works for any $r$, see below. The first 3 polynomials ($i=0$) obviously
satisfy the lemma. If the lemma is satisfied for $Z_{3i+1},Z_{3i+2},Z_{3i+3}$,
from continuity we know the sign of the polynomial in each segment partitioned
by the zeros. For example, $Z_{3i+k}(y)>0$ when $y_{1,3i+k}<y<0$, $k=1,2,3$.
We obtain from the recursion relation

$Z_{3(i+1)+1}(0)>0$ (8)

and

$\begin{split}&Z_{3(i+1)+1}(y_{1,3i+3})\\\
={}&Z_{3i+3}(y_{1,3i+3})+y_{1,3i+3}Z_{3i+1}(y_{1,3i+3})\\\
={}&y_{1,3i+3}Z_{3i+1}(y_{1,3i+3})<0.\end{split}$ (9)

Therefore, there exists a zero of $Z_{3(i+1)+1}$, namely $y_{1,3(i+1)+1}$, in
the open interval $(y_{1,3i+3},0)$ by the intermediate value theorem. Similar
calculations give the other $i+1$ zeros, $y_{m,3(i+1)+1}$, in the intervals
$(y_{m,3i+3},y_{m-1,3i+1})$, where $y_{i+2,3i+3}=-\infty$. Next, the zeros of
$Z_{3(i+1)+2},Z_{3(i+1)+3}$ can be found in the same way, giving inequalities
including $y_{m,3(i+1)+1}<y_{m,3(i+1)+2}<y_{m-1,3i+2}$ and
$y_{m,3(i+1)+2}<y_{m,3(i+1)+3}<y_{m-1,3i+3}$. Combining the results above, we
find

$\begin{split}y_{m,3(i+1)+1}&{}<y_{m,3(i+1)+2}<y_{m,3(i+1)+3}\\\
&{}<y_{m-1,3i+3}<y_{m-1,3(i+1)+1},\end{split}$ (10)

which concludes the proof for $r=2$.

The proof for $r\geq 3$ follows exactly the same spirit and is presented below
for the readers who wish to read a general proof. To facilitate comparison
with the special case of $r=2$ above, we translate the proof to general $r$
verbatim as follows.

The first $(r+1)$ polynomials ($i=0$) obviously satisfy the lemma. If the
lemma is satisfied for $Z_{i(r+1)+j},j=1,\dots,r+1$, from continuity we know
the sign of the polynomial in each segment partitioned by the zeros. For
example, $Z_{i(r+1)+k}(y)>0$ when $y_{1,i(r+1)+k}<y<0$, $k=1,\dots,r+1$. We
obtain from the recursion relation

$Z_{(i+1)(r+1)+1}(0)>0$ (11)

and

$\begin{split}&Z_{(i+1)(r+1)+1}(y_{1,i(r+1)+r+1})\\\
={}&Z_{i(r+1)+r+1}(y_{1,i(r+1)+r+1})\\\
&+y_{1,i(r+1)+r+1}Z_{i(r+1)+1}(y_{1,i(r+1)+r+1})\\\
={}&y_{1,i(r+1)+r+1}Z_{i(r+1)+1}(y_{1,i(r+1)+r+1})<0.\end{split}$ (12)

Therefore, there exists a zero of $Z_{(i+1)(r+1)+1}$, namely
$y_{1,(i+1)(r+1)+1}$, in the open interval $(y_{1,i(r+1)+r+1},0)$ by the
intermediate value theorem. Similar calculations give the other $i+1$ zeros,
$y_{m,(i+1)(r+1)+1}$, in the intervals $(y_{m,i(r+1)+r+1},y_{m-1,i(r+1)+1})$,
where $y_{i+2,i(r+1)+r+1}=-\infty$. Next, the zeros of
$Z_{(i+1)(r+1)+j},j=2,\dots,r+1$ can be found in the same way, giving
inequalities including
$y_{m,(i+1)(r+1)+j}<y_{m,(i+1)(r+1)+j+1}<y_{m-1,i(r+1)+j+1},j=1,\dots,r$.
Combining the results above, we find

$\begin{split}&y_{m,(i+1)(r+1)+1}<y_{m,(i+1)(r+1)+2}<\cdots\\\
<{}&y_{m,(i+1)(r+1)+r+1}<y_{m-1,i(r+1)+r+1}\\\
<{}&y_{m-1,(i+1)(r+1)+1},\end{split}$ (13)

which concludes the proof for general $r$. ∎

This lemma leads directly to Proposition 1, as we have found all the roots of
the partition function. As stated in the lemma, they are all real and
negative.

We note that in two previous papers [49, 50] on integrable systems generalized
from free (para)fermions [51], a result identical to Proposition 1 was
conjectured based on numerical evidence, but a rigorous proof was lacking. The
above proof fills this gap and finds applications in the field of integrable
systems.

##  III Mathematical structure of the distribution of LY zeros

It is possible to formally solve the recursion relation in Eq. (5). We first
extend the initial condition to negative $n$ for convenience. One easily sees
that

$\begin{split}&Z_{n}(y)=1,\ -r\leq n\leq 0,\\\ &Z_{n}(y)=0,\ -2r\leq
n\leq-r-1.\end{split}$ (14)

We then have

$\begin{split}Z_{n}&=(1,0,\cdots,0)\begin{pmatrix}Z_{n}\\\ Z_{n-1}\\\
\vdots\\\ Z_{n-r}\end{pmatrix}\\\ &=(1,0,\cdots,0)\begin{pmatrix}1&&&&y\\\
1&&&&\\\ &1&&&\\\ &&\ddots&&\\\
&&&1&\end{pmatrix}^{n+r}\begin{pmatrix}Z_{-r}\\\ Z_{-r-1}\\\ \vdots\\\
Z_{-2r}\end{pmatrix}\\\ &=(1,0,\cdots,0)\begin{pmatrix}1&&&&y\\\ 1&&&&\\\
&1&&&\\\ &&\ddots&&\\\ &&&1&\end{pmatrix}^{n+r}\begin{pmatrix}1\\\ 0\\\
\vdots\\\ 0\end{pmatrix}.\end{split}$ (15)

To proceed we need to diagonalize

$A=\begin{pmatrix}1&&&&y\\\ 1&&&&\\\ &1&&&\\\ &&\ddots&&\\\
&&&1&\end{pmatrix}.$ (16)

The characteristic polynomial reads

$\lambda^{r+1}-\lambda^{r}-y=0.$ (17)

Assuming all the eigenvalues $\lambda_{i}$ to be different (we will verify
this soon), we see that $A$ is diagonalized by

$U=\begin{pmatrix}1&1&\cdots&1\\\
\lambda_{1}^{-1}&\lambda_{2}^{-1}&\cdots&\lambda_{r+1}^{-1}\\\
\lambda_{1}^{-2}&\lambda_{2}^{-2}&\cdots&\lambda_{r+1}^{-2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\lambda_{1}^{-r}&\lambda_{2}^{-r}&\cdots&\lambda_{r+1}^{-r}\end{pmatrix},$
(18)

which is a Vandermonde matrix of $\lambda_{i}^{-1}$. We arrive at

$\begin{split}Z_{n}&=(U\,\mathrm{diag}\\{\lambda_{1}^{n+r},\cdots,\lambda_{r+1}^{n+r}\\}\,U^{-1})_{1,1}\\\
&=\frac{1}{\prod_{i}\lambda_{i}}\sum_{i}\frac{\lambda_{i}^{n+r+1}}{\prod_{j}^{\prime}(\lambda_{j}^{-1}-\lambda_{i}^{-1})}.\end{split}$
(19)

For $r=1$, or nearest-neighbor blockade, $\lambda_{1,2}=(1\pm\sqrt{1+4y})/2$,

$\begin{split}Z_{n}&=\frac{\lambda_{1}^{n+2}-\lambda_{2}^{n+2}}{\lambda_{1}\lambda_{2}(\lambda_{2}^{-1}-\lambda_{1}^{-1})}=\frac{\lambda_{1}^{n+2}-\lambda_{2}^{n+2}}{\lambda_{1}-\lambda_{2}}\\\
&=\frac{1}{\sqrt{4y+1}}\\\
&\times\bigg{(}\Big{(}\frac{1+\sqrt{4y+1}}{2}\Big{)}^{n+2}-\Big{(}\frac{1-\sqrt{4y+1}}{2}\Big{)}^{n+2}\bigg{)}.\end{split}$
(20)

To find the zeros we set $Z_{n}(y)=0$ and get

$\Big{(}\frac{1+\sqrt{4y+1}}{1-\sqrt{4y+1}}\Big{)}^{n+2}=1.$ (21)

It’s straightforward to find

$y=-\frac{1+\tan^{2}(\phi/2)}{4},$ (22)

where

$\phi=2\pi m/(n+2),m=1,\cdots,\lfloor(n+1)/2\rfloor,$ (23)

which is real, negative and always smaller than $-1/4$.

In order for Eq. (17) to have degenerate roots, the latter must also be roots
of the derivative of the original polynomial. This can only happen when
$\lambda=r/(r+1)$, and correspondingly $y_{0}=-r^{r}/(r+1)^{r+1}$. This $y$ is
also where the number of real roots of Eq. (17) changes, hinting at its
special role. Indeed, we have the following result.

######  Proposition 2.

All the zeros of the partition function as in Proposition 1 are smaller than
$y_{0}$.

###### Proof.

We only need to prove that when $y_{0}\leq y<0$, $Z_{n}(y)$ is always greater
than zero. We prove it by showing a stronger condition,

$\frac{Z_{n+1}}{Z_{n}}>\frac{r}{r+1}.$ (24)

This obviously holds for $Z_{n},n=-r,\cdots,-1$. Assume this to be true for
$Z_{n},n=i,\cdots,i+r-1$, then

$\frac{Z_{i+r}}{Z_{i}}>\Big{(}\frac{r}{r+1}\Big{)}^{r}.$ (25)

The recursion relation gives

$\frac{Z_{i+r+1}}{Z_{i+r}}=1+y\frac{Z_{i}}{Z_{i+r}}>1+y\Big{(}\frac{r+1}{r}\Big{)}^{r}.$
(26)

Then we have

$\frac{Z_{i+r+1}}{Z_{i+r}}>\frac{r}{r+1}\Leftarrow
1+y\Big{(}\frac{r+1}{r}\Big{)}^{r}\geq\frac{r}{r+1}\Leftrightarrow y\geq
y_{0}.$ (27)

∎

We also have the following corollary:

######  Corollary 1.

For $y_{0}\leq y<0$, we have $Z_{n}(y)\to 0$ as $n\to\infty$.

###### Proof.

Using the recursion relation and Proposition 2, we have

$\frac{Z_{n+1}(y)}{Z_{n}(y)}=1+y\frac{Z_{n-r}(y)}{Z_{n}(y)}<1.$ (28)

Thus, $Z_{n}$ decreases monotonically with index $n$. By realizing

$\frac{Z_{n-r}(y)}{Z_{n}(y)}>1,$ (29)

Eq. (28) can be made stronger

$\frac{Z_{n+1}(y)}{Z_{n}(y)}=1+y\frac{Z_{n-r}(y)}{Z_{n}(y)}<1+y.$ (30)

At last, the squeeze theorem gives

$\lim_{n\to\infty}Z_{n}(y)=0,\forall y\in[y_{0},0).$ (31)

∎

##  IV Interpolations between different blockade radii

To better understand Rydberg chain systems, we investigate how the zeros
redistribute as the blockade radius increases. To begin with we focus on the
case $r=1$ to $r=2$, or nearest- to next-nearest-neighbor blockade. For this
purpose we add a next-nearest-neighbor interacting term to the $r=1$
Hamiltonian and the new Hamiltonian reads

$H^{\prime}=-\Delta\sum_{i}n_{i}+V\sum_{i}n_{i}n_{i+2},$ (32)

with the same condition that no two nearest neighbors can be excited
simultaneously. This is the classical ($w=0$) version of the
Fendley–Sengupta–Sachdev model [41]. We note that this model is different from
the classical Rydberg blockade model discussed in previous sections. However,
the same physics should be recovered in the following limits: $V=0$
corresponds to blockade radius $r=1$ and $V=\infty$ corresponds to $r=2$. The
partition function of the new Hamiltonian is then

$Z^{\prime}=\sum_{n}y^{n}\sum_{\\{n_{i}\\}}\exp(-\beta
V\sum_{i}n_{i}n_{i+2}).$ (33)

Figure 2: The zeros of the partition function (33) for $n=36$. The largest
$n/3=12$ of them converge to that of $r=2$ when $V$ goes to $\infty$, while
the others diverge to $-\infty$ exponentially, which are omitted in the lower
panel. Figure 3: The zeros of the partition function (35) for $n=36$. The
largest $n/4=9$ of them converge to that of $r=3$ when $V$ goes to $\infty$,
while the others diverge to $-\infty$ exponentially, which are omitted in the
lower panel.

The Hamiltonian (32) facilitates the understanding of Rydberg chains in the
following two aspects. First, it helps us understand Rydberg chains beyond the
classical Rydberg blockade model. In this aspect, the distribution of LY zeros
of partition function (33) is obtained numerically. We calculate the zeros for
$n=24,\cdots,36$. While for negative $V$ the zeros are in general complex, we
find that when $V$ is positive the zeros are real. Therefore, we conjecture
that the LY zeros remain living on the negative real axis for $V>0$.
Consequently, no phase transitions happen at finite temperature for such
systems with repulsive interactions.

Second, it sheds light on how the LY zeros of classical Rydberg blockade
models with different blockade radii connect with each other. As stated in the
previous sections, for models with blockade radius $r=1$ and $r=2$, the number
of LY zeros are $\lfloor(n+1)/2\rfloor$ and $\lfloor(n+2)/3\rfloor$,
respectively. By numerically studying how the LY zeros evolve during the
interpolation given by Eq. (32), we can gain some insight on how the LY zeros
of Rydberg blockade models with different blockade radii are connected with
each other, and how the number of LY zeros changes as $r$ changes to $r+1$.

Among the $\lfloor(n+1)/2\rfloor$ LY zeros for finite $\beta V$, there is a
sharp difference between the distributions of the largest
$\lfloor(n+2)/3\rfloor$ zeros and the others as $V$ goes to $\infty$: the
former converge to that of the $r=2$ result, while the latter diverge to
$-\infty$ exponentially. The convergence of the larger roots is a ramification
of the fact that as $\beta V$ goes to $\infty$ we should recover the physics
of the model with blockade radius $r=2$. Thus, we conclude that as $r$
increases from 1 to 2, the largest $\lfloor(n+2)/3\rfloor$ zeros of the
blockade model with $r=1$ evolve into each of the LY zeros of the blockade
model with $r=2$. The change in the number of LY zeros of classical Rydberg
blockade models with different blockade radii is attributed to the fact that
the rest of the zeros in the interpolation process eventually diverge to
$-\infty$. Numerically, we find that the dichotomy becomes apparent with
moderate $\beta V\simeq 5$. The results for $n=36$ are shown in Fig. 2.

Due to the unified structures among different blockade radii discussed above,
we expect the interpolations for general $r$ to behave similarly to that
between $r=1$ and $r=2$. Indeed, we numerically check the one between $r=2$
and $r=3$ with an interpolation by a repulsive next-next-nearest-neighbor
interaction

$H^{\prime\prime}=-\Delta\sum_{i}n_{i}+V\sum_{i}n_{i}n_{i+3},$ (34)

with the condition that no two atoms within a radius of two can be excited
simultaneously. The corresponding partition function reads

$Z^{\prime\prime}=\sum_{n}y^{n}\sum_{\\{n_{i}\\}}\exp(-\beta
V\sum_{i}n_{i}n_{i+3}),$ (35)

where $n_{i}$ now sum over configurations with the $r=2$ constraint. The
largest $\lfloor(n+3)/4\rfloor$ zeros approach the $r=3$ results, and the
others diverge to $-\infty$ as $V$ goes to $\infty$, see Fig. 3.

##  V Experimental measurement of the LY zeros

Amenability to individual manipulations and measurements empowers cold Rydberg
atoms as an ideal platform to measure LY zeros experimentally. Here we propose
a setup similar to those in [28, 29] where the zeros can be measured
dynamically.

Figure 4: (a, b) The experimental setup for (a) open boundary condition and
(b) periodic boundary condition. Blue (red) balls denote atoms in the ground
(Rydberg) state, and the green ball denotes the probe atom. Typical time
dependence of $|L(t)|$ as defined in Eq. (38) is shown in (c). Here we choose
$n=24$, $r=1$. The difference between the zeros ($\beta\Delta=-1.37,2.85$) and
non-zeros ($\beta\Delta=-3,2$) is clearly visible.

We first assemble atoms (system, denoted by s) uniformly on a circle, either
to a C-shaped arc for open boundary condition or an O-shaped circle for
periodic boundary condition, with detuning $\Delta$ and inverse temperature
$\beta$, see Fig. 4(a, b). Then we put a probe (denoted by p) atom on the axis
normal to the circle plane, not necessarily in the same plane as the circle to
allow for more general couplings while maintaining a homogeneous one-to-all
interaction,

$H_{\mathrm{sp}}=\lambda n_{\mathrm{p}}\sum_{i}n_{i}=\lambda
n_{\mathrm{p}}n_{\mathrm{s}},$ (36)

where $\lambda$ is a coupling constant and $n_{\mathrm{p}}$ denotes Rydberg
number operator of the probe atom. By preparing the probe atom as a
superposition of the ground state $|0\rangle$ and the Rydberg state
$|1\rangle$, $(|0\rangle+|1\rangle)/\sqrt{2}$, it evolves with a phase factor
proportional to the total Rydberg atom number, weighted by the density matrix
of the system, giving rise to measurable quantities

$\begin{split}\langle\sigma^{x}_{\mathrm{p}}(t)\rangle&=Z(y)^{-1}\sum_{n_{\mathrm{s}}}y^{n_{\mathrm{s}}}\sum_{\\{n_{i}\\}}e^{-\beta
H_{\mathrm{s}}}\cos(\lambda n_{\mathrm{s}}t),\\\
\langle\sigma^{y}_{\mathrm{p}}(t)\rangle&=Z(y)^{-1}\sum_{n_{\mathrm{s}}}y^{n_{\mathrm{s}}}\sum_{\\{n_{i}\\}}e^{-\beta
H_{\mathrm{s}}}\sin(\lambda n_{\mathrm{s}}t).\end{split}$ (37)

Here, $\sigma^{x}_{\mathrm{p}}=|1\rangle\langle 0|+|0\rangle\langle 1|$ and
$\sigma^{y}_{\mathrm{p}}=i(|0\rangle\langle 1|-|1\rangle\langle 0|)$ can be
measured by first performing proper rotations to change to the $x,y$ basis.
Combining the two results gives

$L(t)=\langle\sigma^{x}_{\mathrm{p}}(t)+i\sigma^{y}_{\mathrm{p}}(t)\rangle=Z(e^{\beta\Delta+i\lambda
t})/Z(e^{\beta\Delta}).$ (38)

By tuning $\Delta$ and $t$, one can identify when $L(t)$ goes to zero, which
in turn gives the zeros of the partition function. In particular, the results
we show in the previous sections indicate that $\lambda
t=(2n+1)\pi,n\in\mathbb{N}$, see Fig. 4(c).

##  VI Conclusions

Motivated by the recent advances in cold Rydberg atom experiments, we study LY
zeros in one-dimensional classical Rydberg blockade systems. We find that the
distribution of LY zeros of such models can be obtained analytically. We
rigorously prove that for general blockade radii the LY zeros are real and
negative. Thus, there are no phase transitions at finite temperature. As we
have shown in this paper, the 1d Rydberg chain is one of a small number of
known models, whose distribution of LY zeros can be obtained analytically.
Thus, its mathematical structures may merit further examination. In addition,
we numerically find that these zeros keep living in the negative real axis,
when one turns on a next-nearest-neighbor repulsive interaction to interpolate
between the nearest- and next-nearest-neighbor blockades. We propose a way to
experimentally observe these zeros by coupling the system to a probe atom and
performing dynamical measurements thereon.

## Acknowledgments

We thank Hui Zhai for helpful discussions. Both authors are supported by the
International Postdoctoral Exchange Fellowship Program and the Shuimu Tsinghua
Scholar Program.

## References

  * Yang and Lee [1952] C. N. Yang and T. D. Lee, Statistical theory of equations of state and phase transitions. I. Theory of condensation, Phys. Rev. 87, 404 (1952).
  * Lee and Yang [1952] T. D. Lee and C. N. Yang, Statistical theory of equations of state and phase transitions. II. Lattice gas and Ising model, Phys. Rev. 87, 410 (1952).
  * Asano [1968] T. Asano, Generalization of the Lee-Yang theorem, Prog. Theor. Phys. 40, 1328 (1968).
  * Suzuki [1968a] M. Suzuki, Theorems on the Ising model with general spin and phase transition, J. Math. Phys. 9, 2064 (1968a).
  * Suzuki [1968b] M. Suzuki, Theorems on extended Ising model with applications to dilute ferromagnetism, Prog. Theor. Phys. 40, 1246 (1968b).
  * Griffiths [1969] R. B. Griffiths, Rigorous results for Ising ferromagnets of arbitrary spin, J. Math. Phys. 10, 1559 (1969).
  * Asano [1970] T. Asano, Theorems on the partition functions of the Heisenberg ferromagnets, J. Phys. Soc. Jpn. 29, 350 (1970).
  * Ruelle [1971] D. Ruelle, Extension of the Lee-Yang circle theorem, Phys. Rev. Lett. 26, 303 (1971).
  * Suzuki and Fisher [1971] M. Suzuki and M. E. Fisher, Zeros of the partition function for the Heisenberg, ferroelectric, and general Ising models, J. Math. Phys. 12, 235 (1971).
  * Kurtze and Fisher [1978] D. A. Kurtze and M. E. Fisher, The Yang-Lee edge singularity in spherical models, J. Stat. Phys. 19, 205 (1978).
  * Lieb and Ruelle [1972] E. H. Lieb and D. Ruelle, A property of zeros of the partition function for Ising spin systems, J. Math. Phys. 13, 781 (1972).
  * Heilmann and Lieb [1972] O. J. Heilmann and E. H. Lieb, Theory of monomer-dimer systems, Commun. Math. Phys. 25, 190 (1972).
  * Dobrushin  _et al._ [1985] R. L. Dobrushin, J. Kolafa, and S. B. Shlosman, Phase diagram of the two-dimensional Ising antiferromagnet (computer-assisted proof), Commun. Math. Phys. 102, 89 (1985).
  * Beauzamy [1996] B. Beauzamy, On complex Lee and Yang polynomials, Commun. Math. Phys. 182, 177 (1996).
  * Kim [2004] S.-Y. Kim, Yang-Lee zeros of the antiferromagnetic Ising model, Phys. Rev. Lett. 93, 130604 (2004).
  * Hwang and Kim [2010] C.-O. Hwang and S.-Y. Kim, Yang–Lee zeros of triangular Ising antiferromagnets, Phys. A: Stat. Mech. Appl. 389, 5650 (2010).
  * Lebowitz  _et al._ [2012] J. L. Lebowitz, D. Ruelle, and E. R. Speer, Location of the Lee-Yang zeros and absence of phase transitions in some Ising spin systems, J. Math. Phys. 53, 095211 (2012).
  * Lebowitz and Scaramazza [2016] J. L. Lebowitz and J. A. Scaramazza, A note on Lee–Yang zeros in the negative half-plane, J. Phys.: Condens. Matter 28, 414004 (2016).
  * Heyl  _et al._ [2013] M. Heyl, A. Polkovnikov, and S. Kehrein, Dynamical quantum phase transitions in the transverse-field Ising model, Phys. Rev. Lett. 110, 135704 (2013).
  * Brandner  _et al._ [2017] K. Brandner, V. F. Maisi, J. P. Pekola, J. P. Garrahan, and C. Flindt, Experimental determination of dynamical Lee-Yang zeros, Phys. Rev. Lett. 118, 180601 (2017).
  * Deger and Flindt [2019] A. Deger and C. Flindt, Determination of universal critical exponents using Lee-Yang theory, Phys. Rev. Research 1, 023004 (2019).
  * Deger  _et al._ [2020] A. Deger, F. Brange, and C. Flindt, Lee-Yang theory, high cumulants, and large-deviation statistics of the magnetization in the Ising model, Phys. Rev. B 102, 174418 (2020).
  * Kist  _et al._ [2021] T. Kist, J. L. Lado, and C. Flindt, Lee-Yang theory of criticality in interacting quantum many-body systems, Phys. Rev. Research 3, 033206 (2021).
  * Kurtz [1992] D. C. Kurtz, A sufficient condition for all the roots of a polynomial to be real, Am. Math. Mon. 99, 259 (1992).
  * Borcea and Brändén [2009a] J. Borcea and P. Brändén, The Lee-Yang and Polya-Schur programs. I. Linear operators preserving stability, Invent. Math. 177, 541 (2009a).
  * Borcea and Brändén [2009b] J. Borcea and P. Brändén, The Lee-Yang and Polya-Schur programs. II. Theory of stable polynomials and applications, Commun. Pure Appl. Math. 62, 1595 (2009b).
  * Ruelle [2010] D. Ruelle, Characterization of Lee-Yang polynomials, Ann. Math. 171, 589 (2010).
  * Wei and Liu [2012] B.-B. Wei and R.-B. Liu, Lee-Yang zeros and critical times in decoherence of a probe spin coupled to a bath, Phys. Rev. Lett. 109, 185701 (2012).
  * Peng  _et al._ [2015] X. Peng, H. Zhou, B.-B. Wei, J. Cui, J. Du, and R.-B. Liu, Experimental observation of Lee-Yang zeros, Phys. Rev. Lett. 114, 010601 (2015).
  * 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, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum Kibble–Zurek mechanism and critical dynamics on a programmable Rydberg simulator, Nature 568, 207 (2019).
  * Satzinger  _et al._ [2021] K. J. Satzinger _et al._ , Realizing topologically ordered states on a quantum processor, Science 374, 1237 (2021).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Samajdar  _et al._ [2020] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Complex density wave orders and quantum phase transitions in a model of square-lattice Rydberg atom arrays, Phys. Rev. Lett. 124, 103601 (2020).
  * Kalinowski  _et al._ [2022] M. Kalinowski, R. Samajdar, R. G. Melko, M. D. Lukin, S. Sachdev, and S. Choi, Bulk and boundary quantum phase transitions in a square Rydberg atom array, Phys. Rev. B 105, 174417 (2022).
  * Verresen  _et al._ [2021] R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of toric code topological order from Rydberg blockade, Phys. Rev. X 11, 031005 (2021).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of Rydberg atoms on a kagome lattice, Proc. Natl. Acad. Sci. U.S.A. 118, e2015785118 (2021).
  * Cheng  _et al._ [2021] Y. Cheng, C. Li, and H. Zhai, Variational approach to quantum spin liquid in a Rydberg atom simulator, arXiv:2112.13688 (2021).
  * Giudici  _et al._ [2022] G. Giudici, M. D. Lukin, and H. Pichler, Dynamical preparation of quantum spin liquids in Rydberg atom arrays, Phys. Rev. Lett. 129, 090401 (2022).
  * Fendley  _et al._ [2004] P. Fendley, K. Sengupta, and S. Sachdev, Competing density-wave orders in a one-dimensional hard-boson model, Phys. Rev. B 69, 075106 (2004).
  * Samajdar  _et al._ [2018] R. Samajdar, S. Choi, H. Pichler, M. D. Lukin, and S. Sachdev, Numerical study of the chiral ${\mathbb{Z}}_{3}$ quantum phase transition in one spatial dimension, Phys. Rev. A 98, 023614 (2018).
  * Giudici  _et al._ [2019] G. Giudici, A. Angelone, G. Magnifico, Z. Zeng, G. Giudice, T. Mendes-Santos, and M. Dalmonte, Diagnosing Potts criticality and two-stage melting in one-dimensional hard-core boson models, Phys. Rev. B 99, 094434 (2019).
  * Chepiga and Mila [2019] N. Chepiga and F. Mila, Floating phase versus chiral transition in a 1d hard-boson model, Phys. Rev. Lett. 122, 017205 (2019).
  * Rader and Läuchli [2019] M. Rader and A. M. Läuchli, Floating phases in one-dimensional Rydberg Ising chains, arXiv:1908.02068 (2019).
  * Maceira  _et al._ [2022] I. A. Maceira, N. Chepiga, and F. Mila, Conformal and chiral phase transitions in Rydberg chains, arXiv:2203.01163 (2022).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nat. Phys. 14, 745 (2018).
  * Serbyn  _et al._ [2021] M. Serbyn, D. A. Abanin, and Z. Papić, Quantum many-body scars and weak breaking of ergodicity, Nat. Phys. 17, 675 (2021).
  * Alcaraz and Pimenta [2020a] F. C. Alcaraz and R. A. Pimenta, Free fermionic and parafermionic quantum spin chains with multispin interactions, Phys. Rev. B 102, 121101(R) (2020a).
  * Alcaraz and Pimenta [2020b] F. C. Alcaraz and R. A. Pimenta, Integrable quantum spin chains with free fermionic and parafermionic spectrum, Phys. Rev. B 102, 235170 (2020b).
  * Fendley [2019] P. Fendley, Free fermions in disguise, J. Phys. A: Math. Theor. 52, 335002 (2019).