Numerical study of the roughness of domain walls in the 2-dimensional random bond Ising model using a weighed-loop algorithm

TL;DR

This paper introduces a weighed-loop algorithm for efficiently finding low-energy states in the 2D random bond Ising model, demonstrating improved performance over existing methods and analyzing domain wall roughness across various conditions.

Contribution

The paper presents a novel weighed-loop algorithm applicable to any weighed graph network, optimized for exploring rough energy landscapes in the 2D random bond Ising model.

Findings

The weighed-loop algorithm outperforms Niedermayer's algorithm at low temperatures and high disorder.

High-temperature roughness behavior matches low-temperature expansion predictions.

Temperature, disorder, and system size significantly influence domain wall roughness.

Abstract

We introduce a weighed-loop algorithm that is applicable to any weighed graph network. It is designed to prefer a route of energetically unfavourable bonds in the lattice that can then be flipped without changing the structure inside and outside the enclosed loop. Due to this property there are effectively no energy barriers thus making this algorithm very suitable for finding low energy states in very rough energy landscapes. We apply this algorithm to the random bond Ising model with domain walls and show that the weighed-loop algorithm can outperform Niedermayer's algorithm for low enough temperatures and high enough disorder. We consolidate the high-temperature behaviour of the roughness of a domain-wall with a low-temperature expansion presented in this paper and show agreement with results from our simulations. The effects of temperature, disorder, and system size on the roughness of domain-walls is also investigated.

Tables

Table 1. Table 1 : For both the Niedermayer as well as the weighed-loop algorithm a 𝑎 a has been measured for a variety of parameter values. At each parameter point the greatest value, and thus the algorithm with faster decorrelation, has a light grey background to visualize the trend that the weighed-loop algorithm starts to outperform Niedermayer’s algorithm for low enough temperatures and high enough disorder. The best values for a 𝑎 a come from averages of the best fits to the first 120 120 120 data points obtained from simulations similar as shown in Fig.  3 .

$\Delta$	Niedermayer ($a$)	weighed-loop ($a$)
$H=32$, $\beta =1$
$0.1$	$6.5\left(2\right)\times {10}^{-4}$	$1.35\left(4\right)\times {10}^{-5}$
$0.3$	$1.8\left(6\right)\times {10}^{-5}$	$1.42\left(8\right)\times {10}^{-5}$
$0.5$	$2.1\left(3\right)\times {10}^{-5}$	$1.3\left(2\right)\times {10}^{-5}$
$H=64$, $\beta =1$
$0.1$	$9.4\left(9\right)\times {10}^{-7}$	$8.7\left(6\right)\times {10}^{-7}$
$0.3$	$1.4\left(1\right)\times {10}^{-6}$	$8.9\left(8\right)\times {10}^{-7}$
$0.5$	$1.0\left(1\right)\times {10}^{-6}$	$1.2\left(2\right)\times {10}^{-6}$
$H=128$, $\beta =1$
$0.1$	$2.0\left(2\right)\times {10}^{-7}$	$4.5\left(4\right)\times {10}^{-8}$
$0.3$	$5.3\left(6\right)\times {10}^{-7}$	$8.1\left(7\right)\times {10}^{-8}$
$0.5$	$6.4\left(5\right)\times {10}^{-7}$	$1.1\left(2\right)\times {10}^{-7}$
$H=32$, $\beta =2$
$0.1$	$7.9\left(4\right)\times {10}^{-5}$	$9.4\left(6\right)\times {10}^{-6}$
$0.3$	$1.1\left(2\right)\times {10}^{-6}$	$1.1\left(1\right)\times {10}^{-5}$
$0.5$	$0.5\left(1\right)\times {10}^{-6}$	$1.3\left(3\right)\times {10}^{-5}$
$H=64$, $\beta =2$
$0.1$	$2.6\left(3\right)\times {10}^{-7}$	$8.4\left(9\right)\times {10}^{-7}$
$0.3$	$4.7\left(5\right)\times {10}^{-7}$	$1.2\left(2\right)\times {10}^{-6}$
$0.5$	$1.9\left(2\right)\times {10}^{-7}$	$1.8\left(5\right)\times {10}^{-6}$
$H=128$, $\beta =2$
$0.1$	$0.4\left(2\right)\times {10}^{-6}$	$6.1\left(8\right)\times {10}^{-8}$
$0.3$	$1.9\left(3\right)\times {10}^{-7}$	$1.3\left(2\right)\times {10}^{-7}$
$0.5$	$8.7\left(9\right)\times {10}^{-8}$	$2.5\left(3\right)\times {10}^{-7}$

Equations

$$\mathcal{H}(\sigma )=-\frac{1}{2}\sum _{i\in \Lambda }\sum _{j\in {\mathcal{N}}_i}{J}_{ij}{\sigma }_i{\sigma }_j,$$

$$\mathcal{F}(p)$$

$$A_{\text{high}}(L,\beta J)\sim {\sinh }^{-1}\left(2\beta J+\ln \tanh \beta J\right),$$

$$|\mathcal{F}(p){|}^2$$

$$A_1(p,H,\beta J)=\frac{1}{2H^2}{\csc }^2\left(\frac{\pi p}{H}\right)$$

$$A(p,H,\beta J)=\sum _{d=0}^{\infty }{A}_d(p,H,\beta J)P(d,H,\beta J),$$

$$P(d,H,\beta J)=I_0^{-1}\left(2H{e}^{-2\beta J}\right){\left(\frac{H^d}{d!}\right)}^2{e}^{-4d\beta J},$$

$$A_{\text{low}}(H,\beta J)=\frac{e^{-2\beta J}}{2H}\frac{I_1(2H{e}^{-2\beta J})}{I_0(2H{e}^{-2\beta J})}.$$

$$A(H,\beta J)=A_{\text{high}}(H,\beta J)\frac{I_1\left(4H^2{A}_{\text{high}}(H,\beta J)\right)}{I_0\left(4H^2{A}_{\text{high}}(H,\beta J)\right)},$$

$$A_{\text{high}}(H,\beta J)=\frac{1}{4H}{\sinh }^{-1}\left(2\beta J+\ln \tanh \beta J\right).$$

$$A_{\mu \to \nu }(\lambda )=\min \left(1,e^{-\beta \left(E_{\nu }-E_{\mu }\right)}\frac{g_{\nu }(\lambda )}{g_{\mu }(\lambda )}\right).$$

$$\Pi (\mu \to \nu )=\sum _{\lambda \in {\Theta }_{\mu \to \nu }}{g}_{\mu }(\lambda )A_{\mu \to \nu }(\lambda ),$$

$$\Pi (\nu \to \mu )=\sum _{\lambda \in {\Theta }_{\nu \to \mu }}{g}_{\nu }(\lambda )A_{\nu \to \mu }(\lambda ).$$

$$P_{\text{acc}}=\{\begin{array}{ll}1-e^{-\beta (E_{ij}-E_0)},& \text{mbox}if{E}_{ij}\le E_0\\ 0,& \text{mbox}otherwise,\end{array}$$

$$A(p,H,\beta J,\Delta )=A(H,\beta J,\Delta )\times$$

$${\left[\{c(H,\beta J,\Delta )x^{-7/3}{\}}^b+{\csc }^{2b}(\pi p/H)\right]}^{1/b}.$$

Full text

I Introduction

The behaviour of field and current-driven magnetic domain walls is of great interest due to the promising applications that domain walls have in magnetic memory devices and logic devices Parkin et al. (2008); Allwood et al. (2005). When the magnetic material contains impurities, that can either be created by accident during the manufacturing process or injected intentionally, the motion of the domain walls at low enough driving currents is dominated by the movement of segments of the wall between so-called pinning sites Mougin et al. (2007). A critical field or current separates the creep regime as described above and the flow regime in which the driving force is large enough for linear behaviour Lemerle et al. (1998). Interestingly, this behaviour is found in a plethora of physical systems like propagating crack fronts Bouchaud et al. (1993) and moving vortices in superconductorsBlatter et al. (1994).

In this paper we choose to study the square-lattice Ising model where disorder is introduced in the form of random bonds between nearest neighbouring sites. Although the two-dimensional Ising model provides a highly idealized playground it does contain most of the characteristics of thin films of magnetic material with perpendicular anisotropy (PMA). In both cases the individual spins are either (almost completely) up or down resulting in a very sharp domain wall such that the Ising model relates closely to thin-film PMA’s.

The dynamics in the creep regime and static properties of domain walls are related through the roughness (or wandering) exponent. In case of domain walls in the 2D random-bond Ising model (RBIM) this exponent is $\zeta_{\text{RB}}=2/3$ Huse and Henley (1985); Blatter et al. (1994) and is in agreement with experiments on PMA’s Lemerle et al. (1998); Alvarez et al. (2010); Rodríguez-Rodríguez et al. (2010); Cayssol et al. (2004). In case of random-field disorder this exponent is $\zeta_{\text{RF}}=1$ T. (1997); Dong et al. (2012) in contrast to $\zeta_{\text{T}}=1/2$ when thermal fluctuations dominate the behaviour. When both random-bond disorder and thermal fluctuations are considered a cross-over takes place at a typical length scale characterized by the Larkin length $L_{c}$ Feigel’man et al. (1989) above which the random-bond disorder induced roughness dominates.

Here we set out to investigate the static properties of square-lattice RBIM domain walls by means of Monte Carlo simulations. Due to disorder the energy landscape becomes rough and in most cases results in highly degenerate ground states. Many conventional algorithms result in over-sampling parts of phase space due to the exponentially suppressed probability of leaving a certain meta-stable state. In the last few decades many new update algorithms have emerged Niedermayer (1988); Lee (1993); Houdayer (2001); Zhu et al. (2015) as well as other simulation schemes Bovier et al. (2015); Torlai and Melko (2016); Swendsen and Wang (1986) to overcome these obstacles. Although most loop algorithms follow local update schemes the big advantage is that they are not inhibited by energetic barriers between meta-stable states. We propose a new type of loop algorithm, a weighed-loop algorithm, where the loop is not constructed by a random walk but instead allows for different possible routes to have different weights assigned as to improve the acceptance probability of suggested updates.

The outline of this paper is as follows. First we discuss the model in more detail in section II and combine previously obtained theoretical results for the static properties of domain walls together with our low-disorder low-temperature expansion. In section III we present the weighed-loop algorithm, prove detailed balance and ergodicity, and discuss some other properties. We also compare the weighed-loop algorithm to Niedermayer’s algorithm Niedermayer (1988) and show that there exists a region in the parameter space consisting at low temperature and high disorder where the weighed-loop algorithm is decorrelates faster. In section IV we present the results and the analysis of our simulations on static properties of the domain walls after which we conclude with a discussion and conclusion in section V.

II Model and Theory

In the standard Ising model each lattice point $i$ is occupied by a spin up $\sigma_{i}=+1$ or a spin down $\sigma_{i}=-1$ particle. The Hamiltonian is defined by

[TABLE]

where $\Lambda$ is the entire lattice, $\mathcal{N}_{i}$ are the neighbouring spins of spin $i$ and $J_{ij}$ the bond strength between spins $i$ and $j$. Using a $L\times H$ size rectangular lattice and using boundary conditions with a periodic boundary along the $L$-direction and an anti-periodic boundary along the $H$-direction, the lattice is topological homeomorph to a Klein bottle. Therefore the lattice induces at least one domain wall. Throughout the paper we assume ferromagnetic bonds of varying strengths; $J_{ij}=J\pm J\Delta$ with $0\leq\Delta\leq 1$. The disorder is local and uncorrelated between bonds such that the disorder averages are given by $\langle J_{ij}\rangle_{\text{dis}}=J$ and $\langle J_{ij}J_{i^{\prime}j^{\prime}}\rangle_{\text{dis}}=J^{2}+J^{2}\Delta^{2}\delta_{ii^{\prime}}\delta_{jj^{\prime}}$.

The horizontal displacement of a domain wall at height $h$ for a realized configuration is given by the function $\operatorname{dw}(h)=\sum_{l=1}^{L}\sigma_{lh}$, where $0\leq h<H$ and $(l,h)$ is the Cartesian coordinate for the position of spin $\sigma_{lh}$. Below the critical temperature the probability that the wall is as wide as it is high is negligible. To ensure that the horizontal displacement is calculated correctly the anti-periodic boundary can be shifted such that it does not cross the domain wall. This elementary definition allows for a quick computation of the position of domain walls at the cost of a loss of accuracy due to isolated pockets of reversed spins. By definition overhangs are also not taken into account although these are ignored in derivations for roughness exponents as well. Both isolated pockets of spins and overhangs are energetically expensive excitations and so disappear for low temperatures and disorder.

The system is in thermal contact with a bath at temperature $T$. We define the inverse temperature $\beta:=1/(k_{B}T)$ with $k_{B}$ the Boltzmann constant. We start each simulation with exactly one domain wall. Since extra pairs of domain walls cost macroscopic amounts of energy we are practically ensured to always have exactly one domain wall below the critical temperature $T_{c}$. An example of a typical spin configuration and the resulting horizontal displacement of the domain wall for a system of size $100\times 100$ at $\beta J=0.7$ in the absence of disorder is shown in Fig. 1.

To study the static properties of domain walls in more detail we consider the discrete Fourier transform of the horizontal displacement

[TABLE]

and the structure factor $A(p,H,\beta J,\Delta)=\left\langle|\mathcal{F}(p)|^{2}\right\rangle$ for modes $0<p<H-1$ as a thermal and disorder average over all possible spin configurations. We will show that without disorder $A(p,H,\beta J,\Delta=0)=A(H,\beta J)\csc^{2}{\left(\pi p/H\right)}$ where $A(H,\beta J)$ is a scaling factor. For a system without disorder at high temperatures the scaling factor $A_{\text{high}}(H,\beta J)$ is given by

[TABLE]

which follows from capillary-wave theory Fisher et al. (1982) and results for the surface tension of the 2D Ising model Onsager (1944).

The next step is to make a low-temperature expansion for large domain walls in a system without disorder which will yield a low-temperature scaling factor $A_{\text{low}}(H,\beta J)$. We then combine this with the high-temperature scaling factor $A_{\text{high}}(H,\beta J)$ into $A(H,\beta J)$ that works for all temperatures.

Low-temperature expansion. Consider the ground state of a system without disorder. The domain wall is then given by a straight wall without any perturbations which do not contribute to the structure factor. Let $A_{d}(p,H,\beta J)$ be the average value of $A(p,H,\beta J)$ of all systems with exactly $d$ defects. A defect is defined as a single horizontal displacement over a certain length $b$. Without affecting the non-zero modes we can assume that all defects occur at only one side of the domain wall without a loss of generality. Similarly, since we are eventually interested in the absolute norm of the Fourier components we assume that defects start at position $h=0$. For a domain wall with a single defect the horizontal displacement function is given by $\operatorname{dw}(h)=\theta(h-b)$ with $\theta(x)$ the step function. Using Eq. (2) this results in

[TABLE]

Averaging over all non-trivial values of $b$ and taking the large $H$ limit while $0<p<(L-1)/2$ yields

[TABLE]

At low temperatures and for large systems we can safely assume that the defects are non-interacting such that $A_{d}(p,H,\beta J)=dA_{1}(p,H,\beta J)$. We can then sum over all possible number of defects to get the correct expression for the structure factor

[TABLE]

where $P(d,H,\beta J)$ is the probability of finding $d$ defects.

Each defect increases the energy of the system by $4J$ and so the probability of finding $d$ defects can be rewritten as $P(d,H,\beta J)=Z^{-1}(H,\beta J)g(d,H)\exp\left(-4d\beta J\right)$ with $Z$ the normalizing constant (or partition function) and $g(d,H)$ the number of possible ways $d$ defects can occur in a system of height $H$. We approximate $g(d,H)$ by assuming defects are still uncorrelated and so in the large $H$ limit $g(d,H)\approx\binom{H}{d}^{2}\approx\left(H^{d}/d!\right)^{2}$. Note that this is an overestimation as destructive interference between defects is not taken into account. The probability of finding $d$ defects is then given by

[TABLE]

where $I_{n}(x)$ denotes the Bessel function of the first kind. Using Eqs. (5)-(7) this results in low temperature scaling factor

[TABLE]

Scaling factor. To provide a single expression for the scaling factor we now naively combine the results at high temperatures Fisher et al. (1982) with the low-temperature expansion Eq. (8) such that the scaling factor $A(H,\beta J)$ behaves similarly to the approximations in their respective temperature regimes. First let $a(H,\beta J)$ be the unknown pre-factor of $A_{\text{high}}(H,\beta J)$ in Eq. (3). For low temperatures $\beta_{c}\ll\beta$ we have $A_{\text{high}}(H,\beta J)\approx 2a(H,\beta J)e^{-2\beta J}$. Using the approximation $I_{n}(z)\approx\left(z/2\right)^{n}/\Gamma(n+1)$ for small $z$ results in the approximation $A_{\text{low}}(H,\beta J)\approx e^{-4\beta J}/2$. Similarly just below the critical point $0<\beta J-\beta_{c}J\ll 1$ we have $A_{\text{high}}(H,\beta J)\approx a(H,\beta J)/4\left(\beta J-\beta_{c}J\right)^{-1}$. For large $z$ we have $I_{n}(z)\approx e^{z}/\sqrt{2\pi z}$ such that $A_{\text{low}}(H,\beta J)\approx e^{-2\beta J}/(2H)$. When $a(H,\beta J)=1/(4H)$ the low-temperature expansion of $A_{\text{high}}$ and high-temperature expansion of $A_{\text{low}}$ are connected. An elegant solution that fits both low- and high-temperature expansions of the scaling factor is then given by

[TABLE]

with

[TABLE]

III Weighed-Loop Algorithm

The Ising model has two-fold degenerate ground states that become highly degenerate when disorder is introduced in large enough systems. The energy landscape becomes rough and at low temperatures many algorithms get stuck in metastable states surrounded by energy barriers that become exponentially higher as temperature is lowered. Many single-spin-flip and cluster-growing algorithms suffer from this since they are dependent on all of the bonds around and inside the cluster. Although loop algorithms are local update algorithms the main advantage is that they only depend on the bonds that the loop follows and not on the bonds inside the surrounded cluster. A part of a configuration is shown in Fig. 2 that illustrates this point. For this example let the bonds be ferromagnetic (straight bond) or antiferromagnetic (wiggled bond). Each spin has three or more energetically favourable bonds (green) and so flipping any individual spin always increases energy. The dashed line shows a cluster of spins that, if flipped all together, will lower the energy. The main idea behind the weighed-loop algorithm is that bonds are iteratively added to loop to form a closed chain of bonds with a bias toward bonds that are energetically favourable to swap. Here we will explain the algorithm in more detail, prove detailed balance and ergodicity, and compare the autocorrelation times for domain-walls with those from Niedermayer’s algorithm.

Details algorithm. Let $G$ be the graph given by the lattice $\Lambda$ such that the spins are the faces of $G$. The algorithm first selects a starting node of $G$ at random. Given that the algorithm cannot go backwards each step after the first the algorithm has three possible edges of $G$ to add to the current chain. Each edge corresponds to a bond $J_{ij}$ and the algorithm adds the bond with a relative weight $W(J_{ij})=\exp(-\beta J_{ij}\sigma_{i}\sigma_{j})$. Although any weight can be assigned to any of the bonds we choose these particular values as the relative probability between two bonds that can be followed by the weighed-loop algorithm equal to the relative local Boltzmann weights of the bonds. Available bonds with equal energy have equal probability of being chosen and at lower temperatures bonds that will lower the energy more than other bonds have a higher probability. The algorithm continues adding edges of $G$ until a closed cycle in $G$ is constructed. The spins enclosed by the closed cycle define the cluster to flip. Since the lattice is homeomorph to a Klein bottle it is possible for the closed cycle to enclose the entire lattice. If such a loop is constructed, an additional loop is constructed. If both loops enclose the entire lattice, the spins enclosed by both loops define the cluster to flip. The algorithm will almost always construct a closed loop with an additional tail. Since the weights remain unchanged in the tail after the update it is used as part of the new loop.

Given that flipping the cluster of spins defined by the loop $\lambda$ transitions the system from $\mu$ to $\nu$, we define the acceptance probability of $\lambda$ as

[TABLE]

In the above equation $g_{\mu}(\lambda)$ is the probability of selecting $\lambda$ in $\mu$ which is defined by the product of the individual stochastic choices of the chaining of the bonds.

Properties. To prove that with the acceptance probability (11) the weighed-loop algorithm satisfies the detailed balance equation let $\Theta_{\mu\rightarrow\nu}$ be the set of all loops which transition $\mu$ to $\nu$ and let $\Theta_{\nu\rightarrow\mu}$ be the set of all loops which transition $\nu$ to $\mu$. The transitions probability are

[TABLE]

Each loop in $\Theta_{\mu\rightarrow\nu}$ is also a loop in $\Theta_{\nu\rightarrow\mu}$ hence summing over $\Theta_{\mu\rightarrow\nu}$ yields the same result as summing over $\Theta_{\nu\rightarrow\mu}$. Using Eq. (11) and the equilibrium distribution given by the Boltzmann weights $\pi(\mu):=\exp(-\beta E_{\mu})$ one can show that detailed balance $\pi(\mu)\Pi(\mu\rightarrow\nu)=\pi(\nu)\Pi(\nu\rightarrow\mu)$ holds.

To prove ergodicity of the weighed-loop algorithm it suffices to show it can perform any single spin flip at finite temperature. The starting node of $G$ is chosen at random from a flat distribution for all possibilities. The weight $W(J_{ij})$ for any edge is non-zero and so there is always a finite probability that a loop around a single spin is constructed. The resulting acceptance probability in Eq. (11) is, by construction, non-zero as well and so there is a non-zero probability the proposed loop is accepted and the single spin inside it flipped. Thus the weighed-loop algorithm is ergodic.

Although we made the choice that the relative weights of possible edges only depends on the energy difference one would get after flipping the bonds many other extensions are possible. The weights can be any type of function and it might be interesting to see the effects of a simulations with relative or absolute directional preferences. By doubling the relative weight of the edge in the forwards direction for example, the area enveloped by suggested loops would on average be larger.

Comparison with Niedermayer’s algorithm. We compare the weighed-loop algorithm to Niedermayer’s algorithm Niedermayer (1988). Niedermayer’s algorithm can be viewed as an extension of the Swendsen-Wang Swendsen and Wang (1987) and Wolff Wolff (1989) algorithms that can be applied to glassy spin systems. Niedermayer’s algorithms is almost the same as Wolff’s algorithm except that any spin can be added to the cluster with varying probability. For Niedermayer’s algorithm, with $E_{ij}:=\sigma_{i}\sigma_{J}J_{ij}$, the acceptance probability is given by

[TABLE]

with $E_{0}$ a free to choose parameter. If $E_{0}$ is larger than all possible $E_{ij}$ the acceptance probability for flipping the cluster is always $1$. Unfortunately, there is not much known about the effect of the $E_{0}$ on the correlation times Newman and Barkema (1999). We have set $E_{0}=J+\Delta J$ so that the probability for flipping a selected cluster is always equal to unity.

For both Niedermayer’s algorithm and the weighed-loop algorithm we ran multiple simulations, each exactly $10$ minutes on the same system, on the same $24$ different disorder configurations for varying $H$, $\beta$, and $\Delta$. We measured the autocorrelation time of the squared amplitude of the first mode of the domain wall $C(t):=\left\langle\left(|\mathcal{F}(1,t)|^{2}-\left\langle|\mathcal{F}(1)|^{2}\right\rangle\right)\left(|\mathcal{F}(1,0)|^{2}-\left\langle|\mathcal{F}(1)|^{2}\right\rangle\right)\right\rangle$ for both algorithms. Here $\mathcal{F}(1,t)$ is the first Fourier mode, the mode with the largest correlation time, at time $t$ and brackets denote the disorder and thermal average. For small $t$ we can assume exponential decay $C(t)\sim\exp(-at)$ with equal pre factors for both algorithms. Figure 3 shows typical autocorrelation times for both algorithms evaluated on the same disorder configuration. We measure $a$ for both algorithms for varying parameters $H$, $\beta$, and $\Delta$ where larger values of $a>0$ lead to faster decorrelation. The results are listed in Table 1. We see a clear trend that the weighed-loop algorithm has faster decreasing correlations in the parameter regime where conventional algorithms tend to get stuck (low temperatures and high disorder). Additionally, on larger time scales as can be seen in Fig. 3, we observe that the weighed-loop algorithm also decorrelates exponentially on longer time scales. This might be an indication that the weighed-loop algorithm tends to hop more frequently between far-removed local minima.

IV Results

In this section we present several results obtained using the weighed-loop algorithm. First we compare simulations in the absence of disorder with the theoretical result Eqs. (9,10) obtained from our low-temperature expansion in combination with the high-temperature expression resulting from the surface tension of the 2D Ising model Onsager (1944) and capillary-wave theory Fisher et al. (1982). Next we induce disorder in the system and measure the typical length $L_{c}$ at which the system switches from a roughness dominated by thermal fluctuations to higher length scales at which the roughness is determined by disorder.

In the absence of disorder. First we look at the simplest case $\Delta=0$ for which only thermal averaging is needed. The scaling factor $A(H,\beta J)$ for both $H=64$ and $H=256$ is shown as a function of temperature in Fig. 4(a). The data is in agreement with theory over the complete temperature-range that was used in the simulations. Error margins close to the critical temperature, at low values for $\beta J$, are much higher due to larger thermal fluctuations. This is also the region in which overhangs and local spin pockets start to occur. All these phenomena result in a less accurate approximation of structure of the domain wall. For extremely low temperatures the domain walls are mostly perfectly straight and so the error margins are increased due to the low probability for the wall to be in an excited state. Figure 4(b) shows a data collapse of the scaled mode-dependent factor $A(p,H,\beta J)/A(H,\beta J)$ as a function of $p/H$ to show the universal behaviour in the absence of disorder. The collapse indeed works best for moderate temperatures.

Introducing disorder. The addition of randomness to the system brings with a couple of complications. Besides the usual thermal averaging one also needs to perform disorder averaging resulting in much longer simulation runs than is the case without disorder. Additionally, the energy landscape becomes rough and so many more pockets of local energetically favourable configurations far away removed in phase space need to be explored to measure averages accurately. A trade-off is to be made between simulating large systems with faster dynamics and easier averaging at higher temperatures and less disorder or simulating smaller system in the more problematic regime of the $\Delta$-$\beta$ parameter space.

When disorder is introduced in an infinite system there is a typical length scale, known as the Larkin length $L_{c}$ Feigel’man et al. (1989), above which the roughness of the wall is determined by the disorder instead of thermal fluctuations. The roughness exponent from thermal fluctuations $\zeta_{\text{T}}=1/2$ differs from the exponent for random-bond disorder, for one-dimensional walls, $\zeta_{\text{RB}}=2/3$ Huse and Henley (1985); Blatter et al. (1994). The roughness exponent also appears in the average of the fourier modes. One finds $A(p,H,\beta J,\Delta)\sim p^{-(1+2\zeta)}$ in the low-mode limit, which is in agreement with $A(p,H,\beta J,\Delta)\sim\csc^{2}(\pi p/H)$ and $\zeta_{\text{T}}=1/2$ in the absence of disorder. The Larkin length can be translated into a typical mode $p_{c}/H:=1/(2L_{c})$ at which $\zeta$ changes from thermal to disorder-induced. Since the transition is relatively sharp we can make a fit for all $p$ of the form

[TABLE]

Here $A(H,\beta J,\Delta)$, $b$, and $c(H,\beta J,\Delta)$ are free fit-parameters from which $p_{c}$ can be extracted. For large systems and high enough temperatures we can scale out the height dependence $A(H,\beta J,\Delta)\sim 1/H$ as can be seen from Fig. 4(a).

In Fig. 5 we look at $A(p,H,\beta J,\Delta)$ as a function of $p$ while varying the height $H$ of the system and keeping $\beta J=2.0$ and $\Delta=0.15$ fixed. Simulations were performed for $H=8,16,24,\dots,128$. For very small systems $A(H,\beta J,\Delta)\sim 1/H$ does no longer hold and additionally simulations of such systems does not yield data above the Larkin length $L_{c}$. Therefore, the unweighed fit of the form in Eq. (IV) is performed on data from systems with $H\geq 48$ for which the higher-mode portion already sufficiently collapses. The figure also shows the best unweighed fit as a solid black line as well as a dashed and dotted line to visualize the power-law behaviour in the different regimes. We find a transition at $p_{c}\approx 0.04H$ or equivalently $L_{c}\approx 12.5$ in units of lattice spacing. Clearly one can extract the transition from a series of simulations of systems of slightly varying heights larger than the Larkin length rather well even for systems with disorder. We also find that $A(H,\beta J,\Delta)\sim 1/H$ still holds, in the large-system limit at low enough temperatures, even when disorder is introduced to the system.

Next we fix the height of the system $H=64$ as well as the temperature $\beta=2.0$ while varying the amount of disorder $\Delta=0.00,0.02,0.04,\dots,0.30$. The increase of disorder leads to a decrease in Larkin length and consequently an increase in where the transition $p_{c}/H$ occurs. Additionally, the overal pre factor $A(H,\beta J,\Delta)$ increases and thus an intuitive collapse for all $p/H$ is not possible. Instead Fig. 6(a) displays $A(H,\beta J,\Delta)$ as a function of $p/H$ for all values of $\Delta$ at which simulations were performed. The solid black line goes through $p_{c}/H$ at each value for disorder and thus indicates the crossover point. The pre factor $A(H,\beta J,\Delta)$ is shown as a function of disorder in Fig. 6(b) together with $A(H,\beta J(1-\Delta))$ as defined in Eq. 9 which suggests an exponential increase in roughness as disorder is increased.

We also performed measurements at varying temperatures $\beta=0.5,0.6,\dots,2.8$ at height $H=64$ and $\Delta=0.15$. Figure 7(a) shows $A(p,H,\beta J,\Delta)$ together with a best fit to Eq. IV as a function of $p/H$, together with the transition points denoted by the solid black line. As in the absence of disorder we observe an increase in roughness as temperature is increased. We compare the pre factors $A(H,\beta J,\Delta)$ Fig. 7 similar to Fig. 4(a) but now for systems at equal height at $\Delta=0.0$ and $\Delta=0.15$. At high temperatures the thermal fluctuations overshadow the effects of disorder and so no difference is observed for $A(H,\beta J,\Delta)$ and Eqs. (9,10) still hold. When temperature is lowered, however, disorder plays a more prominent role and $A(H,\beta J,\Delta)$ seems to plateau to the case without disorder.

We find good agreement between the data from our simulations using the weighed-loop algorithm and Eqs. (9,10 for large systems without disorder. The introduction of disorder to the system leads to different behaviour at different length scales. At lower length scales, in the thermal regime, $A(p,H,\beta J,\Delta)\sim\csc^{2}(\pi p/H)\sim p^{-2}$ still holds even for small systems or when there is disorder in the system. As expected we find that the roughness is increased when temperature is increased, larger systems are studied, or disorder is increased. When there is disorder, however, the roughness behaves differently at large length scales $A(p,H,\beta J,\Delta)\sim p^{-7/3}$ which we could also verify with our simulations. The typical length scale $L_{c}$, at which the roughness depends equally on disorder and thermal fluctuations, decreases when disorder is increased or temperature is decreased but is independent on the size of the system.

V Discussion and Conclusion

In this paper we have presented the weighed-loop algorithm that can be implemented for any graph-type model for which the energy only depends on the (one-dimensional) edge of a cluster. Here, we applied it to the square-lattice random bond Ising model (RBIM). More specifically, we have investigated the typical roughness of one-dimensional domain walls in the RBIM. The basic principle of the weighed-loop algorithm is that a weighed directed random walk is performed on a graph. Instead of weighing each possible direction equally we choose to favour directions that are energetically favourable when chosen. This increases the acceptance of proposed updates and thus decreases correlation times especially at lower temperatures or when there is disorder in the system.

The flexibility of the weighed-loop algorithm lies in the function that determines the weights. Throughout this paper we have chosen a function that only depends on the energy contribution of the optional bonds. However, one could also induce a relative or absolute preferred direction by adjusting the weights of optional bonds in different directions. Increasing the odds on straight paths would increase the area of the enclosed cluster for example. It would be interesting to investigate the effects of different weight-functions on correlation times. Additionally, the algorithm can be used in combination with other update schemes and parallel tempering to reduce thermalization.

Next we compared the weighed-loop algorithm to one of many popular algorithms, Niedermayer’s algorithm Niedermayer (1988), used to study spin glassy systems. We compared the autocorrelation functions of the square amplitude of the first Fourier mode for domain walls in the RBIM for the two algorithms. The weighed-loop algorithm outperforms Niedermayer’s algorithm at low temperature high disorder; the part of parameter space at which many issues arise due to the rough energy landscape and the exponentially high energy barriers.

To provide a good description for all temperatures below the critical temperature for the roughness of domain walls in the absence of disorder we conciliated the high-temperature expression with our low-temperature expansion. We show good agreement with data obtained with the weighed-loop algorithm and the description in Eqs. (9,10).

Once disorder is introduced to the system there appears to be a minimal disorder-induced roughness even when temperature tends to zero. We have investigated the effects of the amount of disorder, temperature, and system size independently. We find that the overall roughness of domain-walls increases as temperature, disorder, or system size is increased. Besides the change in overall roughness we also observe different behaviour of the roughness at different length scales conform the known theory. For low modes, or equivalently lengths larger than the Larkin length $L_{c}$ Feigel’man et al. (1989), the roughness exponent $\zeta_{\text{RB}}=2/3$ is different from the exponent at higher modes $\zeta_{\text{T}}=1/2$. The cross-over point $p_{c}/H$, at which the transition from disorder-induced to thermal-induced roughness takes place, depend on the temperature and amount of disorder in the system. It increases as disorder is increased or temperature is lowered. In the future we would like to explore this more quantitatively to accurately describe how $L_{c}$ depends on these parameters, both with simulations and theoretically.

VI Acknowledgements

We thank E. van der Bijl and G. T. Barkema for valuable feedback, useful discussions, and insights at the early stages of this work. RK is grateful to the Institute for Theoretical Physics at Utrecht University and the Lorentz Institute for Theoretical Physics at Leiden University for the hospitality during the course of this work. This work is part of the D-ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).