††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Classical fully packed loop model with attractive interactions on the square
lattice

Bhupen Dabholkar  Laboratoire de Physique Théorique, Université de Toulouse,
CNRS, UPS, France  Xiaoxue Ran  Department of Physics and HKU-UCAS Joint
Institute of Theoretical and Computational Physics,The University of Hong
Kong, Pokfulam Road, Hong Kong SAR, China  Junchen Rong  Institut des Hautes
Études Scientifiques, 91440 Bures-sur-Yvette, France  Zheng Yan  Department of
Physics and HKU-UCAS Joint Institute of Theoretical and Computational
Physics,The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China  G.
J. Sreejith  IISER Pune, Dr Homi Bhabha Road, Pune 411008, India  Zi Yang Meng
<EMAIL_ADDRESS>Department of Physics and HKU-UCAS Joint Institute of
Theoretical and Computational Physics,The University of Hong Kong, Pokfulam
Road, Hong Kong SAR, China  Fabien Alet<EMAIL_ADDRESS>Laboratoire de
Physique Théorique, Université de Toulouse, CNRS, UPS, France

###### Abstract

We study a classical model of fully packed loops on the square lattice, which
interact through attractive loop segment interactions between opposite sides
of plaquettes. This study is motivated by effective models of interacting
quantum matter arising in frustrated magnets or Rydberg atom arrays, for which
loop degrees of freedom appear at low energy. Through the combination of Monte
Carlo simulations and of an effective height field theory, we find that the
critical point known to occur at infinite temperature gives rise to a high-
temperature critical phase with floating exponents. At lower temperature, the
system transitions via a Kosterlitz-Thouless phase transition to a nematic
phase where lattice rotation symmetry is broken. We discuss consequences for
the phase diagram of the quantum loop model on the same lattice.

##  I Introduction

An important notion in the renormalization group theory is the emergence of
effective degrees of freedom at low energies. These new degrees of freedom can
have local structures which take the form of a constraint. For instance, for
degrees of freedom that live on the bonds of a lattice, a gauge-like condition
can emerge which requires that every site of the lattice is touched by a fixed
number of occupied bonds. Related statistical mechanical models such as dimer
or loop models arise as effective theories in many physical situations, such
as in frustrated magnetic systems [1, 2] as, e.g. in spin ice [3], Rydberg
atom arrays [4, 5, 6, 7], models of high-${\rm T}_{c}$ superconductors [8],
adsorption physics [9], quantum Hall effects [10, 11], topological order [12],
deconfined quantum critical points [13, 14, 15, 16, 17, 18, 19, 20, 21, 22],
etc. Loop models also have a long history in statistical physics [23, 24, 25,
13, 26], in relation to Potts models [27], Temperley-Lieb algebras [28],
polymers and O($N$) models [29, 30], Schramm-Loewner evolution [31], or
percolation. These models often assign fixed fugacity for loops [24, 23], but
there are few results when the loop segments interact [32], even though loop
interactions naturally arise in effective models of quantum condensed matter
[8, 33].

In this work, we study a two-dimensional (2D) classical statistical mechanical
model of fully packed loops which attract locally. With the help of a
directed-loop Monte Carlo algorithm [34, 35, 36, 37, 38, 39] and a Coulomb gas
[24] approach formulated in terms of a height-field description of the loop
constraint [40, 41], we obtain evidence for the existence of a finite-
temperature Kosterlitz-Thouless (KT) transition separating a high-temperature
critical phase from a low-temperature nematic phase. Our results have
similarities with those obtained for the classical dimer model with attractive
interactions [42, 35, 43], albeit with specific differences that we highlight.

Aside from their interest in two-dimensional statistical mechanics in
extending previous works on loop models [23, 24, 30, 25, 26, 44], our results
are also relevant for quantum-constrained models. First, the ground-state wave
function at a Rokhsar-Kivelson point [8] (or its generalizations [45, 46]) in
the phase diagram of quantum loop models (QLM) maps to the partition function
of a classical loop model. It is possible to construct extended quantum loop
models (following the prescription in Ref. [45], see an example in Ref. [47]
for a dimer model) whose ground-phase diagram is entirely given by the finite-
temperature phase diagram of a classical (interacting) loop model, as the one
we describe in this work. The second connection is made by realizing that the
phase diagram of the classical model and the methods we use in its inference
can serve to guide us in mapping out the finite-temperature phase diagram [48]
and transitions of the quantum loop model [49, 48, 50, 51, 52, 53, 54] (see
e.g. the finite-temperature phase diagram of the quantum dimer model [55]).
Such quantum-constrained models host a rich set of phases [56, 57, 58, 7, 59,
60, 61, 53, 54, 51, 52, 49, 48, 50] and have recently been shown to be
relevant in the context of Rydberg atom arrays [62, 63, 64, 65, 4, 5, 7, 66],
where the Rydberg blockade effectively implements the loop or dimer
constraint.

The rest of the paper is organized as follows: In Sec. II, we introduce the
classical loop model. Section III introduces the different physical
observables computed in this work. Section IV provides a field theoretical
perspective to the model and phase diagram in the form of a Coulomb gas
analysis. The Monte Carlo simulation results and their analysis are given in
Sec. V, in which, Sec. V.1 presents winding number fluctuations and Sec. V.2
an analysis of the low-temperature order parameter and its susceptibility. In
Sec. V.3, we discuss the behavior of various correlation functions in the
high-temperature phase. These results are analyzed in light of the Coulomb gas
predictions of Sec. IV. We present our conclusions and some perspectives in
Sec. VI. Appendix A describes the Monte Carlo directed loop algorithm used in
our numerical study and Appendix B contains further results on correlation
functions.

##  II Model and methods

Configurations. —  Configurations of the fully packed loop model on a square
lattice require two loop segments (or “dimers”) to touch each site of a square
lattice, and are in one-to-one correspondence with configurations of the six-
vertex model [67, 68, 69, 70, 71, 72]. The ice-rule constraint of the six-
vertex model associates an arrow on each bond and only allows vertices which
have two arrows pointing inwards and two outwards from the lattice site. Under
this constraint, there are six possible vertex configurations on the square
lattice as shown in Fig. 1 (a). The mapping from the six-vertex model to the
loop model on the square lattice is illustrated in Fig. 1 (b). If we place
dimers on two incoming arrows on all sites of a sublattice of the square
lattice, dimers will collectively form fully packed loops as every site is
touched by exactly two dimers (‘loop segments’).

Figure 1: (a) The allowed vertex types for the 6-vertex model. (b)
Correspondence between a 6-vertex configuration and the fully packed loop
configuration on the square lattice. The solid and open circles represent
sites of the A and B sublattices, respectively. Placing dimers on all incoming
arrows of the vertices on the A-sublattice produces a fully-packed loop
configuration.

Energetics —  Loop or vertex models often associate a fugacity with each
closed loop or to each type of vertex respectively, to define the
corresponding partition function [23]. The model that we study here associates
an interaction energy term between proximate parallel loop segments, similar
to the classical interacting dimer models [42, 35]. We consider the following
partition function and energy for an interacting fully-packed loop model on
the square lattice

$\displaystyle Z$ $\displaystyle=$ $\displaystyle\sum_{c}e^{-\beta E_{c}}$
$\displaystyle E_{c}$ $\displaystyle=$ $\displaystyle
V(N(\begin{picture}(22.0,10.0)(-8.0,-7.8)\put(-3.0,-10.0){\line(1,0){12.0}}\put(8.2,-10.0){\rotatebox{90.0}{\rule[-1.0pt]{11.38092pt}{1.99997pt}}}\put(-4.0,-10.0){\rotatebox{90.0}{\rule[-1.0pt]{11.38092pt}{1.99997pt}}}\put(9.0,0.0){\line(-1,0){12.0}}
\put(-3.0,-10.0){\circle*{4.0}}\put(9.0,-10.0){\circle*{4.0}}\put(-3.0,0.0){\circle*{4.0}}\put(9.0,0.0){\circle*{4.0}}
\end{picture})+N(\begin{picture}(22.0,10.0)(-8.0,-7.8)\put(-3.0,-10.0){\rule[-1.0pt]{11.38092pt}{1.99997pt}}\put(9.0,-10.0){\line(0,10){11.0}}\put(-3.0,0.0){\line(0,-10){12.0}}\put(-3.0,0.0){\rule[-1.0pt]{11.38092pt}{1.99997pt}}
\put(-3.0,-10.0){\circle*{4.0}}\put(9.0,-10.0){\circle*{4.0}}\put(-3.0,0.0){\circle*{4.0}}\put(9.0,0.0){\circle*{4.0}}
\end{picture})),$ (1)

where the summation in the partition function $Z$ is over all fully-packed
loop configurations on the square lattice and $\beta=1/T$ is the inverse
temperature. We assign an energy $E_{c}$ to each covering in which
$(N(\begin{picture}(22.0,10.0)(-8.0,-7.8)\put(-3.0,-10.0){\line(1,0){12.0}}\put(8.2,-10.0){\rotatebox{90.0}{\rule[-1.0pt]{11.38092pt}{1.99997pt}}}\put(-4.0,-10.0){\rotatebox{90.0}{\rule[-1.0pt]{11.38092pt}{1.99997pt}}}\put(9.0,0.0){\line(-1,0){12.0}}
\put(-3.0,-10.0){\circle*{4.0}}\put(9.0,-10.0){\circle*{4.0}}\put(-3.0,0.0){\circle*{4.0}}\put(9.0,0.0){\circle*{4.0}}
\end{picture})+N(\begin{picture}(22.0,10.0)(-8.0,-7.8)\put(-3.0,-10.0){\rule[-1.0pt]{11.38092pt}{1.99997pt}}\put(9.0,-10.0){\line(0,10){11.0}}\put(-3.0,0.0){\line(0,-10){12.0}}\put(-3.0,0.0){\rule[-1.0pt]{11.38092pt}{1.99997pt}}
\put(-3.0,-10.0){\circle*{4.0}}\put(9.0,-10.0){\circle*{4.0}}\put(-3.0,0.0){\circle*{4.0}}\put(9.0,0.0){\circle*{4.0}}
\end{picture}))$ counts the number of plaquettes with parallel loop segments.
Note that there is no energy assigned to a plaquette that has more than two
loop segments. Here we set $V=-1$, which corresponds to attractive
interactions between loop segments. We assume periodic boundary conditions for
square lattices of linear size $L$ with $N=L^{2}$ total number of sites. The
model Eq. (II) is the limiting case of the quantum loop model on the square
lattice [49] obtained when kinetic terms vanish in the QLM. To the best of our
knowledge, this fully-packed loop model with aligning interactions has never
been studied earlier.

Limiting cases —  The model admits two simple limits. At infinite temperature,
it is equivalent to the 6-vertex model at the ice-point with equal fugacities
for all vertices in Fig. 1 (a) which is critical with power-law correlators
(see the precise description below). At $T=0$, there are two configurations
which minimize the energy ($E_{0}=-L^{2}$). These are nematic configurations
with $L$ horizontal or vertical loops that wrap around the boundary. The
$\pi/2$ lattice rotation symmetry is broken at $T=0$, and since this model
admits only discrete energies - the first excited states have energies
$E_{1}=-L^{2}+4$ \- we expect a finite-temperature transition into a low-
temperature nematic phase. As will be shown below, this transition is of
Kosterlitz-Thouless (KT) type.

While the two limiting phases (critical and nematic) are easily identified,
one cannot exclude other intervening phases. We will explore the finite-
temperature phase diagram of the model using directed-loop Monte Carlo
simulation [34, 35, 36, 37, 38, 39], which allows for efficient non-local
moves. The precise implementation we use is presented in Appendix A. The
simulations are supplemented by a field-theoretical analysis in terms of a
Coulomb gas description of the system (Sec. IV).

##  III Physical Observables

In this section, we describe the observables measured during the Monte Carlo
simulations to characterize the phases and the transitions.

Winding number fluctuations —  Fully-packed loop configurations on the square
lattice can be associated with two integer winding numbers $W_{x}$ and
$W_{y}$. To compute $W_{y}$ ($W_{x}$), draw a horizontal (vertical) line that
cuts across $L$ lattice bonds oriented in the $y$ ($x$) direction. For a given
configuration, we denote by $N_{o}$ and $N_{e}$ the number of loop segments on
the odd and even bonds that cross this line. The winding numbers are defined
as $N_{e}-N_{o}$. Each winding number $W_{x}$ and $W_{y}$ vary between $-L$
and $L$, and there is at least one fully-packed loop configuration for any
pair $(W_{x},W_{y})$ in this range. Note that the loop constraint ensures that
the winding numbers calculated using different parallel lines are the same.

On account of translation symmetry, the equilibrium average values of
$W_{x},W_{y}$ vanish, but not their fluctuations

$\langle W^{2}\rangle=\frac{1}{2}\langle W^{2}_{x}+W^{2}_{y}\rangle,$ (2)

which have useful physical content and can easily be measured in Monte Carlo
simulations [42, 35, 55, 73, 74]. Low-temperature order parameter —  We can
identify the low-temperature phase through the rotational symmetry breaking,
nematic order parameter [75, 35, 43, 54, 53]

$D=\frac{1}{N}|N_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}-N_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}|,$ (3)

with $N_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}=\sum_{{\bf r}}n_{\hbox{\hskip
1.42271pt{\vrule height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf r})$ and
$N_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}=\sum_{{\bf
r}}n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf r})$,
where $n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf r})$
denotes a horizontal loop segment at the site ${\bf r}$. It is 1 if a loop
segment occupies the edge between ${\bf r}$ and ${\bf r}+(1,0)$ and is 0 if
the edge is empty. $n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf r})$ denotes a vertical
loop segment at site at lattice site ${\bf r}$, and is 1 if a loop segment
occupies the edge between ${\bf r}$ and (${\bf r}+(0,1)$). The order parameter
$D$ is $1$ in the two nematic ground-states, and vanishes ($\langle
D\rangle=0$) at infinite temperature.

We also compute the associated susceptibility [75, 35, 43]:

$\chi_{D}=N(\langle D^{2}\rangle-\langle D\rangle^{2})$ (4)

and monitor its temperature dependence. As shown below, the divergence of
$\chi_{D}$ allows us to determine the transition temperature and the form of
the divergence can be further used to infer the nature of the transition [76,
77, 78, 79].

Loop-segment (dimer) correlators —  We consider the connected correlation
function [80, 81, 82, 35, 41] between loop-segments separated by a vector
${\bf r}=(x,y)$: $C_{\alpha,\beta}({\bf r})=\langle n_{\alpha}({\bf
0})n_{\beta}({\bf r})\rangle-1/4$, where $\alpha,\beta$ can be
${\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}},{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}$. The expectation value
$\langle n_{\alpha}\rangle\langle n_{\beta}\rangle=1/4$ has been subtracted to
get the connected correlator. In the Monte Carlo simulations, we average over
all possible initial positions ${\bf 0}$ of the first loop segment, as well as
all equivalent pairs $\alpha,\beta$. For simplicity, we will focus on the
lattice direction ${\bf r}=(x=r,0)$ and consider three types of loop-segment
correlations, longitudinal, transverse, and crossed, respectively defined as:

$\displaystyle C^{L}(r)$ $\displaystyle=$ $\displaystyle\langle
n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
0})n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(r,0)\rangle-1/4,$ (5) $\displaystyle C^{T}(r)$ $\displaystyle=$
$\displaystyle\langle n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf 0})n_{\hbox{\hskip
1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}(r,0)\rangle-1/4,$ (6)
$\displaystyle C^{C}(r)$ $\displaystyle=$ $\displaystyle\langle
n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
0})n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}(r,0)\rangle-1/4.$ (7)

We will also consider correlators associated to the nematic order parameter

$\langle D_{\bf 0}D_{\bf r}\rangle=\langle(n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}(0)-n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(0))(n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf
r})-n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
r}))\rangle\\\ =C^{L}(r)+C^{T}(r)-2C^{C}(r)$ (8)

where in the last line we again focus on the direction ${\bf r}=(x=r,0)$.

Monomer correlators —  We also measure the monomer-monomer correlator [83,
84], which requires going beyond the definition of the fully-packed loop
configurations space by allowing two test monomers – sites touched by only one
dimer – while the rest of the sites are all touched by two loop segments. No
monomer is included in the fully-packed loop model, and the Monte Carlo
configurations generated by the directed loop algorithm, once a directed loop
is finished, do not contain monomers. However, during the intermediate steps
of the Monte Carlo process, the directed loop algorithm precisely samples the
extended phase space with monomers, allowing to sample the monomer correlator
defined below (see Ref. [42, 35, 36] and Appendix A for details). We define
monomer correlation function

$M({\bf r})=\langle m({\bf 0})m({\bf r})\rangle.$ (9)

where the presence of a monomer at site ${\bf r}$ is denoted by $m({\bf r})=1$
and $m({\bf r})=0$ otherwise. The monomer correlation function $M({\bf r})$ is
estimated as the fraction of such Monte Carlo samples of configurations with
two monomers separated by ${\bf r}$ [84, 36, 43].

##  IV Theoretical framework

Before analyzing the results of the numerical simulations, we first describe
the finite-temperature phase diagram of the model using a field theoretical
analysis. Following its success in two-dimensional models of statistical
mechanics [85, 86, 87, 88, 89, 90, 40, 91, 41, 35, 43, 92], we use a Coulomb
gas description [24], formulated in terms of a height field $h(p)$. This field
lives on plaquettes $p$ of the lattice, and is defined (up to an irrelevant
constant) in the following way: when turning clockwise around A-sublattice
sites of the square lattice, the height increases (decreases) by $1/2$ (i.e.
$h\rightarrow h\pm 1/2$) if one crosses a loop segment (an empty edge). At the
microscopic level, it can be shown that the value of the height field inside a
small patch, can be changed by $1$ without changing the local loop segment
configuration, simply by a change of configuration at far away points (similar
argument as for the dimer model [35]). To see this, consider a region
encircled by a pair of loops separated by a plaquette (). Local changes to
these loops that convert them into a single zig-zag loop () surrounding the
region changes the heights by 1 (holding heights on the exterior fixed)
everywhere inside the region and irrespective of the distance from these
loops.

This indicates that the physical action should be invariant under height
shifts $h\rightarrow h\pm 1$. Promoting the height field to the continuum
$h({\bf r)}$, we expect the effective action to be:

$S=\int d^{2}r\,[g(T)\pi(\nabla h({\bf r}))^{2}+v\cos(4\pi h({\bf r}))].$ (10)

Such a free compact boson model can be used to describe several two-
dimensional statistical physics models, including the Tomonaga–Luttinger
liquid (through bosonization, see e.g [93]) and the XXZ spin chain (see e.g.
[94]). We briefly justify this action below, and discuss its validity
alongside the numerical results in Sec. V.

This action is of the sine-Gordon type [95, 96]. Here $g(T)$ is the Coulomb
gas coupling constant, which depends on microscopic details and on temperature
$T$. At infinite temperature, we have $g(T=\infty)=1/3$ from exact results for
the six-vertex model at the ice point. 111The ice-point corresponds to the XXZ
spin chain at $\Delta=1/2$ in notations where $\Delta>0$ corresponds to
ferromagnetic interactions. The value of $g$ is given by $\Delta=\cos(\pi g)$,
see e.g. Ref. [94, 104].. This action displays the competition between the
first term ($g(T)\pi(\nabla h({\bf r}))^{2}$) which alone describes the
critical phase (rough in the height language) to be encountered at high
temperature, and the $v\cos(4\pi h({\bf r}))$ “vertex” term whose minima
corresponds to the two nematic configurations for which the average height is
constant (flat configurations) and takes values $\bar{h}=\pm 1/4$. $v$ can
also depend on temperature but its exact dependency is not relevant as long as
it remains positive such that the two nematic configurations are always
favored.

In the Coulomb gas language and given the periodicity of the height
$h\rightarrow h+1$ in the microscopic configurations, the later vertex term
can be identified with an electric charge $e=2$ operator. This term is
irrelevant at infinite temperature where $g=1/3$, but becomes relevant when
$g\geq g_{c}=1$ (a general electric charge $e$ operator reads $\exp(i2e\pi
h)$, and has scaling dimension $e^{2}/(2g)$, and thus becomes relevant when
$g\geq e^{2}/4$). As interactions favor the flat nematic phases, we expect $g$
to increase (from its $g(T=\infty)=1/3$ value) as the temperature is lowered.

The Coulomb gas analysis predicts a Kosterlitz-Thouless phase transition [98,
99, 95, 96] from a high-temperature critical phase to the low-temperature
nematic phase, and furthermore provides predictions for several observables.
First, the winding fluctuations can be related to the Coulomb gas constant
[35]:

$\langle W^{2}\rangle=\sum_{n\in\mathbb{Z}}n^{2}e^{-g\pi
n^{2}}/\sum_{n\in\mathbb{Z}}e^{-g\pi n^{2}},$ (11)

as used in Fig. 2 below. This allows in particular to extract the Kosterlitz-
Thouless transition temperature $T_{KT}$ at the predicted critical Coulomb gas
constant $g_{c}=1$.

Next, the leading terms for the dimer/loop segment occupation operator in the
continuum have been identified in Ref. [41] as:

$\displaystyle n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
r}=(x,y))$ $\displaystyle=$
$\displaystyle\frac{1}{2}+(-)^{x+y+1}\nabla_{y}h-\frac{X}{2i}(\exp(2i\pi
h({\bf r}))-\exp(-2i\pi h({\bf r})))$ (12) $\displaystyle n_{\hbox{\hskip
1.42271pt{\vrule height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf
r}=(x,y))$ $\displaystyle=$
$\displaystyle\frac{1}{2}+(-)^{x+y}\nabla_{x}h+\frac{X}{2i}(\exp(2i\pi h({\bf
r}))-\exp(-2i\pi h({\bf r})))$ (13)

The loop segment occupation is thus composed of a gradient part and a vertex
part. The vertex part of the loop segment operator can be expressed in
harmonics of $2\pi h$ (as $h\equiv h+1$) and microscopic $\pi/2$ rotations of
the model give $h\rightarrow-h$ and $h\rightarrow h+1/2$. It can be identified
with an electric charge $e=1$ in the Coulomb gas.

Note that the overall sign in front of the gradient depends on the convention
for the height (odd or even sublattice). The constant $X$ cannot be fixed
easily and we need an external exact solution (see below) – in fact, we expect
it to be re-normalized, that is to change with temperature. This gives the
following predictions for the leading terms of the correlators defined in Eqs.
(5)-(7)

$\displaystyle C^{L}({\bf r}=(x,y))$ $\displaystyle=$ $\displaystyle\langle
n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(0)n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
r})\rangle-1/4=(-)^{x+y}A\frac{x^{2}-y^{2}}{(x^{2}+y^{2})^{2}}+\frac{B}{(x^{2}+y^{2})^{1/2g}}$
(14) $\displaystyle C^{T}({\bf r}=(x,y))$ $\displaystyle=$
$\displaystyle\langle n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}(0)n_{\hbox{\hskip
1.42271pt{\vrule height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf
r})\rangle-1/4=(-)^{x+y}A\frac{y^{2}-x^{2}}{(x^{2}+y^{2})^{2}}+\frac{B}{(x^{2}+y^{2})^{1/2g}}$
(15) $\displaystyle C^{C}({\bf r}=(x,y))$ $\displaystyle=$
$\displaystyle\langle n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(0)n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf
r})\rangle-1/4=(-)^{x+y}A\frac{2xy}{(x^{2}+y^{2})^{2}}-\frac{B}{(x^{2}+y^{2})^{1/2g}}$
(16)

The coefficient $A=\frac{1}{4g\pi^{2}}$ is fixed by the operator product
expansion (12) and the two-point correlation function $\langle
h(x)h(y)\rangle$ known exactly for free compact boson conformal field theory
[100]. We also have $B=X^{2}/2$, however, its dependence on $g$ is not
universal. At $T=\infty$, exact expressions for the XXZ spin chain [94] give
$B\simeq 0.01795$, see Table 1 in Ref. [94] (see also Ref. [82]).

Finally, we note that on the lattice, a monomer creates a dislocation of $\pm
1$ in the height field. The prediction of the monomer correlator decaying as

$M(r)\propto r^{-g}$ (17)

follows [35, 43] from the identification of the monomer operator with the
$m=\pm 1$ magnetic charge operator (the sign depends on the sublattice) with a
scaling dimension $gm^{2}/2$.

This interpretation parallels the one for the interacting classical dimer
model [42, 35, 43] with the following three minor (albeit important for
numerics) distinctions: (i) The infinite temperature value of the Coulomb gas
constant $g=1/3$ renders the vertex contribution (scaling as $r^{-3}$)
subleading with respect to the dipolar contribution (scaling as $r^{-2}$),
which explains why it is often not reported in the polarization fluctuations
for the 6-vertex model [80]. For the dimer problem, we have $g(T=\infty)=1/2$
and both terms contribute equally to the $r^{-2}$ decay of the dimer
correlators [82]. (ii) The critical value of the Coulomb gas constant at the
critical point is $g_{c}=1$ (instead of $g_{c}=4$ for the dimer model),
consistent with the lower degeneracy of the ground-states ($2$ nematic ground-
states instead of $4$ columnar ground-states for the dimer model) and
resulting in a larger value of the anomalous dimension of the low-temperature
order parameter $\eta_{D}={1/g_{c}}=1$ for the loop model (see below) instead
of $\eta_{D}=1/g_{c}=1/4$ for the dimer model at their respective KT
transitions [42]. (iii) Lastly, we found the critical value for the winding
fluctuations $\langle W^{2}\rangle$ is much larger for the loop model, which
allows for a statistically meaningful measurement in the Monte Carlo
simulations. The very small value of $\langle W^{2}\rangle$ for the dimer
model does not allow for an accurate Monte Carlo determination of the critical
point using the value of the winding number fluctuations.

Figure 2: (a) MC results for winding number fluctuations as a function of $T$.
The gray dashed line shows the critical winding number fluctuations $\langle
W^{2}\rangle_{c}=0.07958$, which is obtained from Eq. (11) with $g_{c}=1$ at
the transition point. Inset is a zoom-in for the $1.33\leqslant T\leqslant
1.41$ region. (b) Finite-size scaling for the estimated transition temperature
as a function of system size. The finite-size $T_{KT}(L)$ data points are
obtained from (a). The black curve shows the fit to Eq. (18). The
extrapolation to the thermodynamic limit gives $T_{KT}=1.425(1)$. The inset
shows $\langle W^{2}\rangle$ as a function of $g$ according to the relation in
Eq. (11). Figure 3: (a) The nematic order parameter $\langle D\rangle$ as a
function of temperature for different system sizes. (b) Susceptibility
$\chi_{D}$ of the nematic order parameter defined in Eq. (4). Data collapse is
performed in (c) the critical phase with $T>T_{KT}$. Here we use $\eta_{D}=1$
and $T_{KT}=1.425$, and all the data point nicely collapse onto a single
curve.

##  V MC simulation results

We present our MC simulation results in this section. It contains results for
observables from which we can precisely estimate the critical temperature
$T_{KT}$: the winding number fluctuations (Sec. V.1) and the essential
singularity of the nematic susceptibility in the KT transition (Sec. V.2).
Section V.3 presents results for different correlation functions in the high-
temperature critical phase, confirming the field theoretical analysis
presented in Sec. IV.

###  V.1 Winding number fluctuations

The numerical results for the winding number fluctuations $\langle
W^{2}\rangle$ as a function of temperature $T$ of the classical loop model are
shown in Fig. 2 (a). We simulate system sizes up to $L=128$ for this
measurement. These data directly provide the temperature dependence (albeit on
finite size) of the Coulomb gas constant, which will later be compared with
other estimates of $g(T)$. At the transition point, the analysis of Sec. IV
predicts the critical Coulomb gas constant $g_{c}$ to be $1$, corresponding to
the critical winding number fluctuations $\langle W^{2}\rangle_{c}=0.07958$
[from Eq.( 11), also see inset of Fig. 2(b)]. The predicted critical value
$\langle W^{2}\rangle_{c}$ is shown as the gray dashed horizontal line in Fig.
2 (a) and in its inset.

We estimate the transition temperature $T_{KT}(L)$ for each system size as the
temperature at which the winding number fluctuations cross the critical value,
which is in turn estimated from a linear fit of the data points near $\langle
W^{2}\rangle_{c}$ [Fig. 2 (a)inset]. This estimate has an obvious finite-size
dependence. To determine the transition temperature $T_{KT}$ in the
thermodynamic limit, we use the following finite-size scaling relation for a
KT transition [101, 102, 92]:

$\frac{1}{T_{KT}(L)}=\frac{1}{T_{KT}}+\frac{C}{\log(L/L_{0})^{2}},$ (18)

where $C$ is a constant. By fitting the estimated $T_{KT}(L)$ in Fig. 2 (b)
with Eq. (18), we obtain $T_{KT}=1.425(1)$.

Figure 4: Equal-time loop-segment correlation functions (a) $C^{L}+C^{T}$ as a
function of $r$, (b) the log-log plot for $|C^{L}+C^{T}|$ (absolute value is
used to correct for very small negative values occurring at large $r$, large
$T$ due to statistical fluctuations caused by the finite Monte Carlo
sampling). (c) The crossed correlations $C^{C}$ (Eq. (7)), and (d) The log-log
plot of the monomer correlations $\langle M(r)\rangle$ in Eq. (9). The system
size is $L=256$ for the loop-segment correlators and $L=400$ for monomer
correlations. These data correspond to the high-temperature critical phase,
that is temperatures above the estimated $T_{KT}=1.425$. Gray curves are
power-law fits (in their respective fitting range) according to the scaling
form $B^{\prime}/r^{1/g}+C$ for $C^{L}+C^{T}$,
$(-1)^{r}A^{\prime}/r^{\alpha_{S}}+B^{\prime}/r^{\alpha_{U}}+C$ for $C^{C}$,
and $B^{\prime}/r^{g}+C$ for $M(r)$. In all cases, we add a constant to the
power-law fits to account for a small non-vanishing value of correlators at
large-distance in our finite-size Monte Carlo simulations. Figure 5:  (a) The
equal-time loop-segment correlation functions of $(-)^{r}(C^{L}-C^{T})$ for
$L=256$. A weak staggered part can be observed in this representation
(particularly visible for the lowest temperatures at short distances),
signaling a small uniform component for $(C^{L}-C^{T})$. The gray curves fit
$(C^{L}-C^{T})$ to the form
$(-1)^{r}A^{\prime}/r^{\alpha_{S}}+B^{\prime}/r^{\alpha_{U}}+C$. (b) The two
exponents $\alpha_{S}$ and $\alpha_{U}$ obtained from this fit. Figure 6: The
Coulomb gas constant $g$ obtained from Fig. 4 (b), (c), (d), and $\langle
W^{2}\rangle$. The vertical gray dash line indicates the transition point
$\beta_{KT}=1/T_{KT}\simeq 0.7$. The upper gray dash line in the horizontal
direction denotes the critical value $g_{c}=1$; the middle one at $g=1/2$
corresponds to the case where the vertex and dipolar terms have the same
contribution to scaling, the lowest gray dash lines indicates the infinite
temperature value $g=1/3$.

###  V.2 Nematic order parameter and susceptibility

Our results for the finite-temperature behavior of the nematic order parameter
$\langle D\rangle$ defined in Eq.( 3) are presented in Fig. 3 (a) for
different system sizes. We clearly observe the existence of a nematic phase at
low temperature, where $\langle D\rangle$ takes a finite value, and a high-
temperature phase where $\langle D\rangle$ vanishes relatively slowly as
system sizes increase. The associated nematic susceptibility $\chi_{D}$ [Eq.
(4)], represented in Fig. 3 (b), shows a clear diverging peak (with system
size) in the temperature range where $\langle D\rangle$ starts to vanish.

There is a clear shift in the temperature of the peak in $\chi_{D}$ as system
size varies, as this can be used to perform a data collapse in order to cross-
validate the transition temperature $T_{KT}$ obtained in the previous section.
In the vicinity of the KT transition, we indeed expect that the susceptibility
$\chi_{D}$ obeys the scaling behavior [99]

$\chi_{D}\sim L^{2-\eta_{D}}f\left(L\exp(-\frac{K}{\sqrt{T-T_{KT}}})\right)$
(19)

for $T>T_{KT}$ where $K$ is a constant and $\eta_{D}=1/g_{c}=1$ the anomalous
dimension [42]. Such data collapse has been used in the literature to
determine the Kosterlitz-Thouless transition temperature $T_{KT}$ in many 2D
systems such as the 2D XY model [101], magnetic thin films [102], triangular
lattice transverse field Ising model [103], or for the pairing transition in
various 2D fermionic lattice models [76, 77, 78, 79]. We use the data in the
$T>T_{KT}$ region to rescale the $y$ axis as $\chi_{D}L^{-(2-\eta_{D})}$ and
the $x$ axis as $L\exp(-\frac{K}{\sqrt{T-T_{KT}}})$ as shown in Fig. 3 (c). We
obtain that values $T_{KT}=1.425,\eta_{D}=1$ provide a good data collapse,
resulting in a good agreement with the $T_{KT}$ obtained in Fig. 2 (b).

###  V.3 Correlation functions

The height description of the loop-segment correlations given in Eqs. (14),
(15) and (16) suggests that the correlators have two contributions - from the
vertex and the dipolar part [41, 94, 82, 36]. We show the correlators
calculated from the Monte Carlo methods for $L=256$ system and present the
fits to their expected forms in Appendix B. Below we present a simpler
approach to fitting - by considering combinations of the correlators that
separate out the vertex and dipolar terms.

We first consider the sum $C^{L}+C^{T}=\langle n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(0)n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
r})\rangle+\langle n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}(0)n_{\hbox{\hskip
1.42271pt{\vrule height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf
r})\rangle-\frac{1}{2}$ which should contain only a vertex contribution
$2B/(x^{2}+y^{2})^{1/2g}$. This combination for the direction ${\bf r}=(r,0)$
is shown in Figs. 4(a) and (b) in linear and log scales respectively.
Consistent with the expected form, the combination shows a power-law scaling
with distance $r$ with an exponent (slope in the log plot) that increases with
temperature. The estimated value of $g$ from this combination is discussed
further below.

We then consider the crossed correlators $C^{C}=\langle
n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(0)n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf r})\rangle-1/2$ along the
direction ${\bf r}\equiv(r,0)$; where the correlation is expected to be
dominated by the vertex term according to Eq. 16. Monte Carlo estimates of
$C^{C}(r)$ for different temperatures are shown in Fig. 4 (c). The crossed
correlators show expected power law scaling at large distances but with a
possible oscillatory subleading correction that affects the short distance
correlations, which is visible at higher temperatures. We tentatively
attribute this effect to further subleading terms that do not cancel for $y=0$
and are not included in Eq. 16. This is confirmed by a fit to the form
$(-1)^{r}A^{\prime}/r^{\alpha_{S}}+B^{\prime}/r^{\alpha_{U}}+C$ for $C^{C}(r)$
in Fig. 4 (c), where we find that the amplitude $A^{\prime}$ of the
oscillating term is always small, and almost vanishing as the temperature
lowers towards the critical point (e.g. $A^{\prime}\sim 0.03$ for $T=4$ and
$A^{\prime}\sim 0.002$ for $T=1.3$). The results of this fit for $\alpha_{U}$
allow to estimate $g$ as obtained from Eq. 16 and as presented below.

Next, we study the combination $C^{L}-C^{T}=\langle
n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip
1.77829pt}}}}(0)n_{\hbox{\vbox{\vbox{{\hrule
height=1.42271pt,width=5.69046pt,depth=0.0pt}\vskip 1.77829pt}}}}({\bf
r})\rangle-\langle n_{\hbox{\hskip 1.42271pt{\vrule
height=5.69046pt,width=1.42271pt,depth=0.0pt}}}(0)n_{\hbox{\hskip
1.42271pt{\vrule height=5.69046pt,width=1.42271pt,depth=0.0pt}}}({\bf
r})\rangle$ which, based on Eqs. (14) and (15), is expected to have a purely
staggered dipolar contribution $(-)^{r}/r^{2}$. Our results for
$(-)^{r}(C^{L}-C^{T})({\bf r})$ in the ${\bf r}=(r,0)$ direction are presented
in Fig. 5 (a). We observe that there is a small but non-vanishing uniform
component in the numerical data (which appears as a staggered part in Fig. 5
(a) due to the $(-)^{r}$ factor). To account for this, we fit $C^{L}-C^{T}$ to
a form $(-1)^{r}A^{\prime}/r^{\alpha_{S}}+B^{\prime}/r^{\alpha_{U}}+C$. The
constant $C$ accounts for a non-zero value of this correlator present only at
temperatures close to the phase transition, which we attribute to the finite
sizes used in our Monte Carlo simulations. Here $\alpha_{U}$ is meant to
describe a subleading correction to the vertex part not included in Eqs. (14)
and (15). The estimates for $\alpha_{S},\alpha_{U}$ are presented in Fig. 5
(b), where we find that $\alpha_{S}$ is very close to the predicted value $2$
all along the high-temperature critical phase, and $\alpha_{U}>\alpha_{S}$
confirming the subleading nature of this uniform correction. If we fit the
data fixing $\alpha_{U}$ to be zero, the exponent $\alpha_{S}$ is always
larger than its expected value of $2$.

Finally, the monomer-monomer correlator $M(r)$ in Eq. (9) should decay only
with the vertex contribution $1/r^{g}$ [104, 42, 36, 55]. We present the
monomer correlations $M(r)$ at different temperatures in Fig. 4 (d) for
$L=400$ (within the directed loop algorithm, we can get good statistics for
$M(r)$ for larger systems than for loop-segment correlators). The log-log plot
shows a clear power-law decay above the Kosterlitz-Thouless transition
temperature.

We now collect, in Fig. 6, the estimates of the Coulomb gas constant obtained
from the fits to the correlators $C^{L}+C^{T}$ (from Fig. 4(a) and (b)),
$C^{C}$ (from Fig. 4(c)) and $M(r)$ (from Fig. 4(d)) as well as from the
winding number fluctuations $\langle W^{2}\rangle$ in Fig. 2. Fig. 6 shows the
temperature dependence of $g$ as a function of inverse temperature
$\beta=1/T$. We find that as $\beta$ increases (i.e. as the temperature
decreases), the Coulomb gas constant increases from its infinite temperature
value $g=1/3$, which is consistent with the expectation that attractive
interactions tend to stiffen the loops. As the temperature decreases from the
$T=\infty$ point to finite but high temperature, the dipolar part in Eqs.
(14), (15) and (16) dominates down to a temperature $T\approx 3$, below which
$g>1/2$ and the vertex part takes over down to $T_{KT}\simeq 1.425$
$(\beta_{KT}\simeq 0.7)$ where $g=1$. The various estimates of $g$ are overall
in good agreement (we note the $g$ value obtained from $C^{C}$ is less
accurate due to the subleading oscillations at high temperature) with each
other, and consistent with the theoretical expectations of Sec. IV.

##  VI Discussion and conclusions

In this work, we investigated the finite-temperature phase diagram of a
classical model of fully-packed loops on the square lattice with attractive
local interactions between loop segments. With the help of a directed-loop
Monte Carlo algorithm and a field theoretical analysis based on a height
description of loop configurations, we are able to locate the finite
temperature Kosterlitz-Thouless transition, separating a critical phase at
$T>T_{KT}$ and a nematic phase below $T_{KT}$. We find that in the loop model
the anomalous dimension at the KT transition $\eta_{D}=1$, is four times
larger than that in the classical dimer model [42]. The high-temperature
critical phase is fully characterized by the temperature dependence of the
Coulomb gas constant presented in Fig. 6, which is obtained using several
different concurrent estimates.

An interesting, closely related system to consider would be a similar
classical model, but with repulsive interactions ($V>0$) between fully packed
loops, favoring large-winding sectors. Analogous repulsive interactions in the
dimer model result in a continuous phase transition from a critical to
staggered phase, which has been argued to be in the two-dimensional Ising
universality class [92].

We connect our results to the quantum loop model on the square lattice. From
our analysis, we expect that the QLM on the square lattice should also host a
critical phase at any sufficiently high-temperature parametrized by a Coulomb
gas constant $g(T/t,V/t)$ which depends on temperature and potential energy,
similar to the quantum dimer model [55]. In general, a high-temperature
critical phase can be found in the constrained entropic scaling regime of
several strongly constrained quantum systems, which can extend down to low
temperatures (see Ref. [105] for an extensive discussion). At large negative
ratio of potential to kinetic energy ($V/t\ll 0$), the QLM hosts a nematic
ground-state. From our results, we conclude that the finite-temperature phase
transition to the nematic phase in the QLM should occur as a Kosterlitz-phase
transition that can be described using the same analysis provided here. The
QLM also hosts a plaquette ground-state in a finite range of $-0.35\lesssim
V/t\ll 1$ [53]. We believe that the finite-temperature phase transition to
this plaquette phase should be of KT type too, with an effective action
described by Eq. (10) but with negative $v$, as the two plaquette ground-
states have average height $\bar{h}=0,\frac{1}{2}$. It would be interesting to
find a classical model with a similar phase transition and low-temperature
phase. Finally, we note that the directed loop algorithm that we use can be
directly implemented as a new move [55] within the sweeping cluster algorithm
[106, 107] for the QLM, allowing the study of its finite-temperature phase
diagram fully taking into account the loop constraints and winding
fluctuations.

Rydberg atom arrays form a new type of platform where constraints (due to the
Rydberg blockade) play an important role to determine the ground-state phase
diagram, with a rich variety of phases observed [62, 64, 65, 108, 5, 7]. To
the best of our knowledge, the finite-temperature phase transitions out of
these phases has not been studied experimentally so far. It would be
interesting to see where the finite-temperature critical phase that we find
here could be relevant in some experimental regimes where the fully-packed
constraint is a relevant approximation in Rydberg atom arrays.

Acknowledgments — We acknowledge support from the ANR/RGC Joint Research
Scheme sponsored by Research Grants Council of Hong Kong SAR of China (Project
No. A_HKU703/22) and French National Research Agency (grant
ANR-22-CE30-0042-01). XXR, ZY and ZYM further acknowledge the support from the
Research Grants Council of Hong Kong SAR of China (Project Nos. 17301420,
17301721, AoE/P-701/20, 17309822, HKU C7037-22G), and BD, GJS and FA the
support from the joint PhD program between CNRS and IISER Pune, as well as the
grant NanoX ANR-17-EURE-0009 in the framework of the French “Programme des
Investissements d’Avenir”. The research of JR is supported by the Huawei Young
Talents Program at IHES. We acknowledge the use of HPC resources from CALMIP
(grants 2022-P0677 and 2023-P0677), GENCI (projects A0110500225 and
A0130500225), the HPC2021 system under the Information Technology Services,
the Blackbody high-performance computing system at the Department of Physics,
University of Hong Kong and Param Brahma computing facility at IISER Pune.”.
GJS and BD thank K. Damle for useful discussions as well as TIFR, Mumbai for
hospitality during the completion of this work.

## References

  * Moessner and Sondhi [2001a] R. Moessner and S. L. Sondhi, Phys. Rev. Lett. 86, 1881 (2001a).
  * Moessner and Sondhi [2001b] R. Moessner and S. L. Sondhi, Phys. Rev. B 63, 224401 (2001b).
  * Bramwell and Gingras [2001] S. T. Bramwell and M. J. P. Gingras, Science 294, 1495 (2001), https://www.science.org/doi/pdf/10.1126/science.1064761 .
  * 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, Science 374, 1242 (2021).
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Proceedings of the National Academy of Sciences 118, e2015785118 (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, Nature 595, 227 (2021).
  * Yan  _et al._ [2022] Z. Yan, R. Samajdar, Y.-C. Wang, S. Sachdev, and Z. Y. Meng, Nat. Commun. 13, 5799 (2022).
  * Rokhsar and Kivelson [1988] D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. 61, 2376 (1988).
  * Blunt  _et al._ [2008] M. O. Blunt, J. C. Russell, M. d. C. Giménez-López, J. P. Garrahan, X. Lin, M. Schröder, N. R. Champness, and P. H. Beton, Science 322, 1077 (2008).
  * Gruzberg  _et al._ [1999] I. A. Gruzberg, A. W. W. Ludwig, and N. Read, Phys. Rev. Lett. 82, 4524 (1999).
  * Read and Saleur [2001] N. Read and H. Saleur, Nuclear Physics B 613, 409 (2001).
  * Fendley [2008] P. Fendley, Annals of Physics 323, 3113 (2008).
  * Nahum  _et al._ [2011] A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Phys. Rev. Lett. 107, 110601 (2011).
  * Nahum  _et al._ [2015] A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Phys. Rev. X 5, 041048 (2015).
  * Alet  _et al._ [2006a] F. Alet, G. Misguich, V. Pasquier, R. Moessner, and J. L. Jacobsen, Phys. Rev. Lett. 97, 030403 (2006a).
  * Powell and Chalker [2008] S. Powell and J. T. Chalker, Phys. Rev. Lett. 101, 155702 (2008).
  * Charrier  _et al._ [2008] D. Charrier, F. Alet, and P. Pujol, Phys. Rev. Lett. 101, 167205 (2008).
  * Powell and Chalker [2009] S. Powell and J. T. Chalker, Phys. Rev. B 80, 134413 (2009).
  * Chen  _et al._ [2009] G. Chen, J. Gukelberger, S. Trebst, F. Alet, and L. Balents, Phys. Rev. B 80, 045112 (2009).
  * Sreejith and Powell [2015] G. J. Sreejith and S. Powell, Phys. Rev. B 92, 184413 (2015).
  * Sreejith  _et al._ [2019] G. J. Sreejith, S. Powell, and A. Nahum, Phys. Rev. Lett. 122, 080601 (2019).
  * Sreejith and Powell [2014] G. J. Sreejith and S. Powell, Phys. Rev. B 89, 014404 (2014).
  * Nienhuis [2010] B. Nienhuis, in  _Exact methods in low-dimensional statistical physics and quantum computing_ , edited by J. Jacobsen, S. Ouvry, V. Pasquier, D. Serban, and L. F. Cugliandolo (Oxford University Press, Oxford, 2010) Chap. 6, pp. 159–195.
  * Nienhuis [1987] B. Nienhuis, in  _Phase transitions and critical phenomena Vol. 11_ , edited by C. Domb and J. Lebowitz (Academic Press, 1987) Chap. 1, pp. 1–53.
  * Jacobsen [2009] J. L. Jacobsen, Conformal field theory applied to loop models, in  _Polygons, Polyominoes and Polycubes_, edited by A. J. Guttman (Springer Netherlands, Dordrecht, 2009) pp. 347–424.
  * Nahum  _et al._ [2013] A. Nahum, P. Serna, A. M. Somoza, and M. Ortuño, Phys. Rev. B 87, 184204 (2013).
  * Kondev  _et al._ [1996] J. Kondev, J. de Gier, and B. Nienhuis, Journal of Physics A: Mathematical and General 29, 6489 (1996).
  * Temperley and Lieb [1971] H. N. V. Temperley and E. H. Lieb, Proc. Roy. Soc. Lond. A 322, 251 (1971).
  * de Gennes [1972] P. de Gennes, Physics Letters A 38, 339 (1972).
  * Jacobsen  _et al._ [2003] J. L. Jacobsen, N. Read, and H. Saleur, Phys. Rev. Lett. 90, 090601 (2003).
  * Cardy [2005] J. Cardy, Annals of Physics 318, 81 (2005), special Issue.
  * Jacobsen and Alet [2009] J. L. Jacobsen and F. Alet, Phys. Rev. Lett. 102, 145702 (2009).
  * Schwandt  _et al._ [2010] D. Schwandt, M. Mambrini, and D. Poilblanc, Phys. Rev. B 81, 214413 (2010).
  * Barkema and Newman [1998] G. Barkema and M. Newman, Physical Review E 57, 1155 (1998).
  * Alet  _et al._ [2006b] F. Alet, Y. Ikhlef, J. L. Jacobsen, G. Misguich, and V. Pasquier, Phys. Rev. E 74, 041124 (2006b).
  * Sandvik and Moessner [2006] A. W. Sandvik and R. Moessner, Physical Review B 73, 144504 (2006).
  * Syljuåsen and Sandvik [2002] O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
  * Syljuåsen and Zvonarev [2004] O. F. Syljuåsen and M. Zvonarev, Physical Review E 70, 016118 (2004).
  * Alet and Sørensen [2003] F. Alet and E. S. Sørensen, Phys. Rev. E 68, 026702 (2003).
  * Kondev and Henley [1996] J. Kondev and C. L. Henley, Nuclear Physics B 464, 540 (1996).
  * Moessner  _et al._ [2004] R. Moessner, O. Tchernyshyov, and S. L. Sondhi, Journal of Statistical Physics 116, 755 (2004).
  * Alet  _et al._ [2005] F. Alet, J. L. Jacobsen, G. Misguich, V. Pasquier, F. Mila, and M. Troyer, Phys. Rev. Lett. 94, 235702 (2005).
  * Papanikolaou  _et al._ [2007] S. Papanikolaou, E. Luijten, and E. Fradkin, Physical Review B 76, 134514 (2007).
  * Kundu and Damle [2023] S. Kundu and K. Damle, Flux fractionalization transition in two-dimensional dimer-loop models (2023), arXiv:2305.07012 [cond-mat.stat-mech] .
  * Castelnovo  _et al._ [2005] C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, Annals of Physics 318, 316 (2005).
  * Balasubramanian  _et al._ [2022] S. Balasubramanian, V. Galitski, and A. Vishwanath, Phys. Rev. B 106, 195127 (2022).
  * Castelnovo  _et al._ [2007] C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, Annals of Physics 322, 903 (2007).
  * Henry and Roscilde [2014] L.-P. Henry and T. Roscilde, Phys. Rev. Lett. 113, 027204 (2014).
  * Shannon  _et al._ [2004] N. Shannon, G. Misguich, and K. Penc, Phys. Rev. B 69, 220403 (2004).
  * Syljuåsen and Chakravarty [2006] O. F. Syljuåsen and S. Chakravarty, Phys. Rev. Lett. 96, 147004 (2006).
  * Plat  _et al._ [2015] X. Plat, F. Alet, S. Capponi, and K. Totsuka, Phys. Rev. B 92, 174402 (2015).
  * Roychowdhury  _et al._ [2015] K. Roychowdhury, S. Bhattacharjee, and F. Pollmann, Phys. Rev. B 92, 075141 (2015).
  * Ran  _et al._ [2023] X. Ran, Z. Yan, Y.-C. Wang, J. Rong, Y. Qi, and Z. Y. Meng, Phys. Rev. B 107, 125134 (2023).
  * [54] Z. Yan, X. Ran, Y.-C. Wang, R. Samajdar, J. Rong, S. Sachdev, Y. Qi, and Z. Y. Meng,  arXiv:2205.04472 (2022) [cond-mat.str-el] .
  * Dabholkar  _et al._ [2022] B. Dabholkar, G. J. Sreejith, and F. Alet, Phys. Rev. B 106, 205121 (2022).
  * Yan  _et al._ [2021a] Z. Yan, Z. Zhou, O. F. Syljuåsen, J. Zhang, T. Yuan, J. Lou, and Y. Chen, Phys. Rev. B 103, 094421 (2021a).
  * Yan  _et al._ [2022] Z. Yan, Z. Y. Meng, D. A. Huse, and A. Chan, Phys. Rev. B 106, L041115 (2022).
  * Verresen and Vishwanath [2022] R. Verresen and A. Vishwanath, Phys. Rev. X 12, 041029 (2022).
  * Yan  _et al._ [2021b] Z. Yan, Y.-C. Wang, N. Ma, Y. Qi, and Z. Y. Meng, npj Quantum Mater. 6, 39 (2021b).
  * Pollmann  _et al._ [2011] F. Pollmann, J. J. Betouras, K. Shtengel, and P. Fulde, Phys. Rev. B 83, 155117 (2011).
  * Banerjee  _et al._ [2013] D. Banerjee, F.-J. Jiang, P. Widmer, and U.-J. Wiese, Journal of Statistical Mechanics: Theory and Experiment 2013, P12010 (2013).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Nature Physics 16, 132 (2020).
  * Glaetzle  _et al._ [2014] A. W. Glaetzle, M. Dalmonte, R. Nath, I. Rousochatzakis, R. Moessner, and P. Zoller, Phys. Rev. X 4, 041037 (2014).
  * Celi  _et al._ [2020] A. Celi, B. Vermersch, O. Viyuela, H. Pichler, M. D. Lukin, and P. Zoller, Phys. Rev. X 10, 021057 (2020).
  * Verresen  _et al._ [2021] R. Verresen, M. D. Lukin, and A. Vishwanath, Phys. Rev. X 11, 031005 (2021).
  * Yan  _et al._ [2023] Z. Yan, Y.-C. Wang, R. Samajdar, S. Sachdev, and Z. Y. Meng, Phys. Rev. Lett. 130, 206501 (2023).
  * Lieb and Wu [1972] E. H. Lieb and F. Y. Wu, in  _Phase transitions and critical phenomena: Exact Results_ , Vol. 1, edited by C. Domb and M. S. Green (Academic, London, 1972) Chap. 8, pp. 331–490.
  * Lieb [1967a] E. H. Lieb, Phys. Rev. 162, 162 (1967a).
  * Lieb [1967b] E. H. Lieb, Phys. Rev. Lett. 18, 1046 (1967b).
  * Lieb [1967c] E. H. Lieb, Phys. Rev. Lett. 19, 108 (1967c).
  * Sutherland [1967] B. Sutherland, Phys. Rev. Lett. 19, 103 (1967).
  * Lieb [1967d] E. H. Lieb, Phys. Rev. Lett. 18, 692 (1967d).
  * Pollock and Ceperley [1987] E. L. Pollock and D. M. Ceperley, Phys. Rev. B 36, 8343 (1987).
  * Henelius  _et al._ [1998] P. Henelius, S. M. Girvin, and A. W. Sandvik, Phys. Rev. B 57, 13382 (1998).
  * Leung  _et al._ [1996] P. W. Leung, K. C. Chiu, and K. J. Runge, Phys. Rev. B 54, 12938 (1996).
  * Paiva  _et al._ [2004] T. Paiva, R. R. dos Santos, R. T. Scalettar, and P. J. H. Denteneer, Phys. Rev. B 69, 184501 (2004).
  * Chen  _et al._ [2021] C. Chen, T. Yuan, Y. Qi, and Z. Y. Meng, Phys. Rev. B 103, 165131 (2021).
  * Costa  _et al._ [2018] N. C. Costa, T. Blommel, W.-T. Chiu, G. Batrouni, and R. T. Scalettar, Phys. Rev. Lett. 120, 187003 (2018).
  * Jiang  _et al._ [2022] W. Jiang, Y. Liu, A. Klein, Y. Wang, K. Sun, A. V. Chubukov, and Z. Y. Meng, Nature Communications 13, 2655 (2022).
  * Youngblood  _et al._ [1980] R. Youngblood, J. D. Axe, and B. M. McCoy, Phys. Rev. B 21, 5212 (1980).
  * Sutherland [1968] B. Sutherland, Physics Letters A 26, 532 (1968).
  * Falco [2013] P. Falco, Phys. Rev. E 88, 030103 (2013).
  * Fisher and Stephenson [1963] M. E. Fisher and J. Stephenson, Phys. Rev. 132, 1411 (1963).
  * Krauth and Moessner [2003] W. Krauth and R. Moessner, Phys. Rev. B 67, 064503 (2003).
  * Andrews  _et al._ [1984] G. E. Andrews, R. J. Baxter, and P. J. Forrester, J. Statist. Phys. 35, 193 (1984).
  * Pasquier [1987] V. Pasquier, Nuclear Physics B 285, 162 (1987).
  * Warnaar  _et al._ [1992] S. O. Warnaar, B. Nienhuis, and K. A. Seaton, Phys. Rev. Lett. 69, 710 (1992).
  * Blöte and Nightingale [1993] H. W. J. Blöte and M. P. Nightingale, Phys. Rev. B 47, 15046 (1993).
  * Blöte and Nienhuis [1994] H. W. J. Blöte and B. Nienhuis, Phys. Rev. Lett. 72, 1372 (1994).
  * Kondev and Henley [1995] J. Kondev and C. L. Henley, Phys. Rev. B 52, 6628 (1995).
  * Kondev and Henley [1994] J. Kondev and C. L. Henley, Phys. Rev. Lett. 73, 2786 (1994).
  * Wilkins and Powell [2020] N. Wilkins and S. Powell, Phys. Rev. B 102, 174431 (2020).
  * Sénéchal [2004] D. Sénéchal, in  _Theoretical Methods for Strongly Correlated Electrons_ (Springer, 2004) pp. 139–186.
  * Lukyanov and Terras [2003] S. Lukyanov and V. Terras, Nuclear Physics B 654, 323 (2003).
  * José  _et al._ [1977] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Phys. Rev. B 16, 1217 (1977).
  * Amit  _et al._ [1980] D. J. Amit, Y. Y. Goldschmidt, and S. Grinstein, Journal of Physics A: Mathematical and General 13, 585 (1980).
  * Note [1] The ice-point corresponds to the XXZ spin chain at $\Delta=1/2$ in notations where $\Delta>0$ corresponds to ferromagnetic interactions. The value of $g$ is given by $\Delta=\mathop{cos}\nolimits(\pi g)$, see e.g. Ref. [94, 104].
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Kosterlitz [1974] J. M. Kosterlitz, Journal of Physics C: Solid State Physics 7, 1046 (1974).
  * Ginsparg [1988] P. Ginsparg, Nuclear Physics B 295, 153 (1988).
  * Archambault  _et al._ [1998] P. Archambault, S. T. Bramwell, J.-Y. Fortin, P. C. W. Holdsworth, S. Peysson, and J.-F. Pinton, Journal of Applied Physics 83, 7234 (1998), https://doi.org/10.1063/1.367855 .
  * Atchison  _et al._ [2019] J. Atchison, A. Bhullar, B. Norman, and D. Venus, Phys. Rev. B 99, 125425 (2019).
  * Isakov and Moessner [2003] S. V. Isakov and R. Moessner, Phys. Rev. B 68, 104409 (2003).
  * Hikihara and Furusaki [1998] T. Hikihara and A. Furusaki, Physical Review B 58, R583 (1998).
  * Castelnovo  _et al._ [2006] C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, Phys. Rev. B 73, 144411 (2006).
  * Yan  _et al._ [2019] Z. Yan, Y. Wu, C. Liu, O. F. Syljuåsen, J. Lou, and Y. Chen, Phys. Rev. B 99, 165135 (2019).
  * Yan [2022] Z. Yan, Phys. Rev. B 105, 184432 (2022).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Nature 595, 233 (2021).

##  Appendix A Directed loop algorithm for loop models

Figure S1:  Directed loop Monte Carlo update steps for the fully packed loop
model. The local configuration to the left of the dotted line can transition
into one of the configurations to its right with probabilities determined by
detailed balance.

We estimate the thermal averages of observables in the classical loop model as
averages over Monte Carlo samples generated by a directed loop algorithm [34,
35, 36, 37, 38, 39] tailored for the loop model. The algorithm is summarized
below.

1\. Given a fully packed loop configuration $\mathcal{C}$, we pick with
uniform probability a site $S$; and then choose one of the two occupied edges
around $S$. With a Metropolis probability $p_{0}={\rm
min}(1,e^{\beta(E(\mathcal{C})-E(\mathcal{C}^{\prime}))})$, the dimer on this
edge is replaced by half a dimer (Fig. S1a). The new configuration
$\mathcal{C}^{\prime}$ has monomers on the site $S$ and at the end (monomer
$M$) of the dimer. The monomer $M$ has a binary valued ‘momentum’ internal
degree of freedom that is, initially, directed into the dimer and away from
$S$. In calculating the configurations $E$, it is assumed that the interaction
between half-dimers and parallel dimers is half that of full dimers. Note that
parallel dimers in a plaquette interact only if the other two edges are empty.
With probability $1-p_{0}$ the move is abandoned in this first step itself.

2\. If the monomer $M$ is moving into a dimer, annihilating it in the process
(Fig. S1b), it can, subsequently, create a dimer on one of the two previously
empty edges connected to the site ahead or the monomer can just reverse its
direction. Transition probabilities are chosen to satisfy detailed balance as
described further below.

3\. If the monomer is moving away from the dimer, growing a dimer in the
process (Fig. S1c), it can, subsequently, destroy one of the two dimers
connected to the node ahead or the monomer can just reverse its direction.
Transition probabilities are chosen to satisfy detailed balance as described
further below.

4\. We repeat steps (2) and (3) till the loop closes. If in the current
configuration $\mathcal{C}^{\prime}$, the monomer $M$ sees the starting site
$S$ ahead of it (Fig. S1d), the loop can terminate and produce a fully packed
configuration $\mathcal{C}$ with a Metropolis probability $p_{\rm term}={\rm
min}(1,e^{\beta(E(\mathcal{C}^{\prime})-E(\mathcal{C}))})$. With probability
$1-p_{\rm term}$, the $M$ reverses the direction instead.

The probabilities $p$ in steps (2) and (3) are chosen to satisfy detailed
balance. As shown in Fig. S1b,c the current configuration $\mathcal{C}$ can
transition into $\mathcal{C}^{\prime}_{0}$ with the monomer direction reversed
or two other configurations $\mathcal{C}^{\prime}_{1,2}$. The probability of
transition to $\mathcal{C}^{\prime}_{1,2}$ is given by

$\frac{W_{1,2}}{Z-\min(W_{1},W_{2},W_{0})}$ (S1)

where $W_{i}=e^{-\beta E(\mathcal{C}_{i})}$, $Z=\sum_{i}W_{i}$.

The monomer correlator $M(r)$ is computed as the histogram of displacements
between $S$ and $M$. In the description of the algorithm above we have used
the convention that the monomer sits at the center of the edge. We can instead
choose to place the monomer $M$ at some position $x\in[0,1]$ on the dimer and
associate an interaction energy between parallel full-dimers and the ‘partial’
dimers that is commensurate with $x$. We find that the choice of $x$ adds a
short range correction that does not affect the scaling properties.

Figure S2: The equal-time loop-segment correlation functions of (a)
longitudinal loop-segment correlations $C^{L}$, (b) transverse correlations
$C^{T}$, and (c) $\langle D_{0}D_{r}\rangle$ in the Eqs. (5), (6), and (8) in
the main text. The system size here is $L=256$. (d) The Coulomb gas constant
$g$, as extracted from the various fits presented in previous panels, as a
function of inverse temperature $\beta$.

##  Appendix B Raw data of the correlation functions

We present in Fig. S2 (a) and (b) the correlators $C^{L}(r)$ and $C^{T}(r)$
(measured in the ${\bf r}=(r,0)$ direction) as well as the fits to the
expressions Eq. (5) and Eq. (6). We furthermore present the correlator
$\langle D_{0}D_{r}\rangle$ associated to the order parameter in Fig. S2 (c),
which we fit to a single power-law as its leading contribution should decay as
$r^{-1/g}$. From these fits, we obtain estimates of $g$ represented in the
panel Fig. S2 (d), which are in overall agreement to those obtained from fits
to the adapted linear combinations of correlators presented in the main text
(see Fig. 6), albeit with slightly larger fluctuations at high temperature.

Figure S3:  Fits to obtain the coefficient $A$ in Eq. (14) for the
longitudinal loop-segment correlations $C^{L}$. Here we show the estimated
values of $A$ from $1/4g\pi^{2}$, where $g$ is obtained in three ways - (black
circular markers) from fitting to
$(-1)^{r}A^{\prime}/r^{\alpha_{S}}+B^{\prime}/r^{1/g}+C$, (red squares) from
$\langle W^{2}\rangle$ and (blue triangles $\Delta$) from $M(r)$ (last two as
described in Fig. 6). We also show $A^{\prime}$ obtained directly from the fit
to $(-1)^{r}A^{\prime}/r^{2}+B^{\prime}/r^{\alpha_{U}}+C$. The four curves are
consistent when $\beta$ is large.

Aside from the evaluation of the Coulomb gas constant, we also comment in this
Appendix on the evaluation of the amplitude [denoted $A$ in Eq. (5), (6) and
(7)] of the staggered part of loop segment correlators, which is also expected
to be universal (with $g$). We evaluate this amplitude from measurement of the
longitudinal correlator $C^{L}(r)$ in four different ways (see Fig. S3): (i)
assuming $A=1/(4\pi^{2}g)$ where $g$ is obtained from the scaling
$(-1)^{r}A^{\prime}/r^{\alpha_{S}}+B^{\prime}/r^{1/g}+C$. (ii) and (iii) are
$A=1/(4\pi^{2}g)$ where $g$ is from the winding number fluctuation $\langle
W^{2}\rangle$ and the monomer correlator $M(r)$ respectively. (iv)
$A^{\prime}$ obtained from the scaling form
$(-1)^{r}A^{\prime}/r^{2}+B^{\prime}/r^{\alpha_{U}}+C$.

We note an overall good agreement between all determinations of this
amplitude, as soon as $\beta\geq 0.2$, albeit with some small discrepancy at
$\beta=0$. Overall these data are consistent with the prediction
$A=1/(4\pi^{2}g)$.