Ordering Behavior of the Two-Dimensional Ising Spin Glass with Long-Range Correlated Disorder

The Ising spin glass consist of Ising spins $s_{\bm{m}}\in\{\pm 1\}$ on the sites $\bm{m}\in\Lambda$ of a two-dimensional lattice, i.e. $\Lambda\subset\mathbb{Z}^{2}$. In this study only square systems are considered, such that the spin glass has $L$ spins in each direction and $|\Lambda|=L^{2}$. The Hamiltonian is given by

where the sum runs over all pairs $\mathcal{M}$ of nearest-neighbor spin sites with periodic boundary conditions (BC) in one and free boundary conditions in the other direction. The bonds, $J_{\bm{m},\bm{n}}$, which represent the interaction between two spins, are random in strength and sign but remain constant over time. Hence, one speaks of a quenched disorder where the system is investigated under a fixed realization of the bonds, $\bm{J}$. Here, the bonds originate from a Gaussian random field, $J(\bm{x})$, which is a function of a continuous position $\bm{x}$, and has zero mean, $\langle J(\bm{x})\rangle=0$ and a covariance given by

$\bm{x},\bm{r}\in\mathbb{R}^{2}$, $a\geq 0$ and $r=\lVert\bm{r}\rVert$. The entries of $\bm{J}$ are given by, $J_{\bm{m},\bm{n}}=J((\bm{m}+\bm{n})/2)$, which ensures that the correlation decays in the same manner along both axes. The correlation exponent, $a$, is the only parameter to control the correlation. For $a=0$ one obtains the Ising model of a ferromagnet or antiferromagnet, respectively, depending on the bond realization. When $a\longrightarrow\infty$ the uncorrelated Ising spin glass model with Gaussian disorder is recovered.

To generate the correlated bonds numerically we utilized the Fourier Filtering Method (FFM) Makse et al. (1996); Ahrens and Hartmann (2011). The FFM is a procedure to create stationary correlated random numbers of previously independent random numbers. Because it is based on the convolution theorem it is possible to benefit from the computationally efficient Fast Fourier Transform Algorithm Cooley and Tukey (1965). For its implementation we relied upon the functions of the FFTW library, version 3.3.5 Frigo and Johnson (2005). Figure 1 shows the average bond correlation along the main axes of a system with $|\Lambda|=46^{2}$ spins calculated by the estimator,

Here $\langle...\rangle_{J}$ denotes the average with respect to the disorder. $\mathcal{M^{\prime}}(\bm{r})\subset\mathcal{M}$ contains those bonds $\{\bm{m},\bm{n}\}$ for which the bond $\{\bm{m}+\bm{r},\bm{n}+\bm{r}\}$ is also on the lattice. The fact that $\mathcal{M^{\prime}}(\bm{r})$ does not contain all the bonds $\{\bm{n},\bm{m}\}$ in $\mathcal{M}$ is due to the free boundary conditions in one direction. The point is that for vectors $\bm{r}$ which are not exactly parallel to the free boundary, the bond $\{\bm{m}+\bm{r},\bm{n}+\bm{r}\}$ does not exist for all bonds of $\mathcal{M}$. For the correlations shown in Figure 1 we only investigated the directions parallel and perpendicular to the free boundary.

The nature of the ground state (GS) of the two-dimensional Ising spin glass is an intriguing subject on its own Newman and Stein (2000); Arguin and Damron (2014). Furthermore, GS computations of finite systems are a well established tool Hartmann and Rieger (2002, 2004) to investigate the glassy behavior of the model in the zero-temperature limit Hartmann (2008). The GS is the spin configuration which minimizes Eq. (1) for a given realization of the bonds. In case of two-dimensional planar lattices there exist exact procedures to generate the GS with a polynomial worst-case running time. This is in contrast to the three or higher-dimensional variants which belong to the class of NP-hard problems Barahona (1982).

There is more than one approach to compute the GS of two-dimensional planar spin glasses, such as the algorithm of Bieche et al. Bieche et al. (1980) and that of Barahona et al. Barahona et al. (1982). The key idea of these algorithms is to create a mapping from the original problem defined on the underlying lattice graph of the spin glass onto a related graph which is constructed in such a manner that the GS can be extracted from a minimum-weight perfect matching, which is polynomially computable. In this work we applied an ansatz which includes Kasteleyn-city subgraphs into the mapping process and thus is more efficient Pardella and Liers (2008); Thomas and Middleton (2007) in terms of speed and memory usage than the above mentioned algorithms. This allowed us to investigate systems up to a linear system size of $L=724$ spins in each direction without needing excessive computational resources. For the computation of the minimum-weight perfect matching the Blossom ${IV}$ algorithm Cook and Rohe (1999) implemented by A. Rohe Rohe was utilized.

To study the spin glass at nonzero temperature we create domain wall (DW) excitations in the system. This was done by computing GSs under periodic and antiperiodic BCs also referred to as P-AP Carter et al. (2002). It works as follows. First, a spin glass with periodic BCs in one direction and free BCs in the other direction is generated with quenched disorder, $\bm{J}^{\text{(p)}}$. Then, its GS configuration, $\bm{s}_{\text{gs}}^{(\text{p})}$, is computed. Next, the periodic BCs are replaced by antiperiodic BCs by reversing the sign of one column of bonds parallel to the direction of periodicity, which leads to $\bm{J}^{\text{(ap)}}$. Afterwards, the new GS configuration, $\bm{s}_{\text{gs}}^{(\text{ap})}$, is calculated. The change of the BCs imposes a DW of minimal energy between the two spin configurations $\bm{s}_{\text{gs}}^{(\text{p})}$ and $\bm{s}_{\text{gs}}^{(\text{ap})}$. The energy of the DW is given by

The geometrical structure of a DW can be characterized by the number $\mathcal{D}_{L}$ of bonds which are included in the surface. To avoid including those bonds which are a direct result of the different BCs, the surface is defined to consist of those bonds which fulfill $J_{\bm{m},\bm{n}}^{(\text{p})}s_{\bm{m}}^{(\text{p})}s_{\bm{n}}^{(\text{p})}J_{\bm{m},\bm{n}}^{(\text{ap})}s_{\bm{m}}^{(\text{ap})}s_{\bm{n}}^{(\text{ap})}<0$, where $\bm{m},\bm{n}$ runs over all unordered pairs of nearest neighbor lattice sites Khoshbakht and Weigel (2018).

Next, we look at the influence of bond correlation on the properties of the previously discussed DW excitations. The absolute value of the DW energy is proportional to the coupling strength between block spins in the zero-temperature limit Bray and Moore (1987a). A stable order is possible if the average absolute value of the DW energy increases with the system size. In the uncorrelated case, $a\longrightarrow\infty$, one obtains a power law behavior Bray and Moore (1987a); Hartmann and Young (2001); Hartmann et al. (2002); Khoshbakht and Weigel (2018),

where $\theta$ is the stiffness exponent with its current best estimate $\theta=-0.2793(3)$ Khoshbakht and Weigel (2018). Since $\theta<0$ there is no stable spin-glass phase for temperatures larger than zero. Figure 3 shows the impact of bond correlation on the scaling of $\langle|\Delta E|\rangle_{J}$. For $a\geq 0.9$ one can still observe the pure power-law decay of the uncorrelated model on sufficiently long length scales. The inset demonstrates that, in this region, the stiffness exponent stays constant at a value equal to that of the uncorrelated model. For values of $a\leq 0.9$ the average $\langle|\Delta E|\rangle_{J}$ initially grows with system size, but starts to decrease for larger sizes, i.e. the curves exhibit a peak. The system size at the peak, $L^{*}$, shifts to larger system sizes on decreasing the correlation exponent $a$. This will be analyzed below.

First we show in figure 4 the results for the average DW energy $\langle\Delta E\rangle_{J}$, which behaves in a somewhat similar manner as the average of the absolute value.

If spin-glass order existed in this model for some value of $a$ one would observe an increase $\langle|\Delta E|\rangle_{J}$ in the limit $L\to\infty$, while at the same time $\langle\Delta E\rangle_{J}$ would remain at, or converge to, zero as a function of $L$. The latter indicates the absence of ferromagnetic order, but not necessarily the absence of spin glass order because for a spin glass the change of the boundary conditions from periodic to antiperiodic is symmetric, i.e., could either increase or decrease the GS energy, so the average would be zero even in the case of spin glass order. An increase of both $\langle|\Delta E|\rangle_{J}$ and $\langle\Delta E\rangle_{J}$ for small values of $L$ and $a$ corresponds to a ferro/antiferromagnetic ordering on local length scales, which is visible in Fig. 2 and will be discussed more below. Whether there is a true ordered phase for very small values of $a$ will be discussed next.

To track how the length scale of local order changes as a function of the correlation strength we measure the size at the peak, $L^{*}$, for both the domain-wall energy and the absolute value, as a function of the correlation exponent $a$. Numerically, $L^{*}$ was computed by fitting a parabola in the vicinity of the peaks of $\langle|\Delta E|\rangle_{J}$ and $\langle\Delta E\rangle_{J}$. As the fit in figure 5 demonstrates, the data for $L^{*}(a)$ is well described by an exponential function

A non-zero value of $a_{\text{crit}}$ would indicate that an ordered phase exists for $a<a_{\text{crit}}$. A true spin-glass phase would be possible if $a_{\text{c}rit}$ for $\langle\Delta E\rangle_{J}$ is smaller than $a_{\text{c}rit}$ for $\langle|\Delta E|\rangle_{J}$. We obtained values of $a_{\text{crit}}=0.13(0.10)$ for $\langle\Delta E\rangle_{J}$ and $0.10(0.17)$ for $\langle|\Delta E|\rangle_{J}$ with quality of the fit $Q=0.87$ and $Q=0.58$, respectively. Thus, a zero value for the critical correlation exponent parameter $a_{c}$ seems likely.

To consider this further, we set $a_{\text{crit}}$ to zero and obtain the values of the other fit parameters, which here are $A_{L}=3.63(22)$, $b_{L}=0.39(4)$, $c_{L}=2.10(6)$ ($Q=0.86$) for $\langle\Delta E\rangle_{J}$ and $A_{L}=3.5(4)$, $b_{L}=0.55(6)$, $c_{L}=1.85(8)$ ($Q=0.63$) for $\langle|\Delta E|\rangle_{J}$. Since the qualities of the fits remain almost identical in comparison to $a_{\text{crit}}\neq 0$ the data is considered to be consistent with $a_{\text{crit}}=0$, which would imply that there is no global order when $a>0$, either ferromagnetic or spin glass. The inset shows that the behavior of $L^{*}(a)$ is also compatible with $c_{L}=2$. Fits of this type, with $a_{\text{crit}}=0$ and $c_{L}=2$ fixed, have quality of the fit larger than $0.4$, which is reasonable. Note that an exponential dependence of the “breakup” length scale of ferromagnetic order as a function of disorder strength was also found in the two-dimensional random field Ising model, both at low temperatures Pytte and Fernandez (1985) and in the GS Seppälä et al. (1998).

The behavior of DW surfaces is regarded as one of the essential parameters which describe the properties of random systems Bray and Moore (1987b); Newman and Stein (2007). DWs separate spins in GS and reversed GS. Their surface is defined as those bonds which belong to the DW, and we denote the surface size by $\mathcal{D}$.

In the uncorrelated case, $a\longrightarrow\infty$, the average DW surface size exhibits a power law Kawashima (2000); Melchert and Hartmann (2007),

where $d_{s}$ is the fractal surface dimension, for which the best numerical estimate at present is $d_{s}=1.27319(9)$ Khoshbakht and Weigel (2018). When $a=0$ the system is a ferro/antiferromagnet and $\langle\mathcal{D}\rangle_{J}=L$, implying that $d_{s}=1$, i.e., the surface is not fractal here. Figure 6 shows the scaling of the DW surface size for different values of $a$. In general, it can be seen that the correlation decreases the number of bonds in the DW surface. Even for $a=0.001$, the data on this log-log plot shows a visible a visible deviation from the linear behavior which occurs for $a$ strictly zero.

To characterize this behavior we fitted pure power-laws to the data. For this purpose, we did not do a full fit to all data points, but used the sliding-window approach instead. Here, four values of $\langle\mathcal{D}\rangle_{J}$ which are adjacent in terms of system size were grouped together in one fit window, respectively. The independent variables of such a fit window were given by $(L_{1},L_{2},L_{3},L_{4})$ with $L_{i}<L_{i+1}$, $i=1,...,4$. The dependent variables corresponded to the data, i.e. $\langle\mathcal{D}_{L_{i}}\rangle_{J}$. The smallest independent variable of each fit window is denoted as $L_{1}$.

For each window, we fitted the power law to the data resulting in a value of $d_{s}$. The dependence of $d_{s}$ as a function of $L_{1}$ for different values of $a$ can be found in figure 7. In the uncorrelated case $d_{s}$ decreases as a function of $L_{1}$, whereas for strong correlations $d_{s}$ increases. In any case, for system sizes $L_{1}<32$ and small values $a\leq 1$ we observed some notable dependence of $d_{s}$ on the system size. This motivated us to include corrections to scaling the power law behavior in Eq. (7) by considering Khoshbakht and Weigel (2018); Hartmann and Moore (2003)

By using Eq. (8) the quality of the fit is larger than $0.79$ for all studied values of $a\leq 0.1$ with smallest linear system size of the fit $L_{\text{min}}^{(\text{fit})}=8$. For larger values of $a$, the pure linear fit was always fine, given our statistical accuracy. Figure 8 shows the resulting fractal surface dimension for all considered values of $a$. At $a\approx 0.4$ the fractal surface dimension starts to decline from $d_{s}=1.27319(9)$ Khoshbakht and Weigel (2018), in the uncorrelated case, to smaller values. This is also visible in Fig. 7. Thus, it appears that for small values of $a$ the fractal structure of the cluster changes, although there is no phase transition. Note that, in order to extract $d_{s}(a)$, we also performed fits of the form $d_{s}(L_{1};a)=d_{s}(a)+\kappa_{d}L_{1}^{\gamma_{d}}$ (not shown; $\kappa_{d}$ and $\gamma_{d}$ also depend on $a$) to the sliding window fractal dimensions shown in Fig. 7. The behavior of this extrapolated fractal dimension also exhibits the same notable decrease of $d_{s}$ for $a\leq 0.4$. Of course, it can not be ruled out that on sufficiently large length scales, i.e., for much larger sizes than are currently accessible, the fractal dimension of the uncorrelated model would be recovered again and thus $d_{s}=1.27319(9)$ Khoshbakht and Weigel (2018) for all values $a>0$.

In this section the effects of the bond correlations on spin correlations in the GS will be discussed. Note that at zero temperature there is no thermal disorder and, in our study which uses bonds with a continuous distribution, the ground state is non-degenerate apart from overall spin inversion.

The two-point spin correlation is given by $\langle s_{\bm{m}}^{(\text{gs})}s_{\bm{m}+\bm{r}}^{(\text{gs})}\rangle_{J}$, where $s_{\bm{m}}^{(\text{gs})}$ denotes a spin in the GS configuration. In the uncorrelated model $\langle s_{\bm{m}}^{(\text{gs})}s_{\bm{m}+\bm{r}}^{(\text{gs})}\rangle_{J}=0$ if $\bm{r}\neq 0$. The correlated bonds induce a local ferro/antiferromagnetic order into the GS. When $a=0$ the system is a ferro/antiferromagnet and the order is global. In the ferromagnetic case $s_{\bm{m}}^{(\text{gs})}s_{\bm{m}+\bm{r}}^{(\text{gs})}=1$ and in the antiferromagnetic case $s_{\bm{m}}^{(\text{gs})}s_{\bm{m}+\bm{r}}^{(\text{gs})}$ alternates between plus and minus one. Therefore, the GS spin correlation can be estimated by

and $\bm{r}=(r_{1},r_{2})\in\mathbb{Z}^{2}$. $\Lambda^{\prime}(\bm{r})\subset\Lambda$, similar to Eq. (3), contains those sites $\bm{m}$ for which $\bm{m}+\bm{r}$ is on the lattice, given the free BC in one direction. Note that this definition means that the correlation is measured for each site $\bm{m}\in\Lambda^{\prime}$ on one of the two sublattices of a checkerboard partition of the square lattice, such that the correlation is insensitive to whether the order is ferromagnetic or antiferromagnetic.

In figure 9 one can see the GS correlation for different values of $a$. The data is well described by a scaling form of type

From the perspective of ordering we are especially interested in the correlation length, $\xi_{\text{gs}}$, that provides the distance over which spins are notably correlated. The straightforward method to extract $\xi_{\text{gs}}$ is by fitting the function of Eq. (10) to the data. The problem with such an approach is that neither $\upsilon$ nor $\varphi$ are known. Also finite-size corrections reduce the match between scaling form and actual data. Thus, we used different approaches to obtain $\xi_{\text{g}s}$.

First, we perform a separate fit for each value of $a$ down to small correlations where the error bars start to exceed one quarter of the correlation value. We observed that for $a\geq 2$, the values of the exponents $\upsilon$ and $\varphi$ did not change much, while for smaller values of $a$ the exponents were a bit smaller. Therefore, we fixed the exponents to the (averaged) values seen for $a\geq 2$ and fitted only with respect to $\xi_{\text{gs}}$.

Second, we also performed a multifit, i.e., we fitted the correlation function simultaneously Hartmann (2015) with one value for $\upsilon$, one value for $\varphi$ and values of the correlation length at many values of $a$. The results of this multifit are also shown in Fig. 9. The obtained values for $\xi_{\text{gs}}$ for these two fitting approaches are shown in Fig. 10. As can be seen, the results from fixing the values of the exponents and from using the multifit approach do not differ much.

Before we discuss the behavior of $\xi_{\text{gs}}(a)$, we describe the third approach we have used to estimate the correlation length. Here, we used the integral estimator that was introduced in Ref. Belletti et al. (2008). It presupposes that for $r\leq\xi_{\text{gs}}$ the correlation function is dominated by a power law of type $r^{-\upsilon}$, whereas for $r>\xi_{\text{gs}}$ the correlation is negligible. As a consequence, the integral

is given by $I_{k}\propto(\xi_{\text{gs}})^{{k+1}-\upsilon}$ and thus

This result would be exact and independent of $k$ if the correlation actually was only a power law. For real correlations, the value of $k$ dictates which part of the correlation function contributes most to the integral. Because such a scaling approach is only valid when $r$ is much larger than the lattice constant a high value of $k$ reduces the systematic error of the method. On the other hand, large values of $k$ increase the statistical error of $I_{k}$. Following the recommendation of Belletti et al. (2009) we used $\xi_{\text{gs}}^{(1,2)}$ as a compromise. Note that since the statistical error of the measured correlation grows with distance $r$, one usually defines a cutoff distance up to which the data is directly used for the integral. Similar to Belletti et al. (2009) we specified this cutoff distance as the value of $r$ where $G_{\text{gs}}$ is smaller than three times its error. For values of $r$ larger than this cutoff distance we computed the integral up to the maximal length of $L$ from fits to Eq. (10). The start value of these fits were set to $r_{\text{min}}^{(\text{fit})}\geq 2$. Because it was observed that the correlation decays slightly differently along the directions with free and periodic BCs, the computations of the GS were also done with an independent set of simulations for full free BCs.

Next, the correlation length $\xi_{\text{gs}}^{(1,2)}$ was extracted from the average of the GS correlation, $G_{\text{gs}}$, along the main axes. The statistical error of $\xi_{\text{gs}}^{(1,2)}$ was estimated by bootstrapping Young (2015); Hartmann (2015) and the integrals according to Eq. (11) were computed by utilizing the midpoint integration rule. For small values of $a$, the contribution to $I_{k}$ from the integral beyond the cutoff gets increasingly large. For instance, when $a=1.15$ this contribution made up approximately $4\%$ of the value of $I_{2}$. Hence, to estimate the total error, we added an extra systematic error to the statistical error. This was done by analyzing the value of $\xi_{\text{gs}}^{(1,2)}$ for two other choices of the cutoff distance, i.e., being the distance where the error of the correlation function is two or four times larger than its estimate, respectively. The maximal deviation of these two values from the standard definition, which uses a magnitude of three error bars to define cutoff, was set to be the systematic error. Furthermore, because the statistical error of $I_{k}$ grows large for small values of the correlation exponent $a$ and the system has to be sufficiently large to neglect boundary effects, values of $a<0.8$ were not considered for this approach.

Figure 10 shows $\xi_{\text{gs}}^{(1,2)}$ as a function of the correlation exponent. In contrast to the data for the length scales $L^{\star}$, the plot shows that the behavior of none of the three correlation lengths that we have defined, $\xi_{\text{gs}}\ \text{(fit)},\ \xi_{\text{gs}}\ \text{(multifit)}$ and $\xi_{\text{gs}}^{(1,2)}$, is close to exponential. Nonetheless, for $a\to 0$ the correlation lengths may converge to this behavior. We found that the data for these three definitions of the correlation length could be well fit to a function of the type

The exponential in Eq. (13) is consistent with the previous results for the length scale $L^{*}$ and will dominate for $a\to 0$. The power-law part is chosen to describe the behavior for large values of $a$.

We have fitted the data to this function. For example, for the data obtained from the multifit we again obtained a small value of $a_{\text{crit}}=0.11(6)$. Thus once again we fixed $a_{\text{crit}}\equiv 0$ in order to determine the values of the remaining fit parameters, obtaining $A_{\xi}=4.0(3)$, $d_{\xi}=0.81(3)$, $b_{\xi}=1.58(8)$ and $c_{\xi}=1.19(4)$ with a good quality of the fit. Similar values, in particular for $c_{\xi}$, are found for the other two definitions of $\xi_{\text{gs}}$ that we used. This means that for large values of the correlation exponent $a$ the behavior of the correlation length as function of $a$ seems to be better described by a power law. Also, the behavior of the peak lengths and of the ferromagnetic correlation lengths differ a lot, but it is still possible that for very small values of $a$, an exponential dependence of $\xi_{\text{gs}}$ on $1/a^{2}$ would be recovered. Unfortunately, to investigate this issue much larger system sizes would have to be treated, well beyond current numerical capabilities.