# Sequential Kibble-Zurek dynamics in the anisotropic Ising model of the
Si(001) surface

G. Schaller  Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400,
01328 Dresden, Germany,  F. Queisser  Helmholtz-Zentrum Dresden-Rossendorf,
Bautzner Landstraße 400, 01328 Dresden, Germany,  S.P. Katoorani  Helmholtz-
Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany,
C. Brand  Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1,
47057 Duisburg, Germany,  C. Kohlfürst  Helmholtz-Zentrum Dresden-Rossendorf,
Bautzner Landstraße 400, 01328 Dresden, Germany,  M.R. Freeman  Department of
Physics, University of Alberta, 4-181 Centennial Center for Interdisciplinary
Science Edmonton, Alberta T6G 2E1, Canada,  A. Hucht  Fakultät für Physik,
Universität Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany,  P.
Kratzer  Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1,
47057 Duisburg, Germany,  B. Sothmann  Fakultät für Physik, Universität
Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany,  M. Horn-von Hoegen
Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, 47057
Duisburg, Germany,  R. Schützhold  Helmholtz-Zentrum Dresden-Rossendorf,
Bautzner Landstraße 400, 01328 Dresden, Germany,  Institut für Theoretische
Physik, Technische Universität Dresden, 01062 Dresden, Germany,

###### Abstract

As a simplified description of the non-equilibrium dynamics of buckled dimers
on the Si(001) surface, we consider the anisotropic 2D Ising model and study
the freezing of spatial correlations during a cooling quench across the
critical point. The dependence of the frozen correlation lengths $\xi_{\|}$
and $\xi_{\perp}$ on the cooling rate obtained numerically matches the Kibble-
Zurek scaling quite well. However, we also find that the ratio
$\xi_{\|}/\xi_{\perp}$ of their frozen values deviates significantly from the
ratio in equilibrium. Supported by analytical arguments, we explain this
difference by the fact that the deviation from equilibrium in the weakly
coupled direction occurs earlier than in the strongly coupled direction.

## Introduction

Von Neumann once [1] compared non-equilibrium theory to a theory of non-
elephants – indicating the richness and complexity of this field, which we are
just beginning to understand. In view of the diverging response time near the
critical point, continuous phase transitions are prototypical candidates for
observing non-equilibrium behavior [2, 3]. A prominent example is the Kibble
mechanism describing the formation of topological defects during symmetry-
breaking phase transitions in the early universe [4]. Later Zurek realized
that quite analogous effects should also occur in condensed matter such as
superfluid helium [5]. The Kibble-Zurek mechanism has been studied in numerous
theoretical (e.g. [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]) and experimental
investigations (e.g. [16, 17, 18, 19, 20, 21]). An important point is the
transition from adiabatic evolution to non-equilibrium behavior (such as
freezing) when approaching or traversing the critical point. Apart from the
original idea of creating topological defects, the general mechanism can also
be applied to the frozen domain structure in symmetry-breaking phase
transitions induced by the critical slowing down.

In the following, we shall study the anisotropic Ising model in two spatial
dimensions [22, 23, 24, 25, 26, 27] with special emphasis on possible
differences in the non-equilibrium behavior between the two directions. Apart
from advancing our fundamental understanding, these investigations are also
motivated by the fact that the buckling dynamics of dimers on the Si(001)
surface can be described by the anisotropic 2D Ising model [28, 29, 30, 31,
32, 33, 34, 35, 36, 37]. Here, we consider the transition from the $p(2\times
1)$ to the $c(4\times 2)$ reconstruction at a critical temperature $T_{\rm
crit}\approx 190~{}\rm K$. Since the (001) face of single-crystalline silicon
belongs to the most important surfaces both in technology and science, our
results will also be relevant in this regard. For example, the dependence of
the frozen domain structure on the cooling rate indicates how sufficiently
homogeneous Si(001) surfaces should be prepared.

## Kibble-Zurek scaling

Let us briefly recapitulate the main arguments leading to the standard Kibble-
Zurek scaling. We consider a symmetry-breaking second-order phase transition
at the critical temperature $T_{\rm crit}$. Approaching the critical point
$T_{\rm crit}$, the equilibrium correlation length $\xi^{\rm eq}$ obeys the
universal scaling behavior

$\displaystyle\xi^{\rm eq}\sim\left|\frac{T-T_{\rm crit}}{T_{\rm
crit}}\right|^{-\nu}\equiv|\tau|^{-\nu}$ (1)

with the universal critical exponent $\nu$ and the dimensionless reduced
temperature $\tau$. Similarly, the response or relaxation time $t_{\rm
relax}^{\rm eq}$ (in equilibrium) scales as

$\displaystyle t_{\rm relax}^{\rm eq}\sim|\tau|^{-z\nu}\sim(\xi^{\rm eq})^{z}$
(2)

with the dynamical critical exponent $z$. The divergence of $t_{\rm
relax}^{\rm eq}$ at the critical point is the hallmark of critical slowing
down.

Now the idea is to infer non-equilibrium properties from these equilibrium
values. Let us assume a linear time-dependence for the cooling protocol
$\tau(t)=-\eta t$ with the cooling rate $\eta$ such that, starting at the time
$t_{\rm in}<0$, the critical point is reached at $t=0$. Then we may estimate
the freezing time via $|t_{\rm freeze}|=t_{\rm relax}^{\rm eq}(\tau_{\rm
freeze})$ where $\tau_{\rm freeze}=\tau(t_{\rm freeze})$, after which the
system has no time to equilibrate any more. Insertion into Eq. (2) yields
$|t_{\rm freeze}|\sim\eta^{-z\nu/(z\nu+1)}$ and thus we obtain the frozen
correlation length from Eq. (1)

$\displaystyle\xi_{\rm freeze}=\xi^{\rm eq}(\tau_{\rm
freeze})\sim\eta^{-\nu/(z\nu+1)}\,,$ (3)

which is referred to as Kibble-Zurek scaling [9, 16, 20].

## Anisotropic Ising model

Now let us apply these ideas to the non-equilibrium dynamics of buckled dimers
on the Si(001) surface which form a rectangular lattice. If we describe the
tilt of the dimer at lattice site $i,j\in\mathbb{Z}$ to the left or the right
by the pseudo-spin variable $\sigma_{i,j}=+1$ or $\sigma_{i,j}=-1$,
respectively, the resulting energy landscape corresponds to the anisotropic
Ising model [28, 29, 30, 32, 33, 34, 35].

$\displaystyle E_{\mbox{\boldmath$\scriptstyle\sigma$}}$ $\displaystyle=$
$\displaystyle-
J_{x}\sum_{i,j}\sigma_{i,j}\sigma_{i+1,j}-J_{y}\sum_{i,j}\sigma_{i,j}\sigma_{i,j+1}$
(4) $\displaystyle-
J_{\times}\sum_{i,j}\sigma_{i,j}[\sigma_{i+1,j+1}+\sigma_{i+1,j-1}]$

with a strong anti-ferromagnetic coupling $J_{x}\approx-25~{}{\rm meV}$ in
$x$-direction (i.e., along the dimer rows) and weaker couplings $J_{y}\approx
3.2~{}{\rm meV}$ in $y$-direction (i.e., across the rows) as well as in
diagonal direction $J_{\times}\approx 2.0~{}{\rm meV}$ [36, 37]. The latter
two can be combined into an effective transversal coupling
$J_{\perp}=J_{y}-2J_{\times}\approx-0.8~{}{\rm meV}$. As a result, the Ising
model (4) favors anti-ferromagnetic order both in $x$\- and $y$-direction. For
convenience, we apply a checker-board transformation
$\sigma_{i,j}\to(-1)^{i+j}\sigma_{i,j}$ after which we have ferromagnetic
order since $J_{x}$ and $J_{y}$ change sign.

As explained above, the Kibble-Zurek scaling is derived from the equilibrium
properties. For the Ising model (4), they can be obtained by Onsager theory
[38]. In terms of the longitudinal $J_{\|}=J_{x}$ and transversal $J_{\perp}$
couplings, the critical temperature $T_{\rm crit}=1/(k_{\rm B}\beta_{\rm
crit})$ is determined by the relation $\sinh(2\beta_{\rm
crit}J_{\|})\sinh(2\beta_{\rm crit}J_{\perp})=1$. Thus, in the limit of strong
anisotropy $J_{\|}\gg J_{\perp}$, we obtain the hierarchy of scales
$J_{\|}\gg\beta_{\rm crit}^{-1}\gg J_{\perp}$. Approaching the critical point
from above, the correlation lengths $\xi_{\|}$ and $\xi_{\perp}$ in $x$\- and
$y$-direction (i.e., along the dimer rows and perpendicular to them) both obey
the scaling (1) with the critical exponent $\nu=1$, though with different pre-
factors [22, 24]. Thus, their ratio stays constant and is given by
$\xi_{\perp}/\xi_{\|}=\sinh(2\beta_{\rm crit}J_{\perp})\approx 2\beta_{\rm
crit}J_{\perp}$.

## Rate equations

For the 2D Ising model, the dynamical critical exponent reads
$z=2+\varepsilon$ where $\varepsilon$ is a small and positive number [39, 40,
41, 42, 43, 44, 45]. Thus, the exponent in the Kibble-Zurek scaling relation
(3) is roughly minus one third.

To test this relation, we have to study the non-equilibrium dynamics of the
Ising model (4). To this end, we employ rate equations for the probabilities
$P_{\mbox{\boldmath$\scriptstyle\sigma$}}$ of the configurations $\sigma$ in
the standard form

$\displaystyle\dot{P}_{\mbox{\boldmath$\scriptstyle\sigma$}}=\sum_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}[R_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\to\mbox{\boldmath$\scriptstyle\sigma$}}P_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}-R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}P_{\mbox{\boldmath$\scriptstyle\sigma$}}]\,.$
(5)

Neglecting correlated flips of two or more pseudo-spins (i.e., dimers), we use
single-flip transition rates

$\displaystyle
R_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\to\mbox{\boldmath$\scriptstyle\sigma$}}=\frac{\Gamma\exp\\{-\beta
E_{\rm
B}\\}}{\exp\\{\beta(E_{\mbox{\boldmath$\scriptstyle\sigma$}}-E_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}})\\}+1}\,.$
(6)

The “knocking” frequency $\Gamma\approx 10^{12}/{\rm s}$ and Arrhenius barrier
height $E_{\rm B}\approx 100~{}{\rm meV}$ are obtained from microscopic
considerations [46, 47, 37]. The Glauber factor in the denominator can also be
motivated by microscopic models, e.g., in the form of a reservoir of two-level
systems or via fermionic tunneling. It ensures that the rate is bounded
$R_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\to\mbox{\boldmath$\scriptstyle\sigma$}}<\Gamma\exp\\{-\beta
E_{\rm B}\\}$ and satisfies the detailed balance condition
$R_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\to\mbox{\boldmath$\scriptstyle\sigma$}}/R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}=\exp\\{\beta(E_{\mbox{\boldmath$\scriptstyle\sigma$}}-E_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}})\\}$
which enforces convergence to thermal equilibrium for constant parameters
$\Gamma$ and $\beta$ etc.

## Numerical simulations

Due to the exponential dimensionality of (5) for an $N_{x}\times N_{y}$ spin
lattice, we calculate trajectory solutions. For a given configuration
$\scriptstyle\sigma$, we propagate time by the stochastic waiting time
$\tau_{\mbox{\boldmath$\scriptstyle\sigma$}}$ found by numerically solving
$\ln(1-r)=-\sum_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}\int_{t}^{t+\tau_{\mbox{\boldmath$\scriptstyle\sigma$}}}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}(t^{\prime})dt^{\prime}$,
with uniformly distributed random number $r\in[0,1]$, and perform a jump to a
different state with the conditional probability [48, 49] given by
$P_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}=R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}/[\sum_{\mbox{\boldmath$\scriptstyle\sigma^{\prime\prime}$}\neq\mbox{\boldmath$\scriptstyle\sigma$}}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime\prime}$}}]$.
In the selection of jumps, we take advantage [50] of the fact that the
$N_{x}N_{y}$ different single-spin flip processes can be grouped into 45
classes with identical energy differences entering the rates (6). Eventually,
denoting the fast Fourier transformed spin lattice by
$\tilde{\sigma}_{k_{x}k_{y}}$, the correlation lengths $\xi_{\|}$ and
$\xi_{\perp}$ are then given by the inverse widths of the one-dimensional
lattices $\sum_{k_{y}}|\tilde{\sigma}_{k_{x}k_{y}}|^{2}$ and
$\sum_{k_{x}}|\tilde{\sigma}_{k_{x}k_{y}}|^{2}$, respectively. Averaging over
multiple trajectories (and the resulting $|\tilde{\sigma}_{k_{x}k_{y}}|^{2}$)
can be used to improve the statistics.

Figure 1: Plot of inverse correlation lengths versus temperature (or,
equivalently, time) for a cooling sweep for a $16000\times 2000$ lattice,
averaged over $100$ trajectories. The correlation length in weakly-coupled
direction (yellow curve) departs earlier than the other (brown curve) from the
equilibrium solutions (red squares and black circles). On top, we added snap-
shots of the time evolution of an example configuration at the respective
temperatures as an illustration.

In Fig. 1, we contrast the time-dependent averaged correlation lengths (solid
curves) with equilibrium versions (symbols) for a cooling sweep. Already at
temperatures above $T_{\rm crit}$, the correlation lengths depart from their
equilibrium limits, but furthermore we see that this happens earlier for the
weakly coupled direction.

The final (i.e., frozen) correlation lengths are depicted in Fig. 2 as a
function of the cooling rate $\eta\propto 1/t_{\rm prot}$. For very fast
sweeps, the system cannot follow and basically remains at the initial
equilibrium values. For intermediate-speed sweeps, we find that both final
correlation lengths follow a universal power-law increase, consistent with the
Kibble-Zurek exponent $\nu/(1+z\nu)\approx 1/3$ in Eq. (3), see the fitted
regions in Fig. 2. For very slow sweeps, we find that finite-size effects
start to play a role.

Again, since the frozen correlation lengths $\xi_{\|}^{\rm freeze}$ and
$\xi_{\perp}^{\rm freeze}$ both obey the scaling (3) in this intermediate
region, though with different pre-factors, their ratio $\xi_{\perp}^{\rm
freeze}/\xi_{\|}^{\rm freeze}$ is roughly constant – similar to the
equilibrium case discussed above where $\xi_{\perp}^{\rm eq}/\xi_{\|}^{\rm
eq}=\sinh(2\beta_{\rm crit}J_{\perp})\approx 2\beta_{\rm crit}J_{\perp}$.
However, we find that these two ratios are not the same, but differ by roughly
a factor of two.

Figure 2: Plot of the final (frozen) correlation lengths for different lattice
sizes and averaged over $100$ trajectories $[N_{x}\times N_{y}\times N_{\rm
trj}]$ for different cooling sweeps, where the system is cooled down from
$k_{\rm B}T_{\rm in}=25~{}\rm meV$ to $k_{\rm B}T_{\rm out}=10~{}\rm meV$ (as
in Fig. 1) in various time intervals, i.e., protocol times $t_{\rm prot}$. For
too fast protocols (left), the system can never follow, but Kibble-Zurek
scaling is recovered for intermediate protocol times. Finite-size effects are
visible for slow protocols.

## 1D Ising model

In order to understand the difference between the strongly and the weakly
coupled directions found above, let us first consider the limiting case
$J_{y}\to 0$ and $J_{\times}\to 0$ of the 2D Ising model (4) where each row
$j$ separately forms a 1D Ising model with $J=J_{x}$

$\displaystyle E_{\mbox{\boldmath$\scriptstyle\sigma$}}^{\rm
1D}=-J\sum_{i}\sigma_{i}\sigma_{i+1}\,.$ (7)

Assuming translational invariance, we may derive an exact evolution equation
for the correlator $c_{a}=\langle\sigma_{i}\sigma_{i+a}\rangle$ depending on
distance $a$. Furthermore, let us introduce the dimensionless conformal time
coordinate $\mathfrak{T}$ with adapted step size $d{\mathfrak{T}}/dt=\Gamma
e^{-\beta E_{\rm B}}$ such that

$\displaystyle\partial_{\mathfrak{T}}c_{a}=-2c_{a}+(c_{a+1}+c_{a-1})\tanh(2\beta
J)\,,$ (8)

with the boundary condition $c_{a=0}=1$.

Setting the left-hand side of Eq. (8) to zero yields the well-known
equilibrium solution $c_{a}=[\tanh(\beta J)]^{a}$. Note that the 1D Ising
model (7) does not have a critical point at finite temperature $T_{\rm
crit}>0$, instead the analogue of a critical point occurs at zero temperature
$T_{\rm crit}=0$ where $\xi$ diverges as $\xi\sim e^{2\beta J}$ [23].

In order to understand the non-equilibrium dynamics governed by Eq. (8), let
us consider the continuum limit where $c_{a+1}+c_{a-1}-2c_{a}$ becomes the
second spatial derivative such that we obtain a diffusion-dissipation equation
$\partial_{\mathfrak{T}}c={\mathfrak{D}}\partial_{x}^{2}c-\gamma c$. For large
temperatures, the diffusion coefficient is small
${\mathfrak{D}}\propto\tanh(2\beta J)\approx 2\beta J\ll 1$ and the damping
term $\gamma\approx 2$ dominates. For small temperatures, the damping rate
$\gamma=2-2\tanh(2\beta J)$ is suppressed as $4e^{-4\beta J}$ and the
diffusion term ${\mathfrak{D}}\propto\tanh(2\beta J)\approx 1$ dominates. In
analogy to the response time $t_{\rm relax}$ in Eq. (2), we may introduce a
response or relaxation time ${\mathfrak{T}}_{\rm relax}$ from the inverse
damping rate $1/\gamma$ which then scales as ${\mathfrak{T}}_{\rm relax}\sim
e^{4\beta J}$, i.e., ${\mathfrak{T}}_{\rm relax}\sim\xi^{2}$. Note, however,
that the diffusion coefficient stays finite even for ${\mathfrak{T}}_{\rm
relax}\to\infty$, i.e., diffusion is still possible.

## Freezing in 1D

Since analyzing the non-equilibrium dynamics by means of analytic solutions of
Eq. (8) is still quite involved, let us consider the weighted sum of
correlations ${\mathfrak{C}}=\sum_{a=1}^{\infty}ac_{a}$ which obeys the
simpler evolution equation

$\displaystyle\partial_{\mathfrak{T}}{\mathfrak{C}}=[2\tanh(2\beta
J)-2]{\mathfrak{C}}+\tanh(2\beta J)\,.$ (9)

In order to provide an explicit example and to study the analogue of critical
slowing down, let us assume the simple cooling protocol $\beta(t)=\kappa t$
(i.e., $T(t)\propto 1/t$) starting at infinite temperature at $t=0$ and
cooling down to zero temperature at $t\to\infty$. Then this infinite interval
of laboratory time $t\in(0,\infty)$ is mapped to a finite interval of
conformal time ${\mathfrak{T}}\in(-\Gamma/[\kappa E_{\rm
B}],0)=({\mathfrak{T}}_{\rm in},0)$.

Incidentally, for our values with $E_{\rm B}\approx 4J_{x}$, we may simplify
Eq. (9) even further. For low temperatures $\beta J\gg 1$, the source term
$\tanh(2\beta J)$ can be approximated by unity and the damping rate $\gamma$
behaves as $4e^{-4\beta J}$ which for $E_{\rm B}=4J$ becomes
$4{\mathfrak{T}}/{\mathfrak{T}}_{\rm in}$. As a result, Eq. (9) simplifies to
$\partial_{\mathfrak{T}}{\mathfrak{C}}=-4{\mathfrak{C}}{\mathfrak{T}}/{\mathfrak{T}}_{\rm
in}+1$. The solution to this equation can be given in terms of the error
function, but we may understand its behavior by means of general arguments. In
the limit of slow cooling rates considered here, we have $|{\mathfrak{T}}_{\rm
in}|\gg 1$. Then, starting with ${\mathfrak{C}}=0$ at
${\mathfrak{T}}={\mathfrak{T}}_{\rm in}$, the value of ${\mathfrak{C}}$
quickly approaches its instantaneous equilibrium value ${\mathfrak{C}}_{\rm
eq}={\mathfrak{T}}_{\rm in}/(4{\mathfrak{T}})$. However, once the response
time ${\mathfrak{T}}_{\rm relax}\sim{\mathfrak{T}}_{\rm in}/{\mathfrak{T}}$
becomes too short ${\mathfrak{T}}_{\rm relax}\sim|{\mathfrak{T}}|$, the system
cannot equilibrate anymore and thus the value of ${\mathfrak{C}}$ freezes in
at ${\mathfrak{T}}_{\rm freeze}\sim\sqrt{|{\mathfrak{T}}_{\rm in}|}$ to its
final value ${\mathfrak{C}}_{\rm freeze}\sim\sqrt{|{\mathfrak{T}}_{\rm in}|}$.
Using the asymptotic behavior of the error function, we may also determine the
pre-factor to ${\mathfrak{C}}_{\rm freeze}=\sqrt{|{\mathfrak{T}}_{\rm
in}|/(2\pi)}$. Furthermore, as $\mathfrak{C}$ scales with the square of the
correlation length $\xi_{\|}$, we obtain $\xi_{\|}^{\rm
freeze}\sim|{\mathfrak{T}}_{\rm in}|^{1/4}=(\Gamma/[\kappa E_{\rm B}])^{1/4}$,
which is the analogue to the Kibble-Zurek scaling (3) for this case.

## 2D Ising model

Now we can apply our findings to the anisotropic Ising model (4). Again
assuming translational invariance, the evolution equation for the correlations
$c_{a,b}=\langle\sigma_{i,j}\sigma_{i+a,j+b}\rangle$ reads

$\displaystyle\partial_{\mathfrak{T}}c_{a,b}$ $\displaystyle=$
$\displaystyle-2c_{a,b}+(c_{a+1,b}+c_{a-1,b})\tanh(2\beta J_{x})$ (10)
$\displaystyle+\,{\cal O}(J_{y},J_{\times})\,.$

Here, we used the separation of scales $J_{x},E_{\rm B}\gg J_{y},J_{\times}$
in order to keep the zeroth order in the first line, while all terms which are
suppressed by $J_{y}$ or $J_{\times}$ are in the second line. Note that these
terms actually contain factors of $\tanh(2\beta J_{y})$ and $\tanh(2\beta
J_{\times})$, i.e., they would not grow without bound even for extremely small
temperatures $\beta\sim 1/J_{y,\times}$. At such ultra-low temperatures, the
rates (6) are exponentially suppressed anyway by the barrier $E_{\rm B}$,
i.e., basically no flips would occur anymore.

Let us study the behavior resulting from Eq. (10). Since there are no source
terms in the first line of Eq. (10) for $b\neq 0$, correlations between rows
can only be created by the small terms in the second line, which also contain
the source term $c_{0,0}=1$. Without them, inter-row correlations $c_{a,b\neq
0}$ are only damped or diffused. Since the diffusion in $x$-direction (i.e.,
along the rows) is very fast, let us focus on the remaining slower evolution
by considering the total correlator between rows
${\mathfrak{R}}_{b}=\sum_{a=-\infty}^{\infty}c_{a,b}$

$\displaystyle\partial_{\mathfrak{T}}{\mathfrak{R}}_{b}=[2\tanh(2\beta
J_{x})-2]{\mathfrak{R}}_{b}+\,{\cal O}(J_{y},J_{\times})\,,$ (11)

in analogy to Eq. (9). In the following, we consider the total correlation
between neighboring rows $b=1$. In equilibrium, this quantity scales with
$\xi_{\|}$, provided that we assume a correlation length $\xi_{\perp}$ of
order unity or more. Thus, if $\xi_{\|}$ becomes very large, the system would
need to shuffle more and more correlations from one row to the next in order
to stay close to equilibrium. However, since this can only be done via the
small coupling terms $\,{\cal O}(J_{y},J_{\times})$, this bottleneck limits
the growth of ${\mathfrak{R}}_{b}$ and thus the system departs from
equilibrium at some point.

Since this departure time is determined by the smallness of $\,{\cal
O}(J_{y},J_{\times})$, the smaller these coupling are, the earlier this
departure occurs. Thus, for small $\,{\cal O}(J_{y},J_{\times})$, this
departure from equilibrium in the weakly coupled direction will occur earlier
than in the strongly coupled direction, where the system can still be close to
equilibrium.

## Conclusions

Figure 3: Low temperature STM image of a Si(001) surface taken at 5 K with
$U_{\rm bias}=1.3~{}\rm V$ and $I_{\rm tunnel}=1~{}\rm nA$. Field of view is
$24\times 16~{}\rm nm^{2}$. Areas with $c(4\times 2)$ reconstruction exhibit a
“honeycomb” pattern. Domain boundaries of the frozen domain structure can be
identified by a zig-zag chain of local $p(2\times 2)$ reconstruction. The two
dark spots are missing dimer vacancies and frizzy vertical lines correspond to
active phase boundary changes (”phasons”).

For a cooling quench of the anisotropic 2D Ising model across the critical
point, we studied the freezing of the spatial correlation lengths
$\xi_{\perp}$ and $\xi_{\|}$. We found that their dependence on the cooling
rate obtained numerically matches the Kibble-Zurek scaling quite well.
However, as a distinct signature of the non-equilibrium dynamics in the
presence of the anisotropy, we also observed that the ratio $\xi_{\perp}^{\rm
freeze}/\xi_{\|}^{\rm freeze}$ of their frozen values differs from the ratio
in equilibrium $\xi_{\perp}^{\rm eq}/\xi_{\|}^{\rm eq}$ by roughly a factor of
two. This difference contrasts with the simple picture which underlies the
Kibble-Zurek scaling and assumes that the whole system stays close to
equilibrium before the freezing time $t_{\rm freeze}$ and basically does not
evolve anymore afterwards. Instead, our combination of analytical and
numerical methods shows that the non-equilibrium dynamics in the two
directions is different and cannot be grasped by a single freezing time
$t_{\rm freeze}$.

Finally, let us discuss potential experimental evidence for the formation of a
frozen domain structure for the Si(001) dimerized surface. The surface
exhibits parallel rows of alternately buckled dimers which arrange in a
$c(4\times 2)$ reconstruction indicating the anti-phase correlation between
neighboring dimer rows [51, 52]. Fig. 3 shows a low-temperature scanning
tunneling microscopy image taken at 5 K after preparation of the Si(001)
surface through flash annealing and rapid cool-down to liquid nitrogen
temperatures $T<100~{}\rm K$. The cooling rate was on the order of 1-10 K/s.
Further experimental detaile 3s can be found elsewhere [47]. The STM image was
taken at constant current conditions with positive sample bias, i.e., in Fig.
3 filled orbitals of the Si atoms are displayed in bright. The dimer rows run
vertically from top to bottom. In each row the alternating buckling along the
dimer row can nicely be identified. The anti-phase correlation between
neighboring dimers cause the $c(4\times 2)$ reconstruction which becomes
apparent as “honeycomb” pattern. During the rapid cool-down the regime of
critical slowing down is reached for $T>T_{\rm crit}$, resulting in a frozen
domain structure which is apparent in Fig. 3. The domain boundaries can be
identified as one-dimensional “defects” separating ordered areas with a
$c(4\times 2)$ reconstruction. As expected from our findings described above
and from electron diffraction [36, 37], these ordered domains are extremely
elongated. However, albeit quite intriguing, these observations can only be
interpreted as a “smoking gun” for the non-equilibrium dynamics studied here
and further studies are required to settle this issue.

###### Acknowledgements.

## Acknowledgments

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) through the Collaborative Research Center SFB 1242 “Nonequilibrium
dynamics of condensed matter in the time domain” (Project-ID 278162697).

## References

  * [1] P Bak and M Paczuski.  Complexity, contingency, and criticality.  PNAS, 92:6689–6696,, 1995.
  * [2] S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Shahar.  Continuous quantum phase transitions.  Rev. Mod. Phys., 69:315–333, Jan 1997.
  * [3] Subir Sachdev.  Quantum Phase Transitions.  Cambridge University Press, 2011.
  * [4] T. W. B. Kibble.  Topology of cosmic domains and strings.  Journal of Physics A: Mathematical and General, 9(8):1387, aug 1976\.
  * [5] W. H. Zurek.  Cosmological experiments in superfluid helium?  Nature, 317:505–508, 1985.
  * [6] W. H. Zurek, U. Dorner, and P. Zoller.  Dynamics of a quantum phase transition.  Physical Review Letters, 95:105701, 2005.
  * [7] Bogdan Damski.  The simplest quantum model supporting the Kibble-Zurek mechanism of topological defect production: Landau-Zener transitions from a new perspective.  Physical Review Letters, 95:035701, 2005.
  * [8] L. Cincio, J. Dziarmaga, M. M. Rams, and W. H. Zurek.  Entropy of entanglement and correlations induced by a quench: Dynamics of a quantum phase transition in the quantum Ising model.  Physical Review A, 75:052321, 2007.
  * [9] Amit Dutta, R. R. P. Singh, and Uma Divakaran.  Quenching through Dirac and semi-Dirac points in optical lattices: Kibble-Zurek scaling for anisotropic quantum critical systems.  Europhysics Letters, 89(6):67001, apr 2010.
  * [10] Cheng-Wei Liu, Anatoli Polkovnikov, and Anders W. Sandvik.  Dynamic scaling at classical phase transitions approached through nonequilibrium quenching.  Phys. Rev. B, 89:054307, Feb 2014.
  * [11] Paul M. Chesler, Antonio M. García-García, and Hong Liu.  Defect formation beyond Kibble-Zurek mechanism and holography.  Phys. Rev. X, 5:021015, May 2015.
  * [12] Julian Sonner, Adolfo del Campo, and Wojciech H. Zurek.  Universal far-from-equilibrium dynamics of a holographic superconductor.  Nature Communications, 6(1):7406, 2015.
  * [13] Pietro Silvi, Giovanna Morigi, Tommaso Calarco, and Simone Montangero.  Crossover from classical to quantum kibble-zurek scaling.  Phys. Rev. Lett., 116:225701, Jun 2016.
  * [14] Balázs Dóra, Markus Heyl, and Roderich Moessner.  The Kibble-Zurek mechanism at exceptional points.  Nature Communications, 10:2254, 2019.
  * [15] C.J.O. Reichhardt, A. del Campo, and C. Reichhardt.  Kibble-Zurek mechanism for nonequilibrium phase transitions in driven systems with quenched disorder.  Commun Phys, 5:173, 2022.
  * [16] 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.
  * [17] Giacomo Lamporesi, Simone Donadello, Simone Serafini, Franco Dalfovo, and Gabriele Ferrari.  Spontaneous creation of Kibble–Zurek solitons in a Bose–Einstein condensate.  Nature Physics, 9:656–660, 2013.
  * [18] Sven Deutschländer, Patrick Dillmann, Georg Maret, and Peter Keim.  Kibble–Zurek mechanism in colloidal monolayers.  PNAS, 112(22):6925–6930, 2015.
  * [19] Jérôme Beugnon and Nir Navon.  Exploring the Kibble–Zurek mechanism with homogeneous Bose gases.  Journal of Physics B: Atomic, Molecular and Optical Physics, 50(2):022002, jan 2017.
  * [20] L.-Y. Qiu, H.-Y. Liang, Y.-B. Yang, H.-X. Yang, T. Tian, Y. Xu, and L.-M. Duan.  Observation of generalized Kibble-Zurek mechanism across a first-order quantum phase transition in a spinor condensate.  Science Advances, 6(21):eaba7292, 2020.
  * [21] Kai Du, Xiaochen Fang, Choongjae Won, Chandan De, Fei-Ting Huang, Wenqian Xu, Hoydoo You, Fernando J. Gómez-Ruiz, Adolfo del Campo, and Sang-Wook Cheong.  Kibble–Zurek mechanism of Ising domains.  Nature Physics, 19:1495–1501, 2023.
  * [22] Barry M. McCoy and Tai Tsun Wu.  The Two-Dimensional Ising Model.  Harvard University Press, Harvard, 1973.
  * [23] Rodney J. Baxter.  Exactly Solved Models in Statistical Mechanics.  Academic Press, London, 1989.
  * [24] Hendrik Hobrecht and Alfred Hucht.  Anisotropic scaling of the two-dimensional Ising model I: the torus.  SciPost Phys., 7:26, Aug 2019.
  * [25] Hendrik Hobrecht and Alfred Hucht.  Anisotropic scaling of the two-dimensional Ising model II: surfaces and boundary fields.  SciPost Phys., 8:32, Aug 2020.
  * [26] Alfred Hucht.  The square lattice Ising model on the rectangle iii: Hankel and Toeplitz determinants.  Journal of Physics A: Mathematical and Theoretical, 54(37):375201, sep 2021.
  * [27] Harold J. W. Zandvliet.  Phase diagram of the square 2D Ising lattice with nearest neighbor and next-nearest neighbor interactions.  Phase Transitions, 96(3-4):187–195, 2023.
  * [28] A. Saxena, E.T. Gawlinski, and J.D. Gunton.  Structural Phase Transitions on the Si(100) Surface.  Surf. Sci., 160(2):618–640, 1985.
  * [29] Masakazu Kubota and Yoshitada Murata.  Streak patterns in low-energy electron diffraction on Si(001).  Phys. Rev. B, 49:4810–4814, Feb 1994.
  * [30] Yoshitada Murata and Masakazu Kubota.  Order-disorder transition on Si(001).  Phase Transit., 53(2–4):125–141, 1995.
  * [31] K. Inoue, Y. Morikawa, K. Terakura, and M. Nakayama.  Order-disorder phase transition on the si(001) surface: Critical role of dimer defects.  Phys. Rev. B, 49:14774–14777, May 1994.
  * [32] Yoshimichi Nakamura, Hiroshi Kawai, and Masatoshi Nakayama.  Influence of defects on the order-disorder phase transition of a Si(001) surface.  Phys. Rev. B, 55:10549–10560, Apr 1997.
  * [33] Hiroshi Kawai, Yoshimichi Nakamura, and Masatoshi Nakayama.  Kinetic One-Dimensional Ising System on a Narrow Si(001)SB Terrace.  J. Phys. Soc. Jpn., 68(12):3936–3940, 1999.
  * [34] Devina Pillay, Brett Stewart, Chee Burm Shin, and Gyeong S. Hwang.  Revisit to the Ising model for order–disorder phase transition on Si(001).  Surf. Sci., 554(2):150–158, 2004.
  * [35] Hiroshi Kawai, Osamu Narikiyo, and Kensuke Matsufuji.  Structural Phase Transition between $c(4\times 2)$ and $p(2\times 2)$ Structures on Si(001) Surface under Observation by Scanning Tunneling Microscopy.  J. Phys. Soc. Jpn., 76(3):034602, 2007.
  * [36] Christian Brand, Alfred Hucht, Giriraj Jnawali, Jonas D. Fortmann, Björn Sothmann, Hamid Mehdipour, Peter Kratzer, Ralf Schützhold, and Michael Horn-von Hoegen.  Dimer coupling energies of the Si(001) surface.  Phys. Rev. Lett., 130:126203, Mar 2023.
  * [37] Christian Brand, Alfred Hucht, Hamid Mehdipour, Giriraj Jnawali, Jonas D. Fortmann, Mohammad Tajik, Rüdiger Hild, Björn Sothmann, Peter Kratzer, Ralf Schützhold, and Michael Horn von Hoegen.  Critical behavior of the dimerized Si(001) surface: A continuous order-disorder phase transition in the 2d Ising universality class.  arXiv:2310.10488, 2023.
  * [38] Lars Onsager.  Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition.  Phys. Rev., 65:117–149, Feb 1944.
  * [39] J. K. Williams.  Monte Carlo estimate of the dynamical critical exponent of the 2D kinetic Ising model.  Journal of Physics A: Mathematical and General, 18(1):49, jan 1985\.
  * [40] Nobuyasu Ito, Makoto Taiji, and Masuo Suzuki.  Accurate estimation of the dynamical critical exponent of the two-dimensional kinetic Ising model based on the Ising machine m-TIS.  Journal of the Physical Society of Japan, 56(12):4218–4220, 1987\.
  * [41] B. Dammann and J. D. Reger.  Dynamical critical exponent of the two-dimensional Ising model.  Europhysics Letters, 21(2):157, jan 1993.
  * [42] Christian Münkel, Dieter W. Heermann, Joan Adler, Misha Gofman, and Dietrich Stauffer.  The dynamical critical exponent of the two-, three- and five-dimensional kinetic Ising model.  Physica A: Statistical Mechanics and its Applications, 193(3):540–552, 1993.
  * [43] M. P. Nightingale and H. W. J. Blöte.  Dynamic exponent of the two-dimensional Ising model and Monte Carlo computation of the subdominant eigenvalue of the stochastic matrix.  Phys. Rev. Lett., 76:4548–4551, Jun 1996.
  * [44] B. C. S. Grandi and W. Figueiredo.  Dynamical critical exponent of a nonequilibrium Ising model.  Phys. Rev. E, 54:4722–4725, Nov 1996.
  * [45] M. P. Nightingale and H. W. J. Blöte.  Monte Carlo computation of correlation times of independent relaxation modes at criticality.  Phys. Rev. B, 62:1089–1101, Jul 2000.
  * [46] J. Da¸browski, E. Pehlke, and M. Scheffler.  Calculation of the surface stress anisotropy for the buckled si(001)(1×2) and p(2×2) surfaces.  Phys. Rev. B, 49:4790–4793, Feb 1994.
  * [47] Y. Pennec, M. Horn-von Hoegen, Xiaobin Zhu, D. C. Fortin, and M. R. Freeman.  Dynamics of an Ising Chain under Local Excitation: A Scanning Tunneling Microscopy Study of Si(100) Dimer Rows at 5 K.  Phys. Rev. Lett., 96:026102, Jan 2006.
  * [48] F. M. Bulnes, V. D. Pereyra, and J. L. Riccardo.  Collective surface diffusion: $n$-fold way kinetic Monte Carlo simulation.  Phys. Rev. E, 58:86–92, Jul 1998.
  * [49] Kurt Binder and Dieter W. Heermann.  Rejection-Free Monte Carlo, pages 179–190.  Springer International Publishing, Cham, 2019.
  * [50] Peter Kratzer.  Multiscale Simulation Methods in Molecular Sciences, volume 42 of NIC Series, chapter Monte Carlo and Kinetic Monte Carlo Methods – A Tutorial, pages 51–76.  John von Neumann Institute for Computing (NIC), Jülich Supercomputing Centre, Forschungszentrum Jülich, 52425 Jülich, Germany, 2009\.
  * [51] Robert A. Wolkow.  Direct observation of an increase in buckled dimers on Si(001) at low temperature.  Phys. Rev. Lett., 68(17):2636, 1992.
  * [52] R. G. Zhao and W. S. Yang.  Atomic structure of the Si(001) c(4 ×2) surface.  Phys. Rev. B, 33(10):6780, 1986.
  * [53] Roy J. Glauber.  Time-dependent statistics of the Ising model.  Journal of Mathematical Physics, 4(2):294–307, 12 1963.
  * [54] Yuji Sakai and Koji Hukushima.  Dynamics of one-dimensional Ising model without detailed balance condition.  Journal of the Physical Society of Japan, 82(6):064003, 2013.

##  Appendix A Appendix: 1D Ising model

A closer look at the Glauber-Arrhenius rates (6) shows that for the energy
differences arising for a single-spin flip at site $i$ in the 1D Ising model

$\displaystyle
E_{\mbox{\boldmath$\scriptstyle\sigma$}}-E_{\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}}=2J\sigma_{i}(\sigma_{i-1}+\sigma_{i+1})\,,$
(12)

the rates can be represented by a quadratic function of the spins near that
site [53, 54] with

$\displaystyle\frac{1}{e^{\beta(E_{\mbox{\boldmath$\scriptstyle\sigma$}}-E_{\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}})}+1}=\frac{1}{2}-\frac{\tanh(2\beta
J)}{4}\sigma_{i}(\sigma_{i-1}+\sigma_{i+1})\,,$ (13)

which with $\sigma_{i}^{2}=1$ enables to convert the rate equations (5) into
exact evolution equations for the non-equilibrium expectation values

$\displaystyle\partial_{t}\langle\sigma_{i}\rangle=\Gamma e^{-\beta E_{\rm
B}}\left[\tanh(2\beta
J)\frac{\langle\sigma_{i+1}\rangle+\langle\sigma_{i-1}\rangle}{2}-\langle\sigma_{i}\rangle\right]$
(14)

and similarly for the two-point functions

$\displaystyle\partial_{t}\langle\sigma_{i}\sigma_{j}\rangle=\Gamma e^{-\beta
E_{\rm B}}\left[-2\langle\sigma_{i}\sigma_{j}\rangle+\tanh(2\beta
J)\times\phantom{\frac{1}{2}}\right.$
$\displaystyle\left.\times\frac{\langle\sigma_{i+1}\sigma_{j}\rangle+\langle\sigma_{i-1}\sigma_{j}\rangle+\langle\sigma_{i}\sigma_{j+1}\rangle+\langle\sigma_{i}\sigma_{j-1}\rangle}{2}\right]\,.$
(15)

If we start in the symmetric phase $\langle\sigma_{i}\rangle=0$, the mean
field $\langle\sigma_{i}\rangle$ in Eq. (14) stays zero as expected. However,
due to $\langle\sigma_{i}\sigma_{i}\rangle=\langle\sigma_{i}^{2}\rangle=1$,
the evolution equation (A) for the correlations contains a source term and
thus correlations are generated even if they are absent initially.

##  Appendix B Appendix: 2D Ising model

For the anisotropic 2D Ising model – compare Eq. (4) with $J_{x}\to J_{\|}$,
$J_{y}\to J_{\perp}$, and $J_{\times}\to 0$ – the energy difference upon
flipping a single spin at site $(i,j)$ is

$\displaystyle
E_{\mbox{\boldmath$\scriptstyle\sigma$}}-E_{\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}}$
$\displaystyle=2J_{\|}\sigma_{ij}(\sigma_{i-1,j}+\sigma_{i+1,j})$
$\displaystyle\qquad+2J_{\perp}\sigma_{ij}(\sigma_{i,j-1}+\sigma_{i,j+1})\,,$
(16)

such that the Glauber rates (6) lead to a quartic representation

$\displaystyle\frac{1}{e^{\beta(E_{\mbox{\boldmath$\scriptstyle\sigma$}}-E_{\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}})}+1}$
$\displaystyle=\frac{1}{2}-f_{x}\sigma_{ij}(\sigma_{i-1,j}+\sigma_{i+1,j})$
$\displaystyle\qquad-f_{y}\sigma_{ij}(\sigma_{i,j-1}+\sigma_{i,j+1})$
$\displaystyle\qquad+g_{x}\sigma_{ij}(\sigma_{i-1,j}+\sigma_{i+1,j})\sigma_{i,j-1}\sigma_{i,j+1}$
$\displaystyle\qquad+g_{y}\sigma_{ij}\sigma_{i-1,j}\sigma_{i+1,j}(\sigma_{i,j-1}+\sigma_{i,j+1})\,,$
(17)

with the abbreviations

$\displaystyle f_{x}$ $\displaystyle=\frac{[1+2\cosh(4\beta
J_{\|})+\cosh(4\beta J_{\perp})]\tanh(2\beta J_{\|})}{16\cosh(2\beta
J_{\|}+2\beta J_{\perp})\cosh(2\beta J_{\|}-2\beta J_{\perp})}\,,$
$\displaystyle f_{y}$ $\displaystyle=\frac{[1+2\cosh(4\beta
J_{\perp})+\cosh(4\beta J_{\|})]\tanh(2\beta J_{\perp})}{16\cosh(2\beta
J_{\|}+2\beta J_{\perp})\cosh(2\beta J_{\|}-2\beta J_{\perp})}\,,$
$\displaystyle g_{x}$ $\displaystyle=\frac{\sinh^{2}(2\beta
J_{\perp})\tanh(2\beta J_{\|})}{4\cosh(4\beta J_{\|})+4\cosh(4\beta
J_{\perp})}\,,$ $\displaystyle g_{y}$ $\displaystyle=\frac{\sinh^{2}(2\beta
J_{\|})\tanh(2\beta J_{\perp})}{4\cosh(4\beta J_{\|})+4\cosh(4\beta
J_{\perp})}\,,$ (18)

which shows that will in general not be possible to obtain exact closed
evolution equations, as first-order terms will couple to third-order
correlators, second-order correlators couple to fourth-order terms and so on.
This hierarchy of correlators can however be closed at some order by mean-
field-type approximations, leading in general to nonlinear equations.

If we would keep only two-point correlators and neglect all three-point and
higher correlations, the above equation (10) would again become a diffusion-
dissipation equation (in the continuum limit), but with an anisotropic
diffusion kernel. The damping term $\gamma$ would then vanish at the critical
temperature $T_{\rm crit}$ of the 2D Ising model (4). However, while the
restriction to two-point correlations is exact for the 1D Ising model, it is
only an approximation for the 2D case. Moreover, this approximation is
expected to become quite inaccurate near the critical point.