- Email: [email protected]
The bilayer Ising model and a generalized Husimi tree approximation
Available online at www.sciencedirect.com
Physica A 335 (2004) 563 – 576
www.elsevier.com/locate/physa
The bilayer Ising model and a generalized Husimi tree approximation James L. Monroe Department of Physics, Penn State University, Beaver Campus, 100 University Dr., Monaca, PA 15061, USA Received 10 September 2003
Abstract We study a variety of aspects of a bilayer system of Ising spins including accurate estimates of the critical temperature for ferromagnetic interactions, scaling of the critical temperature when the interlayer interaction goes to zero, and approximations of phase diagrams for the case when antiferromagnetic interlayer interactions are present including location of the tricritical point. This is all done using an extension of a method previously developed by the author to study a variety of lattice spin systems and which is a generalization of the Bethe lattice approach. c 2003 Elsevier B.V. All rights reserved. PACS: 05.50.+q; 64.60.Cn; 75.70.−i Keywords: Ising spins; Thin 8lms; Bilayer system
1. Introduction For nearly four decades models of the magnetic properties of thin 8lms have been of interest both in terms of various aspects of practical applications as well as fundamental theoretical issues. The simplest model of a thin 8lm, and the 8rst to be studied, is the bilayer system consisting of two layers of Ising spins with nearest neighbor, intraplane interactions, and interplane interactions. In 1964, Ballantine [1] using high-temperature series expansion methods obtained an estimate for the critical temperature of such a system where all interactions were ferromagnetic and of the same strength. Since that time there has been continual interest in obtaining better estimates of the critical temperature for such a system through lengthier series expansions [2,3], analytical results based on mean-8eld theory [4], the Bethe approximation [5,6] and a variational E-mail address: [email protected] (J.L. Monroe). c 2003 Elsevier B.V. All rights reserved. 0378-4371/$ - see front matter doi:10.1016/j.physa.2003.12.018
564
J.L. Monroe / Physica A 335 (2004) 563 – 576
cumulant expansion [7] as well as many other numerical methods including Monte Carlo [8], interfacial tensions approximation [9], a transfer-matrix version of mean-8eld [9,10], a constrained variational approach and correlation length equality [11,12], and a corner-transfer-matrix, renormalization group approach [13]. Additionally the behavior of the critical temperature of such a system in the limit of the interplane interaction going to zero and to a lesser extent the behavior when this interaction approaches in8nity has been studied [4,10,14–17]. Finally the general phase diagram has a rich structure especially when either or both the interlayer and intralayer interactions are allowed to be antiferromagnetic. In the following, we use a method already used to obtain results for frustrated systems such as the Ising antiferromagnet on the triangle lattice [18], Ising systems with multisite interactions [19,20], higher-spin Ising and Potts systems [21], and in its most developed form for the Ising antiferromagnet and ferromagnet on the square lattice [22]. Here it will be used in all three areas mentioned above. First, for the case of equal ferromagnetic intraplaner interactions where the phase diagram is simply a line of phase transitions in zero magnetic 8eld ending in at a critical endpoint (the emphasis in this case is to accurately determine the critical temperature). Second, to examine the behavior of the system in the limit of vanishing interplaner interactions both for the case where the intraplaner interactions on both layers are equal and when they are not as well as when not only the interlayer interaction goes to zero but one of the intraplaner interactions does so as well. Third, to approximate the phase diagram for the bilayer system having an antiferromagnetic interplaner interaction along with ferromagnetic intraplaner interactions when there is a magnetic 8eld present. In this case, we have an example of a metamagnet and there is present a tricritical point. Basically the method we will use is a generalization of the Bethe lattice approach. The Bethe lattice can be thought of as being built in a recursive manner by attaching to a central site on which there is an Ising spin, q Ising spins each interacting with the central site through a pair interaction. Then to each of these q spins are attached q − 1 spins each interacting with one of the q spins through a pair interaction. One continues in this manner ad in8nitum to construct the full Bethe lattice. Of course this will not work when one wants to approximate a bilayer system but one can easily imagine two identical Bethe lattices coupled together with an interplaner interaction as an approximation for the bilayer system. This is the system that has been looked at in some detail very recently by Hu et al. [6]. The basic building block in the construction of a Bethe lattice is a pair of nearest neighbor Ising spins and the interaction connecting them. If one wants to approximate a square lattice one can take q = 4. However, a better approximation of the square lattice can be obtained by taking as a basic building block a square on each corner of which is an Ising spin and each Ising spin interacting with its two nearest neighbors. One can take one of these as a central square and consider attaching a similar square to each corner of the central square. One can then take this collection of sites and to each Ising spin not already having four nearest neighbor interactions attach another of the 4-site square building blocks. Continuing this ad in8nitum one constructs what is denoted by graph theorists as an Husimi tree [23]. The behavior of a spin deep inside the tree can be used as an approximation for the spontaneous magnetization
J.L. Monroe / Physica A 335 (2004) 563 – 576
565
of the original square lattice system just as a site deep inside the Bethe lattice can be used as an approximation of the original square lattice system. Since the Husimi tree approximation more correctly takes into account local Fuctuations than the Bethe lattice approximation one 8nds it more correctly approximates the original system. Of course as described for the Bethe lattice approximation of the bilayer system where two Bethe lattices were connected with interplaner interactions the same can be done for the Husimi tree approximation. Furthermore, there is no need to consider only a 4-site square as the basic building block larger building blocks can be used including ones where rather than connections being made at a single site connections are made at a collection of sites. In Ref. [22] this was done using basic building blocks of as many as 60 sites to approximate the square lattice ferromagnetic Ising model as well as the antiferromagnetic case where a magnetic 8eld is present. In the ferromagnetic case the exact value of the critical temperature is known and results as close as 0.003% of the actual value were obtained. In the following section, the Hamiltonian of the bilayer system is presented as well as a more detailed description of the method used here to study the bilayer Ising model. Then Sections 3–5 contain the results of the three diHerent areas mentioned above. We close in Section 6 with some concluding remarks. 2. Basic method In the following, we approximate a variety of properties of a system consisting of two parallel, horizontal, square lattice Ising models, where on the upper layer the Ising spin on the ith site will be denoted by i and on the lower layer by i , and where by Ising spin we mean speci8cally a spin allowed to take on the values ±1. The Hamiltonian of the system is H = −JA i j − J B i j − IJA i i − h A i − h B i ; (1) ¡i; j¿
¡i; j¿
¡i; j¿
i
i
where ¡ i; j ¿ represents a nearest neighbor pair of sites, JA and JB are, respectively, the intralayer interactions on the upper and lower layers of the two-layer system, IJA is the interlayer interaction, and hA and hB are, respectively, the external magnetic 8eld present on the upper and lower layers of the system. As stated above a single layer of the system can be approximated by a Husimi tree approach where one takes as a basic building block a set of four sites, each site on the corner of a square, and on each site an Ising spin. One can start with a central square and attach to it at each corner a new square. Then taking this arrangement one can attach to each site not having four nearest neighbors a new square. One can proceed in this way, in the limit building an in8nite Husimi tree. This process because of the step-by-step nature of the construction can be described by a recursive system of equations, i.e., a dynamical system, which can be used to calculate, for example, the thermal average of a site of the central square. Based on the behavior of this thermal average one can determine a critical temperature. In terms of dynamical systems language one has for this system a one-dimensional dynamical system and for high
566
J.L. Monroe / Physica A 335 (2004) 563 – 576
temperatures the behavior of the system is governed by a single attracting 8xed point which corresponds to zero magnetization of any one of the sites of the central square. For the case where a continuous phase transition occurs as the temperature is lowered this 8xed point becomes repelling and two new attracting 8xed points are created, i.e., a bifurcation occurs. These two new 8xed points correspond to nonzero thermal averages of a site on the central square, one positive and negative, and hence a spontaneous magnetization. Therefore, the temperature at which the single 8xed point becomes repelling and two new attracting 8xed points are created is the critical temperature. See reference [22] for details. Now in analogy when trying to approximate the bilayer system one can imagine constructing a tree-like structure with a basic building block consisting of two parallel, horizontal squares of four sites each and interacting via nearest neighbor intralayer interactions JA and JB . In addition one has an interlayer, vertical interaction, between nearest neighbor sites in the vertical direction. Hence, one has a basic building block which is a cube consisting of eight sites. Now one can construct a system by considering a central cube and attaching to each vertical edge of this central cube a basic building block. One sees now that the vertical interaction in each building block, i.e., the interlayer interaction, should be taken as only one-half of IJA for the interaction strength. This is because when making the attachments the vertical interactions get a contribution from that of the original building block as well as the new building block the sum of which is the full IJA . Since each connection consists of two sites, one on the upper layer and one on the lower layer, the dynamical system which results is a three-dimensional dynamical system. The dimension of the dynamical system, disregarding any symmetry considerations, is always 2n − 1, where n is the number of sites involved in a connection, again see Ref. [22] for details. The mathematical behavior of the resulting three-dimensional dynamical system (here n = 2), i.e., the presence of attracting 8xed points, two-cycles, etc., governs the physical behavior of the system. For the case where JA , JB , and IJA are positive, i.e., ferromagnetic, and the external magnetic 8elds are zero then based on physical grounds one expects for high temperatures a single attracting 8xed point which bifurcates to two new attracting 8xed points at some point as the temperature is lowered as described previously. This is exactly what one 8nds occurring for this case. As with the case of the square lattice Ising model approximated in Ref. [22] one can consider larger building blocks. In fact in Ref. [22] a very systematic sequence of 8ve diHerent building blocks is used to construct various Husimi trees and the results based on them were presented in the approximation of the square lattice Ising model. For approximating, the bilayer system two of each of these basic building blocks can be placed in parallel, the upper one with nearest neighbor interactions JA and the lower one with interactions JB . Then the two can be coupled to each other with interactions IJA for sites not involved in a connection and one-half IJA for sites on the edge of the building block and hence involved in a connection. Of course the dimension of the dynamical system grows rapidly. To use the 8fth of the 8ve building blocks as done in the square lattice case for this bilayer case would mean 10 sites were involved in each connection resulting in a dynamical system of dimension 210 − 1, which is beyond the computer capabilities used here.
J.L. Monroe / Physica A 335 (2004) 563 – 576
567
All computations have been executed on a personal computer using Mathematica to actually generate the dynamical systems governing the various approximations. Because of this for the bilayer system being approximated here we have only been able to obtain three levels of approximation rather than the 8ve for the square lattice. Even so as will be seen in the next section we are able to obtain reasonably accurate estimates of the critical temperature where by “reasonably” we mean estimates having approximately the same accuracy as obtained by presently available series expansion methods. To obtain the best approximations for the critical behavior we have used our systematic sequence of three approximations along with the Bulirsch-Stoer extrapolation procedure [24–27]. Here one wishes to determine the value of Tc of the bilayer system by knowing the corresponding values of Tc based on some sequence of approximations using increasingly larger building blocks. The value of Tc using the approximations will depend on the size of the basic building blocks used and one denotes them as Tc (L) where L is a measure of the system size. One supposes in general that Tc (L) can be written as Tc (L) = Tc; ∞ + a1 L−!1 + a2 L−!2 + a3 L−!3 + · · · ;
(2)
where 0 ¡ !1 ¡ !2 ¡ !3 ¡ · · · and where Tc; ∞ is the value for the system being approximated. The BST algorithm allows one to start with a sequence of values, here three values Tc (L1 ), Tc (L2 ), and Tc (L3 ) and estimate Tc; ∞ . In particular the algorithm allows one to construct a table of extrapolants of this sequence T01 T11 T02
T21 T12
T03 where T21 is the algorithm’s best estimate for Tc; ∞ . The Tqn are computed from n =0 ; T−1
(3)
T0n = Tc (n) ;
(4)
Tmn
=
n+1 Tm−1
+
n+1 (Tm−1
−
n Tm−1 )
Ln Ln+m
−1 ! n+1 n Tm−1 − Tm−1 1 − n+1 −1 ; n+1 Tm−1 − Tm−2
(5)
where m ¿ 1 and where ! is a free parameter. 3. Estimates of the critical temperature The values of Tc (L) for the three generalized Husimi tree systems used here are given in Table 1 for a variety of combinations of JA , JB ; and IJA with the restriction that all three are ferromagnetic. For simplicity we set JA =1 and then vary the values of
568
J.L. Monroe / Physica A 335 (2004) 563 – 576
Table 1 Critical temperature estimates based on the Bethe lattice as well as the three levels of Husimi trees described in the text
JB JB JB JB JB JB JB JB JB JB JB JB JB
= 1, = 4 = 1, = 2 = 1, = 1 = 1, = 0:5 = 1, = 0:1 = 1, = 0:01 = 1, = 0:001 = 1, = 0:0001 = 1, = 0 = 1=2, = 0:1 = 1=2, = 0:01 = 1=2, = 0:001 = 1=2, = 0:0001
Bethe lattice, results
8-site results
24-site results
48-site results
4.86557 4.20205 3.68483 3.33961 2.98923 2.89616 2.88647 2.88550 2.88539
4.784395068 4.175811973 3.677368547 3.322723648 2.915324057 2.786786000 2.772402380 2.770944746 2.770782562 2.776106370 2.770836387 2.770783103 2.770782568
4.512328839 3.957293525 3.491491794 3.153858714 2.742398308 2.588782964 2.569183645 2.567151728 2.566925090
4.396933786 3.864416292 3.412610923 3.083346974 2.675481529 2.509447499 2.485612169 2.483066419 2.482781460
JB , and . The three systems used here are the generalizations of the 8rst three systems used in Ref. [22] to approximate the square lattice, Ising model to the bilayer system as described above. Hence, they consist of eight sites (the cube described in detail above), a 24-site system consisting of two parallel 12-site building blocks (see Ref. [12]) connected by interlayer interactions, and 8nally a 48-site system consisting of two parallel 24-site building blocks (see Ref. [12]) connected by interlayer interactions. In addition to the Husimi tree estimates in Table 1 we give the estimates found from the Bethe lattice approximation [5,6] which in some sense is a zeroth-order approximation if we consider our three estimates as 8rst, second, and third-order estimates. Both Refs. [5,6] give analytical expressions for the critical temperature although the expression derived in Ref. [6] requires JA = JB . Unfortunately the critical temperatures given by the two expressions diHer except at the two extremes of = 0 and → ∞. We have using the dynamical systems approach applied to our three levels of approximation found critical temperature values for the Bethe lattice which match those given by the analytical expressions of Hu et al. in Ref. [6] and therefore use these values as our Bethe lattice results listed in Table 1. We have presented our estimates of Tc to nine places past the decimal as this level of precision is necessary for the use of the BST algorithm which will be used to get our best estimates of the critical temperature. As expected one sees a systematic decrease in the estimated value of Tc as the size of the basic building blocks is increased. The decrease in the estimated value of Tc corresponds to an increase in the accuracy of the approximations as increasing the size of the basic building blocks allows the system to more properly take into account local Fuctuations. As is the case for systems with short range interactions such as the nearest neighbor interactions present in the system under study here the estimates of Tc found using these approaches are typically too high. Using the three estimates based on the Husimi tree approach given in Table I along with the BST extrapolation algorithm we arrive at our best estimates for Tc . These
J.L. Monroe / Physica A 335 (2004) 563 – 576
569
Table 2 Critical temperature estimates based on the BST algorithm and the three Husimi tree results presented in Table 1 along with estimates based on high temperature series expansions [3,4], constrained variational method and correlation length equality [12], transfer matrix, mean-8eld [10], the corner transfer matrix method [13], and Monte Carlo BST estimates JB JB JB JB JB JB JB JB JB JB
= 1, = 1, = 1, = 1, = 1, = 1, = 1, = 1, = 1, = 1,
I=8 =4 =2 =1 = 0:5 = 0:1 = 0:01 = 0:001 = 0:0001 =0
4.095478 3.621238 3.206424 2.902182 2.516101 2.314020 2.274273 2.269701 2.269185
High temp. series exp.
3:623 ± 0:015 3:208 ± 0:007 2:911 ± 0:006 2:509 ± 0:010 2:307 ± :011 2:271 ± 0:005
CVM & CLE 4.4458 4.0968 3.6233 3.2085 2.5150 2.3350
Trans. mat. & MF
3.6227 3.207530 2.901246 2.529991
Corner trans. mat.
3.20762 2.51429 2.33322
Monte Carlo 4.098 3.610 3.226 2.915 2.532 2.299
2.26917
are given in Table 2 along with estimates based on several other approaches. In the use of the BST algorithm one must choose a value for ! as mentioned following Eq. (5). This is a critical issue. In our case we can use the fact that when = 0 one has Onsager’s exact result and using our three estimates for Tc when = 0 we vary ! determining the value which gives to the numerical accuracy used here Onsager’s result. This occurs when ! = 0:940486. We use this value of ! for the cases where = 0. Also, L is a measure of the system size but the proper value to use here is not entirely clear. Two possibilities are readily apparent one is to take L as the square root of the number of sites in a basic building block and the other is to use the number of sites used in each of the connections made in the construction of the tree, i.e., L = 2, 4, or 6. We have chosen to use the later. The most studied case is the one where all three interactions JA , JB ; and IJA are of equal strength. Lipowski and Suzuki [9] 8nd using the transfer matrix, mean-8eld approach that Tc 3:20755(5). The method is similar to the one presented here in two major ways 8rst it is similar in that one can look at larger and larger systems obtaining increasingly accurate results and second one thereby develops a sequence of approximations which can be extrapolated to get a 8nal best estimate. Lipowski [10] extends the earlier work to obtain Tc 3:207530. These results are best supported by the results of Li et al. [13], which diHer from that of Lipowski by less than 0:003%. Our results underestimate the critical temperature by approximately as much as the CVM and CLE results of Ref. [12] over estimate it, with a diHerence of 0:03% with the Lipowski result. The Monte Carlo results from Ref. [17] give a critical temperature estimate Tc 3:226 which is signi8cantly above our estimate as well as the other four estimates presented in Table 2. As one decreases the value of the interlayer interaction the situation in terms of accuracy of the estimates is not at all clear. For example for = 0:1 from Table 2 we see that the results of Ref. [10] are higher than those of all other method of
570
J.L. Monroe / Physica A 335 (2004) 563 – 576
approximation except the Monte Carlo results of Ref. [17] and are in fact outside the error bounds presented with the high temperature series expansion results. The results of Refs. [12,13], and our own give for = 0:1 a best result of Tc 2:515(1). With further reduction of the value to = 0:01 one sees the corner matrix results or the CVM and CLE results are now greater than that obtained by either the estimates of the Husimi tree, BST extrapolation method or series expansions method with the former estimates even outside the error bars of the series expansion results. Unfortunately numerical values for Tc with still further reduction of have not been given in other references. For values of above 1.0 only for = 2 are numerical values for the estimate of Tc generally available. Five results are presented in Table 2 and there is general agreement amongst four of the results with now the Monte Carlo results being signi8cantly diHerent on the low side. In general there is little systematics comparing one method to another so good estimates of the errors are diPcult. 4. Tc as → 0 The behavior of Tc as → 0 has been of interest since 1970 when Abe [14] in terms of a scaling theory derived a value for the shift exponent which describes the deviation of the critical temperature Tc () from the critical temperature in the decoupled limit, i.e., Tc (0), with de8ned through the following: Tc () − Tc (0) ∼ :
(6)
Suzuki [15] followed up Abe’s work with further analysis. Both predict that when the intralayer interactions, JA and JB ; are equal that the shift exponent is equal to 1=, being the critical exponent of the susceptibility of the two-dimensional model. Oitmaa and Enting [4] generalized to the situation where the intralayer interactions are not equal and found for this situation that = 2=. As Lipowski [10] emphasizes this is particularly interesting because we now are coupling two systems which have diHerent critical temperatures. Lipowski also points out that as of 1998, the time of his paper, these scaling predictions had only modest con8rmation. In fact for the case of JA = JB contradictory predictions existed with Horiguchi and Tsushima [16] obtaining = 2. Lipowski’s transfer-matrix results support the = 2 prediction. Using a decimation approach Lipowski goes on to further generalize the situations looked at by considering a case where one of the intralayer interactions goes to zero as the interlayer interaction also goes to zero. This further changes the value of the shift exponent based on his decimation analysis. In particular Lipowski investigates the case where one of the intralayer interactions is set equal to the interlayer interaction and then the limit of these two interactions going to zero is examined. Because of the fact that we can obtain essentially exact critical temperatures for our Husimi tree systems for any given value of JA , JB , and IJA that is, we can obtain numerical results to arbitrary accuracy, we can obtain very accurate results for the behavior of Tc for the general situations looked at thus far in the literature, i.e., where
J.L. Monroe / Physica A 335 (2004) 563 – 576
571
JA = JB , where JA = JB but both remain 8nite while → 0, and where JB = IJA . In addition we two cases, both previously unexplored, one √ √ will consider the following 3 where JB = IJA and one where JB = IJA , where our results will further support the predictions of Lipowski’s decimation scheme. Of course our values for any critical exponents such as will not be those of the two-dimensional Ising model but those of the Husimi tree systems. Critical exponents for such systems take on their classical or mean-8eld values. One can use the approach we have developed here along with the coherent anomaly method of Suzuki [28] to estimate critical exponents of the bilayer Ising model consisting of two couple square lattice systems whose critical temperature we approximated in the previous section but this is not what we are doing in this section. Rather we investigate the behavior of the shift exponent for the Husimi tree systems themselves as various interaction strengths are allowed to go to zero. To do so we rewrite (6) as Tc () − Tc (0) = A
(7)
and with the value of Tc (0) along with Tc () for two values of one can determine and A. We should point out that for the cases JA = JB and JA = JB but both remain 8nite Hu et al. in Ref. [6] show for the bilayer Bethe lattice that in the former case = 1 and in the later = 2. We rather quickly present our results for the case where the intralayer interactions are equal, JA = JB ; as this is the case most intensely studied in the past and for which there is total agreement that = 1=. We simply verify that we obtain similar results using the classical value for . We use the values for the critical temperature listed in Table 1 for = 0:1; = 0:01; : : : ; = 0 to determine for the values of listed there. For the 8-site building block system and speci8cally the values for Tc (0:1); Tc (0:01); and Tc (0) one obtains = 0:9558, whereas with Tc (0:01); Tc (0:001), and Tc (0) one obtains = 0:9948, and 8nally with Tc (0:001), Tc (0:0001); and Tc (0) one obtains = 0:9995. Clearly then = 1=. Using similar combinations of Tc () for the 24-site building block system we obtain 0:9046, 0:9858, and 0:9985, respectively. Finally for the 48-site building block system we obtain 0:8589, 0:9741, and 0:9971. As expected all systems show in general that for this case = 1 although to get results comparable to the smallest system we would need to go to small values. Because of this we restrict ourselves to results for the 8-site building block system for the remainder of this section. For the case where JA = (1=2) JB we must use more accurate Tc values than those presented in Table 1 because of the small diHerence in Tc values when becomes small. Using Tc values accurate to 14 places past the decimal we have for the 8-site building block system = 1:99524; 1:99995, and 2:00000 for the same combinations of values used in the JA = JB case. Unfortunately in this case our results do not distinguish between the predictions of Oitmaa and Enting [4] of = 2= and the predictions of Horiguchi and Tsushima [16] and Lipowski [10] because of the fact that = 1 and hence the two predictions coincide. We now consider the more interesting case where both IJA and JB go to zero. This is our third class of situations. As stated above in Ref. [10] the case where JB = IJA was used by Lipowski as a test of his decimation procedure. This decimation procedure
572
J.L. Monroe / Physica A 335 (2004) 563 – 576
Table 3 Critical temperatures for the Husimi tree system made up of 8 site basic building blocks and for four diHerent interaction strengths for JB = JB = JB =
√
√ 3
= 0:001 2:77078256200 821975608 2:77078257369 130446337 2:77078260466 002927888
= 0:0001 2:77078256164 706435302 2:77078256168 340872204 2:77078256182 809577441
= 0:00001 2:77078256164 670374439 2:77078256164 681812202 2:77078256164 750987991
= 0:000001 2:77078256164 670338382 2:77078256164 670374498 2:77078256164 670705399
predicts = 2 for the situation in the above paragraph but more importantly can be used for a much larger class of situations. Speci8cally Lipowski 8nds by a decimation process an estimation of how the interactions IJA and JB aHect the behavior of the layer with interactions JA . It is this layer whose critical temperature dominates since we are requiring JB ¡ JA . Lipowski’s Eq. (3.4) in Ref. [10] is incorrect and does not follow from his earlier Eq. (3.3). The problem is a typographical error in the second term in both the numerator and denominator where an extra factor of 2 appears in the exponent of both terms. With these factors the series expansion of the right-hand side of his (3.4) do not give the behavior stated in Ref. [10] following Eq. (3.4). The correct equation follows: JB 1 e cosh(2IJA ) + e−JB J = : (8) ln −JB 2 e cosh(2IJA ) + eJB For small the r.h.s. of Eq. (8) goes as J ∼ 2 JA2 . Therefore, the second layer changes the 8rst layer interaction into an interaction of the form JA = 1 + C2 JA2 thereby predicting a shift exponent with a value of 2 as found in the above paragraph. However if one lets JB = IJA then expansion of the r.h.s. of Eq. (8) goes as J ∼ 3 JA3 and one gets a predicted value for the shift exponent of = 3. Lipowski checks this result using his transfer-matrix, mean-8eld method and obtains 3.05. Using our results for the 8-site building block system shown in Table 3 where we look at even smaller values than used in the preceding case and now using even more accurate numerical values for Tc of our Husimi tree systems we 8nd using Tc (0:001); Tc (0:0001); and Tc (0) one obtains = 3:00066, or Tc (0:0001); Tc (0:00001); and Tc (0) then = 3:000063. Again this is based on the Tc values given in Table 3 which are accurate to 21 8gures past the decimal and with Tc (0) = 2:77078256164670338347. However even this is not enough accuracy to use Tc (0:00001); Tc (0:000001); and Tc (0) as there is no diHerence in the value of Tc (0:000001) and Tc (0) until the 20th decimal place. We could easily get even more accurate values for Tc (0:000001) and Tc (0) so that this case can be dealt with accurately but we have not bothered to do so as clearly the 8rst two estimates of the shift exponent already clearly demonstrate that = 3 as predicted by Lipowski’s decimation approach. In addition to the above which duplicates the situations examined in Ref. [10] we have looked at two other variations which have not been looked at before namely
J.L. Monroe / Physica A 335 (2004) 563 – 576
573
√ √ 3 JB = IJA and JB = IJA . In the 8rst of these two cases expansion of the r.h.s. of (8) shows that the value of predicted by the decimation procedure is 52 . We 8nd based on the Tc values given in Table 3 that using Tc (0:001), Tc (0:0001); and Tc (0) one obtains = 2:5161, or Tc (0:0001); Tc (0:00001); and Tc (0) then = 2:5050, and 8nally Tc (0:00001), Tc (0:000001), and Tc (0) one obtains = 2:5016. In the second of these two cases expansion of the r.h.s. of (8) shows that the value of predicted by the decimation procedure is 73 . We 8nd based on the Tc values given in Table 3 that using Tc (0:001), Tc (0:0001); and Tc (0) one obtains = 2:3750, or Tc (0:0001), Tc (0:00001), and Tc (0) then = 2:3520, and for Tc (0:00001), Tc (0:000001), and Tc (0) one obtains =2:3419. In this case we have gone one step further and for Tc (0:000001), Tc (0:0000001), and Tc (0) one obtains = 2:3376 where one needs to go to 24 decimal place accuracy of the critical temperature values. In summary all our results strongly supports the view that when JA =JB then =1=, while when 0 ¡ JB ¡ JA that = 2, and that for the cases where JB and IJA both go to zero the results of the decimation procedure of Lipowski can be used to obtain the value of . 5. Antiferromagnetic interlayer interaction In the case that one has a combination of ferromagnetic and antiferromagnetic interactions one has a metamagnet. One of the reasons these systems have attracted attention is that there is often a tricritical point present in their phase diagrams and hence the phase diagram has much more structure than when all interactions are ferromagnetic or all are antiferromagnetic [29]. The most often studied metamagnet system is the simple cubic lattice where interactions along two of the axes of the lattice are ferromagnetic and those along the third axis are antiferromagnetic. Here we look at a special case where we have only two layers and consider speci8cally the intralayer interactions to be ferromagnetic and the interlayer interactions to be antiferromagnetic. In the past some others, in particular [6], have investigated bilayer systems with this arrangement of interactions but only in the case where no magnetic 8eld is present. However, it is only in the presence of a magnetic 8eld that one 8nds a tricritical point and it is then, as done with the full simple cubic metamagnet, that one wants to locate the tricritical point in the (h; T ) plane. To our knowledge no results exist for the bilayer system with h = 0. With the above mix of ferromagnetic and antiferromagnetic interactions and with h = hA = hB one has for the phase diagram in the (h; T ) plane the following. For h = 0 there exists the Neel point which indicates the presence of a continuous transition between an antiferromagnetic state and disordered state. As h is increased the critical temperature decreases but the phase transition remains continuous until the tricritical point is reached. Continuing to increase h beyond its value at the tricritical point on gets a line of 8rst order transitions all the time further decreasing the critical temperature. At some point as h is increased the critical temperature drops to zero. The exact value of h where this occurs depends on the relative strengths of the ferromagnetic and antiferromagnetic interactions.
574
J.L. Monroe / Physica A 335 (2004) 563 – 576
In terms of the dynamical system perspective the region where there is a continuous phase transition has a 8xed point behavior as described in the preceding sections. At the tricritical point there are three 8xed points, each corresponding to a state of the system, all neutral, i.e., neither attracting or repelling. In the region of the phase diagram where a 8rst order transition occurs one has the following. For high temperatures there is a single, attracting, real valued 8xed point, corresponding to the disordered phase. As the temperature is lowered this 8xed point becomes less attracting, i.e., the largest eigenvalue of the Jacobian of the dynamical system increases in value toward the value of 1 which indicates a neutral 8xed point. However, in contrast to the situation where a continuous phase transition occurs i.e., as the temperature is lowered the 8xed point corresponding to the disordered state becomes neutral and simultaneously two new attracting 8xed points are created, we have a situation where two new real valued, stable or attracting 8xed points are created before the 8xed point corresponding to the disordered state becomes neutral. Which 8xed point the system is drawn to when two or more attracting, real-valued 8xed points are present can be determined in various ways [30] but as shown in Ref. [31] one method, and seemingly the best method, is to consider the system going to the “most stable” real-valued 8xed point where by “most stable” is meant that 8xed point having the smallest, maximum eigenvalue of the Jacobian of the dynamical system. Hence for a given value of h, assuming it to be a value where a 8rst-order transition occurs, to determine the location of T at the transition point one begins by starting at a high T value, i.e., above the transition point, and lowering T while simultaneously observing the presence of all real-valued 8xed points of the system and their corresponding eigenvalues of the Jacobian. When one has located the point in the (h; T ) plane where there are three real-valued 8xed points equally stable this is the transition point. For a T above this point there will be one of these three 8xed points which is the most stable and for a T below this point there will be a pair of the three 8xed points being the most stable. This pair of 8xed points correspond to the state S+ and S− in Fig. 2 of [29]. They correspond to states created by a staggered magnetic 8eld. We consider only the case where the intralayer interactions of each layer are equal to one another. In Table 4, we give several points along the critical line in the (h; T ) plane for both the case where the absolute value of the antiferromagnetic interaction strength is equal to that of the intralayer interactions and where it’s absolute value is 1 2 the intralayer interaction strength. The location of the tricritical point is given in both cases. All points reported in this table were found using the procedures described above. 6. Conclusions In the above we have generalized a method used in the past to approximate the phase diagram of a variety of lattice spin systems but not thin 8lms. We look speci8cally at the bilayer, Ising system with nearest neighbor interactions and have shown for three diHerent areas of interest for this system, i.e., the value of the critical temperature, the scaling behavior of the critical temperature as various interactions go to zero, and the
J.L. Monroe / Physica A 335 (2004) 563 – 576
575
Table 4 Critical points in the (T ,h) plane based on the Husimi tree method using the 8 site basic building block = −1
= −1=2
T
h
T
h
0:5 1:0 1:5 2:0 2:5 3:0 3:086∗ 3:5 3:677
0:99999998 0:99991 0:9976 0:9840 0:9326 0:7913 0:762097∗ 0:4856 0:0000
0:5 1:0 1:5 2:0 2:5 3:0 3:011∗ 3:05 3:3227
0:49999998 0:4999 0:4980 0:4875 0:4472 0:3225 0:319308∗ 0:3063 0:0000
The values with the asterik are the tricritical points.
phase diagram and location of the tricritical point when an interlayer antiferromagnetic interaction and magnetic 8eld is present. In all cases the method is robust enough to give results similar in accuracy to estimates found by quite diHerent methods when these other estimates were available. The method is easy to implement and therefore in the area of the scaling behavior of the critical temperature as various interactions go to zero we looked at a number of new situations with our results supporting the earlier work of Lipowski [10]. Furthermore, in the case where the interlayer interaction is antiferromagnetic we have for the 8rst time approximations for the phase diagram and most importantly the tricritical point which is part of this phase diagram. Of course the method can be extended to consider other thin 8lm systems where either the Ising spin variable is replaced by some other spin variable which takes on a set of discrete values or one has more than two layers. An increase in the number of layers or the number of states the spin variable is allowed to take increases the dimension of the dynamical system governing the behavior of the model and hence requires more computer resources. The work here has been done on a personal computer. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
L.E. Ballantine, Physica 30 (1964) 1231. G.A.T. Allan, Phys. Rev. B 1 (1970) 352. T.W. Capehart, M.E. Fisher, Phys. Rev. B 13 (1976) 5021. J. Oitmaa, I.G. Enting, J. Phys. A: Math. Gen. 8 (1975) 1097. M.L. Lyra, C.R. da Silva, Phys. Rev. B 46 (1992) 3420. C.K. Hu, N.Sh. Izmailian, K.B. Oganesyan, Phys. Rev. E 59 (1999) 6489. D.L. Lin, H. Che, Y. Xia, Phys. Rev. A 46 (1992) 1805. K. Binder, P.C. Hohenberg, Phys. Rev. B 9 (1974) 2194. A. Lipowski, M. Suzuki, Physica A 198 (1993) 227. A. Lipowski, Physica A 250 (1998) 373. L. Angelini, D. Caroppo, M. Pellicoro, M. Villani, J. Phys. A: Math. Gen. 25 (1992) 5423. L. Angelini, D. Caroppo, M. Pellicoro, M. Villani, Physica A 219 (1995) 447.
576 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
J.L. Monroe / Physica A 335 (2004) 563 – 576
Z.B. Li, Z. Shuai, Q. Wang, H.J. Luo, L. Schulke, J. Phys. A: Math. Gen. 34 (2001) 6069. R. Abe, Prog. Theor. Phys. 44 (1970) 339. M. Suzuki, Prog. Theor. Phys. 46 (1971) 1054. T. Horiguchi, N. Tsushima, Physica A 238 (1997) 295. P.L. Hansen, J. Lemmich, H.H. Ipsen, O.G. Mouritsen, J. Stat. Phys. 73 (1993) 723. J.L. Monroe, Physica A 206 (1994) 218. J.L. Monroe, J. Stat. Phys. 65 (1992) 255. J.L. Monroe, J. Stat. Phys. 67 (1992) 1185. J.L. Monroe, Phys. Rev. E 66 (2002) 066129. J.L. Monroe, Phys. Rev. E 64 (2001) 016126. J.W. Essam, M.E. Fisher, Rev. Mod. Phys. 42 (1970) 272. R. Bulirsch, J. Stoer, Numer. Math. 6 (1964) 413. M. Henkel, A. Patkos, J. Phys. A: Math. Gen. 20 (1987) 2199. M. Henkel, G. Schutz, J. Phys. A: Math. Gen. 21 (1988) 2617. J.L. Monroe, Phys. Rev. E 65 (2002) 066116. M. Suzuki, Coherent Anomaly Method, World Scienti8c, Singapore, 1995. I.D. Lawrie, S. Sarbach, In: C. Domb, J.L. Lebowitz, (Eds.), Phase Transitions and Critical Phenomena, Vol. 9, Academic Press, New York, 1984. [30] M. Pretti, J. Stat. Phys. 111 (2003) 993. [31] J.L. Monroe, Phys. Lett. A 188 (1994) 80.
Physica A 335 (2004) 563 – 576
www.elsevier.com/locate/physa
The bilayer Ising model and a generalized Husimi tree approximation James L. Monroe Department of Physics, Penn State University, Beaver Campus, 100 University Dr., Monaca, PA 15061, USA Received 10 September 2003
Abstract We study a variety of aspects of a bilayer system of Ising spins including accurate estimates of the critical temperature for ferromagnetic interactions, scaling of the critical temperature when the interlayer interaction goes to zero, and approximations of phase diagrams for the case when antiferromagnetic interlayer interactions are present including location of the tricritical point. This is all done using an extension of a method previously developed by the author to study a variety of lattice spin systems and which is a generalization of the Bethe lattice approach. c 2003 Elsevier B.V. All rights reserved. PACS: 05.50.+q; 64.60.Cn; 75.70.−i Keywords: Ising spins; Thin 8lms; Bilayer system
1. Introduction For nearly four decades models of the magnetic properties of thin 8lms have been of interest both in terms of various aspects of practical applications as well as fundamental theoretical issues. The simplest model of a thin 8lm, and the 8rst to be studied, is the bilayer system consisting of two layers of Ising spins with nearest neighbor, intraplane interactions, and interplane interactions. In 1964, Ballantine [1] using high-temperature series expansion methods obtained an estimate for the critical temperature of such a system where all interactions were ferromagnetic and of the same strength. Since that time there has been continual interest in obtaining better estimates of the critical temperature for such a system through lengthier series expansions [2,3], analytical results based on mean-8eld theory [4], the Bethe approximation [5,6] and a variational E-mail address: [email protected] (J.L. Monroe). c 2003 Elsevier B.V. All rights reserved. 0378-4371/$ - see front matter doi:10.1016/j.physa.2003.12.018
564
J.L. Monroe / Physica A 335 (2004) 563 – 576
cumulant expansion [7] as well as many other numerical methods including Monte Carlo [8], interfacial tensions approximation [9], a transfer-matrix version of mean-8eld [9,10], a constrained variational approach and correlation length equality [11,12], and a corner-transfer-matrix, renormalization group approach [13]. Additionally the behavior of the critical temperature of such a system in the limit of the interplane interaction going to zero and to a lesser extent the behavior when this interaction approaches in8nity has been studied [4,10,14–17]. Finally the general phase diagram has a rich structure especially when either or both the interlayer and intralayer interactions are allowed to be antiferromagnetic. In the following, we use a method already used to obtain results for frustrated systems such as the Ising antiferromagnet on the triangle lattice [18], Ising systems with multisite interactions [19,20], higher-spin Ising and Potts systems [21], and in its most developed form for the Ising antiferromagnet and ferromagnet on the square lattice [22]. Here it will be used in all three areas mentioned above. First, for the case of equal ferromagnetic intraplaner interactions where the phase diagram is simply a line of phase transitions in zero magnetic 8eld ending in at a critical endpoint (the emphasis in this case is to accurately determine the critical temperature). Second, to examine the behavior of the system in the limit of vanishing interplaner interactions both for the case where the intraplaner interactions on both layers are equal and when they are not as well as when not only the interlayer interaction goes to zero but one of the intraplaner interactions does so as well. Third, to approximate the phase diagram for the bilayer system having an antiferromagnetic interplaner interaction along with ferromagnetic intraplaner interactions when there is a magnetic 8eld present. In this case, we have an example of a metamagnet and there is present a tricritical point. Basically the method we will use is a generalization of the Bethe lattice approach. The Bethe lattice can be thought of as being built in a recursive manner by attaching to a central site on which there is an Ising spin, q Ising spins each interacting with the central site through a pair interaction. Then to each of these q spins are attached q − 1 spins each interacting with one of the q spins through a pair interaction. One continues in this manner ad in8nitum to construct the full Bethe lattice. Of course this will not work when one wants to approximate a bilayer system but one can easily imagine two identical Bethe lattices coupled together with an interplaner interaction as an approximation for the bilayer system. This is the system that has been looked at in some detail very recently by Hu et al. [6]. The basic building block in the construction of a Bethe lattice is a pair of nearest neighbor Ising spins and the interaction connecting them. If one wants to approximate a square lattice one can take q = 4. However, a better approximation of the square lattice can be obtained by taking as a basic building block a square on each corner of which is an Ising spin and each Ising spin interacting with its two nearest neighbors. One can take one of these as a central square and consider attaching a similar square to each corner of the central square. One can then take this collection of sites and to each Ising spin not already having four nearest neighbor interactions attach another of the 4-site square building blocks. Continuing this ad in8nitum one constructs what is denoted by graph theorists as an Husimi tree [23]. The behavior of a spin deep inside the tree can be used as an approximation for the spontaneous magnetization
J.L. Monroe / Physica A 335 (2004) 563 – 576
565
of the original square lattice system just as a site deep inside the Bethe lattice can be used as an approximation of the original square lattice system. Since the Husimi tree approximation more correctly takes into account local Fuctuations than the Bethe lattice approximation one 8nds it more correctly approximates the original system. Of course as described for the Bethe lattice approximation of the bilayer system where two Bethe lattices were connected with interplaner interactions the same can be done for the Husimi tree approximation. Furthermore, there is no need to consider only a 4-site square as the basic building block larger building blocks can be used including ones where rather than connections being made at a single site connections are made at a collection of sites. In Ref. [22] this was done using basic building blocks of as many as 60 sites to approximate the square lattice ferromagnetic Ising model as well as the antiferromagnetic case where a magnetic 8eld is present. In the ferromagnetic case the exact value of the critical temperature is known and results as close as 0.003% of the actual value were obtained. In the following section, the Hamiltonian of the bilayer system is presented as well as a more detailed description of the method used here to study the bilayer Ising model. Then Sections 3–5 contain the results of the three diHerent areas mentioned above. We close in Section 6 with some concluding remarks. 2. Basic method In the following, we approximate a variety of properties of a system consisting of two parallel, horizontal, square lattice Ising models, where on the upper layer the Ising spin on the ith site will be denoted by i and on the lower layer by i , and where by Ising spin we mean speci8cally a spin allowed to take on the values ±1. The Hamiltonian of the system is H = −JA i j − J B i j − IJA i i − h A i − h B i ; (1) ¡i; j¿
¡i; j¿
¡i; j¿
i
i
where ¡ i; j ¿ represents a nearest neighbor pair of sites, JA and JB are, respectively, the intralayer interactions on the upper and lower layers of the two-layer system, IJA is the interlayer interaction, and hA and hB are, respectively, the external magnetic 8eld present on the upper and lower layers of the system. As stated above a single layer of the system can be approximated by a Husimi tree approach where one takes as a basic building block a set of four sites, each site on the corner of a square, and on each site an Ising spin. One can start with a central square and attach to it at each corner a new square. Then taking this arrangement one can attach to each site not having four nearest neighbors a new square. One can proceed in this way, in the limit building an in8nite Husimi tree. This process because of the step-by-step nature of the construction can be described by a recursive system of equations, i.e., a dynamical system, which can be used to calculate, for example, the thermal average of a site of the central square. Based on the behavior of this thermal average one can determine a critical temperature. In terms of dynamical systems language one has for this system a one-dimensional dynamical system and for high
566
J.L. Monroe / Physica A 335 (2004) 563 – 576
temperatures the behavior of the system is governed by a single attracting 8xed point which corresponds to zero magnetization of any one of the sites of the central square. For the case where a continuous phase transition occurs as the temperature is lowered this 8xed point becomes repelling and two new attracting 8xed points are created, i.e., a bifurcation occurs. These two new 8xed points correspond to nonzero thermal averages of a site on the central square, one positive and negative, and hence a spontaneous magnetization. Therefore, the temperature at which the single 8xed point becomes repelling and two new attracting 8xed points are created is the critical temperature. See reference [22] for details. Now in analogy when trying to approximate the bilayer system one can imagine constructing a tree-like structure with a basic building block consisting of two parallel, horizontal squares of four sites each and interacting via nearest neighbor intralayer interactions JA and JB . In addition one has an interlayer, vertical interaction, between nearest neighbor sites in the vertical direction. Hence, one has a basic building block which is a cube consisting of eight sites. Now one can construct a system by considering a central cube and attaching to each vertical edge of this central cube a basic building block. One sees now that the vertical interaction in each building block, i.e., the interlayer interaction, should be taken as only one-half of IJA for the interaction strength. This is because when making the attachments the vertical interactions get a contribution from that of the original building block as well as the new building block the sum of which is the full IJA . Since each connection consists of two sites, one on the upper layer and one on the lower layer, the dynamical system which results is a three-dimensional dynamical system. The dimension of the dynamical system, disregarding any symmetry considerations, is always 2n − 1, where n is the number of sites involved in a connection, again see Ref. [22] for details. The mathematical behavior of the resulting three-dimensional dynamical system (here n = 2), i.e., the presence of attracting 8xed points, two-cycles, etc., governs the physical behavior of the system. For the case where JA , JB , and IJA are positive, i.e., ferromagnetic, and the external magnetic 8elds are zero then based on physical grounds one expects for high temperatures a single attracting 8xed point which bifurcates to two new attracting 8xed points at some point as the temperature is lowered as described previously. This is exactly what one 8nds occurring for this case. As with the case of the square lattice Ising model approximated in Ref. [22] one can consider larger building blocks. In fact in Ref. [22] a very systematic sequence of 8ve diHerent building blocks is used to construct various Husimi trees and the results based on them were presented in the approximation of the square lattice Ising model. For approximating, the bilayer system two of each of these basic building blocks can be placed in parallel, the upper one with nearest neighbor interactions JA and the lower one with interactions JB . Then the two can be coupled to each other with interactions IJA for sites not involved in a connection and one-half IJA for sites on the edge of the building block and hence involved in a connection. Of course the dimension of the dynamical system grows rapidly. To use the 8fth of the 8ve building blocks as done in the square lattice case for this bilayer case would mean 10 sites were involved in each connection resulting in a dynamical system of dimension 210 − 1, which is beyond the computer capabilities used here.
J.L. Monroe / Physica A 335 (2004) 563 – 576
567
All computations have been executed on a personal computer using Mathematica to actually generate the dynamical systems governing the various approximations. Because of this for the bilayer system being approximated here we have only been able to obtain three levels of approximation rather than the 8ve for the square lattice. Even so as will be seen in the next section we are able to obtain reasonably accurate estimates of the critical temperature where by “reasonably” we mean estimates having approximately the same accuracy as obtained by presently available series expansion methods. To obtain the best approximations for the critical behavior we have used our systematic sequence of three approximations along with the Bulirsch-Stoer extrapolation procedure [24–27]. Here one wishes to determine the value of Tc of the bilayer system by knowing the corresponding values of Tc based on some sequence of approximations using increasingly larger building blocks. The value of Tc using the approximations will depend on the size of the basic building blocks used and one denotes them as Tc (L) where L is a measure of the system size. One supposes in general that Tc (L) can be written as Tc (L) = Tc; ∞ + a1 L−!1 + a2 L−!2 + a3 L−!3 + · · · ;
(2)
where 0 ¡ !1 ¡ !2 ¡ !3 ¡ · · · and where Tc; ∞ is the value for the system being approximated. The BST algorithm allows one to start with a sequence of values, here three values Tc (L1 ), Tc (L2 ), and Tc (L3 ) and estimate Tc; ∞ . In particular the algorithm allows one to construct a table of extrapolants of this sequence T01 T11 T02
T21 T12
T03 where T21 is the algorithm’s best estimate for Tc; ∞ . The Tqn are computed from n =0 ; T−1
(3)
T0n = Tc (n) ;
(4)
Tmn
=
n+1 Tm−1
+
n+1 (Tm−1
−
n Tm−1 )
Ln Ln+m
−1 ! n+1 n Tm−1 − Tm−1 1 − n+1 −1 ; n+1 Tm−1 − Tm−2
(5)
where m ¿ 1 and where ! is a free parameter. 3. Estimates of the critical temperature The values of Tc (L) for the three generalized Husimi tree systems used here are given in Table 1 for a variety of combinations of JA , JB ; and IJA with the restriction that all three are ferromagnetic. For simplicity we set JA =1 and then vary the values of
568
J.L. Monroe / Physica A 335 (2004) 563 – 576
Table 1 Critical temperature estimates based on the Bethe lattice as well as the three levels of Husimi trees described in the text
JB JB JB JB JB JB JB JB JB JB JB JB JB
= 1, = 4 = 1, = 2 = 1, = 1 = 1, = 0:5 = 1, = 0:1 = 1, = 0:01 = 1, = 0:001 = 1, = 0:0001 = 1, = 0 = 1=2, = 0:1 = 1=2, = 0:01 = 1=2, = 0:001 = 1=2, = 0:0001
Bethe lattice, results
8-site results
24-site results
48-site results
4.86557 4.20205 3.68483 3.33961 2.98923 2.89616 2.88647 2.88550 2.88539
4.784395068 4.175811973 3.677368547 3.322723648 2.915324057 2.786786000 2.772402380 2.770944746 2.770782562 2.776106370 2.770836387 2.770783103 2.770782568
4.512328839 3.957293525 3.491491794 3.153858714 2.742398308 2.588782964 2.569183645 2.567151728 2.566925090
4.396933786 3.864416292 3.412610923 3.083346974 2.675481529 2.509447499 2.485612169 2.483066419 2.482781460
JB , and . The three systems used here are the generalizations of the 8rst three systems used in Ref. [22] to approximate the square lattice, Ising model to the bilayer system as described above. Hence, they consist of eight sites (the cube described in detail above), a 24-site system consisting of two parallel 12-site building blocks (see Ref. [12]) connected by interlayer interactions, and 8nally a 48-site system consisting of two parallel 24-site building blocks (see Ref. [12]) connected by interlayer interactions. In addition to the Husimi tree estimates in Table 1 we give the estimates found from the Bethe lattice approximation [5,6] which in some sense is a zeroth-order approximation if we consider our three estimates as 8rst, second, and third-order estimates. Both Refs. [5,6] give analytical expressions for the critical temperature although the expression derived in Ref. [6] requires JA = JB . Unfortunately the critical temperatures given by the two expressions diHer except at the two extremes of = 0 and → ∞. We have using the dynamical systems approach applied to our three levels of approximation found critical temperature values for the Bethe lattice which match those given by the analytical expressions of Hu et al. in Ref. [6] and therefore use these values as our Bethe lattice results listed in Table 1. We have presented our estimates of Tc to nine places past the decimal as this level of precision is necessary for the use of the BST algorithm which will be used to get our best estimates of the critical temperature. As expected one sees a systematic decrease in the estimated value of Tc as the size of the basic building blocks is increased. The decrease in the estimated value of Tc corresponds to an increase in the accuracy of the approximations as increasing the size of the basic building blocks allows the system to more properly take into account local Fuctuations. As is the case for systems with short range interactions such as the nearest neighbor interactions present in the system under study here the estimates of Tc found using these approaches are typically too high. Using the three estimates based on the Husimi tree approach given in Table I along with the BST extrapolation algorithm we arrive at our best estimates for Tc . These
J.L. Monroe / Physica A 335 (2004) 563 – 576
569
Table 2 Critical temperature estimates based on the BST algorithm and the three Husimi tree results presented in Table 1 along with estimates based on high temperature series expansions [3,4], constrained variational method and correlation length equality [12], transfer matrix, mean-8eld [10], the corner transfer matrix method [13], and Monte Carlo BST estimates JB JB JB JB JB JB JB JB JB JB
= 1, = 1, = 1, = 1, = 1, = 1, = 1, = 1, = 1, = 1,
I=8 =4 =2 =1 = 0:5 = 0:1 = 0:01 = 0:001 = 0:0001 =0
4.095478 3.621238 3.206424 2.902182 2.516101 2.314020 2.274273 2.269701 2.269185
High temp. series exp.
3:623 ± 0:015 3:208 ± 0:007 2:911 ± 0:006 2:509 ± 0:010 2:307 ± :011 2:271 ± 0:005
CVM & CLE 4.4458 4.0968 3.6233 3.2085 2.5150 2.3350
Trans. mat. & MF
3.6227 3.207530 2.901246 2.529991
Corner trans. mat.
3.20762 2.51429 2.33322
Monte Carlo 4.098 3.610 3.226 2.915 2.532 2.299
2.26917
are given in Table 2 along with estimates based on several other approaches. In the use of the BST algorithm one must choose a value for ! as mentioned following Eq. (5). This is a critical issue. In our case we can use the fact that when = 0 one has Onsager’s exact result and using our three estimates for Tc when = 0 we vary ! determining the value which gives to the numerical accuracy used here Onsager’s result. This occurs when ! = 0:940486. We use this value of ! for the cases where = 0. Also, L is a measure of the system size but the proper value to use here is not entirely clear. Two possibilities are readily apparent one is to take L as the square root of the number of sites in a basic building block and the other is to use the number of sites used in each of the connections made in the construction of the tree, i.e., L = 2, 4, or 6. We have chosen to use the later. The most studied case is the one where all three interactions JA , JB ; and IJA are of equal strength. Lipowski and Suzuki [9] 8nd using the transfer matrix, mean-8eld approach that Tc 3:20755(5). The method is similar to the one presented here in two major ways 8rst it is similar in that one can look at larger and larger systems obtaining increasingly accurate results and second one thereby develops a sequence of approximations which can be extrapolated to get a 8nal best estimate. Lipowski [10] extends the earlier work to obtain Tc 3:207530. These results are best supported by the results of Li et al. [13], which diHer from that of Lipowski by less than 0:003%. Our results underestimate the critical temperature by approximately as much as the CVM and CLE results of Ref. [12] over estimate it, with a diHerence of 0:03% with the Lipowski result. The Monte Carlo results from Ref. [17] give a critical temperature estimate Tc 3:226 which is signi8cantly above our estimate as well as the other four estimates presented in Table 2. As one decreases the value of the interlayer interaction the situation in terms of accuracy of the estimates is not at all clear. For example for = 0:1 from Table 2 we see that the results of Ref. [10] are higher than those of all other method of
570
J.L. Monroe / Physica A 335 (2004) 563 – 576
approximation except the Monte Carlo results of Ref. [17] and are in fact outside the error bounds presented with the high temperature series expansion results. The results of Refs. [12,13], and our own give for = 0:1 a best result of Tc 2:515(1). With further reduction of the value to = 0:01 one sees the corner matrix results or the CVM and CLE results are now greater than that obtained by either the estimates of the Husimi tree, BST extrapolation method or series expansions method with the former estimates even outside the error bars of the series expansion results. Unfortunately numerical values for Tc with still further reduction of have not been given in other references. For values of above 1.0 only for = 2 are numerical values for the estimate of Tc generally available. Five results are presented in Table 2 and there is general agreement amongst four of the results with now the Monte Carlo results being signi8cantly diHerent on the low side. In general there is little systematics comparing one method to another so good estimates of the errors are diPcult. 4. Tc as → 0 The behavior of Tc as → 0 has been of interest since 1970 when Abe [14] in terms of a scaling theory derived a value for the shift exponent which describes the deviation of the critical temperature Tc () from the critical temperature in the decoupled limit, i.e., Tc (0), with de8ned through the following: Tc () − Tc (0) ∼ :
(6)
Suzuki [15] followed up Abe’s work with further analysis. Both predict that when the intralayer interactions, JA and JB ; are equal that the shift exponent is equal to 1=, being the critical exponent of the susceptibility of the two-dimensional model. Oitmaa and Enting [4] generalized to the situation where the intralayer interactions are not equal and found for this situation that = 2=. As Lipowski [10] emphasizes this is particularly interesting because we now are coupling two systems which have diHerent critical temperatures. Lipowski also points out that as of 1998, the time of his paper, these scaling predictions had only modest con8rmation. In fact for the case of JA = JB contradictory predictions existed with Horiguchi and Tsushima [16] obtaining = 2. Lipowski’s transfer-matrix results support the = 2 prediction. Using a decimation approach Lipowski goes on to further generalize the situations looked at by considering a case where one of the intralayer interactions goes to zero as the interlayer interaction also goes to zero. This further changes the value of the shift exponent based on his decimation analysis. In particular Lipowski investigates the case where one of the intralayer interactions is set equal to the interlayer interaction and then the limit of these two interactions going to zero is examined. Because of the fact that we can obtain essentially exact critical temperatures for our Husimi tree systems for any given value of JA , JB , and IJA that is, we can obtain numerical results to arbitrary accuracy, we can obtain very accurate results for the behavior of Tc for the general situations looked at thus far in the literature, i.e., where
J.L. Monroe / Physica A 335 (2004) 563 – 576
571
JA = JB , where JA = JB but both remain 8nite while → 0, and where JB = IJA . In addition we two cases, both previously unexplored, one √ √ will consider the following 3 where JB = IJA and one where JB = IJA , where our results will further support the predictions of Lipowski’s decimation scheme. Of course our values for any critical exponents such as will not be those of the two-dimensional Ising model but those of the Husimi tree systems. Critical exponents for such systems take on their classical or mean-8eld values. One can use the approach we have developed here along with the coherent anomaly method of Suzuki [28] to estimate critical exponents of the bilayer Ising model consisting of two couple square lattice systems whose critical temperature we approximated in the previous section but this is not what we are doing in this section. Rather we investigate the behavior of the shift exponent for the Husimi tree systems themselves as various interaction strengths are allowed to go to zero. To do so we rewrite (6) as Tc () − Tc (0) = A
(7)
and with the value of Tc (0) along with Tc () for two values of one can determine and A. We should point out that for the cases JA = JB and JA = JB but both remain 8nite Hu et al. in Ref. [6] show for the bilayer Bethe lattice that in the former case = 1 and in the later = 2. We rather quickly present our results for the case where the intralayer interactions are equal, JA = JB ; as this is the case most intensely studied in the past and for which there is total agreement that = 1=. We simply verify that we obtain similar results using the classical value for . We use the values for the critical temperature listed in Table 1 for = 0:1; = 0:01; : : : ; = 0 to determine for the values of listed there. For the 8-site building block system and speci8cally the values for Tc (0:1); Tc (0:01); and Tc (0) one obtains = 0:9558, whereas with Tc (0:01); Tc (0:001), and Tc (0) one obtains = 0:9948, and 8nally with Tc (0:001), Tc (0:0001); and Tc (0) one obtains = 0:9995. Clearly then = 1=. Using similar combinations of Tc () for the 24-site building block system we obtain 0:9046, 0:9858, and 0:9985, respectively. Finally for the 48-site building block system we obtain 0:8589, 0:9741, and 0:9971. As expected all systems show in general that for this case = 1 although to get results comparable to the smallest system we would need to go to small values. Because of this we restrict ourselves to results for the 8-site building block system for the remainder of this section. For the case where JA = (1=2) JB we must use more accurate Tc values than those presented in Table 1 because of the small diHerence in Tc values when becomes small. Using Tc values accurate to 14 places past the decimal we have for the 8-site building block system = 1:99524; 1:99995, and 2:00000 for the same combinations of values used in the JA = JB case. Unfortunately in this case our results do not distinguish between the predictions of Oitmaa and Enting [4] of = 2= and the predictions of Horiguchi and Tsushima [16] and Lipowski [10] because of the fact that = 1 and hence the two predictions coincide. We now consider the more interesting case where both IJA and JB go to zero. This is our third class of situations. As stated above in Ref. [10] the case where JB = IJA was used by Lipowski as a test of his decimation procedure. This decimation procedure
572
J.L. Monroe / Physica A 335 (2004) 563 – 576
Table 3 Critical temperatures for the Husimi tree system made up of 8 site basic building blocks and for four diHerent interaction strengths for JB = JB = JB =
√
√ 3
= 0:001 2:77078256200 821975608 2:77078257369 130446337 2:77078260466 002927888
= 0:0001 2:77078256164 706435302 2:77078256168 340872204 2:77078256182 809577441
= 0:00001 2:77078256164 670374439 2:77078256164 681812202 2:77078256164 750987991
= 0:000001 2:77078256164 670338382 2:77078256164 670374498 2:77078256164 670705399
predicts = 2 for the situation in the above paragraph but more importantly can be used for a much larger class of situations. Speci8cally Lipowski 8nds by a decimation process an estimation of how the interactions IJA and JB aHect the behavior of the layer with interactions JA . It is this layer whose critical temperature dominates since we are requiring JB ¡ JA . Lipowski’s Eq. (3.4) in Ref. [10] is incorrect and does not follow from his earlier Eq. (3.3). The problem is a typographical error in the second term in both the numerator and denominator where an extra factor of 2 appears in the exponent of both terms. With these factors the series expansion of the right-hand side of his (3.4) do not give the behavior stated in Ref. [10] following Eq. (3.4). The correct equation follows: JB 1 e cosh(2IJA ) + e−JB J = : (8) ln −JB 2 e cosh(2IJA ) + eJB For small the r.h.s. of Eq. (8) goes as J ∼ 2 JA2 . Therefore, the second layer changes the 8rst layer interaction into an interaction of the form JA = 1 + C2 JA2 thereby predicting a shift exponent with a value of 2 as found in the above paragraph. However if one lets JB = IJA then expansion of the r.h.s. of Eq. (8) goes as J ∼ 3 JA3 and one gets a predicted value for the shift exponent of = 3. Lipowski checks this result using his transfer-matrix, mean-8eld method and obtains 3.05. Using our results for the 8-site building block system shown in Table 3 where we look at even smaller values than used in the preceding case and now using even more accurate numerical values for Tc of our Husimi tree systems we 8nd using Tc (0:001); Tc (0:0001); and Tc (0) one obtains = 3:00066, or Tc (0:0001); Tc (0:00001); and Tc (0) then = 3:000063. Again this is based on the Tc values given in Table 3 which are accurate to 21 8gures past the decimal and with Tc (0) = 2:77078256164670338347. However even this is not enough accuracy to use Tc (0:00001); Tc (0:000001); and Tc (0) as there is no diHerence in the value of Tc (0:000001) and Tc (0) until the 20th decimal place. We could easily get even more accurate values for Tc (0:000001) and Tc (0) so that this case can be dealt with accurately but we have not bothered to do so as clearly the 8rst two estimates of the shift exponent already clearly demonstrate that = 3 as predicted by Lipowski’s decimation approach. In addition to the above which duplicates the situations examined in Ref. [10] we have looked at two other variations which have not been looked at before namely
J.L. Monroe / Physica A 335 (2004) 563 – 576
573
√ √ 3 JB = IJA and JB = IJA . In the 8rst of these two cases expansion of the r.h.s. of (8) shows that the value of predicted by the decimation procedure is 52 . We 8nd based on the Tc values given in Table 3 that using Tc (0:001), Tc (0:0001); and Tc (0) one obtains = 2:5161, or Tc (0:0001); Tc (0:00001); and Tc (0) then = 2:5050, and 8nally Tc (0:00001), Tc (0:000001), and Tc (0) one obtains = 2:5016. In the second of these two cases expansion of the r.h.s. of (8) shows that the value of predicted by the decimation procedure is 73 . We 8nd based on the Tc values given in Table 3 that using Tc (0:001), Tc (0:0001); and Tc (0) one obtains = 2:3750, or Tc (0:0001), Tc (0:00001), and Tc (0) then = 2:3520, and for Tc (0:00001), Tc (0:000001), and Tc (0) one obtains =2:3419. In this case we have gone one step further and for Tc (0:000001), Tc (0:0000001), and Tc (0) one obtains = 2:3376 where one needs to go to 24 decimal place accuracy of the critical temperature values. In summary all our results strongly supports the view that when JA =JB then =1=, while when 0 ¡ JB ¡ JA that = 2, and that for the cases where JB and IJA both go to zero the results of the decimation procedure of Lipowski can be used to obtain the value of . 5. Antiferromagnetic interlayer interaction In the case that one has a combination of ferromagnetic and antiferromagnetic interactions one has a metamagnet. One of the reasons these systems have attracted attention is that there is often a tricritical point present in their phase diagrams and hence the phase diagram has much more structure than when all interactions are ferromagnetic or all are antiferromagnetic [29]. The most often studied metamagnet system is the simple cubic lattice where interactions along two of the axes of the lattice are ferromagnetic and those along the third axis are antiferromagnetic. Here we look at a special case where we have only two layers and consider speci8cally the intralayer interactions to be ferromagnetic and the interlayer interactions to be antiferromagnetic. In the past some others, in particular [6], have investigated bilayer systems with this arrangement of interactions but only in the case where no magnetic 8eld is present. However, it is only in the presence of a magnetic 8eld that one 8nds a tricritical point and it is then, as done with the full simple cubic metamagnet, that one wants to locate the tricritical point in the (h; T ) plane. To our knowledge no results exist for the bilayer system with h = 0. With the above mix of ferromagnetic and antiferromagnetic interactions and with h = hA = hB one has for the phase diagram in the (h; T ) plane the following. For h = 0 there exists the Neel point which indicates the presence of a continuous transition between an antiferromagnetic state and disordered state. As h is increased the critical temperature decreases but the phase transition remains continuous until the tricritical point is reached. Continuing to increase h beyond its value at the tricritical point on gets a line of 8rst order transitions all the time further decreasing the critical temperature. At some point as h is increased the critical temperature drops to zero. The exact value of h where this occurs depends on the relative strengths of the ferromagnetic and antiferromagnetic interactions.
574
J.L. Monroe / Physica A 335 (2004) 563 – 576
In terms of the dynamical system perspective the region where there is a continuous phase transition has a 8xed point behavior as described in the preceding sections. At the tricritical point there are three 8xed points, each corresponding to a state of the system, all neutral, i.e., neither attracting or repelling. In the region of the phase diagram where a 8rst order transition occurs one has the following. For high temperatures there is a single, attracting, real valued 8xed point, corresponding to the disordered phase. As the temperature is lowered this 8xed point becomes less attracting, i.e., the largest eigenvalue of the Jacobian of the dynamical system increases in value toward the value of 1 which indicates a neutral 8xed point. However, in contrast to the situation where a continuous phase transition occurs i.e., as the temperature is lowered the 8xed point corresponding to the disordered state becomes neutral and simultaneously two new attracting 8xed points are created, we have a situation where two new real valued, stable or attracting 8xed points are created before the 8xed point corresponding to the disordered state becomes neutral. Which 8xed point the system is drawn to when two or more attracting, real-valued 8xed points are present can be determined in various ways [30] but as shown in Ref. [31] one method, and seemingly the best method, is to consider the system going to the “most stable” real-valued 8xed point where by “most stable” is meant that 8xed point having the smallest, maximum eigenvalue of the Jacobian of the dynamical system. Hence for a given value of h, assuming it to be a value where a 8rst-order transition occurs, to determine the location of T at the transition point one begins by starting at a high T value, i.e., above the transition point, and lowering T while simultaneously observing the presence of all real-valued 8xed points of the system and their corresponding eigenvalues of the Jacobian. When one has located the point in the (h; T ) plane where there are three real-valued 8xed points equally stable this is the transition point. For a T above this point there will be one of these three 8xed points which is the most stable and for a T below this point there will be a pair of the three 8xed points being the most stable. This pair of 8xed points correspond to the state S+ and S− in Fig. 2 of [29]. They correspond to states created by a staggered magnetic 8eld. We consider only the case where the intralayer interactions of each layer are equal to one another. In Table 4, we give several points along the critical line in the (h; T ) plane for both the case where the absolute value of the antiferromagnetic interaction strength is equal to that of the intralayer interactions and where it’s absolute value is 1 2 the intralayer interaction strength. The location of the tricritical point is given in both cases. All points reported in this table were found using the procedures described above. 6. Conclusions In the above we have generalized a method used in the past to approximate the phase diagram of a variety of lattice spin systems but not thin 8lms. We look speci8cally at the bilayer, Ising system with nearest neighbor interactions and have shown for three diHerent areas of interest for this system, i.e., the value of the critical temperature, the scaling behavior of the critical temperature as various interactions go to zero, and the
J.L. Monroe / Physica A 335 (2004) 563 – 576
575
Table 4 Critical points in the (T ,h) plane based on the Husimi tree method using the 8 site basic building block = −1
= −1=2
T
h
T
h
0:5 1:0 1:5 2:0 2:5 3:0 3:086∗ 3:5 3:677
0:99999998 0:99991 0:9976 0:9840 0:9326 0:7913 0:762097∗ 0:4856 0:0000
0:5 1:0 1:5 2:0 2:5 3:0 3:011∗ 3:05 3:3227
0:49999998 0:4999 0:4980 0:4875 0:4472 0:3225 0:319308∗ 0:3063 0:0000
The values with the asterik are the tricritical points.
phase diagram and location of the tricritical point when an interlayer antiferromagnetic interaction and magnetic 8eld is present. In all cases the method is robust enough to give results similar in accuracy to estimates found by quite diHerent methods when these other estimates were available. The method is easy to implement and therefore in the area of the scaling behavior of the critical temperature as various interactions go to zero we looked at a number of new situations with our results supporting the earlier work of Lipowski [10]. Furthermore, in the case where the interlayer interaction is antiferromagnetic we have for the 8rst time approximations for the phase diagram and most importantly the tricritical point which is part of this phase diagram. Of course the method can be extended to consider other thin 8lm systems where either the Ising spin variable is replaced by some other spin variable which takes on a set of discrete values or one has more than two layers. An increase in the number of layers or the number of states the spin variable is allowed to take increases the dimension of the dynamical system governing the behavior of the model and hence requires more computer resources. The work here has been done on a personal computer. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
L.E. Ballantine, Physica 30 (1964) 1231. G.A.T. Allan, Phys. Rev. B 1 (1970) 352. T.W. Capehart, M.E. Fisher, Phys. Rev. B 13 (1976) 5021. J. Oitmaa, I.G. Enting, J. Phys. A: Math. Gen. 8 (1975) 1097. M.L. Lyra, C.R. da Silva, Phys. Rev. B 46 (1992) 3420. C.K. Hu, N.Sh. Izmailian, K.B. Oganesyan, Phys. Rev. E 59 (1999) 6489. D.L. Lin, H. Che, Y. Xia, Phys. Rev. A 46 (1992) 1805. K. Binder, P.C. Hohenberg, Phys. Rev. B 9 (1974) 2194. A. Lipowski, M. Suzuki, Physica A 198 (1993) 227. A. Lipowski, Physica A 250 (1998) 373. L. Angelini, D. Caroppo, M. Pellicoro, M. Villani, J. Phys. A: Math. Gen. 25 (1992) 5423. L. Angelini, D. Caroppo, M. Pellicoro, M. Villani, Physica A 219 (1995) 447.
576 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
J.L. Monroe / Physica A 335 (2004) 563 – 576
Z.B. Li, Z. Shuai, Q. Wang, H.J. Luo, L. Schulke, J. Phys. A: Math. Gen. 34 (2001) 6069. R. Abe, Prog. Theor. Phys. 44 (1970) 339. M. Suzuki, Prog. Theor. Phys. 46 (1971) 1054. T. Horiguchi, N. Tsushima, Physica A 238 (1997) 295. P.L. Hansen, J. Lemmich, H.H. Ipsen, O.G. Mouritsen, J. Stat. Phys. 73 (1993) 723. J.L. Monroe, Physica A 206 (1994) 218. J.L. Monroe, J. Stat. Phys. 65 (1992) 255. J.L. Monroe, J. Stat. Phys. 67 (1992) 1185. J.L. Monroe, Phys. Rev. E 66 (2002) 066129. J.L. Monroe, Phys. Rev. E 64 (2001) 016126. J.W. Essam, M.E. Fisher, Rev. Mod. Phys. 42 (1970) 272. R. Bulirsch, J. Stoer, Numer. Math. 6 (1964) 413. M. Henkel, A. Patkos, J. Phys. A: Math. Gen. 20 (1987) 2199. M. Henkel, G. Schutz, J. Phys. A: Math. Gen. 21 (1988) 2617. J.L. Monroe, Phys. Rev. E 65 (2002) 066116. M. Suzuki, Coherent Anomaly Method, World Scienti8c, Singapore, 1995. I.D. Lawrie, S. Sarbach, In: C. Domb, J.L. Lebowitz, (Eds.), Phase Transitions and Critical Phenomena, Vol. 9, Academic Press, New York, 1984. [30] M. Pretti, J. Stat. Phys. 111 (2003) 993. [31] J.L. Monroe, Phys. Lett. A 188 (1994) 80.