Distribution of Lee-Yang zeros and Griffiths singularities in the model of spin glasses
We investigate the distribution of zeros of the partition function of the two- and three-dimensional symmetric Ising spin glasses on the complex field plane. We use the method to analytically implement the idea of numerical transfer matrix which provides us with the exact expression of the partition function as a polynomial of fugacity. The results show that zeros are distributed in a wide region in the complex field plane. Nevertheless we observe that zeros on the imaginary axis play dominant roles in the critical behaviour since zeros on the imaginary axis are in closer proximity to the real axis. We estimate the density of zeros on the imaginary axis by an importance-sampling Monte Carlo algorithm, which enables us to sample very rare events. Our result suggests that the density has an essential singularity at the origin. This observation is consistent with the existence of Griffiths singularities in the present systems. This is the first evidence for Griffiths singularities in spin glass systems in equilibrium.
In order to shed a unique light on the problem of finite-dimensional spin glasses, we discuss the distribution of zeros of the partition function of the Ising spin glass in the complex field plane, the Lee-Yang zeros . All Lee-Yang zeros of a ferromagnetic system are proved to lie on the imaginary axis in the complex field plane from the circle theorem . However, locations of Lee-Yang zeros for spin glass models in general do not obey simple rules. It is therefore interesting and important to study the distribution numerically in order to extract useful information on phase transitions in finite-dimensional spin glasses from a point of view different from conventional numerical methods like Monte Carlo simulations. The first study along this line was due to Ozeki and Nishimori . From data for relatively small systems, they found the tendency which may be consistent with the existence of the Griffiths singularity . Our present work is to expand their study to much larger systems and carry out quantitative evaluations of the density of zeros.
In the diluted ferromagnet, the Griffiths singularity means that the free energy is a non-analytic function of the field at the origin in the temperature range , where is the critical temperature and is that of the corresponding non-random system. This temperature range is called the Griffiths phase. This unusual phenomenon in the otherwise-paramagnetic phase is caused by rare events of bond configurations, large-size connected clusters generated with small but non-vanishing probabilities. This fact was connected with an essential singularity in the density of zeros at the origin  because large clusters tend to have zeros close to the origin at temperatures below . The Griffiths singularity in the diluted ferromagnet has also been studied through the behaviour of the inverse susceptibility analytically  and numerically .
On the spin glass problem, Randeria et al. studied the dynamics and suggested the existence of the Griffiths phase also in the spin glass systems . There are few studies of Lee-Yang zeros of spin glass models in equilibrium except the work by Ozeki and Nishimori  although there are some papers [8, 9] on the distribution of zeros in the complex temperature plane, Fisher zeros . In the present paper we evaluate the density of Lee-Yang zeros of the spin glass systems by an importance-sampling Monte Carlo algorithm [6, 11, 12, 13, 14]. Our results give evidence that the density of zeros has an essential singularity, which suggests the existence of the Griffiths singularity in spin glasses.
The outline of this paper is as follows. In section 2 we explain the exact evaluation of the partition function by numerical transfer matrix. Then, we show the distribution of zeros of the Ising model on the complex field plane in section 3. In section 4, we analyze zeros that have direct relevance to the phase transition and compare our investigation with known characteristics of phase diagrams. In section 5, we examine the density of zeros using the importance-sampling Monte Carlo algorithm which enables us to sample very rare events, and we discuss the existence of the Griffiths singularity in the Ising model. We present our conclusion in section 6.
2 Exact partition function
We consider the Ising model in an external magnetic field on the square and the simple cubic lattices with cylindrical boundary conditions, free in one direction and periodic in the others, which are suited for the transfer matrix method because we take the trace of spins one by one in spiral order. The Hamiltonian is
where denotes nearest neighbour pairs and . The interaction is chosen as or with probability . Using the number of states for given and , the partition function is expressed in terms of a two-variable polynomial as
where (fugacity), is the number of sites (assumed to be even), and is the number of bonds on the lattice.
We focus on the -dependence of for a fixed , since we are interested only in the field term for the investigation of Lee-Yang zeros. Thus, we treat the partition function as a single-variable polynomial
To locate Lee-Yang zeros for the present model, we have to evaluate the partition function for a given set of random interactions analytically and solve the equation for .
The first one is to evaluate the partition function at different external fields by using the numerical transfer matrix method. Then, we solve the set of linear equations for . In this way we know the values of all coefficients of the polynomial and therefore we can solve . In this process, since both and are treated as numbers, we have to be careful to keep very high precisions to avoid rounding errors in solving . In addition, this method is not very efficient because it requires a repetition of essentially the same calculations. An advantage is that the required memory size is very small and is easily parallelized on the grid computer. We use this technique for systems larger than (see the Appendix).
The second method, called the canonical transfer matrix, is to numerically give analytic partition functions at a fixed temperature for all values of . This method evaluates (written as in ) for all at once. The required memory size is times larger than in the first method, but this process is done with only one calculation. Thus, the second method is faster than the previous method on a single CPU, but attention needs to be paid to the precision since we treat as a numerical value.
The third method is called the microcanonical transfer matrix , which evaluates perfectly exact partition functions for all values of and using symbolic manipulations on the computer. This algorithm does not involve numerical errors. However, the memory requirement is times larger than in the second method.
In consideration of the balance between the computational time and the required memory size, we mainly use the canonical transfer matrix method in the present work. This method accompanies rounding errors. The necessary precision depends on the ratio of the largest and the smallest coefficients, which increases very rapidly with the system size or with the temperature decrease. We have verified the reliability of this approach for the ferromagnetic Ising model of sizes comparable to the spin glass system we are interested in as explained in the Appendix.
The system size we calculated is up to in two dimensions and in three dimensions, where the last digit is for the free boundary direction and the others are for the periodic boundary conditions. We investigated at least about samples of random interactions for each system size. We solved the polynomial equation by mathematica by specifying the precision to appropriate values.
3 Distribution of zeros
The free energy for a system with quenched randomness is defined as the average over samples for all possible bond configurations. This implies that the partition function of a quenched system is regarded as being given by the product of partition functions of all possible configurations,
where denotes the quenched free energy and denotes the set of random interactions. This fact allows us to plot zeros of all randomly-chosen samples simultaneously on the complex plane of . The temperature will be measured in units of . A schematic diagram for the distribution of zeros is presented in figure 1.
Figure 2 shows the temperature dependence of the distribution of zeros for systems of size in two dimensions and in three dimensions. Throughout this paper the range of the real axis is from to in two dimensions and from to in three dimensions, and the imaginary axis ranges from to .
Generally, zeros of both dimensions approach the real axis with decreasing temperature. Nearest zeros to the origin lie on the imaginary axis for both dimensions and all temperatures. Let us call the location of such zeros the edge (see figure 1). A clear difference between two and three dimensions appears at low temperature. Zeros lying off the imaginary axis but close to the origin (to be called the body) form a wedged-like shape in two dimensions whereas the body looks relatively rounded in three dimensions.
We next discuss the dependence of the distribution on system size. Figures 3 and 4 show the system-size dependence at from system size to in two dimensions (figure 3) and from to in three dimensions (figure 4). The width of the distribution in the transverse direction does not change significantly as the system size increases and is likely to stay almost the same for sufficiently large systems. Both the edge and the body move toward the real axis as the system size increases in both dimensions. Especially, the body (i.e. zeros off the imaginary axis) tends to move more significantly than the edge. In the thermodynamic limit, if the body reaches the real axis away from the origin, it may imply that an ordered phase (that may exist along the real axis at and around the origin) would survive until the field exceeds a finite value corresponding to the point where the zero (the body) reaches the real axis.
To analyze the density more closely, we divide the area of , (two dimensions) or , (three dimensions) into boxes and count the number of zeros in each box. Figures 5 and 6 show the resulting density plots in two and three dimensions, respectively. It is readily seen that the density is very high on (and around) the imaginary axis in all cases. Regions of low density are cleared to white as can be verified by comparison of figure 3 (left-bottom panel) and figure 5 (bottom panel), for example. Zeros in those regions are expected to have no influence on the system properties in the thermodynamic limit. A marked difference between two and three dimensions is that the former develops a very sharp wedge in the body of the density profile as the size increases whereas, in the latter, the body remains rounded at its bottom part as commented above already. In two dimensions, the edge approaches the origin while the body does not show such an approach to the real axis at a faster rate than the edge with size increase. Also in three dimensions, the body seems to approach the real axis, but the difference between two cases in figure 6 is not clear since the difference of size is small.
4 Approach to the real axis
Based on the qualitative observation in the previous section, we analyze our data quantitatively in this section. We assume that the behaviour of the nearest zero to the real axis of each sample determines the phase transition of the model in the thermodynamic limit. We therefore calculate the average location of the nearest zero of each sample and analyze the system-size dependence of this average location. For two dimensions, we choose the temperatures as and (in the presumed Griffiths phase). In three dimensions, the investigated temperatures are (the spin glass transition temperature ) and (in the spin glass phase).
Figures 8 and 8 show the real part (right part of figures) and the imaginary part (left part of figures) of the average location of nearest zeros as functions of the inverse of the linear size in two and three dimensions 111See the Appendix for reasons to choose as the linear size in three dimensions..
We find that the average location of nearest zeros on the imaginary axis is lower than those off the imaginary axis for both dimensions as has been seen in previous figures. We may safely conclude that, if there is a phase transition caused by an approach of zeros to the origin, zeros on the imaginary axis reach the origin before those off the imaginary axis.
In two dimensions, the imaginary part of the average location seems to approach the origin linearly with the inverse of system size. Solid lines in figure 8 show the best fits to linear functions. The average location has clearly a non-vanishing imaginary part in the limit for . The data for is marginal: It is difficult to discern whether or not the average location reaches the origin in the thermodynamic limit. Analysis from a different standpoint will be presented in the next section. The real part of the average location of zeros off the imaginary axis seems to rapidly approach the origin as .
In three dimensions, we tried the fit of the points on the imaginary axis to , because we expect the zeros to reach the origin in the thermodynamic limit due to the existence of the spin glass phase. Solid lines on the left parts in figure 8 represent the best fits to the power function, for and for .
The value for the critical point in three dimensions may not be readily identified with a critical exponent: In the pure system it is known that the edge approaches the origin as . However, in the present random system, we evaluated the average location of the nearest zeros, for which we have no established scaling that relates with critical exponents. In two dimensions our system does not show conventional critical behaviour at , and the exponent in figure 8 would not have direct relevance to critical exponents.
It is not easy to draw a definite conclusion on the behaviour of zeros in the thermodynamic limit in three dimensions since the system size is small. Nevertheless the following story is not inconsistent with our data: At , if the average location of zeros reaches the real axis, it is likely to be only at the origin. At low temperatures, some people may wish to interpret the data that the system-size dependence of average zeros off the imaginary axis suggests the possible stability of the spin glass phase in the external field. If the imaginary part of the average of zeros off the imaginary axis approaches (triangles on the left bottom of figure 8) and the real part approaches a finite value (triangles on the right bottom), the spin glass phase would be stable in the external field. However, the data on the right-bottom panel of figure 8 may instead be extrapolated to , rather than to a finite value, in the thermodynamic limit. If this is indeed the case, it might imply the absence of the AT line  in the three dimensional Ising model. We should anyway be very careful to draw a conclusion from the present data for small systems.
5 Density of zeros on the imaginary axis
In order to analyze the behaviour of zeros in a more detailed way, we study the density of zeros on the imaginary axis in two dimensions. We observe that zeros on the imaginary axis are in close proximity to the real axis in figures 2 to 6. In particular we are interested in whether or not the density shows an essential singularity at the origin as an indicator of the Griffiths singularity. The density of zeros on the imaginary axis will be denoted as , where is defined by for the location of the edge, the nearest zero to the origin on the imaginary axis, of each sample. In order to compare the density of the model with that of the diluted ferromagnet, we also estimated the density for the bond-diluted ferromagnet with the bond probability being 1/2 at , which is believed to show an essential singularity at the origin .
The density of nearest zeros falls very rapidly as shown in figure 9 both for the conventional simple sampling and for the importance sampling (to be explained now). Since we are interested in the behaviour of where its value is extremely small, a special care should be taken to reduce statistical errors in the estimate of . In general, the evaluation of events that occur with exponentially small probabilities is very difficult if one employs a simple-sampling procedure. We therefore used the importance-sampling Monte Carlo algorithm as follows [6, 11, 12, 13, 14].
The basic idea is to generate bond configurations according to the Metropolis algorithm. The new bond configuration is generated from the current configuration with the probability
where is the location of the edge of a given bond configuration and is a guiding function. The guiding function is ideally equal to because equation (6) then guarantees that the stationary distribution of the Markov chain with equation (6) becomes the inverse of the guiding function, . Therefore, configurations with smaller are generated with higher probabilities. The idea is that the new bond configuration with a small probability is encouraged to be chosen. In order to proceed without the prior knowledge of , we first choose an appropriate function as and try an incremental improvement as follows.
After an update, the nearest zero is calculated from the new bond configuration , and the histogram is incremented at this as . Then the histogram will reach the stationary distribution, multiplied by the true distribution ,
Thus, we can estimate as
In our calculations, the new bond configuration is generated by flipping a single bond out of the current bond configuration . The initial guiding function is generated from the simple-sampling algorithm for about independent bond configurations. We repeat the above process to update until becomes nearly flat for the desired range, and finally we suppose that is close enough to . The estimated autocorrelation time of our data is about a few Monte Carlo steps. However, we employ the data taken at every step because we observed no visible differences in the log scale by the change of steps to take data.
We have obtained the data for of the model and the diluted ferromagnet. Figure 10 shows the estimated density of the edge for in two dimensions at for the model and for the diluted ferromagnet. These temperatures are below the transition points of the corresponding pure ferromagnets and therefore the systems lie in the Griffiths phase if any. With an essential singularity at the origin in mind, we analyze the data by the following function in consideration of the finite-size effect
where is an adjustable parameter expected to vanish in thermodynamic limit 222A fit to a power, , gave -dependent , an inconsistent result.. The best fits to equation (9) have been obtained by choosing for the model and for the diluted ferromagnet as shown in figure 11 333Error estimates in these exponents are chosen to be on the safe side such that the error bars in figure 11 are well included by the margin of several times.. For the infinite system, is observed to vanish and the density becomes
This function has an essential singularity at . If , it is the same function as suggested by Bray and Huifang for the diluted ferromagnet . It is remarkable that the analytical result of  has been confirmed numerically and, in addition, a similar result has been observed for the model. This fact suggests that a common physics is likely to underlie the behaviour, which supports the existence of the Griffiths phase also in the model although connected clusters are not trivially defined in the spin glass system.
In this paper we have evaluated the zeros of the partition function of the symmetric model in two and three dimensions by the transfer matrix method. The results revealed several outstanding features of the distribution of zeros in the complex field plane.
On the qualitative aspects, we found that the distribution has a real part on the complex field plane. Nearest zeros to the origin lie on the imaginary axis. This edge moves toward the real axis as the system size increases in both dimensions. For the problem of the existence of the Griffiths singularity, we estimated the density of the nearest zero on the imaginary axis for each sample using the importance-sampling Monte Carlo algorithm. The idea of this algorithm is to enhance the probability of events with low probabilities. As a consequence, we could observe the density for both the model and the diluted ferromagnet for very small and extremely small . This function has an essential singularity at : Thus this result for the density of the model is compatible with the existence of the Griffiths singularity similarly to the diluted ferromagnet albeit through a slightly different mechanism as indicated by the difference between for the model and for the diluted ferromagnet. This is the first evidence for a Griffiths singularity in spin glass systems in equilibrium. Further work is in progress to clarify the temperature dependence of the density and other properties using the method of  for example, which may give us useful information on the phase transition.
Appendix A The nearest zero of the pure system
We can confirm the precision of our calculations from the computation of zeros of pure systems with uniform , for which all zeros are on the imaginary axis by the circle theorem. It is important to check the precision in these calculations, because we use the canonical transfer matrix which treats as numerical values in the partition function . In tables 1 and 2 we list the location of the edge, nearest zero, as functions of the system size in two and three dimensions, respectively. These results for have been obtained by the method of linear equations.
To analyze the data, it is useful to remember that the angle is scaled as 
where is defined as
The exponent is also expanded as
We extrapolate to using the BST algorithm with . The data in table 1 yields and in two dimensions (figure 13), very close to equation (11). In three dimensions, it is not trivial to define the linear size , since the systems are not regular hexahedrons. Figure 13 presents the analyses by the fitting function . Then we obtained and . These results agree well with equation (11).
For further analysis we tried to fit the data at to to obtain and . We thus conclude
with negligibly small higher order corrections. These results give us justifications of our precision and of the use of as the linear size in three dimensions.
-  Lee T D and Yang C N 1952 Phys. Rev.87 410
-  Ozeki Y and Nishimori H 1988 J. Phys. Soc. Japan57 1087
-  Griffiths R B 1969 Phys. Rev. Lett.23 17
-  Bray A J and Huifang D 1989 Phys. Rev.B 40 6980
-  Bray A J and Moore M A 1982 J. Phys. C: Solid State Phys.15 L765
-  Hukushima K and Iba Y 2008 J. Phys.: Conference Series 95 012005
-  Randeria M, Sethna J P and Palmer R G 1985 Phys. Rev. Lett.54 1321
-  Bhanot G and Lacki J 1993 J. Stat. Phys. 71 259
-  Damgaard P H and Lacki J 1995 Int. J. Modern Phys. C 6 819
-  Fisher M E 1964 The Nature of Critical Points (Lectures in Theoretical Physics vol 7C) ed Brittin W E and Dunham L G (the University of Colorado Press, Boulder, 1965)
-  Hartmann A K 2002 Phys. Rev.E 65 056102
-  Körner M, Katzgraber H G and Hartmann A K 2006 J. Stat. Mech. P04005
-  Monthus C and Garel T 2006 Phys. Rev.E 74 051109
-  Iba Y and Hukushima K 2007 Preprint arXiv:0709.2578
-  Binder K 1972 Physica 62 508
-  Morgenstern I and Binder K 1979 Phys. Rev. Lett.43 1615
-  Creswick R J 1995 Phys. Rev.E 52 R5735
-  Kawashima N and Rieger H 2004 Frustrated Spin Systems ed Diep H T (Singapore: World Scientific) Chapter 9
-  Creswick R J and Kim S -Y 1997 Phys. Rev.E 56 2418
-  de Almeida J R L and Thouless D J 1978 J. Phys. A: Math. Gen.11 983
-  Janke W, Johnston D A and Kenna R 2004 Nucl. Phys.B 682 618
-  Blöte H W J, Luijiten E and Heringa J R 1995 J. Phys. A: Math. Gen.28 6289
-  Henkel M and Schutz G 1988 J. Phys. A: Math. Gen.21 2617