Monte Carlo investigation of the thermal stability of coherent multilayers

Scripta mater. 43 (2000) 887– 892 www.elsevier.com/locate/scriptamat

MONTE CARLO INVESTIGATION OF THE THERMAL STABILITY OF COHERENT MULTILAYERS A. Ullrich, M. Bobeth and W. Pompe Institut fu¨r Werkstoffwissenschaft, Technische Universita¨t Dresden, Hallwachsstrasse 3, D-01062 Dresden, Germany (Received June 8, 2000) (Accepted June 20, 2000) Keywords: Multilayers; Alloys; Annealing; Thermal stability; Computer simulation

Introduction Nanoscale metallic multilayers possess outstanding physical properties for various promising applications. In order to optimize these properties, as-deposited multilayers are often thermally treated at elevated temperatures. A thermal loading of multilayers can also occur under processing and operating conditions. The present paper addresses the question of the thermal stability of A/B layer stacks with thin individual layers of only a few atomic layers. The components A and B are supposed to be immiscible. Nanoscale Co/Cu multilayers, exhibiting the giant magnetoresistance, are a typical example. Because of the small lattice misfit of about 2%, the fcc-Co layers and the Cu layers are coherent. As a first step to the understanding of the thermal decomposition of such multilayer films, we considered a coherent A/B layer stack with B-layer thickness of 4 to 6 monolayers (ML) and considerably thicker A-layers. To investigate the thermally activated rearrangements of atoms in this structure, Monte Carlo (MC) simulations have been performed. As initial state, an A/B layer stack with ideally planar A/B interfaces was assumed. This assumption has most relevance to single crystal superlattices, but also to polycrystalline multilayers with lateral grain size which is large compared to the multilayer period. The aim of the simulations was to elucidate the initiation mechanism of the thermally induced layer decomposition as well as the subsequent development of the phase morphology.

Monte Carlo Algorithm The present simulations were performed on a fcc lattice occupied by A- and B-atoms as well as one vacancy. The chemical binding between atoms has been described by nearest neighbor pair interactions ⑀AA, ⑀AB and ⑀BB. As is well known, such a model of a binary alloy is equivalent to the Ising model. Phase separation at low temperatures occurs if ⑀: ⫽ ⑀AA ⫹ ⑀BB ⫺ 2⑀AB ⬍ 0, where the critical temperature for the fcc lattice structure is given by kBTc/兩⑀兩 ⫽ 2.45 [1]. To study the thermal rearrangements of A- and B-atoms, the migration of one vacancy was simulated by means of a residence-time algorithm (see e.g. Ref. [2,3]). As a so-called rejection-free algorithm [4], this method exhibits a high computational efficiency also at lower temperatures compared to the conventional Metropolis algorithm [5]. One important issue of this approach is the calculation of the 1359-6462/00/$–see front matter. © 2000 Acta Metallurgica Inc. Published by Elsevier Science Ltd. All rights reserved. PII: S1359-6462(00)00507-8

888

MONTE CARLO INVESTIGATION OF COHERENT MULTILAYERS

Vol. 43, No. 10

transition rates of atoms (adjacent to the vacancy) to the vacant site. The jump rate of the i th atom to the vacancy at lattice site ␣ is described by

␯ ␣i ⫽ ƒ ␣i exp共⫺E ␣a i/k BT兲,

i ⫽ 1, 2,· · ·z

(1)

where ƒ␣i is the attempt frequency and Ea␣i is the activation energy for the atom jump. z is the number of nearest neighbors of the vacancy. The activation energy is given by the difference, E␣a i ⫽ E␣S i ⫺ E1␣i, of the binding energy of the ith atom, ES␣i, at the saddle point position of the jump and of the binding energy before the jump, E1␣i. The dependence of the attempt frequency on the kind of the jumping atom and the atom configuration around the vacancy has been neglected in the present study, i. e. ƒ␣i ⫽ ƒ0 ⫽ const. For the estimation of the activation energy, different approximations have been proposed in the literature. In the works by Martin and coworkers (see e.g. [2]), the saddle point energy E␣S i has been assumed to be configuration independent, referred to as ‘constant saddle point approximation’ (CSPA) in the following. This means that E␣S i ⫽ ES ⫽ const if also the dependence on the kind of the jumping atom is ignored. The value of ES essentially determines the waiting times between vacancy jumps, but it is not important for selecting the jump directions. Thus, we split the activation energy E␣a i ⫽ EaA ⫹ E˜␣a i into a configuration-independent part EaA and a second one E˜a␣i which is normalized so that E˜␣a i ⫽ 0 for the case of pure metal A. Then, the transition rate can be written as

␯ ␣i ⫽ ␯ A exp共⫺E˜ ␣a i/k BT兲

(2)

where ␯A ⫽ ƒ0 exp(⫺EaA/kBT) represents the transition rate of one atom to the vacant site in metal A. From the binding energy of the ith atom before it jumps, E1␣i ⫽ niX⑀XX ⫹ (z ⫺ 1 ⫺ niX)⑀AB, the energy E˜a␣i results as 1 E˜ ␣a i ⫽ 关共 z ⫺ 1 ⫺ n iX兲共 ⑀ ⫹ u兲 ⫹ 共1 ⫹ ⌰ X兲n iXu兴, 2

X ⫽ A, B

(3)

where niX is the number of nearest neighbor atoms of the ith atom of the same type X, u ⫽ ⑀AA ⫺ ⑀BB, ⌰A ⫽ ⫺1 and ⌰B ⫽ 1. To capture roughly the dependence of the saddle point energy E␣S i on the local atomic configuration, Rautiainen and Sutton [3] used the approximation ES␣i ⫽ E0 ⫹ (E1␣a ⫹ E␣2 i)/2 where E␣1 i and E␣2 i are the binding energies before and after the jump of the ith atom and E0 is a configuration-independent part. This approach is referred to as ‘mean saddle point approximation’ (MSPA) in the following. Analogously to (3), one obtains 1 E˜ ␣a i ⫽ 共n aX ⫺ n iX兲共 ⑀ ⫺ ⌰ X u兲 4

(4)

where n␣X is the number of nearest neighbor atoms of the ith atom after the jump to site ␣ which are of the same type X. The mean jump frequency ␯␣ of the vacancy at site ␣ is determined by the transition rates, ␯␣i, for z

冘 ␯␣i. Correspondingly, the mean the jumps of the surrounding atoms to the vacant site: ␯␣ ⫽ i⫽1 residence time of the vacancy at lattice site ␣ is given by ␶␣ ⫽ 1/␯␣. The lattice site k, to which the k⫺1

k

冘 ␯␣i ⬍ r␯␣ ⱕ i⫽1 冘 ␯␣i , k ⫽ 1, 2, · · · z, where vacancy jumps, is selected according to the condition i⫽0 ␯␣0 ⫽ 0 and r is a uniform random number between 0 and 1 [2,3]. The residence times of the vacancy at sites ␣ are random quantities, which can be generated according to ␶ran ␣ ⫽ ⫺␶␣ with r as random

Vol. 43, No. 10

MONTE CARLO INVESTIGATION OF COHERENT MULTILAYERS

889

Figure 1. Hole in a 4 ML-thick B-layer observed in a simulation.

sim number [2]. The sum of all residence times ␶ran ⫽ 冘␣ ␶␣ran . To save ␣ represents the simulation time t computational time, we added in our simulations simply the mean values, i. e. tsim ⫽ 冘␣ ␶␣ .

Simulation Results and Discussion The atomic rearrangements during annealing were simulated within a rhombohedral sample consisting of L3 lattice sites (L ⫽ 102) and periodic boundary conditions were applied. To visualize the phase morphology in the following figures, the atom distribution was periodically continued to fill a cube enclosing the rhomboeder. The multilayer period was fixed to 102 monolayers in ⬍111⬎ direction and the thickness of the B-layers was varied from 4 to 6 monolayers. For simplicity, we considered the symmetric case ⑀AA ⫽ ⑀BB. The simulations were performed at a temperature of T ⫽ ⫺1.08 ⑀/kB, corresponding for example to T ⫽ 691 K for a critical phase separation temperature of Tc ⫽ 1573 K (T/Tc ⫽ 0.44), which roughly applies to the Co-Cu system. During the simulation, the normalized time ˜tsim ⫽ ␯A 冘 ␶␣ was calculated. This time has been related ␣

roughly to the real time by employing data on the formation and migration enthalpies of vacancies in pure metal A. This is of course a crude approximation, since the formation enthalpy depends actually on the composition. As an example, copper has been chosen as metal A. The vacancy concentration is given by c␯ with the excess entropy of mixing Sv␯F ⬇ 2 kB and the formation enthalpy HF␯ ⫽ 1.28 eV 3 [6]. Since the vacancy concentration in our simulations, csim ␯ ⫽ 1/L , was artificially enhanced, the real sim sim sim sim time results by the scaling t ⫽ (c␯ /c␯)t ⫽ (c␯ /c␯ ␯A)t˜ . The atomic jump rate ␯A was estimated M by ␯A ⫽ ƒ0 exp(⫺HM ␯ /kBT) where H␯ is the migration enthalpy, and the effective attempt frequency ƒ0 is of the order of three times the Debye frequency [6]. We used here ƒ0 ⫽ 2 ⫻ 1013 Hz and HM ␯ ⫽ 0.7 eV. Similar to classical MC studies, the simulation time can be characterized also by the number of Monte Carlo steps (MCS) where one MC step is defined as the total number of vacancy jumps divided by the number of lattice sites L3. A series of simulations revealed characteristic features of the morphology development of the layer stack starting with the formation of breakthroughs of A-atoms through B-layers, which have been called ‘holes’ in B-layers (Fig. 1). The atomic configurations in Figs. 1 to 3 are represented by plotting the B-atoms only. Figs. 2 and 3 show snapshots of the decomposition of a B-layer observed in a simulation

890

MONTE CARLO INVESTIGATION OF COHERENT MULTILAYERS

Vol. 43, No. 10

Figure 2. Stages of layer decomposition: (a—118 MCS) hole formation and growth, (b—790 MCS) fiber-like network, and (c—15160 MCS) disperse morphology (top views of 4 ML-thick slices of a multilayer including one original 4 ML-thick B-layer, side length of hexagon 72 a0 with a0 as fcc lattice constant).

using the CSPA (3). The observations in the case of the MSPA (4) are similar. After hole formation, layer decomposition proceeds by the lateral extension and coalescence of holes, resulting in a fiber-like network of phase B embedded in phase A (Fig. 2b). In the course of time, the mesh width of the network enlarges and the fibers decompose due to the Rayleigh instability so that finally a disperse morphology develops (Fig. 2c). For a quantitative characterization of the thermal stability of multilayers, we analyzed the time ts (time of layer stability) when the first hole appears. To detect the first hole during the simulations, a cylindrical volume with height dB and diameter dB was laterally scanned from lattice site to lattice site through the original B-layer after each MC step. During this scan, the concentration of A-atoms within

Figure 3. Cross section view of the decomposition of the B-layer in Fig. 2. The drawings show 1 ML-thick slices of one B-layer at different simulation times in MC steps. The development up to 118 MCS suggests hole formation to be caused by interface roughening.

Vol. 43, No. 10

MONTE CARLO INVESTIGATION OF COHERENT MULTILAYERS

891

Figure 4. Logarithmic plot of the time of the appearance of the first hole versus B-layer thickness in terms of the real time (right) and in terms of MC steps (left). The error bars at the symbols represent the standard deviations of the mean values resulting from 16 MC runs (triangles—CSPA, circles—MSPA).

this volume was calculated. The A-concentration threshold above which a hole is recognized was chosen as chole ⫽ 0.9. A logarithmic plot of the resulting time of layer stability as a function of the A B-layer thickness is shown in Fig. 4. For comparison, the beginning of hole formation is also given in terms of MC steps. Both plots in Fig. 4 reveal significant differences in the results obtained from approximations (3) and (4) for the activation energy of vacancy jumps. The MSPA (4) leads to an almost constant scaling factor between time and MC steps of 0.17 seconds/MCS, whereas for the CSPA it is considerably smaller and varies slightly from 0.079 at 4 ML to 0.062 seconds/MCS at 6 ML. The dependence of ts on the layer thickness is similar within the two approximations. The almost linear plot of log(ts) versus layer thickness for the MSPA suggests an exponential dependence of ts on the layer thickness. The CSPA seems to deviate somewhat from this dependence. The observed value for the time of layer stability ts depends also on the area of the considered layer, it decreases with increasing area. The analysis of this dependence would require relatively high computational efforts. Thus, a deeper understanding of hole formation within the framework of analytical models is desirable. The basic process determining hole formation is presumably the roughening of A/B interfaces due to thermal fluctuations. According to Fig. 3, holes are obviously formed by the contact of randomly appearing depressions of the upper and lower A/B interfaces of one B-layer. Starting from planar interfaces as non-equilibrium initial state, an equilibrium roughness develops during a certain time depending on the sample size. Interface fluctuations at thermal equilibrium have been studied comprehensively (see e.g. [7]). According to these studies, surface roughening occurs at temperatures above the so-called roughening temperature TR ⫽ 2␥˜ a2/␲kB where ␥˜ is the surface stiffness and a is the lattice spacing along the direction perpendicular to the surface. Estimating ␥˜ by the A/B interface energy in ⬍111⬎ direction, ␥111 ⫽ 2公3兩⑀兩/a20 ⫽ 0.236 J/m2, and the spacing a by the lattice plane spacing d111 ⫽ 0.208 nm, one finds TR ⫽ 473 K. Thus, our simulations at 691 K were performed above the roughening temperature, and the following estimation should be appropriate. Following [7], the rough A/B interface is described by the function z(x, y) with vanishing mean value (具z典 ⫽ 0). The standard deviation of the roughness in the small slope approximation (兩ⵜz(x, y)兩 ⬍⬍ 1) is then given by 公具z2典 ⫽ [kBT ln(l⌳)/2␲␥˜ ]1/2 where l2 is the considered interface area and ⌳ is an upper cut off in the wave vectors of interface fluctuations. For a rough estimate, ⌳ ⬀ 1/a0 ⬇ 2.8 nm⫺1 has been chosen. With the further parameters T ⫽ 691 K, ␥˜ ⬇ 0.2 J/m2 and the macroscopic length l ⫽ 0.01 m, one obtains a standard deviation of about 0.4 nm and, for the present simulations with l of the order of 100a0, it is about 0.2 nm. This value of the order of 1 ML supports our presumption that hole formation is caused by interface roughening. The probability for the appearance of sufficiently deep

892

MONTE CARLO INVESTIGATION OF COHERENT MULTILAYERS

Vol. 43, No. 10

interface depressions forming holes in the layer increases with annealing time. To the authors’ knowledge, an analytical derivation of this probability as a function of time and sample size has not been given yet. In view of the large influence of crystal anisotropy on layer stability, this problem is possibly beyond the scope of continuum theories and has to be tackled then by extended MC investigations. Conclusions The present MC simulations reveal that the thermal decomposition of coherent A/B multilayers starts by local breakthroughs of component A through the thinner B-layers. The appearance of breakthroughs is suggested to be caused by thermal roughening of A/B interfaces. Simultaneous lateral growth of breakthroughs leads to the formation of a fiber-like morphology of component B embedded in component A which transforms to a disperse morphology in a later stage. The thermal stability of individual layers has been characterized by the time when the first layer breakthrough occurs. This time strongly depends on the layer thickness, similar to an exponential dependence within the considered range from 4 to 6 monolayers. Acknowledgments This work was supported by the Deutsche Forschungsgemeinschaft, SFB 422. The simulations were performed at the Zentrum fu¨r Hochleistungsrechnen, Technische Universita¨t Dresden. References 1. 2. 3. 4. 5.

G. A. Baker, Phys. Rev. 124, 768 (1961). T. Abinandanan, F. Haider, and G. Martin, Acta Mater. 46, 4243 (1998). T. T. Rautiainen and A. P. Sutton, Phys. Rev. B59, 13681 (1999). A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comp. Phys. 17, 10 (1975). K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical Physics—An Introduction, Springer-Verlag, Berlin (1992). 6. Landolt-Bo¨rnstein, in Numerical Data and Functional Relationships in Science and Technology, vol. 25, Atomic Defects in Metals, ed. H. Ullmaier, Springer-Verlag, Berlin (1991). 7. P. Nozieres, in Solids Far from Equilibrium, ed. C. Godreche, p. 111, Cambridge University Press, Cambridge (1992).