# Backbone and shortest-path exponents of the two-dimensional $Q$-state Potts
model

Sheng Fang  MinJiang Collaborative Center for Theoretical Physics, College of
Physics and Electronic Information Engineering, Minjiang University, Fuzhou
350108, China  Department of Modern Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  Da Ke  MinJiang Collaborative
Center for Theoretical Physics, College of Physics and Electronic Information
Engineering, Minjiang University, Fuzhou 350108, China  Wei Zhong
<EMAIL_ADDRESS>MinJiang Collaborative Center for Theoretical Physics,
College of Physics and Electronic Information Engineering, Minjiang
University, Fuzhou 350108, China  Youjin Deng<EMAIL_ADDRESS>MinJiang
Collaborative Center for Theoretical Physics, College of Physics and
Electronic Information Engineering, Minjiang University, Fuzhou 350108, China
Department of Modern Physics, University of Science and Technology of China,
Hefei, Anhui 230026, China  Shanghai Research Center for Quantum Sciences,
Shanghai 201315, China

###### Abstract

We present a Monte Carlo study of the backbone and the shortest-path exponents
of the two-dimensional $Q$-state Potts model in the Fortuin-Kasteleyn bond
representation. We first use cluster algorithms to simulate the critical Potts
model on the square lattice and obtain the backbone exponents
$d_{\textsc{B}}=1.732\,0(3)$ and $1.794(2)$ for $Q=2,3$, respectively.
However, for large $Q$, the study suffers from serious critical slowing down
and slowly converging finite-size corrections. To overcome these difficulties,
we consider the O$(n)$ loop model on the honeycomb lattice in the densely
packed phase, which is regarded to correspond to the critical Potts model with
$Q=n^{2}$. With a highly efficient cluster algorithm, we determine from
domains enclosed by the loops
$d_{\textsc{B}}=1.643\,39(5),1.732\,27(8),1.793\,8(3),1.838\,4(5),1.875\,3(6)$
for $Q=1,2,3,2\\!+\\!\sqrt{3},4$, respectively, and $d_{\rm
min}=1.094\,5(2),1.067\,5(3),1.047\,5(3),1.032\,2(4)$ for
$Q=2,3,2+\sqrt{3},4$, respectively. Our estimates significantly improve over
the existing results for both $d_{\textsc{B}}$ and $d_{\rm min}$. Finally, by
studying finite-size corrections in backbone-related quantities, we conjecture
an exact formula as a function of $n$ for the leading correction exponent.

###### pacs:

05.50.+q (lattice theory and statistics), 05.70.Jk (critical point phenomena),
64.60.F- (equilibrium properties near critical points, critical exponents)

##  I Introduction

The Potts model [1, 2] is an extension of the celebrated Ising model and plays
an important role in the theory of phase transition and critical phenomena. It
also has broad applications in various fields like condensed-matter physics
[2]. Given a connected graph $G\equiv(V,E)$ with $V$ the vertex (site) set and
$E$ the edge set, each vertex has a spin of an integer value
$\sigma=0,1,\cdots,Q-1$, and two spins at the ends of an edge are coupled as
$-J\delta_{\sigma_{i},\sigma_{j}}$, where $J$ represents the interaction
strength and $\delta$ is the Kronecker delta function. Accordingly, the
partition function can be written as

${\mathcal{Z}}_{\rm Potts}=\sum_{\\{\sigma\\}}\prod_{\langle
ij\rangle}e^{J\delta_{\sigma_{i},\sigma_{j}}},$ (1)

where the inverse temperature $\beta$ has already been set to be $1$, the
summation $\\{\sigma\\}$ is over all possible spin configurations, and
$\langle ij\rangle$ is for all pairings of spins on the edges. For $Q=2$, the
Potts model reduces to the Ising model with the Ising coupling strength
$K=J/2$.

The Potts model can be reformulated in graphical representations, including
the Fortuin-Kasteleyn (FK) bond [3, 4] and $Q$-flow (loop) representations [5,
6]. Both the graphical representations can be obtained from the original spin
representation by high-temperature expansion techniques, and, recently, these
two graphical models also have been directly mapped onto each other with the
introduction of a loop-cluster joint model [7].

The FK bond representation is also known as the random-cluster (RC)
representation, in which each edge is either empty or occupied. Each occupied
bond has a statistical weight (relative to an empty one) as $v=\exp(J)-1$, and
each connected component (also called cluster) has a fugacity $Q$. The
partition function of the RC model then reads as

${\mathcal{Z}}_{\textsc{rc}}=\sum_{{\mathcal{A}}\subseteq
G}v^{{\mathcal{N}}_{b}}Q^{{\mathcal{N}}_{c}}\;,$ (2)

where the summation is over all possible subgraphs of $G$, and
${\mathcal{N}}_{b}$ and ${\mathcal{N}}_{c}$ represent the total numbers of
occupied bonds and clusters, respectively. As an important consequence, the
parameter $Q$ can now take any nonnegative real number. For ferromagnetic
coupling $J\\!>\\!0$, the bond weight $v$ is positive, and each configuration
in Eq. (2) has a probability interpretation. As special cases, the RC model
reduces to the standard bond percolation in the $Q\\!\rightarrow\\!1$ limit,
and the spanning tree or forest for $Q\\!\rightarrow\\!0$ [4]. The FK
representation has an important role in conformal field theory [8] and in
stochastic Loewner evolution [9, 10], leading to much advanced theoretical
progresses for the Potts model. Further, based on passing back and forth
between the spin and FK representations, efficient cluster methods (i.e., the
Swendsen–Wang (SW) [11] and Wolff algorithms [12]) are developed and widely
used.

In addition to physical quantities from the original spin representation, the
RC model has very rich geometric structures associated with FK random
clusters, and a variety of critical exponents are used to characterize these
geometric behaviors, which were first introduced in percolation [13, 14, 15].
The Euclidean diameter $\xi_{1}$ of the largest cluster acts as the
correlation length and diverges as $\xi_{1}\\!\sim\\!t^{-\nu}$, with
$t\\!\equiv\\!(v_{c}-v)/v_{c}$ and $v_{c}$ the critical point, and $\nu$ is
frequently called the correlation-length exponent. The mass $s$ of any
critical cluster has a power-law dependence on its diameter $\xi_{s}$ as
$s\sim\xi_{s}^{d_{\textsc{f}}}$ for $s\gg 1$, with $d_{\textsc{f}}$ the
fractal dimension. In two dimensions (2D), the hull and the external perimeter
of clusters can be further defined and have fractal dimensions $d_{\rm hull}$
and $d_{\textsc{ep}}$, respectively. For a pair of connected sites with
distance $|\mathbf{x}|\gg 1$, the graph distance $\mathcal{S}$, i.e., the
minimum length of all connecting paths between them, algebraically diverges as
$\mathcal{S}\sim|\mathbf{x}|^{d_{\rm min}}$, with $d_{\rm min}$ the shortest-
path exponent. If a voltage difference is applied between these two sites, the
total number $N_{\textsc{b}}$ of bonds carrying nonzero current scales as
$N_{\textsc{b}}\sim|\mathbf{x}|^{d_{\textsc{B}}}$, where $d_{\textsc{B}}$ is
called the backbone exponent. Further, the total number $N_{\rm red}$ of red
bonds, which carry all the current, behaves as $N_{\rm
red}\sim|\mathbf{x}|^{d_{\rm red}}$, with $d_{\rm red}$ the red-bond exponent.

In the past decades, the $Q$-state Potts model has been extensively studied
[16, 17, 18, 19]. In 2D, the Potts model exhibits a second-order phase
transition for $Q\\!\leq\\!4$, and the exact values of most critical exponents
have been identified. In the Coulomb gas theory [18, 19], the Coulomb-gas
coupling strength $g$ relates to $Q$ as

$Q=4\cos^{2}(\pi g/4)\;,\hskip 14.22636ptg\in[2,4],$ (3)

and the thermal and magnetic exponents of the leading and subleading
renormalization fields are known as

$\displaystyle y_{\rm t1}=3-6/g\;,\qquad\qquad\;\>y_{\rm t2}=4-16/g\;,$ (4)
$\displaystyle y_{\rm h1}=\frac{(g+2)(g+6)}{8g}\,,\;\;\;\;\;y_{\rm
h2}=\frac{(g+10)(g-2)}{8g}\;,$

where the leading thermal exponent $y_{\rm t1}$ is the inverse of the
correlation-length exponent $\nu$ as $y_{\rm t1}=1/\nu$. The subleading field,
with $y_{\rm t2}\\!<\\!0$, governs the convergence of leading corrections;
logarithmic corrections arise as $Q$ approaches to 4, for which $y_{\rm
t2}\\!\rightarrow\\!0$.

The fractal dimension $d_{\textsc{f}}$ of FK clusters is identical to the
leading magnetic exponent as $d_{\textsc{f}}=y_{\rm h1}$, and the exact values
of some other geometric exponents are [20]

$\displaystyle d_{\textsc{\rm hull}}=1+2/g\;,\hskip
56.9055ptd_{\textsc{ep}}=1+g/8\;,$ $\displaystyle d_{\rm
red}\>=(4-g)(3g+4)/g\;.$ (5)

However, there is still a set of geometric critical exponents, including
$d_{\textsc{B}}$ and $d_{\rm min}$, whose exact values are unavailable. Monte
Carlo simulations have been used to estimate them [21, 22, 23, 24, 25, 26],
and Table 1 lists some estimates of $d_{\textsc{B}}$ and $d_{\rm min}$.

In this work, we apply the SW and Wolff cluster methods to simulate the
critical Potts model on the square lattice, expecting to obtain better results
for $d_{\textsc{B}}$ and $d_{\rm min}$ than those in Refs. [21, 23, 24, 26,
22]. An improved algorithm is formulated to classify the occupied bonds into
“bridges” and “nonbridges” [23]. From the fractal dimension of bridge-free
clusters, we determine $d_{\textsc{B}}=1.732\,0(3)$ for $Q=2$ and $1.794(2)$
for $Q=3$. The estimate of $d_{\textsc{B}}$ for $Q=2$ is consistent with Ref.
[24], and for $Q=3$ it nearly rules out the result in Ref. [21]. As $Q$
increases, however, the cluster methods suffer from severe critical slowing-
down and finite-size corrections become very strong, preventing us to obtain
high precision of $d_{\textsc{B}}$ for large $Q$.

($Q$,$g$) |  $d_{\textsc{B}}$ |  $d_{\rm min}$
---|---|---
(1,8/3) |  1.643 36(10)[23] |  1.130 77(2)[26]
present |  1.643 39(5) |  -
(2,3) |  1.732 1(4)[24] |  1.094 0(2) [24]
present |  1.732 27(8) |  1.094 5(2)
(3,10/3) |  1.789 5(5)[27] |  1.066 2(30) [22]
present |  1.793 8(3) |  1.067 5(3)
(2+$\sqrt{3}$,11/3) |  - |  -
present |  1.838 4(5) |  1.047 5(3)
(4,4) |  1.873(4)[27] |  -
present |  1.875 3(6) |  1.032 2(4)
Table 1: Best estimates of $d_{\textsc{B}}$ and $d_{\rm min}$ for the two-
dimensional Potts model with $Q=1,2,3,2+\sqrt{3},4$ and $g$ is the Coulomb-gas
coupling strength. These final quoted values and the error margins are
estimated from a variety of fitting results by taking into account statistical
and systematic errors.

We then consider the O$(n)$ loop model [28] on the honeycomb lattice, of which
the configuration is a gas of nonintersecting loops. Let $n$ be the fugacity
of a loop and $x$ be the statistical weight for a loop segment (an occupied
bond); the partition function reads as

$\displaystyle{\mathcal{Z}}_{\rm loop}=\sum_{\rm
loops}x^{{\mathcal{L}}}n^{{\mathcal{N}}_{\ell}}\;,$ (6)

where ${\mathcal{N}}_{\ell}$ is the number of loops and ${\mathcal{L}}$ is the
total length of all the loops. For $n=1$, the loop model reduces to the flow
(loop) representation of the honeycomb-lattice Ising model, with the coupling
strength as $\tanh(K)=x$ [29].

The phase diagram of the O($n$) loop model is shown in Fig. 1. As the bond
weight $x$ is increased, the system undergoes a second-order phase transition
$x_{+}$ from the dilute phase, consisting of small loops, into a critical
densely-packed (DP) phase, where loops have fractal structures. The
universality of the DP phase is governed by the line $x_{-}$ of stable fixed
points. The exact locations of $x_{\pm}(n)$ are known as [30]

$1/x_{\pm}=\sqrt{2\pm\sqrt{2-n}}\;.$ (7)

The two curves meet at $(n\\!=\\!2,x_{\pm}\\!=\\!1/\sqrt{2})$. In the language
of universality, the O($n$) loop model on the critical branch $x_{+}$ (DP
phase) corresponds to the tricritical (critical) $Q$-state Potts model with
$Q=n^{2}$ [30, 31]. An interesting note is that, at $x_{+}$, the O(1) loop
model is typically regarded as the critical Ising model instead of the
tricritical one-state Potts model.

Figure 1: Phase diagram of the O($n$) loop model on the honeycomb lattice. The
system is in dilute phase if bond weight $x\\!<\\!x_{+}$ and in densely packed
(DP) phase if $x\\!>\\!x_{+}$. Vertical blue arrows sketch the directions of
renormalization flows, and illustrate that the whole DP phase
($x_{+}\\!<\\!x\\!<\\!\infty$) is critical and in the same universality for
each $n$.

In other words, despite of the absence of an exact transformation between the
loop model (6) and the Potts model (1), the two systems have the same Coulomb-
gas coupling strength $g$ as long as $Q=n^{2}$. This correspondence is
expected to hold true for any real $0\leq Q\leq 4$, where, for noninteger $Q$,
the Potts model should be formulated in the FK bond representation. The
relation between $g$ and $Q$ is still given by Eq. (3), with $g\in[2,4]$ for
$x_{-}$ and $g\in[4,6]$ for $x_{+}$. Furthermore, Eqs. (4) and (I) for
critical exponents still hold true. However, what physical quantities in the
loop model, described by these equations, need to be identified, and this has
been largely explored in Refs. [32, 33]. The thermal field along the $x$
direction has exponent $y_{\rm t2}$ in Eq. (4), instead of $y_{\rm t1}$. The
loops and the domains (faces) enclosed by the loops have the fractal
dimensions $d_{\rm hull}$ and $d_{\textsc{f}}=y_{\rm h1}$, respectively. The
quantities, of which the scaling is governed by $y_{\rm t1}$, are identified
as the magnetization and the susceptibility of Ising spins living on the
triangular lattice that is dual to the honeycomb lattice.

Highly efficient worm-type and cluster algorithms have been developed for the
O($n$) loop model with $n\geq 1$ [32, 33]. It is shown that critical slowing-
down barely exists at $x_{+}$, and is completely absent in the DP phase, as
qualitatively understood as following. For worm-type and cluster algorithms,
the dynamic exponent $z$ has the so-called Li-Sokal bound [34] $z\geq z_{\rm
min}=\alpha/\nu$, where $\alpha/\nu$ is the scaling exponent for specific
heat, and, further, this lower bound is rather sharp $z\approx z_{\rm min}$ in
2D. For the critical Potts model, one has $z_{\rm min}\\!=\\!2y_{\rm t1}-2$,
leading to $z\\!>\\!0$ for $Q\\!>\\!2$. Moreover, as $Q$ increases, the
critical slowing-down becomes more severe, since $y_{\rm t1}$ is an increasing
function of $Q$. However, the O($n$) loop model has $\alpha/\nu=2y_{\rm
t2}-2$, and thus, $z_{\rm min}(n\geq 1)\leq 0$ along $x_{+}$ and $z_{\rm
min}\\!<\\!-2$ in the whole DP phase.

We aim to determine the backbone exponent $d_{\textsc{B}}$ and the shortest-
path exponent $d_{\rm min}$ by studying the O($n$) loop model along the line
$x_{-}(n)$. For this goal, we assume that the domains of DP loops and the
corresponding FK random clusters not only have the same fractal dimension, but
also exhibit the same scaling for other geometric properties. In comparison
with the critical Potts model, the advantage of studying the loop model along
$x_{-}(n)$ is twofold: The absence of critical slowing-down and the absence of
finite-size corrections from $y_{\rm t2}$ (corrections from other sources can
still exist).

Our final estimates of $d_{\textsc{B}}$ and $d_{\rm min}$ are given in Table
1, and significantly improve over the existing results. The agreement with the
FK bond representation confirms our assumption that the domains of the DP
O($n$) loops have the same geometric structures as the critical FK random
clusters, and we expect that this correspondence can be extended to the
$x_{+}$ branch.

In the analysis of quantities associated with bridge-free (backbone) clusters
for the DP O($n$) loop model, we observe that, despite the absence of
corrections from $y_{\rm t2}$, finite-size corrections are still significant
and become more and more severe as $n$ increases. For the special case of
$n=2$, logarithmic corrections seem to arise. On the basis of our numerical
results, we conjecture an exact formula, as a function of $n$, for the leading
correction exponent.

The remainder of this paper is organized as follows. Algorithms for simulation
and measurement are described in Sec. II, together with a list of sampled
quantities. Section III presents our results for the backbone exponent and
Sec. IV presents results for the shortest-path exponent. A brief discussion is
given in Sec. V.

##  II Algorithms and Observables

###  II.1 Simulation of the Potts model

For integer $Q$, the celebrated SW and Wolff cluster algorithms are used, and,
for real $Q\geq 1$, the Chayes-Machta (CM) method [35] is applied. They can be
understood by the so-called induced-subgraph picture as following Ref. [32].

For the RC model (2) on graph $G=(V,E)$, one first decomposes $Q$ as
$Q=Q_{\alpha}\\!+\\!Q_{\beta}$. Independently for each cluster in a FK bond
configuration, one then chooses color “$\alpha$” with probability
$Q_{\alpha}/Q$ or color “$\beta$” with probability $Q_{\beta}/Q$, and assigns
the chosen color to all sites in the cluster. Consequently, the lattice sites
are partitioned as $V\\!=\\!V_{\alpha}\cup V_{\beta}$, and the edge set $E$
becomes $E=E_{\alpha}\cup E_{\beta}\cup E_{\alpha\beta}$; an edge $e\in
E_{\alpha\beta}$ connects a pair of sites, respectively, in $V_{\alpha}$ and
$V_{\beta}$, and it is not occupied by definition. It can be seen that,
conditioning on this decomposition, the bond configuration is nothing other
than the combination of a $Q_{\alpha}$-state RC model on the induced subgraph
$G_{\alpha}=(V_{\alpha},E_{\alpha})$ and another $Q_{\beta}$-state RC model on
$G_{\beta}$. Now, one can update these induced RC models by any valid
algorithm. One valid update is the identity operation, which is “do nothing”
for the RC model, corresponding to the “inactive” color of Chayes and Machta
[35]. Of course, we must also include at least one nontrivial update.

A simple choice is to set $Q_{\alpha}=1$ for $Q\geq 1$, so that the
corresponding induced model on $G_{\alpha}$ is the standard bond percolation
and can be trivially updated, i.e., each edge $e\in E_{\alpha}$ is
independently occupied with probability $p$. On this basis, a basic variant of
CM algorithm is formulated as:

  1. 1.

Independently for each cluster, choose the active color “$\alpha$” with
probability $1/Q$ and otherwise “$\beta$” (inactive), resulting in a random
partition $V=V_{\alpha}\cup V_{\beta}$ as well as $E=E_{\alpha}\cup
E_{\beta}\cup E_{\alpha\beta}$.

  2. 2.

Independently for each edge $e\in E_{\alpha}$, place an occupied bond with
probability $p=v/(1+v)$ according to Eq. (2), and for $e\in
E\\!\setminus\\!E_{\alpha}$, do nothing.

Here $E\setminus E_{\alpha}$ represents a complementary subset of $E_{\alpha}$
with edges in $E_{\alpha}$ being excluded. A new FK bond configuration is then
obtained after ignoring/discarding the randomly assigned colors.

For $Q\geq 2$, the efficiency of the CM algorithm can be further improved by
increasing the number of active colors to be $m=[Q]$, with $[Q]$ the ground
integer of $Q$. Namely, for each cluster, an active color
$\alpha_{i}\;(i=1,\cdots,m)$ is chosen with probability $1/Q$, leading to $m$
copies of the bond percolation on $G_{\alpha_{i}}$. For integer $Q$, where the
“do-nothing” operation is absent, the celebrated SW algorithm is recovered. A
single-cluster variant of the CM method is also available [36].

###  II.2 Simulation of the O($n$) model

The induced subgraph picture can also provide a versatile platform to design
Monte Carlo algorithms for the O($n$) model (6) with $n\geq 1$. The basic idea
is to set $n=n_{\alpha}+n_{\beta}$ with $n_{\alpha}=1$, and to simulate the
induced O(1) loop model. More precisely, independently for each loop, one
randomly assigns the active color “$\alpha$” with probability $1/n$ and the
inactive color “$\beta$” with probability $1-1/n$. In addition, one assigns
the active color to all those sites that are not on loops, since they have an
implicit weight of “1” (unity), as in Eq. (6). The O(1) loop model on the
induced subgraph can then be simulated by the worm algorithm. A variant of
worm-type algorithm for the O$(n)$ loop model in 2D and 3D can be found in
Refs. [33, 37].

Note that the honeycomb-lattice loops are the boundaries of Ising domains on
the dual triangular lattice, with the one-to-two correspondence between the
loop and the spin configuration of dual Ising spins. The dual coupling
strength $K^{*}$ relates to $x$ as $2K^{*}\\!=\\!-\ln(x)$, which is
ferromagnetic for $x\\!<\\!1$ and antiferromagnetic for $x\\!>\\!1$. The
number of dual Ising domains is simply
${\mathcal{N}}_{d}={\mathcal{N}}_{\ell}+1$, with ${\mathcal{N}}_{\ell}$ the
loop number and, thus, $n$ is also the fugacity of an Ising domain. In other
words, the O$(n)$ loop model can be regarded as a generalized Ising model on
the triangular lattice, of which the partition function is written as

${\mathcal{Z}}_{\rm GIsn}=\sum_{\\{s\\}}n^{{\mathcal{N}}_{d}}\prod_{\langle
ij\rangle}\exp(K^{*}s_{i}s_{j})\;,$ (8)

where the summation is over all the Ising configurations. Figure 2
demonstrates the relation between the loop model and generalized Ising model.
The Ising spinslive on the sites of the triangular lattice,and the spin values
are specified by the colors of the hexagonal faces. The domain walls are drawn
between the boundaries of domains, which corresponds to the loops in the
O$(n)$ loop model on the honeycomb lattice. In Fig. 2, there are
${\mathcal{N}}_{d}=11$ domains and ${\mathcal{N}}_{\ell}=10$ loops.

Figure 2: Illustration of the generalized Ising model on the triangular
lattice and the O($n$) loop model on the honeycomb lattice with periodic
boundary conditions. The hexagons in red (blue) denote the spin-up (spin-down)
sites and loops, i.e., the domain walls between the generalized Ising spins,
are specified by blue lines.

Applying the induced-subgraph picture to Eq. (8), one obtains a standard Ising
model on induced subgraph $G_{\alpha}$ of the triangular lattice, which can be
simulated by the SW or Wolff cluster algorithms. The bond-occupation
probability is $p=1\\!-\\!e^{-2K^{*}}\\!=\\!1-x$ between a neighboring pair of
parallel spins for $x\\!\leq\\!1$ or $p=1\\!-\\!e^{2K^{*}}\\!=\\!1\\!-\\!1/x$
between anti-parallel spins for $x\\!>\\!1$. Nevertheless, some cautions are
needed for the identify operation to the complementary subgraph
$G\\!\setminus\\!G_{\alpha}$, which can no longer be “do-nothing”. Instead,
one should place occupied bonds on all the edges $e\in
E\\!\setminus\\!E_{\alpha}$, independent of the underlying Ising spins, so
that all sites in $V_{\beta}$ and their neighboring sites in $V_{\alpha}$ are
in the same cluster. As a result, the domain topology on
$G\\!\setminus\\!G_{\alpha}$ is kept unchanged, thanks to the Ising symmetry.
For explicitness, the bond occupation probability for $x\\!\leq\\!1$ is listed
as

$p=\left\\{\begin{array}[]{lcl}1-x&&\text{if $s_{i}=s_{j}$ and $i\in
V_{\alpha},j\in V_{\alpha}$;}\\\ 0&&\text{if $s_{i}\neq s_{j}$ and $i\in
V_{\alpha},j\in V_{\alpha}$;}\\\ 1&&\text{otherwise.}\\\ \end{array}\right.$
(9)

For $x\\!>\\!1$, similar procedure can be applied and the only difference is
as follows: Within the active subgraph $G_{\alpha}$, the bond is placed
between antiparallel spins ($s_{i}\neq s_{j}$) with probability $p=1-1/x$.

For the O$(n)$ loop model in the DP phase, both the cluster algorithm and the
worm method exhibit no critical slowing-down, and their efficiency is
approximately the same. Since we are interested in the geometric structures of
Ising domains in this work, simulating the generalized Ising model (8) can
avoid the complication of passing back and forth between the honeycomb and the
triangular lattice.

###  II.3 Identification of nonbridges

Following Ref. [23], we classify the occupied bonds of a FK bond configuration
into bridges and nonbridges, and study the resulting bridge-free clusters to
determine the backbone dimension $d_{\textsc{B}}$. A bridge is an occupied
bond whose deletion would break a cluster into two, and removing all the
bridges would produce a bridge-free configuration, which consists of isolated
sites and blobs. In other terminology, blobs are also called bridge-free
clusters and biconnected clusters, in which any two sites are connected by at
least two independent paths. For the Ising domains on the triangular lattice,
occupied bonds have to be first placed within the domains before the bond
classification.

Two distinct methods to identify nonbridges are introduced in Refs. [23] and
[38], of which the former makes use of the planarity of 2D lattices. In the
latter, clusters are grown by breadth-first procedure with help of a treelike
data structure; whenever a newly visited bond is to close a loop, one
backtracks the two branches of tree until the joint site, and all passing
bonds are marked as nonbridges. A problem is that a nonbridge may be visited
for many times. Since the number of nonbridges is of system size order, which
is observed in Refs. [23, 24], a higher efficient algorithm is desired to
identify nonbridges.

Figure 3: Identification procedure for nonbridges. (a), Grow the depth-first
oriented walk untill it is a to-be-closed loop, as labeled by the walk
$[0],[1],[2],\cdots,[9]$. (b), Backtrack to the root site ($[2]$), construct a
flat tree (red arrows), and mark all the passing bonds as nonbridges (thick
blue bonds). (c), Apply the last-in-first-out rule so that the old elements
$[9],\cdots,[6]$ are removed and and new elements, $[6]$ and $[7]$, are added
to the walk. (d), Add a new loop and replace the existing root by the new one.
The final configuration consists of four bridge-free clusters, three of which
are isolated sites (those with red circles).

Hereby, we employ a depth-first growth with a chainlike data structure
${\mathcal{W}}[\ell]$ and introduce a treelike graph ${\mathcal{D}}[i]$ with
some minimal depth. The chainlike array ${\mathcal{W}}[\ell]$ is to keep trace
of the depth-first growth: Let ${\mathcal{W}}[\ell]$ denote the newly added
site, ${\mathcal{W}}[\ell\\!+\\!1]$ stores the next to-be-added site. It is
required that, at one step, only one of the unvisited sites connecting to
${\mathcal{W}}[\ell]$ is added, and thus the chain, ${\mathcal{W}}[0]$,
${\mathcal{W}}[1],\cdots,{\mathcal{W}}[\ell\\!+\\!1]$, is nothing other than
an oriented walk. Further, the “last-in-first-out” rule is applied: If all
neighboring sites of ${\mathcal{W}}[\ell]$ have been visited, then
${\mathcal{W}}[\ell]$ is removed and the incoming site
${\mathcal{W}}[\ell\\!-\\!1]$ becomes the walk frontier. Array
${\mathcal{D}}[i]$ is a directed graph consisting of many trees, and each tree
is to represent a bridge-free cluster (including an isolated site). The parent
site of site $i$ is given by ${\mathcal{D}}[i]$, and the root of a tree is the
site with the parent being itself. Initially, ${\mathcal{D}}[i]=i$ is set for
all sites $i\in V$.

The procedure is sketched in Fig. 3. Starting with an initial site, the depth-
first walk is performed until it is a to-be-closed loop, labeled by
$[0],[1],\cdots,[9]$ in Fig. 3(a). Then, the root is found starting from the
joint site, [2], and the loop is identified by backtracking ${\mathcal{D}}[i]$
(if applicable) or ${\mathcal{W}}[\ell]$ to the root. Meanwhile, all the
passing bonds are marked as nonbridges (thick blue bonds) and all the passing
sites take the root as their parent [Fig. 3(b)]. Afterwards, the depth-first
walk is continued with the last-in-first-out rule, as illustrated by Fig.
3(c). When a new loop is added, the existing root is taken place by the new
one if they are different [Fig. 3(d)]. In this way, it is kept that the parent
${\mathcal{D}}[i]$ of site $i$ is always visited no later than site $i$
itself, so that the loop identification by backtracking ${\mathcal{D}}[i]$
works properly. When all bonds in the cluster are visited, the depth-first
walk is done, and all the nonbridges are successfully identified. As a by-
product, all bridge-free clusters are constructed so that if one can further
keep track of the sizes of trees, the associated quantities can be calculated
immediately.

In the backtracking procedure for identifying a loop, ${\mathcal{D}}(i)$ has a
higher priority than ${\mathcal{W}}[\ell]$: If ${\mathcal{D}}(i)\neq i$, then
the backtracking is applied to ${\mathcal{D}}[i]$ instead of
${\mathcal{W}}[\ell\\!-\\!1]$. This is to avoid visiting a nonbridge for many
times. In this work, we observe that trees ${\mathcal{D}}(i)$ are very flat
and their typical depth is of a few units. In practice, the trees can be
further flatten: When a loop is to be formed, one tries to continue the walk
by visiting another neighboring site of ${\mathcal{W}}[\ell]$. In other words,
a self-avoiding walk is performed as long as possible before the backtracking
procedure. Therefore, the present method has a nearly optimal efficiency for
identifying nonbridges.

Sampled quantities.–On the basis of bridge-free clusters, we sample the size
of the largest cluster as $B_{1}=\langle{\mathcal{B}}_{\textsc{1}}\rangle$ and
calculate the second moment of cluster sizes as
$B_{2}=L^{-2}\langle\sum{\mathcal{B}}^{2}\rangle$, where the sum is over all
bridge-free clusters ${\mathcal{B}}$ and $\langle\cdot\rangle$ represents the
ensemble average. For the O$(n)$ loop model, we also distinguish the spin-up
domains from the spin-down ones, and calculate the averaged size
$B_{\textsc{a}}=\langle({\mathcal{B}}_{1\uparrow}+{\mathcal{B}}_{1\downarrow})/2\rangle$.

To determine the shortest-path exponent $d_{\rm min}$, we use the breadth-
first algorithm to grow clusters; the number ${\mathcal{S}}$ of layers after
finishing the cluster growth corresponds to the maximum graph distance between
the seed site and any other site in the cluster. We then sample the longest
graph distance among all the clusters as
$S_{1}=\langle{\mathcal{S}}_{1}\rangle$.

At criticality, these sampled quantities are expected to asymptotically
diverge as a power-law scaling of linear system size $L$. More precisely, they
behave as

$\displaystyle S_{1}$ $\displaystyle\sim$ $\displaystyle L^{d_{\rm
min}}\;,\hskip 42.67912pt\mbox{and}$ $\displaystyle B_{1}$ $\displaystyle\sim$
$\displaystyle B_{\textsc{a}}\sim L^{d_{\textsc{B}}}\;,\hskip
22.76219ptB_{2}\sim L^{2d_{\textsc{B}}-d}\;.$ (10)

##  III Results for backbone exponent

We simulate the critical $Q$-state Potts model on the square lattice and the
DP O$(n)$ loop model with system size $L$ ranging from 4 to 8192. Periodic
boundary conditions are implemented. For the Potts model, we take $Q=2,3$ and
the exact critical point $v_{c}=\sqrt{Q}$. For the loop model, the simulation
is at the branch $x_{-}(n)$ with $n=1,\sqrt{2},\sqrt{3},\sqrt{2+\sqrt{3}},2$,
corresponding to $Q=1,2,3,2+\sqrt{3},4$, respectively. We generate more than
$10^{8}$ independent samples for each system size $L<256$, and no fewer than
$10^{7}$ independent samples for $L\geq 256$.

We perform least-squares fits to the finite-size scaling ansatz

${\mathcal{O}}=L^{y_{{\mathcal{O}}}}(a+b_{1}L^{y_{1}}+b_{2}L^{y_{2}})+c_{0}\;,$
(11)

where ${\mathcal{O}}=S_{1},B_{1},B_{\textsc{a}},B_{2}$ and the corresponding
exponent $y_{{\mathcal{O}}}$ is given by Eq. (10). The term with $c_{0}$
describes the background contribution of the system, and those with $b_{1}$
and $b_{2}$ account for leading and subleading finite-size corrections,
respectively. In comparison with the leading scaling term with amplitude $a$,
all the other three terms play a role of corrections. As a precaution against
correction-to-scaling terms that we miss including in the fitting ansatz, we
impose a lower cutoff $L\geq L_{\rm m}$ on the data points admitted in the
fits, and systematically study the effect on the residuals (denoted by
$\chi^{2}$) by increasing $L_{\rm m}$. In general, the preferred fit for any
given ansatz corresponds to the smallest $L_{\rm m}$ for which the goodness of
the fit is reasonable and for which subsequent increases in $L_{\rm m}$ do not
cause the $\chi^{2}$ value to drop by vastly more than one unit per degree of
freedom (DF). In practice, by “reasonable” we mean that $\chi^{2}/{\rm
DF}\approx 1$. The systematic error is obtained by comparing estimates from
various reasonable fitting ansatz.

$Q$ | $L_{\rm m}$ | $d_{\textsc{B}}$ | $a$ | $b_{1}$ | $y_{1}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---|---
2 | 8 | 1.732 17(9) | 0.7995(6) | 0.1168(9) | -0.66(1) | 4.0/7
| 16 | 1.732 0(1) | 0.8003(8) | 0.121(4) | -0.68(2) | 2.1/6
| 32 | 1.732 0(2) | 0.800(1) | 0.12(1) | -0.68(5) | 2.1/5
| 16 | 1.732 23(8) | 0.7990(5) | 0.104(3) | -5/8 | 2.3/6
| 32 | 1.732 1(1) | 0.7995(8) | 0.098(7) | -5/8 | 1.5/5
3 | 16 | 1.794 0(6) | 0.638(7) | 0.415(2) | -0.44(1) | 10.0/6
| 32 | 1.793 6(9) | 0.64(1) | 0.420(9) | -0.46(3) | 9.8/5
| 16 | 1.793 3(1) | 0.647(1) | 0.415(2) | -0.46 | 12.0/7
| 32 | 1.793 5(2) | 0.645(2) | 0.421(4) | -0.46 | 9.8/6
| 64 | 1.793 2(3) | 0.648(3) | 0.411(10) | -0.46 | 8.6/5
Table 2: Fitting results of $B_{1}$ for the Ising model and of $B_{2}$ for the
$Q=3$ Potts model.

###  III.1 Estimate from the Potts model

The $Q\\!\rightarrow\\!1$ Potts model is the standard bond percolation, and
the geometric structures of critical percolation clusters on the square
lattice have been extensively studied in Ref. [23], which gives
$d_{\textsc{B}}=1.643\,36(10)$.

Figure 4: Estimated backbone exponent $d_{\textsc{B}}$ for the Ising model (a)
and the $Q=3$ Potts model (b). The approximate linearity of red lines
indicates that the leading correction exponent is $y_{1}\approx-0.6$ for $Q=2$
and $-0.4$ for $Q=3$. The upward and downward bending of other curves reflects
the reliability of the finally quoted central values and their error bars.
Figure 5: Backbone exponent $d_{\textsc{B}}$ for the DP O($n$) loop model with
$n=\sqrt{2}$ (a) and $n=\sqrt{3}$ (b). The upward or downward bending of
curves with the $d_{\textsc{B}}$ value deviating from the estimated value
illustrates the reliability of the quoted error bars. By comparing to Fig. 4,
it is seen that the results from the loop model have higher precision.

The Coulomb-gas coupling of the Ising model ($Q=2$) is $g=3$, and, according
to Eq. (4), the subleading thermal and magnetic scaling fields have exponents
$y_{\rm t2}=-4/3$ and $y_{\rm h2}=13/24$, respectively. However, the
subleading magnetic field with $y_{\rm h2}$ is considered to be redundant due
the $Z_{2}$ symmetry of the Ising spins, and the subleading thermal field with
$y_{\rm t2}$ is believed to play no role in the scaling of thermodynamic
quantities in the spin representation. Indeed, finite-size corrections with
exponent $y_{\rm t2}=-4/3$ are not observed in physical observables; for
quantities like energy, magnetization, specific heat, and susceptibility, the
leading finite-size corrections are found to be governed by exponent $-2$.
Nevertheless, it remains open whether corrections with exponent $-4/3$ would
arise in some geometric quantities.

According to the least-squares criterion, the $B_{1}$ data for $Q=2$ are
fitted to Eq. (11). At first, we set $b_{2}=0$ and leave $d_{\textsc{B}}$,
$a$, $b_{1}$, $y_{1}$ and $c_{0}$ as free parameters, but no stable fitting
results can be obtained. By further fixing $c_{0}=0$, which is effectively a
correction term, we obtain $d_{\textsc{B}}=1.732\,0(2)$ and $y_{1}=-0.68(5)$,
of which the correction exponent is much larger than $-2$ or $-4/3$. The
results are shown in Table 2. Therefore, for geometric quantities associated
with bridge-free clusters, the leading finite-size corrections probably arise
from some other resource. In the framework of conformal field theory, an
exponent $-5/8$ appears at several places in the Kac table for the Ising model
[39]. With $y_{1}=-5/8$ being fixed, the fit yields
$d_{\textsc{B}}=1.732\,23(8)$ for $L_{\rm m}=16$ (Table 2). Following the same
procedure, the fit of the $B_{2}$ data gives $d_{\textsc{B}}=1.732\,3(5)$,
which is consistent with $d_{\textsc{B}}$ from $B_{1}$ but has a slightly
larger error bar. Comparing the fitting results of $B_{1}$ and $B_{2}$ for
different $L_{\rm m}$, we take the final estimate as
$d_{\textsc{B}}=1.732\,0(3)$. We choose not to seriously take into account the
smaller error bar from the fits with $y_{1}=-5/8$ being fixed, since it is
only a crude guess and the exact value of $y_{1}$ is not available. In Fig.
4(a), the $B_{1}/L^{d_{\textsc{B}}}$ data are plotted versus $L^{-5/8}$, with
$d_{\textsc{B}}$ chosen to be the central value of the estimate as well as the
central value plus or minus three error bars. The approximate linearity of the
red line and the upward and downward bending curves reflect the reliability of
the final estimate $d_{\textsc{B}}=1.732\,0$ and of the quoted error margin
$0.000\,3$.

The leading correction exponent for the $Q=3$ Potts model, with the Coulomb-
gas coupling $g=10/3$, is $y_{\rm t2}=-4/5$. By fixing $b_{2}=0$ for the
subleading correction term and $c_{0}=0$ for the background contribution, we
fit the $B_{2}$ data to Eq. (11) and obtain $d_{\textsc{B}}=1.793\,6(10)$ and
$y_{1}=-0.46(3)$. The value of $y_{1}$ is much larger than $y_{\rm t2}=-4/5$,
again suggesting that the leading finite-size corrections arise from some
other source. The data are also well described when the term with $c_{0}$ is
included and $y_{1}=-0.46$ is fixed in Eq. (11), as shown in Table 2.
Corrections due to the subleading thermal field with $y_{\rm t2}=-4/5$ should
also exist, and, thus, there simultaneously exist corrections from different
sources. Nevertheless, to identify different correction terms is very
challenging for numerical analysis, particularly because our data are limited.
The fits of $B_{1}$ yield $d_{\textsc{B}}=1.792(2)$. To be conservative, we
quote our final estimate as $d_{\textsc{B}}=1.794(2)$, of which the
reliability is demonstrated by Fig. 4(b).

We did not simulate the Potts model for larger value of $Q$, due to severe
critical slowing-down and finite-size corrections.

###  III.2 Estimate from the loop model

For the O($n$) loop model, the backbone dimension $d_{\textsc{B}}$ can be
obtained from the finite-size scaling of quantities including
$B_{1},B_{\textsc{a}}$, and $S_{2}$, all of which give results consistent with
each other. Nevertheless, the estimates from $B_{2}$ has somewhat better
precision than those $B_{1}$ and $B_{\textsc{a}}$, probably because $B_{2}$
includes the sizes of all Ising domains in its definition. Table 3 summarizes
the fitting results from $B_{2}$ for
$n=1,\sqrt{2},\sqrt{3},\sqrt{2+\sqrt{3}}$, corresponding to
$Q=1,2,3,2+\sqrt{3}$. The fitting details are given below.

$Q$ | $L_{\rm m}$ | $d_{\textsc{B}}$ | $a$ | $b_{1}$ | $y_{1}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---|---
1 | 16 | 1.643 39(1) | 1.4922(2) | -1.529(9) | -1.082(3) | 4.2/6
| 32 | 1.643 40(2) | 1.4918(4) | -1.57(4) | -1.091(9) | 3.0/5
| 64 | 1.643 39(2) | 1.4921(6) | -1.5(1) | -1.08(2) | 2.8/4
| 32 | 1.643 35(1) | 1.4932(3) | -1.07(1) | -1 | 4.2/5
| 64 | 1.643 37(2) | 1.4928(4) | -1.02(3) | -1 | 2.1/4
2 | 32 | 1.732 27(3) | 1.1140(5) | -0.40(3) | -1.00(3) | 5.4/4
| 64 | 1.732 27(5) | 1.1139(8) | -0.4(1) | -1.01(8) | 5.3/3
| 32 | 1.732 27(1) | 1.1140(2) | -0.399(3) | -1 | 5.4/5
| 64 | 1.732 27(2) | 1.1140(3) | -0.400(8) | -1 | 5.4/4
3 | 16 | 1.793 9(1) | 0.898(3) | 0.105(4) | -0.42(3) | 5.4/5
| 32 | 1.793 8(2) | 0.900(4) | 0.11(1) | -0.46(6) | 4.9/4
| 32 | 1.793 66(3) | 0.9024(5) | 0.126(3) | -1/2 | 5.4/5
| 64 | 1.793 70(5) | 0.9018(8) | 0.130(5) | -1/2 | 4.4/4
2 | 8 | 1.838 3(2) | 0.757(3) | 0.210(2 | -0.304(6) | 4.5/6
$+$ | 16 | 1.838 4(2) | 0.755(5) | 0.211(2) | -0.30(1) | 4.3/5
$\sqrt{3}$ | 8 | 1.838 38(3) | 0.7546(4) | 0.2114(7) | -0.3 | 5.0/7
| 16 | 1.838 36(4) | 0.7549(5) | 0.211(1) | -0.3 | 4.3/6
| 32 | 1.838 32(5) | 0.7555(8) | 0.209(2) | -0.3 | 3.5/5

Table 3: Fitting results of $d_{\textsc{B}}$ from $B_{2}$ for the DP O$(n)$
loop model with $n=1,\sqrt{2},\sqrt{3},\sqrt{2+\sqrt{3}}$.

####  III.2.1 $n\\!=\\!1$

From Eq. (7) for the line of stable fixed points, one has the bond weight
$x\\!=\\!1$ for $n\\!=\\!1$, meaning that the coupling strength for the dual
Ising model $K^{*}\\!=\\!-\frac{1}{2}\ln x\\!=\\!0$ and thus the Ising spins
on different lattice sites are independent. With up (down) spins being
interpreted as occupied (empty) sites, the DP O(1) loop model at the branch
$x_{-}$ is just the site percolation on the triangular lattice. This is also
reflected in the induced-subgraph cluster algorithm, as formulated in Sec. II,
where the bond occupation probability $p=1-x=0$ and each site is randomly
occupied or empty with probability $50\%$.

Figure 6: Estimated backbone exponent $d_{\textsc{B}}$ from the domains of the
DP O($n$) loop model with $n=1$ (a) and $n=\sqrt{2+\sqrt{3}}$ (b). The O(1)
loop model, which has the bond weight $x=1$, is actually the site percolation
on the triangular lattice. The upward or downward bending of those curves,
with the $d_{\textsc{B}}$ value deviating from the estimated value,
illustrates the reliability of the quoted error bars.

In the least-squares fit of the second moment $B_{2}$ of bridge-free clusters
by Eq. (11), if $b_{2}=0$ and $b_{1}$, $y_{1}$ and $c_{0}$ are left as free
parameters, then we have $d_{\textsc{B}}=1.643\,38(4)$ but the correction
amplitude $b_{1}$ is consistent with zero within its error bar. With further
fixing $c_{0}=0$, we obtain $d_{\textsc{B}}=1.643\,39(2)$ and
$y_{1}=-1.09(1)$. The $B_{2}$ data can be also well described by Eq. (11) with
$b_{1}=b_{2}=0$ and $c_{0}$ as free parameter, and the backbone dimension is
$d_{\textsc{B}}=1.643\,35(5)$. The fitting results are listed in Table 3. By
considering all different fits, we take the final estimate as
$d_{\textsc{B}}=1.643\,39(5)$, of which the reliability is illustrated in Fig.
6. The fits for $B_{1}$ and $B_{\textsc{a}}$ give
$d_{\textsc{B}}=1.643\,37(9)$ and $d_{\textsc{B}}=1.643\,38(5)$, respectively,
which are consistent with that from $B_{2}$ but have slightly larger error
bars.

####  III.2.2 $n\\!=\\!\sqrt{2}$

As described in Sec. II, the current study of the honeycomb-lattice O($n$)
loop model is to simulate the generalized Ising model Eq. (8) on the
triangular lattice by the induced-subgraph cluster method. To determine the
backbone dimension $d_{\textsc{B}}$ for $Q=2$, one can either simulate the
standard Ising model at the critical branch $1/x_{+}=\sqrt{3}$ and study the
FK random clusters constructed in the cluster simulation or simulate the
O($\sqrt{2}$) loop model in the DP phase and study the Ising domains. Namely,
within the framework of the loop model, there exists a correspondence between
the FK random clusters constructed at the critical branch $x_{+}$ with $n=1$
and the Ising domains in the DP phase with $n^{*}=\sqrt{2}$.

We argue that such a geometric correspondence can be extended to general $n$.
Given the honeycomb-lattice loop model with parameters $(n,x)$, one can always
construct on the triangular lattice the Ising domains and FK-like random
clusters by placing occupied bonds with probability $1-x$ within each domain.
Along the critical branch $x_{+}(n)$, we assume the following duality
relation: Both the fractal dimensions of domains and of FK-like random
clusters, $d_{\textsc{face}}$ and $d_{\textsc{fk}}$, are given by $y_{\rm h1}$
in Eq. (4); nevertheless, the Coulomb-gas coupling $g$ for $d_{\textsc{face}}$
is in $[4,6]$ and calculated from $n$ as $n=-2\cos(\pi g/4)$, while $g_{*}$
for $d_{\textsc{fk}}$ is in $[2,4]$ and relates to $g$ as $g_{*}=16/g$. As a
consequence of this duality relation, we further argue that the FK-like random
clusters for $(n,x_{+})$ have the same geometric structures as the domains of
the DP O($n_{*}$) loop model with $n_{*}=-2\cos(\pi g_{*}/4)$, similarly to
the correspondence between the critical and tricritical $Q$-state Potts model
observed in Ref. [27]. For the marginal case $n=n_{*}=2$, this means that,
apart from multiplicative logarithmic corrections, the FK-like random clusters
and the domains have the same fractal dimension. Finally, it is noted that,
for $n<2$, the FK random clusters in the DP phase are too small to percolate.
A simple example is for $(n\\!=\\!1,x_{-}\\!=\\!1)$, where the bond
probability is $p\\!=\\!1-x=0$ and each FK-like cluster is just a single site.

For the critical Ising model, the induced-subgraph cluster method reduces to
the SW cluster algorithm, which still suffers from some critical slowing-down
due to the logarithmic divergence of the specific heat. Hereby, we determine
the backbone dimension $d_{\textsc{B}}$ from the O$(\sqrt{2})$ loop model at
the branch $x_{-}$, for which critical slowing-down is completely absent.

We perform the least-squares fits for the $B_{2}$ data to Eq. (11). By setting
$b_{2}\\!=\\!0$ and leaving the other correction parameters free, we obtain
$d_{\textsc{B}}\\!=\\!1.732\,19(5)$, and $c_{0}$ is consistent with zero
within its error bar. By further fixing $c_{0}\\!=\\!0$, we have
$d_{\textsc{B}}\\!=\\!1.732\,27(5)$ and $y_{1}\\!=\\!-1.01(8)$. The error bar
of $d_{\textsc{B}}$ can be further suppressed if $y_{1}\\!=\\!-1$ is taken.
Nevertheless, this smaller error bar cannot be taken into account seriously,
since we are not aware from what sources such finite-size corrections arise.
The fitting results are given in Table 3, and the final estimate is quoted as
$d_{\textsc{B}}=1.732\,27(8)$, of which the reliability for the central value
and the quoted error bar is illustrated in Fig. 5. The fits for $B_{1}$ and
$B_{\textsc{a}}$ give $d_{\textsc{B}}=1.732\,2(1)$ and
$d_{\textsc{B}}=1.732\,0(3)$, respectively, in excellent agreement with that
from $B_{2}$.

####  III.2.3 $n\\!=\\!\sqrt{3}$

The same fitting procedure is applied to the data of
$B_{1},B_{\textsc{a}},B_{2}$. It is now necessary to simultaneously include
two correction terms (e.g., with $b_{1}$ and $c_{0}$). For $n\\!=\\!\sqrt{3}$,
we obtain $d_{\textsc{B}}=1.793\,8(3)$ from $B_{2}$, $1.793\,5(7)$ from
$B_{1}$ and $1.794\,0(7)$ from $B_{\textsc{a}}$. The leading correction
exponent is estimated as $y_{1}=-0.42(3)$, which is in good agreement with
$y_{1}=-0.46(3)$ from the $Q=3$ Potts model. Both estimates of $y_{1}$ are
significantly larger than $y_{\rm t2}=-2/5$. This strongly suggests that the
leading corrections associated with backbone clusters are from some scaling
field with exponent near $\approx-0.4$. The fitting results are listed in
Table 3, the final estimate is taken as $d_{\textsc{B}}=1.793\,8(3)$ and shown
Fig. 5(b).

####  III.2.4 $n\\!=\\!\sqrt{2+\sqrt{3}}$

The Coulomb-gas coupling is $g=11/3$, leading to $y_{\rm
t2}\\!=\\!-4/11\\!\approx\\!-0.364$. With $b_{2}=0$, we obtain exponents
$d_{\textsc{B}}=1.838\,2(2)$ and $y_{1}=-0.304(6)$ for $L_{\rm m}=8$, and the
correction amplitudes are also well determined as $b_{1}=0.210(2)$ and
$c_{0}=-0.254(8)$. Although the estimated value of $y_{1}$ is larger than
$y_{\rm t2}$ if the error bar of $y_{1}$ is taken into account, it is
difficult to conclude $y_{1}\neq y_{\rm t2}$. Nevertheless, it is clear that,
as $n$ increases, finite-size corrections become larger and larger, albeit
they are smaller than those for the corresponding Potts model. As
$n\rightarrow 2$, we expect that correction exponent $y_{1}\rightarrow 0$ and
multiplicative logarithmic corrections would arise. The final estimate is
taken as $d_{\textsc{B}}=1.838\,4(5)$ and shown in Fig. 6(b).

###  III.3 Logarithmic corrections for $n=2$

The two branches of the O($n$) loop model meet at $x_{\pm}=1/\sqrt{2}$ for
$n=2$, and the thermal field associated with the bond weight $x$ becomes
marginal (i.e., $y_{\rm t2}=0$). In accordance with the phase diagram in Fig.
1, this thermal field is marginally irrelevant for $x>x_{\pm}$ and marginally
relevant for $x<x_{\pm}$, and, at $x_{\pm}$, its amplitude vanishes. It is
expected that logarithmic corrections, irrespective of multiplicative or
additive forms, are absent in thermodynamic quantities, including energy
density, specific heat, and magnetic susceptibility, as well as in geometric
quantities associated with the sizes of loops or domains.

$L_{\rm m}$ | $y_{D_{2}}$ | $a$ | $b_{1}$ | $y_{1}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---
4 | 1.750 000(8) | 0.89188(4) | 0.079(4) | -2.00(4) | 4.0/8
8 | 1.749 993(10) | 0.89192(5) | 0.12(5) | -2.2(2) | 2.7/7
4 | 1.750 000(6) | 0.89188(3) | 0.0787(6) | -2 | 4.0/9
8 | 1.750 001(7) | 0.89187(3) | 0.079(2) | -2 | 4.0/8
Table 4: Fitting results of the second moment $D_{2}$ of domain sizes for the
O$(2)$ loop model at $x_{\pm}=1/\sqrt{2}$.  Figure 7: Absence of
multiplicative and additive logarithmic corrections in geometric quantities
associated with domain sizes of the O$(2)$ loop model at $x_{\pm}=1/\sqrt{2}$,
as illustrated by the rescaled quantity $D_{2}/L^{7/4}$. The horizontal axis
is the log plot of linear size $L$. All the finite-size data of
$D_{2}/L^{7/4}$ are more or less consistent with the horizontal red dashed
line at $0.89188$, strongly indicating that corrections of logarithmic types
are absent and those of algebraic forms are very small.

To confirm this expectation, we sample the second moment of domain sizes as
$D_{2}=L^{-2}\langle\sum{\mathcal{D}}^{2}\rangle$, where the summation is over
all the dual Ising domains ${\mathcal{D}}$. The domains have a fractal
dimension $d_{\textsc{f}}=y_{\rm h1}$ given by Eq. (4), and, for $n=2$, it is
expected $D_{2}\sim L^{2d_{\textsc{f}}-2}=L^{7/4}$ with $d_{\textsc{f}}=15/8$.
By fixing $b_{2}=c_{0}=0$ so that only one single correction term is included
in Eq. (11), we observe that the reasonably good fits can already be obtained
with the minimum size as small as $L_{\rm m}=4$. The estimated exponent,
$d_{\textsc{f}}=1.875\,00(1)$, is excellent agreement with the exact value
$15/8$. The correction exponent $y_{1}=-2.00(4)$ is well consistent with $-2$,
and, furthermore, the correction amplitude $b_{1}=0.079(4)$ is very small; see
Table 4 for details. For clarity, the rescaled data of $D_{2}/L^{7/4}$ are
plotted versus size $L$ in Fig. 7, and it can be seen that, as long as for
$L\geq 128$, finite-size corrections are buried in statistical noise. The
rapid convergence of $D_{2}/L^{7/4}$ demonstrates the absence of
multiplicative or additive logarithmic corrections for $n=2$. It also strongly
supports that, along the line $x_{-}(n)$ of stable fixed points, finite-size
corrections with $y_{\rm t2}$ should not occur in quantities associated with
loops and domains.

${\mathcal{O}}$ | $L_{\rm m}$ | $d_{\textsc{B}}$ | $a$ | $b_{1}$ | $y_{1}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---|---
$B_{2}$ | 16 | 1.867 1(3) | 0.733(6) | 0.230(3) | -0.29(1) | 3.7/5
| 32 | 1.866 7(4) | 0.740(9) | 0.230(2) | -0.30(2) | 2.8/4
| 64 | 1.866 9(8) | 0.74(2) | 0.229(3) | -0.29(4) | 2.8/3
$B_{\textsc{a}}$ | 32 | 1.869(1) | 0.522(9) | 0.051(8) | -0.20(4) | 1.5/4
| 64 | 1.868(1) | 0.531(10) | 0.046(5) | -0.25(8) | 0.9/3
$B_{1}$ | 16 | 1.866 4(6) | 0.780(6) | 0.149(3) | -0.30(2) | 3.7/5
| 32 | 1.866 2(9) | 0.782(10) | 0.149(2) | -0.31(4) | 3.6/4

Table 5: Fitting results for backbone exponent $d_{\textsc{B}}$ for the O$(2)$
loop model without logarithmic corrections–i.e., Eq. (11).

For the backbone exponent $d_{\textsc{B}}$, we first try to fit the data to
Eq. (11), without logarithmic corrections being included, and the results are
shown Table 5. Taking $B_{2}$ for an example, if we fix $b_{2}=0$ and leave
other parameters free, then no stable results can be obtained. When $c_{0}=0$
is further fixed, we have $d_{\textsc{B}}=1.866\,9(8)$. Similar analysis for
$B_{\textsc{a}}$ and $B_{1}$ gives $d_{\textsc{B}}=1.868(2)$ and $1.866(2)$,
respectively. As shown in Table 5, all the estimates are consistent with each
other. The estimated value of $d_{\textsc{B}}$ is slightly smaller than the
fractal dimension $d_{\textsc{f}}=15/8$.

In the earlier study [21], it was numerically observed that the backbone
exponent $d_{\textsc{B}}$ of the tricritical Potts model reduces to the
fractal dimension $d_{\textsc{f}}$ of FK random clusters, which can be
explained by the fact that the red-bond scaling field is irrelevant (i.e.,
$d_{\rm red}<0$ for $g>4$). For the O(2) loop model with $g=4$, the red-bond
scaling field is marginal $d_{\rm red}=0$, and thus,
$d_{\textsc{B}}=d_{\textsc{f}}=15/8$ might still be expected. Nevertheless,
the fitting results of $d_{\textsc{B}}$ in Table 5 deviate from the predicted
value systematically when the error bars are taken into account. A plausible
reason is due to logarithmic corrections, which may still occur in backbone
clusters. This is also implicitly reflected by the small correction exponent
$y_{1}\\!\approx\\!-0.2$, as shown in Table 5, which might actually correspond
to logarithmic corrections.

${\mathcal{O}}$ | $L_{\rm m}$ | $d_{\textsc{B}}$ | $a$ | $b$ | $\hat{y}_{{\mathcal{O}}}$ | $d_{0}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---|---|---
$B_{2}$ | 64 | 1.875 | 1.171(3) | 0.06(3) | -0.246(2) | 1.928(4) | 2.3/3
| 128 | 1.875 | 1.171(8) | 0.07(9) | -0.246(4) | 1.929(9) | 0.5/2
| 128 | 1.875 6(3) | 1.163(2) | -0.12(8) | -0.25 | 1.66(2) | 2.1/2
$B_{\rm A}$ | 16 | 1.875 | 0.694(4) | -0.042(2) | -0.126(2) | 4.2(1) | 5.1/6
| 32 | 1.875 | 0.700(7) | -0.038(5) | -0.129(3) | 4.4(3) | 4.3/5
| 16 | 1.874 93(8) | 0.6919(9) | -0.042(2) | -0.125 | 4.13(7) | 5.0/6
| 32 | 1.874 8(1) | 0.693(2) | -0.039(4) | -0.125 | 4.2(1) | 4.2/5
$B_{1}$ | 16 | 1.875 | 0.960(6) | -0.01(1) | -0.120(2) | 0.7(1) | 4.3/6
| 32 | 1.875 | 0.956(10) | -0.03(3) | -0.118(3) | 0.6(2) | 4.0/5
| 64 | 1.875 | 0.96(2) | -0.01(6) | -0.120(6) | 0.7(4) | 3.9/4
| 16 | 1.875 3(1) | 0.971(2) | -0.006(8) | -0.125 | 0.88(6) | 4.3/6
| 32 | 1.875 4(2) | 0.970(3) | -0.01(2) | -0.125 | 0.8(1) | 4.1/5
| 64 | 1.875 3(3) | 0.972(5) | 0.01(4) | -0.125 | 0.9(2) | 3.9/4
Table 6: Fitting results of $d_{\textsc{B}}$ for the O$(2)$ loop model with
the logarithmic ansatz (12) and $y_{1}=-1$.

We assume that quantities on the basis of backbone clusters follow the
following finite-szie scaling ansatz

${\mathcal{O}}=L^{y_{{\mathcal{O}}}}[\ln
L+d_{0}]^{\hat{y}_{{\mathcal{O}}}}(a+bL^{y_{1}})\;,$ (12)

where the constant $d_{0}$ is commonly used in the fitting with multiplicative
logarithmic corrections. In comparison with the general scaling ansatz for
logarithmic corrections, we have ignored in Eq. (12) additive logarithmic
corrections which appear in powers of $1/\ln L$. This is because, in the fit,
(i) the effects of parameter $d_{0}$ and of those additive logarithmic
corrections can interfere with each other, and (ii) our data are limited and
can already be reasonably described without multiplicative logarithmic
corrections.

Since no prediction exists for the logarithmic exponent
$\hat{y}_{{\mathcal{O}}}$ for backbone clusters and it is difficult to get an
estimate to $y_{{\mathcal{O}}}$ and $\hat{y}_{{\mathcal{O}}}$ simultaneously,
we first fix $y_{{\mathcal{O}}}$ with expected value and get the estimate of
$\hat{y}_{{\mathcal{O}}}$. For $B_{2}$, we first fix
$y_{{\mathcal{O}}}=2d_{\textsc{B}}-2=7/4$ and leave all other parameters free,
which leads to no stable results. Then we fix $y_{1}=-1$ and obtain
$\hat{y}_{{\mathcal{O}}}=-0.246(2)$ for $L_{\rm m}=64$. More fits have been
done with $y_{1}$ fixed to various values in the interval $[-2,-1)$, which
give the estimate $\hat{y}_{{\mathcal{O}}}=0.24(1)$. Later, we fix $y_{1}=-1$
and $\hat{y}_{{\mathcal{O}}}$ to various values in the range
$-0.23\leq\hat{y}_{{\mathcal{O}}}\leq-0.25$, and we obtain the estimate
$d_{\textsc{B}}=1.875\,3(6)$ by covering all the fitting results. Following
the same procedure, we obtain the estimate $d_{\textsc{B}}=1.874\,7(7)$ and
$\hat{y}_{{\mathcal{O}}}=0.135(15)$ for $B_{\rm A}$,
$d_{\textsc{B}}=1.875\,2(7)$ and $\hat{y}_{{\mathcal{O}}}=0.122(8)$ for
$B_{1}$. The detail of the fitting results are summarized in Table. 6. Thus,
we expect $(y_{{\mathcal{O}}},\hat{y}_{{\mathcal{O}}})$ is $(15/8,1/8)$ for
$B_{\textsc{a}}$ and $B_{1}$, and is $(7/4,1/4)$ for $B_{2}$. This is further
demonstrated in Fig. 8, where the approximately straight lines in the log-log
plot correspond to logarithmic exponents $\hat{y}_{{\mathcal{O}}}=-1/4$ and
$-1/8$.

In the four-state Potts model, the multiplicative logarithmic correction
exponent is $-1/16$ for the magnetization [40], and thus, the size of the
largest FK cluster also has logarithmic exponent $-1/16$, half of our estimate
$-1/8$ for backbone clusters. This is an interesting observation, but it might
be due to the absence of additive logarithmic corrections in Eq. (12).
Unfortunately, it is very challenging to obtain conclusive evidence from fits
simultaneously with multiplicative and additive logarithmic corrections.
Further theoretical insights are needed for the exact value of the logarithmic
exponent for backbone clusters at $n=2$.

Figure 8: Log-log plot of rescaled backbone quantities versus $(\ln L+d_{0})$
for the O$(2)$ loop model. The rescaled quantities are $B_{2}/L^{7/4}$,
$B_{1}/L^{15/8}$, and $B_{\textsc{a}}/L^{15/8}$. For being concise,
$B_{\textsc{a}}$ is multiplied by a constant such that the data for
$B_{\textsc{a}}$ and $B_{1}$ approximately collapse on top of each other. The
existence of multiplicative logarithmic corrections and the corresponding
exponents are clearly illustrated.

###  III.4 Conjectured formula for leading correction exponent in backbone
clusters

The fitting results of $d_{\textsc{B}}$ from the $Q$-state Potts and the DP
O($n$) loop model suggest that, in geometric observables associated with
backbone clusters, the leading correction exponent is
$y_{1}\\!\approx\\!-1.0,-0.6,-0.4,-0.2$ for $Q=1,2,3,2+\sqrt{3}$, respectively
[the estimate $y_{1}\\!\approx\\!-1$ for the O($\sqrt{2}$) loop model might be
due to the small amplitude of the correction term]. For $Q=4$ or $n=2$, there
exists a marginal scaling field, leading to multiplicative and additive
logarithmic corrections. These corrections cannot be accounted by the scaling
field with $y_{\rm t2}$. It is tempting to identify the exact value of $y_{1}$
from the Kac table in the conformal field theory [41].

The Kac table is characterized by a pair of integers $i,j\geq 0$. The value of
a scaling dimension $X_{i,j}$, which relates to the corresponding exponent $y$
as $y=2-X$, is expressed as [42, 39]

$X_{i,j}=\frac{[i(m+1)-jm]^{2}-1}{2m(m+1)}\;,$ (13)

where the conformal anomaly $m$ can be calculated from the Coulomb-gas
coupling as $m=g/(4-g)$ for $g\in[2,4]$ and $m=4/(g-4)$ for $g\in[4,6]$. For
the critical Potts model, this gives $m=2,3,5,\infty$ for $Q\\!=\\!1,2,3,4$,
respectively. The exactly known exponents, in Eqs. (4) and (I) can be mostly
identified in the Kac table with integer $(i,j)$. An exception is the leading
magnetic dimension $X_{\rm h1}$, of which the indices
$i\\!=\\!j\\!=\\!(m\\!+\\!1)/2$ are half integers for $m=2$ (percolation).

On the basis of the estimated values of $y_{1}$, we conjecture that the
corresponding indices are $(i,j)=(2,0)$, leading to

$\displaystyle y_{1}$ $\displaystyle=$ $\displaystyle-(4m+3)/2m(m+1)\;,$ (14)
$\displaystyle=$ $\displaystyle-(4-g)(g+12)/8g\;.$

This predicts $y_{1}\\!=\\!-11/12,-5/8,-23/60,-47/263$ for $Q\\!=\\!1,2$,
$3,2\\!+\\!\sqrt{3}$, respectively, which are surprisingly close to the
numerical results of $y_{1}$. Moreover, Eq. (14) gives
$y_{1}\\!\rightarrow\\!0$ as $g\\!\rightarrow\\!4$, and thus can explain the
appearance of logarithmic corrections for $n=2$.

$Q$ | $L_{\rm m}$ | $d_{\textsc{B}}$ | $a$ | $b_{1}$ | $b_{2}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---|---
$1$ | 64 | 1.643 34(2) | 1.4935(5) | -0.67(3) | -5.8(8) | 2.0/4
| 128 | 1.643 33(3) | 1.4937(8) | -0.69(7) | -5(4) | 2.0/3
$2$ | 16 | 1.732 20(2) | 1.1153(3) | -0.051(3) | -0.57(1) | 8.3/6
| 32 | 1.732 17(3) | 1.1159(6) | -0.059(6) | -0.53(3) | 6.3/5
| 64 | 1.732 20(5) | 1.1154(9) | -0.05(1) | -0.6(1) | 5.6/4
$3$ | 64 | 1.794 1(1) | 0.894(2) | 0.108(9) | -0.04(2) | 6.8/4
| 128 | 1.793 8(2) | 0.899(3) | 0.07(2) | 0.06(6) | 3.2/3
| 256 | 1.794 0(3) | 0.895(6) | 0.10(4) | -0.0(1) | 2.5/2
$2$ | 32 | 1.840 0(3) | 0.714(5) | 0.14(1) | 0.102(7) | 1.6/5
$+$ | 64 | 1.840 0(4) | 0.715(9) | 0.14(2) | 0.10(1) | 1.6/4
$\sqrt{3}$ | 128 | 1.839 3(7) | 0.73(2) | 0.10(4) | 0.13(3) | 0.3/3
Table 7: Fitting results of $d_{\textsc{B}}$ from $B_{2}$ for the DP O($n$)
loop model with $n=1,\sqrt{2},\sqrt{3},\sqrt{2+\sqrt{3}}$. The leading
correction exponent $y_{1}$ is fixed at the value given by Eq. (14) and the
subleading correction exponent is taken as $y_{2}=2y_{1}$.

For the critical $Q$-state Potts model, the external-perimeter fractal
dimension $d_{\textsc{ep}}$ has indices as $(i,j)=(1,0)$ in the Kac table,
which is neighboring to $(i,j)=(2,0)$ conjectured for $y_{1}$. It is noted
that $d_{\textsc{ep}}$ is also the fractal dimension of external perimeters
and hulls of backbone clusters [43]. This observation indirectly supports that
the conjectured formula (14) might be reasonable for backbone clusters.

On this basis, we reanalyze the $B_{2}$ data for the DP O$(n)$ loop model by
Eq. (11) with $y_{1}$ being fixed at the predicted value. Another correction
term with $b_{2}$ or $c_{0}$ is also included, and, for simplicity, we first
fix the subleading correction exponent as $y_{2}=-2$. For $n=1$, the fit with
$c_{0}=0$ gives $d_{\textsc{B}}=1.643\,33(3)$, and the fit with $b_{2}=0$
leads to $d_{\textsc{B}}=1.643\,35(5)$, consistent with each other. The same
procedure has been done for $n=\sqrt{2}$ and $\sqrt{3}$, and the results are
consistent with the previous estimates. For $n=\sqrt{2+\sqrt{3}}$, the
residual $\chi^{2}$ is big until $L_{\rm m}=128$, and the result is
$d_{\textsc{B}}=1.841\,7(5)$ for $c_{0}=0$. When fixing $b_{2}=0$ and leaving
$c_{0}$ free, we obtain $d_{\textsc{B}}=1.841\,5(5)$. Both estimates deviate
from the previous estimate $1.838\,5(4)$. This is probably due to that
subleading finite-size corrections cannot be described by exponent $y_{2}=-2$.

We then set $y_{2}=2y_{1}$ and fix $c_{0}=0$. The final estimates
$d_{\textsc{B}}=1.643\,33(3)$, $1.732\,20(6)$, $1.793\,9(4)$, $1.839\,5(9)$
are obtained for $n=1$, $\sqrt{2}$, $\sqrt{3}$, and $\sqrt{2+\sqrt{3}}$,
respectively, consistent with previous estimates. The fitting results with
this choice are summarized in Table 7. For $n=\sqrt{2}$, we have
$y_{2}=2y_{1}=-5/4$, close to $-1$. The amplitude $b_{1}$ is much smaller than
$b_{2}$, and thus, for small system sizes, the subleading corrections can
surpass the leading ones. This might explain that, in the previous fits, the
estimated correction exponent is $y_{1}\approx-1$ in Table 3. For small $n$,
the precision of $d_{\textsc{B}}$ with this choice is slightly improved.

$Q$ | $L_{\rm m}$ | $d_{\rm min}$ | $a$ | $b_{1}$ | $y_{1}$ | $\chi^{2}/{\rm DF}$
---|---|---|---|---|---|---
| 32 | 1.094 42(7) | 1.1221(6) | -2.01(3) | -0.943(5) | 2.3/4
| 64 | 1.094 5(1) | 1.1212(9) | -2.10(8) | -0.96(1) | 1.1/3
$2$ | 64 | 1.094 54(3) | 1.1209(3) | -2.132(7) | -0.96 | 1.2/4
| 128 | 1.094 53(5) | 1.1210(4) | -2.14(2) | -0.96 | 1.1/3
| 32 | 1.067 45(9) | 1.0833(8) | -1.82(3) | -0.911(6) | 2.7/4
| 64 | 1.067 5(1) | 1.082(1) | -1.9(1) | -0.92(2) | 2.1/3
$3$ | 32 | 1.067 58(3) | 1.0822(2) | -1.865(3) | -0.92 | 5.0/5
| 64 | 1.067 52(5) | 1.0826(4) | -1.878(9) | -0.92 | 2.1/4
| 128 | 1.067 52(7) | 1.0826(6) | -1.88(2) | -0.92 | 2.1/3
| 32 | 1.047 4(1) | 1.0332(9) | -1.71(3) | -0.907(8) | 4.2/4
$2$ | 64 | 1.047 5(2) | 1.033(2) | -1.7(1) | -0.91(2) | 4.1/3
$+$ | 32 | 1.047 49(4) | 1.0328(3) | -1.724(4) | -0.91 | 4.4/5
$\sqrt{3}$ | 64 | 1.047 47(6) | 1.0330(4) | -1.73(1) | -0.91 | 4.1/4
| 128 | 1.047 49(9) | 1.0328(7) | -1.72(2) | -0.91 | 4.0/3
| 32 | 1.031 9(1) | 0.981(1) | -1.46(3) | -0.874(9) | 4.1/4
$4$ | 64 | 1.032 1(2) | 0.980(2) | -1.6(1) | -0.90(2) | 2.6/3
| 64 | 1.032 15(7) | 0.9795(5) | -1.59(1) | -0.9 | 2.6/4
| 128 | 1.032 2(1) | 0.9794(7) | -1.58(2) | -0.9 | 2.6/3
Table 8: Fitting results of the shortest-path exponent $d_{\rm min}$ from
$S_{1}$ for the DP O($n$) loop model with
$n=\sqrt{2},\sqrt{3},\sqrt{2+\sqrt{3}},2$.

##  IV Results for shortest-path exponent

Since the shortest-path exponent $d_{\rm min}$ of percolation ($Q\\!\to\\!1$)
is already estimated to have a high precision [26], we hereby present the
results from the domains of the DP O($n$) loop model with
$n\\!=\\!\sqrt{2},\sqrt{3},\sqrt{2+\sqrt{3}}$, and 2. We perform the least-
squares fits to the $S_{1}$ data by the ansatz (11). For $n=\sqrt{2}$, we
first leave all parameters free but cannot get stable results. Then, we fix
$c_{0}=b_{2}=0$ and obtain the estimate $d_{\rm min}=1.094\,48(13)$ and
$y_{1}\approx-0.96$, close to $-1$. If we fix $y_{2}=-2$ and leave $b_{2}$
free, then the fit gives the consistent estimate as $d_{\rm min}=1.094\,5(1)$
and $y_{1}=-0.97(2)$. In addition, fixing $y_{1}=-1$ and leaving $c_{0}$ free
gives the estimate $d_{\rm min}=1.0945\,5(1)$. We take the final estimate
$d_{\rm min}=1.094\,5(2)$ by considering the systematic error from different
fits. Following the same procedure for $n=\sqrt{3},\sqrt{2+\sqrt{3}},2$, we
obtain $d_{\rm min}=1.067\,5(3),1.047\,5(3),1.032\,2(4)$, respectively. The
details of fits are summarized in Table 8.

Figure 9: Estimated shortest-path exponent $d_{\rm min}$ from the domains of
the DP O($n$) loop model with $n=\sqrt{2}$ (a) and $n=\sqrt{3}$ (b). The
upward or downward bending of those curves, with the $d_{\rm min}$ value
deviating from the estimated value, illustrates the reliability of the quoted
error bars.

As shown in Table 8, we find that the leading correction exponents $y_{1}$ is
always close to $-1$, irrespective of the value of $n$, and, thus, that the
precision of $d_{\rm min}$ remains approximately unchanged as $n$ increases.
Unlike for the backbone exponent $d_{\textsc{B}}$, finite-size corrections do
not become stronger as $n$ increases. In particular, no logarithmic
corrections seem to arise for $n=2$. In Figs. 9 and 10, we plot
$S_{1}/L^{d_{\rm min}}$ versus $L^{-1}$. The approximate linearity of the red
line and the upward and downward bending curves reflect the reliability of the
estimates and the quoted error bars.

Figure 10: Estimated shortest-path exponent $d_{\rm min}$ from the domains of
the DP O($n$) loop model with $n=\sqrt{2+\sqrt{3}}$ (a) and $n=2$ (b). The
upward or downward bending of those curves, with the $d_{\rm min}$ value
deviating from the estimated value, illustrates the reliability of the quoted
error bars.

##  V Discussion

By cluster Monte Carlo methods, we carry out extensive simulations of the
$Q$-state Potts model on the square lattice and of the O$(n)$ loop model on
the honeycomb lattice, up to linear size $L=8192$. The induced-subgraph
picture adds a valuable perspective to understand the celebrated Swendsen-Wang
algorithm and further provides a versatile platform to formulate efficient
cluster or worm-type algorithms for loop models.

The backbone exponent $d_{\textsc{B}}$ and the shortest-path exponent $d_{\rm
min}$ are determined for the Fortuin-Kasteleyn random clusters of the critical
Potts model and for the domains of the DP O($n$) loop model. According to
Coulomb gas theory, the DP O$(n)$ loop model corresponds to the critical
$Q$-state Potts model with $Q=n^{2}$. The excellent agreement of the
$d_{\textsc{B}}$ and $d_{\rm min}$ results between the two systems strongly
supports that, in addition to the overall fractal structure, the domains of
the DP loop model and the corresponding FK random clusters share the same
scaling for other geometric properties.

In comparison with the $Q$-state Potts model, the study of the DP O$(n)$ loop
model benefits from the advantages that the absence of critical slowing down
and the absence of finite-size corrections with exponent $y_{\rm t2}$ in Eq.
(4). Also, an improved algorithm for identifying backbones is formulated, by
adopting a depth-first cluster growth procedure and introducing an auxiliary
array of treelike structure.

As a consequence of these advantages, the estimates of $d_{\textsc{B}}$ and
$d_{\rm min}$ achieve a precision better than the existing results. It can be
seen that, as $Q$ increases, the backbone exponent $d_{\textsc{B}}$ increases
and the shortest-path exponent $d_{\rm min}$ decreases, and, for $Q=4$, the
backbone exponent and the fractal dimension coincide (apart from logarithmic
corrections). This reflects that the critical FK random clusters become more
and more compact as $Q$ increases. Since the exact values of $d_{\textsc{B}}$
and $d_{\rm min}$ are still unknown for the two-dimensional Potts model, our
high-precision results can provide a solid testing ground for future
theoretical explorations.

Finite-size corrections in the DP O$(n)$ loop model are also studied along the
line of stable fixed points. It is confirmed that corrections with exponent
$y_{\rm t2}$ are absent. Particularly, for the marginal case $n=2$,
multiplicative and additive logarithmic corrections are not observed in
quantities associated with the domain sizes and in the shortest path. The
leading correction exponent $y_{1}$ for the shortest path is observed to be
roughly around $-1$, independent of $n$. For backbone clusters, however,
finite-size corrections become larger and larger as $n$ increases, and
logarithmic corrections arise for $n=2$. To account for these corrections,
which are beyond the description of exponent $y_{t2}$, we conjecture an exact
formula, Eq. (14), for the correction exponent $y_{1}$. The predicted value of
$y_{1}$ is consistent with our Monte Carlo data, and can explain the emergence
of logarithmic corrections for $n=2$. Fixing the $y_{1}$ value can lead to
smaller error bars of $d_{\textsc{B}}$, which are, however, not taken in our
final estimate. It remains as an interesting question to explore whether our
conjecture holds true.

##  VI Acknowledgments

Y.D. acknowledges support by the National Natural Science Foundation of China
(under Grant No. 11625522), the Science and Technology Committee of Shanghai
(under grant No. 20DZ2210100), and the National Key R&D Program of China
(under Grant No. 2018YFA0306501). W.Z. acknowledges support by the National
Natural Science Foundation of China Youth Fund (under Grant No. 12105133) and
the Fujian Provincial Natural Science Foundation of China (under Grant No.
2021J011030).

## References

  * Potts [1952] R. B. Potts, Some generalized order-disorder transformations, in  _Mathematical Proceedings of the Cambridge Philosophical Society_ , Vol. 48 (Cambridge University Press, Cambridge, UK, 1952), pp. 106–109.
  * Wu [1982] F. Y. Wu, The Potts model, Reviews of Modern Physics 54, 235 (1982).
  * Kasteleyn and Fortuin [1969] P.W. Kasteleyn and C.M. Fortuin, Phase transitions in lattice systems with random local properties, Physical Society of Japan Journal Supplement 26, 11 (1969).
  * Grimmett [2006] G. Grimmett,  _The Random-Cluster Model_ , edited by H. Kesten (Springer, Berlin, 2006).
  * Essam and Tsallis [1986] J. W. Essam and C. Tsallis, The Potts model and flows I. The pair correlation function, Journal of Physics A: Mathematical and General 19, 409 (1986).
  * Wu [1988] F. Y. Wu, Potts model and graph theory, Journal of Statistical Physics 52, 99 (1988).
  * Zhang  _et al._ [2020] L. Zhang, M. Michel, E. M. Elçi, and Y. Deng, Loop-cluster coupling and algorithm for classical statistical models, Physical Review Letters 125, 200603 (2020).
  * Francesco  _et al._ [2012] P. Francesco, P. Mathieu, and D. Sénéchal,  _Conformal Field Theory_ (Springer Science & Business Media, New York, 2012).
  * Kager and Nienhuis [2004] W. Kager and B. Nienhuis, A guide to stochastic Löwner evolution and its applications, Journal of Statistical Physics 115, 1149 (2004).
  * Cardy [2005] J. Cardy, SLE for theoretical physicists, Annals of Physics 318, 81 (2005).
  * Swendsen and Wang [1987] R. H. Swendsen and J. S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Physical Review Letters 58, 86 (1987).
  * Wolff [1989] U. Wolff, Collective Monte Carlo updating for spin systems, Physical Review Letters 62, 361 (1989).
  * Herrmann and Stanley [1984] H. J. Herrmann and H. E. Stanley, Building blocks of percolation clusters: Volatile fractals, Physical Review Letters 53, 1121 (1984).
  * Herrmann  _et al._ [1984] H. J. Herrmann, D. Hong, and H. Stanley, Backbone and elastic backbone of percolation clusters obtained by the new method of ‘burning’, Journal of Physics A: Mathematical and General 17, L261 (1984).
  * Stauffer and Aharony [2018] D. Stauffer and A. Aharony,  _Introduction to Percolation Theory_ (CRC Press, Boca Raton, FL, 2018).
  * Baxter [2016] R. J. Baxter,  _Exactly solved models in statistical mechanics_ (Elsevier, 2016).
  * Wu [2009] F. Y. Wu,  _Exactly Solved Models: A Journey In Statistical Mechanics-Selected Papers With Commentaries (1963–2008)_ (World Scientific, 2009).
  * Nienhuis [1984] B. Nienhuis, Critical behavior of two-dimensional spin models and charge asymmetry in the Coulomb gas, Journal of Statistical Physics 34, 731 (1984).
  * Nienhuis [1987] B. Nienhuis, Coulomb gas formulation of two-dimensional phase transitions, in  _Phase transitions and Critical Phenomena_ , Vol. 11, edited by C. Domb and J. L. Lebowitz (Academic Press, London, 1987) pp. 1–51.
  * Coniglio [1989] A. Coniglio, Fractal structure of Ising and Potts clusters: Exact results, Physical Review Letters 62, 3054 (1989).
  * Deng  _et al._ [2004a] Y. Deng, H. W. J. Blöte, and B. Nienhuis, Backbone exponents of the two-dimensional $q$-state Potts model: A Monte Carlo investigation, Physical Review E 69, 026114 (2004a).
  * Deng  _et al._ [2010] Y. Deng, W. Zhang, T. M. Garoni, A. D. Sokal, and A. Sportiello, Some geometric critical exponents for percolation and the random-cluster model, Physical Review E 81, 020102(R) (2010).
  * Xu  _et al._ [2014] X. Xu, J. Wang, Z. Zhou, T. M. Garoni, and Y. Deng, Geometric structure of percolation clusters, Physical Review E 89, 012120 (2014).
  * Hou  _et al._ [2019] P. Hou, S. Fang, J. Wang, H. Hu, and Y. Deng, Geometric properties of the Fortuin-Kasteleyn representation of the Ising model, Physical Review E 99, 042150 (2019).
  * Elçi  _et al._ [2016] E. M. Elçi, M. Weigel, and N. G. Fytas, Bridges in the random-cluster model, Nuclear Physics B 903, 19 (2016).
  * Zhou  _et al._ [2012] Z. Zhou, J. Yang, Y. Deng, and R. M. Ziff, Shortest-path fractal dimension for percolation in two and three dimensions, Physical Review E 86, 061101 (2012).
  * Deng  _et al._ [2004b] Y. Deng, H. W. J. Blöte, and B. Nienhuis, Geometric properties of two-dimensional critical and tricritical Potts models, Physical Review E 69, 026123 (2004b).
  * Domany  _et al._ [1981] E. Domany, D. Mukamel, B. Nienhuis, and A. Schwimmer, Duality relations and equivalences for models with O($n$) and cubic symmetry, Nuclear Physics B 190, 279 (1981).
  * Thompson [2015] C. J. Thompson,  _Mathematical statistical mechanics_ (Princeton University Press, 2015).
  * Nienhuis [1982] B. Nienhuis, Exact critical point and critical exponents of O($n$) models in two dimensions, Physical Review Letters 49, 1062 (1982).
  * Nienhuis [1991] B. Nienhuis, Locus of the tricritical transition in a two-dimensional $q$-state Potts model, Physica A: Statistical Mechanics and Its Applications 177, 109 (1991).
  * Deng  _et al._ [2007] Y. Deng, T. M. Garoni, W. Guo, H. W. J. Blöte, and A. D. Sokal, Cluster simulations of loop models on two-dimensional lattices, Physical Review Letters 98, 120601 (2007).
  * Liu  _et al._ [2011] Q. Liu, Y. Deng, and T. M. Garoni, Worm Monte Carlo study of the honeycomb-lattice loop model, Nuclear Physics B 846, 283 (2011).
  * Li and Sokal [1989] X.-J. Li and A. D. Sokal, Rigorous lower bound on the dynamic critical exponents of the Swendsen-Wang algorithm, Physical Review Letters 63, 827 (1989).
  * Chayes and Machta [1998] L. Chayes and J. Machta, Graphical representations and cluster algorithms II, Physica A: Statistical Mechanics and Its Applications 254, 477 (1998).
  * Deng  _et al._ [2009] Y. Deng, X. Qian, and H. W. J. Blöte, Single-cluster dynamics for the random-cluster model, Physical Review E 80, 036707 (2009).
  * Liu  _et al._ [2012] Q. Liu, Y. Deng, T. M. Garoni, and H. W. Blöte, The O($n$) loop model on a three-dimensional lattice, Nuclear Physics B 859, 107 (2012).
  * Huang  _et al._ [2018] W. Huang, P. Hou, J. Wang, R. M. Ziff, and Y. Deng, Critical percolation clusters in seven dimensions and on a complete graph, Physical Review E 97, 022107 (2018).
  * Cardy [1987] J. Cardy, Conformal invariance, in  _Phase transitions and Critical Phenomena_ , Vol. 11, edited by C. Domb and J. L. Lebowitz (Academic Press, London, 1987) pp. 55–123.
  * Salas and Sokal [1997] J. Salas and A. D. Sokal, Logarithmic corrections and finite-size scaling in the two-dimensional 4-state Potts model, Journal of Statistical Physics 88, 567 (1997).
  * Henkel [1999] M. Henkel,  _Conformal invariance and critical phenomena_ (Springer Science & Business Media, 1999).
  * Friedan  _et al._ [1984] D. Friedan, Z. Qiu, and S. Shenker, Conformal invariance, unitarity, and critical exponents in two dimensions, Physical Review Letters 52, 1575 (1984).
  * Grossman and Aharony [1986] T. Grossman and A. Aharony, Structure and perimeters of percolation clusters, Journal of Physics A: Mathematical and General 19, L745 (1986).