# Annealing Dynamics of Regular Rotor Networks: Universality and Its Breakdown

András Grabarits Department of Physics and Materials Science, University of
Luxembourg, L-1511 Luxembourg, Luxembourg  Gaetano Sammartino Department of
Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg,
Luxembourg  Adolfo del Campo Department of Physics and Materials Science,
University of Luxembourg, L-1511 Luxembourg, Luxembourg  Donostia
International Physics Center, E-20018 San Sebastián, Spain

###### Abstract

The spin-vector Monte Carlo model is widely used as a benchmark for the
classicality of quantum annealers but severely restricts the time evolution.
The spin-vector Langevin (SVL) model has been proposed and tested as an
alternative, closely reproducing the real-time dynamics of physical quantum
annealers such as D-Wave machines in the dissipative regime. We investigate
the SVL annealing dynamics of classical $\mathrm{O}(2)$ rotors on regular
graphs, identifying universal features in the nonequilibrium dynamics when
changing the range of interactions and the topology of the graph. Regular
graphs with low connectance or edge density exhibit universal scaling dynamics
consistent with the Kibble-Zurek mechanism, which leads to a power-law
dependence of the density of defects and the residual energy as a function of
the annealing time. As the interaction range is increased, the power-law
scaling is suppressed, and an exponential scaling with the annealing time sets
in. Our results establish a universal breakdown of the Kibble-Zurek mechanism
in classical systems characterized by long-range interactions, in sharp
contrast with previous findings in the quantum domain.

††preprint: APS/123-QED

##  I Introduction

Quantum Annealing (QA) is a heuristic algorithm that exploits quantum effects
to solve optimization problems [1, 2, 3]. QA relies on the gradual
transformation of a quantum system from an initial and easily preparable state
to a target state representing the solution of the optimization problem.
Specifically, the solution is encoded as the ground state of the problem
Hamiltonian. Due to the presence of noise and decoherence, current quantum
annealing devices are best described as programmable many-body open quantum
systems [4, 5, 6, 7]. A unitary time evolution provides an accurate
description for short annealing times [8, 9] and shallow circuits in a digital
implementation. These devices offer an exciting platform to explore non-
equilibrium physics [9]. One central paradigm in this context is the Kibble-
Zurek (KZ) mechanism, describing the dynamics of classical continuous and
quantum phase transitions. It predicts that when a system is driven across a
critical point in a time scale $\tau_{Q}$, adiabaticity breaks down, leading
to the formation of topological defects at a density which varies as a
universal power law with the quench time [10, 11, 12]. In many-body spin
systems of relevance to QA, defects generally represent errors in the
optimization process that limit the preparation of low-energy states. The KZ
mechanism has been used as a benchmark in quantum simulators [13, 14, 15, 16,
17] and annealing devices [18, 19, 7, 20, 8, 9]. In this context, efforts to
identify and establish quantum signatures of the dynamics have focused on
ruling out the behavior of classical models embedded in a thermal bath [1, 21,
22, 23, 24, 25, 26, 27, 20, 28]. In addition, the KZ mechanism holds in
classical spin systems with time-continuous simulated annealing and Glauber
dynamics [29, 30, 31, 32], as well as with Monte Carlo updates [33, 7].

The conventional derivation of the KZ scaling, applicable to classical
systems, relies on spatial locality [34, 14, 35, 36], and yet, the KZ
mechanism has been shown to hold in quantum systems with long-range
interactions decaying as a power law. This is the case for nonintegrable
systems, such as the transverse-field Ising model with power-law interactions
[37, 38, 39] and related systems realizable with Rydberg atom quantum
simulators [16, 40], as well as in systems with a free-fermionic structure,
such as the long-range Kitaev model [41]. By contrast, fully connected systems
exhibit a power-law scaling with the quench time, which is generally
incompatible with KZ mechanism and better explained in terms of Landau-Zener
crossings and quasi-adiabatic approximations [42, 43, 44, 45].

The coupling of a quantum critical system to a bath need not destroy the KZ
behavior and can preserve it by simply renormalizing the equilibrium critical
exponents [46, 7]. However, environmental effects can give rise to anti-
Kibble-Zurek behavior, whereby the density of kinks increases with the
annealing time, as first predicted in short-range systems [47] and observed in
several D-Wave annealing devices [18, 19, 8]. The anti-Kibble-Zurek scaling
has been further established in the long-range Kitaev model [48, 49, 50] and
fully connected systems [51].

In the classical domain, early studies in ion trap systems showed that power-
law Coulomb interactions preserve the KZ scaling expected in short-range
systems [52, 53, 54], with deviations reported in experiments in small samples
[55, 56]. Likewise, the critical dynamics of the dipolar spin ice exhibits KZ
scaling [57, 32]. Recent experiments in systems with dipolar interactions
reported the universal KZ scaling when excitations are topological defects and
found an enhanced suppression with the driving time in the case of
nontopological excitations [58]. Deviations from the KZ scaling are expected
in the presence of coarsening dynamics [59, 60], as has been described in
systems with power-law interactions when crossing equilibrium [61, 62] and
nonequilibrium [63, 64] phase transitions.

Power-law interactions decaying with the distance are natural in physical
systems and can be studied in a variety of simulators using ultracold gases
and Rydberg atoms [65, 66]. Complex couplings involving nonlocal interactions
arise in computational problems, statistical physics, and the study of
networks. Recent experimental progress has demonstrated the realization of
programmable interactions, giving rise to a synthetic emergent geometry
different from the physical one [67].

In spite of this progress, the validity of the KZ mechanism as a function of
the interaction range in classical systems remains to be studied and is the
focus of our work. We consider the annealing dynamics of a regular rotor
network as a function of the connectance or edge density, i.e., the ratio
between the number of edges and the maximal possible value associated with a
fully connected graph. While a universal KZ scaling is observed in networks
with small edge density, our results establish the breakdown of the KZ scaling
as the connectance is increased, giving rise to the onset of an exponential
suppression of the density of excitations with the annealing time.

The theoretical framework for the classical Langevin annealing dynamics of
$O(2)$ rotor networks is introduced in Sec. II, which also discusses the
equilibrium relaxation time and magnetization. This sets the ground for the
study of KZ scaling in graphs with low connectance in Sec. III and its
breakdown at high edge density in Sec. IV. Universal signatures beyond the KZ
mechanism characterizing fluctuations of the defect and energy densities are
discussed in Sec. V. Section VI closes with a summary and outlook.

##  II  Time-evolution and equilibirum properties

This section summarizes the details of the classical stochastic dynamics of
rotor networks and the associated equilibrium statistical properties.

###  II.1 The Spin Vector Langevin model

The QA algorithm relies on the dynamics of a system generated by a Hamiltonian
of the form

$H(t)=A(t)H_{0}+B(t)H_{P},$ (1)

where $A(t)$ and $B(t)$ define the annealing protocol interpolating between a
trivial Hamiltonian $H_{0}$ and the problem Hamiltonian $H_{P}$. These
functions are chosen to satisfy the boundary conditions $A(0)=B(\tau_{Q})=1$
and $B(0)=A(\tau_{Q})=0$, where $\tau_{Q}$ is the duration of the process.
Here, $H_{0}$ and $H_{P}$ are usually Ising type Hamiltonians of the form

$\displaystyle H_{0}=-\sum_{i\in V}\sigma^{x}_{i},$ (2) $\displaystyle
H_{P}=-\sum_{(i,j)\in E}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}-\sum_{i\in
V}h_{i}\sigma_{i}^{z},$ (3)

with $\sigma^{z}_{i}$ and $\sigma^{x}_{i}$ denoting the $z$ and $x$ Pauli
matrices on the $i$-th site. The QA algorithm aims to find the solution to
combinatorial optimization problems, which are encoded in the ground state of
$H_{P}$. The classical hardness of these problems is characterized by both the
underlying graph topology, $\mathcal{G}(V,E)$, and the spin-spin interactions
$J_{ij}$ [68]. Here, $V$ is the set of vertices, and $E$ denotes the set of
edges accounting for the presence of interactions among the vertices. In the
course of the QA algorithm, the system starts from the easy-to-prepare ground
state of $H_{0}$, and under slow adiabatic dynamics, it ends up in the final
ground state. However, in quantum systems, the classical hardness of the given
optimization problem leads to small energy gaps at quantum critical points. As
a result, this leads to time scales similar to those required by classical
algorithms for successful QA processes.

Even though the understanding of the crossing of such critical points in
finite-range models is still in its infancy, valuable insights are offered by
the stochastic Langevin dynamics of classical systems [52, 53, 18, 7, 20, 8].
To this end, we explore the finite-range extension of the ring topology of
classical rotors, where the interactions are described by circulant graphs.
The corresponding spins are represented by classical spin vectors with
trigonometric functions replacing the Pauli operators as
$\sigma_{i}^{z}\rightarrow\sin\theta_{i}$ and
$\sigma_{i}^{x}\rightarrow\cos\theta_{i}$. The angles $\theta_{i}$ describe
classical planar $O(2)$ rotors characterized by the $N$-dimensional vector
$\bm{\theta}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$. Note that this
mapping has been justified to reproduce most of the low-energy physics in the
corresponding quantum spin systems with tunable power-law interactions, and
common benchmarks of classicality rely on it [21, 22, 23, 24, 26, 7, 27, 20,
28]. The corresponding finite-range interacting Hamiltonian is given by

$H(\bm{\theta},t)=-\frac{J(t)}{2}\sum_{i,j=1}^{N}A_{ij}\sin\theta_{i}\sin\theta_{j}-h(t)\sum_{i=1}^{N}\cos\theta_{i},$
(4)

where $A_{ij}$ denotes the elements of the adjacency matrix describing the
circulant graphs $\mathcal{G}(V,E)=C_{i_{N}}(r)$. Here, $N$ denotes the number
of vertices of the graph, i.e., the number of rotors,
$\sum_{i<j}\,A_{ij}=\lvert V\rvert=N$. The integer $r$ sets the interaction
range between spins,

$A_{ij}=\begin{cases}1&\text{if}\;\,|i-j|\leq r,\\\
0&\text{otherwise.}\end{cases}$ (5)

In Eq. (4), the first term takes the role of the $H_{P}$ problem Hamiltonian,
while the second is that of the $H_{0}$ transverse field Hamiltonian.
Circulant graphs provide the simplest generalization of the ring topology by
extending the interaction range to the $r$ nearest neighbors. From the point
of view of graph theory, the distance between two vertices corresponds to the
number of connecting edges in the shortest path linking them. This shrinks the
distances between vertices, leading to a non-trivial geometrical structure.
However, a one-dimensional interpretation provides an accurate description of
the non-equilibrium dynamics near the critical point.

Due to the homogeneity of the system, these graphs are $2r$-regular. For
simplicity, we will refer to $r$ as the regularity in the following. To
characterize these graphs independently of the system size, the graph density,
or connectance, is introduced: $c=\frac{2r}{N-1}$. The connectance
interpolates between the empty graph ($c=0$) and the complete graph $K_{N}$
with $c=1$ [69], as depicted in Fig. 1.

The corresponding dynamics can be implemented in different ways. The Spin-
Vector Monte Carlo (SVMC) model is a convenient method to implement classical
stochastic dynamics, capturing time evolution via discrete Monte Carlo steps
[22, 23]. While it remains a common benchmark for annealing devices [24, 26,
7, 27, 28], it is not suited to describe the real-time dynamics unambiguously.
Its evolution is highly restricted, being parameterized solely by temperature
and lacking the notion of dissipation. The Spin Vector Langevin (SVL) model
replaces the Monte Carlo steps with continuous-time Langevin dynamics that is
flexible enough to account for dissipative effects by including a friction
coefficient subject to the fluctuation-dissipation theorem [20]. It also
provides a remarkable improvement over the SVMC method giving access to the
real-time evolution of the system.

The SVL evolution is encoded in the set of stochastic coupled differential
equations

$\begin{split}&m_{i}\ddot{\theta_{i}}+\gamma_{i}\dot{\theta_{i}}+\frac{\partial
H(\bm{\theta},t)}{\partial\theta_{i}}+\xi_{i}(t)=0\,,\;\;i=1,2,\ldots,N\,\\\
&\frac{\partial
H(\bm{\theta},t)}{\partial\theta_{i}}=-J(t)\sum_{j=1}^{N}A_{ij}\cos\theta_{i}\sin\theta_{j}+h(t)\sin\theta_{i},\end{split}$
(6)

where $\xi_{i}(t)$ represents independent and identical Wiener processes with
zero mean, describing white noise on the $i$-th vertex. This noise generates
classical fluctuations modeling how the external thermal bath affects the
dynamics [7, 21, 22]. The effective mass of the rotors $m_{i}$ weights the
inertia term, and $\gamma_{i}$ denotes the local damping constant. In the
following, a homogeneous mass and damping are chosen, setting $m_{i}=m$ and
$\gamma_{i}=\gamma$. Then, the noise autocorrelation is
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=2m\gamma
k_{B}T\delta_{ij}\delta(t-t^{\prime})$, as required by the fluctuation-
dissipation theorem [70] for a thermal bath at temperature $T$. In this work,
natural units are adopted for simplicity, setting $k_{B}=1,\,m=1$, and for the
temperature, $T=0.001$ is chosen in dimensionless units.

The annealing schedule specifies the time dependence of the Hamiltonian (4),
and in real QA machines, it is determined by the balance of engineering
optimization and the physical limitations of the hardware. We choose a linear
protocol

$\quad J(t)=J_{0}\frac{t}{\tau_{Q}}\quad\text{and}\quad
h(t)=h_{0}\left(1-\frac{t}{\tau_{Q}}\right),$ (7)

where $J_{0}$ and $h_{0}$ are arbitrary constants controlling the strength of
$H_{P}$ and $H_{0}$, respectively. As for the numerical implementation, the
Langevin equations are integrated with the help of the multi-dimensional
explicit order $2.0$ weak scheme of Ref. [71, 20].

Figure 1:  Schematic representation of $O(2)$ networks associated with
circulant graphs of order $N=10$ and increasing connectance from left to
right, interpolating between the cycle graph $C_{N}$ and the complete graph
$K_{N}$.

###  II.2 Equilibrium critical exponents and order parameter

As in the ring topology with $r=1$, regular rotor networks undergo a
paramagnetic-ferromagnetic-like phase transition at values of $J(t)$ and
$h(t)$ for which the two competing terms in the Hamiltonian are of the same
order of magnitude for arbitrary $N$ and $r$. The system starts from the
paramagnetic phase, where the angles are approximately paramagnetically
aligned $\theta_{i}\ll 1$. Here, the Hamiltonian can be approximated as
$H(\bm{\theta})\approx-\frac{J(t)}{2}\sum_{i,j=1}^{N}A_{ij}\theta_{i}\theta_{j}-h(t)N$.
The corresponding linearized Langevin equations read

$m\ddot{\theta_{i}}+\gamma\dot{\theta_{i}}+h(t)\theta_{i}+J(t)\sum_{\begin{subarray}{c}j=-r\\\
j\neq 0\end{subarray}}^{r}\theta_{i+j}+\xi_{i}(t)=0\,,\;\forall i.$ (8)

From Eq. (8), both the critical point and the critical exponent of the
relaxation time can be determined in the different damping regimes [72, 52,
53, 20]. In the overdamped regime, using the approximation
$\gamma\dot{\theta}\gg\ddot{\theta}$, one finds

$\tau\sim\left|\theta/\dot{\theta}\right|\sim\frac{\gamma}{h(t)-2rJ(t)}\,,$
(9)

setting the critical point at $h(t)=2rJ(t)$. To make the corresponding
critical time independent of the connectance, $J_{0}=1$ and $h_{0}=2r$ are
chosen, yielding $t_{c}=\tau_{Q}/2$. In the underdamped case, the inertial
term $\ddot{\theta}$ dominates and the equilibrium relaxation time reads

$\tau\simeq\left|\theta/\ddot{\theta}\right|^{1/2}\simeq\left|\frac{m}{h(t)-2rJ(t)}\right|^{1/2}.$
(10)

These equations imply that the regularity does not affect the relaxation time
beyond shifting the critical point. The corresponding critical exponents obey
$z\nu=1$ in the overdamped regime and $z\nu=1/2$ in the underdamped regime. In
both cases, the natural choice for the control parameter is
$\epsilon(t)=h(t)-2rJ(t)$, so that the transition takes place at
$\epsilon(t_{c})=0$. Shifting the beginning of the annealing schedule by
$\tau_{Q}/2$ leads to $t_{c}\equiv 0$, convenient for further analysis. This
way, the annealing schedule extends over the interval of
$[-\tau_{Q}/2,\tau_{Q}/2]$, with the linearized control parameter given by
$\epsilon(t)=2r\,t/\tau_{Q}$ around the critical point. Consistently, near
$t_{c}$, $H_{P}$ and $H_{0}$ vary over the same energy scale, independently of
$N$ and $r$. The quadratic part scales as
$\sum_{i,j}A_{ij}\sin\theta_{i}\sin\theta_{j}\sim Nr$, while $H_{0}$ acquires
the same order of magnitude thanks to the additional scaling of $h_{0}$ in
$r$, $h_{0}\sum_{i}\cos\theta_{i}\sim Nr$.

Figure 2: Equilibirum correlation length for the overdamped and underdamped
regimes, captured by the spatial and ensemble averages of the correlators of
the $\sin\theta_{i}$ values. Panels show the scaling collapse for $(a)$ the
overdamped and $(b)$ underdamped regimes via rescaling the distance by the
regularity $r$, as motivated by the universal KZ limit. The inset displays the
divergence of the correlation length as a function of the proximity to the
critical point, quantified by $\epsilon$. The corresponding critical exponent
is independent of $r$, up to numerical precision. The results were obtained by
averaging over $10^{4}$ trajectories and with quench time set to
$\tau_{Q}=2560$ to ensure adiabatic time-evolution. Figure 3:  Evolution of
the average magnetization $\mathcal{M}_{z}$ as the order parameter for
different values of $c$ in the $(a)$ small, $(b)$ intermediate, and $(c)$
fully connected limits. The equilibrium dependence on the control parameter
was extracted from time evolutions in the adiabatic limit, with
$\tau_{Q}=2560$ for $N=400$ rotors and averaged over $10^{4}$ trajectories.
The critical properties observed in the chain topology are invariant with
respect to $c$, maintaining the behavior around the critical point.

Next, the equilibrium correlation length is investigated via the angle-angle
correlation function defined as

$G(d)=\overline{\langle\sin\theta_{i}\sin\theta_{i+d}\rangle_{\mathrm{traj}}}\approx\overline{\langle\theta_{i}\theta_{i+d}\rangle_{\mathrm{traj}}},$
(11)

where the last approximation is valid near the critical point. Here, the
expressions $\langle\cdots\rangle_{\mathrm{traj}}$ and $\overline{\cdots}$
stand for the average over the stochastic trajectories and the spatial
coordinates, respectively. As shown in App. A, this allows for a compact
analytical treatment in terms of a Hamiltonian quadratic in the angle
variables $\theta_{i}$, yielding for the correlator

$G(d)\sim r^{-2}e^{-d/\xi_{\mathrm{eq}}},\quad\xi_{\mathrm{eq}}(r)\sim
r^{3/2}\epsilon^{-1/2}.$ (12)

This result also implies the universal scaling of the correlation function as

$r^{2}G(d)\sim e^{-d/\xi_{\mathrm{eq}}}.$ (13)

In contrast to the SMVC model, the SVL method provides access to the numerical
investigation of equilibrium properties. This is achieved through an
adiabatically slow time evolution, allowing the system to reach thermal
equilibrium at temperature $T$ at all time instances via noise-induced
fluctuations. To this end, the correlation function is computed by averaging
over $3\times 10^{4}$ trajectories for $\tau_{Q}=2560$. For both damping
regimes, the rescaled correlators exhibit an accurate scaling collapse as a
function of $d/\xi_{\mathrm{eq}}$, while the correlation length grows
approximately as $\xi_{\mathrm{eq}}\sim\lvert\epsilon\rvert^{-1/2}$ near the
critical point, indicating $\nu=1/2$. The results of this numerical analysis
are displayed in Fig. 2.

Finally, the effect of the finite-range interaction on the order parameter is
put to the test. In analogy with the $r=1$ case, the classical version of the
magnetization is introduced as the order parameter. This measure accurately
quantifies the distance of the final rotor configuration from the
ferromagnetic ground state,

$\mathcal{M}_{z}(t)=\frac{1}{N}\sum_{i=1}^{N}\left\langle|\sin\theta_{i}(t)|\right\rangle.$
(14)

The magnetization is determined for the same value of $\tau_{Q}=2560$ and
number of trajectories. Remarkably, the sudden increase of $\mathcal{M}_{z}$
around $t_{c}=0$ follows the same shape regardless of $r$ and $N$. The results
are displayed in Fig. 3 for the underdamped case in the regimes of small,
intermediate, and large values of $c$.

In short, the SVL model captures the equilibrium properties of classical
ferromagnetic $O(2)$ rotor network, indicating a continuous phase transition
in which the critical properties are robust against variations of the
interaction range, as long as the fundamental symmetries are not altered
compared to the chain case.

##  III  Universal Kibble-Zurek scaling regime and its breakdown

This section investigates how the excess energy density and the number of
near-ferromagnetic domains are related to the quench time $\tau_{Q}$ and the
connectance $c$.

###  III.1 The Kibble-Zurek mechanism

The KZ mechanism predicts a universal power-law scaling of the density of
excitations or topological defects generated while crossing a phase transition
point [10, 11, 12]. The critical slowing down prevents the system from
adjusting to the driving-induced changes, leading to the formation of
topological defects in the final state. To estimate the defect density, the KZ
mechanism utilizes equilibrium scaling theory, characterizing the power-law
divergence of the correlation length $\xi$ and the relaxation time $\tau$,

$\xi=\xi_{0}|\epsilon|^{-\nu}\quad\text{and}\quad\tau=\tau_{0}|\epsilon|^{-z\nu},$
(15)

in the proximity of the transition. Here, $\nu$ and $z$ denote respectively
the correlation-length critical exponent and the dynamic critical exponent,
with $\tau_{0}$ and $\xi_{0}$ fixing the model-dependent time and length
scales. Matching the driving rate and the relaxation time provides the freeze-
out time-scale $\hat{t}$. Inside the corresponding freeze-out regime,
$t\in[-\hat{t},\hat{t}]$ dynamics ceases to be adiabatic. This sets the
associated freeze-out length scale $\hat{\xi}$, determining the typical length
scale of the domains

$\begin{split}&\epsilon(\hat{t}\,)/\dot{\epsilon}(\hat{t}\,)=\tau_{0}\lvert\epsilon(\hat{t}\,)\rvert^{-z\nu}\Rightarrow\hat{t}=\tau^{\frac{1}{z\nu+1}}_{0}\tau^{\frac{z\nu}{z\nu+1}}_{Q}\,,\\\
&\hat{\xi}=\xi_{0}\lvert\epsilon(\hat{t}\,)\rvert^{-\nu}.\end{split}$ (16)

Correspondingly, defect-free domains will extend over regions of size
$\hat{\xi}^{d-D}$ leading to the average number of defects $\langle n\rangle$
scaling as

$\langle
n\rangle\propto\hat{\xi}^{-(d-D)}\propto\xi^{-(d-D)}_{0}\left(\frac{\tau_{0}}{\tau_{Q}}\right)^{\frac{(d-D)\nu}{z\nu+1}},$
(17)

valid for defects extending along $D$ dimensions, when the spatial dimension
of the system is $d>D$ [14].

###  III.2 Connectance dependence of the Kibble-Zurek scaling

As in long-range quantum Ising models [37], defects arising in circulant rotor
networks preserve their point-like or zero-dimensional character and are thus
$\mathbb{Z}_{2}$ kinks with $D=0$. Within the domains, the rotors are
approximately aligned ferromagnetically with $\theta_{i}\approx\pm\pi$ and
suddenly change their orientation at the location of a kink. This allows for a
one-dimensional representation of the defects, the density of which can be
estimated via the projected rotor values as

$n_{1}=\frac{\mathcal{N}_{1}}{N}=\frac{1}{2N}\sum_{i=1}^{N-1}[1-\text{sgn}(\sin\theta_{i})\text{sgn}(\sin\theta_{i+1})].$
(18)

According to the numerical results, this “one-dimensional” defect density
follows precisely the same power law as in the chain topology, sufficiently
below the fully connected limit $c\lesssim 0.15$,

$\displaystyle n_{1}\sim\tau^{-1/4}_{Q},\quad\text{ overdamped},$ (19)
$\displaystyle n_{1}\sim\tau^{-1/3}_{Q},\quad\text{ underdamped}.$ (20)

Even though the circulant graph topology does not modify the power law, it
plays a crucial role in the universal properties of the defect density beyond
rescaling the driving rate as $\sim r/\tau_{Q}$. Tentatively, substituting the
control parameter $\epsilon(t)=2r\,t/\tau_{Q}$ in the KZ rate equations (16)
yields the following freeze-out scales and the density of defects,

$\begin{split}&\hat{t}=\tau^{\frac{1}{z\nu+1}}_{0}\left(\frac{\tau_{Q}}{r}\right)^{\frac{z\nu}{z\nu+1}}\,,\\\
&\hat{\xi}=\xi_{0}|\epsilon\big{(}\hat{t}\big{)}|^{-\nu}=\xi_{0}\left(\frac{\tau_{Q}}{\tau_{0}r}\right)^{\frac{\nu}{z\nu+1}}\,,\\\
&n_{1}\propto\hat{\xi}^{-(d-D)}\propto\xi^{-(d-D)}_{0}\left(\frac{\tau_{0}r}{\tau_{Q}}\right)^{\frac{(d-D)\nu}{z\nu+1}}.\end{split}$
(21)

The identification of defects as kinks and the one-dimensional topology of the
resulting ferromagnetic domains imply that $d=1$ and $D=0$. The critical
exponents $\nu,\,z$ and the microscopic constants $\tau_{0},\xi_{0}$ may
generally depend on the regularity in the KZ power-law regime. As shown in Eq.
(9) and Eq. (10) and in App. A both $z\nu$ and $\nu$ are independent of $r$.
In particular, $z=1,\,2$ a for the overdamped and underdamped regimes,
respectively, while $\nu=1/2,\,1/2$. As for the model-dependent constants, the
regularity only enters via $\xi_{0}\sim r^{3/2}$, while the relaxation time is
independent of the graph topology, $\tau_{0}(r)\sim r^{0}$, see also Eqs. (9)
and (10).

Thus, the KZ exponents preserve their one-dimensional values,
$\alpha_{\mathrm{KZ}}=1/4$ and $\alpha_{\mathrm{KZ}}=1/3$ for the overdamped
and underdamped regimes, respectively. The KZ scalings, however, will not
follow a universal function of $c$ due to the explicit $r$ dependence of the
correlation length and the control parameter near the critical point. This
leads to the results

$\displaystyle\hat{\xi}\propto(r^{5}\tau_{Q})^{1/4}\Rightarrow
n_{1}\propto(r^{5}\tau_{Q})^{-1/4},\,\text{ overdamped},$ (22)
$\displaystyle\hat{\xi}\propto N(r^{\frac{7}{2}}\tau_{Q})^{1/3}\Rightarrow
n_{1}\propto(r^{\frac{7}{2}}\tau_{Q})^{-1/3},\,\text{ underdamped}.$ (23)

As presented in Fig. 4a for the overdamped case and in Fig. 9b in App. B for
the underdamped regime, the curves of $n_{1}$ follow the analytical
approximations within numerical precision and exhibiting a precise scaling
collapse for system sizes, $N=100,200,400$.

In addition, the scaling of the freeze-out time scale also acquires a
regularity dependence,

$\hat{t}\propto\begin{cases}&\left(\frac{\tau_{Q}}{r}\right)^{1/2},\quad\text{overdamped},\\\
&\left(\frac{\tau_{Q}}{r}\right)^{1/3},\quad\text{underdamped}.\end{cases}$
(24)

###  III.3  Universal energy density beyond nearest-neighbors

Next, we investigate the relationship between $n_{1}$ and excess energy
production. The circulant graph topology disrupts the simple one-to-one
correspondence between the number of near-ferromagnetic domains and the excess
energy. All rotor configurations within the interaction range contribute to
the latter, which is defined as follows

$\begin{split}&\Delta
E(t)=H(\bm{\theta},t)-E_{\mathrm{min}}(\bm{\theta},t),\\\
&\rho_{E}=\frac{\Delta
E(\tau_{Q})}{Nr}=1-\frac{1}{4}\sum_{i,j=1}^{N}A_{ij}\sin\theta_{i}\sin\theta_{j}.\end{split}$
(25)

Figure 4: Averages of the defect and excess-energy densities as functions of
$\tau_{Q}$ in the overdamped regime with the connectance varying from the
chain topology towards the fully connected limit for system sizes
$N=100,\,200,\,400$. Panels $(a)-(c)$: One-dimensional defect density. For
values $c\lesssim 0.15$, the same power law of the chain topology is observed
with a universal regularity rescaling. For intermediate and densely connected
graphs, a crossover is displayed towards the adiabatic regime, exhibiting a
universal exponential shape. Panels $(d)-(f)$: The graph defect density
$n_{2}$ behaves similarly, exhibiting power laws of $r$ different from the KZ
regime but identical exponential behavior in the adiabatic regimes. Panels
$(g)-(i)$: The excess energy density closely matches the behavior of graph
defect density with high precision, including additional features in the
absence of defects that stem from the dynamics of small angle deviations
around the ferromagnetic direction. In all regimes, the curves were averaged
over $10^{3}$ trajectories.

Here, the first equation accounts for the excess energy during the evolution,
while the second is for its density at the end of the schedule. The
normalization in the second row makes the density scale invariant with respect
to both $N$ and $r$. Here, the minimum of the instantaneous ground state
energy $E_{\mathrm{min}}(\bm{\theta},t)=\min_{\\{\theta_{i}\\}_{i\in
V}}H(\bm{\theta},t)$ is defined as the global minimum of the Hamiltonian at
time $t$. By translational invariance, all angles are identical, and the
search for the instantaneous energy minimum simplifies to the minimization of
a single-variable function. For convenience, the schedule is considered in the
interval $[0,\tau_{Q}]$ with the critical time given by $t_{c}=\tau_{Q}/2$,

$\begin{split}\partial_{\theta}H(\theta,t)&=rN\partial_{\theta}\left[-\frac{t}{\tau_{Q}}\sin^{2}\theta-2\left(1-\frac{t}{\tau_{Q}}\right)\cos\theta\right]=0\\\
&\Rightarrow\frac{t}{\tau_{Q}}\sin\theta\cos\theta=\left(1-\frac{t}{\tau_{Q}}\right)\sin\theta.\end{split}$
(26)

For $t\leq\tau_{Q}/2$ this equation has the trivial solution $\theta=0$ and
with minimum energy
$E_{\mathrm{min}}(\bm{\theta},t)=-2rN\left(1-t/\tau_{Q}\right)$, while for
$t>\tau_{Q}/2$ the solution is given by $\cos\theta=\frac{\tau_{Q}-t}{t}$.
Thus, the minimum energy yields

$\begin{split}E_{\mathrm{min}}(\bm{\theta},t)=-rN\biggl{[}\frac{t}{\tau_{Q}}\biggl{[}1&-\left(\frac{\tau_{Q}-t}{t}\right)^{2}~{}\biggr{]}\\\
&+2\frac{\tau_{Q}}{t}\left(1-\frac{t}{\tau_{Q}}\right)^{2}~{}\biggr{]}.\end{split}$
(27)

This expression provides a convenient way for the numerical computation of Eq.
(25). The excess energy is set by the actual rotor configuration with
contributions attributed to the kinks as well as angle deviations within the
domains. However, these angle fluctuations provide negligibly small
contributions compared to the energy scale of the kinks, and $\rho_{E}$ to
leading order can be expressed in terms of $n_{1}$. For values of $\tau_{Q}$
where the freeze-out length sufficiently exceeds the interaction range
($\hat{\xi}>r$), every defect carries an energy contribution of $\sim r^{2}$.
This originates from the antiparallel rotors on opposite sides of the defect.
To this end, $n_{1}$ is also generalized by counting all rotor-rotor
fluctuations within the interaction range, introducing the graph defect
density

$n_{2}=\frac{\mathcal{N}_{2}}{Nr}=\frac{1}{4Nr}\sum_{i,j=1}^{N}A_{ij}[1-\text{sgn}(\sin\theta_{i})\text{sgn}(\sin\theta_{j})].$
(28)

Within the domain of the KZ scaling, the graph defect density agrees well with
the numerically exact excess energy density. Due to the normalization with the
connectivity, $2rN$, both $\rho_{E}$ and $n_{2}$ are also proportional to
$n_{1}$ by a factor of $r$. Remarkably, $\rho_{E}$ and $n_{2}$ follow the same
universal function of $\tau_{Q}/\hat{t}(r)$ independently of both $r,\,N$ and
the respective damping regimes up to the leading order,

$\begin{split}&n_{2}\propto
rn_{1}\propto\left(\frac{\tau_{Q}}{\hat{t}(r)}\right)^{-1/2},\\\
&\rho_{E}=\frac{\Delta E}{Nr}\propto
n_{2}\propto\begin{cases}&(r\,\tau_{Q})^{-1/4},\quad\text{overdamped},\\\
&(\sqrt{r}\,\tau_{Q})^{-1/3},\,\text{underdamped},\end{cases}\end{split}$ (29)

as displayed in Fig. 4d,g for the overdamped case for different values of $N$
and $c$ (see also Fig. 9d,g in App. B for the underdamped case).

##  IV Breakdown of Kibble-Zurek mechanism and the fully connected limit

In this section, the breakdown of the universal KZ scaling regime is
investigated in terms of the regularity $r$, connectance $c$, and quench time
$\tau_{Q}$. Subsequently, the universal signatures of annealing dynamics near
the fully connected limit are explored. Finally, universal properties of
$\rho_{E}(t)$ are analyzed at intermediate times,
$t\in[-\tau_{Q}/2,\tau_{Q}/2]$ in both the KZ and the adiabatic regimes.

###  IV.1 Fast quenches

As in the TFIM and the $\phi^{4}$ classical stochastic model [73, 74], the
universal power-law scaling of both $n_{1}$ and $\rho_{E}$ breaks down for
rapid annealing processes. For small $\tau_{Q}$, the condition $\hat{\xi}<r$
coincides with the fast quench breakdown of the KZ scaling. In particular, the
$\tau_{Q}/r$ rescaling in the control parameter allows the system to stay
inside the freeze-out regime for the whole process. The fast quench breakdown
can be captured by matching the freeze-out time with the equilibrium
relaxation time, $\hat{t}=\tau(\tau_{Q})$ [73] at the end of the process. It
sets the condition for $\tau_{Q}$ and the lower threshold of the freeze-out
length scale

$\begin{split}&\left(\frac{\tau_{Q}}{r}\right)^{\frac{z\nu}{z\nu+1}}=\tau_{0}r^{-z\nu}\Rightarrow\tau^{\mathrm{fast}}_{Q}\propto
r^{-z\nu},\\\ &\Rightarrow\hat{\xi}\propto
r^{3/2}\left(\frac{\tau_{Q}}{r}\right)^{\frac{\nu}{z\nu+1}}\propto
r^{3/2-\nu}=r,\end{split}$ (30)

using $\nu(r)=1/2$ in the last step. Thus, the fast quench analysis via
matching the relaxation and freeze-out time scales provides the same result as
using the condition of $\hat{\xi}=r$,

$\tau^{\mathrm{fast}}_{Q}\sim
r^{-z\nu}\sim\begin{cases}&r^{-1},\,\,\,\,\,\quad\text{overdamped},\\\
&r^{-1/2},\quad\text{underdamped}.\end{cases}$ (31)

Furthermore, this approach is validated by the fact that by increasing either
$r$ or $\tau_{Q}$, the freeze-out time scale remains inside the domain of
validity of the linear approximation, $\hat{t}/\tau_{Q}\ll 1$. The fast quench
breakdown is exhibited in Fig. 4a,d,g for $n_{1},\,n_{2}$ and $\rho_{E}$,
respectively. As $r$ is increased, these measures terminate earlier at
constant plateaus attributed to the fast quench regime (see also App. B for
the underdamped case). However, for larger values of $c$, these plateaus
disappear in the limit $N\rightarrow\infty$.

###  IV.2 Fully connected limit and adiabatic regime

As discussed in Sec. III.3, considering the similarities between the one-
dimensional rotor chain and the TFIM model in both damping regimes, one could
expect the complete rotor network with graph $K_{N}$ to parallel the behavior
of the Lipkin-Meshkov-Glick (LMG) model [42]. However, the results clearly
show that this is not the case.

The KZ regime only persists until the onset of adiabaticity identified by the
condition $n_{1}=\frac{1}{N}$. This sets the upper threshold of of $\tau_{Q}$
in terms of $r$ and $N$

$\displaystyle
r^{-3/2}\left(\frac{\tau_{Q}}{r}\right)^{-\frac{\nu}{z\nu+1}}\propto\frac{1}{N}$
(32) $\displaystyle\Rightarrow\tau^{\mathrm{ad}}_{Q}\propto
N^{4}r^{-5}\quad\text{and}\quad\propto N^{3}r^{-7/2},$ (33)

for the overdamped and the underdamped regimes, respectively. Note that this
condition is equivalent to matching the freeze-out length scale with the
system size, $\hat{\xi}=N$.

Thus, the universal KZ power-law regime completely disappears around the
threshold value of the connectance for which the fast quench time scale (30)
and the adiabatic one in (32) become of the same order in the limit
$N,r\rightarrow\infty$. This leads to the breakdown scale of

$\tau^{\mathrm{fast}}_{Q}=\tau^{\mathrm{ad}}_{Q}\Rightarrow r\sim N\rightarrow
c\sim\mathrm{O}(1),$ (34)

for both the underdamped and overdamped regimes. The breakdown scales suggest
that the larger $r$ is, the earlier the KZ regime sets in. However, as $N$
increases, the adiabatic regime takes over the dynamics even faster,
eventually washing out the power-law KZ scaling. Remarkably, the corresponding
threshold connectance remains finite in the thermodynamic limit, with
$c\approx 0.15$ for both damping regimes. When the connectance
$c\in[0.15,0.5]$, the KZ regime disappears, and the defect densities exhibit
finite fast quench plateaus followed by an intermediate non-universal sharply
decaying regime. The corresponding shrinking and vanishing of the universal KZ
regime are displayed in Fig. 4b,c for $n_{1}$ and Fig. 4e,f,h,i for $n_{2}$
and $\rho_{E}$ in the overdmaped regime (see also Fig. 9 in App. B for the
underdamped regimes).

Beyond this scale, $c\gtrsim 0.5$, the defect and energy densities turn into a
sharp, exponential decay. This regime exhibits new distinct universal behavior
taking over the dynamics, with the fast quench plateaus disappearing with
increasing $N$. In both cases, the universal regime also becomes independent
of $N$.

The exponential decay of $n_{1}$ and $n_{2}$ is shown in Fig. 4c,f,i (see also
Fig. 9c,f,i in App. B), for quench times, where the last effect pair can still
be detected numerically for the used number of stochastic trajectories. These
curves become universal as a function of the rescaled annealing time
$r^{\frac{z\nu}{z\nu+1}}\tau_{Q}$,

$\begin{split}&\log n_{1}\sim\log
n_{2}\sim-0.25\,r^{1/2}\tau_{Q},\quad\text{overdamped},\\\ &\log n_{1}\sim\log
n_{2}\sim-0.7\,r^{1/3}\tau_{Q},\quad\text{underdamped}.\end{split}$ (35)

Although the decaying part exhibits non-universal behavior in the intermediate
regime $c\in[0.15,0.5]$, the best scaling collapse for moderate values of
$\tau_{Q}\lesssim 10$ is obtained by the same rescaling as in Eq. (35). This
is shown in Fig. 4b,d in the overdamped regime, for $n_{1}$ and $n_{2}$,
respectively (see also App. B for the underdamped case).

By contrast, the excess energy can take finite values also for much larger
values of $\tau_{Q}$ due to small angle deviations compared to the
$\theta=\pi/2$ ferromagnetic direction. Within this defect-free regime,
$\rho_{E}$ is governed by these small angle fluctuations. Its behavior thus
differs from that of $n_{1}$ and $n_{2}$. As shown in Fig. 4i for the
overdamped regime (see also App. B for the underdamped case), the exponential
decay of the excess energy becomes independent of $r$ up to exponential
accuracy and follows the curve given by the best numerical fit,

$\displaystyle\log\rho_{E}\sim\,0.8\tau_{Q},\quad\text{overdamped},$ (36)
$\displaystyle\log\rho_{E}\sim\,0.06\tau_{Q},\,\,\,\text{underdamped}.$

As for the intermediate regime $c\in[0.15,0.5]$, the best scaling collapse for
moderate values of $\tau_{Q}$ is achieved by a slightly different rescaling of
$\tau_{Q}\rightarrow\sqrt{r}\,\tau_{Q}$, as observed in Fig. 4h for the
overdamped regime. Quite remarkably, in the underdamped regime, a perfect
scaling collapse is observed independently of $r$, as shown in App. B.

###  IV.3 Universal dynamics and scaling collapse

So far, the defect and excess energy densities have been investigated at the
end of the time evolution as a function of $r$ and $\tau_{Q}$. To provide a
more complete picture of the emerging universal properties, dynamical
characteristics at intermediate times during the quench are also put to the
test.

Figure 5: Time evolution of the excess energy density near the critical point.
The universal character of the dynamics is demonstrated for different values
of $r,\,\tau_{Q}$ and $N$ by displaying $\rho_{E}$ as a function of
$t/\hat{t}$ in the KZ regime both for $(a)$ the overdamped and $(b)$ the
underdamped regime. Blue symbols were obtained with $N=201$ while the red
curves with $N=401$. In both cases, $10^{4}$ trajectories were employed.

The time-evolution of $\rho_{E}(t)$ near the critical point is in good
agreement with the predictions of the KZ mechanism. The growth of
$\rho_{E}(t)$ becomes independent of $r,\,\tau_{Q}$ and $N$ for connectances
inside the KZ regime following a single universal curve as a function of
$t/\hat{t}$ up to numerical precision. This is demonstrated in Fig. 5a,b for
the overdamped and underdamped regimes, respectively.

Moreover, additional universal signatures emerge in both the KZ and adiabatic
regimes with $r$ independent time-evolution along the whole process. As shown
in Fig. 6a,b, the curves of $\rho_{E}(t)$ exhibit precise scaling collapse for
different values of $N$. This universality only holds inside the KZ and
adiabatic regimes, as also shown by the deviations of the curves with
$N=200,\,r=12$ reaching the boundary of the KZ regime and $N=400,\,r=25$ near
the onset of adiabaticity. However, for fixed $\tau_{Q}$, the time-evolution
exhibits non-universal behavior with respect to $r$ even when plotted for the
rescaled time, $t/\tau_{Q}$. As displayed in 6a,b by decreasing either $r$ or
$\tau_{Q}$, $\rho_{E}(t)$ exhibits stronger non-adiabatic effects with larger
final values and more pronounced peaks around the critical point. In this
case, the final values follow either KZ power laws or decay exponentially in
the adiabatic regime, depending on $r$ and $\tau_{Q}$.

We note that varying the annealing time while keeping fixed the ratio
$\tau_{Q}/\hat{t}(r)$ in the KZ regime, or the product of
$r^{\frac{z\nu}{1+z\nu}}\tau_{Q}$ in the adiabatic regime, unveils a precise
scaling collapse of the date for different values of $N$ and $r$. This is
demonstrated in Fig. 6c in the KZ regime and in Fig. 6d in the adiabatic case.
As with the dependence of $n_{1}$ and $\rho_{E}$ on $\tau_{Q}$, the
nonequilibrium universal dynamics is brought out upon separate rescaling by
$N$ and $r$, rather than by only $c$. The proper rescaling can take different
forms depending on the investigated quantities, the damping regimes, and the
power-law or the exponential behavior.

Figure 6:  Time-evolution of the energy density $\rho_{E}(t)$ as a function of
$t/\tau_{Q}$ for the overdamped regime for different values of $N,\,\tau_{Q}$
and $r$. Panels $(a),\,(b)$: Fixed quench times for pairs of $r$ and $N$
exhibiting size-independent behavior but becoming non-universal with respect
to $r$ both in $(a)$ the KZ and $(b)$ the adiabatic regimes. Panel $(c)$:
Varying $\tau_{Q}$ with fixed $r\tau_{Q}=5120$. Good scaling collapse is found
in the KZ regime for different values of $N$. Panel $(d)$: Adiabatic regime of
$c\gtrsim 0.5$ and $\tau_{Q}\gg\hat{t}$ featuring similar scaling collapse in
the adiabatic regime. In all panels, averages involve $10^{4}$ trajectories,
with blue symbols for $N=201$, and red lines for $N=401$. Figure 7:  Cumulant
ratios as a function of $c$ for $n_{1},\,n_{2}$ and $\rho_{E}$, for fixed
values of $\tau_{Q}$ and for $N=100,\,200,\,400$. Panels $(a)-(c)$: The ratio
$\kappa_{2}/\kappa_{1}$ behaves similarly and exhibits the predicted universal
linear growth for $\rho_{E}$ and $n_{2}$, while it remains constant for
$n_{1}$. For relatively small values of $c$, a crossover is observed between
the sub- and super Poissonian behavior within the KZ regime, while the same
happens for larger values of $c$ in the case of the $n_{1}$ around the
breakdown of the power-law regime. These cumulant ratios are estimated from
$5\times 10^{4}$ trajectories. Panels $(d)-(f)$: The ratio
$\kappa_{3}/\kappa_{1}$ exhibits a universal quadratic growth with $c$ inside
the KZ regime for $\rho_{E}$ and $n_{2}$, while it follows the predicted
constant value for $n_{1}$. Cumulants are extracted from $10^{4}$
trajectories.

##  V Beyond the Kibble-Zurek mechanism: fluctuations in the defect and energy
densities

In this section, we investigate the second and third cumulants of the defect
densities and excess energy. Beyond the KZ scaling of the average defect
densities, higher-order cumulants are of utmost importance regarding the full
energy distribution. In the quantum case, it was shown for the TFIM that the
probability distribution to find $n$ kinks in the system at the end of the
quench, $P(n)=\mathbb{E}[\delta(\mathcal{N}-n)]$, follows a Poisson binomial
distribution [75, 76], as demonstrated in annealing devices [7, 27, 8]. In
fully connected quantum systems, a negative binomial distribution of
fractional index was found instead [45]. In the classical domain, a binomial
distribution was predicted for point-like defects [75, 77, 31, 20, 78, 79],
with a Poisson binomial distribution occurring in the more general setting
[77, 80]. Fourier transforming $P(n)$ yields the characteristic function
$\hat{P}(\phi)=\mathbb{E}[e^{i\phi n}]$, whose logarithm gives the cumulant
generating function $\ln{\mathbb{E}[e^{i\phi
n}]}=\sum_{p=1}^{\infty}\frac{\kappa_{p}}{p!}(i\phi)^{p}$, where $\kappa_{p}$
is the $p$-th order cumulant. It has been argued that all cumulants share the
same power-law scaling exhibited by the mean, $\kappa_{q}\propto\kappa_{1}$
for every $q$. As a result, the ratio $\kappa_{q}/\kappa_{1}$ is constant and
independent of $\tau_{Q}$ [75, 77]. For instance, for the TFIM
$\kappa_{2}/\kappa_{1}=2-\sqrt{2}\approx 0.578$ [75, 17], as verified in
D-Wave devices [7]. As for classical systems, the same binomial relation of
the first three cumulants has been shown for the chain topology [20]. In what
follows, we focus on the second and third cumulants, $\kappa_{2}$ and
$\kappa_{3}$, beyond the average density ($\kappa_{1}$).

Figure 8: Second cumulants of the excess energy density as a function of
$\tau_{Q}$ for system size $N=100,\,200,\,400$ and connectances, $c\in[0,1]$.
Panels $(a)-(c)$: cumulant ratio $\kappa_{2}(\rho_{E})/\rho_{E}$ in the
overdamped and underdamped cases, for different values of $c$ in the KZ,
intermediate and adiabatic regimes. In the KZ regime, the cumulant ratio is
independent of $\tau_{Q}$ and increases universally with $c^{2}$. In the
adiabatic regime, the ratio becomes proportional to the average decaying
exponentially, as also highlighted with the additional fittings. Panels
$(d)-(f)$: Second cumulant in the underdamped regime as a function of
$\tau_{Q}$ in the same three regimes of $c$. Inside the KZ regime, the same
power law is displayed as for the averages, but with a different universal
rescaling with $r$. In the intermediate regime $c\in[0.15,0.5]$, a perfect
scaling collapse is found as a function of $\tau_{Q}$. The adiabatic limit
exhibits an exponential decay, matching approximately the square of the
average. Panels $(g)-(i)$: Similar features with the corresponding power-laws
in the overdamped regime with a reasonably faster adiabatic decay and less
universal behavior for intermediate connectance values. An ensemble of
$10^{4}$ trajectories is used.

Their numerical estimate is performed with the same number of $5\times 10^{4}$
trajectories. As shown in Fig. 7c,f the $\kappa_{3}/\kappa_{1}$ and
$\kappa_{2}/\kappa_{1}$ cumulant ratios of the one-dimensional defects follow
precisely the same value as that of the $r=1$ case and of the TFIM. This
agreement breaks down at slightly smaller values of the connectance, $c\approx
0.1$, than the boundary of the universal KZ scaling of the averages.
Correspondingly, $\kappa_{2}(n_{1})$ and $\kappa_{3}(n_{1})$ also become
universal functions of the $r\tau_{Q}$ and $\sqrt{r}\tau_{Q}$ for the
overdamped and underdamped regime, respectively (for further details see App.
C and App. D).

Remarkably, in the fast quench regime, the proportionality is maintained with
the same constant value as in the power law regime, as proposed in Ref. [80].
This indicates that the shape of the distribution is the same under fast
quenches and in the KZ regime. However, this is not true at intermediate
quenches as signaled by the non-universal peaks in a small interval of
$\tau_{Q}$ before entering the KZ regime (see Apps. C and D). As demonstrated
in Apps. C and D, in the adiabatic regime, the cumulants of $n_{1}$ and
$n_{2}$ follow the same exponential law, given by the best numerical fit,

$\begin{split}&\log\kappa_{2}(n_{1})\sim\,-0.25r^{1/2}\tau_{Q},\,\log\kappa_{3}(n_{1})\sim\,-0.75r^{1/2}\tau_{Q},\\\
&\log\kappa_{2}(n_{1})\sim\,1.4r^{1/3}\tau_{Q},\,\log\kappa_{3}(n_{1})\sim\,2.1r^{1/3}\tau_{Q}.\end{split}$
(37)

This implies that up to numerical precision, the proportional nature of the
cumulants turns into a power law relation up to exponential accuracy,

$\kappa_{2}(n_{1})\sim n^{2}_{1},\quad\kappa_{3}(n_{1})\sim n^{3}_{1}.$ (38)

This is in good agreement with the discussion of Sec. IV.2 on the universal
signatures in the intermediate regime. The exponential fits in Eq. (37)
further demonstrate the universal rescaling
$r^{\frac{z\nu}{1+z\nu}}\,\tau_{Q}$ in the adiabatic limit.

In the case of $\rho_{E}$, the cumulant ratios behave similarly to those of
$n_{1}$, as observed in Fig. 8a and in App. D. In the KZ regime, the ratio
$\kappa_{2}(\rho_{E})/\kappa_{1}(\rho_{E})$ approximately follows constant
plateaus, which get shifted as $r$ is increased. This happens because, as for
the averages, the excess energy depends on the actual rotor configuration and
includes contributions from all rotors, including those at the sides of each
kink. Due to the graph theoretical normalization $1/Nr$, this effect boosts
$\kappa_{2}(\rho_{E})$ and $\kappa_{3}(\rho_{E})$ with $r^{2}$ and $r^{3}$,
respectively. This yields the following universal power laws,

$\displaystyle\kappa_{2}(\rho_{E})\sim\left(r^{-3}\tau_{Q}\right)^{-1/4},\,\kappa_{3}(\rho_{E})\sim\left(r^{-7}\tau_{Q}\right)^{-1/4},$
(39)
$\displaystyle\kappa_{2}(\rho_{E})\sim\left(r^{-5/2}\tau_{Q}\right)^{-1/3},\,\kappa_{3}(\rho_{E})\sim\left(r^{-11/2}\tau_{Q}\right)^{-1/3}.$
(40)

As shown in Fig. 7a,b and Fig. 7d,e, the corresponding $\kappa_{3}/\kappa_{1}$
and $\kappa_{2}/\kappa_{1}$ cumulant ratios acquire scales $\sim r$ and $\sim
r^{2}$, respectively leading to

$\displaystyle\frac{\kappa_{2}(n_{2})}{n_{2}}\sim\frac{\kappa_{2}(\rho_{E})}{\rho_{E}}\sim\frac{\frac{n_{1}}{Nr^{2}}r^{4}}{rn_{1}}\propto
c,$ (41)
$\displaystyle\frac{\kappa_{3}(n_{2})}{n_{2}}\sim\frac{\kappa_{3}(\rho_{E})}{\rho_{E}}\sim\frac{\frac{n_{1}}{N^{2}r^{3}}r^{6}}{rn_{1}}\propto
c^{2},$ (42)

where the proportionality of the $n_{1}$ cumulants was explicitly written out.
Note that the ratios exhibit precise scaling collapse as a function of $c$ in
contrast to the universal $r$ dependence of the averages. Thus, energy
cumulants naturally exhibit a crossover from the sub- to the super-Poissonan
regime as $c$ is increased. This also implies that the corresponding
distributions are no longer concentrated around the average and fluctuate
asymmetrically.

In contrast to $n_{1}$, for fast quenches, the proportionality breaks down for
the cumulant ratios of both $n_{2}$ and $\rho_{E}$, as displayed for
$\kappa_{2}(\rho_{E})$ in Fig. 8a,b,c for the small $c\lesssim 0.15$,
intermediate $c\in[0.15,0.5]$, and large $c\gtrsim 0.5$ connectance regimes,
respectively. In the adiabatic regimes, the same exponential dependence is
found as for the averages, as demonstrated in Fig. 8b,c for
$\kappa_{2}(\rho_{E})$, also highlighted by additional exponential fittings
(see also App. D for $\kappa_{3}(\rho_{E})$). However, the best numerical fits
yield the following $r$ independent, universal exponential behavior,

$\displaystyle\log\kappa_{2}(\rho_{E})\sim\,1.6\tau_{Q},\quad\log\kappa_{3}(\rho_{E})\sim\,2.4\tau_{Q},$
(43)
$\displaystyle\log\kappa_{2}(\rho_{E})\sim\,0.12\tau_{Q},\quad\log\kappa_{3}(\rho_{E})\sim\,0.18\tau_{Q},$
(44)

up to exponential accuracy in the overdamped and underdamped, respectively.
Note that the $r$ independent behavior of the decay is dominated by small
angel fluctuations rather than by the occurrence probability of the last pair
of defects.

By contrast, the cumulants of $n_{2}$ exhibit the same characteristics as
those of $n_{1}$ up to exponential accuracy, as in this regime, defect
densities are dominated by the vanishing occurrence probability of a single
pair of defects. Accordingly, the proportionality between the cumulants of
$\rho_{E}$ and $n_{2}$ follow similar power law relations as in the case of
$n_{1}$,

$\displaystyle\kappa_{2}(\rho_{E})\sim\rho^{2}_{E}~{},$
$\displaystyle\kappa_{3}(\rho_{E})\sim\rho^{3}_{E},$ (45)
$\displaystyle\kappa_{2}(n_{2})\sim n^{2}_{2}~{},$
$\displaystyle\kappa_{3}(n_{2})\sim n^{3}_{2}.$

Notably, keeping $r\sim N,\,c\lesssim 0.1$ enforces the same power law
relations as in Eq. (45).

Finally, let us also provide details about the surviving universal features in
the intermediate regime. As shown in Fig. 8a, a perfect scaling collapse is
found as a function of $\tau_{Q}$ in the underdamped regime. In the overdamped
regime, universal signatures only survive up to not too large annealing times,
$\tau_{Q}\lesssim 10$, with the rescaling of $\sqrt{r}\,\tau_{Q}$, as shown in
Fig. 8e.

##  VI Conclusions

We have investigated the annealing dynamics of finite-range interacting
$\mathrm{O}(2)$ rotor networks on circulant regular graphs using numerical
simulations and approximate analytical methods. This model provides an ideal
testbed to benchmark quantum annealing devices and test the Kibble-Zurek
mechanism in finite-range Ising models interacting with a thermal bath. Using
the SVL model, both the full-time evolution and equilibrium properties can be
accessed. The latter is achieved by performing adiabatically slow time
evolutions. In particular, we showed that for regular networks, the critical
exponents and the average magnetization near the critical point remained
invariant with respect to variations of the connectance. The only dependence
was involved in the correlation length, which varied by a multiplicative
factor.

Numerical results revealed that the KZ mechanism accurately described the
universal power-law scaling of the density of point-like defects separating
one-dimensional near-ferromagnetic domains. However, this measure could not
account for the excess energy generated during the process. Therefore, the
generalized graph defect density was introduced, which matched the excess
energy density with high precision. For these quantities, the KZ power laws
provided an accurate description in terms of the equilibrium critical
exponents, however, with different regularity dependencies compared to the
one-dimensional defects.

The robustness of the universal KZ regime was explored as a function of the
annealing time and the connectance. The defect and excess energy densities
exhibited a transition to the universal fast quench regime below the
corresponding regularity-dependent time scale. At the other extreme, adiabatic
evolution dominated the dynamics above a system size- and regularity-dependent
time scale. These threshold time scales identify a finite connectance for the
breakdown of the KZ regime. Above the adiabatic connectance scale, the defect
and excess energy densities decayed exponentially. For both the one-
dimensional and graph defects, similar universal properties were observed,
albeit with regularity dependencies distinct from those in the KZ regime. In
the case of the excess energy, the adiabatic regime was dominated by small-
angle deviations around the ground state. This resulted in a regularity-
independent universal exponential decay in defect-free regimes of $\tau_{Q}$.

Finally, we analyzed the statistical distribution of the kink and energy
densities by exploring the second and third cumulants. For the one-dimensional
defect density, constant cumulant ratios were found in the KZ regime, with the
same ratios as in the TFIM. Remarkably, this behavior persisted in the fast-
quench regime as well. In the adiabatic regime, the ratios followed a power-
law relation, which matched that of the graph defect density up to exponential
accuracy. In the case of the graph defect and excess energy densities, the
proportionality between the cumulants acquired additional linear and quadratic
regularity dependencies for the second and third cumulants, respectively. In
contrast to the one-dimensional defects, this behavior broke down for fast
quenches, while in the adiabatic regime, the cumulants followed the same
power-law relation. Specifically, in this regime, the graph defect cumulants
matched those of the one-dimensional defects up to exponential accuracy.
However, in the case of the excess energy, different regularity-independent
universal behavior was found. Similar to the averages, this behavior arose
from the dynamics of small-angle deviations around the ferromagnetic
directions in defect-free regimes of $\tau_{Q}$.

In short, we have established the interplay of universality and its breakdown
in the annealing dynamics of classical rotor networks with finite-range
interactions. Our findings not only contribute to the understanding of
nonequilibrium statistical mechanics of complex networks but also offer a
benchmark for the performance of quantum simulators and annealing devices. As
a step in this direction, during the writing of this manuscript, the preprint
[81] reported the use of SVL on cyclic and random graphs as a benchmark in
complex optimization by the Vector Ising Spin Annealer.

###### Acknowledgements.

This research was funded by the Luxembourg National Research Fund (FNR), grant
reference 17132054. For the purpose of open access, the authors have applied a
Creative Commons Attribution 4.0 International (CC BY 4.0) license to any
Author Accepted Manuscript version arising from this submission.

##  Appendix A Connectance dependence of the correlation length and time near
the critical point

In this appendix, we present an approximate analytical derivation of the
spatial correlation function, unveiling its dependence on the connectance.
Near the critical point, in the small angle limit, the Hamiltonian in Eq. (4)
can be expanded as

$H(t)\approx-\frac{J(t)}{2}\sum_{i,j=1}^{N}A_{ij}\theta_{i}\theta_{j}-h(t)\sum_{i=1}^{N}\theta^{2}_{i}.$
(46)

Figure 9: Averages of defect and excess energy densities as a function of
$\tau_{Q}$ in the overdamped regime with the connectance varying from the
chain topology towards the fully connected limit for system sizes
$N=100,\,200,\,400$. Panels $(a)-(c)$: One-dimensional defect density for
values $c\lesssim 0.15$. The same power law is observed as in the chain
topology with a universal regularity rescaling. For the intermediate and large
connectances, a crossover is displayed towards the adiabatic regime,
exhibiting a universal exponential shape. Panels $(d)-(f)$: Similar
observations for $n_{2}$ exhibiting the predicted universal KZ power laws of
$r$ for $c\lesssim 0.15$. Different universal rescaling emerges in the
intermediate and adiabatic regimes. In the former, universal behavior breaks
down for $\tau_{Q}\gtrsim 10$. Panels $(g)-(i)$: Excess energy density,
matching up to high precision the behavior of the finite range defect
densities. In contrast to $n_{2}$, for $c\gtrsim 0.15$, a perfect scaling
collapse is observed as a function of $\tau_{Q}$. In all regimes, curves
involve averages over $10^{3}$ trajectories.

The correlation function is captured by directly considering the angle-angle
equal-time two-point function, which is also expanded to leading order as

$\begin{split}G(d,t)&=\overline{\langle\sin\theta_{n}(t)\sin\theta_{n+d}(t)\rangle_{\mathrm{traj}}}\\\
&\approx\overline{\langle\theta_{n}(t)\theta_{n+d}(t)\rangle_{\mathrm{traj}}}\\\
&=\frac{1}{N}\sum_{n=1}^{N}\langle\theta_{n}(t)\theta_{n+d}(t)\rangle_{\mathrm{traj}},\end{split}$
(47)

with $\overline{\cdots}$ denoting the spatial average over the rotors
exploiting the spatial translational invariance of the system. The average
over the random stochastic trajectories is replaced by the equilibrium
canonical averaging with the Boltzmann weight $e^{-\beta H(t)}$ as dictated by
the fluctuation-dissipation theorem. Here, the inverse temperature
$\beta=1/T=1000/J$ is given by the applied value of $T=0.001$ and where
$k_{B}=1$ was set for convenience. To extract the equilibrium correlation
length, the angles are expanded in Fourier series as

$\theta_{n}=\frac{1}{N}\sum_{k}\,e^{-2\pi ikn}\theta_{k},\quad
k=\frac{2\pi}{N},\,\frac{4\pi}{N},\,\dots,\,2\pi.$ (48)

This leads to the following Hamiltonian

$\begin{split}H(t)&=-\frac{1}{2}J(t)\sum_{k,k\prime}\theta_{k}\theta_{k^{\prime}}\,\sum_{n=1}^{N}\,e^{2i\pi(k+k^{\prime})n}\sum_{m=-r,\,m\neq
0}^{r}\,e^{-2i\pi
km}+h(t)\frac{1}{N^{2}}\sum_{n=1}^{N}\,e^{-2i\pi(k+k^{\prime})n}\sum_{k,k^{\prime}}\theta_{k}\theta_{k^{\prime}}\\\
&=-NJ(t)\sum_{k}\,\theta_{k}\theta_{-k}\,\sum_{m=1}^{r}\cos\left(2\pi
kn\right)+Nh(t)\sum_{k}\,\theta_{k}\theta_{-k}=\sum_{k}\,A_{k}\lvert\theta_{k}\rvert^{2},\end{split}$
(49)

where the identity $\frac{1}{N}\sum_{n=1}^{N}\,e^{-2i\pi
kn}=\delta_{k,-k^{\prime}}$ was used. Furthermore, as the angles are real
variables, the products were rewritten as
$\theta_{k}\theta_{-k}=\lvert\theta_{k}\rvert^{2}$. Finally, we also
introduced the overall constant of the quadratic Hamiltonian as
$A_{k}=\frac{1}{N}\left[-2J\sum_{m=1}^{r}\cos\left(2\pi km\right)+h\right]$.
Here, for the sake of convenience, the time arguments were dropped everywhere.

Next, the Fourier decomposition is also applied to the correlation function

$G(d)\approx\left\langle\frac{1}{N}\sum_{n=1}^{N}\theta_{n}\theta_{n+d}\right\rangle=\frac{1}{N^{2}}\sum_{k}\,e^{-2i\pi
dk}\langle\lvert\theta_{k}\rvert^{2}\rangle,$ (50)

where the same steps were performed as above in Eq. (49). The single $k$
averages are simply evaluated by Gaussian integrations with variance $\beta
A_{k}$

$\langle\lvert\theta_{k}\rvert^{2}\rangle=\frac{1}{\beta
A_{k}}=\frac{2TN}{h-2J\sum_{m=1}^{r}\cos\left(2\pi km\right)}.$ (51)

For large enough values of $d$ capturing the long-distance fall-off of the
correlations, the sum can be approximated by an integral with the new variable
of $y=k/N$,

$G(d)\approx\frac{T}{\pi}\int_{0}^{\infty}\mathrm{d}y\frac{e^{iy\,d}}{h-2J\sum_{m=1}^{r}\cos\left(\frac{2\pi
km}{N}\right)}.$ (52)

Expanding to leading order, the denominator of the integrand determines the
decay of the correlations

$\begin{split}G(d)&\approx\frac{T}{\pi}\int_{0}^{\infty}\mathrm{d}y\frac{e^{iy\,d}}{h-2rJ+\frac{1}{12}r(r+1)(2r+1)y^{2}}\approx\frac{T}{2\pi}\int_{0}^{\infty}\mathrm{d}y\frac{e^{iy\,d}}{h-2rJ+\frac{r^{3}}{6}y^{2}}\\\
&\sim
r^{-3}\int_{0}^{\infty}\mathrm{d}y\frac{e^{iy\,d}}{6\frac{h-2rJ}{r^{3}}+y^{2}}\propto
r^{-3}\xi\,e^{-d/\xi}\propto r^{-2}e^{-d/\xi},\end{split}$ (53)

with the correlation length given by the pole of the denominator close to the
critical point, where the order parameter can be linearized,
$\epsilon=h-2rJ\sim 2r\frac{t}{\tau_{Q}}$,

$\xi_{\mathrm{eq}}(r)\propto r^{3/2}\epsilon^{-1/2}.$ (54)

Note that the $r^{-2}$ constant factor emerged as the product of $r^{-3}$ in
Eq. (53) and the overall $r$ dependence of the correlation length stemming
from the regular part of the denominator, $\xi\sim r$.

Figure 10: Second cumulant of $n_{1}$ in the overdamped and underdamped
regimes for the KZ, intermediate, and adiabatic connectance regimes. Panels
$(a)-(c)$: Ratios of the first and second cumulants for the overdamped and
underdamped regimes for the three connectance regimes, respectively. The ratio
remains constant in the KZ and fast quench regimes. In the intermediate and
adiabatic regimes, the proportionality breaks down, and an approximate power
law relation is observed. Panels $(d)-(f)$: Second cumulants of $n_{1}$ as a
function of $\tau_{Q}$ in the underdamped regime $(a)$ showing the same power-
law and crossover behavior in the KZ regime while exhibiting universal
exponential decay in $(b)$ the intermediate and $(c)$ adiabatic ones. The
universal behavior is governed by a different $r$ dependence and only survives
up to $\tau_{Q}\lesssim 10$ in the intermediate regime. Panels $(g)-(i)$:
Similar findings in the overdamped regime. Simulations involve an ensemble of
$10^{4}$ trajectories.  Figure 11: Third cumulants of the excess energy
density as a function of $\tau_{Q}$ for system size $N=100,\,200,\,400$, and
connectance $c\in[0,1]$. Panels $(a)-(c)$: Cumulant ratios of
$\kappa_{3}(\rho_{E})/\rho_{E}$ both for the overdamped and underdamped cases,
for different values of $c$ in the KZ, intermediate and adiabatic regimes. In
the KZ regime, the cumulant ratios are independent of $\tau_{Q}$ and increase
universally with $c^{2}$. In the adiabatic regime, the ratio becomes
proportional to the square of the average decaying exponentially, as also
highlighted with the additional fittings. Panels $(d)-(f)$: Third cumulants in
the underdamped regime as a function of $\tau_{Q}$ in the three connectance
regimes. Inside the KZ regime, the same power law is displayed as for the
averages, but with a different universal rescaling with $r$. Remarkably, in
the $c\in[0.15,0.5]$ intermediate regime, perfect scaling collapse is found as
a function of $\tau_{Q}$. The adiabatic limit exhibits an exponential decay,
matching approximately the square of the average. Panels $(g)-(i)$: Similar
features with the corresponding power-laws in the overdamped regime,
exhibiting a faster adiabatic decay and less universal behavior for
intermediate connectances. The reported data is extracted from an ensemble of
$5\times 10^{4}$ trajectories.

As demonstrated in Fig. 2 in the main text, near the critical point,
$t_{c}\approx\tau/100$, the fitted exponential decay of the correlator
increases linearly with the connectance, up to high precision, as observed in
the insets for both the overdamped and underdamped regimes, verifying the
above result. Intuitively, the above picture can also be understood as the
correlation length being linearly scaled with the interaction range compared
to the one-dimensional case. Thus, near the critical point,
$\xi_{\mathrm{eq}}\sim r\left(t/\tau_{Q}\right)^{-1/2}\sim
r^{3/2}\epsilon^{-1/2}$. Here, the $3/2$ power balances for the additional
connectance dependence of the order parameter.

##  Appendix B Average defect and excess energy density in the underdamped
regime

In this appendix, we show the numerical results for the average defect density
in the underdamped regime. As shown in Fig. 9a,d,g, $n_{1},\,n_{2}$ and
$\rho_{E}$ follow precisely the predicted power laws as a function of
$\tau_{Q}$ and $r$, (22) in the KZ regime with $c\lesssim 0.15$,

$\displaystyle n_{1}\propto\left(r^{7/2}\tau_{Q}\right)^{-1/3},$ (55)
$\displaystyle\rho_{E}\sim n_{2}\propto\left(r^{1/2}\tau_{Q}\right)^{-1/3}.$
(56)

In the intermediate, $c\in[0.15,0.5]$ regime, universal signature only survive
up to $\tau_{Q}\lesssim 10$ as a function of $r^{1/3}\tau_{Q}$, as shown in
Fig. 9e,h, respectively for $n_{1}$ and $n_{2}$. In the case of the excess
energy, however, a remarkable scaling collapse is observed independently of
$r$, as shown in Fig. 9h. In the adiabatic limit, $c\gtrsim 0.5$, $n_{1}$ and
$n_{2}$ follow the same universal curve up to exponential accuracy as a
function of $r^{1/3}\tau_{Q}$, while $\rho_{E}$ becomes independent of $r$ as
governed by the dynamics of small angle deviations in defect-free regimes of
$\tau_{Q}$,

$\displaystyle n_{1}\sim n_{2}\sim e^{-0.7r^{1/2}\tau_{Q}},$ (57)
$\displaystyle\rho_{E}\sim e^{-0.06r^{1/4}\tau_{Q}}.$ (58)

##  Appendix C Second cumulant of the one-dimensional defects

In this appendix, a further numerical demonstration is provided for the
$\tau_{Q}$ and $r$ dependence of the second cumulants of the one-dimensional
defect density, $\kappa_{2}(n_{1})$, and the corresponding ratio
$\kappa_{2}(n_{1})/n_{1}$.

As demonstrated in Fig. 10a in the KZ regime, the proportionality is preserved
independently of $r$ with the same value as in the TFIM,
$\kappa_{2}(n_{1})/n_{1}\approx 0.578$. Remarkably, this universal behavior
survives for fast quenches as well. Additionally, between the KZ and fast
quench regimes, $r$-dependent non-universal peaks appear. In the intermediate
($c\in[0.15,0.5]$) and adiabatic ($c\gtrsim 0.5$) regimes, a sharp decrease is
displayed after similar non-universal peaks, as shown in Fig. 10b,c. This is
in good agreement with the power law relation between $\kappa_{2}(n_{1})$ and
$n_{1}$.

The dependence on $\tau_{Q}$ in the KZ regime is exhibited in Fig. 10d,g for
the overdamped and underdamped regimes, respectively, following the universal
power-laws predicted by the KZ mechanism in Eq. (22). In the intermediate
regime, universal scaling collapse can only be achieved for $\tau_{Q}\lesssim
10$ by the same regularity rescalings as in the adiabatic regime,
$\tau_{Q}\rightarrow r^{\frac{z\nu}{z\nu+1}}\tau_{Q}$. In the latter case,
curves with different values of $N$ and $r$ follow the same universal
exponential decay. These characteristics are displayed in Fig. 10e,f for the
overdamped case and in Fig. 10h,i for the underdamped case.

##  Appendix D Third cumulants

In this appendix, the third cumulants are investigated as a function of
$r,\,c$ and $\tau_{Q}$. As demonstrated in Fig. 11a, the excess energy
exhibits the same features in the KZ regime as for $\kappa_{2}(\rho_{E})$ in
Sec. V. The cumulant ratios follow approximately constant lines that get
shifted with $r$ beyond the fast quench regime. In this latter limit,
$\kappa_{3}(\rho_{E})/\rho_{E}$ converges to zero. As demonstrated in Fig.
11b,c, in the intermediate and adiabatic regimes of the connectance, the
ratios exhibit an exponential decay, as also indicated by the additional
fittings. This is in agreement with the power law relation between
$\kappa_{3}(\rho_{E})$ and $\rho_{E}$.

Furthermore, Fig. 11d,g shows the dependence on $\tau_{Q}$ following precisely
the universal power laws in the KZ regime in Eq. (43) for the overdamped and
underdamped regimes, respectively. In the intermediate regime with
$c\in[0.15,.5]$, universality breaks down around $\tau_{Q}\approx 10$, while
it survives in the underdamped regime, as shown in Fig. 11e and Fig. 11h,
respectively. In the adiabatic regime with $c\gtrsim 0.5$, a universal,
regularity-independent exponential decay is observed in agreement with the
results on $\kappa_{2}(\rho_{E})$ reported in Sec. IV.2.

In the case of $n_{1}$, similar features are found as in the case of
$\kappa_{2}(n_{1})$ in App. C. In the domain of the KZ power law scaling, the
$\kappa_{3}(n_{1})/n_{1}$ ratio follows precisely the same constant value as
in the TFIM independently of $r$. Similar to the $\kappa_{2}(n_{1})/n_{1}$
ratios, the proportionality also survives in the fast quench regime. A small
interval of non-universal peaks also appears before the KZ regime, as shown in
Fig. 12a. In the intermediate and adiabatic regimes, the ratios exhibit
exponential decay in agreement with the power-law relation, as shown in Fig.
12b,c.

As for the $\tau_{Q}$ dependence of $\kappa_{3}(n_{1})$, the numerical results
follow the same universal regularity and power law dependences in the KZ
regime as for the average given in Eq. (22). These features are demonstrated
in Fig. 12d,g. In the intermediate and adiabatic regimes, the universal
exponential decay is displayed in agreement with Eq. (35), as demonstrated in
Fig. 12e,f and Fig. 12h,i for the overdamped and underdamped regimes,
respectively.As in the case of $n_{1}$ and $\kappa_{2}(n_{1})$, the universal
scaling collapse in the intermediate regime as a function of
$r^{\frac{z\nu}{z\nu+1}}\tau_{Q}$ only survives up to $\tau_{Q}\lesssim 10$.

Figure 12: Third cumulant of $n_{1}$ in the overdamped and underdamped cases
for the KZ, intermediate, and adiabatic connectance regimes. Panels $(a)-(c)$:
Ratios of $\kappa_{3}(n_{1})/n_{1}$ for the overdamped and underdamped
regimes, for the three connectance regimes, respectively. The ratio remains
constant in the KZ and fast quench regimes. In the intermediate and adiabatic
regimes, the proportionality breaks down, and an approximate power law
relation is observed. Panels $(d)-(f)$: Third cumulants of $n_{1}$ as a
function of $\tau_{Q}$ in the underdamped regime showing the same power-law
and crossover behavior $(a)$ in the KZ regime, while exhibiting universal
exponential decay $(b)$ in the intermediate, $(c)$ and adiabatic regimes. The
universal behavior is governed by a different $r$ dependence and only survives
up to $\tau_{Q}\lesssim 10$ in the intermediate regime. Panels $(g)-(i)$:
Similar findings in the overdamped regime. The reported data is built from
$5\times 10^{4}$ trajectories.

## References

  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse ising model, Phys. Rev. E 58, 5355 (1998).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018).
  * Hauke  _et al._ [2020] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Perspectives of quantum annealing: methods and implementations, Reports on Progress in Physics 83, 054401 (2020).
  * Boixo  _et al._ [2013] S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A. Lidar, Experimental signature of programmable quantum annealing, Nature Communications 4, 2067 (2013).
  * Boixo  _et al._ [2016] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y. Smirnov, M. Mohseni, and H. Neven, Computational multiqubit tunnelling in programmable quantum annealers, Nature Communications 7, 10327 (2016).
  * Mishra  _et al._ [2018] A. Mishra, T. Albash, and D. A. Lidar, Finite temperature quantum annealing solving exponentially small gap problem with non-monotonic success probability, Nature Communications 9, 2917 (2018).
  * Bando  _et al._ [2020] Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J. Gómez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo, and H. Nishimori, Probing the universality of topological defect formation in a quantum annealer: Kibble-zurek mechanism and beyond, Phys. Rev. Research 2, 033369 (2020).
  * King  _et al._ [2022a] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Coherent quantum annealing in a programmable 2,000 qubit ising chain, Nature Physics 18, 1324 (2022a).
  * King  _et al._ [2024] A. D. King, A. Nocera, M. M. Rams, J. Dziarmaga, R. Wiersema, W. Bernoudy, J. Raymond, N. Kaushal, N. Heinsdorf, R. Harris, K. Boothby, F. Altomare, A. J. Berkley, M. Boschnak, K. Chern, H. Christiani, S. Cibere, J. Connor, M. H. Dehn, R. Deshpande, S. Ejtemaee, P. Farré, K. Hamer, E. Hoskinson, S. Huang, M. W. Johnson, S. Kortas, E. Ladizinsky, T. Lai, T. Lanting, R. Li, A. J. R. MacDonald, G. Marsden, C. C. McGeoch, R. Molavi, R. Neufeld, M. Norouzpour, T. Oh, J. Pasvolsky, P. Poitras, G. Poulin-Lamarre, T. Prescott, M. Reis, C. Rich, M. Samani, B. Sheldan, A. Smirnov, E. Sterpka, B. T. Clavera, N. Tsai, M. Volkmann, A. Whiticar, J. D. Whittaker, W. Wilkinson, J. Yao, T. J. Yi, A. W. Sandvik, G. Alvarez, R. G. Melko, J. Carrasquilla, M. Franz, and M. H. Amin, Computational supremacy in quantum simulation (2024), arXiv:2403.00910 [quant-ph] .
  * Kibble [1976] T. W. B. Kibble, Topology of cosmic domains and strings, Journal of Physics A: Mathematical and General 9, 1387 (1976).
  * Kibble [1980] T. Kibble, Some implications of a cosmological phase transition, Physics Reports 67, 183 (1980).
  * Zurek [1985a] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985a).
  * Dziarmaga [2010] J. Dziarmaga, Dynamics of a quantum phase transition and relaxation to a steady state, Advances in Physics 59, 1063–1189 (2010).
  * del Campo and Zurek [2014] A. del Campo and W. H. Zurek, Universality of phase transition dynamics: Topological defects from symmetry breaking, International Journal of Modern Physics A 29, 1430018 (2014).
  * Cui  _et al._ [2016] J.-M. Cui, Y.-F. Huang, Z. Wang, D.-Y. Cao, J. Wang, W.-M. Lv, L. Luo, A. del Campo, Y.-J. Han, C.-F. Li, and G.-C. Guo, Experimental trapped-ion quantum simulation of the kibble-zurek dynamics in momentum space, Scientific Reports 6, 33381 (2016).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, P. Zoller, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum kibble–zurek mechanism and critical dynamics on a programmable rydberg simulator, Nature 568, 207 (2019).
  * Cui  _et al._ [2020a] J.-M. Cui, F. J. Gómez-Ruiz, Y.-F. Huang, C.-F. Li, G.-C. Guo, and A. del Campo, Experimentally testing quantum critical dynamics beyond the kibble–zurek mechanism, Communications Physics 3, 44 (2020a).
  * Gardas  _et al._ [2018] B. Gardas, J. Dziarmaga, W. H. Zurek, and M. Zwolak, Defects in Quantum Computers, Sci. Rep. 8, 4539 (2018).
  * Weinberg  _et al._ [2020] P. Weinberg, M. Tylutki, J. M. Rönkkö, J. Westerholm, J. A. Åström, P. Manninen, P. Törmä, and A. W. Sandvik, Scaling and diabatic effects in quantum annealing with a d-wave device, Phys. Rev. Lett. 124, 090502 (2020).
  * Subires  _et al._ [2022] D. Subires, F. J. Gómez-Ruiz, A. Ruiz-García, D. Alonso, and A. del Campo, Benchmarking quantum annealing dynamics: The spin-vector langevin model, Phys. Rev. Res. 4, 023104 (2022).
  * Smolin and Smith [2013] J. A. Smolin and G. Smith, Classical signature of quantum annealing (2013), arXiv:1305.4904 [quant-ph] .
  * Wang  _et al._ [2013] L. Wang, T. F. Rønnow, S. Boixo, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Comment on: “Classical signature of quantum annealing”, arXiv e-prints , arXiv:1305.5837 (2013), arXiv:1305.5837 [quant-ph] .
  * Shin  _et al._ [2014] S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani, How ”quantum” is the d-wave machine? (2014), arXiv:1401.7087 [quant-ph] .
  * Albash  _et al._ [2015a] T. Albash, W. Vinci, A. Mishra, P. A. Warburton, and D. A. Lidar, Consistency tests of classical and quantum models for a quantum annealer, Phys. Rev. A 91, 042314 (2015a).
  * Albash  _et al._ [2015b] T. Albash, T. F. Rønnow, M. Troyer, and D. A. Lidar, Reexamining classical and quantum models for the d-wave one processor, The European Physical Journal Special Topics 224, 111 (2015b).
  * Albash and Marshall [2021] T. Albash and J. Marshall, Comparing relaxation mechanisms in quantum and classical transverse-field annealing, Phys. Rev. Appl. 15, 014029 (2021).
  * Bando and Nishimori [2021] Y. Bando and H. Nishimori, Simulated quantum annealing as a simulator of nonequilibrium quantum dynamics, Phys. Rev. A 104, 022607 (2021).
  * Rajak  _et al._ [2023] A. Rajak, S. Suzuki, A. Dutta, and B. K. Chakrabarti, Quantum annealing: an overview, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 381, 20210417 (2023).
  * Krapivsky [2010] P. L. Krapivsky, Slow cooling of an ising ferromagnet, Journal of Statistical Mechanics: Theory and Experiment 2010, P02014 (2010).
  * Jeong  _et al._ [2020] K. Jeong, B. Kim, and S. J. Lee, Nonequilibrium kinetics of excess defect generation and dynamic scaling in the ising spin chain under slow cooling, Phys. Rev. E 102, 012114 (2020).
  * Mayo  _et al._ [2021] J. J. Mayo, Z. Fan, G.-W. Chern, and A. del Campo, Distribution of kinks in an ising ferromagnet after annealing and the generalized kibble-zurek mechanism, Phys. Rev. Res. 3, 033150 (2021).
  * Fan  _et al._ [2023] Z. Fan, A. del Campo, and G.-W. Chern, Kibble-zurek mechanism for nonequilibrium generation of magnetic monopoles in spin ices (2023), arXiv:2307.05267 [cond-mat.stat-mech] .
  * Suzuki [2011] S. Suzuki, Kibble-zurek mechanism in simulated annealing and quantum annealing, Journal of Physics: Conference Series 302, 012046 (2011).
  * Zurek [1985b] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985b).
  * del Campo  _et al._ [2013] A. del Campo, T. W. B. Kibble, and W. H. Zurek, Causality and non-equilibrium second-order phase transitions in inhomogeneous systems, J. Phys.: Cond. Matt. 25, 404210 (2013).
  * del Campo  _et al._ [2022] A. del Campo, F. J. Gómez-Ruiz, and H.-Q. Zhang, Locality of spontaneous symmetry breaking and universal spacing distribution of topological defects formed across a phase transition, Phys. Rev. B 106, L140101 (2022).
  * Jaschke  _et al._ [2017] D. Jaschke, K. Maeda, J. D. Whalen, M. L. Wall, and L. D. Carr, Critical phenomena and kibble–zurek scaling in the long-range quantum ising chain, New Journal of Physics 19, 033032 (2017).
  * Sadhukhan  _et al._ [2020] D. Sadhukhan, A. Sinha, A. Francuz, J. Stefaniak, M. M. Rams, J. Dziarmaga, and W. H. Zurek, Sonic horizons and causality in phase transition dynamics, Phys. Rev. B 101, 144429 (2020).
  * Li  _et al._ [2022] B. W. Li, Y. K. Wu, Q. X. Mei, R. Yao, W. Q. Lian, M. L. Cai, Y. Wang, B. X. Qi, L. Yao, L. He, Z. C. Zhou, and L. M. Duan, Probing critical behavior of long-range transverse-field ising model through quantum kibble-zurek mechanism (2022), arXiv:2208.03060 [quant-ph] .
  * Chepiga and Mila [2021] N. Chepiga and F. Mila, Kibble-zurek exponent and chiral transition of the period-4 phase of rydberg chains, Nature Communications 12, 414 (2021).
  * Dutta and Dutta [2017] A. Dutta and A. Dutta, Probing the role of long-range interactions in the dynamics of a long-range kitaev chain, Phys. Rev. B 96, 125113 (2017).
  * Caneva  _et al._ [2008] T. Caneva, R. Fazio, and G. E. Santoro, Adiabatic quantum dynamics of the lipkin-meshkov-glick model, Phys. Rev. B 78, 104426 (2008).
  * Acevedo  _et al._ [2014] O. L. Acevedo, L. Quiroga, F. J. Rodríguez, and N. F. Johnson, New dynamical scaling universality for quantum networks across adiabatic quantum phase transitions, Phys. Rev. Lett. 112, 030403 (2014).
  * Defenu  _et al._ [2018] N. Defenu, T. Enss, M. Kastner, and G. Morigi, Dynamical critical scaling of long-range interacting quantum magnets, Phys. Rev. Lett. 121, 240403 (2018).
  * Gherardini  _et al._ [2023] S. Gherardini, L. Buffoni, and N. Defenu, Universal defects statistics with strong long-range interactions (2023), arXiv:2305.11771 [quant-ph] .
  * Silvi  _et al._ [2016] P. Silvi, G. Morigi, T. Calarco, and S. Montangero, Crossover from classical to quantum kibble-zurek scaling, Phys. Rev. Lett. 116, 225701 (2016).
  * Dutta  _et al._ [2016] A. Dutta, A. Rahmani, and A. del Campo, Anti-Kibble-Zurek Behavior in Crossing the Quantum Critical Point of a Thermally Isolated System Driven by a Noisy Control Field, Phys. Rev. Lett. 117, 080402 (2016).
  * Singh and Gangadharaiah [2021] M. Singh and S. Gangadharaiah, Driven quantum spin chain in the presence of noise: Anti-kibble-zurek behavior, Phys. Rev. B 104, 064313 (2021).
  * King  _et al._ [2022b] E. C. King, M. Kastner, and J. N. Kriel, Long-range kitaev chain in a thermal bath: Analytic techniques for time-dependent systems and environments (2022b), arXiv:2204.07595 [quant-ph] .
  * Singh  _et al._ [2023] M. Singh, S. Dhara, and S. Gangadharaiah, Driven one-dimensional noisy kitaev chain, Phys. Rev. B 107, 014303 (2023).
  * Puebla  _et al._ [2020] R. Puebla, A. Smirne, S. F. Huelga, and M. B. Plenio, Universal anti-kibble-zurek scaling in fully connected systems, Phys. Rev. Lett. 124, 230602 (2020).
  * del Campo  _et al._ [2010] A. del Campo, G. De Chiara, G. Morigi, M. B. Plenio, and A. Retzker, Structural defects in ion chains by quenching the external potential: The inhomogeneous kibble-zurek mechanism, Phys. Rev. Lett. 105, 075701 (2010).
  * Chiara  _et al._ [2010] G. D. Chiara, A. del Campo, G. Morigi, M. B. Plenio, and A. Retzker, Spontaneous nucleation of structural defects in inhomogeneous ion chains, New Journal of Physics 12, 115003 (2010).
  * Nigmatullin  _et al._ [2016] R. Nigmatullin, A. del Campo, G. De Chiara, G. Morigi, M. B. Plenio, and A. Retzker, Formation of helical ion chains, Phys. Rev. B 93, 014106 (2016).
  * Ulm  _et al._ [2013] S. Ulm, J. Roßnagel, G. Jacob, C. Degünther, S. T. Dawkins, U. G. Poschinger, R. Nigmatullin, A. Retzker, M. B. Plenio, F. Schmidt-Kaler, and K. Singer, Observation of the kibble–zurek scaling law for defect formation in ion crystals, Nature Communications 4, 2290 (2013).
  * Pyka  _et al._ [2013] K. Pyka, J. Keller, H. L. Partner, R. Nigmatullin, T. Burgermeister, D. M. Meier, K. Kuhlmann, A. Retzker, M. B. Plenio, W. H. Zurek, A. del Campo, and T. E. Mehlstäubler, Topological defect formation and spontaneous symmetry breaking in ion coulomb crystals, Nature Communications 4, 2291 (2013).
  * Hamp  _et al._ [2015] J. Hamp, A. Chandran, R. Moessner, and C. Castelnovo, Emergent coulombic criticality and kibble-zurek scaling in a topological magnet, Phys. Rev. B 92, 075142 (2015).
  * Du  _et al._ [2023] K. Du, X. Fang, C. Won, C. De, F.-T. Huang, W. Xu, H. You, F. J. Gómez-Ruiz, A. del Campo, and S.-W. Cheong, Kibble–zurek mechanism of ising domains, Nature Physics 19, 1495 (2023).
  * Biroli  _et al._ [2010] G. Biroli, L. F. Cugliandolo, and A. Sicilia, Kibble-zurek mechanism and infinitely slow annealing through critical points, Phys. Rev. E 81, 050101 (2010).
  * Ma  _et al._ [2024] T.-C. Ma, H.-Q. Shi, H.-Q. Zhang, and A. del Campo, Universal critical holography and domain wall formation (2024), arXiv:2406.05167 [cond-mat.stat-mech] .
  * Libál  _et al._ [2020] A. Libál, A. del Campo, C. Nisoli, C. Reichhardt, and C. J. O. Reichhardt, Quenched dynamics of artificial colloidal spin ice, Phys. Rev. Res. 2, 033433 (2020).
  * Samajdar and Huse [2024] R. Samajdar and D. A. Huse, Quantum and classical coarsening and their interplay with the kibble-zurek mechanism (2024), arXiv:2401.15144 [quant-ph] .
  * Reichhardt  _et al._ [2022] C. J. O. Reichhardt, A. del Campo, and C. Reichhardt, Kibble-zurek mechanism for nonequilibrium phase transitions in driven systems with quenched disorder, Communications Physics 5, 173 (2022).
  * Reichhardt and Reichhardt [2023] C. Reichhardt and C. J. O. Reichhardt, The kibble-zurek scenario and coarsening across nonequilibrium phase transitions in driven vortices and skyrmions (2023), arXiv:2307.02685 [cond-mat.mes-hall] .
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Chomaz  _et al._ [2022] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Dipolar physics: a review of experiments with magnetic quantum gases, Reports on Progress in Physics 86, 026401 (2022).
  * Periwal  _et al._ [2021] A. Periwal, E. S. Cooper, P. Kunkel, J. F. Wienand, E. J. Davis, and M. Schleier-Smith, Programmable interactions and emergent geometry in an array of atom clouds, Nature 600, 630 (2021).
  * Lucas [2014] A. Lucas, Ising formulations of many np problems, Frontiers in Physics 2 (2014).
  * Newman [2018] M. Newman,  _Networks_ (Oxford University Press, 2018).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Reports on Progress in Physics 29, 255 (1966).
  * Kloeden and Platen [2011] P. Kloeden and E. Platen,  _Numerical Solution of Stochastic Differential Equations_ , Stochastic Modelling and Applied Probability (Springer Berlin Heidelberg, 2011).
  * Laguna and Zurek [1998] P. Laguna and W. H. Zurek, Critical dynamics of symmetry breaking: Quenches, dissipation, and cosmology, Phys. Rev. D 58, 085021 (1998).
  * Zeng  _et al._ [2023] H.-B. Zeng, C.-Y. Xia, and A. del Campo, Universal breakdown of kibble-zurek scaling in fast quenches across a phase transition, Phys. Rev. Lett. 130, 060402 (2023).
  * Kou and Li [2023] H.-C. Kou and P. Li, Varying quench dynamics in the transverse ising chain: The kibble-zurek, saturated, and presaturated regimes, Phys. Rev. B 108, 214307 (2023).
  * del Campo [2018] A. del Campo, Universal statistics of topological defects formed in a quantum phase transition, Phys. Rev. Lett. 121, 200601 (2018).
  * Cui  _et al._ [2020b] J.-M. Cui, F. J. Gómez-Ruiz, Y.-F. Huang, C.-F. Li, G.-C. Guo, and A. del Campo, Experimentally testing quantum critical dynamics beyond the kibble-zurek mechanism, Communications Physics 3, 44 (2020b).
  * Gómez-Ruiz  _et al._ [2020] F. J. Gómez-Ruiz, J. J. Mayo, and A. del Campo, Full counting statistics of topological defects after crossing a phase transition, Phys. Rev. Lett. 124, 240602 (2020).
  * Gómez-Ruiz  _et al._ [2022] F. J. Gómez-Ruiz, D. Subires, and A. del Campo, Role of boundary conditions in the full counting statistics of topological defects after crossing a continuous phase transition, Phys. Rev. B 106, 134302 (2022).
  * Thudiyangal and del Campo [2024] M. Thudiyangal and A. del Campo, Universal vortex statistics and stochastic geometry of bose-einstein condensation (2024), arXiv:2401.09525 [cond-mat.quant-gas] .
  * Xia  _et al._ [2024] C.-Y. Xia, H.-B. Zeng, A. Grabarits, and A. del Campo, Kibble-zurek mechanism and beyond: Lessons from a holographic superfluid disk (2024), arXiv:2406.09433 [cond-mat.stat-mech] .
  * Cummins and Berloff [2024] J. S. Cummins and N. G. Berloff, Vector ising spin annealer for minimizing ising hamiltonians (2024), arXiv:2403.16608 [quant-ph] .