Critical behavior and phase diagrams of a spin-1 Blume-Capel model with random crystal field interactions: An effective field theory analysis
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Critical behavior and phase diagrams of a spin-1 Blume-Capel model with randomcrystal field interactions: An effective field theory analysis
Yusuf Y¨uksel, ∗ ¨Umit Akıncı, and Hamza Polat † Department of Physics, Dokuz Eyl¨ul University, TR-35160 Izmir, Turkey (Dated: October 10, 2018)A spin-1 Blume-Capel model with dilute and random crystal fields is examined for honeycomb andsquare lattices by introducing an effective-field approximation that takes into account the correla-tions between different spins that emerge when expanding the identities. For dilute crystal fields, wehave given a detailed exploration of the global phase diagrams of the system in k B T c /J − D/J planewith the second and first order transitions, as well as tricritical points. We have also investigatedthe effect of the random crystal field distribution characterized by two crystal field parameters
D/J and △ /J on the phase diagrams of the system. The system exhibits clear distinctions in qualitativemanner with coordination number q for random crystal fields with △ /J, D/J = 0. We have alsofound that, under certain conditions, the system may exhibit a number of interesting and unusualphenomena, such as reentrant behavior of first and second order, as well as a double reentrance withthree successive phase transitions. PACS numbers
Contents
I. Introduction II. Formulation III. Results and Discussion
IV. Conclusions Acknowledgements A. Fundamental correlation functions of the system for a square lattice B. The complete set of twenty one linear equations of a honeycomb lattice References I. INTRODUCTION
Spin-1 Blume-Capel (BC) model is one of the most extensively studied models in statistical mechanics and con-densed matter physics. The model exhibits a variety of multicritical phenomena such as a phase diagram with orderedferromagnetic and disordered paramagnetic phases separated by a transition line that changes from a continuous phasetransition to a first-order transition at a tricritical point. On the other hand, as an extension of the model, BC modelwith a random crystal field represents the critical behavior of He − He mixtures in a random media, i.e., aerogelwhere S = 0 and S = ± He and He atoms, respectively . From the theoretical point of view,BC model with a random crystal field (RCF) has been studied by a variety of techniques such as cluster variationalmethod (CVM) , Bethe lattice approximation (BLA) , effective field theory (EFT) , finite cluster approximation(FCA) , mean field theory (MFT) , Monte Carlo (MC) simulations , pair approximation (PA) , and renor-malization group (RG) method . Among these studies, EFT and MFT have been widely used to investigate thethermal and magnetic properties of BC model with a RCF distribution. For example, Kaneyoshi and Mielnicki inves-tigated the phase diagram of the system for a honeycomb lattice by using EFT with correlations and they found someimportant differences from the results obtained by the standard MFT. Similarly, in a recent paper, Yan and Deng considered the same model within the framework of EFT, and they derived the expressions of magnetizations forhoneycomb and square lattices. On the other hand, in several studies based on MFT, the authors paid attentionto the effects of crystal field dilution on the phase diagrams of the system and they observed that the system mayexhibit a reentrant behavior, as well as first order phase transitions.However, EFT and MFT studies mentioned above have some unsatisfactory results. Namely, the results obtainedby EFT are limited to second-order phase transitions and tricritical points, and a detailed description of first-ordertransitions has not been reported. Other than this, it is well known that magnetic systems with dilute crystal fieldsexhibit qualitatively similar characteristics when compared to site dilution problem of magnetic atoms. From thispoint of view, for a BC model with diluted crystal fields, MFT predicts that the phase transition temperature of thesystem will remain at a finite value until zero concentration is reached. In this context, we believe that BC modelwith RCF still deserves particular attention for investigating the proper phase diagrams, especially the first-ordertransition lines that include reentrant phase transition regions. Conventional EFT approximations include spin-spincorrelations resulting from the usage of the Van der Waerden identities, and provide results that are superior to thoseobtained within the traditional MFT. However, these conventional EFT approximations are not sufficient enough toimprove the results, due to the usage of a decoupling approximation (DA) that neglects the correlations betweendifferent spins that emerge when expanding the identities. Therefore, taking these correlations into consideration willimprove the results of conventional EFT approximations. In order to overcome this point, recently we proposed anapproximation that takes into account the correlations between different spins in the cluster of a considered lattice .Namely, an advantage of the approximation method proposed by this study is that no decoupling procedure is usedfor the higher-order correlation functions.In this paper, we intent to investigate the effects of RCF distributions on the phase diagrams of spin-1 BC modelon 2D lattices, namely honeycomb ( q = 3) and square ( q = 4) lattices. For this purpose, the paper is organized asfollows: In Sec. II we briefly present the formulations. The results and discussions are presented in Sec. III, andfinally Sec. IV contains our conclusions. II. FORMULATION
In this section, we give the formulation of the present study for a 2D lattice which has N identical spins arranged.We define a cluster on the lattice which consists of a central spin labeled S , and q perimeter spins being the nearestneighbors of the central spin. The cluster consists of ( q + 1) spins being independent from the spin operator ˆ S . Thenearest-neighbor spins are in an effective field produced by the outer spins, which can be determined by the conditionthat the thermal average of the central spin is equal to that of its nearest-neighbor spins. The Hamiltonian describingour model is H = − J X h i,j i S zi S zj − X i D i ( S zi ) , (1)where the first term is a summation over the nearest-neighbor spins with S zi = ± , D i on the secondsummation represents a random crystal field, distributed according to a given probability distribution. In this paper,we primarily deal with two kinds of probability distributions, namely, a quenched diluted crystal field distributionand a double peaked delta distribution which are given by Eqs. (2) and (3), respectively as follows P ( D i ) = pδ ( D i − D ) + (1 − p ) δ ( D i ) , (2) P ( D i ) = 12 { δ [ D i − ( D − △ )] + δ [ D i − ( D + △ )] } . (3)where p denotes the concentration of the spins on the lattice which are influenced by a crystal field D .We can construct the mathematical background of our model by using the approximated spin correlation identities by taking into account random configurational averages hh{ f i } S zi ii r = (cid:28)(cid:28) { f i } Tr i ( S zi ) exp ( − βH i )Tr i exp ( − βH i ) (cid:29)(cid:29) r , (4) hh{ f i } ( S zi ) ii r = (cid:28)(cid:28) { f i } Tr i ( S zi ) exp ( − βH i )Tr i exp ( − βH i ) (cid:29)(cid:29) r , (5)where β = 1 /k B T , { f i } is an arbitrary function which is independent of the spin variable S i and the inner h ... i andthe outer h ... i r brackets represents the thermal and random configurational averages, respectively.In order to apply the differential operator technique , we should separate the Hamiltonian (1) into two parts as H = H i + H ′ . Here, the effective Hamiltonian H i includes all the contributions associated with the site i , and theother part H ′ does not depend on the site i . − H i = ES zi + D i ( S zi ) , (6)where E = J P j S zj is the local field on the site i . If we use the matrix representations of the operators S zi and ( S zi ) for the spin-1 system then we can obtain the matrix form of Eq. (6) − H i = E + D − E + D . (7)Hereafter, we apply the differential operator technique in Eqs. (4) and (5) with { f i } = 1. From Eq. (4) we obtainthe following spin identity with thermal and configurational averages of a central spin for a lattice with a coordinationnumber q as hh S z ii r = ** q Y j =1 (cid:2) S zj sinh( J ∇ ) + ( S zj ) { cosh( J ∇ ) − } (cid:3)++ r F ( x ) | x =0 . (8)The function F ( x ) in Eq. (8) is defined by F ( x ) = Z dD i P ( D i ) f ( x, D i ) , (9)where f ( x, D i ) = 1 P n =1 exp( βλ n ) X n =1 h ϕ n | S zi | ϕ n i exp( βλ n ) , (10)= 2 sinh( βx )2 cosh( βx ) + e − βD i . In Eq. (10), λ n denotes the eigenvalues of − H i matrix in Eq. (7), and ϕ n represents the eigenvectors correspondingto the eigenvalues λ n of − H i matrix. With the help of Eq. (10), and by using the distribution functions defined inEqs. (2) and (3), the function F ( x ) in Eq. (9) can be easily calculated by numerical integration. Hereafter, we willfocus our attention on the construction of the correlation functions, as well as magnetization and quadrupole momentidentities of a honeycomb lattice with q = 3. A brief formulation of the fundamental spin identities for a square latticewith q = 4 can be found in Appendix A.By expanding the right-hand side of Eq. (8) for a honeycomb lattice with q = 3, we get the longitudinal magneti-zation as m z = hh S z ii r = l + 3 k hh S ii r + 3( l − l ) hh S ii r + 3 l hh S S ii r +6( k − k ) hh S S ii r + 3( l − l + l ) hh S S ii r + k hh S S S ii r + 3( l − l ) hh S S S ii r +3( k − k + k ) hh S S S ii r +( − l + 3 l − l + l ) hh S S S ii r , (11)We note that, for the sake of simplicity, the superscript z is omitted from the correlation functions on the right-handside of Eq. (11). The coefficients in Eq. (11) are defined as follows: l = F (0) ,l = cosh( J ∇ ) F ( x ) | x =0 , k = sinh( J ∇ ) F ( x ) | x =0 ,l = sinh ( J ∇ ) F ( x ) | x =0 , k = cosh( J ∇ )sinh( J ∇ ) F ( x ) | x =0 ,l = cosh ( J ∇ ) F ( x ) | x =0 , k = sinh ( J ∇ ) F ( x ) | x =0 ,l = cosh( J ∇ )sinh ( J ∇ ) F ( x ) | x =0 , k = cosh ( J ∇ )sinh( J ∇ ) F ( x ) | x =0 ,l = cosh ( J ∇ ) F ( x ) | x =0 , (12)Next, the average value of the perimeter spin in the system can be written as follows, and it is found as m = hh S zδ ii r = hh S z sinh( J ∇ ) + ( S z ) { cosh( J ∇ ) − }ii r F ( x + γ ) | x =0 , (13) hh S ii r = a (cid:0) − hh ( S z ) ii r (cid:1) + a hh S z ii r + a hh ( S z ) ii r , (14)with the coefficients a = F ( γ ) ,a = sinh( J ∇ ) F ( x + γ ) | x =0 ,a = cosh( J ∇ ) F ( x + γ ) | x =0 , (15)where γ = ( q − A is the effective field produced by the ( q −
1) spins outside the system and A is an unknownparameter to be determined self-consistently. In the effective-field approximation, the number of independent spinvariables describes the considered system. This number is given by the relation ν = hh ( S zi ) S ii r . As an example forthe spin-1 system, 2 S = 2 which means that we have to introduce the additional parameters hh ( S z ) ii r and hh ( S zδ ) ii r resulting from the usage of the Van der Waerden identity for the spin-1 Ising system. With the help of Eq. (5),quadrupolar moment of the central spin can be obtained as follows hh ( S z ) ii r = ** q Y j =1 (cid:2) S zj sinh( J ∇ ) + ( S zj ) { cosh( J ∇ ) − } (cid:3)++ r G ( x ) | x =0 , (16)where the function G ( x ) is defined as G ( x ) = Z dD i P ( D i ) g ( x, D i ) . (17)Definition of the function g ( x, D i ) in Eq. (17) is given as follows and the expression in Eq. (18) can be evaluated byusing the eigenvalues and corresponding eigenvectors of the effective Hamiltonian matrix in Eq. (7). g ( x, D i ) = 1 P n =1 exp( βλ n ) X n =1 h ϕ n | ( S zi ) | ϕ n i exp( βλ n ) , (18)= 2 cosh( βx )2 cosh( βx ) + e − βD i . Hence, we get the quadrupolar moment by expanding the right-hand side of Eq. (16) hh ( S z ) ii r = r + 3 n hh S ii r + 3( r − r ) hh S ii r + 3 r hh S S ii r +6( n − n ) hh S S ii r + 3( r − r + r ) hh S S ii r + n hh S S S ii r +3( r − r ) hh S S S ii r + 3( n − n + n ) hh S S S ii r +( − r + 3 r − r + r ) hh S S S ii r , (19)with r = G (0) ,r = cosh( J ∇ ) G ( x ) | x =0 , n = sinh( J ∇ ) G ( x ) | x =0 ,r = sinh ( J ∇ ) G ( x ) | x =0 , n = cosh( J ∇ )sinh( J ∇ ) G ( x ) | x =0 ,r = cosh ( J ∇ ) G ( x ) | x =0 , n = sinh ( J ∇ ) G ( x ) | x =0 ,r = cosh( J ∇ )sinh ( J ∇ ) G ( x ) | x =0 , n = cosh ( J ∇ )sinh( J ∇ ) G ( x ) | x =0 ,r = cosh ( J ∇ ) G ( x ) | x =0 , (20)Corresponding to Eq. (13), hh ( S zδ ) ii r = hh S z sinh( J ∇ ) + ( S z ) { cosh( J ∇ ) − }ii r G ( x + γ ) , (21) hh S ii r = b (cid:0) − hh ( S z ) ii r (cid:1) + b hh S z ii r + b hh ( S z ) ii r . (22)where b = G ( γ ) ,b = sinh( J ∇ ) G ( x + γ ) | x =0 ,b = cosh( J ∇ ) i G ( x + γ ) | x =0 . (23)Eqs. (11), (14), (19) and (22) are the fundamental spin identities of the system. When the right-hand sides of Eqs.(8) and (16) are expanded, the multispin correlation functions appear. The simplest approximation, and one of themost frequently adopted is to decouple these correlations according to (cid:10)(cid:10) S zi ( S zj ) ...S zl (cid:11)(cid:11) r ∼ = hh S zi ii r (cid:10)(cid:10) ( S zj ) (cid:11)(cid:11) r ... hh S zl ii r , (24)for i = j = ... = l . The main difference of the method used in this study from the other approximations inthe literature emerges in comparison with any decoupling approximation (DA) when expanding the right-hand sidesof Eqs. (8) and (16). In other words, one advantage of the approximation method used in this study is that nouncontrolled decoupling procedure is used for the higher-order correlation functions.For spin-1 Ising system with q = 3, taking Eqs. (11), (14), (19) and (22) as a basis, we derive a set of linear equationsof the spin identities. At this point, we assume that (i) the correlations depend only on the distance between thespins, (ii) the average values of a central spin and its nearest-neighbor spin (it is labeled as the perimeter spin) areequal to each other with the fact that, in the matrix representations of spin operator ˆ S , the spin-1 system has theproperties ( S zδ ) = S zδ and ( S zδ ) = ( S zδ ) . Thus, the number of the set of linear equations obtained for the spin-1Ising system with q = 3 reduces to twenty one, and the complete set is given in Appendix B.If Eq. (B1) is written in the form of a 21 ×
21 matrix and solved in terms of the variables x i [( i = 1 , , ..., e.g., x = hh S z ii r , x = hh S S ii r , ... )] of the linear equations, all of the spin correlation functions, as well as magnetizations andquadrupolar moments can be easily determined as functions of the temperature and Hamiltonian parameters. Sincethe thermal and configurational average of the central spin is equal to that of its nearest-neighbor spins within thepresent method then the unknown parameter A can be numerically determined by the relation hh S z ii r = hh S ii r or x = x . (25)By solving Eq. (25) numerically at a given fixed set of Hamiltonian parameters we obtain the parameter A . Thenwe use the numerical values of A to obtain the spin correlation functions which can be found from Eq. (B1). Notethat A = 0 is always the root of Eq. (25) corresponding to the disordered state of the system. The nonzero root of A in Eq. (25) corresponds to the long-range ordered state of the system. Once the spin identities have been evaluatedthen we can give the numerical results for the thermal and magnetic properties of the system. Since the effectivefield γ is very small in the vicinity of k B T c /J , we can obtain the critical temperature for the fixed set of Hamiltonianparameters by solving Eq. (25) in the limit of γ → III. RESULTS AND DISCUSSION
In this section, we will discuss the effect of the crystal field distributions defined in Eqs. (2) and (3) on the globalphase diagrams of the system where the second and first order transitions are shown by solid and dashed curves,respectively with tricritical points (shown by hollow circles) for honeycomb ( q = 3) and square ( q = 4) lattices. Also,in order to clarify the type of the transitions in the system, we will give the temperature dependence of the orderparameter. A. Phase diagrams of the system with dilute crystal field
In this section, we illustrate the phase diagrams and magnetization curves of the system with a dilute crystal fielddistribution defined in Eq.(2) where crystal field D is turned on, or turned off with probabilities p and (1 − p ) on thelattice sites, respectively. In Figs. (1a) and (1c), we plot the phase diagrams of the system in ( k B T c /J − D/J ) planefor honeycomb and square lattices with coordination numbers q = 3 and q = 4, respectively. As seen in Figs. (1a)and (1c), phase diagrams of the system can be divided into three parts with different concentration values p . Forthe curves in the first group with p < p ∗ , the system always exhibits a second order phase transition with a finitecritical temperature k B T c /J which extent to D/J → −∞ . If the concentration p reaches its critical value p ∗ then thecritical temperature depresses to zero. Physical reason underlying this behavior can be explained as follows; when weselect sufficiently large negative crystal field values (i . e .D/J → −∞ ), all of the spins in the system tend to align in S = 0 state. As p increases starting from zero, the ratio of spins which aligned in S = 0 state increases, and therefore,magnetization weakens, and accordingly, critical temperature of the system decreases. According to our numericalresults, the critical concentration value is obtained as p ∗ = 0 . q = 3 and p ∗ = 0 . q = 4. In the secondgroup of the phase diagrams in Figs. (1a) and (1c), the system exhibits a reentrant behavior of second order, whereasthe curves in the third group, exhibit a reentrant behavior of first order with a tricritical point at which a first ordertransition line turns into a second order transition line. Besides, the curves which exhibit a reentrant behavior offirst (or second) order, depress to zero at three successive values of crystal field D/J = − . , − . , − .
0. Moreover,in
D/J → ∞ limit, the system behaves like spin-1 / p = 1 .
0. In the case of p = 0, the ratio of spins that behavelike S = ± p increases. Therefore, for 0 ≤ p ≤ .
0, all transition lines have finite critical temperatureswhich increase with increasing p values for D/J → ∞ . At this point, we also note that if we select
D/J = 0 in Eq.(2), all lattice sites expose to a crystal field D i /J = 0 independent from p . Hence, all transition lines intersect eachother on the point ( D/J, k B T c /J ) = (0 , . q = 3, and ( D/J, k B T c /J ) = (0 , . q = 4. Meanwhile,previous studies based on EFT are not capable of obtaining first order transition lines of the system. From this pointof view, we see that our method improves the results of the other EFT works and we take the conventional EFTmethod one step forward by investigating the global phase diagrams, especially the first-order transition lines thatinclude reentrant phase transition regions. TABLE I: Critical concentration p ∗ obtained by several methodsand the present work for honeycomb ( q = 3) and square ( q = 4)lattices.Lattice EFT-I EFT-II MFT PA Present Work q = 3 0.484 0.492 1.0 0.5 0.3795 q = 4 0.604 0.610 1.0 0.667 0.5875 On the other hand, Figs. (1b) and (1d) shows the phase boundary in ( k B T c /J − p ) plane which separates theferromagnetic and paramagnetic phases with D/J → −∞ . According to this figure, critical temperature k B T c /J ofsystem decreases gradually, and ferromagnetic region gets narrower as p increases, and k B T c /J value depresses tozero at p = p ∗ . Such a behavior is an expected fact in dilution problems. Numerical value of critical concentration p ∗ for honeycomb ( q = 3) and square ( q = 4) lattices is given in Table I, and compared with the other works in theliterature. As seen in Table I, numerical values of p c for q = 3 and q = 4 are new results in literature. Furthermore,MFT predicts that the system always has a finite critical temperature and exists in a ferromagnetic state at lowertemperatures in D/J → −∞ limit, except that p = 1 .
0. This artificial result can be regarded as a failure of the MFT.In Fig.(2), we plot the temperature dependencies of magnetization curves corresponding to the phase diagramsdepicted in Fig. (1) for q = 3. As seen in Fig. (2), as p increases then critical temperature k B T c /J values decrease for D/J <
0, except the reentrant phase transition temperatures which occur at low temperatures. On the other hand,effect of increasing p values on the shape of magnetization curves depends on value of D/J . Namely, in Figs. (2a) and(2b) we see that ground state saturation values of magnetization curves decreases as p increases for D/J = − . − .
1. Moreover, for
D/J = − .
1, magnetization curves of the system exhibit a broad maximum at low temperaturesfor p = 0 .
37, and a reentrant behavior of second order for p = 0 .
4. If we select
D/J = − . p = 0 , . , . p > .
3. If p increasesfurther, a reentrant behavior of second order appears for p = 0 .
53, and we see a broad maximum at low temperaturesfor p = 0 .
517 and 0 . p decreases. This broad maximum behavior of magnetization curvesoriginates from the increase in the number of spins directed parallel to the z-direction with increasing temperature, dueto the thermal agitation. For D/J = − . m = 1 at the ground stateand reentrant behavior disappears. If we increase D/J further, for example for
D/J = − . m = 1 and thesystem always undergoes a second order phase transition from a ferromagnetic phase to a paramagnetic phase withincreasing temperature, which can be seen in Fig. (2f) with D/J = 10 .
0. As a common property of the curves inFig. (2), we see that effect of p on the saturation values, as well as temperature dependence of magnetization curvesstrictly depend on the strength of D/J . Hence, according to us, the presence of dilute crystal fields on the systemshould produce a competition effect on the phase diagrams of the system. We also note that, although it has not beenshown in the present work, magnetization curves for q = 4 corresponding to the phase diagrams depicted in Fig. (1c)exhibit qualitatively similar behavior with those of Fig. (2) with q = 3.As seen in Fig. (1), for a dilute crystal field distribution defined in Eq. (2), the global phase diagrams which areplotted in ( k B T c /J − D/J ) plane, as well as the phase boundaries in ( k B T c /J − p ) plane for q = 3 exhibit qualitativelysimilar characteristics when compared with those for q = 4. Hence, in order to examine the phase diagrams which areplotted in ( k B T c /J − D/J ) plane in Figs. (1a) and (1c) in detail, we plot the evolution of the global phase diagramsin Fig. (3) only for q = 4. From this point of view, Fig. (3a) shows how the phase diagrams in Fig. (1c) evolvewhen the concentration p changes from 0.5 to 0.6. As seen in Fig. (3a), we observe a second order phase transitionline with a finite critical temperature k B T c /J which extent to D/J → −∞ for p = 0 . p increases, namelyfor p = 0 .
580 and 0 . − . < D/J < − .
0, as wellas a high temperature phase boundary which extents to
D/J → −∞ . If p increases further, such as for p = 0 . .
585 and 0 . − . < D/J < − .
0, and the phase diagrams exhibit a bulge on the right hand side of ( k B T c /J − D/J ) plane, whereasanother transition line emerges within the range of −∞ < D/J < − .
0, which disappears for p > . p changes from 0.6 to 0.7 can be seen in Fig.(3b). As seen in this figure, the curves for p = 0 .
60, 0 .
62, 0 .
64, 0 .
66 exhibit a reentrant behavior of second order,while for p = 0 .
68 reentrance disappears and for p = 0 .
70 and 0 .
71 double reentrance with three successive secondorder phase transitions occurs in a very narrow region of
D/J . On the other hand, increasing values of p generatesfirst order phase transitions with tricritical points, as well as reentrant behavior of first order. This phenomena isillustrated in Fig. (3c). From Fig. (3c), we see that, the second order transition temperatures decrease as absolutevalue of D/J increases, and turn into first order transition lines at tricritical points. Evidently, the phase diagramschange abruptly for p ≥ . p = 0 . p = 0 . p ≤ . p ≥ . D/J = − . D/J = − .
0, respectively. In order to investigate the phase transition features of the systemfurther, we should continue increasing the value of p . In Fig. (3d), we see that the curves for p = 0 .
72, 0 .
74, 0 . .
78 exhibit a reentrant behavior of first order, whereas the curves with p = 0 .
80, 0 .
82, and 0 .
84 exhibit doublereentrance with two first order and a second order transition temperature.
B. Phase diagrams of the system with random crystal field
Next, in order to investigate the effect of the random crystal fields defined in Eq.(3) on the thermal and magneticproperties of the system, we represent the phase diagrams and corresponding magnetization curves for honeycomb( q = 3) and square lattices ( q = 4) throughout Figs. (4) and (8).We note that random crystal field distribution given in Eq. (3) with D/J = 0 corresponds to a bimodal distributionfunction, while for △ /J = 0, we obtain pure BC model with homogenous crystal field D/J . In Fig. (4), phase diagramsof the system corresponding to the bimodal distribution function are shown in ( k B T c /J − △ /J ) plane. For a bimodaldistribution, the phase diagrams have symmetric shape with respect to △ /J which comes from the fact that p = 1 / q . Namely, for q = 3, transition temperature decreaseswith increasing △ /J and exhibits double reentrance with three second order phase transition temperatures, thenfalls to zero at △ /J = 3 . △ /J increases then the transition temperature of the system for q = 4 decreases and remains at a finite value for △ /J → ∞ which means that ferromagnetic exchange interactions for q = 3 are insufficient for the system to keep itsferromagnetic order for △ /J > .
0, while for q = 4 these interactions are dominant in the system, and the presenceof a disorder in the crystal fields cannot destruct the ferromagnetic order.At the same time, in order to see the effect of the random crystal fields with △ /J, D/J = 0 on the phase diagramsand magnetization curves of the system for q = 3 and 4, we plot the phase diagrams in ( k B T c /J − D/J ) plane inFig. (5) and variation of the corresponding magnetization curves with temperature in Figs. (6) and (7), respectively.At first sight, it is obvious that the phase diagrams in Fig.(5) represent evident differences in qualitative mannerwith coordination number q . That is, as seen in Fig. (5a), the curve corresponding to △ /J = 0 represents thephase diagram of pure BC model for a honeycomb lattice which exhibits a reentrant behavior of first order with firstand second order transition lines, as well as a tricritical point. From Fig. (5a), we see that as △ /J increases thenthe tricritical point and first order transitions disappear, and the first order reentrance turns into double reentrancewith three transition temperatures of second order, and phase transition lines shift to positive crystal field directionwithout changing their shapes. On the other hand, the situation is very different for a square lattice. Namely, asseen in Figs. (5a) and (5b), △ /J = 0 curves for q = 3 and q = 4 are qualitatively identical to each other. However,as seen in Fig. (5b), for △ /J = 0, first order transition lines and tricritical points do not disappear from the systemfor q = 4, but shift to negative crystal field values. Besides, the system does not exhibit double reentrance for q = 4.Furthermore, for △ /J ≥ . △ /J ≥ D/J varies. Namely, critical temperature k B T c /J in Fig. (5a) reduces to zero at D/J = △ /J − .
0. On the otherhand, first order transition temperatures in Fig. (5b), reduce to zero at
D/J = −△ /J .It is important to note that these observations are consistent with the results shown in Figs. (1a) and (1c). Inother words, the distribution function given in Eq. (2) with p = 0 . D/J = 2 D /J is identical to Eq. (3) for △ /J = ± D /J and D/J = D /J . For example, according to Eq. (2), if we select D /J = 4 . p = 0 .
5, it meansthat half of the spins on the lattice sites expose to a crystal field
D/J = 0, while a crystal field given by
D/J = 8 . △ /J = ± D /J and D/J = D /J by using Eq. (3),we generate the same distribution again. Hence, we expect to get the same results in Figs. (1) and (5) for these systemparameters. For instance, for D /J = 4 . D/J = 8 .
0, and the system exhibits a ferromagneticorder in the ground state, which can also be seen in Fig. (5a) with a critical temperature k B T c /J = 1 . q = 4, and for the whole temperature region on the phase diagrams. Therefore, the state(para-or ferro), as well as thermal and magnetic properties of a selected ( k B T /J, D/J ) point with respect to p = 0 . k B T /J, D/ J ) in Figs. (5a) and (5b) with respect tothe curve △ /J = ± D/J , respectively. Moreover, the qualitative differences between Figs. (5a) and (5b) mentionedabove are strongly related to the percolation threshold value of the lattice. Namely, distribution function Eq.(3) isvalid only for p = 0 .
5. However, as seen in Table I, we obtain p c < . q = 3, and p c > . q = 4.In Fig. (6), we examine the temperature dependence of magnetization curves for q = 3, corresponding to the phasediagrams shown in Fig. (5a) with D/J = − .
0. Fig. (6a), shows how the temperature dependence of magnetizationcurves evolve when △ /J changes. According to Fig. (6a), magnetization curves saturate at a partially ordered stateat low temperatures. Besides, for △ /J = 1 . , .
74 and 1 .
80, the system undergoes three successive phase transitionsof second order, which confirms the existence of double reentrance. Similarly, Fig. (6b) shows how the shape of themagnetization curves change as
D/J changes for constant △ /J = 6 .
0. As seen in Fig. (6b), magnetization curvesexhibit a second order phase transition from a ferromagnetic (fully ordered) to a paramagnetic phase at certain valuesof crystal field, namely at
D/J = 4 . , . .
4, whereas for
D/J = 3 . , . , . . D/J = 3 . , . . D/J = 3 .
3, we observe doublereentrance. Fig. (7) shows the magnetization curves for q = 4, corresponding to the phase diagrams shown in Fig.(5b). In Fig. (7a), we see that magnetization curves exhibit a second order phase transition from a paramagneticphase to a fully ordered ferromagnetic phase for △ /J = 0 and 5 .
0. On the other hand, the curves corresponding to
D/J = 5 .
5, 5 . .
0, saturate at a partially ordered state at low temperatures, and exhibit a broad maximum withincreasing temperature which depresses gradually as △ /J increases, then fall rapidly at a second order phase transitiontemperature. The broad maximum behavior observed in these curves disappears for △ /J = 8 .
0. Additionally, Fig.(7b) represents the magnetization versus temperature curves for q = 4 with with D/J = − .
0. In Fig. (7b), itis clearly evident that, at low temperatures, the system saturates at a partially ordered phase for △ /J = 4 .
0, 5 . . .
0, while for △ /J = 3 . .
8, a reentrant behavior of first order occurs. Again we see that, there is acompetition between ferromagnetic exchange interactions and disorder effects in crystal fields which determines thesaturation values and temperature dependence of magnetization curves of the system.Finally, dependence of magnetization of the system on the crystal field △ /J for fixed temperature values k B T /J =0 . , . , . . D/J = 2 . q = 3 and 4, respectively. We see that at sufficientlylow temperatures such as k B T /J = 0 .
01, magnetization curves exhibit three phases for q = 3. A first order transitionis characterized by a gap in this figure. On the left panel in Fig. (8), which is plotted for q = 3, we observe twosuccessive first order phase transitions for k B T /J = 0 .
01. The first one is from the fully ordered ferromagnetic phase( m = 1 .
0) to the partly ordered phase ( m = 0 . m = 0 . q = 4.Hence, we observe two phases. Namely, for k B T /J = 0 .
01, a first order phasetransition from a fully ordered phase ( m = 1 .
0) to a partly ordered phase ( m = 0 . m = 0 . q = 3 and 4. Since in Figs. (5a) and (5b), all phase diagrams exhibit similarbehavior as D/J varies, behavior of magnetization curves in Fig. (8) should be the same for different
D/J values.Namely, the left panel of Fig. (8) can be regarded as the variation of magnetization with △ /J within the range D/J + 1 . < △ /J < D/J + 4 . q = 3. for a selected value of D/J . It is possible to observe the similar behavioron the right panel in Fig. (8). This general behavior is an expected result, since the behavior of the system dependson the relation between △ /J and D/J parameters for a given temperature, not their values independently.
IV. CONCLUSIONS
In this work, we have studied the phase diagrams of a spin-1 Blume-Capel model with diluted and random crystalfield interactions on two dimensional lattices. We have introduced an effective-field approximation that takes intoaccount the correlations between different spins in the cluster of a considered lattice and examined the phase diagramsas well as magnetization curves of the system for different types of crystal field distributions, namely, dilute crystalfields and a double peaked delta distribution, given by Eqs. (2) and (3), respectively.For dilute crystal fields, we have given a detailed exploration of the global phase diagrams of the system in k B T c /J − D/J plane with the second and first order transitions, as well as tricritical points. We have also shown that the systemwith dilute crystal fields exhibits a percolation threshold value p c which can not be predicted by standard MFA. Inaddition, we have observed multi-reentrant phase transitions for specific set of system parameters.On the other hand, we have investigated the effect of the random crystal field distribution characterized by twocrystal field parameters D/J and △ /J on the phase diagrams of the system. As a limited case, we have alsofocused on a bimodal distribution with D/J = 0. Particulary, we have reported the following observations for abimodal distribution: It has been found that the phase diagrams have symmetric shape with respect to △ /J whichcomes from the fact that p = 1 /
2. The transition temperatures are of second order, and the system exhibit differentcharacteristic features depending on the coordination number q . Besides, we have realized that the system may exhibitclear distinctions in qualitative manner with coordination number q for random crystal fields with △ /J, D/J = 0.Moreover, we have discussed a competition effect which arises from the presence of dilution, as well as random crystalfields, and we have observed that saturation values of the magnetization curves are strongly related to these effects.As a result, we can conclude that all of the points mentioned above show that our method improves the conventionalEFT methods based on decoupling approximation. Therefore, we hope that the results obtained in this work may bebeneficial from both theoretical and experimental points of view. Acknowledgements
One of the authors (Y.Y.) would like to thank the Scientific and Technological Research Council of Turkey(T ¨UB˙ITAK) for partial financial support. This work has been completed at Dokuz Eyl¨ul University, GraduateSchool of Natural and Applied Sciences. Partial financial support from SRF (Scientific Research Fund) of DokuzEyl¨ul University (2009.KB.FEN.077) (H.P.) is also acknowledged.0
Appendix A: Fundamental correlation functions of the system for a square lattice
Magnetization of the central spin for a square lattice is given as follows hh S z ii = µ + 4 c hh S ii r + 4( µ − µ ) hh S ii r +6 µ hh S S ii r + 12( c − c ) hh S S ii r +6( µ − µ + µ ) hh S S ii r + 4 c hh S S S ii r +12( µ − µ ) hh S S S ii r + 12( c − c + c ) hh S S S ii r +4( µ − µ + 3 µ − µ ) hh S S S ii r + µ hh S S S S ii r +4( c − c ) hh S S S S ii r + 6( µ − µ + µ ) hh S S S S ii r +4( c − c + 3 c − c ) hh S S S S ii r +( µ − µ + 6 µ − µ + µ ) hh S S S S ii r , (A1)where the coefficients are given by µ = F (0) ,µ = sinh ( J ∇ ) F ( x ) | x =0 , c = sinh( J ∇ ) F ( x ) x =0 ,µ = cosh( J ∇ ) F ( x ) | x =0 , c = sinh( J ∇ ) cosh( J ∇ ) F ( x ) x =0 ,µ = cosh ( J ∇ ) F ( x ) | x =0 , c = sinh ( J ∇ ) F ( x ) x =0 ,µ = sinh ( J ∇ ) cosh( J ∇ ) F ( x ) | x =0 , c = cosh ( J ∇ ) sinh( J ∇ ) F ( x ) x =0 ,µ = cosh ( J ∇ ) F ( x ) | x =0 , c = sinh ( J ∇ ) cosh( J ∇ ) F ( x ) x =0 ,µ = sinh ( J ∇ ) cosh ( J ∇ ) F ( x ) x =0 , c = cosh ( J ∇ ) sinh( J ∇ ) F ( x ) x =0 ,µ = cosh ( J ∇ ) F ( x ) | x =0 ,µ = sinh ( J ∇ ) F ( x ) x =0 . Quadrupolar moment corresponding to equation (B1) defined as (cid:10) h ( S z ) (cid:11) i = ρ + 4 η hh S ii r + 4( ρ − ρ ) hh S ii r +6 ρ hh S S ii r + 12( η − η ) hh S S ii r +6( ρ − ρ + ρ ) hh S S ii r + 4 η hh S S S ii r +12( ρ − ρ ) hh S S S ii r + 12( η − η + η ) hh S S S ii r +4( ρ − ρ + 3 ρ − ρ ) hh S S S ii r + ρ hh S S S S ii r +4( η − η ) hh S S S S ii r + 6( ρ − ρ + ρ ) hh S S S S ii r +4( η − η + 3 η − η ) hh S S S S ii r +( ρ − ρ + 6 ρ − ρ + ρ ) hh S S S S ii r , (A2)where ρ = G (0) ,ρ = sinh ( J ∇ ) G ( x ) | x =0 , η = sinh( J ∇ ) G ( x ) x =0 ,ρ = cosh( J ∇ ) G ( x ) | x =0 , η = sinh( J ∇ ) cosh( J ∇ ) G ( x ) x =0 ,ρ = cosh ( J ∇ ) G ( x ) | x =0 , η = sinh ( J ∇ ) G ( x ) x =0 ,ρ = sinh ( J ∇ ) cosh( J ∇ ) G ( x ) | x =0 , η = cosh ( J ∇ ) sinh( J ∇ ) G ( x ) x =0 ,ρ = cosh ( J ∇ ) G ( x ) | x =0 , η = sinh ( J ∇ ) cosh( J ∇ ) G ( x ) x =0 ,ρ = sinh ( J ∇ ) cosh ( J ∇ ) G ( x ) x =0 , η = cosh ( J ∇ ) sinh( J ∇ ) G ( x ) x =0 ,ρ = cosh ( J ∇ ) G ( x ) | x =0 ,ρ = sinh ( J ∇ ) G ( x ) x =0 . Finally, perimeter spin identities are as follows hh S ii r = α (1 − hh ( S ) ii r ) + α hh S ii r + α hh ( S ) ii r , (A3) hh S ii r = ω + ω hh S ii r + ( ω − ω ) hh S ii r (A4)1with the coefficients α = F ( γ ) ω = G ( γ ) | x =0 α = sinh( J ∇ ) F ( x + γ ) ω = sinh( J ∇ ) G ( x + γ ) | x =0 α = cosh( J ∇ ) F ( x + γ ) ω = cosh( J ∇ ) G ( x + γ ) | x =0 where γ = ( q − A with q = 4, and the functions F ( x ) and G ( x ) are defined in equations (9) and (17). Appendix B: The complete set of twenty one linear equations of a honeycomb lattice hh S z ii r = l + 3 k hh S ii r + 3( l − l ) hh S ii r + 3 l hh S S ii r +6( k − k ) hh S S ii r + 3( l − l + l ) hh S S ii r + k hh S S S ii r + 3( l − l ) hh S S S ii r +3( k − k + k ) hh S S S ii r +( − l + 3 l − l + l ) hh S S S ii r hh S S ii r = (3 l − l ) hh S ii r + 3 k hh S ii r + 3( l − l + l + l ) hh S S ii r +6( k − k ) hh S S ii r + k hh S S S ii r +( − l + 3 l − l − l + 3 l + l ) hh S S S ii r +3( k − k + k ) hh S S S ii r hh S S S ii r = ( l − l + 3 l + 3 l ) hh S S ii r + (6 k − k ) hh S S ii r +( − l + 3 l − l − l + 3 l + l ) hh S S S ii r +(3 k − k + k + 3 k ) hh S S S ii r hh S ii r = a (1 − hh ( S z ) ii r ) + a hh S z ii r + a hh ( S z ) ii r hh S S ii r = a hh S ii r + a hh S S ii r + ( a − a ) hh S S ii r hh S S S ii r = a hh S S ii r + a hh S S S ii r + ( a − a ) hh S S S ii r hh S ii r = b (1 − hh ( S z ) ii r ) + b hh S z ii r + b hh ( S z ) ii r hh S S ii r = b hh S ii r + b hh S S ii r + ( b − b ) hh S S ii r hh S S ii r = b hh S ii r + b hh S S ii r + ( b − b ) hh S S ii r hh S S ii r = b hh S ii r + b hh S ii r hh S S S ii r = b hh S S ii r + b hh S S ii r hh S S S ii r = b hh S S ii r + b hh S S ii r hh S S S ii r = b hh S S ii r + b hh S S S ii r + ( b − b ) hh S S S ii r hh S S S ii r = b hh S S ii r + b hh S S S ii r + ( b − b ) hh S S S ii r hh S S S ii r = b hh S S ii r + b hh S S S ii r + ( b − b ) hh S S S ii r hh S ii r = r + 3 n hh S ii r + 3( r − r ) hh S ii r + 3 r hh S S ii r +6( n − n ) hh S S ii r + 3( r − r + r ) hh S S ii r + n hh S S S ii r +3( r − r ) hh S S S ii r + 3( n − n + n ) hh S S S ii r +( − r + 3 r − r + r ) hh S S S ii r hh S S ii r = (3 r − r ) hh S ii r + 3 n hh S ii r + (3 r + 3 r − r + 3 r ) hh S S ii r +6( n − n ) hh S S ii r + n hh S S S ii r +( − r + 3 r − r − r + 3 r + r ) hh S S S ii r +3( n − n + n ) hh S S S ii r hh S S ii r = (3 r − r ) hh S ii r + 3 n hh S ii r + (3 r + 3 r − r + 3 r ) hh S S ii r +6( n − n ) hh S S ii r + (3 n − n + n + 3 n ) hh S S S ii r +( − r + 3 r − r − r + 3 r + r ) hh S S S ii r hh S S S ii r = ( r − r + 3 r + 3 r ) hh S S ii r + ( − n + 6 n ) hh S S ii r − r + 3 r − r − r + 3 r + r ) hh S S S ii r +(3 n − n + n + 3 n ) hh S S S ii r hh S S S ii r = ( r − r + 3 r + 3 r ) hh S S ii r + ( − n + 6 n ) hh S S ii r ( − r + 3 r − r − r + 3 r + r ) hh S S S ii r +(3 n − n + n + 3 n ) hh S S S ii r hh S S S ii r = ( r − r + 3 r + 3 r ) hh S S ii r + ( − n + 6 n ) hh S S ii r +(3 n − n + n + 3 n ) hh S S S ii r +( − r + 3 r − r − r + 3 r + r ) hh S S S ii r . (B1)3 ∗ Also at Dokuz Eyl¨ul University, Graduate School of Natural and Applied Sciences, Turkey † Electronic address: [email protected] M. Blume, Phys. Rev. , 517 (1966). H. W. Capel, Physica , 966 (1966). A. Maritan, M. Cieplak, M. R. Swift, F. Toigo, J.R. Banavar, Phys. Rev. Lett. , 221 (1992). C. Buzano, A. Maritan, A. Pelizzola, J. Phys. Condens. Matter , 327 (1994). E. Albayrak, Physica A , 1529 (2011). T. Kaneyoshi, J. Phys. C , L557 (1986). T. Kaneyoshi, J. Phys. C , L679 (1988). T. Kaneyoshi and J. Mielnicki, J. Phys. Condens. Matter , 8773 (1990). T. Kaneyoshi, Phys. Status Solidi B , 313 (1992). S. L. Yan and L. L. Deng, Physica A , 301 (2002). A. Benyoussef and H. Ez-Zahraouy, J. Phys. Condens. Matter , 3411 (1994). V. Ilkovic, Phys. Status Solidi B , K7 (1995). A. Benyoussef, T. Biaz, M. Saber, and M. Touzani, J. Phys. C , 5349 (1987). M. E. S. Borelli and C. E. I Carneiro, Physica A , 249 (1996). C. E. I Carneiro, V. B. Henriques, and S. R. Salinas, J. Phys. Condens. Matter , 3687 (1989). N. Boccara, A. El Kenz, and M. Saber, J. Phys. Condens. Matter , 5721 (1989). C. E. I Carneiro, V. B. Henriques, and S. R. Salinas, J. Phys. A Math. Gen. , 3383 (1990). L. Bahmad, A. Benyoussef, and A. El Kenz, J. Magn. Magn. Mater , 397 (2008). I. Puha and H. T. Diep, J. Magn. Magn. Mater. , 85 (2001). D. P Lara and J. A. Plascak, Physica A , 443 (1998). N. S. Branco and B. M. Boechat, Phys. Rev. B , 11673 (1997). ¨U. Akıncı, Y. Y¨uksel, and H. Polat, Phys. Rev. E , 061103 (2011). F. C. S´aBarreto, I. P. Fittipaldi, B. Zeks, Ferroelectrics , 1103 (1981). R. Honmura and T. Kaneyoshi, J. Phys. C , 3979 (1979). T. Kaneyoshi, Acta Phys. Pol. A , 703 (1993). I. Tamura, T. Kaneyoshi, Prog. Theor. Phys. , 1892 (1981). -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 50.00.20.40.60.81.01.21.4 k B T c / J D/J p=00.10.20.3 0.40.50.6 0.70.80.91.0 q=3(a) k B T c / J p q=3ferromagnetic region paramagnetic region(b) p*D/J=-10.0 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 50.00.51.01.52.02.5 k B T c / J D/J p=00.10.20.30.40.5 0.6 0.7 0.80.91.0 q=4(c) (d) k B T c / J p q=4ferromagnetic region paramagnetic regionD/J=-10.0p* FIG. 1: (a) Phase diagrams of the system for q = 3 in a ( k B T c /J − D/J ) plane corresponding to dilute crystal field distributiondefined in Eq. (2). The solid and dashed lines correspond to second- and first-order phase transitions, respectively. The opencircles denote the tricritical points, and the numbers on each curve represent the value of concentration p . (b) Phase diagramsof the system for q = 3 in a ( k B T c /J − p ) plane with a selected value of the crystal field D/J = − .
0. (c) Phase diagrams ofthe system for q = 4 in a ( k B T c /J − D/J ) plane corresponding to dilute crystal field distribution defined in Eq. (2). The solidand dashed lines correspond to second- and first-order phase transition, respectively. The open circles refer to the tricriticalpoints, and the numbers on each curve represent the value of concentration p . (d) Phase diagrams of the system for q = 4 in a( k B T c /J − p ) plane with a selected value of the crystal field D/J = − . D/J=-10q=3 p=00.1
D/J=-3.1q=3 (a) (b)(c) (d)(e) (f)
D/J=-2.5q=3 m m D/J=-2.0q=3
D/J=-1.5q=3 k B T/Jk B T/Jk B T/Jk B T/J m k B T/J
D/J=10.0q=3 k B T/J p=1.0
FIG. 2: (Color online) Temperature dependence of magnetization corresponding to Fig. (1a) with some selected values ofcrystal field. (a)
D/J = − .
0, (b)
D/J = − .
1, (c)
D/J = − .
5, (d)
D/J = − .
0, (e)
D/J = − .
5, and (f)
D/J = 10 .
0. Thenumbers on each curve denote the value of concentration p . The solid and dashed lines correspond to second- and first-orderphase transitions, respectively. -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.00.00.10.20.30.40.50.60.70.8 f dcc bgf eb k B T c / J D/J p=0.575 (a) 0.580 (b) 0.583 (c) 0.584 (d) 0.585 (e) 0.587 (f) 0.588 (g) q=4(a) a -3.5 -3.0 -2.5 -2.0 -1.5 -1.00.00.20.40.60.81.01.21.41.61.8 p=0.71 k B T c / J D/J q=40.70.680.660.640.620.60(b) -3.0 -2.8 -2.6 -2.4 -2.2 -2.00.00.20.40.60.81.01.21.4 k B T c / J D/J p=0.719 0.717 0.71640.7160.712q=4(c) -2.6 -2.4 -2.2 -2.0 -1.8 -1.60.00.20.40.60.81.01.21.41.6 k B T c / J D/J p=0.720.74 q=4
FIG. 3: (Color online) Evolution of the phase diagrams corresponding to Fig. (1c). The numbers on each curve denote thevalue of concentration p . The solid and dashed lines correspond to second- and first-order phase transitions, respectively. Theopen circles indicate the tricritical points. k B T c / J /J D/J=0q=3 /J D/J=0q=4
FIG. 4: Phase diagrams of the system in a ( k B T c /J − △ /J ) plane for a bimodal crystal field distribution corresponding to Eq.(3) with D/J = 0 .
0. Left and right-hand side panels are plotted for q = 3 and q = 4, respectively. -2 0 2 4 6 8 100.00.20.40.60.81.01.21.41.6 (a) 8.06.04.03.0 k B T c / J D/J q=3/J=02.0 -10 -8 -6 -4 -2 0 2 4 6 8 100.00.51.01.52.02.5 (b) 8.06.04.03.0 k B T c / J D/J q=4/J=02.0
FIG. 5: (Color online) Phase diagrams of the system in a ( k B T c /J − D/J ) plane corresponding to random crystal fielddistribution defined in Eq. (3) for (a) q = 3, (b) q = 4. The numbers on each curve denote the value of △ /J . The opencircles represent the tricritical points, and the solid and dashed lines correspond to second- and first-order phase transitions,respectively. m k B T/J
D/J=-1.0q=3/J=1.61.641.661.68 1.74 1.80 (a) m k B T/J /J=6.0q=3D/J=4.4(b)
FIG. 6: (Color online) (a) Temperature dependence of magnetization curves corresponding to Fig. (5a) for q = 3 (a) with D/J = − . △ /J , and (b) with △ /J = 6 . D/J . m k B T/J
D/J=2.0q=4(a) /J=05.0 5.55.86.08.0 /J=3.8D/J=-4.0q=4 m k B T/J (b) 3.7 4.0 7.06.05.0
FIG. 7: (Color online) (a) Temperature dependence of magnetization curves for q = 4 corresponding to Fig. (5b) for (a) D/J = 2 . D/J = − . △ /J . ddcb (a) 0.01 (b) 0.05 (c) 0.1 (d) 0.2 mm /J D/J=2.0q=3a da (a) 0.01 (b) 0.05 (c) 0.1 (d) 0.2D/J=2.0q=4 /J b c FIG. 8: (Color online) Variation of the magnetization curves as a function of △ /J with D/J = 2 .
0. Left and right panels areplotted for q = 3 and 4, respectively. Each curve is plotted for different temperature values, namely k B T /J = 0 . , . , . ..