- Email: [email protected]
Thermodynamic and critical properties of dilute XY magnets: Monte Carlo study
Solid State Communications 149 (2009) 1000–1003
Contents lists available at ScienceDirect
Solid State Communications journal homepage: www.elsevier.com/locate/ssc
Thermodynamic and critical properties of dilute XY magnets: Monte Carlo study Yun-Zhou Sun, Lin Yi ∗ , Yi-Hua Gao Department of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
article
info
Article history: Received 26 August 2008 Received in revised form 3 March 2009 Accepted 10 April 2009 by S. Miyashita Available online 21 April 2009 PACS: 07.05.Tp 21.60.Ka 02.70.Uu 75.10.Hk
abstract The thermodynamic and critical properties of two dimensional dilute XY spin magnets on a triangular lattice are studied by using a hybrid Monte Carlo method. The critical temperatures of different nonmagnetic impurity densities are obtained by two types of methods, including finite-size scaling of plane magnetic susceptibility and helicity modulus. It is found that the dilution plays an important effect on the critical temperature and thermodynamic equilibrium properties. The critical temperature decreases with the increase of magnetic impurity density ρ and the BKT phase transition vanishes at the site percolation threshold pc = 0.5. © 2009 Elsevier Ltd. All rights reserved.
Keywords: D. Thermodynamic properties D. BKT phase transition E. Mont Carlo simulation E. XY model
1. Introduction Over the years, the critical behavior and thermodynamic properties of many ideal statistical lattice spin models, such as planar rotator model, XY model and easy-plane Heisenberg model, have been deeply and extensively studied to describe a variety of physical situations, for example, solid solutions, magnetical ordered systems and adsorption phenomena. In some real physical systems, the occurrence of defects or impurities, if these are not magnetic, may affect the critical properties due to the missing bonds between spins. Recently, the study of nonmagnetic impurity in two dimensional spin models has been the subject of much interest [1–5]. It is well known that the so called Berezinskii–Kosterlitz–Thouless (BKT) phase transition [6,7] governed by the condensation of topological defects is found to be common in these spin models, Monte Carlo (MC) simulation showed that the BKT critical temperature decreases with the increase of non-magnetic impurity density ρ . However, MC simulation for planar rotator model with a Metropolis single spin update [1] showed that the critical temperature falls to zero at the non-magnetic impurity density ρ above the site percolation threshold on square lattice. This phenomena is also found in the dilute system of Josephson junctions [8]. On the other hand,
through extensive MC simulation, the critical temperature was found to drop to zero at the site percolation threshold [2,4]. The previous works were all focused on the discussion of spin system on the square lattice where the site percolation threshold is pc ≈ 0.59. As we all known, BKT phase transition also exists in the XY model on triangular lattice, where the site percolation threshold is pc = 0.5 [9]. On a triangular lattice, Wolff cluster algorithm has been used successfully to study the relation between percolation temperature, defined as the temperature at which spanning clusters start to appear in a large system, and the BKT critical temperature for planar rotator model [10]. On the other hand, BKT phase transition of the two dimensional XY magnets has very interesting features that make the question of the influence of disorder deduced by dilute magnets worth studying. Therefore, as a comparison, it is also interesting to study the dilute XY magnets on a triangular lattice by extensive MC simulation, especially in the low temperature and high non-magnetic density that the single spin Metropolis update is not efficient. In this work, we use an efficient hybrid MC method to discuss the effects of dilute spin magnets on the thermodynamic and critical properties. 2. Methods and results
∗
Corresponding author. Tel.: +86 27 87557705; fax: +86 27 87557705. E-mail address: [email protected] (L. Yi).
0038-1098/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ssc.2009.04.010
The Hamiltonian of the dilute XY magnets is described by [4]
Y.-Z. Sun et al. / Solid State Communications 149 (2009) 1000–1003
H = −J
X
σi σj sin θi sin θj cos(ϕi − ϕj ),
1001
(1)
hi , j i y
where spin vector SEi = (Six , Si , Siz ) = (sin θi cos ϕi , sin θi sin ϕi , cos θi ). J is the ferromagnetic coupling constant, and hi, ji indicates the nearest neighbor sites. σ is taken to be 1 or 0 depending on whether the spin magnet is occupied or not. Cluster algorithms, such as Swendsen-Wang, Wolff algorithm, are proved to be efficient to decrease critical slowing down, even in the present of anisotropic Dzyaloshinskii-Moriya interaction [11]. Here a hybrid MC approach, including one Wolff update [12] of planar components of the spins followed by one standard Metropolis update [13] and four over-relaxation updates [14,15], which change the configuration but keep the energy unchanged, was applied to this model on a triangular lattice. In an over-relaxation update, a spin is reflected across the effective potential of its P Ej = J j σj (SEjx + SEjy + SEjz ) with nearest neighbors, defined as B the reflection effected by SEi → 2
Ej SEi ·B
kBEj k2
Ej − SEi . The Wolff cluster B
update is similar to the application to a pure XY system, except that spins at vacant sites are all set to zero. In the actual simulation, a random bond can be constructed by choosing a random direction, projecting spins on that direction, and then assigning unit to spins of parallel or antiparallel projections to it. The over-relaxation and Wolff updates are efficient to reduce critical slowing down especially at low temperature where the single spin Metropolis update is inefficient. From an initial spin configuration by random selecting occupied site with probability 1 − ρ , the simulation was performed with periodic boundary conditions for system size N = L × L, where the sizes of lattice were considered as L = 20, 30, 40, 60, 80 and 100. Then the size of occupied site is Nmag = N (1 − ρ). All the measurements were carried out by increasing the temperature from a selected initial temperature. 1 ∼ 2 × 104 MC steps were discarded for thermal equilibration and about 4 ∼ 5×105 MC steps were used to get thermal averages at each temperature. In order to avoid a correlation, measurements were taken every 2 ∼ 6 MC steps. Two different methods, including finite-size scaling of inplane susceptibility and helicity modulus, were used to get the critical temperature TC for different ρ . Some thermodynamics are also obtained during the simulation. The specific heat is defined as CV = [hE 2 i − hE i2 ]/(Nmag kB T 2 ),
(2)
where E is the energy, kB is the Boltzmann constant. According P E = (Mx , My , Mz ) = i σi SEi , to the definition of magnetization M
E . The susceptibility and we define the order parameter m = M
susceptibility component can be written as [16]
χ = [hm2 i − hmi2 ]/Nmag kB T , χ α = Mα2 − hMα i2 /Nmag kB T .
(3) (4)
With the average of χ x and χ y , we can get the in-plane susceptibility
χ 0 = (χ x + χ y )/2.
(5)
For the dilute XY magnets, defined by Eq. (1), the expression of helicity modulus, Υ , on the triangular lattice can be defined as [17–19]
Υ (T ) = − √
hH i 3Nmag
−√
2J 2
×
#2 + X (ˆeij · xˆ )σi σj sin θi sin θj sin(ϕi − ϕj ) hi,ji
Here eˆ ij is the unite vector pointing from site j to site i. xˆ is a selected basis vector in one coordinate. In the presenting figures, the error bars that are not shown indicate the statistical errors are smaller than the symbols. For convenience, the temperature, the specific heat and the susceptibility are taken as the unites of J, J /kB and (g µB )2 /J respectively, where g is the Lande factor. Fig. 1 shows the result of the specific heat as a function of temperature at ρ = 0.05 for different lattice sizes. As shown in the figure, the maximum of specific heat is close to the low temperature as the increase of lattice sizes until stopping at the same temperature for large lattice sizes, such as L = 60, 80 and 100 in our simulation. The situation that the specific heat is independent of lattice sizes for large enough lattice is a typical indication of BKT transition character [20]. Then the maximum of specific heat is dislocated from BKT transition and the phase transition temperature can not be obtained by the information of specific heat. However, from the data of susceptibility, we can get some more useful information about critical properties. Here we take impurity density ρ = 0.05 as an example. Fig. 2 is a plot of susceptibility as a function of temperature. The location of the maximum is quickly near to the critical temperature TC with the increasing lattice sizes. From the maximum value of χ of different system sizes, according to the finite-size scaling relation χ ∝ L2−η , here we assert the relation still valid in dilute case, and then we can get the value of critical exponent η [21]. If the BKT phase transition exists, renormalization group theory shows that η equals 0.25 at the temperature of phase transition point. Our result of linear fit at ρ = 0.05 is η = 0.233 ± 0.005 which is comparable with the theoretical value 0.25. Then the phase transition in the case of dilute spin is still a BKT phase transition type, the same as that in pure model. Since the exisitence of BKT phase transition is verified, the calculation of critical temperature will be necessary. From the pseudo-critical temperature TC (L), obtained from the temperature of maximum of χ , the critical temperature, TC , can be obtained by the following finite-size relation [22,23] TC (L) = TC + π 2 /(4c (ln L)2 ).
2 3kB TNmag
*"
Fig. 1. Specific heat as a function of temperature at ρ = 0.05 for different lattice sizes.
.
(6)
(7)
The estimated critical temperature is TC = 1.023 ± 0.004 at ρ = 0.05 by this method. Though the application of this method to get TC is very efficient, but it is difficult to get TC (L) for the fluctuations at peak, as well as more MC steps and more times of simulation
1002
Y.-Z. Sun et al. / Solid State Communications 149 (2009) 1000–1003
Fig. 2. Susceptibility as a function of temperature at ρ = 0.05 for different lattice sizes.
Fig. 3. Application of the finite-size scaling of in-plane susceptibility to estimate critical temperature at ρ = 0.05.
are needed to get statistical average, especially at the high nonmagnetic density and low temperature. Therefore, we use another two efficient methods to get TC in the following discussions. As mentioned above, near and below TC , the finite-size scaling relation is still valid for in-plane susceptibility: χ 0 ∝ L2−η . Asserting the value of η is 0.25 at the critical point, then from the crossing point of χ 0 ∝ L2−η vs. T for different system sizes, we can get TC . This method is proved very efficient in many References [4,24,25]. As an example, Fig. 3 shows the results of this method to get TC at ρ = 0.05. The estimated TC is 1.014 ± 0.004 from the inset of this figure. This result is a little higher than the value obtained from TC (L). This is reasonable because here η is taken as a value in theory. Another method to estimate TC is based on the calculation of helicity modulus. Renormalization group theory [7] shows that there is a universal relation between the helicity modulus and the critical temperature. TC can be estimated from the intersection between Υ (T ) and the straight line Υ = 2kB T /π . For large enough system, the intersection will be nearer to the critical temperature, and Υ (T ) will have a deeper drop in the critical region. Due to the limit of system sizes in our simulation, an over-estimate of TC will be obtained by this method. As an example, Fig. 4 shows
Fig. 4. Helicity modulus as a function of temperature for several lattice sizes at ρ = 0.05. The inset shows the view near the estimated critical temperature.
the results of different lattice size at ρ = 0.05. From the largest lattice size L = 100, the critical temperature is estimated at TC = 1.030 ± 0.002, a little higher than the data from the inplane susceptibility χ 0 , but comparable with the result from TC (L). It is noted that the fluctuations or error bars are obvious for small system due to the effects of dilute magnets. For small systems, due to the initial distributions of spins are selected randomly, the effects of missing bonds between spins, which amplify the finite size effects, on the thermodynamics are more obvious than that for large systems, just as shown in the figures above. While for large enough system, the fluctuations due to impurity disorder will be smoothed out for low non-magnetic density. Then in our simulations, it is not necessary to take an average when ρ ≤ 0.25. However, in our practical application for this model, the statistical errors are obvious with the increase of ρ . In this work, the data of ρ > 0.25 are the results by statistical averages over four groups obtained by different random numbers. The helicity modulus for different ρ at L = 100 are shown in Fig. 5(a) and (b). Due to the critical slowing down in the low temperature, the statistical errors are obvious when ρ is near the threshold limit, where there is no spanning cluster to appear. The helicity modulus almost disappears at ρ = 0.50, just as shown in Fig. 5. Fig. 6 shows the last results of critical temperatures from two methods for different ρ . We find that TC decreases with the increase of ρ . On the other hand, from the plot, TC has an intuitionistic linear relation with ρ at both sides of ρ = 0.45. By linear fit of TC from Υ and χ 0 below ρ = 0.45, we find that the critical temperature falls to zero at ρ = 0.505 and ρ = 0.492 separately, almost the same as the site percolation threshold pc = 0.5. It is verified that the BKT phase transition vanishes at the triangular lattice percolation limit, the same result as this model on square lattice. 3. Conclusions In conclusion, the BKT critical temperatures of dilute XY magnets on a triangular lattice were obtained separately by finitesize scaling of in-plane susceptibility and renormalization group relation of helicity modulus by a hybrid MC algorithm. We found that the critical temperature decreases with the increase of ρ , and falls to zero at the site percolation limit. Meanwhile, by a linear fit of TC when ρ is near the percolation limit, we found that critical temperature vanishes at the site percolation threshold
Y.-Z. Sun et al. / Solid State Communications 149 (2009) 1000–1003
a
1003
b
Fig. 5. Helicity modulus as a function of temperature for lattice size L = 100 at different non-magnetic impurity density ρ . The overall trend and the results near percolation threshold are shown in (a) and (b) separately.
Acknowledgements This work was supported by the National Natural Science Foundation of China under Grants No. 10774053 and the Natural Science Foundation of Hubei Province under Grants No. 2007ABA035. References
Fig. 6. Critical temperature estimated by the two methods for different ρ . The straight line is the result of linear fit from helicity modulus near percolation limit.
pc = 0.5 and then the phase transition disappears. The BKT phase transition is also verified by peak of specific heat and critical exponents η, as well as the trend of helicity modulus. On appearance of dilute magnets, the lattice is less and less connected and the system becomes disorder. Though the disorder deduced by impurity becomes lager and larger when impurity density is near the percolation limit, the percolation spanning cluster still exists and the disorder is not enough to prevent the long range ordering (BKT) until impurity density exceeds the site percolation limit.
[1] S.A. Leonel, P.Z. Coura, A.R. Pereira, L.A.S. Mól, B.B. Costa, Phys. Rev. B 71 (2003) 094423. [2] B. Berche, A.I. Farińas-Sánchez, Yu. Holovatch, R. Paredes V, Eur. Phys. J. B 36 (2003) 91. [3] G.M. Wysin, Phys. Rev. B 71 (2005) 094423. [4] G.M. Wyain, A.R. Pereira, I.A. Marques, S.A. Leonel, P.Z. Coura, Phys. Rev. B 72 (2005) 094418. [5] L.M. Castro, A.S.T. Pires, J.A. Plascak, J. Magn. Magn. Mater. 248 (2002) 62. [6] V.L. Berezinskii, Sov. Phys. JEPT 34 (1972) 610. [7] J.M. Kosterlitz, D.J. Thouless, J. Phys. C 6 (1973) 1181. [8] Y.E. Lozovick, L.M. Pomirchi, Phys. Solid State 35 (1993) 1248. [9] D. Stauffer, A. Aharony, Introduction to Percolation Theory, Taylor and Francis, London, 1994. [10] P.W. Leung, C.L. Henley, Phys. Rev. B 43 (1991) 752. [11] Y.Z. Sun, L. Yi, X.G. Zhao, H.P. Liu, Solid State Commun. 144 (2007) 61. [12] U. Wolff, Phys. Rev. Lett. 62 (1989) 361. [13] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [14] F.R. Brown, T.J. Woch, Phys. Rev. Lett. 58 (1987) 2394. [15] M. Creutz, Phys. Rev. D 36 (1987) 515. [16] E. Rastelli, S. Regina, A. Tassi, Phys. Rev. B 70 (2004) 174447. [17] D.H. Lee, J.D. Joannopoulos, J.W. Negele, D.P. Landau, Phys. Rev. B 33 (1986) 450. [18] T. Ohta, D. Jasnow, Phys. Rev. B 20 (1979) 139. [19] S. Miyashita, H. Shiba, J. Phys. Soc. Jpn. 53 (1984) 1145. [20] B.V. Cosa, A.S. Pires, Phys. Rev. B 64 (2001) 092407. [21] S. Lee, K.C. Lee, Phys. Rev. B 57 (1998) 8472. [22] S.T. Bramwell, P.C.W. Holdsworth, Phys. Rev. B 49 (1994) 8811. [23] S.G. Chung, Phys. Rev. B 60 (1999) 11761. [24] L.A.S. Mól, A.R. Pereira, W.A. Moura-Melo, Phys. Rev. B 67 (2003) 132403. [25] L.A.S. Mól, A.R. Pereira, H. Chamati, S. Romano, Eur. Phys. J. B 50 (2006) 541.
Contents lists available at ScienceDirect
Solid State Communications journal homepage: www.elsevier.com/locate/ssc
Thermodynamic and critical properties of dilute XY magnets: Monte Carlo study Yun-Zhou Sun, Lin Yi ∗ , Yi-Hua Gao Department of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
article
info
Article history: Received 26 August 2008 Received in revised form 3 March 2009 Accepted 10 April 2009 by S. Miyashita Available online 21 April 2009 PACS: 07.05.Tp 21.60.Ka 02.70.Uu 75.10.Hk
abstract The thermodynamic and critical properties of two dimensional dilute XY spin magnets on a triangular lattice are studied by using a hybrid Monte Carlo method. The critical temperatures of different nonmagnetic impurity densities are obtained by two types of methods, including finite-size scaling of plane magnetic susceptibility and helicity modulus. It is found that the dilution plays an important effect on the critical temperature and thermodynamic equilibrium properties. The critical temperature decreases with the increase of magnetic impurity density ρ and the BKT phase transition vanishes at the site percolation threshold pc = 0.5. © 2009 Elsevier Ltd. All rights reserved.
Keywords: D. Thermodynamic properties D. BKT phase transition E. Mont Carlo simulation E. XY model
1. Introduction Over the years, the critical behavior and thermodynamic properties of many ideal statistical lattice spin models, such as planar rotator model, XY model and easy-plane Heisenberg model, have been deeply and extensively studied to describe a variety of physical situations, for example, solid solutions, magnetical ordered systems and adsorption phenomena. In some real physical systems, the occurrence of defects or impurities, if these are not magnetic, may affect the critical properties due to the missing bonds between spins. Recently, the study of nonmagnetic impurity in two dimensional spin models has been the subject of much interest [1–5]. It is well known that the so called Berezinskii–Kosterlitz–Thouless (BKT) phase transition [6,7] governed by the condensation of topological defects is found to be common in these spin models, Monte Carlo (MC) simulation showed that the BKT critical temperature decreases with the increase of non-magnetic impurity density ρ . However, MC simulation for planar rotator model with a Metropolis single spin update [1] showed that the critical temperature falls to zero at the non-magnetic impurity density ρ above the site percolation threshold on square lattice. This phenomena is also found in the dilute system of Josephson junctions [8]. On the other hand,
through extensive MC simulation, the critical temperature was found to drop to zero at the site percolation threshold [2,4]. The previous works were all focused on the discussion of spin system on the square lattice where the site percolation threshold is pc ≈ 0.59. As we all known, BKT phase transition also exists in the XY model on triangular lattice, where the site percolation threshold is pc = 0.5 [9]. On a triangular lattice, Wolff cluster algorithm has been used successfully to study the relation between percolation temperature, defined as the temperature at which spanning clusters start to appear in a large system, and the BKT critical temperature for planar rotator model [10]. On the other hand, BKT phase transition of the two dimensional XY magnets has very interesting features that make the question of the influence of disorder deduced by dilute magnets worth studying. Therefore, as a comparison, it is also interesting to study the dilute XY magnets on a triangular lattice by extensive MC simulation, especially in the low temperature and high non-magnetic density that the single spin Metropolis update is not efficient. In this work, we use an efficient hybrid MC method to discuss the effects of dilute spin magnets on the thermodynamic and critical properties. 2. Methods and results
∗
Corresponding author. Tel.: +86 27 87557705; fax: +86 27 87557705. E-mail address: [email protected] (L. Yi).
0038-1098/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ssc.2009.04.010
The Hamiltonian of the dilute XY magnets is described by [4]
Y.-Z. Sun et al. / Solid State Communications 149 (2009) 1000–1003
H = −J
X
σi σj sin θi sin θj cos(ϕi − ϕj ),
1001
(1)
hi , j i y
where spin vector SEi = (Six , Si , Siz ) = (sin θi cos ϕi , sin θi sin ϕi , cos θi ). J is the ferromagnetic coupling constant, and hi, ji indicates the nearest neighbor sites. σ is taken to be 1 or 0 depending on whether the spin magnet is occupied or not. Cluster algorithms, such as Swendsen-Wang, Wolff algorithm, are proved to be efficient to decrease critical slowing down, even in the present of anisotropic Dzyaloshinskii-Moriya interaction [11]. Here a hybrid MC approach, including one Wolff update [12] of planar components of the spins followed by one standard Metropolis update [13] and four over-relaxation updates [14,15], which change the configuration but keep the energy unchanged, was applied to this model on a triangular lattice. In an over-relaxation update, a spin is reflected across the effective potential of its P Ej = J j σj (SEjx + SEjy + SEjz ) with nearest neighbors, defined as B the reflection effected by SEi → 2
Ej SEi ·B
kBEj k2
Ej − SEi . The Wolff cluster B
update is similar to the application to a pure XY system, except that spins at vacant sites are all set to zero. In the actual simulation, a random bond can be constructed by choosing a random direction, projecting spins on that direction, and then assigning unit to spins of parallel or antiparallel projections to it. The over-relaxation and Wolff updates are efficient to reduce critical slowing down especially at low temperature where the single spin Metropolis update is inefficient. From an initial spin configuration by random selecting occupied site with probability 1 − ρ , the simulation was performed with periodic boundary conditions for system size N = L × L, where the sizes of lattice were considered as L = 20, 30, 40, 60, 80 and 100. Then the size of occupied site is Nmag = N (1 − ρ). All the measurements were carried out by increasing the temperature from a selected initial temperature. 1 ∼ 2 × 104 MC steps were discarded for thermal equilibration and about 4 ∼ 5×105 MC steps were used to get thermal averages at each temperature. In order to avoid a correlation, measurements were taken every 2 ∼ 6 MC steps. Two different methods, including finite-size scaling of inplane susceptibility and helicity modulus, were used to get the critical temperature TC for different ρ . Some thermodynamics are also obtained during the simulation. The specific heat is defined as CV = [hE 2 i − hE i2 ]/(Nmag kB T 2 ),
(2)
where E is the energy, kB is the Boltzmann constant. According P E = (Mx , My , Mz ) = i σi SEi , to the definition of magnetization M
E . The susceptibility and we define the order parameter m = M
susceptibility component can be written as [16]
χ = [hm2 i − hmi2 ]/Nmag kB T , χ α = Mα2 − hMα i2 /Nmag kB T .
(3) (4)
With the average of χ x and χ y , we can get the in-plane susceptibility
χ 0 = (χ x + χ y )/2.
(5)
For the dilute XY magnets, defined by Eq. (1), the expression of helicity modulus, Υ , on the triangular lattice can be defined as [17–19]
Υ (T ) = − √
hH i 3Nmag
−√
2J 2
×
#2 + X (ˆeij · xˆ )σi σj sin θi sin θj sin(ϕi − ϕj ) hi,ji
Here eˆ ij is the unite vector pointing from site j to site i. xˆ is a selected basis vector in one coordinate. In the presenting figures, the error bars that are not shown indicate the statistical errors are smaller than the symbols. For convenience, the temperature, the specific heat and the susceptibility are taken as the unites of J, J /kB and (g µB )2 /J respectively, where g is the Lande factor. Fig. 1 shows the result of the specific heat as a function of temperature at ρ = 0.05 for different lattice sizes. As shown in the figure, the maximum of specific heat is close to the low temperature as the increase of lattice sizes until stopping at the same temperature for large lattice sizes, such as L = 60, 80 and 100 in our simulation. The situation that the specific heat is independent of lattice sizes for large enough lattice is a typical indication of BKT transition character [20]. Then the maximum of specific heat is dislocated from BKT transition and the phase transition temperature can not be obtained by the information of specific heat. However, from the data of susceptibility, we can get some more useful information about critical properties. Here we take impurity density ρ = 0.05 as an example. Fig. 2 is a plot of susceptibility as a function of temperature. The location of the maximum is quickly near to the critical temperature TC with the increasing lattice sizes. From the maximum value of χ of different system sizes, according to the finite-size scaling relation χ ∝ L2−η , here we assert the relation still valid in dilute case, and then we can get the value of critical exponent η [21]. If the BKT phase transition exists, renormalization group theory shows that η equals 0.25 at the temperature of phase transition point. Our result of linear fit at ρ = 0.05 is η = 0.233 ± 0.005 which is comparable with the theoretical value 0.25. Then the phase transition in the case of dilute spin is still a BKT phase transition type, the same as that in pure model. Since the exisitence of BKT phase transition is verified, the calculation of critical temperature will be necessary. From the pseudo-critical temperature TC (L), obtained from the temperature of maximum of χ , the critical temperature, TC , can be obtained by the following finite-size relation [22,23] TC (L) = TC + π 2 /(4c (ln L)2 ).
2 3kB TNmag
*"
Fig. 1. Specific heat as a function of temperature at ρ = 0.05 for different lattice sizes.
.
(6)
(7)
The estimated critical temperature is TC = 1.023 ± 0.004 at ρ = 0.05 by this method. Though the application of this method to get TC is very efficient, but it is difficult to get TC (L) for the fluctuations at peak, as well as more MC steps and more times of simulation
1002
Y.-Z. Sun et al. / Solid State Communications 149 (2009) 1000–1003
Fig. 2. Susceptibility as a function of temperature at ρ = 0.05 for different lattice sizes.
Fig. 3. Application of the finite-size scaling of in-plane susceptibility to estimate critical temperature at ρ = 0.05.
are needed to get statistical average, especially at the high nonmagnetic density and low temperature. Therefore, we use another two efficient methods to get TC in the following discussions. As mentioned above, near and below TC , the finite-size scaling relation is still valid for in-plane susceptibility: χ 0 ∝ L2−η . Asserting the value of η is 0.25 at the critical point, then from the crossing point of χ 0 ∝ L2−η vs. T for different system sizes, we can get TC . This method is proved very efficient in many References [4,24,25]. As an example, Fig. 3 shows the results of this method to get TC at ρ = 0.05. The estimated TC is 1.014 ± 0.004 from the inset of this figure. This result is a little higher than the value obtained from TC (L). This is reasonable because here η is taken as a value in theory. Another method to estimate TC is based on the calculation of helicity modulus. Renormalization group theory [7] shows that there is a universal relation between the helicity modulus and the critical temperature. TC can be estimated from the intersection between Υ (T ) and the straight line Υ = 2kB T /π . For large enough system, the intersection will be nearer to the critical temperature, and Υ (T ) will have a deeper drop in the critical region. Due to the limit of system sizes in our simulation, an over-estimate of TC will be obtained by this method. As an example, Fig. 4 shows
Fig. 4. Helicity modulus as a function of temperature for several lattice sizes at ρ = 0.05. The inset shows the view near the estimated critical temperature.
the results of different lattice size at ρ = 0.05. From the largest lattice size L = 100, the critical temperature is estimated at TC = 1.030 ± 0.002, a little higher than the data from the inplane susceptibility χ 0 , but comparable with the result from TC (L). It is noted that the fluctuations or error bars are obvious for small system due to the effects of dilute magnets. For small systems, due to the initial distributions of spins are selected randomly, the effects of missing bonds between spins, which amplify the finite size effects, on the thermodynamics are more obvious than that for large systems, just as shown in the figures above. While for large enough system, the fluctuations due to impurity disorder will be smoothed out for low non-magnetic density. Then in our simulations, it is not necessary to take an average when ρ ≤ 0.25. However, in our practical application for this model, the statistical errors are obvious with the increase of ρ . In this work, the data of ρ > 0.25 are the results by statistical averages over four groups obtained by different random numbers. The helicity modulus for different ρ at L = 100 are shown in Fig. 5(a) and (b). Due to the critical slowing down in the low temperature, the statistical errors are obvious when ρ is near the threshold limit, where there is no spanning cluster to appear. The helicity modulus almost disappears at ρ = 0.50, just as shown in Fig. 5. Fig. 6 shows the last results of critical temperatures from two methods for different ρ . We find that TC decreases with the increase of ρ . On the other hand, from the plot, TC has an intuitionistic linear relation with ρ at both sides of ρ = 0.45. By linear fit of TC from Υ and χ 0 below ρ = 0.45, we find that the critical temperature falls to zero at ρ = 0.505 and ρ = 0.492 separately, almost the same as the site percolation threshold pc = 0.5. It is verified that the BKT phase transition vanishes at the triangular lattice percolation limit, the same result as this model on square lattice. 3. Conclusions In conclusion, the BKT critical temperatures of dilute XY magnets on a triangular lattice were obtained separately by finitesize scaling of in-plane susceptibility and renormalization group relation of helicity modulus by a hybrid MC algorithm. We found that the critical temperature decreases with the increase of ρ , and falls to zero at the site percolation limit. Meanwhile, by a linear fit of TC when ρ is near the percolation limit, we found that critical temperature vanishes at the site percolation threshold
Y.-Z. Sun et al. / Solid State Communications 149 (2009) 1000–1003
a
1003
b
Fig. 5. Helicity modulus as a function of temperature for lattice size L = 100 at different non-magnetic impurity density ρ . The overall trend and the results near percolation threshold are shown in (a) and (b) separately.
Acknowledgements This work was supported by the National Natural Science Foundation of China under Grants No. 10774053 and the Natural Science Foundation of Hubei Province under Grants No. 2007ABA035. References
Fig. 6. Critical temperature estimated by the two methods for different ρ . The straight line is the result of linear fit from helicity modulus near percolation limit.
pc = 0.5 and then the phase transition disappears. The BKT phase transition is also verified by peak of specific heat and critical exponents η, as well as the trend of helicity modulus. On appearance of dilute magnets, the lattice is less and less connected and the system becomes disorder. Though the disorder deduced by impurity becomes lager and larger when impurity density is near the percolation limit, the percolation spanning cluster still exists and the disorder is not enough to prevent the long range ordering (BKT) until impurity density exceeds the site percolation limit.
[1] S.A. Leonel, P.Z. Coura, A.R. Pereira, L.A.S. Mól, B.B. Costa, Phys. Rev. B 71 (2003) 094423. [2] B. Berche, A.I. Farińas-Sánchez, Yu. Holovatch, R. Paredes V, Eur. Phys. J. B 36 (2003) 91. [3] G.M. Wysin, Phys. Rev. B 71 (2005) 094423. [4] G.M. Wyain, A.R. Pereira, I.A. Marques, S.A. Leonel, P.Z. Coura, Phys. Rev. B 72 (2005) 094418. [5] L.M. Castro, A.S.T. Pires, J.A. Plascak, J. Magn. Magn. Mater. 248 (2002) 62. [6] V.L. Berezinskii, Sov. Phys. JEPT 34 (1972) 610. [7] J.M. Kosterlitz, D.J. Thouless, J. Phys. C 6 (1973) 1181. [8] Y.E. Lozovick, L.M. Pomirchi, Phys. Solid State 35 (1993) 1248. [9] D. Stauffer, A. Aharony, Introduction to Percolation Theory, Taylor and Francis, London, 1994. [10] P.W. Leung, C.L. Henley, Phys. Rev. B 43 (1991) 752. [11] Y.Z. Sun, L. Yi, X.G. Zhao, H.P. Liu, Solid State Commun. 144 (2007) 61. [12] U. Wolff, Phys. Rev. Lett. 62 (1989) 361. [13] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [14] F.R. Brown, T.J. Woch, Phys. Rev. Lett. 58 (1987) 2394. [15] M. Creutz, Phys. Rev. D 36 (1987) 515. [16] E. Rastelli, S. Regina, A. Tassi, Phys. Rev. B 70 (2004) 174447. [17] D.H. Lee, J.D. Joannopoulos, J.W. Negele, D.P. Landau, Phys. Rev. B 33 (1986) 450. [18] T. Ohta, D. Jasnow, Phys. Rev. B 20 (1979) 139. [19] S. Miyashita, H. Shiba, J. Phys. Soc. Jpn. 53 (1984) 1145. [20] B.V. Cosa, A.S. Pires, Phys. Rev. B 64 (2001) 092407. [21] S. Lee, K.C. Lee, Phys. Rev. B 57 (1998) 8472. [22] S.T. Bramwell, P.C.W. Holdsworth, Phys. Rev. B 49 (1994) 8811. [23] S.G. Chung, Phys. Rev. B 60 (1999) 11761. [24] L.A.S. Mól, A.R. Pereira, W.A. Moura-Melo, Phys. Rev. B 67 (2003) 132403. [25] L.A.S. Mól, A.R. Pereira, H. Chamati, S. Romano, Eur. Phys. J. B 50 (2006) 541.