Approaching equilibrium and the distribution of clusters
We investigate the approach to stable and metastable equilibrium in Ising models using a cluster representation. The distribution of nucleation times is determined using the Metropolis algorithm and the corresponding model using Langevin dynamics. We find that the nucleation rate is suppressed at early times even after global variables such as the magnetization and energy have apparently reached their time independent values. The mean number of clusters whose size is comparable to the size of the nucleating droplet becomes time independent at about the same time that the nucleation rate reaches its constant value. We also find subtle structural differences between the nucleating droplets formed before and after apparent metastable equilibrium has been established.
Understanding nucleation is important in fields as diverse as materials science, biological physics, and meteorology (1); (2); (3); (4); (5); (6); (7); (8); (9). Fundamental progress was made when Gibbs assumed that the nucleating droplet can be considered to be a fluctuation about metastable equilibrium, and hence the probability of a nucleating droplet is independent of time (10). Langer (11) has shown that the probability of a nucleating droplet can be related to the analytic continuation of the stable state free energy in the limit that the metastable state lifetime approaches infinity. Hence the assumption by Gibbs is valid in this limit. It has also been shown that the Gibbs assumption is correct in systems for which the interaction range (12); (13).
For metastable states with finite lifetimes equilibrium is never reached because a large enough fluctuation would initiate the transformation to the stable state. However, if the probability of such a fluctuation is sufficiently small, it is possible that systems investigated by simulations and experiments can be well approximated as being in equilibrium. Hence, for metastable lifetimes that are very long, we expect the Gibbs assumption to be a good approximation.
In practice, nucleation is not usually observed when the lifetime of the metastable state is very long. Processes such as alloy formation, decay of the false vacuum, and protein crystallization generally occur during a continuous quench of a control parameter such as the temperature. It is natural to ask if the nucleation process that is observed occurs when the system can be reasonably approximated by one in metastable equilibrium. If so, the nucleation rate will be independent of time.
It is usually assumed that metastable equilibrium is a good approximation when the mean value of the order parameter and various global quantities are no longer changing with time. As an example, we consider the nearest-neighbor Ising model on a square lattice and equilibrate the system at temperature in a magnetic field . The relatively small value of the linear dimension was chosen in order to avoid nucleation occurring too quickly. At time the sign of the magnetic field is reversed. In Fig. b we plot the evolution of the magnetization and the energy per spin using the Metropolis algorithm. The solid lines are the fits to an exponential function with the relaxation time . In the following we will measure the time in terms of Monte Carlo steps per spin. A major goal of our work is to address the question, “Can the system be treated as being in metastable equilibrium for ?”
If the nucleation rate is independent of time, the probability of a nucleating droplet occurring at time after the change of magnetic field is an exponentially decreasing function of time. To understand this dependence we divide the time into intervals and write the probability that the system nucleates in a time interval as , where the nucleation rate is a constant. The probability that nucleation occurs in the time interval is given by
If we assume that is small and write , we can write
where is the probability that the system nucleates at a time between and after the change of the magnetic field. In the following we ask if the nucleation rate and the mean values of the order parameter and other thermodynamic quantities become independent of time at approximately the same time after a quench or is the approach to metastable equilibrium more complicated?
In Sec. II we determine the probability distribution of the nucleation times and find that the nucleation rate becomes a constant only after a time that is much longer than the relaxation time of and . In Sec. III we study the microscopic behavior of the system and determine the relaxation time for , the mean number of clusters of size , to approach its equilibrium value (14). Our main result is that is an increasing function of , and the time required for to reach its equilibrium value is the same order of magnitude as for values of comparable to the nucleating droplet. That is, the time for the number of clusters that are the size of the nucleating droplet to reach its equilibrium value is considerably longer than the time for the mean value of the order parameter to become independent of time within the accuracy that we can determine.
In Secs. IV and V we show that there are subtle differences between the structure of the nucleating droplets which occur before and after metastable equilibrium appears to have been achieved. This difference suggests the possibility of finding even greater differences in the nucleating droplets in systems of physical and technological importance. We summarize and discuss our results in Sec. VI. In the Appendix we study the evolution of the clusters after a quench to the critical temperature of the Ising model and again find that that the clusters equilibrate in size order, with the smaller clusters equilibrating first. Hence in principle, an infinite system will never equilibrate. How close to equilibrium a system needs to be and on what spatial scale so that it can be treated by equilibrium methods depends on the physical process of interest.
Ii Distribution of nucleation times
We simulate the Ising model on a square lattice with interaction range with the Hamiltonian
where is the external field. The notation in the first sum means that the distance between spins and is within the interaction range . We studied both nearest-neighbor () and long-range interactions (). The interaction strength is scaled as , where is the number of interaction neighbors per spin. The external field and the temperature are measured in terms of . All of our simulations are at temperature , where is the critical temperature. For the critical temperature is . For the mean field result is a good approximation to the exact value of the critical temperature (21). As discussed in Sec. I the system is equilibrated in a magnetic field . The time corresponds to the time immediately after the magnetic field is reversed.
The clusters in the Ising model are defined rigorously by a mapping of the Ising critical point onto the percolation transition of a properly chosen percolation model (22); (23); (9). Two parallel spins that are within the interaction range are connected only if there is a bond between them. The bonds are assigned with the probability for and near the spinodal, where is the density of the stable spins, and is the inverse temperature. Spins that are connected by bonds form a cluster.
Because the intervention method (15) of identifying the nucleating droplet is time consuming (see Sec. IV), we use a simpler criterion in this section to estimate the nucleation time. We monitor the size of the largest cluster (averaged over 20 bond realizations) and estimate the nucleation time as the time when the largest cluster first reaches a threshold size . The threshold size is chosen so that the largest cluster begins to grow rapidly once its size is greater than or equal to . Because is larger than the actual size of the nucleating droplet, the nucleation time that we estimate by this criterion will be 1 to 2 Monte Carlo steps per spin later than the nucleation time determined by the intervention method. Although the distribution function is shifted to slightly later times, the nucleation rate is found to be insensitive to the choice of the threshold.
Figure b shows for and , where is the probability that nucleation has occurred between time and . The results for were averaged over 5000 runs. The mean size of the nucleating droplet is estimated to be approximately 25 spins for this value of . Note that is an increasing function of for early times, reaches a maximum at , and fits to the expected exponential form for . The fact that falls below the expected exponential for indicates that the nucleation rate is reduced from its equilibrium value and that the system is not in metastable equilibrium. Similar nonequilibrium effects have been observed in Ising-like (16); (17) and continuous systems (18). We conclude that the time for the nucleation rate to become independent of the time after the change of magnetic field is much longer than the relaxation time of the magnetization and energy. We will refer to nucleation that occurs before metastable equilibrium has been reached as transient nucleation.
In order to see if the same qualitative behavior holds near the pseudospinodal, we simulated the long-range Ising model with and . In the mean-field limit the spinodal field is at (for ). A plot of for this system is shown in Fig. a and is seen to have the same qualitative behavior as in Fig. b for ; the relaxation time . In Fig. b the distribution of nucleation times is shown, and we see that does not decay exponentially until . According to Ref. (19), should become comparable to in the limit because the free energy is described only by the magnetization in the mean-field limit. We find that the difference between and is smaller for than for , consistent with Ref. (19).
Iii Relaxation of clusters to metastable equilibrium
Given that there is a significant time delay between the relaxation of the magnetization and the energy and the equilibration of the system as measured by the nucleation rate, it is interesting to monitor the time-dependence of the cluster-size distribution after the reverse of the magnetic field. After the change the system gradually relaxes to metastable equilibrium by forming clusters of spins in the stable direction. How long is required for the number of clusters of size to reach equilibrium? In particular, we are interested in the time required for clusters that are comparable in size to the nucleating droplet.
We first consider and monitor the number of clusters of size at time . To obtain good statistics we chose and averaged over 5000 runs. Figure 4 shows the evolution of , which can be fitted to the exponential form:
We find that for . By doing similar fits for a range of , we find that the time for the mean number of clusters of size to become time independent increases linearly with over the range of that we can simulate (see Fig. b). The extrapolated value of corresponding to the mean size of the nucleating droplet ( spins by direct simulation) is . That is, it takes a time of for the mean number of clusters whose size is the order of the nucleating droplets to become time independent. The time is much longer than the relaxation time of the macroscopic quantities and and is comparable to the time for the nucleation rate to become independent of time.
Because the number of clusters in the nucleating droplet is relatively small for except very close to coexistence (small ), we also consider a long-range Ising model with and (as in Fig. b). The relaxation time of the clusters near the pseudospinodal fits to a power law with (see Fig. b). We know of no theoretical explanation for the qualitatively different dependence f the relaxation time on near coexistence () and near the spinodal (). If we extrapolate to , the approximate size of the nucleating droplet, we find that the equilibration time for clusters of the size of the nucleating droplet is , which is comparable to the time for the nucleation rate to become independent of time.
To determine if our results are affected by finite size effects, we compared the equilibration time of the clusters for lattices with linear dimension and . The equilibration times of the clusters were found to be unaffected.
Iv Structure of the nucleating droplet
Because nucleation can occur both before and after the system is in metastable equilibrium, we ask if there are any structural differences between the nucleating droplets formed in these two cases. To answer this question, we determine the nature of the nucleating droplets for the one-dimensional (1D) Ising model where we can make (and hence the size of the nucleating droplets) large enough so that the structure of the nucleating droplets is well defined. In the following we take , , and . The relaxation time for is , and the time for the distribution of nucleation times to reach equilibrium is .
We use the intervention method to identify nucleation (15). To implement this method, we choose a time at which a nucleating droplet might exist and make many copies of the system. Each copy is restarted using a different random number seed. The idea is to determine if the largest cluster in each of the copies grows in approximately the same place at about the same time. If the percentage of copies that grow is greater than 50%, the nucleating droplet is already in the growth phase; if it is less than 50%, the time chosen is earlier than nucleation. We used a total of 20 trials to make this determination.
Our procedure is to observe the system for a time after the intervention and determine if the size of the largest cluster exceeds the threshold size at approximately the same location. To ensure that the largest cluster at is the same cluster as the original one, we require that the center of mass of the largest cluster be within a distance of the largest cluster in the original configuration. If these conditions are satisfied, the nucleating droplet is said to grow. We choose , , and . (In comparison, the size of the nucleating droplet for the particular run that we will discuss is spins.)
There is some ambiguity in our identification of the nucleation time because the saddle point parameter is large but finite (9). This ambiguity manifests itself in the somewhat arbitrary choices of the parameters , , and . We tried different values for , , and and found that our results depend more strongly on the value of the parameter than on the values of and . If we take , the nucleating droplets almost always occur one to two Monte Carlo steps per spin later than for . The reason is that the linear size of the nucleating droplet is typically 6 to , and its center of mass might shift more than during the time . If such a shift occurs, a cluster that would be said to grow for would not be counted as such because it did not satisfy the center of mass criterion. This shift causes an overestimate of the time of the nucleating droplet. A reasonable choice of is 20% to 40% of the linear size of the nucleating droplet. The choice of parameters is particularly important here because the rate of growth of the transient nucleating droplets is slower than the growth rate of droplets formed after metastable equilibrium has been reached. Hence, we have to identify the nucleating droplet as carefully as possible.
Because nucleation studies are computationally intensive, we used a novel algorithm for simulating Ising models with a uniform long-range interaction (20). The algorithm uses a hierarchical data structure to store the magnetization at many length scales, and can find the energy cost of flipping a spin in time , rather than the usual time , where is the spatial dimension.
Figure 6 shows the fraction of copies for which the largest cluster grows as a function of the intervention time. For this particular run the nucleating droplet is found to occur at .
We simulated 100 systems in which nucleation occurred before global quantities such as became independent of time, , and 100 systems for which nucleation occurred after the nucleation rate became time independent (). We found that the mean size of the nucleating droplet for is with a standard deviation of in comparison to the mean size of the nucleating droplet for of and . That is, the nucleating droplets formed before metastable equilibrium has been reached are somewhat smaller.
We introduce the cluster profile to characterize the shape of the largest cluster at the time of nucleation. For a particular bond realization a spin that is in the stable direction might or might not be a part the largest cluster due to the probabilistic nature of the bonds. For this reason bond averaging is implemented by placing 100 independent sets of bonds between spins with probability in the stable direction. The clusters are identified for each set of bonds, and the probability that spin is in the largest cluster is determined. The values of for the spins in a particular bin are then averaged using a bin width equal to . This mean value of is associated with . Note that the spins that point in the unstable direction are omitted in this procedure. The mean cluster profile is found by translating the peak position of each droplet to the origin.
Figure a shows the mean cluster profile formed after metastable equilibrium has been established (). The position is measured in units of . For comparison we fit to the form
with , and by construction. In Fig. b we show a comparison of to the Gaussian form with and . Note that Eq. (5) gives a better fit than a Gaussian, which underestimates the peak at and the wings. Although Unger and Klein (12) derived Eq. (5) for the magnetization saddle point profile, we see that this form also provides a good description of the cluster profile.
A comparison of the cluster profiles formed before and after metastable equilibrium is shown in Fig. 8. Although both profiles are consistent with the form in Eq. (5), the transient nucleating droplets are more compact, in agreement with the predictions in Ref. (19).
We also directly coarse grained the spins at the time of nucleation to obtain the density profile of the coarse-grained magnetization (see Fig. a). The agreement between the simulation and analytical results (24) are impressive, especially considering that the analytical form is valid only in the limit . The same qualitative differences between the nucleating droplets that occur before and after metastable equilibrium is found (see Fig. b), although the magnetization density profile is much noisier than that based on the cluster analysis.
V Langevin simulations
It is interesting to compare the results for the Ising model and the Langevin dynamics of the model. One advantage of studying the Langevin dynamics of the theory is that it enables the efficient simulation of systems with a very large interaction range . If all lengths are scaled by a large value of , the effective magnitude of the noise decreases, making faster simulations possible.
The coarse grained Hamiltonian analogous to the 1D ferromagnetic Ising model with long-range interactions in an external field can be expressed as
where is the coarse-grained magnetization. A dynamics consistent with this Hamiltonian is given by,
where is the mobility and represents zero-mean Gaussian noise with .
For nucleation near the spinodal the potential has a metastable well only for . The magnitude of and at the spinodal are given by and , and are found by setting . The distance from the spinodal is characterized by the parameter . For , the bottom of the metastable well is near , specifically .
The stationary solutions of the dynamics are found by setting . Besides the two uniform solutions corresponding to the minima in , there is a single nonuniform solution which approximates the nucleating droplet profile when the nucleation barrier is large. When , the profile of the nucleating droplet is described by Eq. (5) with , , and (12).
where is replaced by its central difference approximation. Numerical stability requires that , but it is often desirable to choose even smaller for accuracy.
As for the Ising simulations, we first prepare an equilibrated system with in the stable well corresponding to the direction of the external field . At the external field is reversed so that the system relaxes to metastable equilibrium. We choose , , , , and . The scaled length of the system is chosen to be . We choose to be large so that, on length scales of , the metastable fluctuates near its equilibrium value .
After nucleation occurs will rapidly grow toward the stable well. To determine the distribution of nucleation times, we assume that when the value of the field in any bin reaches , nucleation has occurred. This relatively crude criterion is sufficient for determining the distribution of nucleation times if we assume that the time difference between the nucleation event and its later detection takes a consistent value between runs.
Figure 10 compares the distribution of 50,000 nucleation times for systems with and with and . The distribution shows the same qualitative behavior as found in the Metropolis simulations of the Ising model (see Fig. b). For example, the distribution of nucleation times is not exponential for early times after the quench. As expected, the nucleation rate decreases as increases. Smaller values of and give similar results for the distribution.
To find the droplet profiles, we need to identify the time of nucleation more precisely. The intervention criterion, which was applied in Sec. IV, is one possible method. In the Langevin context we can employ a simpler criterion: nucleation is considered to have occurred if decays to the saddle-point profile (given by Eq. (5) for ) when is evolved using noiseless dynamics (26); (19). For fixed these two criteria agree in the limit, but can give different results for finite (27).
In Fig. 11 we plot the average of 1,000 density profiles of the nucleating droplets formed after metastable equilibrium has been established for and . Note that there are noticeable deviations of the averaged profiles from the theoretical prediction in Eq. (5), but the deviation is less for . The deviation is due to the fact that the bottom of the free energy well in the metastable state is skewed; a similar deviation was also observed in the Ising model. We also note that the individual nucleating droplets look much different from their average. It is expected that as increases, the profiles of the individual nucleating droplets will converge to the form given by Eq. (5).
In Fig. 12 we compare the average of 1,000 density profiles of nucleating droplets before and after metastable equilibrium has been established. As for the Ising model, there are subtle differences consistent with the predictions of Ref. (19). The transient droplets have slightly lower background magnetization and compensate by being denser and more compact.
Although the time-independence of the mean values of macroscopic quantities such as the magnetization and the energy is often used as an indicator of metastable equilibrium, we find that the observed relaxation time of the clusters is much longer for sizes comparable to the nucleating droplet. This longer relaxation time explains the measured non-constant nucleation rate even when global quantities such as the magnetization appear to be stationary. By identifying the nucleating droplets in the one-dimensional long-range Ising model and the Langevin equation, we find structural differences between the nucleating droplets which occur before and after metastable equilibrium has been reached. Our results suggest that using global quantities as indicators for metastable equilibrium may not be appropriate in general, and distinguishing between equilibrium and transient nucleation is important in studying the structure of nucleating droplets. Further studies of transient nucleation in continuous models of more realistic systems would be of interesting and practical importance.
Finally, we note a subtle implication of our results. For a system to be truly in equilibrium would require that the mean number of clusters of all sizes be independent of time. The larger the cluster, the longer the time that would be required for the mean number to become time independent. Hence, the bigger the system, the longer the time that would be required for the system to reach equilibrium. Given that the system is never truly in metastable equilibrium so that the ideas of Gibbs, Langer, and others are never exactly applicable, when is the system close enough to equilibrium so that any possible simulation or experiment cannot detect the difference? We have found that the magnetization and energy are not sufficient indicators for nucleation and that the answer depends on the process being studied. For nucleation the equilibration of the number of clusters whose size is comparable to the size of the nucleating droplet is the relevant indicator.
Appendix A Relaxation of clusters at the critical temperature
Accurate determinations of the dynamical critical exponent have been found from the relaxation of the magnetization and energy at the critical temperature. In the following we take a closer look at the relaxation of the Ising model by studying the approach to equilibrium of the distribution of clusters of various sizes.
We consider the Ising model on a square lattice with . The system is initially equilibrated at either zero temperature (all spins up) or at , and then instantaneously quenched to the critical temperature . The Metropolis algorithm is used.
As a check on our results we first determine starting from . Scaling arguments suggest that approaches its equilibrium value as (28)
where the static critical exponents are and for finite and and in the mean-field limit. The fit of our results in Fig. b to Eq. (9) yields the estimate for and for , which are consistent with previous results (29); (30). Note that no time scale is associated with the evolution of .
We next determined , the number of clusters of size at time after the temperature quench. Because all the spins are up at , the number of (down) clusters of size begins at zero and increases to its (apparent) equilibrium value . The value of the latter depends on the size of the system.
Figure 14 shows the evolution of clusters of size for one run. Because we know of no argument for the time dependence of except in the mean-field limit (30), we have to rely on empirical fits. We find that the time-dependence of can be fitted to the sum of two exponentials,
where , , , and are parameters to be fitted with .
Figure a shows the relaxation time as a function of for at starting from . Note that the bigger the cluster, the longer it takes to reach its equilibrium distribution. That is, small clusters form first, and larger clusters are formed by the merging of smaller ones. The -dependence of can be approximately fitted to a power law with the exponent 0.4.
To prepare a configuration at , the system is randomized with approximately half of the spins up and half of the spins down. The temperature is instantaneously changed to . As before, we focus on the relaxation of down spin clusters. In contrast to the case, the evolution of the clusters falls into three classes (see Fig. b). For small clusters (), monotonically decreases to its equilibrium value. This behavior occurs because the initial random configuration has an abundance of small clusters so that lowering the temperature causes the small clusters to merge to form bigger ones. For intermediate size clusters (), first increases and then decreases to its equilibrium value. The initial growth is due to the rapid coalescence of smaller clusters to form intermediate ones. After there are enough intermediate clusters, they slowly coalesce to form bigger clusters, which causes the decrease. For clusters with , slowly increases to its equilibrium value. The range of sizes for these different classes of behavior depends on the system size. In all three cases can be fitted to the sum of two exponentials. One of the two coefficients is negative for for which overshoots its equilibrium value. The relaxation time is plotted in Fig. b as a function of .
Acknowledgements.We thank Aaron O. Schweiger for very useful discussions. Bill Klein acknowledges the support of Department of Energy grant # DE-FG02-95ER14498 and Kipton Barros was supported in part by the National Science Foundation grant # DGE-0221680. Hui Wang was supported in part by NSF grant # DUE-0442581. The simulations at Clark University were done with the partial support of NSF grant # DBI-0320875.
- Dimo Kashchiev, Nucleation: Basic Theory with Applications (Butterworths-Heinemann, Oxford, 2000).
- N. E. Chayen, E. Saridakis, R. El-Bahar, and Y. J. Nemirovsky, Mol. Biol. 312 (4), 591 (2001).
- N. Wang, Y. H. Tang, Y. F. Zhang, and C. S. Lee, Phys. Rev. B 58, R16 024 (1998).
- S. Auer and D. Frenkel, Nature 413, 711 (2001).
- N. Delgehyr, J. Sillibourne, and M. Bornens, J. Cell Science 118, 1565 (2005).
- E. Pechkova and C. Nicolini, J. Cell Biochem. 85 (2), 243 (2002).
- A. R. Fersht, Proc. Natl. Acad. Sci. U. S. A. 92, 10869 (1995).
- D. Stauffer, A. Coniglio, and D. W. Heermann, Phys. Rev. Lett. 49, 1299 (1982).
- W. Klein, H. Gould, N. Gulbahce, J. B. Rundle, and K. Tiampo, Phys. Rev. E 75, 031114 (2007).
- J. D. Gunton, M. san Miguel, and P. Sahni, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic Press, New York, 1983), Vol. 8.
- J. S. Langer, Ann. Phys. (NY) 41, 108 (1967).
- C. Unger and W. Klein, Phys. Rev. B 29, 2698 (1984).
- More precisely, the Gibbs assumption is correct in systems for which the interaction range and which are not too close to the spinodal. See Ref. (9) for more details.
- We assume that when the nucleation rate and mean number of clusters of a given size become apparently independent of time that they have reached their equilibrium values.
- L. Monette, W. Klein, and M. Zuckermann, J. Stat. Phys. 66, 117 (1992).
- D. W. Heermann and C. Cordeiro, Int. J. Mod. Phys. 13, 1419 (2003).
- K. Brendel, G. T. Barkema, and H. van Beijeren, Phys. Rev. B 71, 031601 (2005).
- H. Huitema, J. van der Eerden, J. Janssen, and H. Human, Phys. Rev. B 62, 14690 (2000).
- A. O. Schweiger, K. Barros, and W. Klein, Phys. Rev. E 75, 031102 (2007).
- K. Barros, manuscript in preparation.
- E. Luijten, H. W. J. Blöte, and K. Binder, Phys. Rev. E 54, 4626 (1996).
- W. Klein in Computer Simulation Studies in Condensed Matter Physics III, edited by D. P. Landau, K. K. Mon, and H. B. Schuttler (Springer-Verlag, Berlin, Heidelberg, 1991).
- A. Coniglio and W. Klein, J. Phys. A 13, 2775 (1980).
- The density profile of the nucleating droplet of the Ising model has been calculated analytically in the limit (K. Barros, unpublished). The result is consistent with the form in Eq. (5) with , , and , where is the magnitude of at the spinodal. From this analytical solution, the calculated parameters are found to be , , which are very close to the values fitted to the simulation data, , , .
- J. G. Gaines, in Stochastic Partial Differential Equations, edited by A. M. Etheridge (Cambridge University Press, Cambridge, 1995), pp. 55–71.
- A. Roy, J. M. Rickman, J. D. Gunton, and K. R. Elder, Phys. Rev. E 57, 2610 (1998).
- Consider a fluctuation that decays to the metastable phase under noiseless dynamics. To perform the intervention method we make many copies of the configuration and examine the percentage that grow after a given waiting time. Although the expected drift is a decay to the metastable phase, every copy has time to sample a path in configuration space. It is possible that during this waiting time the majority of copies discover and grow toward the stable phase, contradicting the result from the zero-noise criterion. However, for the sampling path will be dominated by the drift term and the two nucleation criteria agree.
- M. Suzuki. Phys. Lett. A 58, 435 (1976) and M. Suzuki. Prog. Theor. Phys. 58, 1142 (1977). See also A. Linke, D. W. Heermann, P. Altevogt, and M. Siegert, Physica A 222, 205 (1995).
- M. Nightingale and H. Blöte, Phys. Rev. B 62, 1089 (2000).
- L. Colonna-Romano, A. I. Mel’cuk, H. Gould, and W. Klein, Physica A 209, 396 (1994).