# Phase transition in fluctuations of interacting spins at infinite
temperature

V. N. Gorshkov  National Technical University of Ukraine, 37 Prospect
Peremogy, Kyiv 03056, Ukraine  N. A. Sinitsyn  T-4, Los Alamos National
Laboratory, Los Alamos NM 87545  D. Mozyrsky<EMAIL_ADDRESS>T-4, Los
Alamos National Laboratory, Los Alamos NM 87545

(August 28, 2024)

###### Abstract

The high temperature limit of interacting spins is usually not associated with
ordering or critical phenomena. Nevertheless, spontaneous fluctuations of a
local spin polarization at equilibrium have nontrivial dynamics even in this
limit. Here, we demonstrate that the spin noise power spectrum of these
fluctuations can undergo discontinuous changes as a function of an external
magnetic field. As a simple illustration, we consider a model of Ising-like
long range spin-spin interactions with a transverse magnetic field as a
control parameter. This system undergoes a phase transition associated with
disappearance of the noise power peak responsible for the most detrimental
decoherence effect of the interactions.

A phase transition is a discontinuous change of some measurable characteristic
of a many-body system. Traditionally, phase transitions in condensed matter
have been associated with sharp changes of a long range order at equilibrium
when a control parameter changes across a critical value. The phase
transitions are observed as divergence of linear response characteristics,
such as the magnetic susceptibility near the critical point.

Any long-range ordering, however, is possible only at a sufficiently low
temperature, $T$. Moreover, in the large temperature limit, the measurable
linear response characteristics are vanishing because the changes of
microstate probabilities in the equilibrium Gibbs distribution are suppressed
by the factor $1/(k_{B}T)$.

In this letter, contrarily to such expectations, we argue that interacting
spin systems may show the basic features of a phase transition even in the
large temperature limit at thermodynamic equilibrium. This critical phenomenon
should be accessible for experimental studies by spin noise spectroscopy (SNS)
SNS , which probes dynamics of mesoscopic spin fluctuations at equilibrium.
The basic measurable characteristic for the SNS is the noise power spectrum.
It is formally related to the dynamic magnetic susceptibility, as articulated
by the Fluctuation Dissipation Theorem SNS . However, unlike the
susceptibility, the power spectrum is not suppressed and has nontrivial
features even at the infinite temperature.

For demonstration, we suggest Metal–Organic Frameworks – the organic materials
that incorporate metallic ions with uncompensated spins in a regular array MOF
. The atomic spins can be placed at a sufficient distance from each other to
remove exchange interactions. Thus, such spins form a macroscopic array of
qubits with considerable coherence time. At cryogenic temperatures, $T<4K$,
the phonon-related mechanisms of decoherence in molecular spin qubits are
suppressed. The main sources of the decoherence are then hyperfine, of order
of $\sim 10-10^{2}$G, and magnetic dipole fields from the neighboring spin
qubits ($\sim 10^{2}-10^{3}$G) mol-nano . The hyperfine effects can be
suppressed further by isotopic purification purified , but the random long-
range dipole fields represent an unsolved problem for quantum control. Thus,
large qubit systems induce collective decoherence effects that are absent on
the level of a single qubit.

Here, we explore the coherent behavior of long-range interacting spins in a
moderate external magnetic field. The most unfortunate effect of the dipole
spin-spin interactions is the relatively slow collective relaxation of the
local spontaneous spin field fluctuations. The latter mis-align the net
effective fields that act on individual spins from the external field axis.
Therefore, there is no possibility to define a single rotation axis for the
spins.

The main finding of our theory is that this most detrimental decoherence
effect disappears at relatively week, of an order of the typical dipole field
fluctuation, value of the external field. This is manifested as a
discontinuous disappearance of one of the peaks in the noise power spectrum at
a critical value of the external field. Our theory demonstrates other
similarities with the conventional phase transitions, such as divergence of
the spin-lattice relaxation rate near a critical point.

As a theoretical model we consider the Hamiltonian with long-range Ising-like
spin-spin interactions in a controllable transverse magnetic field, $f$:

$H=\sum_{i<j}J_{ij}\sigma^{z}_{i}\sigma^{z}_{j}+f\sum_{i}\sigma^{x}_{i}\,,$
(1)

where $\sigma_{z}^{i}$ and $\sigma_{x}^{i}$ are the $z$ and $x$ components for
the $i$’th spin, $J_{ij}$ are coupling constants ab .

The temperatures $T\sim 1$K are much larger, in energy scale, than the dipole
fields. Hence, there is no spin ordering at the thermodynamic equilibrium, and
the equilibrium spin density matrix is proportional to a unit matrix. The
local spontaneous spin fluctuations create an uncompensated local spin
polarization, whose dynamics is characterized by an ensemble averaged spin
auto-correlation function,

$C(t)\equiv\frac{1}{2}\langle\\{\sigma^{z}_{i}(0),\sigma^{z}_{i}(t)\\}\rangle,$
(2)

where $\\{\ldots\\}$ stands for anticommutator and the averaging is over the
measurement results of the spin and over all possible initial spin
configurations in the system at $t=0$. This correlator is directly related to
the noise power spectrum SNS-exp ; SNS-exp-rev ; SNS :

$g(\omega)=\int dtC(t)e^{i\omega t},\quad\omega>0.$ (3)

We will assume that $J_{ij}$ is long range and so each spin in the sum in Eq.
(1) interacts with a large number of spins in its neighborhood. Therefore the
effective magnetic field, acting on the $i$-th spin along z-direction,

$h_{i}=\sum_{j\neq i}J_{ij}\sigma^{z}_{j}\,,$ (4)

can be viewed as a slow (semi)classical variable. This approximation was
justified for the large temperature limit in QD1 ; QD2 . The spins then follow
the Larmor equation of motion:

${\dot{\sigma}}_{i}={\bf\Omega}_{i}\times{\bf\sigma}_{i}\,,$ (5)

where ${\bf\Omega}_{i}$ is a 3-component vector,
${\bf\Omega}_{i}=(2f,0,2h_{i})$; $h_{i}$ can be time-dependent and should be
determined self-consistently.

If fields $h_{i}$ were constant (in time), the spin correlator would be the
sum of a constant and an oscillating component:

$C^{0}(t)=\frac{h^{2}_{i}}{h^{2}_{i}+f^{2}}+\frac{f^{2}}{h^{2}_{i}+f^{2}}\cos\Big{(}2t\sqrt{h^{2}_{i}+f^{2}}\Big{)}\,.$
(6)

The Fourier transform (3) would have then two delta-peaks: at $\omega=0$ and
at $\omega=2\sqrt{h^{2}_{i}+f^{2}}$. The spectral weight
$h^{2}_{i}/(h^{2}_{i}+f^{2})$ of the $\omega=0$ peak represents the fraction
of spin polarization that does not experience the relaxation, whereas
$f^{2}/(h_{i}^{2}+f^{2})$ is associated with the Larmor precession around
${\bm{\Omega}}_{i}$.

Since $h_{i}$ are different for different $i$, one should average Eq. (6) over
$h_{i}$. The width of the Larmor peak would then correspond to the standard
inhomogeneous broadening. The averaging over $h$-fields, however, would only
modify the overall spectral weight of the zero-frequency peak, without any
broadening to it.

In Fig. 1 we show the results of our numerical simulations for a spin chain
with long-range interactions in Eq. (6) suppl . Figure 1 confirms the presence
of both noise power peaks for relatively small $f$. However, for various
choices of the coupling distribution, inevitably, at the external field above
some critical value, the zero-frequency peak disappears completely.

Figure 1:  The numerically generated (see supplementary file suppl  for
details) spectrum of the spin auto-correlation function in the spin chain at
different transverse fields $f$. Curve 1: $f=0.5$; Curve 2: $f=0.75$; Curve 3:
$f=1.0$; Curve 4: $f=1.4$; Curve 5: $f=1.6$; Curve 6: $f=1.8$ ; Curve 7:
$f=2.0$. For all curves $\sum_{j}J_{ij}^{2}=W^{2}/4$, where $W=3.62$.

This observation was not expected because within the simple picture in which
the fields $h_{i}$ are static, the zero-frequency peak does not disappear even
in very large external fields. For $f\rightarrow\infty$, the ratio of the
weights of the delta-peaks from (6) decays as a power law $\sim
h_{i}^{2}/f^{2}$. Hence, it seems that in order to just suppress the
correlated lattice effects, the external field $f$ must be considerably larger
than the equilibrium root-mean-square of $h_{i}$ at $f=0$. However, Fig. 1
shows that the zero-frequency peak disappears at some critical $f=f_{c}$,
which is comparable to the typical $|h_{i}|$.

This resembles the behavior of a magnetization near the critical point of a
ferromagnetic phase transition. The role of the order parameter in our case is
played by the weight of the zero-frequency peak, which is the noise power
integrated over the frequency range where this peak has noticeably nonzero
values. This weight is finite for $f<f_{c}$ but is zero otherwise.

The analogy with conventional phase transitions can be enhanced by considering
behavior of this peak’s spectral weight for $f$ near $f_{c}$. The field
$h_{i}$ in Eq. (6) has a quasi-static and a quickly fluctuating, due to fast
spin rotations, components. The quasi-static one, $h$, varies from one spin to
another randomly. If $J_{ij}$ are long range, then the distribution function
$P(h)$ for the $h$-fields can be approximated by a Gaussian,

$P(h)=\frac{1}{\sqrt{2\pi h_{0}^{2}}}\ e^{-\frac{h^{2}}{2h_{0}^{2}}}\,,$ (7)

where $h_{0}^{2}$ is the variance of the quasi-static interaction field. On
the other hand, from the definition of $C(t)$ and $h_{i}$, Eqs. (2, 4), and
the interpretation of the first term in (6) as the quasi-static part of the
field contribution to the noise power, we find that $\langle h_{0}^{2}\rangle$
is expressible via the weight of the zero-frequency peak:

$h_{0}^{2}=\sum_{j,\,j\neq i}J_{ij}^{2}\,\frac{h^{2}_{j}}{h^{2}_{j}+f^{2}}\,,$

where we used the fact that there are no same-time correlations between
different spins at high temperatures. This leads to a self-consistency
equation on $h_{0}$:

$\frac{W^{2}}{4}\,\int dhP(h)\,\frac{h^{2}}{h^{2}+f^{2}}=h_{0}^{2},$ (8)

where $W^{2}\equiv 4\sum_{j\neq i}J_{ij}^{2}$.

Figure 2:  Numerical solution of Eq. (8) for $W=3.62$.

Equation (8) can be readily solved numerically. It has always a trivial
solution $h_{0}=0$, but for sufficiently small $f$ it also has a non-zero
solution which is shown in Fig. 2. For $f\rightarrow 0$, $h_{0}\rightarrow
W/2$, which is the typical spin-spin field fluctuation at $f=0$. As $f$
increases, the distribution $P(h)$ narrows until its width $h_{0}$ is $0$.
This point corresponds to the critical field: $f_{c}=W/2$. For values of
$f>f_{c}$, only the trivial solution exists, whereas at slightly smaller $f$
(than $f_{c}$), the non-zero solution is

$h_{0}=f_{c}\sqrt{\delta/3},\quad\delta=(f_{c}^{2}-f^{2})/f_{c}^{2}\,.$ (9)

This (square root) behavior is consistent with the presence of a 2nd order
phase transition, with $h_{0}$ being the order parameter. The self-consistency
equation yields $f_{c}\simeq 1.81$. This result is in a very good agreement
with numerical simulations: Comparing such estimates with Fig. 1, the
transition takes place between curves 6 and 7, which correspond to the values
of $f=1.8$ and $f=2.0$, respectively.

The disappearance of the zero-frequency peak does not mean that the spin
interaction effects are not present. For $f>f_{c}$, they are still manifested
in the broadening of the Larmor peak. However, this effect is very different
from the mis-aligning of quasi-static fields. The spin interaction effects
resemble now the effect of a fast noise that produces Lorentzian-like
broadening of the Lamor precessions around the external field direction.

Broadening of the zero-frequency peak: The delta-peak at $\omega=0$ associated
with the first term in Eq. (6) has zero-width, which is not the case with the
curves in Fig. 1. This broadening is the consequence of the dynamics of
$h_{i}$. Thus, the phase transition and the broadening of the delta-peak have
the same origin.

Not only the zero-frequency peak looses power with growing $f$ but also it
flattens as $f\rightarrow f_{c}$ (Fig. 1). This suggests that the width of
this peak has critical behavior near the phase transition as well. To confirm
this, we rewrite Eq. (5) as an equation that contains only $\sigma^{z}_{i}$
components of the spins. For that we solve the equations containing the time
derivatives of $\sigma^{x}_{i}$ and $\sigma^{y}_{i}$ in terms of $h_{i}$ and
$\sigma^{z}_{i}$ and substitute these solutions (for $\sigma^{x}_{i}$ and
$\sigma^{y}_{i}$) into the remaining equation for $\sigma^{z}_{i}$. The result
is

${\dot{\sigma}^{z}_{i}}=-4f^{2}\int_{0}^{t}dt^{\prime}\sigma^{z}_{i}(t^{\prime})\,\cos\Big{[}2\int_{t^{\prime}}^{t}d\tau
h_{i}(\tau)\Big{]}\,,$ (10)

where $h_{i}$ is given by Eq. (4).

We separate the field $h_{i}(\tau)$ in Eq. (10) into low frequency (slow) and
high frequency (fast) components:

$h(\tau)=h_{s}+h_{f}(\tau)\,,$ (11)

with $h_{s}$ being nearly a constant (yet, a random number distributed
according to the distribution $P(h_{s})$, e.g. Eq. (7)) and $h_{f}(\omega)$
having non-zero Fourier components only at $|\omega|\geq 0$. That is, the
power spectrum of the latter,
$g_{f}(\omega)=\langle|h_{f}(\omega)|^{2}\rangle$, can be associated with the
part of the spectrum (e.g. Fig. 1) that includes the peak at finite (i.e.
Larmor) frequencies. This part corresponds, approximately, to the Fourier
transform of the second term in Eq. (2), averaged over the distribution
$P(h_{s})$ suppl .

To obtain the low frequency spin dynamics, we apply the same frequency
separation to $\sigma^{z}_{i}$, as we did for $h$ in (11). Then, when writing
Eq. (10) for a slow spin (in the following denoted by $\sigma_{s}^{z}$), we
can average over the rapidly varying field $h_{f}$. By virtue of the central
limit theorem the latter can be assumed to be a Gaussian variable (in the
limit of long range coupling), and therefore, upon averaging over $h_{f}$ and
assuming that $h_{s}$ is a constant, the cosine in Eq. (10) transforms into a
function of $t-t^{\prime}$; (see suppl  for explicit details). The resulting
equation can be readily solved by the Laplace transform and one obtains the
fluctuation spectrum of slow spins,
$g_{s}(\omega,h_{s})=\langle|\sigma_{s}^{z}(\omega,h_{s})|^{2}\rangle$. In the
supplemental materials we show that spectrum $g_{s}(\omega,h_{s})$ is a
Lorentzian, with the width being a function of the slow field $h_{s}$ suppl .

The spectrum of the spin autocorrelation function at sufficiently small
frequencies can be obtained by averaging $g_{s}(\omega,h_{s})$ with respect to
the distribution of slow fields $h_{s}$ according to Eq. (7), and accounting
for that the slow spins constitute a fraction $h_{s}^{2}/(h_{s}^{2}+f^{2})$ of
all spins, e.g. the first term in Eq. (6),

$g(\omega)=\int
dh_{s}P(h_{s})\,\frac{h_{s}^{2}}{f^{2}+h_{s}^{2}}\,g_{s}(\omega,h_{s})\,.$
(12)

Noise spectrum at small $\omega$ near the transition point: Since at the
transition point (i.e., for $f\rightarrow f_{c}$) the order parameter $h_{0}$
vanishes (i.e., $P(h_{s})$ approaches delta-function), we can evaluate the
spectrum $g(\omega)$ by setting $h_{s}=0$ in $g_{s}(\omega,h_{s})$. Evaluation
of the Lorentzian $g_{s}(\omega,0)$ is presented in the supplemental materials
suppl , and we quote the final result

$g(\omega)=\frac{\delta}{3}\,\frac{2\sqrt{2\pi}f_{c}}{\delta^{2}\omega^{2}+2\pi
f_{c}^{2}}\,.$ (13)

Equation (13) clearly has properties consistent with those for the numerical
power spectra presented in Fig. 1: as $f$ approaches $f_{c}$, the integrated
weight of the spectrum in the interval between $0$ and a frequency of the
order of $O(f)$ is of the order of $\delta\ll 1$. Furthermore, as $f$
approaches $f_{c}$, the spectrum of the slow spins flattens and, at the
critical point, becomes equal to $0$ at the final interval at sufficiently low
frequencies. Barring a contradiction of terms, the relaxation rate,
$\sqrt{2\pi}f_{c}/\delta$, diverges for slow spins. This is formally opposite
to the divergent correlation length in case of conventional 2nd order phase
transitions.

Noise spectrum at small $\omega$ for $f\ll f_{c}$: At sufficiently small $f$,
the typical values of $h_{s}$ greatly exceed $f$ and therefore the relative
fraction of slow spins, i e. the fraction in the integrand in Eq. (12) can be
replaced by $1$. The width $\Gamma$ of the Lorentzian $g_{s}(\omega,h_{s})$ is
evaluated in the supplemental materials suppl : It turns out to be a strong
(Gaussian) function of $h_{s}$. Then, changing integration variable in Eq.
(12) from $h_{s}$ to $\Gamma(h_{s})$, with logarithmic accuracy we obtain

$g(\omega)\sim\int_{0}^{\Gamma(0)}d\Gamma\,\frac{\Gamma^{\alpha}}{\omega^{2}+\Gamma^{2}}\,,$
(14)

where $\alpha=(W^{2}-4h_{0}^{2})/4h_{0}^{2}$ and $\Gamma(0)\sim
f^{3/2}/W^{1/2}$; see suppl  for the details of derivation.

We see that for $\alpha\leq 1$ and $|\omega|\ll\Gamma(0)$, Eq. (14) yields
$g(\omega)\sim|\omega|^{\alpha-1}$. This power law dependence of the
fluctuation spectrum is consistent with the shapes of the spectra obtained by
numerical simulations in Fig. 1. Indeed, for $f\ll W$,
$\alpha\simeq\sqrt{2\pi}f/W$, and therefore the power low dependence (with the
exponent $\alpha-1$) suggests that the spectra with smaller $f$ tend to
$\infty$ as $|\omega|\rightarrow 0$ more steeply than those with greater $f$.
This is exactly the case with the curves 1 and 2 in Fig. 1.

Apart from the phase transition, at $\alpha=1$ the system exhibits a
crossover: For $\alpha<1$, the fluctuation spectrum is divergent at zero
frequency, which indicates that some spins/spin domains (with sufficiently
high $h_{s}$) do not evolve and remain frozen indefinitely. For $\alpha>1$,
Eq. (14) suggests that such picture breaks down (so that all spins have finite
relaxation time). This happens when $h_{0}=W/\sqrt{8}$, which corresponds to
$f\simeq 0.9$ (for $W=3.62$), e.g. Fig. 2. This value is clearly in agreement
with what we observe in Fig. 1: While at $\omega=0$ the spectrum for the curve
2 with $f=0.75$ is divergent, it is finite for curve 3 with $f=1$.

The origins for the phase transition and the crossover can be understood if we
consider a model with infinite range coupling ($J_{ij}=J\sim
1/\sqrt{N}=const$) LMG . This model is integrable and, although missing the
effects of disorder, provides an opportunity for a detailed, nearly exact
analysis of the critical behavior. In the supplemental materials suppl  we
show that the dynamics of the effective field $h$ (i.e. magnetization) in this
model corresponds to the dynamics of a particle moving in a Ginzburg-Landau-
like potential. The phase transition can be related to the
emergence/disappearance of two potential minima/wells, each corresponding to a
fixed point for $h$. On the other hand, the crossover to the state with frozen
magnetization is related to the existence of a classical trajectory that
connects between these two wells. That is, for sufficiently small values of
the transverse field $f$, the “particle” is bound to travel within a single
well, while for $f$ higher than a certain threshold value, it can transition
between the two wells.

In conclusion, direct measurements of the spin fluctuations can reveal unusual
dynamic phase transitions that have no analogs in the conventional theory of
critical phenomena. Our phase transition is associated with the disappearance
of the low frequency fluctuations at certain (critical) value of the applied
magnetic field. Also, at somewhat lower values of the field the system
exhibits a crossover, which can be associated with the emergence/disappearance
of very slowly relaxing states (i.e. sufficiently large spin domains) related
to the existence of the dynamically stable fixed points for system’s
magnetization. suppl .

The study of similar phase transitions in many-body fluctuations should be of
practical importance. In addition to the qubit control applications, such
critical points reveal nonperturbative, including topological, changes in the
energy functional even though the system does not follow its mimimum. It
should be insightful to relate such phenomena to the properties of the ordered
phases that emerge in the same systems at finite temperatures, and thus
predict the phase diagrams using only easy-to-detect high-temperature noise
characteristics.

###### Acknowledgements.

Acknowledgements: We thank A. Balatsky, R. Malla and I. Martin for valuable
discussions. The work at LANL is supported by LDRD program.

## References

  * (1) N. A. Sinitsyn, and Y. V. Pershin, “The Theory of Spin Noise Spectroscopy: A Review”, Rep. Prog. Phys. 79, 106501 (2016).
  * (2) C.-J. Yu, S. von Kugelgen, M. D. Krzyaniak, W. Ji, W. R. Dichtel, M. R. Wasielewski, D. E. Freedman, “Spin and Phonon Design in Modular Arrays of Molecular Qubits”, Chem. Mater. 32, 23, 10200 (2020).
  * (3) N. V. Prokof’ev and P. C. E. Stamp, “Quantum Relaxation of Ensembles of Nanomagnets”, Journal of Low Temperature Physics 113, 1147 (1998).
  * (4) W. Wernsdorfer,“Quantum dynamics in molecular nanomagnets”, Comptes Rendus Chimie 11, 1086 (2008).
  * (5) Similar spin Hamiltonian has been considered in connection with time crystals; see e.g., W. W. Ho, S. Choi, M. D. Lukin, D. A. Abanin, “Critical Time Crystals in Dipolar Systems”, Phys. Rev. Lett. 119, 010602 (2017).
  * (6) S. A. Crooker, D. G. Rickel, A. V. Balatsky, and D. L. Smith, “Spectroscopy of spontaneous spin noise as a probe of spin dynamics and magnetic resonance”, Nature 431, 49-52 (2004).
  * (7) V. S. Zapasskii, “Spin-noise spectroscopy: from proof of principle to applications”, Advances in Optics and Photonics 5, 131 (2013).
  * (8) K. A. Al-Hassanieh, V. V. Dobrovitski, E. Dagotto, and B. N. Harmon, “Numerical Modeling of the Central Spin Problem Using the Spin-Coherent-State P-Representation”, Phys. Rev. Lett. 97, 037204 (2006).
  * (9) A. Fischer, I. Kleinjohann, N. A. Sinitsyn, and F. B. Anders, “Cross-correlation spectra in interacting quantum dot systems”, Phys. Rev. B 105, 035303 (2022).
  * (10) Supplemental material, with additional references LL1 ; LMG , describes our numerics and defines different tested distributions of $J_{ij}$. It also provides a discussion of the phase transition in exactly solvable Lipkin-Meshkov-Glick model, and describes self-consistent calculations of the power spectrum in the entire range of frequencies.
  * (11) N. A. Sinitsyn, V. V. Dobrovitski, S. Urazhdin, A. Saxena, “Geometric Control Over the Motion of Magnetic Domain Walls”, Phys. Rev. B 77, 212405 (2008).
  * (12) Y. Huang, T. Li, Z. Yin,“Symmetry-breaking dynamics of the finite-size Lipkin-Meshkov-Glick model near ground state”, Phys. Rev. A 97, 012115 (2018).

##  I Supplementary material for “Phase transition in fluctuations of
interacting spins at infinite temperature”

##  II Numerical simulations

We have carried out numerical simulations for the system of $1000$ classical
spins on 1-dimensional ring (i.e., with periodic boundary conditions) by
preparing spins with random initial $\sigma^{z}_{i}(0)=\pm 1$ at equal weights
(e.g. high temperature limit), and with
$\sigma^{x}_{i}(0)=\sigma^{y}_{i}(0)=0$. We found that the choice of $J_{ij}$
does not change the behaviour of the correlation function qualitatively as
long as the spin-spin coupling is sufficiently long range, e.g. Figs. 1-3 of
this supplementary. We also calculated the auto-correlation functions for the
effective fields $h_{i}(t)$, which show behaviour very similar to the spin
auto-correlation functions, signalling that correlations between different
spins are relatively weak, (see, e.g., the insets in Figs. 1-3). Furthermore,
we observe in Figs. 1-3 that at small frequencies the behaviour of the spin
auto-correlation functions for different choices of spin-spin couplings (i.e.
red, blue and green curves) are very close even at the quantitative level.
This is in agreement with out theoretical expectations that the small
frequency behaviour of the fluctuation spectra is controlled by the quantity
$\sum_{j}J_{ij}^{2}$, rather than by the specific functional choice of
couplings between spins.

In Fig. 1 of the main text, we present the Fourier transform of the
correlation function, $g(\omega)$, at different values of $f$. For relatively
small $f$, e.g. curves 1-5, the correlation function exhibits two peaks: one
at zero frequency, and another at, roughly, at Larmor frequency $2f$. While
the nature of the latter peak is self-evident, the origin and the relevant
physics of the former, zero-frequency peak is less obvious. We emphasize that
the aforementioned behaviour occur for rather arbitrary long range couplings.

We argue (see also the main text) that the peak at $\omega=0$ is related to
the existence of random domains of spins initially having the same z-component
(i.e., $\pm 1$) and thus having relatively low relaxation times. Indeed, spins
in such domains would “see” sufficiently strong effective fields $h_{i}$ and,
therefore, would be less susceptible to transitions due to the transverse
field $f$. That is, when $h_{i}\gg f$, the individual spins in such domains
would tend to align along $z$ direction rather than rotate along total
magnetic field (of magnitude $2\sqrt{h_{i}^{2}+f^{2}}$; see more discussion
below). In Fig. 4 we present “screenshots” for the time evolution of spin’s
z-components showing the existence of such “metastable” domains that confirm
this argument.

Figure 3:  Numerical simulations: Spectrum of spin auto-correlation function
corresponding at $f=1.0$ for different couplings; Inset: Spectrum of the
effective field auto-correlation function $C_{h}(t)=(1/2)\langle
h_{i}(0)h_{i}(t)+h_{i}(t)h_{i}(0)\rangle$, with
$h_{i}=\sum_{j}J_{ij}\sigma^{z}_{i}$. Here $g_{h}(\omega)=\int d\omega
e^{i\omega t}\,C_{h}(t)$. For red curve $J_{ij}=1/|i-j|$, for blue curve
$J_{ij}=a/|i-j|^{2}$ and for the green curve $J_{ij}=\exp[(1-|i-j|)/b]$. The
values of coefficients $a$ and $b$ are chosen to satisfy condition
$\sum_{j}J_{ij}^{2}=W^{2}/4$, where $W=3.62$ for all three curves. Figure 4:
Numerical simulations: Same as in Fig. 1, but for $f=1.8$. Figure 5:
Numerical simulations: Same as Figs. 1 and 2, but for $f=2.0$. Figure 6:
Snapshots for simulated time-evolution of local spin $z$-polarization in a 1D
spin chain with long-range interactions for $f=0.5$ and $J_{ij}=1/|i-j|$. The
horizontal axis is the spin index. The relaxation time is strongly dependent
on the spin-fluctuation size. The relaxation time for large domains is
anomalously long and exceeds the entire accessible simulation time. Green
arrow points to such a stable spin domain.

At higher transverse fields $f$, between curves 6 and 7 (for $f=1.8$ and
$f=2.0$, respectively), the peak at $\omega=0$ disappears, which indicates
that a certain transition is taking place. This transition is related to the
existence of local, nearly constant in time, but random in space effective
fields (along z-direction) that persists below some critical value of the
transverse field $f$, and vanish above this critical value. Therefore we
believe that the system exhibits a dynamical phase transition between the two
regimes, with the root-mean-squared average of such fields being an order
parameter for such transition.

##  III Interpretations of the dynamic phase transition: Lipkin-Meshkov-Glick
model

The Lipkin-Meshkov-Glick (LMG) model LMG  is a simplified version of our
model, which can be understood in great detail. It is obtained from our model
after an assumption of all-to-all equal interactions of $N$ spins:

$H=(J/2)\left(\sum_{k=1}^{N}\sigma_{k}^{z}\right)^{2}+f\sum_{k=1}^{N}\sigma_{k}^{x}$
(15)

The spin projections

$S_{x,z}\equiv\frac{1}{2}\sum_{k=1}^{N}\sigma^{x,z}_{k}$ (16)

commute with the net spin operator $S^{2}$. Therefore, The Hamiltonian (15)
splits into invariant sectors with all possible values of $S=0,\ldots,N/2$,
which have different multiplicity. The Hamiltonian for the sector with spin
$S$, up to a constant, is given by the LMG model:

$H_{LMG}=2JS_{z}^{2}+2fS_{x}.$ (17)

In the large temperature limit, all sectors with different $S$ will contribute
to the noise power spectrum. However, for $N\gg 1$, the main contribution will
be produced by the sectors with typical spin fluctuation size

$\sqrt{\langle S_{z}^{2}\rangle}=\sqrt{N}/2\gg 1\,,$ (18)

where we have used Eq. (16) and the fact that different individual spins are
uncorrelated at the same time.

Let us introduce spherical coordinates $(\theta,\phi)$ that parameterize the
large-spin state, so that

$S_{z}=S\cos\theta,\quad S_{x}=S\sin\theta\cos\phi.$ (19)

Comparing Eq. (19) with Eq. (18), we notice that the value of $S$ can be
related with the total number of spins $N$ through the relation

$S^{2}\times\frac{1}{2}\simeq\frac{N}{4}\,,$ (20)

so that the typical values of $S$ are of the order of $\sqrt{N/2}$. That is,
in the limit of large $N$, the calculation of the correlation function for the
z-component of spin magnetization, $\langle S_{z}(0)S_{z}(t)\rangle$, formally
reduces to distribution of $S$’s according to Gaussian distribution $P(S)$
with variance $\sqrt{N/2}$, finding trajectories $S_{z}(t)$ for a fixed value
of $S$ (corresponding to the initial conditions $S_{x}(0)=S_{y}(0)=0$ and
$S_{z}(0)=S$), and then averaging the quantity $S_{z}(t,S)$ with respect to
$P(S)$.

Eq. (20) requires a comment. One could have argued that the coefficient $1/2$
in the left hand side of Eq. (20) should instead be $1/3$, which would
correspond to the angles $\theta$ and $\phi$ being uniformly distributed.
Indeed, under such assumption the average of $S_{z}^{2}$ is equal to $1/3$:

$\langle S_{z}^{2}\rangle=\frac{1}{4\pi}\int d\,\Omega
S^{2}\cos^{2}\theta=\frac{S^{2}}{3}\,,$

where $d\Omega=d\phi d\cos\theta$. This, however, is not the case in the
present calculation: Here we are interested in the correlation function of the
spin magnetization, and, therefore, the initial magnetization is constrained
by the condition $S_{x}=S_{y}=0$, with $S_{z}$ being random. Then, the
appropriate averaging corresponds to treating angle $\theta$ as a random
number, uniformly distributed in the interval between $0$ and $\pi$. This
gives factor $1/2$ in the left hand side of Eq. (20).

The classical spin dynamics is found by minimizing the action with an
effective Lagrangian LL1

$L=S\dot{\phi}\cos\theta-H_{LMG}(\theta,\phi),$

where in our case

$H_{LMG}(\theta,\phi)=2JS^{2}\cos^{2}\theta+2Sf\sin\theta\cos\phi.$

The equations of motion read:

${\dot{\theta}}=-2f\sin\phi\,,$ (21)
$\sin\theta{\dot{\phi}}=4JS\cos\theta\sin\theta-2f\cos\theta\cos\phi\,.$ (22)

After straightforward manipulation Eqs. (21,s22) can be cast in the following
form:

$\frac{d}{dt}(\sin\theta\cos\phi)=-4JS\cos\theta\sin\theta\sin\phi\,,$ (23)
$\frac{d}{dt}(\sin\theta\sin\phi)=4JS\cos\theta\sin\theta\cos\phi-2f\cos\theta\,.$
(24)

Using Eq. (21), we rewrite Eq. (23) as

$\frac{d}{dt}\big{(}f\sin\theta\cos\phi+JS\cos^{2}\theta\big{)}=0\,,$ (25)

which is nothing but the energy conservation, i.e.,
$H_{LMG}(\theta,\phi)=const$ along a trajectory. Using the initial conditions,
$\cos\theta=1$ and $\sin\theta=0$ at $t=0$, Eq. (25) can be integrated to give
the relation between $\phi$ and $\theta$:

$f\cos\phi=JS\sin\theta\,.$ (26)

Eliminating $\sin\theta$ and $\cos\theta$ in Eq. (24) and using Eqs. (21, 26),
Eq. (24) can be rewritten as:

$\frac{d^{2}\cos\theta}{dt^{2}}=4\cos\theta\Big{[}2(JS)^{2}(1-\cos^{2}\theta)-f^{2}\Big{]}\,.$
(27)

Equation (27) has the desired form: it is a closed equation for the
z-component of the spin magnetization $S_{z}=Scos\theta$.

Denoting $\cos\theta=z$, we see that Eq. (27) corresponds to that for particle
moving in external (1D) potential. The total energy of this “particle” can be
written as

$\frac{{\dot{z}}^{2}}{2}+V(z)=0\,.$ (28)

with

$V(z)=-2Az^{2}+Bz^{4}+2A-B\,.$ (29)

where

$A=2(JS)^{2}-f^{2}\,,\quad B=2(JS)^{2}\,.$ (30)

In Eqs. (28-30) we have again used the initial conditions $z(0)=\pm 1$ and
${\dot{z}}(0)=0$.

The behaviour of the “potential eenergy” $V(z)$ is presented schematically in
Figure 5. For $f<JS\sqrt{2}$, the potential is a double well, while for
$f>JS\sqrt{2}$ the potential energy has a single minimum at $z=0$. Therefore
point $f_{c}=JS\sqrt{2}$ corresponds to the phase transition discussed in the
main text. Indeed, we note that the parameter $W=2\sqrt{\sum_{j}J_{ij}^{2}}$
corresponds in the LMG model, with the typical spin size (18), to the
combination $W=2J\sqrt{N}\approx 2\sqrt{2}JS$, and hence to the critical
condition in LMG model given by

$f_{c}=W/2\,,$

which reproduces the result of our estimates of the critical point in the main
text.

Figure 7:  Effective “Ginzburg-Landau” potentials for different phases.

Another important point corresponds to the value of $f$ being equal to $JS$,
e.g. the red curve versus the blue curve in Fig. 5. Indeed, for $f$ smaller
than $JS$, the “particle” starting at $z=\pm 1$ with zero initial velocity,
cannot penetrate through the barrier and reach the other well. Such a
transition is impossible (red curve in Fig. 5). For the blue curve such
transition is possible. Therefore, for the dynamics in the potential given by
the red curve, the quantity $z(0)z(t)$ does not change its sign, which leads
to the divergence of the integral $\int dte^{i\omega t}z(0)z(t)$ at
$\omega=0$. This corresponds to results of the numerical simulation as well as
to the theoretical results of the main text: for small enough $f$, the
spectrum of the (z-component) spin fluctuations is divergent at $\omega=0$.
For the dynamics along the blue curve $z(t)$ changes its sign and therefore
the integral $\int dtz(0)z(t)$ is no longer divergent. This explains the
crossover behaviour exhibited by the fluctuation spectrum in Fig. 1 and Eq.
(14) of the main text.

Finally, we must note that in the limit $N\rightarrow\infty$ the fluctuation
spectrum is not going to exhibit a true phase transition (i.e., fully
vanishing spectral density at and around $\omega=0$ for $f>W/2$), but a
crossover. This is related to the fact that as $N\rightarrow\infty$ there are
trajectories corresponding to the values of $S$ that violate the condition
$f>\sqrt{2}JS$. For such trajectories the spectral weight at $\omega=0$ is
finite (or even infinite). The probabilities for such trajectories for large
enough $f$ are exponentially small, but finite from the formal standpoint
(recall that $S$ is a random variable obeying Gaussian distribution). This is
related to the infinite interaction range in the LMG model, which allows for
its integrability (e.g. “energy conservation” for the effective coordinate
$z$). We conjecture that the finite interaction range (between spins) would
effectively keep $S$ finite, thus leading to a true (and finite) critical
value for $f$, as observed in the numerical simulations; see, for example, the
inset in Figure 1 of the main text. A detailed study of of this issue,
however, lies out of the scope of this paper.

##  IV Noise Spectrum: time-dependent mean-field analysis

.

In this section we present the missing details of the calculation of the
spectrum of the spin auto-correlation function outlined in the main text.

Upon averaging over the rapidly varying fields $h_{f}(t)$, Eq. (10) of the
main text can be written as

${\dot{\sigma}^{z}_{s}}=-4f^{2}\int_{0}^{t}dt^{\prime}\sigma^{z}_{s}(t^{\prime})F(t-t^{\prime})\,,$
(31)

with

$F(\tau)=\cos(2h_{s}\tau)\exp\Big{[}-4\int\frac{d\omega}{2\pi}\,g_{f}(\omega)\frac{1-\cos(\omega\tau)}{\omega^{2}}\Big{]}.$
(32)

Note that the spectrum $g_{f}(\omega)$ can be associated with the Fourier
transform of the averaged Larmor part of the spectrum for the correlation
function $C^{0}(t)$ in Eq. (6) of the main text,

$g_{f}^{0}(\omega)=\frac{W^{2}}{4}\int dhP(h)\frac{\pi
f^{2}}{f^{2}+h^{2}}\Big{[}\delta(\omega-\omega_{h})+\delta(\omega+\omega_{h})\Big{]}\,,$
(33)

where $\omega_{h}=2\sqrt{f^{2}+h^{2}}$. The superscript in $g_{f}^{0}(\omega)$
indicates that as far as the actual frequency dependence of the spectrum is
concerned, this equation is a poor approximation: Indeed, due to the time
dependence of $h$, the lines at $\pm\omega_{h}$ (i.e., the delta-functions)
are broadened, which will clearly affect the shape of $g_{h}(\omega)$.
Furthermore, $g_{f}^{0}(\omega)$ does not contain contributions from the
cross-correlations between different spins. However, the variance of $h_{f}$,

$S_{f}=\int(d\omega/2\pi)g_{f}(\omega)\,,$

depends neither on the widths of the lines, nor on the spin cross-
correlations, - the latter is equal to zero since different spins are
uncorrelated at same times at infinite temperature. Therefore, while the shape
of the spectrum given by Eq. (33) is a too crude an approximation, the
estimate for its variance based on Eq. (33) is quite reasonable. This we
obtain that

$S_{f}=\frac{W^{2}}{4}\int
dhP(h)\frac{f^{2}}{f^{2}+h^{2}}=\frac{W^{2}}{4}-h_{0}^{2}\,,$ (34)

where $h_{0}$ is given by Eq. (8) of the main text. In the following we will
see that knowledge of the variance $S_{f}$ is sufficient for our purposes.

Solving Eq. (31) (supplemented by the initial condition $\sigma_{s}^{z}(0)=\pm
1$) by the method of Laplace transform, one readily obtains the power spectrum
for the slow spins,

$g_{s}(\omega,h_{s})=\frac{2\Gamma(\omega,h_{s})}{\big{[}\omega-\Omega(\omega,h_{s})\big{]}^{2}+\Gamma^{2}(\omega,h_{s})},$
(35)

where $\Gamma$ and $\Omega$ are the real and imaginary parts of
$2f^{2}\int_{0}^{\infty}d\tau e^{-i\omega\tau}F(\tau)$, respectively. This
equation enters in Eq. (12) of the main text.

Noise spectrum at small $\omega$ near the transition point:

In the limit of vanishing $h_{s}$ and small $\omega$, functions
$\Gamma(\omega,h_{s})$ and $\Omega(\omega,h_{s})$ in Eq. (35) can be
approximated as $\Gamma(0,0)$ and
$\omega\,\partial\Omega(0,0)/\partial\omega$, i.e..,

$\Gamma\simeq
2f^{2}\int_{0}^{\infty}d\tau\exp\Big{[}-4\int\frac{d\omega}{2\pi}\,g_{f}(\omega)\frac{1-\cos(\omega\tau)}{\omega^{2}}\Big{]}\,,$
(36) $\Omega(\omega)\simeq 2\omega
f^{2}\int_{0}^{\infty}d\tau\tau\exp\Big{[}-4\int\frac{d\omega}{2\pi}\,g_{f}(\omega)\frac{1-\cos(\omega\tau)}{\omega^{2}}\Big{]}\,.$
(37)

In evaluating these integrals we apply the saddle point method. While in the
vicinity of the phase transition (i.e. for $f\sim W/2$) the accuracy of the
method is poorly controlled (i.e., within $O(1)$, since the exponent does not
contain a large parameter), such an approximation is consistent with the
assumption that the distribution $P(h_{s})$ is Gaussian (e.g. Eq. (7) of the
main text). Indeed, the latter assumption is based on the central limit
theorem, which also utilizes the saddle point method reif .

The saddle point (in the integral over $\tau$ in Eq. (36)) is at
$\tau^{\ast}=0$. Expanding the function in the exponent in Eq. (36) around the
saddle point, the resulting integral can be written as

$\Gamma\simeq 2f^{2}\int_{0}^{\infty}d\tau\exp(-2S_{f}\,\tau^{2})\,,$ (38)

where the variance $S_{f}$ is given by Eq. (34), $S_{f}\simeq
W^{2}/4=f_{c}^{2}$. This yields

$\Gamma\simeq\sqrt{2\pi}f_{c}\,.$

In calculation of $\Omega(\omega,0)$ we apply the same approach. Then Eq. (37)
can be approximated as

$\Omega(\omega)\simeq 2\omega
f^{2}\int_{0}^{\infty}d\tau\tau\exp(-2S_{f}\,\tau^{2})\,,$ (39)

which yields

$\Omega(\omega)\simeq(f/f_{c})^{2}\omega\,.$

Substituting $\Gamma$ and $\Omega(\omega)$ into Eq. (35), we obtain the
expression for the spectrum of slow spins,

$g_{s}(\omega)=\frac{2\sqrt{2\pi}f_{c}}{\delta^{2}\omega^{2}+2\pi
f_{c}^{2}}\,,$ (40)

with $\delta$ given by Eq. (9) of the main text.

Substituting this expression into Eq. (12) of the main text and replacing the
factor $h_{s}^{2}/(f^{2}+h_{s}^{2})$ by $h_{s}^{2}/f_{c}^{2}$ (recall that
$h_{s}\ll f$ in the vicinity of the transition point), we obtain

$g(\omega)=\frac{2\sqrt{2\pi}f_{c}}{\delta^{2}\omega^{2}+2\pi f_{c}^{2}}\int
dh_{s}P(h_{s})\,\frac{h_{s}^{2}}{f_{c}^{2}}\,,$ (41)

with $P(h_{s})$ given by Eq. (7) of the main text. Performing the remaining
Gaussian integral, we obtain Eq. (13) of the main text.

Noise spectrum at small $\omega$ for $f\ll f_{c}$:

At sufficiently small $f$, the typical values of $h_{s}$ greatly exceed $f$
and therefore the relative fraction of slow spins, i e. the fraction in the
integrand in Eq. (12) in the main text, can be replaced by $1$. Furthermore,
in the limit of small $f$ and at small enough $\omega$ (see below), we can
neglect $\Omega(\omega,h_{s})$ in Eq. (35) and set $\omega=0$ in the
expression for $\Gamma(\omega,h_{s})$. Then, Eqs. (12) of the main text can be
simplified as

$g(\omega)=\int dh_{s}\,\frac{e^{-\frac{h_{s}^{2}}{2h_{0}^{2}}}}{\sqrt{2\pi
h_{0}^{2}}}\,\frac{2\Gamma(h_{s})}{\omega^{2}+\Gamma^{2}(h_{s})}\,,$ (42)

where $\Gamma(h_{s})$ is equal to

$2f^{2}\int_{0}^{\infty}d\tau\cos{(2h_{s}\tau)}\exp\Big{[}-4\int\frac{d\omega}{2\pi}\,g_{f}(\omega)\frac{1-\cos(\omega\tau)}{\omega^{2}}\Big{]}\,.$
(43)

Writing cosine as a sum of two exponents of complex arguments, we can evaluate
the integrals over $\tau$ by the saddle point method. Note that in the limit
$f\ll W$ the exponent in the integrand function has a large parameter ($W/f$)
and so the use of the saddle point method is well justified. The saddle points
are given by the following equation,

$\pm
ih_{s}=2\int\frac{d\omega}{2\pi}\,g_{f}(\omega)\,\frac{\sin(\omega\tau^{\ast})}{\omega}\,.$
(44)

Assuming that the argument of sine in the above equation is small compared to
$1$ (in magnitude), we find the saddle points at

$\tau^{\ast}\simeq\pm ih_{s}/(2S_{f})\,,$ (45)

with the variance $S_{f}$ given by Eq. (34). In what follows we will see that
the relevant values of $h_{s}$ are of the order of $\sqrt{Wf}$; see Eq. (46)
below. Also, from Eq. (34) and Eq. (8) of the main text, one can readily
verify that $S_{f}\simeq\sqrt{\pi/8}\,Wf$ for $f\ll W$, and so
$\tau^{\ast}\sim i/\sqrt{Wf}$. Furthermore, from Eq. (33) as well as from the
numerical noise spectra presented in Figure 1 of the main text, one can see
that $g_{f}(\omega)$ decays at frequencies of the order of a few $f$. Then the
quantity $\omega\tau^{\ast}$ in $\sin(\omega\tau^{\ast})$ in Eq. (44) is
effectively bounded by $\sqrt{f/W}\ll 1$, and, therefore, the approximate
solution for the saddle point in Eq. (45) is well justified. Then , after a
straightforward calculation we obtain

$\Gamma(h_{s})=f^{2}\sqrt{\frac{2\pi}{S_{f}}}\exp\Big{(}-\frac{h_{s}^{2}}{2S_{f}}\Big{)}\,.$
(46)

From Eq. (46 we see that the relevant values for $h_{s}$ are indeed bounded by
$\sqrt{S_{f}}\sim\sqrt{Wf}$ for $f\ll W$ and, therefore, our calculation for
$\Gamma(h_{s})$ is internally consistent.

Substituting Eq. (46) into Eq. (42) one obtains the desired expression for the
low frequency part of the fluctuation spectrum. It is, however, useful to
change the integration variable in the resulting integral from $h_{s}$ to
$\Gamma$. Then one obtains

$g(\omega)\sim\int_{0}^{\Gamma(0)}d\Gamma\,\Big{|}\frac{dh_{s}(\Gamma)}{d\ln\Gamma}\Big{|}\,\frac{\Gamma^{\alpha}}{\omega^{2}+\Gamma^{2}}\,,$
(47)

where we have dropped the proportionality constant to simplify notations. Here

$\alpha=\frac{W^{2}}{4h_{0}^{2}}-1\,.$ (48)

Note that for $f\ll W$, $\alpha\simeq\sqrt{2\pi}f/W$, and so $\alpha\ll 1$.

Since $h_{s}$ depends on $\Gamma$ as “square root of logarithm”, within
logarithmic accuracy one can assume that the factor $|dh_{s}/dln\Gamma|$ in
Eq. (47) is a constant. Therefore Eq. (47) reduces to Eq. (14) of the main
text.

## References

  * (1) Yi Huang, Tongcang Li, and Zhang-qi Yin,“Symmetry-breaking dynamics of the finite-size Lipkin-Meshkov-Glick model near ground state”, Phys. Rev. A 97, 012115 (2018).
  * (2) N. A. Sinitsyn, V. V. Dobrovitski, S. Urazhdin, A. Saxena, “Geometric Control Over the Motion of Magnetic Domain Walls”, Phys. Rev. B 77, 212405 (2008).
  * (3) F. Reif, “Fundamentals of statistical and thermal physics”, McGraw-Hill (1965).