edited by JOHN W. MOORE
Simulation of Two-Dimensional, Jahn-Teller Phase-Transitions in Solids W. J. A. Maaskant and R. A. G. de Graaff Gorlaeus Laboratories, State University of Leiden, 2300 RA Leiden, The Netherlands
Since the discovery in 1957 by Dunitz and Orgel (1) of phase transitions in transition metal compounds originating from a so-called cooperative Jahn-Teller effect, many examples have been found. These phase transitions occur often around room temperature. The cooperative Jahn-Teller effects are important, since through the deformation of the crystals the exchange paths of magnetic interactions are determined. This paper describes a simulation of two-dimensional phase transitions that mimics part of what occurs in a certain class of cooperative Jahn-Teller compounds. First an introduction is given t o the Jahn-Teller effect, then the types of crystals and the cooperative effects are discussed. Finally the assumptions in relation to the simulation are given. In 1937 Jahn and Teller (2) published their theorem stating that for nonlinear molecules or complexes with degenerate electronic states there are always vibrational modes present that can lift the degeneracy. There are different kinds of Jahn-Teller effects, and, moreover, these may differ appreciably in energetic consequences. In this paper the discussion will be restricted to the so-called E @ 6 effect, which occurs in octahedrallv coordinated Cr2+or Cu2+ions. The capital E refers to thelectronic ground state of these inns. The symbol ( denotes the twofold degenerate E-vibrational mode, umhich occurs in regular octahedra. In Figure 1 the two inde~endentvibrational modesaredrawn.Theseare called QQand Qp. The relative lengths of the arrows on the ligands of the octahedra are given by the following formulae:
picted most easily by showing the potential energy surface in the space of the vibrational modes. Using the linear approximation the Jahn-Teller octahedron has minimum energy for p = po. the radius of the bottom of a trough in the potential energy surface (Fig. 2). Due to nonlinear effects the bottom of this potential energy valley isnot flat but shows wells near the angles 6 O,2n/3, and - 2 ~ 1 3 (Fig. 3). For configurations near these wells the octahedron is elongated along one of its fourfold axes. I t can easily be checked that the elongations along the z , I,and y axes are, respectively,
Figure 1. Normal m d e s Q and C$ of a regular octahedron.
Often it is convenient to designate these modes by their radial forms
and
These vibrational normal modes lower the ground state by lifting the degeneracy. The energy stabilization is in general a function of the cation, the anions, or ligands constituting the octahedron and also of the packing in the crystal. For Cu2+and Cr2+with halide ions as ligands a reasonable estimate of this lowering of the energy is 2000 em-' (see, e.g., Reinen and Friebel (3)).Therefore these octahedra are not regular but are strongly distorted. These distortions are de966
Journal of Chemical Education
Figure 2. Potential energy surface lnthe I&. Teller effect.
Orplane fwthe linear E @ 4 Jahn-
where the nearest-neiehbor interaction is eenerallv the strongest. An ordered state with parallel alignGents is E a ~ ~ e d ferrodistortive. Antiferrodistortive states are also encountered. General reviews on the Jahn-Teller effect are given, e.g., by Englman ( 6 ) ,Reinen and Friebel (3), and Bersuker The fact that these octahedra are preferentially elongated and not compressed is not very well understood, but it is borne out by experiment (3).The barriers between the wells in the potential energy surface are (for the compounds under consideration) of the order of 100-300 em-', corresponding to kT for temperatures of 150450 K. This means that in these compounds phase transitions are to be expected around room temperature, connected with orientational order of the elongation directions. A so-called three-states Potts model of a Jahn-Teller octahedron has been proposed by Hock, Schroder, and Thomas ( 4 ) , which is useful in discussing cooperative phenomena. Each octahedron is assumed to be elongated, along the x , y, or a axis. Steep wells correspond to these elongations, but the free energy of the octahedron is lowered for rising temperatures due to a gain in self-entropy, which is caused by reorientation of the directions of elongation. If there were only two states instead of three, the correspondence with the familiar Ising model (Newell and Montroll (5))would have been exact. In crystals containing an appreciable number of JahnTeller ions of this type, an interaction between those octahedra tends to direct each others elongation. This may occur via strain, where the parameters of the original unit cell are altered. I t may also occur by way of neighbor interaction,
rn. , ~
,
~
The crystals to which the computer simulation applies are hexagonal perovskites, with the formula ABX3. Here A = Cs+ or Rh+, B = C U ~or + Cr2+, and X = C1-, Br-, or I-. In Figure 4 the crystal structure of the ideal hexagonal perovskite is given. Clearly visible are the columns of face-sharing octahedra, which contain a t their centers the Jahn-Teller active ions (B), which replace Mg. The six vertices of these octahedra are formed bv the X ions. The X ions form. together wirh the A ions, so-called close-packed layers. In the planc pwpendicular to the direction of the n ~ l u m n s(the c axis), the R ions form a hexagonal pattern. Each R ion has, therefore, six nearest neighbur B ims. At high trmperatures these crystills hare an averngwl structure of the idea1 hexapP&/mmc. When onal perovskite type with the space aroup . the temperature &lowered, phase transitions t o lower symmetry occur, due to the interactions of the Jahn-Teller octahedra. Crama and Maaskant (8) give a comprehensive account of the different possible orderings that are theoretically expected and that are found in nature. The computer program to he discussed refers only to phase transitions in a plane perpendicular to the c axis of the ldescribed crystal type. Effects of strain will not be treated; i t is assumed that only nearest neighbors interact. Crama and Maaskant (8)derived the number of interaction parameters using group theory. For the simulation, however, a simplified interaction is assumed; the interaction energy is -V for parallel nearest neighbors and 1'2 V for differently oriented neighbors. The interaction is proportional t o the "scalar" product of the linear combinations of the normal coordinates (eqs 1 4 ) . For example, Q,Q, = 1and QxQ, = -'I%. Computational Details
1 Q3
Figure 3. Warping in the lower pdential energy suriace by nonlinear JahnTeller coupling.
During the simulation a hexagonal array of 144 sites is shown, arranged in 12 rows and 12 columns. The sites correspond to the Jahn-Teller centers of a plane of the crystal structure described in the previous section. Each site displays an X, a Y, or a Z, representing the elongation direction of that particular site. For the sites a t the edges, periodic boundary conditions have been imposed. Therefore all sites are equivalent and have six nearest neighbors a t identical distances. The method of calculation consists of determining the Boltzmann factors for the three possible orientations on a site, the interaction energy being -V between identical neighbors and %V between different centers, as explained previously. Suppose a certain site has n, X-oriented, n, Yoriented, and n, Z-oriented neighbors (n, n, n, = 6). Then, according to Boltzmann, the probability of having an X a t the central site is P , = exp I(-n, + %n, + %n,)VlkTJ (6)
+ +
Here k is the Boltzmann constant and T the absolute temperature. The prohabilities of having a Y or a Z, p, and p,, can be obtained by cyclical permutation from p,. Usually the sum of these prohabilities is normalized to 1. In order to determine the actual orientation of the central site, a random number between 0 and 1 is chosen. If this number is smaller thanp,, the orientation is X. When this number isp, or between p, and (p, p,) the orientation will he Y. The orientation will he Z when the number is equal to or larger than (p, + p,). This calculation is done for subsequent sites in a row and then for the next row etc. On completion, the results in memory are moved to the textscreen. Then the additional rows and columns for the periodic boundary conditions are filled and the calculation is started again. I t is not difficult to see that only a limited number of
+
u
.w
Figure 4. me ideal hexagonal psrovskite crystal structure of CsMgCb. One column has been omitted tor reasons of clarily
Volume 63 Number 11 November 1986
967
possible configurations of the nearest neighbors exist for a given site. When n,, n, and n, are given in this order, the following sequences are found: (6001, (510), (501), (4201, (4111, (402), (330), (321). (3121, and (2221, and, exceptforthe last combination. the cvclical oermuted ones. Since for a particular value of the interaction constant this calculation is repeated 144 X 300 times, the procedure is speeded u p considerably by calculating the Boltzmann factors for these confieurations onlv once and storing the results in a tahle. calcula~onswith integers are much faster than floating-point calculations. Therefore, instead of normalizing the total probability on 1, the probabilities were multiplied by 16777215 (256 X 256 X 256 -1) and the respective values were stored as an unsigned three-byte number in the table. The nromam ~ . . consists of the actual simulation, which is repeated 300 times, and of a stage in which a plot isshown of the correlation roefficient (vide infral, which is calculated for each interaction constant. The plot is on the high-resolution screen and requires some time t o he read. During that time the calculation of the Boltzmann factors for the next interaction value is done. In this type of program, which uses a so-called Monte Carlo method, the quality of the random numher generator (in reality a pseudo-random-number generator) is of the utmost importance since it determines the finer details of the behavior of the system near phase transitions. Since in this case the program serves for demonstration purposes, a second requirement is that i t should work fast. Therefore, the generator works as follows: A full page (256 hytes) is filled with random numhers, each between 0 and 255, created with a Basic random number generator. Since the probabilities, the Boltzmann factors, are approximated by unsigned integers of three hytes, the random numbers must be expressed similarlv. The required random numbers are composed of three su&essive bytes from the page mentioned. After having used a particular triplet of bytes, the next set will start with the middle byte of the previous set. The first byte and the second of the used triplet is added modulo 256, and the result olaced in the first bvte. thus ~rovidinea revised page of random hytes. Since for one interaction constant 300 X 144 random numbers are selected, the Dace of hytes is used very often. The quality of the distrib&n of numbers was checked by counting the numbers between 0 and 15, between 16 and 31, etc. Though at the start the distribution was relatively flat, after a certain number of passages through the pa& i t became less so. However, the procedure used did not affect the results in a noticeable way and an improvement, which certainly would take more time, seemed not to he required. Each time a new interaction constant is given, the simulation starts with the same random set of X's, Y's, and Z's. This has the advantage that the new equilibrium situation is found by the calculation itself. Since this is observable, i t is more convincine. Near the arrav of orientations, the number of iterations, t h i number of X's, Ys, and Z's, and the interaction constant are given. For each value of the interaction parameter 300 iterations are carried out. After the first 200 it is assumed equilibrium has been reached. During the last 100 passages the correlation (local order) between nearestneighbor sites is calculated. Two identically oriented neighboring sites give a value 01 1; two that are dillerent g k e - l 2 . The numhers nre added and in the Basir program (floatingooint mlrulation) di\.idrd b\. the numl~erof pairs of neiehbors for the while array, during 100 cycles. The average value of this local order. the correlation coefficient. has a maximum value of +1, corresponding to total alignment of the orientations, and a minimum value of -'I%, which corresponds to an antiferrodistortive ordering. This correlation coefficient is plotted on the high-resolution screen as a function of the interaction parameter. In the next section this correlation coefficient will be further discussed.
e or cover .
966
Journal of chemical Education
Discussion The program starts with a menu. The default values provide for a change of the interaction constant, which is in reality (VIkT), from the value -2.0 in steps of 0.1 to the value of +2.0. The whole of the simulation requires about 40 min, which in practice may be too long. For this reason the menu offers other possibilities. The initialization takes about 1min, which also may last too long in a lecture. Therefore, on completion of the initialization process, the word "start" is shown on the screen. On pressing the return key the simulation starts immediately. As mentioned, the int~mctionconstant, which implicitly contains the temperature T, has a positivr and n negative ranye, which urt. scanned in steps from negative t o positive d u e s (Fie. , " 51. , This is unrealistic. sinre in orartice the interaction constant is either positive or negative. So to begin with the IVlkT) is neeative, and it can be imaeined that for . . rising temperature iis influence is graduall; diminished. This is imitated bv increasing the value of the effective interaction constant toward 0. The point 0 would then correspond to infinite temperature. On the positive side a further increase of (VIkT) from 0 to 2.0, corresponds to the falling of the temperature. The results are clearly demonstrated on the textscreen and in the plot (Fig. 5) produced during the process. For large negative values of (VIkT), neighboring sites tend to he different in their orientations. The process converges to an ordered structure. which is antiferrodistortive (Fig. 6a). I t should he pointed'out here that the choice of the numbers of rows and columns, 12. is not arbitrary. If a numher unequal to a multiple of 3 is 'chosen, problems would arise due to mismatchiue. The correiation coefficient, shown on the plot is exactly -'I?. for the most negative values of (VIkT). When the effective interaction constant increase>,a certain disorder is developed, which is seen on the text wwrn !I., changing orientation axes and by the total numbers of different orienta-
Figure 5. The correlation coefficient(vertically)against the effective interaction constant (horirantally).
Y
L
X
Y
Z
Y
X
z x y z x y x y z x y z y
z
x
y
z
Y
Y
x
Y
z x y r x y x
y
z
x 0
y
z
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
Y Y
b
Figure 6. Part of the antiferrodistortiveordering (a)and of the ferrca'islortive ordering (b).
tions. Also the correlation coefficient increases on the olot. '['his hecomes dmmnrir at the phase trnnsition, after which, the nrrrelation coefficient inrreoseaauicklv toward 0. For an effective interaction constant with the vaiue 0, the correlation coefficient is also 0, as is to he expected. Thereafter the correlation in~reasesfurther, and suddenly much steeper when a second phase transition sets in. This is a transition toward a fully aligned ordering (Fig. 6h). The correlation coefficient then mounts up to +l.I t should be noted that i t is comoletelvirrelevant which orientation is found: all X. all Y. or ali 2. 1" particular it is intereiting ro see during the'fairli steep riseofthe n~rrelariuncoefficienthou, the domains with different orientations compete for domination. As was stated in the former section, it is assumed that after 200 simulations the system is in equilibrium. However, occasionally two domains are equally strong and are maintained up to the end of the simulations. This is the phenomenon of twinning, which occurs naturally in the type of crystals under consideration. The case in point results in an anomalous value for the correlation coefficient. A remark should he made on the differences in hehavior near the ferro- and the antifimodistorrive phase trnnsitions. The rise of the correlation coefficient is much steener for the ferro- than for the antiferrodistortive transitions. The reason lies in the number of possible confieurationsas discussed in the former section. suppose a site shbuld have an X in the end. Then only one configuration, (600), will have the lowest ~
energy in the ferrodistortive case. For the antidistortive ordering, however, when the central site should have an X orientation, there are several configurations with equal probability: (0,6,0), (0,5,1), (0,4,2), (0,3,3), (0,2,4), (0,1,5), and (0.0.6). In summary a demonstration of a simulation of a threestates Potts model, a variant of the Ising model, for a twodimensional structure has been described. The program runs on an Apple .. II+/e comDuter. The simulation is related to phase transitions in a class of inorganic compounds, showing a cooperative Jahn-Teller effect. The .oroeram. . . . ~. a r t.l v written in Rasic, but mostly in Assemhly language is available on a 9,-in. disk from Prvjecr SI?HAPHIM. Acknowledgment The authors gratefully acknowledge interesting discussions with W. G. Haije. Llterature Cited
~
-.
Jahn~TellerEffect in Molecules and Cryst&"; Wiley-I
7. Borsuker. I. B. "The Jahn-Teller Elfact and Vibronie In+raetions in Madern Chemisfry"; Plenum: New York,1981. 8. Crsma, W. J.; Maaskant. W. J. A. Physics 1983.121B.219.
Volume 63 Number 11 November 1986
969
.
