Phys Chem Minerals (1992) 19:196-202

PHUlCI I CHEMISTRY NMIIi[IIAL6 9 Springer-Verlag 1992

The Energetics of Interaction Between Oxygen Vacancies in Sillimanite: A Model for the MuUite Structure S. Padlewski 1, V. Heine 1, and G.D. Price 2 1 Cavendish Laboratory, Madingley Road, Cambridge CB30HE, UK 2 Department of Geological Sciences,University College London, Gower Street, London WC1E 6BT, UK Received November 25, 1991 / Accepted April 25, 1992 Abstract. A microscopic model is introduced to discuss the modulated structure of mullite. The oxygen vacancies of this aluminosilicate are known to play a central role. In particular, a single vacancy strongly orders its surrounding A1/Si tetrahedral sites. It is shown in this work that if two oxygen vacancies approach too closely to one another, their A1/Si "dressing" overlap. This situation results in repulsive interaction. The field of interaction between the vacancies is estimated with the use of an atomistic computer simulation. We use a Bragg-Williams type of theory to dicuss the ordering pattern of the vacancies. Due to "frustration" between the two dominant repulsive interactions, our model predicts a modulated phase transition in agreement with observation.

1 Introduction Mullite is a mineral of composition A12(A12+ zxSi2 2x)O10 where the compositional index x lies in the range 0.17
-

x

2 Ni4 + + O 2 - --* 2 A13 + + vacancy,

We have discussed in a recent work the local aluminium/silicon distribution around a single oxygen vacancy. We found that an extremely strong ordering scheme is established around it (Padlewski et al. 1991). We called such A1/Si ordering the "dressing" of the vacancy (Fig. 2). On formation, the vacancy releases its two adjacent cations (the so-called T* atoms). They repel one another to T* sites by adjacent oxygen atoms at the so-called O* sites (Fig. 1), that consequently become coordinated by three cations (Fig. 1). We found that these so-called 3-clusters contain exclusively aluminium ions for coulombic reasons and are accompanied in the adjacent ab-planes by S i - O c - S i groups above and below the vacancy (Fig. 2). The vacancy with its "dressing" of ordered A1/Si atoms is considered in this study as a modular unit. This approach is quite natural because the ordering inside such a complex is energetically very favourable. Moreover the whole A1/Si "dressing" carries the appropriate stoichiometry, which is constrained by the substitutional reaction (1.1). The oxygen vacancy and its A1/Si "dressing" is an object which has a large and extended spatial extent (Fig. 2). When two vacancies approach too closely, their

(1.1)

where Oc is the oxygen site where the vacancy can reside (Fig. 1). The average (i.e. disordered) structure of mullite consists of chains of A10 6 octahedra running along the c-direction which do not participate directly in the ordering process (Angel and Prewitt 1986). The chains are cross linked via bonds of type T - O ~ - T , where T represents a tetrahedral site occupied by silicon or aluminium.

IlL Fig. L Structure around an oxygen vacancy placed in the centre (shaded) of a supercell of sillimanite. The two adjacent tetrahedrally coordinated T atoms migrate to T* sites near the O* atoms which are then in the middle of three T atoms forming a 3-cluster

197 central plane

adjacent planes

T

I Al'-~r*

---~'3

bT

T

T

T

T

Si

\ br

AI*

[

a

AI*

AI* bt

AI

AI*

AI

Y

T

Al* C

a

d

v

Fig. 2. A1/Si ordering around a single oxygen vacancy (shaded). The two AI* atoms would have been linked if the central O0 site had been occupied by an oxygen atom instead of the vacancy

T 1IF

T

A

s,]

D

AI*I

B

Si*

Si

AI*

Fig. 3. Vacancies (full squares) and their "dressing" regions (shaded). A and B: isolated vacancies. C, D and E: vacancies whose dressings overlap

Fig. 4a-e. Energetically highly unsatisfactory situations arising when two vacancies approach too closely. In (a) a cluster with four ions is formed. This situation is unfavourable because of the excess charge. In (b) and (e) one of the four At* atoms has no O* atom to attach itself to. In (d) the vacancy between the silicon atoms results in two energetically expensive (T, T, Si*) clusters and finally in (e) the unlikely (Si, Si, AI*) duster is formed. The numbering of the sites follows that of Fig. 2

mutual "dressings" overlap and result in interaction (Fig. 3). The interaction is due to local strain and Coulomb forces, but the dominant effect is due to the problems with the AI* atoms liberated by the vacancies. Five situations are shown in Fig. 4. Consider first the vacancy shown shaded in Fig. 4(b), with its two Al* atoms. If another vacancy is placed in position 2 in Fig. 4(b) (the numbering of positions follows that of Fig. 2), then one of its AI* atoms has no O* a t o m to attach itself to where the central vacancy is. Instead it has to wander off to some less favourable site out of the ab-plane: in any case it is an energetically unfavourable situation. The absence of an appropriate O* a t o m to attach to also applies in Fig. 4(c). Conversely in Fig. 4(a), one AI* a t o m from a second vacancy in position 1 would join the AI* a t o m from the first (shaded) vacancy to form a 4-cluster, i.e. with four A13+ ions around one (O*) 2- ion, which would have a high net charge and very high C o u l o m b

energy as shown by Padlewski et al. (1991). Figure 4(e) shows a second vacancy in an ab layer adjacent to that of the first vacancy: its AI* ion would form a 3-cluster but this would contain two Si 4+ ions which again is energetically unfavourable because of the large total charge (Padlewski et al. 1991). The important result of the above discussion is that there is a strong repulsion for a second vacancy in position 1, 2, 3, 4 or 5 in Fig. 2 with respect to the central vacancy. The strong repulsion is represented schematically in Fig. 3 by the overlap of the striped "dressings" around the vacancies. In Sect. 3 the detailed interaction between vacancies will be worked out, but the strong repulsion just described is the dominant part. As a guide to the detailed treatment to follow in Sect. 4, we now sketch how the repulsive interaction between vacancies results in an incommensurate structure. We write the energy as the function

a

198 1 2

2•p,

(1.2)

J1

J1

3

where P~ is the probability of finding a vacancy on Oc site i and E({Pi}) is the energy per unit cell (a.b.c) of the disordered sillimanite structure containing four T atoms per unit cell. The compositional index x fixes the average concentration of the vacancies as

p = l x.

-ira

,,

(1.3)

The sum in (1.2) runs over the N sites Oo of the crystal and we need to introduce the factor 2IN for proper normalisation (two Oc sites per unit cell). In expression (1.2), the vacancy is formally considered as a defect so that the whole energy is an expansion around the defect-free sillimanite material (x = 0) (Holland and Carpenter 1986). The term Es is the energy per unit cell of sillimanite (P = 0) structure, while Er is the cost of the substitutional exchange reaction (1.1) and is formally the energy to introduce a vacancy into the structure. Finally the term C~j (Fig. 6) represents the interaction between two vacancies placed on the Oc sites i and j respectively, which we determine by computer simulation in Sect. 3. The ordering of the model (1.2) near the transition temperature can be approached by expanding the probability distribution P~about the average:

P~=_P+ A~.

(1.4)

Here Az>0 denotes an excess of vacancies and A~<0 an excess of oxygen. A full analysis will be carried out in Sect. 4, considering all directions of the modulation as well as commensurate superlattices. For the moment let us assume that A i varies in the a-direction, so that we can characterise Az on a whole bc-plane of vacancy sites by a planar variable A, for plane n. This variable characterises the vacancy density including its variation with position on bc plane n. We shall see in Section 4 that there is a simple modulation on each plane but the only point relevant here is that its form is the same on each plane apart from sign which is included in A,. The energy of the system can then be written 12 E({A,}) = c o n s t a n t - 7 - - y' ~ Jr A, A,+r g

n

(1.5)

r

The new coupling terms Jr with r = 1, 2 (truncated after second neighbour) are linear combinaisons of the relevant Czj from equation (1.2). Note that consecutive planes of Oc sites are separated by a distance a/2 (Fig. 2) so that a structure modulated in the a-direction with wave vector q can be represented by

A. = A cos(q na/2). Inserting (1.6) into (1.5) yields

(1.6)

J2 Fig. 5. Signs of possible A, on three consecutive planes. Starting with a positive sign on the left layer, the ./1 and Ja (each being negative) tend to produce opposite signs on the right hand layer

The coupling terms 31 and J2 will both be shown to be very repulsive, i.e. negative, and of comparable magnitude, for the reasons already discussed in connection with Fig. 4. Under such circumstances it is easy to show (by differentiating E(q) with respect to q) that the energy (1.7) has a minimum at some incommensurate value q, given by

(~)

cos q.

-

IJll

4J J21"

It is common (Yeomans 1988; Selke 1988) to interpret this result pictorially in the manner of Fig. 5 representing the signs of some Ar on consecutive planes. Suppose dr is positive on the left hand plane. Then J1 being negative wants (energetically) the sign of A to be opposite on neighbouring planes, which by two hops would result in a positive sign on the right hand layer of Fig. 5. However J2 is also negative, so that the energy is lowered if A has opposite signs on next nearest neighbour planes. That would tend to give a negative sign on the right hand plane of Fig. 5. Such a situation is known as frustration. It shows that any simple ordering pattern is not very satisfactory, and a deeper analysis then gives a modulated structure, at least at the level of (1.7). Equation (1.7) and Fig. 5 demonstrates in outline the origin of the incommensurate structure of mullite. In fact using the values of Cij obtained from computer modelling in Sect. 3, we shall in Sect. 4 obtain the modulation wave vector qlc = 89c* _+0.69 a*.

(1.9)

The c*/2 in (1.9) arises quite simply: Fig. 2 already suggests that vacancies like to be spaced at a distance 2c apart in the c-direction so that they can "share" the pair of Si atoms on the intervening a b-plane. By subtracting/adding the reciprocal lattice vector a* we finally obtain from (1.9) the "observable" wave vector inside the first Brillouin zone as would be reported from a diffraction measurement Qlc = 89c* +_0.31 a*,

E (q) = constant-- ~ -

(1.8)

which has the right order of magnitude.

(1.10)

199 There is another interpretation of the origin of the incommensurate structure, and in particular why the modulation is along the a-direction. McConnell and Heine (1984; also Heine and MacConnell 1984) have shown that any incommensurate modulation can be expressed in term of two component difference structures C1 and C2 which are modulated by sine and cosine functions respectively and which interact energetically. The point is that, given the direction of the modulation wave vector q~c, there is a unique symmetry relation between C 1 and C2, which has been discussed for mullite by McConnell and Heine (1985). The symmetry of the A1/Si ordering around a vacancy (Fig. 2) sets the symmetry of one of the components; let that be C2 to be consistent with McConnell and Heine (1985). If qlc is along a*, then the symmetry of C1 implied by this C2 is precisely that of the sillimanite ordering pattern. For a low concentration of vacancies we can think of the C 1 material between the C2 regions as just being sillimanite (with some disordered vacancies). Such an overall ordering pattern clearly makes sense intuitively. The important point is that C1 would not turn out to have the sillimanite symmetry if q~c were along a different, high symmetry direction such as b* or c*. The above analysis of (1.6) to (1.10) is only valid near the ordering transition temperature To, which can also be estimated from C~j: it turns out to be well above the melting point of mullite which suggests why the IC phase of mullite is so strongly ordered and stable. Our analysis also does not allows any variation of Qlc with x, contrary to observation. These are matters requiring more detailed theory and calculations which will be taken up in further work.

2 Computer Simulation Atomistic computer simulation has already proved to be of use in many fields of solid state physics and mineralogy. In many cases, it is the only tool available which enables microscopic processes to be inferred. We have recently used this approach (Padlewski et al. 1991) in order to determine the local aluminium/silicon ordering scheme around a single oxygen vacancy in sillimanite. This could not be deduced clearly by experiment because the X-ray scattering factor of A1 and Si are very similar, and also because some difficulties arises when one tries to analyse diffraction pattern of a system showing important substitutional disorder (Welberry 1985). As in our previous paper, we use in this work a static lattice energy calculation (Leslie 1984). In this all ions of a supercell (one or several unit cells of sillimanite) are placed on their approximately correct positions. The supercell is reproduced indefinitely by periodic boundary conditions and all ions interact with each other via two and three body potentials. For the two body potential, we use a Buckingham potential which includes Van der Waals attraction and an exponential repulsion due to electron cloud overlap. A cut off is introduced beyond the second nearest neighbour distance. A formal charge model is used for all ions and the Coulomb energy is

calculated with the Ewald technique. We also use a shell model (Catlow et al. 1982) for the oxygen sites where "core" and "shell" interact harmonically. Angle dependent three body interactions are used to model the covalent bonding of the tetrahedrally coordinated aluminium and silicon to their nearest neighbour oxygen (Lewis and Catlow 1985). The empirical potentials used in this study are the same as that of Bertram et al. (1990). In this study, we used sillimanite supercells of different size, up to eight unit cells. Here a "unit cell" of dimension (a.b.c) refers to the cell of disordered sillimanite, i.e. containing four T sites and two Oo sites in one layer. We constructed supercells with different geometries, for example (2a.2b.2c), (4a.b.2c) and (a.2b.2c). Into these supercells we introduced, according to the reaction exchange (1.1), one or two oxygen vacancies. In each case, the aluminium and silicon were arranged around the vacancy(ies) according to the ordering scheme indicated in Fig. 2. Note that even with one vacancy per supercell we obtain information about vacancy-vacancy interaction because the vacancy in one supercell interacts with those in neighbouring supercells. Away from their zone of influence the A1 and Si ions were arranged in the sillimanite ordering scheme which is known to be the lowest energy arrangement in the absence of vacancies (Bertram et al. 1990). However when two vacancies are close one to the other, their "dressings" overlap (Fig. 4) and the best A1/Si arrangement around them was found.

3 Field of Interaction Between the Oxygen Vacancies In this section, we discuss how the computer simulation enables us to estimate the coupling terms C~j which occur in the energy parametrisation (1.2). Note that the interactions Cij in (1.2) include all direct and indirect effects, e.g. the atomic displacements and A1/Si ordering in the region where the dressings of the vacancies overlap. Since the dressing around each vacancy is relatively short ranged, the interaction energies C~j will tend to zero when the dressing no longer overlap. Thus we only consider the Cij defined in Fig. 6. Some of them do not result in overlapping effect, i.e. C3ao, C13o or C4oo. They were found to have a smaller interaction energy. More distant ones have been set to zero. One could in principle simply do a least-squares fit of the eleven coupling terms considered in Fig. 6 to the calculated energies of a large number of vacancy configurations. However since the number of configurations needed to be kept small because of computer time, we designed them in order to obtain specific combinations of the C~i such that each C~ can be determined. For the following we only calculated the total energy of fifteen different configurations. For example in order to obtain Coo~ we used a supercell of size (2a.b.2c) and another one of size (4a.b.2c). The first one is presented explicitly in Fig. 7 where we also indicate the first neighbour vacancies which are reproduced by the periodic boundary conditions. A similar picture provides the A1/ Si distribution of the other supercell (4a.b.2e). By using the energy expression (1.2) and also the coupling terms

200 Table 1. Numerical value obtained for the coupling between the vacancies. The four first terms result from overlapping effect and are all repulsive as expected Coupling eV

Cllo

C22o

Cool

Clll

Czoo

C~oo

- 1.02

-0.36

-0.72

-0.17

+0.20

+0.04

Coupling eV

Co2o +0.22

Co4o +0.05

C31o -0.02

Coo2 +0.19

C13o +0.08

;05 v a

a+b

Fig. 6. Field of interaction around a central oxygen vacancy (shaded) that are considered in our modelling. The interactions in bold character result from the overlapping effectsof Fig. 4

+

I-

1

ls i A A AI

AI

central

plane

AI

AI

AI

|si Is' L _ o "l / Si

adjacent

/ si c~-I

planes

Si

a

Fig. 7. A1/Si ordering in a block of 2a.b.2e supercells with one vacancy per supercell (square box). In bold, the T and Oc sites belonging to the central supercell. The others are reproduced by periodic boundary conditions. Some neighbouring vacancies which can interact with the central one (shaded), are represented

defined in Fig. 6, the energies of the two supercells take the forms E1 = 4 E s + E v - Co zo - Co4o - Coo2 - C4oo,

(2.1)

E2 = 8 E s + E v -

(2.2)

C020 -- C040 -

Coo 2 .

The computer simulation enables us to calculate the energy of the both supercells and we respectively find E1 = -2214.67 eV and E 2 = -4526.35 eV. We also need to calculate the energy of the defect-free sillimanite structure, for which we find Es = -577.93 eV per unit cell (a.b.c). The coupling term Coo 4 is directly deduced from (2.1) and (2.2): C 4 o o = E2 -- 2 E s - -

E1 = + 0.04 eV.

(2.3)

The interaction in this case is positive, i.e. attractive, and also very weak because there is no overlapping. We generalised this procedure for other coupling terms, and used fifteen supercells of different size and vacancy(ies)

distribution to obtain the eleven exchange interactions of Fig. 6 in a similar manner. The results are presented in Table 1. Note that the energies of the interaction are significantly higher for those with direct overlap of the main dressings shown in Fig. 4. This result justifies the cut-off that we used for Cu: we find Clto/C4oo,,~25 for example. For symmetry reasons (i.e. the space group of mullite Pbam does not contain the element my) the exchange interactions C22o and C2-2o (not shown in Fig. 6) for example are not formally equivalent. In these cases, we only consider their average effect. The first term Cllo is the most repulsive. As already mentioned it is the average of Cllo which corresponds to the situation of Fig. 4(b), and C~_~0 which corresponds to a very similar situation when a vacancy is placed in position 3 of Fig. 2 (see Fig. 4(c)). In each case, the T* atom from the second vacancy does not find an Oo atom (because the first vacancy is there) to attach itself to. Close inspection shows that after the relaxation the cation is pushed toward the Oc site adjacent to the vacancy, on the next a b-plane. This situation results in a distortion among the columns of A106 octahedra (Fig. 1). The next term C22o of Table 1 is also repulsive corresponding to the geometry of Fig. 4(a). As mentioned in Sect. 1, the resultant 4-cluster of two T and two T* ions around one Oo atom is energetically unfavourable for coulombic reason (Padlewski et al. 1991). The term Coos is strongly repulsive while Coo2 is attractive. This can be explained quite simply: in Fig. 1 the S i - O o - S i group is situated about the Oo (4) site because such ordering promotes local charge balance. In the simplest picture, we can attribute an excess negative ( - ) charge to the vacancy site and excess positive charge ( + ) to the S i - Q - S i bond. If we place two vacancies along the c-direction, the configuration giving a " +-+ " sequence, i.e. two consecutive vacancies along the c-direction surrounded by silicon bonds at the extremities, is less favourable than the " + - + - " sequence with alternating vacancies and silicon bonds. This picture is in agreement with the interactions Coo~ < 0 and C o o 2 > 0 of Table 1. Finally, the term C~a~ occurs when a vacancy is placed on the site Or in Fig. 4(e). This situation creates a 3-cluster of type (Si,Si,AI*) which is known (Padlewski et al. 1991) to be somewhat energetically unfavourable. Note also that the coupling along the b-direction and a-direction are not the same (orthorombic structure). The first direction appers slightly more favourable for aligning the vacancies along (Co2o > C2o0) than the other.

201

4 OrderingofVacanciesandIsing-TypeModel

As already mentioned in Sect. 1, it is not a proper treatment of the free energy and does not explain the variation of Q~c with the composition index x, a matter that will be taken up in a subsequent more thorough analysis. In the previous section, the analysis of Table 1 led us to infer that the vacancies tend to lie along lines in the b-direction (Co2o, C040>0), with the density of vacancies on one line depending on the composition x and the temperature. These lines are likely to be spaced a distance 2c (C0oa<0 and Coo2>0). Hence each b c atomic layer can appear in two different ways according to the phase of the modulation, as shown in Fig. 8. We can describe both types of layer by writing

In this section, we discuss the ordering scheme that we would expect in mullite from interactions C u discussed above. We can rewrite the energy (1.2) in term of the variables Ai that are defined in (1.4). The concentration of vacancies is linked to P by (1.3) so that the variables Ai need to fulfil

~Ai=0.

(4.1)

apart

i

The total energy per unit cell (a.b.e) takes the form C o j -1~ 2 cuaiAj.

E({A,})=Es+2EvP-p2~ j

i

Ai

j

on layer n = A.

COS(~. ri),

(4.3)

(4.2) with A. being positive or negative respectively for the two types of layer (see Fig. 8). Note that our notation applies equally to the atomic layers containing Oo sites with X = 0 and with X = a/2 in the unit cell. Equation (4.3) is the precise definition of A. used in the energy expression (1.5). The corresponding Jr in (1.5) are given by

What kind of ordering can emerge from the above model? We follow below the simplest type of heuristic analysis to establish the occurrence of a modulated structure.

z

z

more vacancies; A i>0

2c

2c

less

vao~tncies

A i>O

less vacancies A i<0

c

more vacancies

more vacancies A i>0

0

less vacancies A i
Ar negative

positive

Fig. 8. Two types of bc-layers of ordered vacancies. These layers interact with each other with repulsive -/1 and J2 coupling terms so that a modulated structure emerges

I

I

I

1.2 eV (4.4)

The main contribution to d~ is given by the situation depicted in Fig. 4(b, c) and is therefore very repulsive as we discussed in Sect. 3. The next term J2 is also very repulsive because of the unlikely 4-cluster of Fig. 4(a). Note that d3 and J4 do not involve any direct overlapping of the vacancy dressings and are for that reason an order of magnitude smaller.

~y

~y Ar

J1 = 2 C a l o + 2 C 1 3 0 - 4 C l 1 1 = J: = 2 C220 + C2o0 = - 0.52 eV J3 =2C31o = - 0 . 0 4 eV 9J4 = C 4 0 o = q - 0.04 eV.

Ai
J

In

d

'

I

,

v C) I v

T

T I

0.2

,

.1

I

0.4

I

06.

qa

I

I

0.8

Fig. 9. Fourier transform d(q) of the coupling terms between the vacancies along the line qa a*+0.5c*. The maximum fixes the wave vector of the modulated structure at the transition temperature (here about q~c~0.5c * + 0.69 a*)

202 We have therefore justified the energy expression (1.5) and fleshed out the discussion (1.5)-(1.8) in Sect. 1 of the origin of the incommensurate modulated structure. The structure is a modulation in the a* direction of the whole pattern in the bc atomic layers shown in Fig. 8, which is in good agreement with the structure observed by Cameron (1977). Note that the wave vector e*/2 from (4.3) combines with the modulation wave vector q~ (1.6), (1.8) in the a* direction to give a modulation with the total wave vector qic given in (1.9). An alternative way of presenting this result is to consider quite generally

A~=A cos(q.ri+ O),

(4.5)

which encompasses both incommensurate modulation and single superlattice ordering. Here q is a general wave vector in the reciprocal space and O is a phase. The energy (4.2) with (4.5) becomes E (q) = constant - 1 j (q) A2

(4.6)

where J (q) = ~ Coj cos (q (r i - rj)). (4.7) J Clearly to minimise the energy (4.6) with respect to q, we want to maximise J(q) in (4.7). In fact we find that the maximum of J(q) occurs along the line q, a* +e*/2 at the value of qtc already given in (1.9). We show in Fig. 9 the value of J(q) along this line. Note that it is just the quantity in square brackets in (1.7) and that the maximum at an incommensurate value of qa, i.e. not at the end point q , = 0 or 1, arises because Jz is negative and is substantial in magnitude (to be precise I321> JJlJ/4 as follows from (1.8)). In fact, the wave vector (1.9) is not only maximum along the direction q,a* + e*/2, but all over the Brillouin zone. In further work we shall show that Tc over the whole range of composition where mullite exists (Cameron 1977; Angel and Prewitt 1986, 1987; Angel et al. 1991) is well above the melting point. For example, we find that at x = 0.3 the transition temperature is about T~ ~ 5000 ~ C. Despite overestimate due to mean field theory (a factor of about 3/2 higher than the real Tc) such a temperature is much higher than the melting point Tm~1850~ C of mullite (Holland and Carpenter 1986). This is why the disordered structure is never observed in reality. Finally the wave vector Qlc is experimentally found to depend on the composition (Cameron 1977) which our simple analysis above does not explain. Our discussion does not mention temperature, entropy or free energy: it is purely heuristic. However a proper theory will show that the treatment here is valid for sinusoidal modulation just below To, but that to obtain the variation of Q~c with x it is necessary to consider the squaring up of modulation at lower temperature.

5 Conclusion

Computer simulation has enabled us to determine the complete field of interactions between two oxygen vacan-

cies in sillimanite in order to model the mullite structure. It is known that a single vacancy strongly orders its surrounding A1/Si tetrahedral sites which we called its dressing. If two vacancies appraoch one another, their dressings overlap. In particular if they appraoch very closely, there is a very strong repulsion because of mutual interference with the T* atoms as discussed with regard to Fig. 4. F r o m the picture of the vacancy-vacancy interactions, we have discussed the ordering pattern of vacancies which minimises the total energy. The solution turns out to be a modulated distribution of b c-layers (see Fig. 8) along the a-direction. The simple analysis here is only valid at the transition temperature but already gives quite good agreement with experiment. We also find T~ well above the melting point which explains the stability of mullite.

Acknowledgements.We wish to thank Prof. J.D.C. McConnell, Dr. R. Angel, Dr. I.L. Jones and Dr. M. Dove for useful and helpful discussions. We acknowledge the NERC grant GR3/6970 for providing computer facilities. One of us (SP) also wishes to thank the French Ministry of Research and Industry for financial support.

References

Angel RJ, Prewitt CT (1986) Crystal structure of Mullite: A reexamination of the average structure. Am Mineral 71:1476 1482 Angel RJ, Prewitt CT (1987)The incommensurate structure of Mullite by Patterson synthesis. Acta Cristallogr B 43:116-126 Angel RJ, McMullan RK, Prewitt CT (1991) Substructure and superstructure of Mullite by neutron diffraction. Am Mineral (submitted) Bertram UC, Heine V, Jones JL, Price GD (1990) Computer modelling of A1/Si ordering in sillimanite. Phys Chem Minerals 17:326-333 Catlow CRA, Dixon M, Mackrodt WC (1982) Interionic potential in ionic Solid, vol 166. Springer, Berlin Heidelberg New York, pp 130-161 Cameron WE (1977) Mullite: A substituted alumina. Am Mineral 62:747-755 Heine V, McConnell JDC (1984) The origin of incommensurate structures in insulators. J Phys C 17:1199-1220 Holland TJB, Carpenter MA (1986) Aluminium/silicon disordering and melting in siltimanite at high pressure. Nature 320:151-153 Leslie M (1984) Daresbury Laboratory International Report. DL/ SCI/TM 36T Lewis GV, Catlow CRA (1985) Potential model for ionic oxides. J Phys C 18:1149-1161 McConnell JDC, Heine V (1985) Incommensurate structure and stability of Mullite. Phys Rev B 31:6140-6142 McConnell JDC, Heine V (1984) An aid to the structural analysis of incommensurate phases. Acta Crystallogr A 40:473-483 Padlewski S, Heine V, Price GD (1991) Atomic ordering around the oxygen vacany in sillimanite: A model for the mullite structure. Phys Chem Minerals (in press) Selke W (1988) The ANNNI model Theoretical analysis and experimental application. Phys Rep 170:214-264 Welberry TR (1985) Diffuse X-ray scattering and models of disorder. Rep Prog Phys 48:1543-1593 Yeomans J (1988) The theory and application of axial Ising models. Solid State Phys 41:151-200