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

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. We observe a crossover from 1D to 2D behavior. For rapid cooling, we find effectively 1D behavior in the strongly coupled direction, for which we provide an exact analytic solution of the non-equilibrium dynamics. For slower cooling rates, we start to see 2D behavior where our numerical simulations show an approach to the usual Kibble-Zurek scaling in 2D.

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, 16, 17, 18, 19, 20, 21, 22]) and experimental investigations (e.g. [23, 24, 25, 26, 27, 28, 29, 30]). 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 this Letter, we study the anisotropic Ising model in two spatial dimensions [31, 32, 33, 34, 35, 36, 37, 38] with special emphasis on possible differences 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 [39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. 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.

Figure 1: Low temperature scanning tunneling microscope (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, whereas domain boundaries can be identified by a zig-zag chain of local $p(2\times 2)$ reconstruction, with dimer buckling and resulting domains (vertically running rows) indicated in the inset. The two dark spots are frozen horizontal domain boundaries, while frizzy vertical lines correspond to active phase boundary changes, i.e., mobile ”phasons” [49].

Experimental observations

Let us start by presenting experimental evidence for the formation of frozen domain structures on 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 [50, 51]. Fig. 1 shows a low-temperature scanning tunneling microscopy (STM) image taken at 5 K after preparation of the Si(001) surface through flash annealing and rapid cool-down across $T_{\rm crit}$ to liquid nitrogen temperatures $T<100~{}\rm K$. The cooling rate was on the order of 1-10 K/s. Further experimental details can be found elsewhere [49]. The STM image was taken at constant current conditions with positive sample bias, i.e., in Fig. 1 filled orbitals of the Si atoms are displayed in bright. Along each vertically running row, the alternating buckling can be nicely 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 subsequent to sample preparation, the regime of critical slowing down [52, 53] is reached near $T_{\rm crit}$, where the system falls out of equilibrium, resulting in a frozen domain structure which is apparent in Fig. 1. The domain boundaries can be identified as one-dimensional “defects” separating ordered areas with a $c(4\times 2)$ reconstruction. As also confirmed by electron diffraction [47, 48], these ordered domains are extremely elongated, i.e., the frozen correlation length $\xi_{\|}$ along the dimer rows, i.e., in vertical direction, is much larger than the horizontal one $\xi_{\perp}$ across the rows.

Anisotropic Ising model

In order to understand the behavior of the correlation lengths mentioned above, we need to model 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 for the pseudo-spin configuration $\sigma$ corresponds to the anisotropic Ising model [39, 40, 41, 43, 44, 45, 46]

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

$$-J_{\rm d}\sum_{i,j}\sigma_{i,j}[\sigma_{i+1,j+1}+\sigma_{i+1,j-1}]\,.$$

Combining microscopic considerations with experimental data, a strong anti-ferromagnetic coupling along dimer rows $J_{x}\approx-25~{}{\rm meV}$ and weaker couplings across rows $J_{y}\approx 3.2~{}{\rm meV}$ and in diagonal direction $J_{\rm d}\approx 2.0~{}{\rm meV}$ have been found [47, 48]. The latter two can be combined into an effective transversal coupling $J_{\perp}=J_{y}-2J_{\rm d}\approx-0.8~{}{\rm meV}$. As a result, the actual surface favors anti-ferromagnetic order in both directions. For convenience however, we apply a checker-board transformation $\sigma_{i,j}\to(-1)^{i+j}\sigma_{i,j}$ after which both $J_{\|}=J_{x}$ and $J_{\perp}$ are positive and the transformed model favors ferromagnetic order.

Rate equations

We study the non-equilibrium dynamics of the Ising model (1) via standard rate equations for the probabilities $P_{\mbox{\boldmath$\scriptstyle\sigma$}}$ of the pseudo-spin configurations $\sigma$

$$\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$}}]\,.$$ (2)

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

$$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}\,.$$ (3)

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 [54, 49, 48]. The Glauber factor in the denominator can also be motivated by microscopic models for surface-bulk interactions, e.g., in the form of a reservoir of two-level systems [55] or via fermionic tunneling [56]. 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 evolution towards thermal equilibrium for constant parameters $\Gamma$ and $\beta$ etc.

1D Ising model

Now we are in the position to study the non-equilibrium dynamics of the Ising model (1) during a cooling quench. As the first step, we consider a very rapid cooling rate such that the system basically has no time for an exchange between the dimer rows, i.e., in the weakly coupled direction. Then, as also confirmed by numerical simulations (see Fig. 2), we may consider the limiting case $J_{y}\to 0$ and $J_{\rm d}\to 0$ of the 2D Ising model (1) such that each row $j$ separately forms a 1D Ising model with $J=J_{x}=J_{\|}$

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

where after the checker-board transformation $J>0$. 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$. With the dimensionless conformal time coordinate $d{\mathfrak{T}}/dt=\Gamma e^{-\beta(t)E_{\rm B}}$ we find (with the boundary condition $c_{a=0}=1$) [supplement]See supplemental material.

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

Setting the left-hand side to zero yields the well-known equilibrium solution $c_{a}=[\tanh(\beta J)]^{|a|}$. The 1D Ising model (4) 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 the correlation length $\xi$ diverges as $\xi\sim e^{2\beta J}$ [33].

In order to understand the non-equilibrium dynamics governed by Eq. (5), 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. Thus, 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. (5) is still quite involved, let us consider the weighted sum of correlations ${\mathfrak{C}}=\sum_{a=1}^{\infty}ac_{a}$ which obeys

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

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}}(t)={\mathfrak{T}}_{\rm in}e^{-\kappa E_{\rm B}t}$ with ${\mathfrak{T}}_{\rm in}=-\Gamma/[\kappa E_{\rm B}]$.

Incidentally, for our values with $E_{\rm B}\approx 4J_{x}$, we may simplify Eq. (6) 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 behaves as $\gamma\approx 4e^{-4\beta J}$ which for $E_{\rm B}=4J$ becomes $4{\mathfrak{T}}/{\mathfrak{T}}_{\rm in}$. As a result, Eq. (6) 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. For $\kappa\ll\Gamma/E_{\rm B}$, 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{\pi|{\mathfrak{T}}_{\rm in}|/8}$. 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 for this case.

Figure 2: Freezing of the weighted sum of correlations ${\mathfrak{C}}$ for cooling quenches $\beta(t)=\kappa t$. The black solid curve represents the analytic solution of Eq. (6) for the 1D Ising model (4) for $\kappa=10^{8}~{}\rm s^{-1}meV^{-1}$ and the dashed blue curve for $\kappa=10^{9}~{}\rm s^{-1}meV^{-1}$ while the symbols correspond to full numerical simulations of the 2D Ising model (1) on a $1600\times 200$ lattice averaged over $100$ trajectories, see below. Correlations between rows remain negligible for these parameters (not shown). The dotted red curve depicts the equilibrium value in 1D and the horizontal dashed-dotted lines represent the approximation ${\mathfrak{C}}_{\rm freeze}\approx\sqrt{\pi|{\mathfrak{T}}_{\rm in}|/8}$ for the two cases with ${\mathfrak{T}}_{\rm in}=-10$ and ${\mathfrak{T}}_{\rm in}=-100$, respectively. The vertical dotted line depicts the critical temperature in 2D.

2D Ising model

As the next step, let us study slower cooling rates, where we start to see two-dimensional behavior. For simplicity, we consider the two-dimensional Ising model in terms of the effective transversal coupling strength $J_{\perp}$ introduced above

$$E_{\mbox{\boldmath$\scriptstyle\sigma$}}^{\rm 2D}=-J_{\|}\sum_{i,j}\sigma_{i,j}\sigma_{i+1,j}-J_{\perp}\sum_{i,j}\sigma_{i,j}\sigma_{i,j+1}\,.$$ (7)

The equilibrium properties of this model can be obtained by Onsager theory [58]. Its 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 $|T-T_{\rm crit}|^{-\nu}$ with the critical exponent $\nu=1$, though with different pre-factors [32, 35]. In particular, their ratio stays constant and is given by $\xi_{\perp}^{\rm eq}/\xi_{\|}^{\rm eq}=\sinh(2\beta_{\rm crit}J_{\perp})\approx 2\beta_{\rm crit}J_{\perp}\approx 0.1$.

Now considering the correlator $c_{a,b}=\langle\sigma_{i,j}\sigma_{i+a,j+b}\rangle$ and its evolution equation analogous to Eq. (5), we find that the terms on the right-hand side containing $J_{\perp}$ do also involve four-point correlators such as $\langle\sigma_{i,j}\sigma_{k-1,\ell}\sigma_{k+1,\ell}\sigma_{k,\ell\pm 1}\rangle$. To close this set of equations approximately, we may apply a perturbative expansion scheme and neglect terms of order $J_{\perp}^{2}$ [supplement]. Then, since the four-point correlators have already small pre-factor $\sim J_{\perp}$ out front, we may approximate them by their zeroth order $\langle\sigma_{i,j}\sigma_{k-1,\ell}\sigma_{k+1,\ell}\sigma_{k,\ell\pm 1}\rangle\approx\langle\sigma_{i,j}\sigma_{k,\ell\pm 1}\rangle\langle\sigma_{k-1,\ell}\sigma_{k+1,\ell}\rangle$ where the short-range correlations $\langle\sigma_{k-1,\ell}\sigma_{k+1,\ell}\rangle$ within a row $\ell$ can be approximated by their equilibrium value $\langle\sigma_{k-1,\ell}\sigma_{k+1,\ell}\rangle\approx\tanh^{2}(\beta J_{\|})$ in 1D (as discussed above). After that, we obtain an approximate diffusion-dissipation equation in 2D

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

$$+\beta J_{\perp}^{\rm eff}(c_{a,b+1}+c_{a,b-1})+\,{\cal O}(J_{\perp}^{2})\,,$$

but with a strongly reduced transversal coupling strength $J_{\perp}^{\rm eff}=2J_{\perp}/\cosh(2\beta J_{\|})$. As an intuitive picture, the pre-existent strong in-row correlation requires many spin flips for cross-row equilibration and thus renders this process very hard. In the continuum limit, the ratio between the effective diffusion coupling strengths ${\mathfrak{D}}_{\|}\propto\tanh(2\beta J_{\|})$ and ${\mathfrak{D}}_{\perp}\propto\beta J_{\perp}^{\rm eff}$ in the two directions in Eq. (8) determines the ratio of the equilibrium correlation lengths which is consistent with the results above $\xi_{\perp}^{\rm eq}/\xi_{\|}^{\rm eq}=\sqrt{{\mathfrak{D}}_{\perp}/{\mathfrak{D}}_{\|}}=\sqrt{2\beta J_{\perp}/\sinh(2\beta J_{\|})}\approx 2\beta_{\rm crit}J_{\perp}$ close to the critical point.

In the non-equilibrium case, Eq. (8) allows us to understand the crossover from 1D to 2D. Inserting our values, we find $\beta_{\rm crit}J_{\perp}^{\rm eff}|_{\rm crit}\approx 0.01$, which explains why a time interval of $|{\mathfrak{T}}_{\rm in}|=10$ (dashed blue curve in Fig. 2) is too short to generate correlations between the rows. Only for slower sweep rates, such as $|{\mathfrak{T}}_{\rm in}|=100$ (solid black curve in Fig. 2), these correlations start to build up (not shown), even though they are still very weak.

Numerical simulations

Finally, let us compare the analytical approximation schemes discussed above to a full numerical simulation of Eqns. (1-3) with time-dependent $\beta(t)$ [supplement]. Due to the exponential dimensionality of (2) for an $N_{x}\times N_{y}$ spin lattice, we calculate trajectory solutions. For a given configuration $\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 [59, 60] 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}$}}]$ at time $t+\tau_{\mbox{\boldmath$\scriptstyle\sigma$}}$. In the selection of jumps, we use [61] that the $N_{x}N_{y}$ different single-spin flip processes can be grouped into 45 classes with identical energy differences entering the rates (3). 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 3: Plot of inverse correlation lengths versus temperature (or, equivalently, time) for a linear cooling sweep $T(t)=-\eta t$ with $\eta=174\times 10^{6}~{}\rm K/s$ (i.e., $\Gamma t_{\rm prot}=10^{6}$ in Fig. 4) for a $16000\times 2000$ lattice, averaged over $100$ trajectories. The correlation length in weakly-coupled direction (orange curve) departs earlier than the other (brown curve) from the equilibrium solutions (red squares and black circles). On top, we added snap-shots ($67\times 50$ pseudo-spins) of the time evolution of an example configuration at the respective temperatures as an illustration.

In Fig. 3, we contrast the time-dependent averaged correlation lengths (solid curves) with equilibrium versions (symbols) for a cooling sweep (starting after an equilibration phase). 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. 4 as a function of the cooling rate $\eta\propto 1/t_{\rm prot}$. For extremely fast sweeps $\Gamma t_{\rm prot}<10^{3}$, the system cannot follow and basically remains at the initial equilibrium values. However, already for $\Gamma t_{\rm prot}=\,{\cal O}(10^{4})$, we start observing a growth of $\xi_{\|}$ while $\xi_{\perp}$ is still negligible, i.e., effectively 1D behavior. For intermediate-speed sweeps, such as $\Gamma t_{\rm prot}=\,{\cal O}(10^{6})$, 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$, see the fitted regions in Fig. 4. For very slow sweeps, we find that finite-size effects start to play a role.

Figure 4: 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. 3) 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.

Conclusions

As a model for the buckling dynamics of dimers on Si(001) surfaces, we consider the anisotropic Ising model in 2D and study the freezing of spatial correlations (which determine the final domain structure) when cooling through the critical point. Depending on the cooling rate, we find a crossover from effectively 1D behavior for rapid cooling to 2D behavior for slower cooling rates. Similar crossovers can also be observed for symmetry-breaking fields [62]. Exact analytic solution of the 1D case with $\beta(t)=\kappa t$ yields a scaling of the frozen correlation lengths $\xi^{\rm freeze}\propto\kappa^{1/4}$. For the 2D case with $T(t)=-\eta t$ at intermediate cooling rates $\eta$, we numerically found a scaling with $\eta^{\nu/(1+z\nu)}\approx\eta^{1/3}$ which is consistent with the Kibble-Zurek scaling in 2D. These theoretical predictions strongly encourage experimental investigations of the frozen domain structure. The variation of the cooling rate in this anisotropic 2D Ising system will provide a well-controlled experimental approach to Kibble-Zurek dynamics in condensed matter.

Appendix A Appendix: 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

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

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

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

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 [52, 53].

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. (10) yields $|t_{\rm freeze}|\sim\eta^{-z\nu/(z\nu+1)}$ and thus we obtain the frozen correlation length from Eq. (9)

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

which is referred to as Kibble-Zurek scaling [9, 23, 28].

Appendix B Appendix: 1D Ising model

A closer look at the Glauber-Arrhenius rates (3) shows that for the energy differences arising for a single-spin flip at site $i$ in the 1D Ising model ($\sigma_{i}^{\prime}=-\sigma_{i}$)

$$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 [63, 64] with

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

For the expectation value of a spin at site $i$ this implies via Eq. (2)

$$\partial_{t}\langle\sigma_{i}\rangle =\sum\limits_{\mbox{\boldmath$\scriptstyle\sigma$}}\sigma_{i}\dot{P}_{\mbox{\boldmath$\scriptstyle\sigma$}}=\sum\limits_{\mbox{\boldmath$\scriptstyle\sigma$},\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}}(\sigma_{i}^{\prime}-\sigma_{i})R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}}P_{\mbox{\boldmath$\scriptstyle\sigma$}}$$ (14)

$$=-2\sum\limits_{\mbox{\boldmath$\scriptstyle\sigma$}}\sigma_{i}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to F_{i}\mbox{\boldmath$\scriptstyle\sigma$}}P_{\mbox{\boldmath$\scriptstyle\sigma$}}$$

$$=-\Gamma e^{-\beta E_{\rm B}}\times$$

$$\qquad\times\sum\limits_{\mbox{\boldmath$\scriptstyle\sigma$}}\left[\sigma_{i}-\frac{\tanh(2\beta J)}{2}(\sigma_{i-1}+\sigma_{i+1})\right]P_{\mbox{\boldmath$\scriptstyle\sigma$}}$$

$$=\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]\,,$$

where $F_{i}$ denotes the flip operator for the $i$th spin and we used furthermore $\sigma_{i}^{2}=1$. With this, we have obtained exact evolution equations for the non-equilibrium expectation values.

Similarly, we can write the two-point functions as $\partial_{t}\langle\sigma_{i}\sigma_{j}\rangle=-2\sum_{\mbox{\boldmath$\scriptstyle\sigma$}}\sigma_{i}\sigma_{j}(R_{\mbox{\boldmath$\scriptstyle\sigma$}\to F_{i}\mbox{\boldmath$\scriptstyle\sigma$}}+R_{\mbox{\boldmath$\scriptstyle\sigma$}\to F_{j}\mbox{\boldmath$\scriptstyle\sigma$}})P_{\mbox{\boldmath$\scriptstyle\sigma$}}$, which eventually leads to

$$\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.$$

$$\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 (B) for the correlations contains a source term and thus correlations are generated even if they are absent initially.

Taking $j=i+a$ and assuming translational invariance, we obtain from (B) an equation for the correlator

$$\partial_{t}c_{a}=\Gamma e^{-\beta E_{\rm B}}\left[-2c_{a}+\tanh(2\beta J)(c_{a-1}+c_{a+1})\right]\,,$$ (16)

which upon transformation to the conformal time coordinate becomes Eq. (5) in the main text.

Appendix C Appendix: 2D Ising model

For the anisotropic 2D Ising model – compare Eq. (1) with $J_{x}\to J_{\|}$, $J_{y}\to J_{\perp}$, and $J_{\rm d}\to 0$ – the energy difference upon flipping a single spin at site $(i,j)$ (i.e., $\sigma_{i,j}^{\prime}=-\sigma_{i,j}$) is

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

$$\quad-2J_{\perp}\sigma_{i,j}(\sigma_{i,j-1}+\sigma_{i,j+1})\,,$$ (17)

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

$$\frac{1}{e^{\beta(E_{\mbox{\boldmath$\scriptstyle\sigma$}^{\prime}}-E_{\mbox{\boldmath$\scriptstyle\sigma$}})}+1} =\frac{1}{2}-f_{\|}\sigma_{i,j}(\sigma_{i-1,j}+\sigma_{i+1,j})$$

$$\quad-f_{\perp}\sigma_{i,j}(\sigma_{i,j-1}+\sigma_{i,j+1})$$

$$\quad+g_{\|}\sigma_{i,j}(\sigma_{i-1,j}+\sigma_{i+1,j})\sigma_{i,j-1}\sigma_{i,j+1}$$

$$\quad+g_{\perp}\sigma_{i,j}\sigma_{i-1,j}\sigma_{i+1,j}(\sigma_{i,j-1}+\sigma_{i,j+1})\,,$$ (18)

with the abbreviations

$$f_{\|} =\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})}\,,$$

$$f_{\perp} =\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})}\,,$$

$$g_{\|} =\frac{\sinh^{2}(2\beta J_{\perp})\tanh(2\beta J_{\|})}{4\cosh(4\beta J_{\|})+4\cosh(4\beta J_{\perp})}\,,$$

$$g_{\perp} =\frac{\sinh^{2}(2\beta J_{\|})\tanh(2\beta J_{\perp})}{4\cosh(4\beta J_{\|})+4\cosh(4\beta J_{\perp})}\,,$$ (19)

which shows that it will in general not be possible to obtain exact closed evolution equations, as e.g. first-order terms couple to three-point functions ($F_{ij}$ implements a flip of spin at $ij$)

$$\partial_{t}\langle\sigma_{i,j}\rangle =-2\sum_{\mbox{\boldmath$\scriptstyle\sigma$}}\sigma_{i,j}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to F_{ij}\mbox{\boldmath$\scriptstyle\sigma$}}P_{\mbox{\boldmath$\scriptstyle\sigma$}}$$

$$=\Gamma e^{-\beta E_{\rm B}}\Big{[}-\langle\sigma_{i,j}\rangle$$

$$\qquad+2f_{\|}\langle\sigma_{i-1,j}+\sigma_{i+1,j}\rangle+2f_{\perp}\langle\sigma_{i,j-1}+\sigma_{i,j+1}\rangle$$

$$\qquad-2g_{\|}\langle(\sigma_{i-1,j}+\sigma_{i+1,j})\sigma_{i,j-1}\sigma_{i,j+1}\rangle$$

$$\qquad-2g_{\perp}\langle\sigma_{i-1,j}\sigma_{i+1,j}(\sigma_{i,j-1}+\sigma_{i,j+1})\rangle\Big{]}\,,$$ (20)

and analogously two-point functions couple to four-point functions

$$\partial_{t}\langle\sigma_{i,j}\sigma_{k,\ell}\rangle =-2\sum_{\mbox{\boldmath$\scriptstyle\sigma$}}\sigma_{i,j}\sigma_{k,\ell}(R_{\mbox{\boldmath$\scriptstyle\sigma$}\to F_{ij}\mbox{\boldmath$\scriptstyle\sigma$}}+R_{\mbox{\boldmath$\scriptstyle\sigma$}\to F_{k\ell}\mbox{\boldmath$\scriptstyle\sigma$}})P_{\mbox{\boldmath$\scriptstyle\sigma$}}$$

$$=\Gamma e^{-\beta E_{\rm B}}\Big{[}-2\langle\sigma_{i,j}\sigma_{k,\ell}\rangle$$

$$\quad+2f_{\|}\langle\sigma_{k,\ell}(\sigma_{i-1,j}+\sigma_{i+1,j})\rangle$$

$$\quad+2f_{\|}\langle\sigma_{i,j}(\sigma_{k-1,\ell}+\sigma_{k+1,\ell})\rangle$$

$$\quad+2f_{\perp}\langle\sigma_{k,\ell}(\sigma_{i,j-1}+\sigma_{i,j+1})\rangle$$

$$\quad+2f_{\perp}\langle\sigma_{i,j}(\sigma_{k,\ell-1}+\sigma_{k,\ell+1})\rangle$$

$$\quad-2g_{\|}\langle\sigma_{k,\ell}(\sigma_{i-1,j}+\sigma_{i+1,j})\sigma_{i,j-1}\sigma_{i,j+1}\rangle$$

$$\quad-2g_{\|}\langle\sigma_{i,j}(\sigma_{k-1,\ell}+\sigma_{k+1,\ell})\sigma_{k,\ell-1}\sigma_{k,\ell+1}\rangle$$

$$\quad-2g_{\perp}\langle\sigma_{k,\ell}\sigma_{i-1,j}\sigma_{i+1,j}(\sigma_{i,j-1}+\sigma_{i,j+1})\rangle$$

$$\quad-2g_{\perp}\langle\sigma_{i,j}\sigma_{k-1,\ell}\sigma_{k+1,\ell}(\sigma_{k,\ell-1}+\sigma_{k,\ell+1})\rangle\Big{]}\,,$$ (21)

and so on. The above hierarchy can be closed at some order by mean-field-type approximations, leading in general to nonlinear equations, which however is expected to become inaccurate near the critical point.

If we would keep only two-point correlators and neglect all three-point and higher correlations, the above equation 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 a mean-field approximation to the critical temperature $T_{\rm crit}$ of the 2D Ising model (1). 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.

In the strongly anisotropic scenario relevant for our considerations, noting that $g_{\perp}=\,{\cal O}\{J_{\perp}\}$ and considering the strong correlations along rows, we can approximate terms like $\langle\sigma_{k,\ell}\sigma_{i-1,j}\sigma_{i+1,j}\sigma_{i,j\pm 1}\rangle\approx\langle\sigma_{k,\ell}\sigma_{i,j\pm 1}\rangle\langle\sigma_{i-1,j}\sigma_{i+1,j}\rangle\approx\langle\sigma_{k,\ell}\sigma_{i,j\pm 1}\rangle\tanh^{2}(\beta J_{\|})$. Taking $k=i+a$ and $\ell=j+b$, and exploiting translational invariance, we then obtain for the correlator (further using that $g_{\|}=\,{\cal O}\{J_{\perp}^{2}\}$),

$$\partial_{t}c_{a,b} =\Gamma e^{-\beta E_{\rm B}}\Big{[}-2c_{a,b}+4f_{\|}(c_{a-1,b}+c_{a+1,b})$$

$$\qquad+4[f_{\perp}-g_{\perp}\tanh^{2}(\beta J_{\|})](c_{a,b-1}+c_{a,b+1})\Big{]}$$

$$\qquad+\,{\cal O}\{J_{\perp}^{2}\}\,.$$ (22)

Furthermore inserting $4f_{\|}=\tanh(2\beta J_{\|})+\,{\cal O}\{J_{\perp}^{2}\}$ and $4[f_{\perp}-g_{\perp}\tanh^{2}(\beta J_{\|})]=2\beta J_{\perp}/\cosh(2\beta J_{\|})+\,{\cal O}\{J_{\perp}^{2}\}$ and employing the mapping to the conformal time as before, we recover Eq. (8) in the main text. An approximation to the critical temperature is then obtained from solving $1=4f_{\|}+4[f_{\perp}-g_{\perp}\tanh^{2}(\beta_{\rm crit}J_{\|})]$.

If instead we had neglected all four-point correlators completely, we would obtain the above equation with $g_{\perp}\to 0$, which with $4f_{\perp}=\beta J_{\perp}[1+\cosh^{-2}(2\beta J_{\|})]+\,{\cal O}\{J_{\perp}^{2}\}$ would lead to a larger diffusion constant along the weakly-coupled direction. Also the equation for the approximate critical point $1=4f_{\|}+4f_{\perp}$ would yield a less accurate result in this case.

Appendix D Appendix: Numerical Methods

Stochastic propagation

To solve a rate equation of the form

$$\dot{\mbox{\boldmath$P$}}={\cal L}(t)\mbox{\boldmath$P$}\,,$$ (23)

with probability vector $P$ and time-dependent rate matrix ${\cal L}(t)={\cal L}_{0}(t)+{\cal L}_{1}(t)$ (with diagonal part ${\cal L}_{0}(t)$ and an off-diagonal part ${\cal L}_{1}(t)$) via trajectories, we employ a Dyson-type series for the full propagator

$${\cal P}(t,0) ={\cal P}_{0}(t,0)+\int_{0}^{t}dt_{1}{\cal P}_{0}(t,t_{1}){\cal L}_{1}(t_{1}){\cal P}_{0}(t_{1},0)$$

$$\qquad+\int_{0}^{t}dt_{2}\int_{0}^{t_{2}}dt_{1}{\cal P}_{0}(t,t_{2}){\cal L}_{1}(t_{2}){\cal P}_{0}(t_{2},t_{1})\times$$

$$\qquad\qquad\times{\cal L}_{1}(t_{1}){\cal P}_{0}(t_{1},0)+\ldots$$ (24)

As ${\cal L}_{0}(t)$ is a diagonal matrix, it commutes with itself at different times, and we can compute the associated propagator simply via

$${\cal P}_{0}(t_{2},t_{1})=\exp\left\{\int_{t_{1}}^{t_{2}}{\cal L}_{0}(t^{\prime})dt^{\prime}\right\}\,.$$ (25)

Therefore, assuming that at time $t$ we are in state $\sigma$, the associated probability to remain in that state up to time $t+\tau$ is

$$P_{\rm stay}^{\mbox{\boldmath$\scriptstyle\sigma$}}(t,\tau)=\exp\left\{-\sum_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\neq\mbox{\boldmath$\scriptstyle\sigma$}}\int_{t}^{t+\tau}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}(t^{\prime})dt^{\prime}\right\}\,.$$ (26)

The waiting time distribution for the first jump to a different state is then

$$w^{\mbox{\boldmath$\scriptstyle\sigma$}}(t,\tau)=\frac{d}{d\tau}\left[1-P_{\rm stay}^{\mbox{\boldmath$\scriptstyle\sigma$}}(t,\tau)\right]\,.$$ (27)

To obtain accordingly distributed waiting times from uniformly distributed random numbers $r\in[0,1]$ we thus equate those with the cumulative probability of jumping $r=1-\exp\left\{-\sum_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\neq\mbox{\boldmath$\scriptstyle\sigma$}}\int_{t}^{t+\tau}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}(t^{\prime})dt^{\prime}\right\}$, such that we solve the equation

$$\ln(1-r)=-\sum_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\neq\mbox{\boldmath$\scriptstyle\sigma$}}\int_{t}^{t+\tau}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}(t^{\prime})dt^{\prime}$$ (28)

for the waiting time $\tau$. This can be achieved by Newton-Raphson iteration, where a suitable initial guess can be obtained by assuming that rates are roughly constant over the waiting time

$$\tau_{0}=-\frac{\ln(1-r)}{\sum\limits_{\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}\neq\mbox{\boldmath$\scriptstyle\sigma$}}R_{\mbox{\boldmath$\scriptstyle\sigma$}\to\mbox{\boldmath$\scriptstyle\sigma^{\prime}$}}(t)}\,.$$ (29)

Extraction of correlation lengths

Denoting the (fast) Fourier-transformed grid for an $N_{x}\times N_{y}$ lattice as

$$\tilde{\sigma}_{k\ell} =\sum_{ij}\sigma_{ij}e^{-2\pi{\rm i}(i-1)(k-1)/N_{x}}e^{-2\pi{\rm i}(j-1)(\ell-1)/N_{y}}\,,$$ (30)

it follows that the quantities

$$\sum_{k}{\left|\tilde{\sigma}_{k\ell}\right|}^{2} =N_{x}\sum_{q}\left[\sum_{ij}\sigma_{ij}\sigma_{i,j+q}\right]e^{+2\pi{\rm i}q(\ell-1)/N_{y}}\,,$$

$$\sum_{\ell}{\left|\tilde{\sigma}_{k\ell}\right|}^{2} =N_{y}\sum_{q}\left[\sum_{ij}\sigma_{ij}\sigma_{i+q,j}\right]e^{+2\pi{\rm i}q(k-1)/N_{x}}\,,$$ (31)

are Fourier transforms of the convolutions $\sum_{ij}\sigma_{ij}\sigma_{i,j+q}$ and $\sum_{ij}\sigma_{ij}\sigma_{i+q,j}$. By ensemble-averaging (runs can and have been performed in parallel), we thus have

$$\sum_{k}{\langle{\left|\tilde{\sigma}_{k\ell}\right|}^{2}\rangle} =N_{x}\sum_{q}\left[\sum_{ij}{\langle\sigma_{ij}\sigma_{i,j+q}\rangle}\right]e^{+2\pi{\rm i}q(\ell-1)/N_{y}}\,,$$

$$\sum_{\ell}{\langle{\left|\tilde{\sigma}_{k\ell}\right|}^{2}\rangle} =N_{y}\sum_{q}\left[\sum_{ij}{\langle\sigma_{ij}\sigma_{i+q,j}\rangle}\right]e^{+2\pi{\rm i}q(k-1)/N_{x}}\,.$$ (32)

Accordingly, when the expectation values decay like ${\langle\sigma_{ij}\sigma_{i+q,j}\rangle}\propto e^{-{\left|q\right|}/\xi_{\|}}$ and ${\langle\sigma_{ij}\sigma_{i,j+q}\rangle}\propto e^{-{\left|q\right|}/\xi_{\perp}}$ with correlation lengths $\xi_{\|}$ and $\xi_{\perp}$, it follows that we can obtain the latter from the inverse widths of $\sum_{\ell}{\left|\tilde{\sigma}_{k\ell}\right|}^{2}$ and $\sum_{k}{\left|\tilde{\sigma}_{k\ell}\right|}^{2}$, respectively. The subtraction of an average magnetization has no impact on our simulations. For the simulations of fast quenches of the form $\beta(t)=\kappa t$ starting at infinite temperatures, it should be noted that the initially flat Fourier transforms $\sum_{\ell}{\langle{\left|\tilde{\sigma}_{k\ell}\right|}^{2}\rangle}_{0}\to N_{x}N_{y}^{2}$ and $\sum_{k}{\langle{\left|\tilde{\sigma}_{k\ell}\right|}^{2}\rangle}_{0}\to N_{x}^{2}N_{y}$ will for the quenches shown in Fig. 2 in the main text still have considerable tails, such that for the numerically stable computation of the weighted correlators we used a cutoff that was given by the 10-fold correlation length

$$\mathfrak{C} =\sum_{a=1}^{\lceil 10\xi_{\|}\rceil}\frac{a}{N_{x}N_{y}}\sum_{ij}{\langle\sigma_{i,j}\sigma_{i+a,j}\rangle}$$

$$=\sum_{a=1}^{\lceil 10\xi_{\|}\rceil}\frac{a}{N_{x}^{2}N_{y}^{2}}\sum_{k\ell}{\langle{\left|\tilde{\sigma}_{k\ell}\right|}^{2}\rangle}e^{-2\pi{\rm i}a(k-1)/N_{x}}\,,$$ (33)

and analogously for the perpendicular direction. Variation of the cutoff had little impact on our results as long as it remained significantly smaller than the lattice size.