# Numerical investigation of the thermodynamic limit for ground states in models with quenched disorder

## Abstract

Numerical ground state calculations are used to study four models with quenched disorder in finite samples with free boundary conditions. Extrapolation to the infinite volume limit indicates that the configurations in “windows” of fixed size converge to a unique configuration, up to global symmetries. The scaling of this convergence is consistent with calculations based on the fractal dimension of domain walls. These results provide strong evidence for the “two-state” picture of the low temperature behavior of these models. Convergence in three-dimensional systems can require relatively large windows.

The structure of the thermodynamic set of states of a system in statistical mechanics is studied formally through the infinite volume limits of correlation functions . If a nested sequence of systems with given Hamiltonian and boundary conditions has spin correlation functions that converge in the infinite volume limit, a thermodynamic state can be defined. For example, in a ferromagnet with fixed, positive fields at the boundary, the single-spin correlation function converges to a positive value, defining an “up” state. For disordered spin systems, the question of the number of thermodynamic states is a subtle one . Whether there are many thermodynamic states in some sense  or a small number of states related by simple global symmetries  (e.g., two spin-flip related states in an Ising spin glass) has been a most controversial point for low-dimensional systems. Part of this debate has been over what are the most useful methods for determining the structure of thermodynamic states, spin overlaps $`P(q)`$  or correlation functions in subsystems  and it is unclear whether Monte Carlo simulations at finite temperature can be used to study large enough systems .

This letter describes the results of numerical computations which address the structure of states in disordered systems in the thermodynamic limit, at zero temperature. Two two-dimensional models, an Ising spin glass and a charge density wave (CDW) model (also referred to here as an elastic medium model) and two three-dimensional models, a CDW model and a dimer matching model that is equivalent to non-intersecting lines in a random medium (similar to vortex lines in type-II superconductors), were studied. The ground states were computed for a sequence of free boundary conditions and the configurations in a fixed finite subsystem (or “window”) were compared. This study is a particular instance of the numerical approach suggested by Newman and Stein , who have presented detailed arguments that the existence of many states, as in the Parisi solution  of the mean field spin glass, gives rise to “chaotic size dependence” . The principle result derived from the simulations presented here is that the window configurations converge to a single fixed configuration with probability one. These computations strongly support the picture of a small number of ground states related by global symmetries, consistent with the droplet model . The details of the convergence to a single fixed configuration as the boundary grows has a scaling behavior which is well-described by a simple picture of domain walls.

The 2d spin glass model (SG) studied has spins $`s_i=\pm 1`$ defined at lattice points $`i`$, with Edwards-Anderson Hamiltonian  $`H_{\mathrm{SG}}=_{ij}J_{ij}s_is_j`$, where $`J_{ij}`$ is chosen independently from a Gaussian distribution for all nearest neighbor bonds $`ij`$. This model is believed to be paramagnetic at finite temperature, but is a spin glass at $`T=0`$; minimal energy large scale excitations of size $`L`$ have an energy $`E(L)L^{\theta _{\mathrm{SG}}}`$ with $`\theta _{\mathrm{SG}}0.27`$ . The discretized CDW or elastic medium model in two dimensions (E2) studied here is also equivalent to a disordered substrate model or vortex lines in two dimensions pinned by quenched disorder . The configurations in this model are defined by complete dimer coverings of a hexagonal lattice, with the Hamiltonian being the sum over covered dimers $`d`$ of dimer weights $`w_d`$, $`H_{E2}=_dw_d`$, where the $`w_d`$ are chosen for each bond from a uniform distribution. In mean field replica calculations, matching problems are found to have replica symmetric solutions . A mapping of the dimer model to a discrete height representation $`h`$ can be made ; the variable $`h`$ corresponds to the scalar phase displacements in CDW models. This model is believed  to have a finite temperature phase transition, with the height-height correlations $`h(r)h(0)`$ behaving as $`\mathrm{ln}(r)`$ in the high-$`T`$ phase and as $`\mathrm{ln}^2(r)`$ in the low-$`T`$ phase. In this model, $`\theta _{\mathrm{E2}}=0`$ ($`E(L)\mathrm{const}`$.) The model E2 can be extended to three dimensions in two distinct ways; both are both studied here. One extension is that of dimer covering (matching) on a cubic lattice (M3), which can be mapped to a set of vortex lines with hard-core repulsion . It has a Hamiltonian identical to that for E2, with the covering dimers a subset of the edges in a simple cubic lattice. The other 3-D model is the three dimensional CDW or elastic medium model (E3) ; in the continuum limit, $`\theta _{E3}=1`$ (consistent with numerics in Ref. .) The low temperature phase of the elastic medium models have been studied using both  renormalization group and replica symmetry breaking techniques, which are usually, though not exclusively, interpreted physically as describing systems with few states or many states, respectively.

These models of disordered systems were studied using polynomial-time combinatorial optimization algorithms . The spin glass was studied on a triangular lattice, using the method developed by Barahona , rather than the string method which is often used . The minimum weighted matching algorithm  used for the implementation of Barahona’s algorithm was the algorithm described in Ref. . Calculations were made for at least $`10^3`$ samples of up to $`512^2`$ spins. The model E2 can be mapped to a bipartite matching problem  and was solved using the algorithm of Ref.  for at least $`10^3`$ samples of sizes up to $`1024^2`$ sites. The same algorithm was used for model M3, 3-D matching, with up to $`128^3`$ sites with at least $`10^3`$ samples, while the push-relabel maximum flow algorithm as implemented in Ref.  was used to study the 3-D elastic medium model E3 (up to $`64^3`$ sites with at least $`10^3`$ samples). The algorithms used determine ground states up to global symmetry transformations. For example, in the spin glass, unsatisfied bonds (bonds with $`J_{ij}s_is_j<0`$) are calculated, rather than $`s_i`$. Configurations related by symmetries are considered identical here, so that a “two-state” picture for spin glasses naturally appears as a single state in the computations .

The effect of system size was studied extensively for free boundary conditions. The disorder realizations for each sample $`S_L^\alpha `$ of linear size $`L`$ (with $`L^d`$ spins or sites) were generated so that $`S_L^\alpha `$ was a subsystem of a given infinite volume sample $`\alpha `$. Two finite samples $`S_L^\alpha `$ and $`S_L^{}^\alpha `$ had the same quenched disorder in their intersection. Each finite sample was centered at an origin $`C`$, so that a sequence of samples $`S_L^\alpha ,S_L^{}^\alpha ,S_{L^{\prime \prime }}^\alpha ,\mathrm{}`$ with $`L<L^{}<L^{\prime \prime }<\mathrm{}`$ gives a nested set of square or cubic samples centered about $`P`$. The $`L\mathrm{}`$ limit could then be numerically studied for a number of infinite samples $`\alpha `$. The free boundary conditions were assumed to be typical for the models SG and M3. In the elastic models, free boundary conditions give ground states with lower energy than boundary conditions that would introduce a uniform strain in the elastic models in the infinite volume limit; such uniform strain states are not considered here.

The configuration differences for samples of different sizes $`L<L^{}`$ were computed by comparing the exact ground states in the volumes of size $`L^d`$ where $`S_L^\alpha `$ and $`S_L^{}^\alpha `$ overlapped. Spin glass ground states in two dimensions were compared by finding the differences in unsatisfied bonds. An example of such a ground state comparison by bond overlap is shown in Fig. 1. For the models with dimer matchings (E2 and M3), the configurations are compared by finding the symmetric difference of the dimer sets in the common volume. The natural comparison for the height configurations for the model E3 is to determine where the gradients of the heights in the intersection volume differ.

The primary quantity of interest that was computed was the (sampled) probability that a change in boundary conditions resulted in any change in the ground state configuration in a window of size $`w`$ centered at $`C`$. The probability $`P(L^{},L,w)`$ is defined as the probability that the configuration in the window region changes as the system size is increased from $`L`$ to $`L^{}`$, that is, that the ground state configuration for $`S_L^{}^\alpha `$ differs from that for $`S_L^\alpha `$ in the volume of size $`w`$ centered at $`C`$. This quantity was estimated by sampling over a large number of samples $`\alpha `$ for various $`L^{}`$, $`L`$, and $`w`$ . This measurement is sensitive to all gauge invariant spin correlation functions in the window volume.

A plot of the data for $`P(L^{},L,w)`$, as a function of $`w`$ for various $`L^{}`$ and $`L`$, is shown in Fig. 2 for the spin glass. Assuming scale invariance, $`P`$ should be a function of the two ratios $`L^{}/L`$ and $`L/w`$. The data is consistent with this hypothesis, for large values of $`w`$ and $`L`$. For fixed $`L/w`$, $`P(L^{},L,w)`$ approaches a constant for large $`w`$ or large $`L`$. Note that to within error estimates, $`P(L^{},L,w)`$ is independent of $`L^{}`$ for $`L^{}/L=2,4,8`$: the probability of change in a finite window is approximately independent of the magnitude of expansion in the boundary, for $`L^{}2L`$ ($`P`$ does decrease noticeably as $`L^{}L`$.) In addition, for fixed $`w`$, $`P(L^{},L,w)`$ decreases approximately as a power law in $`L`$ (by the even vertical spacing of the data points for $`P1`$.) The data strongly suggest, by extrapolation to larger values of $`L`$ for fixed $`w`$, that the probability of changing the configuration in a window of size w goes to zero for $`L^{}/L\mathrm{}`$ as $`L\mathrm{}`$, implying convergence to a unique thermodynamic ground state (up to global symmetries).

The data can be explained by simple assumptions about the convergence of the configurations as $`L\mathrm{}`$ and the properties of domain walls or defect lines. For the spin glass model, induced domain walls are lines where the bonds change from satisfied to unsatisfied or vice versa. In models which are represented by a matching (E2 and M3) defects are also line objects and are composed of bonds where the dimer covering changes. In model E3, the induced walls are surfaces where the height gradient changes. Defect lines have fractal dimension $`d_f^{\mathrm{SG}}=1.27(1)`$ for model SG and $`d_f^{\mathrm{E2}}=1.25(1)`$ for model E2 . For the 3-D elastic medium, a shift in boundary conditions introduces a domain wall of dimension $`d_f^{\mathrm{E3}}=2.60(5)`$ , while localized string defects were computed during the course of this work to have fractal dimension of $`d_f^{\mathrm{M3}}=1.65(4)`$ in the 3-D matching model. If the fractal dimension of the defects is large enough ($`d_f>d/2`$) that no more than $`O(1)`$ defects of size $`L`$ can co-exist in the volume $`L^d`$, the expected number of defects of linear size $`L`$ introduced upon expansion to size $`L^{}`$ is bounded above by a constant. Whether boundary changes do induce a number of defects that saturate this bound is less clear a priori. For the models where $`\theta 0`$, finite changes at the boundary are likely to induce as many defects as possible, as the large scale defect cost is comparable to the cost of local changes at the boundary. The probability that a line or surface will intersect a window of size $`w`$ is then the ratio of the number of volumes of size $`w^d`$ that intersect the defect to the number of areas of size $`w^d`$ in the area $`L^d`$, giving the form

$$P(L^{},L,w)=c(L^{}/L)(L/w)^\kappa ,$$
(1)   
for large $`L/w`$, with $`\kappa =dd_f`$ by the supposition of a single dominant defect and, by these numerical results, the coefficient function $`c(L^{}/L)`$ quickly converges to a constant value for $`L^{}/L2`$. This form can be checked by plotting $`P`$ as a function of $`L/w`$ and comparing with a line of slope $`d_fd`$, as shown in Fig. 3. The match between this prediction and the data in $`d=2`$ is quite good; a two-parameter fit (varying $`c(\mathrm{})`$ and $`\kappa `$) gives exponents that agree with $`\kappa =dd_f`$ to within $`0.05`$ for models SG and E2. Differences of this order are within statistical fluctuations and apparent finite size effects. In addition, numerical study of a number of configurations for three values of $`L`$ (e.g., $`S_L^\alpha `$, $`S_L^{}^\alpha `$, $`S_{(L^{})^2/L}^\alpha `$) for the $`d=2`$ spin glass suggests that the location of $`L`$-scale defects in a volume is nearly independent of $`L^{}`$, giving more support to the conclusion that there is convergence to a unique state in these models.

The 3-D results also indicate convergence to a single state, as $`P(L/w)0`$ for $`L/w\mathrm{}`$. The quantitative fits are also consistent with a defect picture, but have a larger uncertainty. For the 3-D elastic medium, the data are consistent with $`P(L/w)^{d_f^{\mathrm{E3}}d}`$ for fixed $`w`$, as shown in Fig. 3(c), though larger sample sizes would be useful. For the problem $`M3`$, the behavior $`P(L/w)^{1.35}`$, with $`\kappa =dd_f`$, can be fit to the largest $`L/w`$ values, though only over a small range. Note that $`P>0.9`$ for $`L/w4`$ and $`P>0.5`$ for $`L/w=8`$. Under an expansion $`L^{}/L=2`$, the configuration in a window usually changes for $`L/w<8`$. Such change in small systems mimics the predictions of a many-states picture.

In summary, the infinite-volume limit for four model disordered systems was studied numerically by computing ground state configurations in fixed volumes embedded in systems of successively larger sizes. Strong evidence was found for convergence to a unique state (up to global symmetries), even in cases where $`\theta 0`$. The convergence to a unique state in $`d=2`$ can be understood in detail by estimating the chance of a defect wall intersecting a given area upon a boundary change. The 3d model results are more qualitative: while it appears that the system converges to a unique state, the ratio of scales ($`L^{}/L`$, $`L/w`$) required is larger, so that systems of size $`L>50`$ are needed. Polynomial ground state algorithms are not available for the 3d spin glass and this system is not directly addressed here, but these results suggest that one should be cautious in interpreting finite temperature Monte Carlo results  and ground state calculations  in small systems.

I would like to thank Daniel Fisher for stimulating discussions and comments on a draft of this paper. At the completion of this paper, I became aware of related results for the spin glass in $`d=2`$ by Palassini and Young . This work was supported in part by the National Science Foundation (DMR-9702242) and by the Alfred P. Sloan Foundation.