Monte Carlo investigation of the tricritical point stability in a three-dimensional Ising metamagnet 3 M.Zˇukoviˇc∗ and T.Idogaki 1 0 2 Department of Applied Quantum Physics, Graduate School of Engineering, n a J Kyushu University, Fukuoka 812-8581, Japan 9 ] h c e m Abstract. We use Monte Carlo simulations to study multicritical properties of an Ising - t a metamagnet in an external field. According to the mean field theory predictions, a three- t s . dimensional layered metamagnet is expected to display a tricritical point decomposition to t a m a critical endpoint and a bicritical endpoint, when a ratio between intralayer ferromagnetic - d and interlayer antiferromagnetic couplings becomes sufficiently small. Our simulations show n o no evidence of such a decomposition and produce a tricritical behaviour even for a coupling c [ ratio as small as R = 0.01. 1 v 6 1 8 1 . 1 PACS codes: 75.10.Hk; 75.30.Kz; 75.40.Cx; 75.40.Mg; 75.50.Ee. 0 3 Keywords: Ising metamagnet; Monte Carlo simulation; Multicritical behaviour; Tricritical 1 : v point decomposition. i X r a ∗Corresponding author. Permanent address: Department of Applied Quantum Physics, Graduate School of Engi- neering, Kyushu University, Fukuoka 812-8581, Japan Tel.: +81-92-642-3811; Fax: +81-92-633-6958 1 1.Introduction Ising metamagnets, systems with ferromagnetic and antiferromagnetic couplings simulta- neously present, have attracted much interest because it is possible to induce novel kinds of critical behaviour by forcing competition between these couplings, in particular by applying a magnetic field. They are generally believed to exhibit a tricritical point (TCP) which separates a second-order phase transition at high temperatures and low fields from a first- order phase transition at low temperatures and high fields (Fig.1a). Although this picture has in principle been confirmed experimentally1, some ambiguities concerning the tricritical behaviour still remain. In particular, the mean-field theory2 (MFT) predicts splitting of the TCP into a critical endpoint (CE) and a bicritical endpoint (BCE) (as shown in Fig.1b), if the ratio between total intrasublattice ferromagnetic and total intersublattice antiferromag- netic couplings is sufficiently small. Such a way, for example, if we increased the field at any temperature between T and T , the system would first undergo a first-order transition CE BCE from an antiferromagnetic phase to a generally different antiferromagnetic phase, and then a second-order transition to a paramagnetic phase. Since, as we know, the MFT neglects fluctuations, which could destroy this “middle” phase, some more sophisticated methods have been employed to give an answer to the question of whether the kind of phase dia- gram shown in Fig.1b can really exist or not. However, no definite conclusions have been drawn yet. So far, a two-dimensional next-nearest neighbor model (nnn−model) for Ising antiferromagnet has been studied by a Monte Carlo renormalization group3 (MCRG) with no indications of the TCP splitting, however, transfer matrix techniques4 did not reproduce (possibly because of the limited strip widths which could be studied) tricritical behaviour at very small ratios R. The same model has been also studied in three dimensions, which is more interesting case from the experimental point of view as well as from such a respect that fluctuations are smaller than in two dimensions and hence, the MFT predictions are more likely to be correct. Monte Carlo (MC) and MCRG methods applied to a variety of Hamiltonians for both meta−model and nnn−model of an Ising antiferromagnet showed that the exponent behaviour was consistent with the MFT, however, could not detect any change in the phase diagram itself5. However, more recent MC simulations on nnn−model produced only tricritical behaviour for R ≥ 0.056, i.e. the result contradicting to the MFT predictions. 2 Clear distinction in critical behaviour between two- and three-dimensional models was observed inthe antiferromagnetic Blume-Capel model, which, according to the MFT, should also show decomposition of the TCP. While MC simulations produced only tricritical be- haviour in two dimensions7, they provided clear evidence of the decomposition into a CE and a BCE in a three-dimensional model8. In light of the previous studies on Ising meta- magnets, which seemed to favor a non-decomposition of the TCP, quite perplexing results have been recently obtained by bothexperimental9-13 and Monte Carlo14-16 studies onone of the typical Ising metamagnets - FeBr . Although they failed to confirm the decomposition 2 of the TCP, they quite convincingly showed the existence of anomalies of the magnetization and the specific heat, which could be associated with the decomposition. These anomalies were attributed to the effectively weak ferromagnetic intralayer interaction and the high interlayer coordination - features present in the real compound. In this paper, we perform high precision MC simulations and provide convincing evidence of the stability of the tricritical point in the three-dimensional meta − model of an Ising metamagnet for R ≥ 0.01. Such a way, we show that merely weak, or even extremely weak, ferromagneticintralayer interactionshouldnotcausedecompositionoftheTCP, aspredicted by the MFT. Hence, in light of the previous results obtained for FeBr , the high interlayer 2 coordination, which is present in the real compound but absent in our model, seems to be of crucial importance in producing of such a decomposition, if it really takes place. 2.Model and simulation technique The considered system is the spin-1 Ising metamagnet described by the Hamiltonian 2 H = −J1 X sisj −J2 X sisk −hXsi , (1) <i,j> <i,k> i wheres = ±1isanIsingspin, < i,j >and< i,k >denotethesumovernearestneighborsin i theplaneandinadjacent planes, respectively, andhisanexternal magneticfield. Wechoose J > 0 and J < 0 so that each of the planes is ferromagnetic, but antiferromagnetically 1 2 coupled to adjacent planes. According to the MFT scenario, it is only the ratio R = z J /z |J |, where z and z are 1 1 2 2 1 2 numbers of nearest neighbors of the site i in the plane and adjacent planes, respectively, which determines the phase diagram. While for R > 3 it should look like the one in Fig.1a, 5 3 i.e. it should display a TCP with tricritical exponents α = 1,β = 1,γ = 1,λ = 1, for t 2 t 4 t t 2 0 < R < 3 the TCP is expected to split up into a CE and a BCE, as shown in Fig.1b. Then, 5 in the latter case, both the CE and the BCE should probably keep usual three-dimensional critical exponents: α ≈ 0.11,β ≈ 0.32,γ ≈ 1.24. At R = 3 the MFT predicts a four-order 5 critical point with yet different set of exponents. We have performed MC simulations on simple cubic lattice samples of linear sizes ranging from L = 16 to L= 40, assuming the periodic boundary condition throughout. We used an antiferromagnetic initial spin configuration at low temperatures and the field not exceeding the critical value h (p) = pz |J |, and a ferromagnetic one at high temperatures. As we c 2 2 moved in the temperature-field (T,h) space, we used the last spin configuration as an input for calculation at the next point. Spin updating followed a Metropolis dynamics. Averages were calculated using at most 25,000 Monte Carlo steps per spin (MCS/s) after equilibrating over another 5,000 to 10,000 MCS/s. Since we focused on the tricritical region, which for very small ratios R lies at very low temperatures, in order to prevent huge fluctuations at field-heating andfield-coolingprocesses (thepathofmeasurement wouldbevirtuallyparallel to the phase boundary), we only performed (h ↑)+(h ↓) loops, i.e. raised and lowered the field at fixed temperature and measured : the direct and staggered magnetizations m and m , respectively s m = hMi/N, where M = |Xsi +Xsj| , (2) i∈A j∈B ms = hMsi/N, where Ms = |Xsi −Xsj| , (3) i∈A j∈B where N is a total number of sites, and A, B denote sublattices made up of ”spin-up” and ”spin-down” planes, respectively, and the corresponding direct and staggered susceptibilities per site χ and χ , respectively s (hM2i−hMi2) χ = , (4) Nk T B (hM2i−hM i2) χ = s s . (5) s Nk T B These quantities were used to determine the nature as well as a location of the transition. The simulations were performed on the vector supercomputer FUJITSU VPP700/56. 3.Results and discussion 4 Westudied phase transitionsofthesystem forfourdifferent valuesofR: 0.5, 0.2, 0.05and 0.01, with particular attention paid to the region of the expected TCP. Fig.2 depicts phase boundaries in the vicinity of the TCPs for the respective values of R. At sufficiently low temperatures, both the direct and staggered magnetizations displayed besides jumps also pronounced hysteresis formed by the (h ↑)+(h ↓) cycles - features which are characteristic for a first-order transition. In Fig.2 the upper and lower branches in the low temperature region (blank circles) are determined by the jumps of the m curve in the (h ↑) and (h ↓) s processes, respectively, andoutlinetheregionofthecoexistence oftheantiferromagneticand ferromagnetic phases. Increasing temperature makes the hysteresis shrink, and eventually disappear at a point which gave us a first guess for the location of the TCP (arrows). Note, however, that if this is not done for sufficiently large L, the above mentioned phenomena, ac- companying first-order transitions, could be suppressed by the finite-size effects, and hence, lattices of small L are not suitable for such a kind of investigation. Upon further increase in temperature, the (h ↑) and (h ↓) magnetization curves collapse onto a single smooth curve, signifying a second-order transition (filled circles). Its location was determined by the lo- cation of the staggered susceptibility peak, extrapolated for L brought to infinity. As we can see, none of the cases shown in Fig.2 shows any indications of the TCP decomposition. On the other hand, according the MFT the splitting should take place in each of the cases and in very noticeable scale (e.g., for R = 0.05 the separation between the CE and BCE temperatures is estimated to 50%). Note that, in contrast to the results obtained for the three-dimensional nnn − model6, the varying value of R has virtually no influence on the slope of the transition in the tricritical region (however, it is not possible to see it right away from Fig.2 because of different scales). At this stage, judging by the obtained results, we could draw a preliminary conclusion that there is no decomposition of the TCP in this model, at least for R ≥ 0.01. In order to confirm this claim, we picked up the case of R = 0.01 (it is the most likely candidate for the decomposition, if there is any) and investigated it more closely. In particular, we ran additional simulations for lattices of larger L and finally extrapolated to L → ∞. This data presented on a fine scale allowed us to localize the TCP with fairly high accuracy to (k T /|J |,h /|J |) = (0.0150±0.0001,1.999138±0.000002). Then we examined behaviour B t 2 t 2 of some physical quantities as they approach the TCP. Fig.3 depicts a log-log plot of the 5 direct and staggered susceptibilities vs field near the TCP. While a good fit to a power-low behaviour can be observed in a fairly narrow region quite near the TCP, a more distant region shows a deviation from such a behaviour, probably as a consequence of logarithmic corrections predicted for tricritical point in three dimensions17. Nevertheless, slopes of the both lines seem to asymptotically approach the magnitude close to 1 and 1 for χ and χ, 2 st respectively. This means, however, that both the staggered and direct susceptibilities take on exponents which are very close to the tricritical ones γ = 1 and λ = 1 rather than t t 2 those which characterize an usual critical behaviour. These results just confirm the previous conclusion about non-decomposition of the TCP and put it on firmer ground. In Fig.4 we depict both the tricritical temperature and the tricritical field vs R. While the tricritical temperature shows a linear dependence over whole range of R, which for a certain smaller range of R was also observed in a diluted nnn−model18, the tricritical field displays a slight curvature. As expected, for R → 0, the tricritical temperature moves to zero, while the tricritical field approaches the exact value for the zero-temperature critical field h (0)/|J | = z (= 2). c 2 2 4.Conclusions We have investigated the possibility of the decomposition of the tricritical point in a three-dimensional layered Ising metamagnet. Since the mean-field theory predicts such a decomposition for small values of the ratio R, we mainly focused on that region. However, we observed only tricritical behaviour with no signs of the decomposition for R ≥ 0.01. In the case of the smallest value of R = 0.01, where, according to the MFT, the TCP is most likely to decompose, we managed to locate the TCP with a precision of approximately 1% in temperature and 0.0001% in the field. Even for such a small R the analysis of the critical exponents clearly showed tricritical behaviour. Therefore, we conclude that it is very unlikely that the TCP decomposes for any value of R, although some very small possibility of the decomposition for R < 0.01 still remains. Hence, recently found anomalies in FeBr 2 (note that it has relatively high value of R11,14) make us believe that the high interlayer co- ordination (interlayer superexchange paths present in the real material), supposedly causing the anomalies by inducing local thermal excitations of the second antiferromagnetic phase 6 for weak intralayer couplings J , could play a key role in the possible TCP decomposition. 1 1 E. Stryjewski and N. Giordano, Adv. Phys. 26, 487, (1977). 2 K. Motizuki, J. Phys. Soc. Japan 14, 759, (1959). 3 D. P. Landau and R. H. Swendsen, Phys. Rev. B 33, 7700, (1986). 4 H. J. Herrmann, Phys. Lett. 100A, 256, (1984). 5 H. J. Herrmann, E. B. Rasmussen and D. P. Landau, J. Appl. Phys. 53, 7994, (1982). 6 H. J. Herrmann and D. P. Landau, Phys. Rev. B 48, 239, (1993). 7 J. D. Kimel, S. Black, P. Carter and Y.-L. Wang, Phys. Rev. B 35, 3347, (1987). 8 Y.-L. Wang and J. D. Kimel, J. Appl. Phys. 69, 6176, (1991). 9 M. M. P. de Azevedo, Ch. Binek, J. Kushauer, W. Kleemann and D. Bertrand, J. Magn. Magn. Mater. 140-144, 1557, (1995). 10 J.Pelloth, R.A. Brand, S.Takele, M.M. P. deAzevedo, W. Kleemann, Ch.Binek, J.Kushauer and D. Bertrand, Phys. Rev. B 52, 15372, (1995). 11 H. Aruga Katori, K. Katsumata and M. Katori, Phys. Rev. B 54, R9620, (1996). 12 K.Katsumata, H.ArugaKatori, S.M.Shapiroand G.Shirane, Phys. Rev. B 55, 11466, (1997). 13 O. Petracic, Ch. Binek and W. Kleemann, J. Appl. Phys. 81, 4145, (1997). 14 W. Selke and S. Dasgupta, J. Magn. Magn. Mater. 147, L245, (1995). 15 W. Selke, Z. Phys. B 101, 145, (1996). 16 M. Pleimling and W. Selke, Phys. Rev. B 56, 8855, (1997). 17 F. J. Wegner and E. K. Riedel, Phys. Rev. B 7, 248, (1973). 18 S.Galam, P.Azaria and H.T.Diep, J. Phys.: Condens. Matter 1, 5473, (1989). 7 h TCP h t T T T t N (a) h CE h CE h BCE BCE T T T CE BCE N T (b) FIG. 1: Schematic phase diagrams of a metamagnet displaying two different kinds of behaviour: the second-order and first-order transition lines (solid and dashed, respectively) either (a) meet at the tricritical point (TCP) or, (b) as predicted by the MFT, can end up by a critical endpoint (CE) and a bicritical endpoint (BCE). 8 1.99 2.00 (b) 1.98 (a) 1.97 1.99 J|2 1.96 J|2 h/| h/| 1.95 1.98 1.94 1.93 1.97 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.20 0.24 0.28 0.32 0.36 k T/|J | k T/|J | B 2 B 2 2.001 2.0002 2.000 2.0000 (d) (c) 1.999 1.9998 1.998 1.9996 h/|J|2 1.997 h/|J|21.9994 1.996 1.9992 1.995 1.9990 1.994 1.9988 1.993 1.9986 0.008 0.010 0.012 0.014 0.016 0.018 0.04 0.05 0.06 0.07 0.08 0.09 k T/|J | k T/|J | B 2 B 2 FIG. 2: Phase diagrams near the TCP in the T −h plane for four different values of R: (a) 0.5, (b) 0.2, (c) 0.05 and (d) 0.01. The TCP (arrow) separates second-order transitions (filled circles) from first-order transitions (blank circles). 9 103 c c st st 2 c 10 slope=1/2 , slope=1 c 10 10-5 10-4 10-3 10-2 (h-h)/|J | t 2 FIG.3: Log-logplotofthedirectandstaggeredsusceptibilitiesχandχ ,respectively,asafunction st of (h − h )/|J | for k T/|J | = 0.0150, h /|J | = 1.999138 and R = 0.01. The crosses and the t 2 B 2 t 2 diamonds represent data for two different lattice sizes L = 30 and L = 40, respectively, and the solid lines of slopes 1/2 and 1 (expected tricritical exponents for χ and χ , respectively) serve as st guides to the eye. 10