Cation adatom diffusion on the NiO(0 0 1) surface by molecular dynamics simulation

Surface Science 486 (2001) 46±54

www.elsevier.nl/locate/susc

Cation adatom di€usion on the NiO(0 0 1) surface by molecular dynamics simulation T.E. Karakasidis a,*, D.G. Papageorgiou b, G.A. Evangelakis c a

Department of Civil Engineering, University of Thessaly, Pedion Areos, GR-38334 Volos, Greece b Department of Materials Science, University of Ioannina, 45110 Ioannina, Greece c Department of Physics, Solid State Division, University of Ioannina, 45110 Ioannina, Greece Received 25 September 2000; accepted for publication 29 March 2001

Abstract We present results concerning the di€usion of Ni2‡ adatom on the (0 0 1) surface of NiO obtained by molecular dynamics simulations based on a rigid ion potential model. A wide temperature region was covered ranging from 0.29Tm up to 0.85Tm , Tm being the melting point of the model system. Two possible adatom positions were found on the surface in accordance with static calculations. From the detailed analysis of the ionic trajectories it came out that the adatom di€uses on the surface via hopping and exchange mechanisms. Both processes exhibit Arrhenius behavior from where we deduced the corresponding migration energies. In addition, we found two distinct temperature regions re¯ecting di€erent energetic requirements for hopping di€usion. This is due to the spontaneous creation of anionic adatom that combine with the cationic adatom at high temperatures. Moreover, we found that although the frequency rates for hopping and exchange are comparable, the exchange mechanisms participate to the di€usion coecient more than the hopping process by as much as an order of magnitude at high temperatures. Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Molecular dynamics; Di€usion and migration; Surface di€usion; Nickel oxides; Single crystal surfaces; Adatom

1. Introduction Metal particles and adlayers on oxide surfaces are related to many technologically important applications such as metal±oxide contacts in microelectronic and photovoltaic devices, coatings for corrosion passivation, gas sensors and oxidesupported transition metal catalysts [1±3]. The knowledge of adatom equilibrium positions on the

* Corresponding author. Tel.: +30-421-74054; fax: +30-42174169/62660. E-mail address: [email protected] (T.E. Karakasidis).

surface as well as their di€usive behavior are quantities necessary for the understanding of the ®rst stages of interface formation or the development of thick metal layers on oxides [4,5]. Several experiments are devoted to the deposition of metals on metal oxides. Reviews of such experiments can be found in Refs. [1±3]. The majority of these experiments focus on the deposition of di€erent metals than the cation present in the oxide. In general, these studies concern the crystallographic relation between the metal adparticles and the substrate, while results referring to the adatom position at the interface or its di€usion behavior are rather limited. Theoretically, there

0039-6028/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 9 - 6 0 2 8 ( 0 1 ) 0 1 0 6 3 - 9

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

are ab initio studies concerning the structure and the binding energies of metal adsorbates on oxide surfaces [6±11]. From some of these studies it is concluded that the on-top O2 site is the most favorable for the cation adatom, followed by the hollow fourfold position in half distance between two O2 surface ions, while the cation on-top position exhibits greater energetic requirements. In addition, there are several static calculations in the aim of ®nding the cationic coverage dependence of some basic properties of the interface, looking at the heat of segregation [12±15]. Although static calculations are very useful, they are limited to zero temperature computations. Atomistic methods like molecular dynamics and Monte Carlo simulations are well suited for the study of adsorbate dynamics, since they take into account temperature e€ects and thus provide insight to the di€usion mechanisms and other dynamical properties. Although there are a lot of studies, mainly experimental, concerning the hetero-di€usion on metal oxides like MgO [16,17], CoO [18], NiO [19±22] or Al2 O3 [17,21] and Mg2 Al2 O4 [17], to our knowledge a detailed work on the di€usion of cation adatoms on oxide surfaces is rather lacking. Such a work could enlighten our understanding of growth, oxidation and other processes [4,5]. The aim of the present investigation is the study of the di€usive behavior of the cation adatom on the NiO(0 0 1) surface. This choice lies on the simplicity of this face and also on the availability of experimental [23±25] and simulation [26] results showing its stability. The paper is organized as follows: we start with the computational details before reporting results on the cationic adatom position on the surface. Then we proceed with the detailed presentation of the di€usion mechanisms observed, to end up with some concluding remarks. 2. Model and computational details The simulations were carried out using slab geometry in which the simulation cell was partially ®lled and included a region of vacuum. The simulation box containing the slab was a parallelepi-

47

ped with edges parallel to the x: [1 0 0], y: [0 1 0] and z: [0 0 1] directions. The slab consisted of 1728 ions arranged on 12 planes at each direction. In the z direction there was an empty space of four times the dimension of the slab. The use of periodic boundary conditions resulted in a system with two free (0 0 1) surfaces perpendicular to the z direction on the two opposite sides of a slab of in®nite extent in x and y directions. The simulations were carried out in the constant temperature canonical ensemble using the Nose scheme [27]. The equations of motion were integrated by means of the Verlet algorithm and a time step of 10 15 s. For the ionic interactions we adopted a rigid ion potential developed for NiO [28] that already gave satisfactory results to studies of structural [29,30] and di€usive properties [31] of P the 5(3 1 0)[0 0 1] grain boundary of NiO and also of the dynamical [32] and di€usive properties [33] of the NiO(0 0 1) surface. This potential has the following analytic form:   qa qb r Cab Vab ˆ ‡ Aab exp ; …1† qab r r6 the ®rst term refers to the Coulomb force; the second one represents a repulsive term, while the last one stands for the Van-der-Waals attractive contributions. Rigid ion potentials are limited by the fact that the electronic polarization is not taken fully into account, but their advantage is that they permit to perform quite long simulations. Adopting such a potential, we assume that NiO is completely ionic, the cations and anions having formal charges Ni2‡ and O2 , respectively. In addition, this model is derived from bulk properties; therefore its utilization in surface simulations makes the additional assumption that the modi®ed environment of ions on the surface does not a€ect the potential. In Eq. (1), qa;b stand for the charges of the ions a and b, r for their distance, while A, C and q are parameters that have been adjusted to experimental quantities. The values of the parameters are presented in Table 1. The short range interactions were truncated at a cuto€ distance 1.67a, a being the lattice parameter. The Coulombic contributions have been evaluated using the Ewald method [34]. Technical details

48

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

Table 1 Parameters of the short range rigid-ion potential for NiO [28] Term

Potential functional

Potential parameters

V‡‡ V‡ V

0 A‡ exp… r=q‡ † A exp… r=q †

V

Fifth-degree polynomial (1:5 < r < 2:75) Third-degree polynomial (2:75 < r < 2:9) C =r6 …2:9 < r)

± 556.43exp… r=0:3553† 22764.3exp… r=0:149† (r < 1:5) ±

V V

± ±

Energy is in electronvolts and length in Angstroms.

concerning Ewald parameters can be found in Ref. [33]. With this potential, the melting point for the bulk system was estimated at Tm ˆ 3500 K compared to the experimental value of 2257 K. The elevated melting point is a characteristic of rigid ion potentials. Others using a similar model for the NiO [26] found a melting point of the same magnitude. It is known that nickel oxide is slightly nonstoichiometric with a de®ciency in nickel and that the majority of the defects are cation vacancies. These vacancies are doubly or singly charged, their relative percentage being a function of the temperature and the partial oxygen pressure [35]. In the present study, we consider a stoichiometric surface and doubly charged Ni adatom. Of course, Ni adatoms can have other charges and can be also associated with charged vacancies or even more complex defects. Several such defects are known to exist for the case of the bulk from simulation results for the bulk and several grain boundaries [36±41]. Such systems require di€erent treatment permitting charge transfer between the ions and are dicult to be treated by classical molecular dynamics. However, although these simpli®cations may be important, if we are interested in bulk properties, they can serve as a rather reasonable approach when dealing with deposition of Ni on the NiO surface [26]. Since the spontaneous creation of an adatom is rather scarce within the time scale of our simulations, we created one by taking a cation from one surface and putting it on the other one. This corresponds to a cationic adatom concentration on the surface of about

1.4%, resulting in an over-saturated surface. Nevertheless, since we deal with single atom properties, we do not expect that the di€usion mechanisms and the corresponding migration energies would be a€ected. We equilibrated the system for about 30 ps, before starting the production runs, the duration of which varied between 245 and 500 ps, depending on the temperature. We have veri®ed that the presence of the adatom on one surface and the vacancy on the other do not a€ect the di€usion results by performing some runs with only one adatom on one surface. The frequency of a given di€usion mechanism is the total number of di€usion events, N, divided by the simulation length Trun and normalized by the average number of adatoms hNa i: Cˆ

1 N hNa i Trun

…2†

The di€usion coecient can, then, be obtained using the expression: 1 D ˆ Cd 2 4

…3†

where C is the frequency rate of the particular mechanism and d the corresponding distance the adatom covered.

3. Results 3.1. Adatom position 3.1.1. Zero temperature calculations Our ®rst concern was to ®nd the equilibrium adatom positions on the surface. These calculations were performed with one extra Ni2‡ adatom on one surface of the slab. We found two possible adatom positions. The lower energy position was found on top of a surface oxygen (Fig. 1(a) and (b)), in agreement with the literature [7±11]. It worth's noting, that the underlying oxygen ion is raised from the surface layer, at a distance nearly one third that of the nickel adatom from the surface, while the neighboring ions on the surface undergo slight distortions.

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

49

Fig. 1. Relaxed nickel adatom positions on the NiO(0 0 1) surface: (a) on-top site top view, (b) on-top site side view, (c) hollow site top view, (d) hollow site side view. Black and gray circles stand for cations and anions respectively, while white particle for the cation adatom.

The other adatom position, with energy of 0.58 eV higher than the on-top position, is the hollow site, located in the diagonal linking two oxygen ions of the surface layer. This result is also in accordance with the literature [7,9±11]. As we can see in Fig. 1(c) and (d) these oxygen ions do not occupy perfect surface positions, but are displaced in the in-plane directions and raised from the surface by about one third the adatom±surface distance, in contrast to the previous case (the on-top O site). 3.1.2. High temperature results Bearing these results in mind we performed a detailed adatom trajectory analysis at various temperatures. We found that the cationic adatom spends more time at the hollow sites than on the on-top sites and an important amount of time at transitional positions as it is di€using in a rather important way. This result, that is not expected by energetic considerations (static calculations showed that the on-top position is energetically favored), can be understood taking into account the dynamics of the system. At ®nite temperature, the atoms vibrate giving rise to conformations that facilitate the passage of the adatom from one

hollow site to a neighboring one: as it can be seen in Fig. 1(c), around one oxygen there are four equivalent hollow positions and only one on-top site. Passing from a hollow site to another hollow site is energetically easier than from an on-top position to another on-top site. In addition, we observed that close to the nickel adatom there are oxygen ions raised from the surface, a result consistent with static calculations and in agreement with the ®ndings of Kantorovich et al. [6]. The e€ect can be seen in Fig. 2, where we present a snapshot with the adatom viewed from the top and from the side at T ˆ 1000 K (0:29Tm ). At higher temperatures the situation is di€erent: an oxygen ion usually becomes adatom that combines with the nickel adatom. In Fig. 3, we present a side view of the surface at T ˆ 2800 K (0:80Tm ). We can see the nickel adatom and close to it the oxygen adatom forming a NiO admolecule. This di€erent behavior of the surface oxygen atoms is very important not only for the structural and vibrational properties of the surface, but also for the di€usive behavior of the cationic adatom. Indeed, we have to expect that the Ni adatom's di€usion must be in¯uenced by the presence of an

50

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

(a)

Fig. 4. Arrhenius diagram of the concentration of nickel ( ) and oxygen () adatoms.

(b)

Fig. 2. Snapshot showing a nickel adatom position at T ˆ 1000 K (black circles are nickel ions, white circle is the nickel adatom and gray circles are oxygen ions): (a) top view, (b) side view. Note the two surface oxygen ions close to the adatom that are raised from their normal surface positions.

Fig. 3. Surface layer side view showing a nickel adatom at T ˆ 2800 K. (Black circles are nickel ions and gray circles are oxygen ions.) Note the nickel and oxygen ions that are at adatom positions forming a NiO admolecule.

anion in its neighborhood. This e€ect is discussed in detail in the di€usion Section 3.3. The observed creation of oxygen adatoms could be related to segregation phenomena. 3.2. Adatom lifetime and concentration We calculated the concentration of Ni2‡ and O adatoms, Fig. 4. We can see that the Ni2‡ 2

concentration is virtually constant and equal to the value imposed by the simulation. However, the concentration of O2 adatom appears as being thermally activated with activation energy of 1.32 eV. As we can see, at high temperatures the concentration of O2 adatoms becomes comparable to that of the Ni2‡ adatom. Such a behavior could be explained by the formation of Ni2‡ and O2 pairs at higher temperatures. This ®nding suggests that the creation of oxygen adatoms is induced by the presence of nickel adatoms and can be related to segregation and wetting phenomena. We also calculated the average lifetime of Ni2‡ and O2 adatoms as a function of temperature and the obtained results are reported in Fig. 5. As we can see, at low temperatures up to 2000 K (0.57Tm ), the lifetime of Ni2‡ adatoms is very large, then as the temperature increases the lifetime decreases by as much as two order of magnitudes to reach a limiting value of about 10 ps. The lifetime of the O2 adatom is very small at low temperatures (the creation of oxygen adatoms is not very frequent), to increase practically to the same value as the Ni2‡ adatom, above 2000 K. These values have to be considered only as an indication since they could be a€ected by the ®xed cationic concentration used in our simulations. We have to mention here that the nickel adatom as well as the oxygen adatom are not always the same as there are exchanges that take place with surface ions.

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

Fig. 5. Lifetime of nickel ( ) and oxygen () adatoms as a function of temperature.

3.3. Di€usion mechanisms The detailed trajectory analysis revealed two basic cationic di€usion mechanisms: hopping and exchange. 3.3.1. Hopping di€usion mechanism As we have already mentioned in the previous section, the nickel adatom spends most of the time on hollow sites. In fact we observed that it circulates around the surface oxygen that is localized in the center of a cationic square. As it is already discussed, there are four possible hollow sites around the central oxygen and one on-top site. We found that the cationic adatom spends an important amount of time passing from one such hollow site to another hollow site and less on an on-top site. Therefore, we considered that the Ni2‡ adatom is localized in the cationic square. Consequently, a hopping event happens when the adatom passes from a cationic square to a neighboring one. We measured the corresponding jump frequencies and the results are presented in Fig. 6. As we can see, at about T ˆ 2000 K (0:57Tm ) there is a change of the activation energy, from 0.09 eV at low temperatures, to 0.62 eV at higher temperatures. This change in the energetic requirements of the same di€usion mechanism can be explained considering that at high temperatures a surface oxygen ion becomes adatom and combines with the cationic adatom. The presence of the oxygen

51

Fig. 6. Arrhenius diagrams for hopping (j, ) and exchange ( ) frequencies of nickel adatom. Black and white squares refer to the low and high temperature hopping frequency rates.

adatom close to the nickel adatom apparently slows down the di€usion of the nickel adatom due to the extra electrostatic lateral interactions. The low temperature behavior is consistent with the vibrational results for the cation adatom on a (0 0 1) surface obtained by molecular dynamics simulation [32] at room temperature. From that study it was found that the cationic adatom is responsible for a phonon mode at low frequencies in the in-plane direction and another one at high energy in the perpendicular to the surface direction. 3.3.2. Exchange mechanism In this di€usion mechanism, besides the adatom one or more surface cations participate also. In the simple exchange mechanism the adatom goes down to the surface on the position of surface cation which in his turn becomes adatom. In the case of the double exchange process two surface cations are involved. A representative example of double exchange event is presented in Fig. 7 in a succession of snapshots. The nickel adatom, labeled A, goes down to the surface to the position of a surface cation B. Atom B moves simultaneously towards the position of another surface cation on position C, that is ®nally popped out in adatom position. We found that the double exchange events appear to be more frequent than the simple exchanges although their frequencies are of the same order of magnitude. At higher temperatures multiple exchanges incorporating more than

52

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

two surface ions have been observed, but their frequency was much smaller than that of simple and double exchanges. The total number of exchanges (regardless the number of ions participating in the mechanism) was counted for all temperatures and the results are presented in the Arrhenius diagram of Fig. 6, from where we deduce an apparent activation energy of 1.1 eV. As we can see in this ®gure, at low temperatures the hopping mechanism is the dominant di€usion process, while at higher temperatures the frequencies of the two mechanisms become comparable with the exchange mechanism being slightly more important. However, the contribution of the mechanisms observed to the di€usion coecient is not the same as the di€usion lengths are di€erent. The average length contribution pto the di€usion coecient for p hopping is about a 2=2, for simple p exchange a 2=2 and for double exchange a 2, a being the lattice constant. Using Eq. (3) we calculated the di€usion coecients for hopping and exchange respectively. Given that the di€usion length of the double exchange is twice as large as that of the simple exchange, its contribution to the di€usion coecient will be four times larger than that of simple exchange. Consequently, the double exchange mechanism is the dominant exchange type di€usion process and it has been considered for the estimation of the exchange di€usion coef®cient. The obtained results for hopping and exStarting Configuration

change are presented in Fig. 8. As we can see hopping di€usion is dominant at low temperatures (T < 0:57Tm ), while the exchange mechanism becomes more important than hopping at higher temperatures. Similar behavior is found in hetero-di€usion experiments. To our knowledge there are no direct measurements of cation adatom di€usion on the NiO(0 0 1) surface as a function of temperature. However, from measurements concerning the

Fig. 8. Arrhenius diagrams for the di€usion coecients of the cationic adatom. Squares and circles stand for exchange and hopping di€usion events, while black and white square denote the low and high temperature results.

After 4 ps

C C

A B

A

B

After 5.5ps

After 7.5 ps

C

C

B

B

A

A

Fig. 7. Double exchange mechanism at T ˆ 2800 K. The nickel adatom (labeled A) takes the position of the nickel surface ion (labeled B) that replaces the surface nickel ion labeled C. Finally, the C particle is popped out in adatom position.

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

di€usion of Cr3‡ and Co2‡ carried on MgO(0 0 1) surface [22], that has similar structure with the NiO(0 0 1) surface, a change of the di€usion activation energy was observed above a certain temperature attributed to another mechanism becoming dominant above this temperature. Such a change in the slope was not observed for the NiO(1 1 0) surface [21,22] at least for the temperature range studied. However, the NiO(1 1 0) surface exhibits a vicinal character and therefore a direct comparison with the NiO(0 0 1) face is not straightforward. It should be emphasized that the important contribution of the exchange di€usion mechanism, if present also in the cases of hetero-deposition, could play an important role to segregation and wetting of oxide surfaces. 4. Conclusions In this paper we present results obtained by molecular dynamics simulation in a wide temperature range (0.29Tm ±0.85Tm ) concerning the diffusion of Ni2‡ adatom on a NiO(0 0 1) surface. We found that the adatom occupies two di€erent positions, one on top of a surface oxygen ion and another one in the diagonal linking two nearest surface oxygen ions, with a marked preference for the latter. The oxygen ions that are nearest neighbors of the cation adatom are raised from the surface by about one third the adatom±surface distance. As temperature rises there is an increased possibility that one of these oxygen ions becomes adatom. Di€usion of the cationic adatom takes place through two mechanisms, simple hopping and double exchange. Hopping is dominant at low temperatures, while exchange becomes important at higher temperatures. Hopping requires activation energy of 0.09 eV for temperatures up to 0:57Tm and 0.62 eV for higher temperatures. This increase in the activation energy is due to the presence of oxygen adatoms close to the nickel adatom at high temperatures. The exchange type di€usion mechanisms require activation energy of 1.1 eV. In addition, although the frequency rates for hopping and exchange are comparable, the exchange mechanisms contribute to the di€usion

53

coecient more than the hopping process by as much as an order of magnitude at high temperatures.

References [1] C.T. Campbell, Surf. Sci. Rep. 27 (1997) 1. [2] C.R. Henry, Surf. Sci. Rep. 31 (1998) 235. [3] V.E. Henrich, P.A. Cox, The Surface Science of Metal Oxides, Cambridge University Press, 1994. [4] A. Atkinson, Rev. Modern Phys. 57 (1985) 437. [5] S. Mitra, Mater. Sci. Engng A 174 (1994) 103. [6] L. Kantorovich, A. Shluger, J. Gunster, J. Stultz, S. Krischok, D.W. Goodman, P. Stracke, V. Kempter, Nuclear Instrum. Meth. Phys. Res. B 157 (1999) 162. [7] J. Goniakowski, Phys. Rev. B 57 (1998) 1935. [8] R. Wu, A.J. Freeman, Phys. Rev. B 51 (1995) 5408. [9] Y. Li, D.C. Langreth, M.R. Pederson, Phys. Rev. B 52 (1995) 6067. [10] J.R. Smith, T. Hong, D.J. Srolovitz, Phys. Rev. Lett. 72 (1994) 4021. [11] C. Li, R. Wu, A.J. Freeman, C.L. Fu, Phys. Rev. B 48 (1993) 8317. [12] E.A. Colburn, W.C. Mackrodt, P.W. Tasker, in: G. Simkovich, V.S. Stubican (Eds.), Transport in Nostoichiometric Compounds, NATO ASI B129, 1984, p. 397. [13] P.W. Tasker, E.A. Colbourn, W.C. Mackrodt, J. Am. Ceram. Soc. 68 (1985) 74. [14] P. Masri, P.W. Tasker, J.P. Hoare, J.H. Harding, Surf. Sci. 173 (1986) 439. [15] W.C. Mackrodt, P.W. Tasker, J. Am. Ceram. Soc. 72 (1989) 1576. [16] C.M. Lin, V.S. Stubican, J. Am. Ceram. Soc. 70 (1987) C73. [17] V.S. Stubican, G. Huzinec, D. Damjanovic, J. Am. Ceram. Soc. 68 (1985) 181. [18] P.S. Maiya, W.K. Chen, N.L. Peterson, Metal. Trans. 1 (1970) 801. [19] C.M. Lin, V.S. Stubican, J. Am. Ceram. Soc. 73 (1990) 587. [20] G. Dhalenne, A. Revcolevschi, C. Monty, in: G. Simkovich, V.S. Stubican (Eds.), NATO ASI B129, 1985, p. 371. [21] E.P. Macey, V.S. Stubican, J. Am. Ceram. Soc. 76 (1993) 557. [22] V.S. Stubican, Phil. Mag. A 68 (1993) 809. [23] M.R. Castell, P.L. Wincott, N.G. Condon, C. Muggelberg, G. Thornton, S.L. Dudarev, A.P. Sutton, G.A.D. Briggs, Phys. Rev. B 55 (1997) 7859. [24] M.R. Castell, S.L. Dudarev, G.A.D. Briggs, A.P. Sutton, Phys. Rev. B 59 (1999) 7342. [25] N. Floquet, L.-C. Dufour, Surf. Sci. 126 (1983) 543. [26] P.M. Oliver, G.W. Watson, S.C. Parker, Phys. Rev. B 52 (1995) 5323. [27] S. Nose, J. Chem. Phys. 81 (1984) 511.

54

T.E. Karakasidis et al. / Surface Science 486 (2001) 46±54

[28] C. Massobrio, M. Meyer, J. Phys.: Condens. Matter. 3 (1991) 279. [29] M. Meyer, C. Waldburger, Mater. Sci. Forum 126±128 (1993) 229. [30] T. Karakasidis, M. Meyer, Modelling Simulation Mater. Sci. Engng. 8 (2000) 117. [31] T. Karakasidis, M. Meyer, Phys. Rev. B 55 (1997) 13853. [32] T. Karakasidis, D.G. Papageorgiou, G.A. Evangelakis, Appl. Surf. Sci. 162±163 (2000) 233. [33] T. Karakasidis, G.A. Evangelakis, Surf. Sci. 436 (1999) 193.

[34] P.P. Ewald, Ann. der. Physik 64 (1921) 253. [35] R. Farhi, G. Petot-Ervas, J. Phys. Chem. Solids 39 (1978) 1175. [36] C.R.A. Catlow, W.C. Mackrodt, M.J. Norgett, A.M. Stoneham, Phil. Mag. 35 (1977) 177. [37] W.C. Mackrodt, Solid State Ionics 12 (1984) 175. [38] D.M. Du€y, P.W. Tasker, Phil. Mag. 50 (1984) 143. [39] D.M. Du€y, P.W. Tasker, Physica 131B (1985) 46. [40] D.M. Du€y, P.W. Tasker, Phil. Mag. 54 (1986) 759. [41] J. Rabier, J. Soullard, M.P. Puls, Phil. Mag. 61 (1990) 99.