What do Monte Carlo simulations tell us about compressible Ising models?

What do Monte Carlo simulations tell us about compressible Ising models?

D. P. Landau

Center for Simulational Physics, The University of Georgia, Athens, GA 30622, U.S.A.

ABSTRACT

Monte Carlo simulations have begun to illuminate the nature of phase transitions and universality classes for compressible Ising models. A comprehensive analysis of a Landau-Ginsburg Wilson hamiltonian for systems with elastic degrees of freedom predicts that there should be four cases with different behavior, depending upon symmetries and thermodynamic constraints. We shall describe the results of careful Monte Carlo simulations for a simple compressible Ising model that can be easily modified to correspond to each of the four cases.

Keywords: Compressible Ising model; Phase transition; Monte Carlo; Finite size scaling

I. INTRODUCTION

The study of critical phenomena in magnetic systems is a mature endeavor with quite substantial results from both theory and experiment. In addition, high resolution Monte Carlo studies of Ising models have provided quite precise values for critical temperatures and exponents for several different three dimensional Ising models. In many cases, however, the experimental data deviate from theoretical predictions close to T_c. One possible explanation is that real materials have elastic, not rigid, lattices. High resolution measurements of the specific heat of DAG (dysprosium aluminum garnet) [1] show that quite close to the critical point the apparently diverging peak becomes rounded, see Fig. 1, in spite of the fact that the single crystals that were used were of exceedingly high quality. Is the elastic nature of the lattice the culprit? The question of what happens when the "lattice" is allowed to be compressible has long been the subject of theoretical scrutiny (the theoretical background is reviewed elsewhere [2-4]). Of course, the equivalence between the Ising model and the binary alloy model adds to the general interest in the problem.

One simple example is the Si/Ge alloy for which the covalent interactions give rise to strongly directional nearest neighbor bonds that dominate the behavior. (Ising spins s_i = ±1 are then equivalent to Si or Ge atoms, respectively; and hence, the concentration of Ge atoms plays the role of the magnetization in the corresponding magnetic system. In the same vein, the chemical potential is equivalent to the magnetic field in the magnetic interpretation.) The appropriate model contains both two-body and three-body interactions so that both the bond lengths and bond angles are somewhat constrained. The particular values of the constants were determined by fitting to the properties of Si/Ge alloys, but the goal here is actually to determine the behavior of a generic compressible Ising model. Of course, it remains to be seen if this type of model fully encompasses the range of physical behavior characteristic of "compressible Ising models" or if the problem is more subtle.

A comprehensive analysis of a Landau-Ginsburg Wilson hamiltonian for systems with elastic degrees of freedom determined [3,4] that there were four distinct cases that would exhibit quite different behavior. These depend upon symmetries (e.g. the coupling between the elastic and magnetic degrees of freedom) as well as the thermodynamic constraints and are listed in Table I.

With the initiation of extensive Monte Carlo simulations of distortable Ising nets with elastic couplings the "modern", comprehensive approach to the study of problems in physics becomes achievable (see Fig. 2). Thus, the purpose of this manuscript is to present an overview of the state of our understanding of compressible Ising models resulting from more than a decade of simulations performed in the Center for Simulational Physics.

II. MODEL AND METHOD

In these studies the Ising model-binary alloy equivalence is used, and for the Ising ferromagnet Si/Ge alloys correspond to a physical example. "Spins" are placed on a distortable diamond net with 8 sites per unit cell. The simulation cell contains L ×L ×L unit cells and the boundary conditions were periodic. Multiple hamiltonians were considered but all contained distance dependent two-body and three-body terms. The initial choice was the Keating [7] interatomic potential, whose "stiffness" parameters E and A were determined by the macroscopic elastic constants of the crystal and which has been used to describe the structural properties of mixed systems [2].

with

(Kronecker's d), where µ_A and µ_B are chemical potentials

where the e(S_i,S_j) are chemical binding energies

where E(S_i,S_j) are bond btiffnesses", R₀(S_i,S_j) ideal lattice neighbor distances, and

where A_i,j,k(S_i,S_j,S_k) are "angular stiffnesses". In these equations, i-j denotes a bond between nearest neighbors i and j, while i-j-k denotes the angle of the nearest-neighbor bonds i-j and j-k with vertex at site j. The vector _ij = _i - _j, where _i is the position of site i. Additional degrees of freedom are the linear sizes of the simulation cell L_x, L_y and L_z which fluctuate if the pressure is held fixed, e.g. P = 0. These simulations were repeated and extended [8] using the Stillinger-Weber potential [9]. Also containing two-body and three-body interactions, this potential has the decided advantage that a simple Stillinger-Weber alloy expands upon heating whereas a Keating model shows unphysical shrinkage with increasing temperature.

The Monte Carlo sampling [10] was implemented as follows: For a single particle of type S_i at position _i, a new type is randomly chosen as is a slightly displaced new position , while keeping the other particles and the simulation box dimensions fixed. This random trial move is accepted or rejected according to the usual Metropolis criterion. After all particles have been considered, new box dimensions , and are randomly chosen. The Metropolis acceptance criterion uses

where D is the energy change associated with this global distortion of the system and the latter term describes the change in translational entropy due to a volume change. (Since the total number of particles was constant, although the number of Si and Ge were not individually fixed, the simulations were performed in the semi-grand canonical ensemble.) Constant volume simulations were performed by simply turning off volume changing moves. In such cases the fixed volume was set at about 3/4 of the way between the volumes for pure Si and pure Ge at T = 0. For some simulations parallel tempering [11,12] was implemented to overcome thermal slugishness. Intially runs of only 10⁴ MCS were used to determine parameter ranges of interest, but runs of length 10⁵ - 10⁷ MCS or more were used for more serious data taking.

Because of the Ising model - binary alloy model equivalence, either the magnetic field H or the chemical potential difference Dm can be used for presenting results. Similarly, either the magnetization M or the Ge concentration c_Ge provide a measure of the order parameter. However, it should be noted that the elastic degrees of freedom result in a non-equivalence of the constant-M and constant-H ensembles. This has to do with the fact that the usual grand-canonical particle bath has the properties of a fluid, while it would have to have the properties of a solid (i.e., in particular, apply coherency stresses) in order to ensure equivalence. This difference was first noted by Vandeworp and Newman [13].

Monte Carlo data were obtained by sweeping the chemical potential at low temperature or by sweeping the temperature at fixed chemical potential for higher temperatures. Because of pronounced hysteresis at low temperatures, thermodynamic integration was used to find the intersection of the free energies and hence the location of the transition. At higher temperatures histogram reweighting [14] was used to locate the transition at, or near, the critical point. Of course, in a finite system transitions are rounded and shifted [15], so finite size scaling was then used to extrapolate to the thermodynamic limit.

III. RESULTS

For the compressible ferromagnet at constant pressure, the transition at low temperatures was strongly first order and free energy integrations were needed to locate the transition. At higher temperatures the hysteresis was much smaller and a "law of equal areas" could be used to estimate the location of the transition. An example is given in Fig. 3, where the data are shown using binary alloy variables. At intermediate temperatures the hysteresis shown in the inset can only be seen using a fine scale observation. Since the spin-up / spin-down symmetry of the rigid Ising model is no longer valid, the search for the critical point necessitated scanning the two-dimensional chemical potential-temperature space. With the aid of histogram reweighting, the critical point was located and excellent finite size scaling was obtained using mean-field exponents. Clear, systematic deviations were found for other choices of exponents so they could be excluded. The 4th order cumulants of the order parameter were also found to cross at a different value than the rigid Ising one (U₄^* ~ 0.47), but the crossing point matched the mean field value quite well. See Fig. 4. The conclusion that the critical behavior was mean field-like was later found in the analysis of the Landau-Ginsburg-Wilson Hamiltonian mentioned earlier.

Qualitatively similar results were obtained for the compressible Ising ferromagnet when the Stillinger-Weber potential is substituted for the Keating potential [8]. Although the scales of interest for the temperature and chemical potential change, the asymptotic behavior of the normalized thermodynamic quantities as the critical point is approached does not. Finite scaling of the reduced fourth order cumulant, shown in Fig. 5, demonstrates how well mean-field exponents work in this situation as well.

For the constant volume case, the analysis of the Landau-Ginsburg-Wilson Hamiltonian [3] suggests that two transitions should occur at low temperatures: one from a Si-rich state to one that is about equal in Si and Ge concentration, followed by a second transition to a state that was almost pure Ge. In contrast, however, the Monte Carlo data obtained by sweeping chemical potential at low temperatures revealed a relatively smooth increase of the Ge-concentration with varying chemical potential and no apparent transition. However, a more careful, high resolution study showed that the low temperature data indeed indicated very slight, smeared out hysteresis loops for the range of chemical potential between ~ 0.45-0.49. Qualitatively similar hysteresis could be seen for both "high" and "low" values of chemical potential. It was reproducible, was rather insensitive to the length of the runs, and was matched by relatively similar hysteresis in the internal energy as well. The hysteresis slowly diminished as the temperature increased, and at substantially higher temperature it disappeared completely for the concentration within the range of sizes that could be studied. Even so, an examination of the distribution of nearest neighbors of each species showed that a double peaked distribution remained below some chemical potential dependent temperature. It thus appeared as though there was a phase transition separating the low temperature, medium concentration and the disordered states over a wide range of chemical potential (see Fig. 6) but the nature of the transition was perplexing.

Ising model critical exponents, for the compressible Ising antiferromagnet at constant volume [19].

These results could only be understood by augmenting the numerical data with "snapshots" of the system. The visualization program AViz, developed at the Technion [17,18], was used to produce many views of the system for different conditions. Typical snapshots of the system (see e.g. Fig. 7) showed, unexpectedly, that Si rich regions formed slabs of approximately fixed thickness. Instead of the thickness growing with enhanced concentration, additional slabs would simply form. This behavior became increasingly distinct for larger systems. Particularly long time scales were associated with the formation or dissolution of these slabs as the temperature or chemical potential were swept, and this property gave rise to the rather strange "smeared" hysteresis loops. Forming the planes took longer and thus produced the asymmetry in the hysteresis. At higher temperatures the formation of the slabs seemed to indicate the transition to an "ordered" state. Moreover, the presence of the slabs introduces a special length scale into the problem that had not been included in the LGW Hamiltonian study. The resultant phase diagram, depicted in Fig. 6, shows that instead of the two predicted phase lines a single, closed phase boundary is present. The range of coexistence along the 1st order phase boundary, as determined from the hysteresis, is quite narrow, and for the lattice sizes that were accessible it seemed to disappear completely at high temperature. The phase boundary appears to be first order everywhere except at the single point where the derivative dT/dH = 0. Although it is now clear what went wrong in the analysis of this case in Ref. 3 (incorrect assumption of a finite interfacial tension), the physical mechanism which produces the observed slab formation still needs a clearer explanation.

For Monte Carlo simulations of compressible Ising antiferromagnets, essentially the same model (distortable diamond net with a Stillinger-Weber potential) was used. This was done so that as few differences as possible were introduced relative to the ferromagnetic model; however, it was necessary to change the sign and magnitude of the interaction between unlike species to produce an antiferromagnetic ground state. (Of course, in this case the resultant alloy is no longer related in any way to Si/Ge and must be regarded merely as a model of purely theoretical interest.) The phase boundaries were first estimated for modest size lattices, and then more careful runs were made for specific values of chemical potential or temperature. At both constant volume and constant pressure the resultant phase boundaries separating the ordered (antiferromagnetic) and disordered phases were 2nd order. Careful finite size analyses of the critical behavior revealed rigid Ising exponents. For example, Fig. 8 shows a finite size scaling plot, made with rigid

Quite similar behavior was found for the compressible Ising antiferromagnet at constant pressure [4]. Furthermore, as shown in Fig. 9, cumulant crossings for both cases occurred at the value found previously for a rigid Ising model. Since no indication of crossover to any other kind of behavior was seen, substantially larger systems would probably be required to see some other possible kind of asymptotic behavior, if it exists! In such a case the computational resources needed for extending studies of the present model would be orders of magnitude beyond the capability of current computer resources. Thus, this remains a challenge for future generations of students.

IV. CONCLUSIONS

Monte Carlo studies of Ising models on a distortable diamond net have provided a test of the "global" nature of theoretical predictions. For the ferromagnet at constant pressure the critical behavior is mean-field-like. For the constant volume forromagnet a closed first order line was found instead of the two predicted first order lines, terminating in critical points. In the "ordered" phase, the less favorable species forms slabs of approximately fixed thickness, the number of which increases with increasing concentration. For both constant pressure and constant volume the compressible Ising antiferromagnet showed a closed, Ising like 2nd order phase boundary separating the ordered and disordered states. Of course, these studies were restricted to a single type of compressible Ising model, and the general classifications of such models may be more complex. Much remains to be done.

Acknowledgments

We thank B. Düenweg, M. Laradji, F. Tavazza, L. Cannavacciuolo, J. Adler, and X. Zhu for extremely fruitful collaborations. This work was supported by NSF grant No DMR-0341874.

[1] D. P. Landau, B. E. Keen, B. Schneider, and W. P. Wolf, Phys. Rev.3, 2310 (1971).

[2] B. Dünweg and D. P. Landau, Phys. Rev. B 48, 14182 (1993).

[3] B. Dünweg, Computersimulationen xu Phasenübergängen und kritischen Phänomenen, Habilitationsschrift, U. Mainz (2000).

[4] X. Zhu, F. Tavazza, D. P. Landau, and B. Dünweg, Phys. Rev. B 72, 104102 (2005).

[5] This prediction relied on the assumption (Ref. 3) that the interfacial tension between the Si-rich and Ge-rich phases is finite. This has turned out to be incorrect. There is, instead, a complete suppression of capillary waves [6], and also a phase behavior which differs substantially from that assumed in Ref. 3.

[6] B. J. Schulz, B. Dünweg, K. Binder, and M. Müller, Phys. Rev. Lett. 95, 096101 (2005).

[7] P. N. Keating, Phys. Rev. 145, 627 (1966).

[8] M. Laradji, D. P. Landau, and B. Dünweg, Phys. Rev. B 51, 4898 (1995).

[9] F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).

[10] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, 2nd Edition (Cambridge U. Press, Cambridge, 2005).

[11] K. Hukushima and K. Nemoto, Journal Phys. Soc. of Japan 65, 1604 (1996).

[12] R.H. Swendsen and J.S. Wang, Phys. Rev. Lett. 57, 2607 (1986).

[13] E. M. Vandeworp and K. E. Newman, Phys. Rev. B 52, 4086 (1995); E. M. Vandeworp and K. E. Newman, Phys. Rev. B 55, 14222 (1997)

[14] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).

[15] V. Privman (Ed.) Finite Size Scaling and Numerical Simulation of Statistical Systems, World Scientific, Singapore, 1990.

[16] F. Tavazza, D. P. Landau, and J. Adler, Phys. Rev. B 70, 184103 (2004).

[17] J. Adler,A. Hashibon, and G. Wagner, Recent Developments in Computer Simulation Studies in Condensed Matter Physics, XIV, ed. D. P. Landau, S. P. Lewis, and H.-B. Schüttler, p.160 (Springer, Heidelberg, 2002).

[18] J. Adler, Comput. in Sci. and Engineering, 5, 61 (2003)

[19] L. Cannavacciuolo and D. P. Landau, Phys. Rev. B 71, 134104 (2005).

Received on 16 September, 2006

Publication Dates

• Publication in this collection
  16 Oct 2006
• Date of issue
  Sept 2006

History

• Received
  16 Sept 2006