Computer Simulation of the Critical Behavior of 3D Disordered Ising Model
Abstract
The critical behavior of the disordered ferromagnetic Ising model is studied numerically by the Monte Carlo method in a wide range of variation of concentration of nonmagnetic impurity atoms. The temperature dependences of correlation length and magnetic susceptibility are determined for samples with various spin concentrations and various linear sizes. The finitesize scaling technique is used for obtaining scaling functions for these quantities, which exhibit a universal behavior in the critical region; the critical temperatures and static critical exponents are also determined using scaling corrections. On the basis of variation of the scaling functions and values of critical exponents upon a change in the concentration, the conclusion is drawn concerning the existence of two universal classes of the critical behavior of the diluted Ising model with different characteristics for weakly and strongly disordered systems.
pacs:
05.40.a, 64.60.Fr, 75.10.HkI Introduction
Analysis of the critical behavior of disordered systems with quenched structural defects is of considerable theoretical and experimental interest. Most real solids contain quenched structural defects, whose presence can affect the characteristics of the system and may strongly modify the behavior of the systems during phase transitions. This leads to new complex phenomena in structurally disordered systems, which are associated with the effects of an anomalously strong interaction of fluctuations of a number of thermodynamic quantities, when any perturbation introduced by structural defects (even in small concentration) may strongly change the state of the system. The description of such systems requires the development of special analytic and numerical methods.
The following two questions arise when the effect of structural disorder on secondorder phase transitions is investigated: (i) do the critical exponents of a homogeneous magnet change upon its dilution by nonmagnetic impurity atoms? and (ii) if these exponents change, are the new critical exponents universal (i.e., independent of the structural defect concentration up to the percolation threshold)? The answer to the first question was obtained in (1), where it was shown that the critical exponents of systems with quenched structural defects change as compared to their homogeneous analogs if the critical exponent of the heat capacity of a homogeneous system is positive. This criterion is satisfied only by 3D systems whose critical behavior can be described by the Ising model. A large number of publications are devoted to the study of the critical behavior of diluted Isinglike magnets by the renormgroup methods, the numerical Monte Carlo methods, and experimentally (see review (2)). An affirmative answer has been obtained to the question concerning the existence of a new universal class of the critical behavior, which is formed by diluted Isinglike magnet. It remains unclear, however, whether the asymptotic values of critical exponents are independent of the rate of dilution of the system, how the crossover effects change these values, and whether two or more regimes of the critical behavior exist for weakly and strongly disordered systems; these questions are the subjects of heated discussions.
This study is devoted to numerical analysis of the critical behavior of a diluted 3D Ising model in a wide range of concentration of quenched point defects. The fundamental importance of the results of this study is due to stringent requirements to simulation conditions imposed in the course of investigations; the wide range of linear dimensions of lattices () analyzed in this work; the chosen temperature range of simulation close to the critical temperature with , which makes it possible to single out the asymptotic values of characteristics; high statistics used for averaging of thermodynamic and correlation functions over various impurity configurations; the application of finitesize scaling technique (3) for processing the result of simulation, which makes it possible to obtain scaling function for thermodynamic functions apart from their asymptotic values; and application of corrections to scaling for determining the asymptotic values of critical exponents.
Ii Computer simulation technique and results
We consider a model of a disordered spin system in the form of a cubic lattice with linear size under certain boundary conditions. The microscopic Hamiltonian of the disordered Ising model can be written in the form
(1) 
where is the shortrange exchange interaction between spins fixed at the lattice sites and assuming values of . Nonmagnetic impurity atoms form empty sites. In this case, occupation numbers assume the value or and are described by the distribution function
(2) 
with , where  is the concentration of the impurity atom. The impurity is uniformly distributed over the entire system, and its position is fixed in simulation for an individual impurity configuration. We consider here disordered systems with spin concentrations , , , and .
To suppress the effect of critical slowing down and correlation of various spin configurations, we used the singlecluster Wolf algorithm, which is most effective in this respect (4); (5). A Monte Carlo step per spin (MCS) was assumed to correspond to 1020 rotations of a Wolf cluster depending on the linear size of the lattice being simulated, the spin concentration of the system, and the closeness of the temperature to the critical point. The stabilization of thermodynamic equilibrium required Monte Carlo steps, and steps were allotted to statistical averaging of quantities being simulated for a given impurity configuration. To determine the average values of thermodynamic and correlation functions, averaging over various impurity configurations was carried out along with statistical averaging (averaging was carried out over samples for , over samples for , and over samples for and ).
In the course of simulation of various spin systems, correlation length and susceptibility were carried out on lattices with a linear size in accordance with the following relations:
(3) 
where , , and
(4) 
where are the coordinates of the th lattice site; angle brackets indicate statistical averaging over Monte Carlo steps, and the bar indicates averaging over impurity configurations. The temperature dependences and were determined in the temperature interval for samples with and a linear size in the range of . For samples with the remaining spin concentrations, temperatures were chosen in the interval of for values of ranging from to . In computer simulation, the value of for each temperature was limited by the lattice size for which the correlation length and susceptibility of the system attained their asymptotic values.
In accordance with the results obtained in (4); (5) and the results of our investigations, the chosen simulation conditions ensure equilibrium values for measurable thermodynamic quantities for all lattice sizes and spin concentrations studied here since the autocorrelation times for magnetization and energy turn out to be not longer than ten Monte Carlo steps per spin even for chosen temperatures closest to the critical temperature (with allowance for the number of turns of the Wolf cluster taken as a step).
Iii Method of FiniteSize Scaling
It is known that the secondorder phase transition considered here can be manifested only in the thermodynamic limit, when the volume of the system and the number of particles in it tend to infinity. To determine the asymptotic values of thermodynamic quantities exhibiting an anomalous behavior near the critical temperature from their values determined on finite lattices, the concepts of the scaling theory concerning the generalized uniformity of thermodynamic functions in the critical region relative to scale transformations of the system are widely used. These concepts formed the basis of various methods of finitesize scaling. Here, we apply the method proposed in (3) and tested by the authors in analysis of the results of simulation of the critical behavior of 2D and 3D pure Ising models.
The idea of this method (3) is that, in accordance with the scaling theory, the size dependence of a certain thermodynamic quantity defined on a finite lattice in zero magnetic field can be presented in the critical region in the form
(5) 
where is the critical exponent for the thermodynamic quantity . Taking into account the fact that the correlation length in the critical region behaves as , we can write
(6) 
Then expression (5) can be written in the form
(7) 
where the relation between scaling functions and is defined in the form of the relation
(8) 
If correlation length plays the role of quantity , Eq. (7) defines as a function of only one variable . This leads to a relation that makes it possible to find the asymptotic value of any thermodynamic quantity in terms of directly measurable values of and the scaling function of ,
(9) 
where function is defined by the expression
(10) 
Scaling function , defined in the interval , where is the value of the argument independent of in the critical region, must satisfy the following asymptotic conditions: and .
To satisfy the asymptotic conditions, we chose, analogously to (3), the scaling function for susceptibility and correlation length in the form of a polynomial dependence of , as well as of :
(11)  
(12) 
with coefficients selected for each temperature using the least squares method.
Here, we implement the following scheme of finitesize scaling.

For an arbitrary value of in the critical temperature range, the values of and are measured for lattices with increasing size .

The thermodynamic value of quantity is determined as the value of , which is found to be independent of within the error of measurements.

The results of measurements for are processed by the least squares method to determine the corresponding functional form for scaling function .

The procedure is repeated for other values of in the range of .

Averaged scaling function is determined on the basis of functions , determined for various temperatures for a fixed spin concentration of the samples.

The temperature dependence is determined for asymptotic values of the thermodynamic quantity by substituting and into relation (9).
Figures 14 show by way of example the scaling functions for correlation length and susceptibility , obtained for systems with spin concentrations and at various temperatures using the polynomial approximation in variable (Figs. 1, 4) and in variable (Figs. 2, 4). It can be seen from the figures that the scaling functions show a tendency towards a single universal curve for each spin concentration in the entire range of variation of scaling variable .
Figures 5 and 6 show the averaged scaling functions for the correlation length and susceptibility for various spin concentrations , which were obtained using the polynomial approximation in (Fig. 6) and in (Fig. 6). The averaged scaling functions demonstrate a tendency indicating the possible existence of two classes of universal critical behavior for the diluted Ising model with different modes of behavior for weakly () and strongly () disordered systems.
Table 1 contains asymptotic values of and obtained using averaged scaling functions for various temperatures and spin concentrations. The errors in values of and take into account statistical errors in the measured values of and , as well as approximation errors.
p=0.95  T  4.265  4.275  4.280  4.285  4.295  4.315  4.335 

pol  62.44(15)  21.25(5)  17.06(4)  14.41(3)  11.40(2)  8.37(2)  6.76(2)  
exp  62.31(15)  21.19(6)  17.02(4)  14.38(3)  11.38(2)  8.36(2)  6.76(2)  
pol  14467(80)  1748(16)  1130(4)  819(3)  515(2)  282(2)  187(1)  
exp  14359(81)  1724(12)  1126(5)  813(4)  513(2)  281(2)  187(1)  
p=0.80  T  3.51  3.52  3.53  3.54  3.55  3.57  
pol  26.16(9)  16.50(4)  12.51(2)  10.31(2)  8.79(3)  7.01(3)  
exp  26.11(9)  16.46(4)  12.49(3)  10.30(2)  8.76(3)  7.00(3)  
pol  2612(17)  1060(5)  618(2)  424(2)  312(2)  201(2)  
exp  2603(18)  1055(5)  615(2)  423(2)  310(2)  200(1)  
p=0.60  T  2.430  2.435  2.440  2.445  2.450  2.460  
pol  46.03(15)  29.37(7)  22.49(8)  18.33(5)  15.70(4)  12.41(4)  
exp  45.86(13)  29.29(7)  22.40(8)  18.27(5)  15.65(5)  12.37(4)  
pol  7943(55)  3289(17)  1953(15)  1308(7)  967(5)  611(4)  
exp  7881(47)  3268(15)  1937(14)  1298(7)  961(5)  608(4)  
p=0.50  T  1.851  1.854  1.857  1.861  1.865  1.874  
pol  49.55(46)  36.38(12)  29.57(9)  23.86(6)  20.23(4)  15.38(3)  
exp  49.04(38)  36.34(16)  29.48(11)  23.81(7)  20.17(4)  15.35(3)  
pol  9456(356)  5036(36)  3365(26)  2209(12)  1603(6)  939(4)  
exp  9169(165)  5030(43)  3349(25)  2199(12)  1592(6)  934(3) 
Iv Calculation of Critical Characteristics
Asymptotic critical exponent of a thermodynamic quantity is described by the expression
(13) 
where and are the critical amplitudes above and below the critical temperature, respectively. A power law of the type (13) is accurate only in the limit . To calculate critical exponents in the intermediate nonasymptotic regime, we must introduce additional correcting terms to power law (13). In accordance with the Wegner expansion (6) we have
(14) 
where are nonuniversal amplitudes and is the critical exponent of the correction to scaling. Here, in the calculation of the characteristics of disordered systems, we use the first correction to the asymptotic behavior for the correlation length and susceptibility:
(15) 
(16) 
and calculated critical exponents , and , as well as the critical temperatures using the least squares method for the best approximation of the data presented in Table 1 by expressions (15) and (16). Table 3 contains the values of critical parameters obtained for various spin concentrations using the initial data corresponding to various approximations for scaling functions, as well as their values averaged over approximations. It can be seen that the critical exponents form two groups with close values to within experimental error. The first group corresponds to and , i.e., to weakly disordered systems with spin concentrations , larger than the impurity percolation threshold ( for cubic systems), while the second group with and , corresponds to strongly disordered systems with , where is the spin percolation threshold ( for cubic systems), ; in the latter case, two mutually penetrating (spin and impurity) clusters exist in the system. Fractal effects of these two penetrating clusters may be responsible for the change in the type of critical behavior for strongly disordered systems. We can consider that the averaged values of critical exponents , and for weakly disordered systems and , and for strongly disordered systems are the final results of our investigations. It should be noted that the values of the critical exponents obtained for weakly disordered systems correlate with the values of , and (), obtained in (7) by the renormalizations group methods in the sixloop approximation, which are valid only for systems with low concentrations of impurities.
p=0.95  pol  0.6883(15)  1.3339(25)  0.141(52)  0.152(50)  4.26264(4)  4.26269(3) 

exp  0.6935(26)  1.3430(33)  0.113(64)  0.142(54)  4.26265(5)  4.26270(3)  
aver  0.6909(33)  1.3385(54)  0.137(56)  4.26267(4)  
p=0.80  pol  0.6960(29)  1.3473(30)  0.180(107)  0.193(74)  3.49937(21)  3.49954(14) 
exp  0.6947(28)  1.3421(30)  0.147(94)  0.192(71)  3.49940(21)  3.49961(14)  
aver  0.6956(29)  1.3447(40)  0.178(87)  3.49948(18)  
p=0.60  pol  0.7272(37)  1.4253(34)  0.221(147)  0.201(63)  2.42409(11)  2.42404(6) 
exp  0.7233(24)  1.4054(43)  0.184(92)  0.192(109)  2.42414(8)  2.42423(7)  
aver  0.7253(36)  1.4154(107)  0.199(103)  2.42413(9)  
p=0.50  pol  0.7372(25)  1.4299(26)  0.164(159)  0.195(74)  1.84503(7)  1.84512(3) 
exp  0.7368(26)  1.4266(30)  0.242(96)  0.226(66)  1.84503(7)  1.84519(3)  
aver  0.7370(33)  1.4283(33)  0.207(100)  1.84509(6) 
Authors  

Birgeneau et al.,  p=0.60  
1983, (8)  p=0.50  
Belanger et al.,  p=0.46  
1986, (9)  
Slanic et al.,  p=0.93  
1998, (10)  
Slanic et al.,  p=0.93  
1999, (11)  
Mitchell et al.,  p=0.75  
1986, (12)  p=0.50 
The above values of critical exponents and are also in good agreement with the available results of experiments with diluted Isinglike magnets (Table 3).
Authors  

Wang et al., 1990, (13)  300  
Heuer, 1993,(14)  60  
Wiseman et al.,1998, (15)  64  
80  
Ballesteros et al., 1998, (16)  128  
Calabrese et al., 2003, (17)  256  0.8  
Berche et al., 2005, (18)  96  0.85  
Murtazaev et al., 2004, (19)  60  
Table 4 contains the latest results of Monte Carlo simulation of the critical behavior of the diluted Ising model obtained by various authors. Each work cited in the table has its merits due to the application of various techniques for processing the results of simulation, as well as disadvantages associated with the small size of experimental lattices, which does not ensure reliable asymptotic values of the quantities being measured, or with insufficient statistics of averaging over various impurity configurations for obtaining reliable results, or with disregard of the effect of nonasymptotic corrections to scaling in the calculation of critical exponents (the inclusion of these corrections is especially important for samples with spin concentrations and and strongly disordered systems). The results obtained in (14); (19) can be treated as supporting our conclusions; these ideas were formulated in our earlier publications (20) on computer simulation of critical dynamics of the disordered Ising model. In fact, in spite of the attractiveness of the idea about a single universal critical behavior with the asymptotic values of critical exponents independent of the spin concentration, which was supported by the authors of (16), the results obtained in (16) did not make it possible to adequately explain the results obtained for samples with using the universal critical exponent of scaling correction for all systems. At the same time, the nonasymptotic values of critical exponents obtained in (16) demonstrated explicit dependence on and, assuming that is not unique, led to two sets of asymptotic critical exponents for weakly and strongly disordered systems. The results of the remaining studies carried out on samples with a single spin concentration coincide as a rule with our results to within experimental error, although some mismatching associated in all probability with the abovementioned drawbacks also takes place.
V Conclusions
The results of our investigations lead to the following conclusions.
 (i)

Scaling functions and values of critical exponents for the correlation length and susceptibility demonstrate the existence of two classes of universal critical behavior for the diluted Ising model with various characteristics for weakly and strongly disordered systems.
 (ii)

The values of critical exponents obtained for weakly disordered systems are in good agreement with the results of the fieldtheoretical description to within statistical errors of simulation and the numerical approximations used.
 (iii)

The results of simulation match the results of experimental studies of the critical behavior of diluted Isinglike magnets.
Acknowledgements.
The authors thank D.P. Landau, M. Novotny, V. Yanke, and N. Ito for fruitful discussions of the results of this study during the 3rd International Seminar on Computational Physics in Hangzhow (China, November 2006). This study was supported by the Russian Foundation for Basic Research (project nos. 040217524 and 040239000) and partly by the Program of the President of the Russian Federation (grant no. MK8738.2006.2). The authors are grateful to the Interdepartmental Supercomputer Center of the Russian Academy of Sciences for providing computational resources.Footnotes
 Email: prudnikv@univer.omsk.su
References
 A. B. Harris, J. Phys. C 7, 1671 (1974).
 R. Folk, Yu. Holovatch, and T. Yavorskii, Usp. Fiz. Nauk 173, 175 (2003) [Phys. Usp. 46, 169 (2003)].
 J. K. Kim, A. J. de Souza and D.P Landau, Phys. Rev. E 54, 2291 (1996).
 M. Hennecke and U. Heyken, J. Stat. Phys. 72, 829 (1993).
 D. Ivaneyko, J. Ilnytskyi, B. Berche et al., eprint condmat/0603521 (2006).
 F. J. Wegner, Phys. Rev. B 5, 4529 (1972).
 A. Pelissetto and E. Vicari, Phys. Rev. B 62, 6393 (2000).
 R. J. Birgeneau, R. A. Cowly, G. Shirane, et al., Phys. Rev. B 27, 6747 (1983).
 D. P. Belanger, A. R. King and V. Jaccarino, Phys. Rev. B 34, 452 (1986).
 D. P. Belanger, Z. Slanic and J. A. FernandezBaca, J. Magn. Magn. Mater 177181, 171 (1998).
 Z. Slanic, D. P. Belanger and J. A. FernandezBaca, Phys. Rev. Lett. 82, 426 (1999).
 P. W. Mitchell, R. A. Cowely, H. Yoshizawa et al., Phys. Rev. B 34, 4719 (1986).
 J. S. Wang, M. Wohlert, H. Muhlenbein et al., Physica A 166, 173 (1990).
 H. O. Heuer, J. Phys. A 26, L333 (1993).
 S. Wiseman and E. Domany, Phys. Rev. Lett. 81, 22 (1998); Phys. Rev. E 58, 2938 (1998).
 H.G. Ballesteros, L.A. Fernandez, V. MartinMayor et al., Phys. Rev. B 58, 2740 (1998).
 P. Calabrese, V. MartinMayor, A. Pelissetto et al., Phys. Rev. E 68, 036136 (2003).
 D. Ivaneyko, J. Ilnytskyi, B.Berche et al., Condens. Matter Phys. 8, 149 (2005).
 A. K. Murtazaev, I. K. Kamilov, and A. B. Babaev, Zh. Éksp. Teor. Fiz. 126, 1377 (2004) [JETP 99, 1201 (2004)].
 V. V. Prudnikov and A. N. Vakilov, Pis’ma Zh. Éksp. Teor. Fiz. 55, 709 (1992) [JETP Lett. 55, 741 (1992)]; Zh. Éksp. Teor. Fiz. 103, 962 (1993) [JETP 76, 469 (1993)].