- Email: [email protected]
Computer simulation of Ostwald ripening for ion beam synthesis of buried layers
Nuclear Instruments and Methods in Physics Research B 178 (2001) 135±137
www.elsevier.nl/locate/nimb
Computer simulation of Ostwald ripening for ion beam synthesis of buried layers T. Pohl
a,*
, C. Hammerl a, B. Rauschenbach b, U. R ude
c
a Experimentalphysik IV, Institut fur Physik, Universitatstr. 1, 86159 Augsburg, Germany Institut f ur Ober¯ achenmodi®zierung Leipzig und Institut fur Experimentelle Physik II, Universitat Leipzig, 04318 Leipzig, Germany c Lehrstuhl f ur Informatik 10 (Systemsimulation), Institut fur Informatik, Universitat Erlangen, Nurnberg, 91058 Erlangen, Germany
b
Abstract Ion beam synthesis of buried layers can be realized by a two-stage process of ion implantation and a post-implantation thermal treatment. During ion implantation precipitates are formed after exceeding the solid solution state and grow with increasing ¯uence. Thermal treatments can stimulate an onward growth of the implantation-induced precipitates up to the state of coalescence, where closed buried layers can be formed. In order to investigate such Ostwald ripening processes, computer simulations based on a newly developed model were carried out. The model consists of both the diusion behaviour and the precipitate evolution. Within one simulation time step the numerical calculation of diusion and precipitate changes is alternately performed until a stationary numerical solution is achieved. The simulation of the diusion part of the ripening process is realized by a multigrid algorithm, which leads to near real time calculations. The obtained simulation results are compared to the predictions of the theory of Lifshitz, Slyozov and Wagner (LSW) and the advanced model of Voorhees and Glicksman. Ó 2001 Published by Elsevier Science B.V. Keywords: Ostwald ripening; Diusion; Precipitates; Computer simulation; Ion beam synthesis
1. Introduction In the past several models for simulating the Ostwald ripening process were developed. As a main disadvantage the diusion was either totally neglected [1,2], only roughly approximated [3,4] or limited to one or two dimensions [5]. With the
increase in computing power and the development of highly ecient numerical methods it is now possible to simulate Ostwald ripening accurately considering the in¯uence of diusion in three dimensions and large simulation volumes.
2. Model equations * Corresponding author. Tel.: +49-821-598-3412; fax: +49821-598-3425. E-mail address: [email protected] (T. Pohl).
Virtually all models concerning the Ostwald ripening have the distinction of a stationary precipitated phase and a dissolved phase in common.
0168-583X/01/$ - see front matter Ó 2001 Published by Elsevier Science B.V. PII: S 0 1 6 8 - 5 8 3 X ( 0 1 ) 0 0 5 1 1 - 0
136
T. Pohl et al. / Nucl. Instr. and Meth. in Phys. Res. B 178 (2001) 135±137
In the presented model the latter is described by the diusion equation extended by a source term: oc r; t DDc r; t Q r; t; ot
1
where D is the diusion coecient of the solute, c r; t the concentration of solute atoms at the position r at time t and Q r; t the change in concentration due to absorption or emission of solute atoms by the precipitates. Since D is constant, in particular independent of the local concentration, interdiusion eects are neglected. The behaviour of the precipitates is described by a reaction equation based on the Gibbs± Thomson equation [6]: dNi t kNi t2=3 c ri ; t dt
cGT Ni ;
2
where Ni t speci®es the amount of particles in the precipitate indexed with i, k is a rate constant, c ri ; t is the concentration of solute atoms in the vicinity of the precipitate positioned at ri and cGT Ni is the Gibbs±Thomson concentration depending on the precipitate's size. It is assumed that the precipitates are spherical with a ®xed centre and that the lattice distortion energy can be ignored. To combine these two equations, it is necessary to replace Q r; t in Eq. (1) by Q r; t
P X i1
d jri
rj
dNi t : dt
The constant P is the amount of precipitates involved. The multiplication with the delta function leads to a localization of the absorbed respectively emitted particles. In general, the resulting equation system is far too complex to be solved analytically. So the implemented simulation software uses a multigrid algorithm for the diusion equation and a Runge± Kutta algorithm for the reaction equation to get a numerical solution (for details see [7]). 3. Comparison with the results of the LSW theory At ®rst, the in¯uence of the total precipitate volume fraction on the time-independent scaled
Fig. 1. Comparison of the time-independent radius distribution with the results of the LSW theory. The numbers show the total initial precipitate volume fraction.
radius distribution was studied and compared to the results of the LSW theory. The simulated system is a closed ensemble of SiO2 precipitates randomly positioned in Si bulk material at 1300°C. The simulation volume has a cubic shape with a length of 3.5 lm. As a starting radius distribution a Gauss function with a mean radius of 45 nm and a width of 3 nm was used. The three lines in Fig. 1 show the resulting timeindependent radius distribution for dierent initial total precipitate volume fractions. With increasing total precipitate fraction the distributions become broader and more symmetrical as is predicted by the model of Voorhees and Glicksman [4]. For R=R < 1 the simulation results agree very well with the LSW theory. The asymmetric peaks for lower volume fractions (2.1% and 0.5%) result from statistical ¯uctuations due to the smaller number of precipitates. To measure the accuracy of the simulation the initial amount of particles was compared to the amount after 3 104 s. The violation of mass conservation was about 1 ppm. 4. Simulation of the evolution of a buried layer In order to investigate the synthesis of a buried layer after ion implantation the evolution of precipitates during thermal treatment was simulated (see Table 1 for system parameters). The system
T. Pohl et al. / Nucl. Instr. and Meth. in Phys. Res. B 178 (2001) 135±137
137
Table 1 Simulation parameters for the MoO2 system Parameter
Description
Units
Value
Reference
r c1 D
Surface tension Equilibrium solubility Diusion coecient
Jm 2 cm 3 cm2 s
2.25 5:66 1018 0:03 exp 1:347 eV=kB T
[8] [9] [9]
consists of 2000 MoO2 precipitates in Mo bulk material at 500°C. The precipitates are randomly positioned up to a maximum depth of 500 nm. To model the Gaussian implantation pro®le the radius of the precipitates depend on their depth: 10 nm radius at the edge of the implantation pro®le; 30 nm in the middle. The solute atom concentration was homogeneously initialized with the equilibrium solubility c1 . In lateral direction Neumann boundary conditions are implemented which are comparable to quasi-periodic boundary conditions. The remaining two boundaries representing the interfaces to vacuum and to the bulk material, respectively, are ®xed to the equilibrium solubility by Dirichlet boundary conditions. Fig. 2 shows the precipitate volume fraction depending on the depth. The small precipitates at the edge of the implantation pro®le are completely dissolved which leads to an increase in solute atom concentration in these areas as can be seen in Fig. 3. csuper is the solute atom concentration exceeding the equilibrium solubility c1 . This screening of the inner precipitates against the lower edge concentrations results in the growing of the inner precipitates. A continuous evolution would coalesce precipitates in the centre of the implantation pro®le and a buried layer could be established.
Fig. 2. Evolution of the precipitate volume fraction.
1
Fig. 3. Evolution of the normalized solute atom concentration.
5. Conclusion Partially a very good agreement with results of the LSW theory and the model of Voorhees and Glicksman could be shown. The implementation of the new model is especially useful for the simulation of areas close to a surface typical for the ion beam synthesis scenario. References [1] I.M. Lifshitz, V.V. Slyozov, J. Phys. Chem. Solids 19 (1±2) (1961) 35. [2] C. Wagner, Z. Elektrochem. 65 (7±8) (1961) 581. [3] P.W. Voorhees, M.E. Glicksman, Acta Metall. 32 (11) (1984) 2001. [4] P.W. Voorhees, M.E. Glicksman, Acta Metall. 32 (11) (1984) 2013. [5] T. K upper, N. Masbaum, Acta Metall. 42 (6) (1994) 1847. [6] G. Martin, in: Solid State Phase Transformations in Metals and Alloys, Les editions de physique, 1978. [7] T. Pohl, U. R ude, B. Rauschenbach, to be published. [8] E. Lax, J. d'Ans, Taschenbuch f ur Chemiker und Physiker, Springer, Berlin, 1998. [9] H. Landolt, K.-H. Hellwege, B. Predel, Landolt±B ornstein, Numerical Data and Functional Relationships in Science and Technology, New series, Springer, Berlin.
www.elsevier.nl/locate/nimb
Computer simulation of Ostwald ripening for ion beam synthesis of buried layers T. Pohl
a,*
, C. Hammerl a, B. Rauschenbach b, U. R ude
c
a Experimentalphysik IV, Institut fur Physik, Universitatstr. 1, 86159 Augsburg, Germany Institut f ur Ober¯ achenmodi®zierung Leipzig und Institut fur Experimentelle Physik II, Universitat Leipzig, 04318 Leipzig, Germany c Lehrstuhl f ur Informatik 10 (Systemsimulation), Institut fur Informatik, Universitat Erlangen, Nurnberg, 91058 Erlangen, Germany
b
Abstract Ion beam synthesis of buried layers can be realized by a two-stage process of ion implantation and a post-implantation thermal treatment. During ion implantation precipitates are formed after exceeding the solid solution state and grow with increasing ¯uence. Thermal treatments can stimulate an onward growth of the implantation-induced precipitates up to the state of coalescence, where closed buried layers can be formed. In order to investigate such Ostwald ripening processes, computer simulations based on a newly developed model were carried out. The model consists of both the diusion behaviour and the precipitate evolution. Within one simulation time step the numerical calculation of diusion and precipitate changes is alternately performed until a stationary numerical solution is achieved. The simulation of the diusion part of the ripening process is realized by a multigrid algorithm, which leads to near real time calculations. The obtained simulation results are compared to the predictions of the theory of Lifshitz, Slyozov and Wagner (LSW) and the advanced model of Voorhees and Glicksman. Ó 2001 Published by Elsevier Science B.V. Keywords: Ostwald ripening; Diusion; Precipitates; Computer simulation; Ion beam synthesis
1. Introduction In the past several models for simulating the Ostwald ripening process were developed. As a main disadvantage the diusion was either totally neglected [1,2], only roughly approximated [3,4] or limited to one or two dimensions [5]. With the
increase in computing power and the development of highly ecient numerical methods it is now possible to simulate Ostwald ripening accurately considering the in¯uence of diusion in three dimensions and large simulation volumes.
2. Model equations * Corresponding author. Tel.: +49-821-598-3412; fax: +49821-598-3425. E-mail address: [email protected] (T. Pohl).
Virtually all models concerning the Ostwald ripening have the distinction of a stationary precipitated phase and a dissolved phase in common.
0168-583X/01/$ - see front matter Ó 2001 Published by Elsevier Science B.V. PII: S 0 1 6 8 - 5 8 3 X ( 0 1 ) 0 0 5 1 1 - 0
136
T. Pohl et al. / Nucl. Instr. and Meth. in Phys. Res. B 178 (2001) 135±137
In the presented model the latter is described by the diusion equation extended by a source term: oc r; t DDc r; t Q r; t; ot
1
where D is the diusion coecient of the solute, c r; t the concentration of solute atoms at the position r at time t and Q r; t the change in concentration due to absorption or emission of solute atoms by the precipitates. Since D is constant, in particular independent of the local concentration, interdiusion eects are neglected. The behaviour of the precipitates is described by a reaction equation based on the Gibbs± Thomson equation [6]: dNi t kNi t2=3 c ri ; t dt
cGT Ni ;
2
where Ni t speci®es the amount of particles in the precipitate indexed with i, k is a rate constant, c ri ; t is the concentration of solute atoms in the vicinity of the precipitate positioned at ri and cGT Ni is the Gibbs±Thomson concentration depending on the precipitate's size. It is assumed that the precipitates are spherical with a ®xed centre and that the lattice distortion energy can be ignored. To combine these two equations, it is necessary to replace Q r; t in Eq. (1) by Q r; t
P X i1
d jri
rj
dNi t : dt
The constant P is the amount of precipitates involved. The multiplication with the delta function leads to a localization of the absorbed respectively emitted particles. In general, the resulting equation system is far too complex to be solved analytically. So the implemented simulation software uses a multigrid algorithm for the diusion equation and a Runge± Kutta algorithm for the reaction equation to get a numerical solution (for details see [7]). 3. Comparison with the results of the LSW theory At ®rst, the in¯uence of the total precipitate volume fraction on the time-independent scaled
Fig. 1. Comparison of the time-independent radius distribution with the results of the LSW theory. The numbers show the total initial precipitate volume fraction.
radius distribution was studied and compared to the results of the LSW theory. The simulated system is a closed ensemble of SiO2 precipitates randomly positioned in Si bulk material at 1300°C. The simulation volume has a cubic shape with a length of 3.5 lm. As a starting radius distribution a Gauss function with a mean radius of 45 nm and a width of 3 nm was used. The three lines in Fig. 1 show the resulting timeindependent radius distribution for dierent initial total precipitate volume fractions. With increasing total precipitate fraction the distributions become broader and more symmetrical as is predicted by the model of Voorhees and Glicksman [4]. For R=R < 1 the simulation results agree very well with the LSW theory. The asymmetric peaks for lower volume fractions (2.1% and 0.5%) result from statistical ¯uctuations due to the smaller number of precipitates. To measure the accuracy of the simulation the initial amount of particles was compared to the amount after 3 104 s. The violation of mass conservation was about 1 ppm. 4. Simulation of the evolution of a buried layer In order to investigate the synthesis of a buried layer after ion implantation the evolution of precipitates during thermal treatment was simulated (see Table 1 for system parameters). The system
T. Pohl et al. / Nucl. Instr. and Meth. in Phys. Res. B 178 (2001) 135±137
137
Table 1 Simulation parameters for the MoO2 system Parameter
Description
Units
Value
Reference
r c1 D
Surface tension Equilibrium solubility Diusion coecient
Jm 2 cm 3 cm2 s
2.25 5:66 1018 0:03 exp 1:347 eV=kB T
[8] [9] [9]
consists of 2000 MoO2 precipitates in Mo bulk material at 500°C. The precipitates are randomly positioned up to a maximum depth of 500 nm. To model the Gaussian implantation pro®le the radius of the precipitates depend on their depth: 10 nm radius at the edge of the implantation pro®le; 30 nm in the middle. The solute atom concentration was homogeneously initialized with the equilibrium solubility c1 . In lateral direction Neumann boundary conditions are implemented which are comparable to quasi-periodic boundary conditions. The remaining two boundaries representing the interfaces to vacuum and to the bulk material, respectively, are ®xed to the equilibrium solubility by Dirichlet boundary conditions. Fig. 2 shows the precipitate volume fraction depending on the depth. The small precipitates at the edge of the implantation pro®le are completely dissolved which leads to an increase in solute atom concentration in these areas as can be seen in Fig. 3. csuper is the solute atom concentration exceeding the equilibrium solubility c1 . This screening of the inner precipitates against the lower edge concentrations results in the growing of the inner precipitates. A continuous evolution would coalesce precipitates in the centre of the implantation pro®le and a buried layer could be established.
Fig. 2. Evolution of the precipitate volume fraction.
1
Fig. 3. Evolution of the normalized solute atom concentration.
5. Conclusion Partially a very good agreement with results of the LSW theory and the model of Voorhees and Glicksman could be shown. The implementation of the new model is especially useful for the simulation of areas close to a surface typical for the ion beam synthesis scenario. References [1] I.M. Lifshitz, V.V. Slyozov, J. Phys. Chem. Solids 19 (1±2) (1961) 35. [2] C. Wagner, Z. Elektrochem. 65 (7±8) (1961) 581. [3] P.W. Voorhees, M.E. Glicksman, Acta Metall. 32 (11) (1984) 2001. [4] P.W. Voorhees, M.E. Glicksman, Acta Metall. 32 (11) (1984) 2013. [5] T. K upper, N. Masbaum, Acta Metall. 42 (6) (1994) 1847. [6] G. Martin, in: Solid State Phase Transformations in Metals and Alloys, Les editions de physique, 1978. [7] T. Pohl, U. R ude, B. Rauschenbach, to be published. [8] E. Lax, J. d'Ans, Taschenbuch f ur Chemiker und Physiker, Springer, Berlin, 1998. [9] H. Landolt, K.-H. Hellwege, B. Predel, Landolt±B ornstein, Numerical Data and Functional Relationships in Science and Technology, New series, Springer, Berlin.