Molecular Dynamics Study of Explosive Crystallization of SiGe and

Erik J. Albenze, Michael O. Thompson, and Paulette Clancy ... Joseph M. Zaug , Sorin Bastea , Jonathan C. Crowhurst , Michael R. Armstrong , and Nick ...
1 downloads 0 Views 627KB Size
Download PDF
5628

Ind. Eng. Chem. Res. 2006, 45, 5628-5639

Molecular Dynamics Study of Explosive Crystallization of SiGe and Boron-Doped SiGe Alloys Erik J. Albenze,† Michael O. Thompson,‡ and Paulette Clancy† Department of Chemical and Biomolecular Engineering and Department of Materials Science and Engineering, Cornell UniVersity, Ithaca, New York 14853

Explosive crystallization of SiGe alloys and boron-doped SiGe alloys is studied using a special molecular dynamics construction developed by the authors. For the case of Si1-xGex alloys, the simulations indicate that explosive crystallization should occur over the entire compositional range of Si1-xGex alloys at velocities high enough to prevent Ge segregation. The results show qualitative agreement with experimental results for Si-rich systems: Over a range of 0-25% germanium, the simulations predicted a nonlinear decrease in velocity that roughly matches the slope and extent of the rapid drop in interface velocity observed experimentally. An analogous investigation of boron-doped Si1-xGex alloys showed that the limit of incorporation of boron during explosive crystallization is strongly dependent on the amount of germanium present and that explosive crystallization is capable of limiting transient enhanced diffusion throughout the crystallization process. Introduction Increasing performance demands from electronic devices dictates that the size of the device features must continually shrink,1 in turn requiring that source and drain regions are increasingly heavily doped. This constitutes an obstacle due to the difficulty in maintaining control over the location of dopant atoms within specified junction depths, e.g., to prevent “channeling” following implantation of dopant atoms. Additionally, the ion implantation process by which the samples are typically doped causes damage (including amorphization) to the crystalline matrix which must be annealed by rapid heating to restore the crystallinity of the sample. In the case of p-type dopants such as boron, the annealing process facilitates unwanted boron diffusion in a process known as transient enhanced diffusion (TED),2-7 which can also enhance clustering of the boron atoms. This creates an additional problem since clusters of boron atoms are not electronically active. One method which has been shown to decrease TED is to alloy Si with Ge which has been shown to reduce the diffusivity of the B atoms.8,9 As a result, it would be desirable to have a better understanding of the crystallization processes of heavily doped silicon-germanium alloys. All of the concerns mentioned above have their roots in atomic-level diffusional processes, lending impetus to a computational study of fast crystallization processes with atomic-level detail. This is the main aim of this paper. Silicon and germanium are capable of forming two distinct and well-characterized solid phases, one crystalline and one amorphous. Both of these phases have local tetrahedral coordination. While the crystalline phase possesses a diamond structure, the amorphous phase has little or no long-ranged order. The metallic liquid phase of silicon has a coordination number of approximately 6.4, defined as the average number of nearest neighbors per atom; in contrast, both solid phases are semiconductors with an average coordination number around 4.0. The melting points of the crystalline phases of Si and Ge are 1685 and 1210 K, respectively. Melting points for the amorphous phases are considerably lower, 1450 K for Si and 975 K for * To whom correspondence should be sent. Tel.: (607) 255-6331. Fax: (607) 255-9166. E-mail: [email protected]. † Department of Chemical Biomolecular Engineering. ‡ Department of Materials Science and Engineering.

Ge.10-12 Thus the melting points of the two phases for both Si and Ge differ by over 200 K. Importantly, this difference in melting points allows a liquid phase to exist at a temperature that is higher than the melting point of the amorphous material but lower than that of the crystalline phase. This situation allows silicon and germanium to belong to a class of materials which are capable of exhibiting a phenomenon known as explosiVe crystallization in which an amorphous phase is converted to crystalline material at high speed (15 m/s for silicon).10,13-20 For Si and Ge, it is now generally accepted that this self-sustaining process occurs through movement of two solid/liquid interfaces separated by a thin liquid layer10,13,14,20 (for further discussion see refs 21 and 22). Driving this reaction are differences in the melting temperatures and enthalpies of the two phases. Once an initial liquid layer has been created (e.g., by the action of a laser), crystallization occurs rapidly. In the case of silicon, atoms crystallizing at the crystal/liquid interface release heat which is thermally conducted through the liquid to the liquid/amorphous boundary. The subsequent melting of amorphous Si withdraws heat from the system, encouraging further crystallization, causing the process to autopropagate. This is essentially a chain reaction process and, experimentally, it self-selects the fastest direction in which to propagate (for silicon, this is the [100] direction). The overall speed of explosive crystallization is dictated by the kinetics of liquid phase crystallization as expressed in the interface response function, briefly described below. Interface response functions describe the relationship between the solid/liquid interface temperature and the solidification velocity of the moving interface. When a crystal seed is in contact with a liquid phase and the temperature is less than the crystal melting point, the crystal/liquid interface moves such that the liquid phase is converted into crystalline material. The speed at which this interface moves will be dictated by the temperature at the interface. At 0 K, there is no energy to allow the interface to move. At the crystal melting point (T ) 1685 K for Si and 1210 K for Ge), the interface will also not move since there is no thermodynamic driving force. At some temperature between the melting point of the crystal and 0 K, the velocity of the moving solid/liquid interface must pass through a maximum. This relation between the temperature and the interface velocity is known as the interface response function

10.1021/ie051361w CCC: $33.50 © 2006 American Chemical Society Published on Web 03/02/2006

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5629

(IRF). The IRF is likely to depend on the orientation of crystal growth since each growth orientation will have different energetic barriers to crystallization and different maximum velocities. It was generally assumed, with little experimental justification at the atomic scale, that an amorphous phase will have the same general kinetic characteristics; namely that a moving amorphous/liquid interface is identifiably present and thus the kinetics of this motion can be expressed according to an interface response function in a manner analogous to that representing the motion of a crystal/liquid interface. However, the legitimacy of this assumption has been recently challenged by computational results which suggest that amorphous silicon melts by a two-mode melting mechanism.21 A more in-depth discussion of the relation between the IRF and explosive crystallization is given in refs 21 and 22. Until recently, little was known about the explosive crystallization process on an atomic level since it occurs much too fast for microscopic analysis using experimental approaches. However, we have recently characterized explosive crystallization process at an atomic scale through a specially constructed molecular dynamics technique.21-22 This work showed that the Stillinger-Weber potential formalism presents the best semiempirical description of the kinetics of silicon and germanium condensed phase transitions but that the original parametrization of Stillinger and Weber had difficulty describing the amorphous phase of silicon.21 A new parametrization for both silicon and germanium21 was more successful, and this parametrization will be used in the first part of this work. To date, there has been only one experimental study on explosive crystallization of Si1-xGex alloys reported in the literature.23 Yu used a lateral explosive crystallization technique to process alloys where x ranged from 1 to 0.95 and a vertical explosive crystallization technique to process alloys where x ranged from 0 to 0.35. The results of this study demonstrated a considerably nonlinear dependence of velocity on composition where the velocity decreased dramatically as the percentage of the minority component increased for both Ge-rich and Si-rich alloys. In the composition range 0.90 e x < 0.95, explosive crystallization was not self-sustaining. No experiments were reported between 0.35 < x < 0.90. Yu and co-workers24 also performed a complementary simulation study of the crystallization of Si1-xGex alloys from the melt. This study involved simulating the effect of a laser-induced melt and subsequent fast crystallization sample of Si1-xGex alloys, including a determination of the crystallization velocity. Using the original parametrization of the Stillinger-Weber potential, they provided maximum velocity results in agreement with their experimental results, showing a nonlinear dependence of the velocity with Ge composition. The goal of this research is to study explosive crystallization processes for SiGe alloys and boron-doped SiGe alloys. We will study the effect of alloying and the presence of boron dopants on the explosive crystallization process, emphasizing the ability of silicon-germanium alloys to crystallize across the entire compositional range, a study which has not been made experimentally so far. We will present results concerning the ability of the process to crystallize amorphous material without changing the dopant distribution in the system, a technologically relevant design criterion. Simulation Details Explosive crystallization will be simulated using a molecular dynamics method we developed for Si.22 The initial configuration of the system is shown in Figure 1a. It consists of 10 240

Figure 1. (a) Schematic of the simulation cell showing three coexisting phases (crystal-liquid-amorphous), the placement of the heat baths, and static regions at either end. (b) Starting configuration for explosive crystallization samples initiated with a ‘laser’ in which no liquid layer is initially present. The arrow denotes the region affected by the ‘laser’. (c) Initial setup for samples used to determine interface response functions for undoped Si1-xGex alloys.

atoms arranged in a rectangular box with approximate dimensions of 440 Å × 22 Å × 22 Å (an aspect ratio of 20:1). This box was shown previously to be large enough to capture the dynamics of the process while maximizing computational speed;22 in addition, simulation cells containing 6656 atoms were tested and proved too small to establish a steady state. The initial sample consists of a crystalline section (66 Å) adjoining a liquid region (22 Å) and an amorphous section (352 Å). The amorphous material possessed a coordination number of 4.11 and an average bond angle of 108° ( 15°, consistent with ref 25. Samples of Si1-xGex alloys were constructed to represent the entire compositional range. Initially, all atoms are designated as being silicon. A random number generator then selected atoms to be relabeled as germanium up to the desired composition, and the density of the sample was adjusted to alleviate the strain induced from this compositional change. Experimentally, the explosive crystallization process begins when an amorphous sample is subjected to intense energy, such as by a laser. Thus, starting with the configuration shown in Figure 1a is equivalent to observing the process after initiation has taken place. We will investigate below whether this has any bearing on the subsequent crystallization process. Periodic boundary conditions were used in the y- and z-directions only. In the x-direction, the first four atomic layers at each end of the cell were frozen (to minimize edge effects) while the next eight layers were used as heat baths to control heat loss by rescaling the velocities of the atoms in this region every time step. The remainder of the cell was unconstrained using Brown and Clarke’s integration scheme26 for the equations of motion.

5630

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

For SiGe alloys, the interatomic potential models are based on the Stillinger-Weber (SW) potential model27 which consists of a term for pairwise interactions and one that represents threebody interactions. Two different parameter sets are discussed in this paper: One set is based on the ‘original’ SW potentials, i.e., the original Si-Si potential27 and the Yu et al. parametrization for Ge-Ge and Si-Ge interactions.24 The second parameter set (referred to as ‘modified’ SW) is taken from ref 21. For the modified SW potentials, the Si-Ge cross terms were determined using traditional ad hoc mixing rules that employ an arithmetic mean for the length parameter (σ) and a geometric mean for the energy parameters ( and λ). This approach for Si-Ge models was also used by Yu et al.24 Using two different parametrizations also provides a means to evaluate sensitivity to parameter set. For SiGe alloys containing boron, it was only practical to employ the original SW parameter set for Si-Si and Ge-Ge interactions since a complete set of models exists for a system containing boron within a SiGe matrix: The B-B interaction parameters used here are those of Rasband and Clancy,28 the Si-B interaction parameters are those of Noya et al.,29 and the Ge-B interaction parameters are due to Wang and Clancy.30 No such parameter set is available for the modified SW parameters. For boron-doped alloys, two sets of simulations were performed. The first set was performed as described above while the second set simulated the action of a virtual laser, modeled as a sudden energy input to the system, to initiate melting in the box; see Figure 1b. The second set of simulations did not initially contain a mediating liquid layer; which was instead generated by the simulated laser. Using the virtual laser to initiate melting should provide a more realistic representation of the experimental technique where the liquid layer is formed by the action of a nanosecond pulsed laser. The action of a laser was modeled by increasing the temperature of the first 66 Å of amorphous material adjacent to the crystalline material to 1550 K for a time of 30 ps after which time the region was modeled as unconstrained. These conditions were chosen to give sufficient time to fully form a liquid layer (10 ps proved sufficient) of comparable size to the steady-state liquid layer (27-117 Å for pure samples as shown in ref 22) at a temperature that will not melt the crystalline material. Using a higher intensity ‘laser’ (say, with identical time and melt width values but a temperature of 2500 K) initiated unwanted melting in the crystalline material adjacent to the amorphous layers and was thus rejected. In each case, whether laser-initiated or with a preexisting liquid layer, the samples were evolved until one of three conditions occurred: a steady state was reached, crystal growth was quenched by the amorphous side heat bath, or the sample completely crystallized. The definition of ‘steady state’ was taken to be the point at which the width of the liquid layer remained constant for at least 250 ps. For a study that focuses on the kinetic properties of moving interfaces, it is essential to be able to determine the location of the interface and follow its progress through the sample as a function of time. The primary method used to determine the location of the interface was the angular order parameter of Uttormark et al.31 which designates an atom as being either ‘solid-like’ or ‘liquid-like’. To be solid-like, an atom has four nearest neighbors with the nearest neighbor cutoff being chosen as the first minimum in the radial distribution function and a value for A < 0.4 determined by the following equation

A)

∑i (cos θi + (1/3))2

where the sum is taken over the triplet sets of the four nearest neighbors and θi is the angle made by the triplet set (for the tetrahedral angle, cos θ ) - (1/3) so A will equal zero for perfect tetrahedral structure). To track interface motion, the sample is divided up into slices 1-2 atomic layers thick in the direction of crystallization (the z-direction). Each atom can then be assigned as either ‘solid-like’ or ‘liquid-like’, allowing a calculation of the fraction of solid-like atoms in each slice and hence the determination as to the prevalent phase in a given slice: characteristic values were obtained for bulk samples of crystal (where the fraction of solid-like atoms is >98%), amorphous (60-85% solid-like), and liquid phases (1 K/Å) may affect the diffusion in the liquid layer thereby affecting crystal growth velocity. The density difference mentioned above also helps to explain the result that velocities predicted using the original SW models show higher values if derived from the interface response functions, whereas the velocities predicted using the modified SW models show higher values if derived from the explosively

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5633

Figure 7. Interface velocity as a function of Si composition using modified SW potentials. Open squares: IRF-determined maximum velocities; filled squares: maximum velocities from explosive crystallization runs. Crosses: experimental data;23 plus signs represent the maximum velocity as determined by the interface response function simulations of ref 24. The dashed line represents a linear velocity-composition relationship, showing the nonlinearity of the data.

crystallized samples. We suspect that this difference can be explained in terms of the differences in liquid densities which also change sign when switching from the original SW potentials to the modified SW potentials. In summary, crystallization velocities derived from the explosive crystallization runs do not match those from the IRF-derived values due to the volume constraints imposed on the liquid layer by the presence of the two solid phases. We believe that the velocity values predicted by the interface response functions are the more accurate values since the presence of the vacuum at the edge of the IRF samples allows the liquid region more freedom to equilibrate at its correct density while the explosively crystallized samples are fixed at the density set by the two adjoining solid phases. This contention is supported by the fact that equilibrium densities for the IRFderived samples are closer to equilibrium values for both the modified and original SW potentials. Comparison is also made to simulation data by Yu et al.24 who derived maximum velocity data from interface response functions. As can be seen in Figure 7, Yu’s data are significantly different to the results obtained here for xGe g 0.25. We believe that the explanation for this discrepancy arises from the manner in which Yu’s simulations were carried out: In Yu’s case, the initially crystalline sample was melted with a simulated laser. Once the melting is complete, and the ‘laser’ turned off, crystal regrew at the expense of the melt, but the temperature and velocity of the crystal-liquid interface are not immediately at steady state. Yu and co-workers obtained an interface response function for a given alloy composition using the transient temperature and velocity prior to reaching steady state, essentially capturing the velocity-temperature data ‘on the fly’ as the interface temperature changed with time during the crystallization process. In contrast, our simulation data only took data points once the temperature profile and velocity had remained constant for at least 300 ps. In addition, Yu was limited computationally to much smaller system sizes: while Yu’s systems were large enough to crystallize a distance of ∼30 Å of alloy, the results provided in this work crystallized a distance of over 100 Å. We consider the method used here to obtain IRF-derived velocity-time data to be superior to Yu’s. The

closer proximity of Yu’s results to the experimental data is, we believe, largely fortuitous. Comparison of the crystallization velocities predicted by explosive crystallization simulations to the experimental results of Yu23 is shown in Figures 6 and 7. The experimental data show that alloying has a dramatic effect on the maximum crystallization velocity. The velocity plummets when the amount of Ge in a Si-rich sample is increased from 0 to 35%, or when the amount of Si in a Ge-rich sample increases from 0 to 5%. No discernible explosive crystallization could be observed for SiGe samples containing between 90 and 95% germanium.23 Our computational results show a similar decrease in crystallization velocity as Ge is added to a Si-rich sample (from 65 to 100% Si), but show a relatively constant velocity over the composition range 0-50% Si, within the uncertainty in the velocity data. Discussion of EC and IRF Results. The most striking difference between the computational and experimental results is that the simulations predict that explosive crystallization can occur over the entire compositional range, whereas Yu’s experiments were unable to observe explosive crystallization above 5% Si. This does not necessarily mean that explosive crystallization cannot occur at these conditions, but it is very difficult to observe experimentally. For example, Yu suggests the possibility of success with higher substrate temperatures. However, at higher substrate temperatures, the competing process of solid phase epitaxy can also take place with the result that the explosive crystallization can only travel a short distance before it reaches crystalline material formed by solid phase epitaxy, thereby halting the process. In summary, only the results for simulations of pure germanium and Si1-xGex samples containing less than 25% germanium are qualitatively comparable to the experimental data. The qualitative discrepancy between experimental and computational results for samples with more than 25% Ge, which are particularly noticeable for Ge-rich samples, may be caused by the parameter set used for Ge, the mixing rule used for SiGe interactions, or, more broadly, point to inadequacies of the SW potential model itself that cannot be fixed by any reparametrization. We will now consider these possible sources of error in turn. The Ge potential used here is based largely on the rather simplistic model used by Yu et al.24 in which the only changes made to the original SW model for Si were to rescale parameters to reproduce experimental values of the cohesive energy and lattice constant while leaving the shape of the potential unaltered. The reparametrization of this model re-weighted the respective contributions of the two- and three-body terms, but left the shape of the individual terms unaltered. Wang and Clancy showed that the SW pair potential using Yu’s germanium parametrization is significantly more short-ranged than that of a more accurate density functional theory calculation.30 The interaction length may play a bigger role in simulations of interfaces such as these, where the interfacial atoms will ‘feel’ the influence of atoms in the adjoining phase. The same simple mixing rule used by Yu and co-workers was used here to obtain the Si-Ge cross terms for the modified potentials. Unfortunately, there is no way to know the effect of the interaction length or how to set the mixing rule without comparison to ab initio calculations. Finally, while we believe the modified SW potentials to be an improvement over the original SW potentials, they still remain a relatively simple empirical model. Again, the consideration of alternative, more accurate, potential models is currently not possible. Even though tight-binding parameters exist for Si and Ge32-35 and hence a tight-binding molecular dynamics simula-

5634

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Figure 8. An example of the slowing of the propagating interface due to the presence of Ge (top). This sample is a snapshot of an explosive crystallization sample of Si0.95Ge0.05. The thick black line represents the interface. The presence of the germanium atoms to the left of the positions labeled A and B retard the interface, leading to its nonplanarity. Also presented is an example of the effect on the shape of the interface of adding Ge to SiGe, here for an explosively crystallized sample of Si0.78Ge0.22 (bottom). The increased percentage of Ge causes the interface, represented by the thick black line, to again become planar.

tion of the explosive crystallization process is theoretically possible, it is practically impossible due to the computational effort necessary to simulate systems containing enough material (around 10 000 atoms) to produce artifact-free data. Results using the modified SW models and, to a lesser extent, those using the original SW model, mimic the nonlinear decrease in crystallization velocity observed experimentally as the percent of germanium range from 0 to 25%. In fact, the velocity change over the range of 0-20% Ge for the modified SW potentials of 6 m/s compares favorably with the experimental change of 7 m/s. The nonlinear nature of this composition dependence arises, we believe, because the advancing interface slows while it incorporates the significantly larger Ge atoms onto lattice sites. This phenomenon results in a locally nonplanar growth front in the vicinity of Ge atoms illustrated in Figure 8. With increasing Ge concentration, this effect is exacerbated. Once the composition of Ge reaches 20-25%, the nonplanarity of the interface is greatly diminished since there are germanium atoms along the entire interface, also shown in Figure 8. A similar interface slowing was observed by Yu and Clancy in simulations of SiGe crystallization from a melt.36 This may also help to explain why the velocity-composition profile is flat for germanium-rich alloys. With germanium atoms positioned along the length of the interface, the entire interface is slowed and further addition of germanium does not cause a significant change in velocity. Width of the Mediating Liquid Layer. It was also possible to predict the extent of the width of the mediating liquid as a function of alloy composition. This information has not been determined experimentally. The width of the liquid layer stayed relatively constant throughout the composition range of 0-25% germanium most accurately predicted by the simulations and, indeed, over the entire composition range. The liquid is 65 ( 10 Å wide for the original SW potentials, and 51 ( 9 Å for the modified SW potentials. As discussed in ref 22, the width of

the liquid layer depends primarily on the amorphous-side heat bath. We expect the width of the liquid layer to be relatively uncoupled from the crystallization velocity, assuming that the liquid layer width has a non-negligible thickness. Considering the stability of the process, it is not surprising that liquid layer widths for SiGe alloys are comparable to those obtained for pure components. Segregation During Crystallization. We can also calculate the segregation coefficient, the preference of one component to reside in the solid or liquid phase. This is defined as the mole fraction of the component in the solid phase divided by that of the liquid phase. The segregation coefficients for germanium in Si1-xGex alloys were determined by dividing the mole fraction of Ge in a 5.5 Å slice (roughly two atomic layers) immediately before the interfacial slice by the mole fraction of Ge in a 5.5 Å slice immediately after the interface. Use of smaller slices gave poor statistics for the 5% germanium case (where one atomic layer contains only three Ge atoms). The roughness of the interface for Si-rich samples could make evaluating the segregation coefficients difficult; however, since the interfacial slice was not included in the calculation this effect should be minimized. At these high growth velocities (>7 m/s), the segregation coefficient is essentially unity, implying that no segregation is observed for any composition in agreement with similar calculations by Yu and Clancy.57 However, the interface roughening observed during the crystallization of the Si-rich samples (see Figure 8) suggests that the interface prefers to bind silicon atoms, i.e., implying that segregation is occurring. There is, in fact, some experimental evidence to support this idea: Experimental values of the segregation coefficients for SiGe alloys containing less than 4% Ge have been determined to be ∼0.85 for a velocity near 4 m/s,37 suggesting that segregation may indeed occur at the velocities used here. Explosive Crystallization of Boron-Doped Si1-xGex Alloys The computational ‘doping’ process used here was achieved using a random number generator to select atoms from among the 10 240-atom sample to be relabeled as boron atoms until the desired percentage of boron was reached. We considered the effect of boron doping on the explosive crystallization of two alloys, Si0.91Ge0.09 and Si0.78Ge0.22, chosen to be compositions shown by Yu’s experiments23 to be capable of exhibiting explosive crystallization. The pure components and the two selected SiGe alloys were doped with increasing amounts of boron, up to 8 at. % to determine the doping limit for which crystallization still occurred. As we shall show, not all the systems we investigated were capable of supporting explosive crystallization. Results for successful explosive crystallization runs are shown in Figures 9 and 10, where it can be seen that the limit of boron doping varies considerably with Ge composition. For pure Si, the explosive crystallization process was predicted to occur for concentrations up to 4% boron. In contrast, explosive crystallization could only be observed for doping levels up to 1% boron for pure Ge and Si0.91Ge0.09 samples. The doping limit for the Si0.78Ge0.22 alloy was the lowest value; only 0.5% boron was capable of supporting explosive crystallization. Higher concentrations of boron than the limits given in the preceding paragraph led to one of two scenarios. In one scenario, explosive crystallization began, only to be followed by recession of the crystal front with the eventual conversion of the entire sample to liquid; the amorphous front also continued to recede throughout the process. In the second scenario, crystal growth

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5635

Figure 9. Interface velocity as a function of boron composition for systems modeled using original SW potentials. Open data points represent laserinitiated simulations; filled data points represent simulations which initially possessed a liquid layer. Vertical lines connect data at the same composition. Dashed lines represent a linear fit to data for a given Ge composition. Symbols designate the following compositions: diamonds, Si; squares, Ge; triangles, Si0.91Ge0.09; circles, Si0.78Ge0.22.

Figure 11. Effect of boron atoms to disrupt the advancement of a planar growth front. The three pictures are snapshots taken 25 ps apart for a Si system containing 1% boron, modeled using original SW potentials. From top to bottom: (a) the advancing crystal front slows near a boron atom. (b) The interstitial boron atom continues to slow crystallization. (c) The crystallization front regains its roughly planar symmetry with the boron atom incorporated into the lattice as an interstitial.

Figure 10. Liquid layer width as a function of boron composition for systems modeled using the original potentials. Key as in Figure 9.

would not commence at all. For all the cases we investigated, if the original set of temperature conditions, namely the heatbath temperatures and initial kinetic energy, did not result in a successful explosive crystallization run, new temperature conditions were set until the sample was successfully explosively crystallized or until the sample was deemed to be unable to explosively crystallize under the current computational setup. A successful explosive crystallization run was defined as exhibiting steady-state crystallization for at least 250 ps. As shown in the results depicted in Figures 9 and 10, increasing the amount of boron in the sample makes crystallization increasingly difficult. In cases where the crystal grows and then recedes, the initial conditions were clearly amenable to growth at the start of the simulation. However, after the crystal/liquid front has propagated a short distance, at most 20 atomic layers, there is a gradual build-up in the heat generated in the vicinity of the growing interface and, at some point, the crystal/liquid front is forced to recede. These results imply that increasing the amount of boron in the sample destabilizes the crystallization process. In cases where growth could not be sustained, the crystal front does not have sufficient energy to advance into the liquid and the liquid layer takes on glassy characteristics.

Effect of B-Doping on Velocity and Liquid Layer Width. As we have seen above, doping SiGe samples with boron has a significant effect on the overall velocity of the crystallization process. It also has a significant effect on the width of the liquid layer, increasing as the extent of boron is increased. The results for all the samples we investigated show a scattered but approximately linear dependence of the velocity on boron concentration, as shown in Figure 9. There is a drop of approximately 2 m/s in crystallization velocity for every 0.5% increase in boron concentration for all of the samples studied except for the Si0.78Ge0.22 samples. For this last case, the dependence of velocity on concentration is significantly greater, 4.9 m/s drop in velocity per 0.5% increase in boron concentration; 1% boron-doping for this alloy could not support explosive crystallization. Since we have only one data point at this composition, its results should carry less weight. Atomic representations of the crystallization process were analyzed to determine the reason for the drop in velocity with increasing boron concentration. As shown in Figure 11, the advancing crystallization front is slowed in the vicinity of the boron impurity atom causing a nonplanar crystallization front. As the front passes, the boron atom is incorporated into the lattice at an interstitial position due to the speed of the passing crystallization front. Subsequently, the boron atom annihilates with an adjacent vacancy and moves to a lattice site. Since the boron impurity atoms do not join the crystallizing front as rapidly as the silicon and germanium atoms, there is an overall slowing of the process with the presence of boron. The trapping

5636

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Table 2. Percent of Boron Atoms Occupying Interstitial Sites for SiGe Systems for Boron Doping up to 1% Following Explosive Crystallization

sample

percentage of B atoms on interstitial sites for a 0.5% B sample

percentage of B atoms on interstitial sites for a 1% B sample

Si Si0.9Ge0.1 Si0.75Ge0. 25 Ge

0 21 24 20

0 32 no successful EC 0

Original SW potentials.

of boron on interstitial sites does not appear to be related to boron clustering since there is generally only one boron atom in the vicinity. The boron does not appear to diffuse through the solid and simply find a vacancy; rather, the boron eventually occupies the lattice site (annihilating the vacancy) that it should have settled upon as the front passed. It should be noted that other mechanisms are possible for boron incorporation into the lattice. The effect of the extent of boron doping on the width of the liquid layer was also studied with the results presented in Figure 10. For the pure Si samples, the concentration of boron in the system has very little effect on the liquid layer width until the concentration of boron is high enough (2% boron) to cause a few boron atoms to occupy interstitial sites. In contrast, samples containing germanium showed a strong, roughly linear, increase in the width of the liquid layer as the amount of boron present is increased. The trend indicates that more germanium in the lattice gives a steeper increase in the width of the liquid layer as the amount of boron is increased but, given the scatter in the few data points, it is hard to state this conclusively. The borondoped pure Ge sample shows the largest increase in liquid layer width from 79 Å for pure Ge to 212 Å for a 1% boron-doped sample. Changes in the liquid layer width for the Si0.78Ge0.22 and Si0.91Ge0.09 alloy samples were smaller: This ranged from 79 to 126 Å, for 0-0.5% boron doping for Si0.78Ge0.22, and 64-106 Å for 0-1% boron doping for Si0.91Ge0.09. These results can be understood when considered in conjunction with the observation of the failed runs in which the crystal front receded leaving an entirely liquid sample. As the level of boron doping increases, it became increasingly more difficult to find a set of conditions which would successfully explosively crystallize the sample. Increasing boron doping destabilizes the process resulting in a larger liquid layer width. Final Doping Profile and Activation of Boron Atoms. For industrial applications, any significant accumulation of boron atoms must be considered since control of the final doping profile is crucial. No significant boron redistribution or accumulation was observed for all of the samples studied here except the 1% boron-doped Si0.91Ge0.09 sample, for which two of the four runs showed boron redistribution and clustering which will be discussed further below. Additionally, it is important that boron atoms occupy electronically active lattice sites. Table 2 clearly demonstrates an increased tendency of the boron atom to occupy an interstitial site in SiGe samples compared to pure silicon samples. This tendency has been investigated by Wang and co-workers30 in which extensive ab initio calculations were performed to determine the energetics of interstitial defects involving boron in Si, Ge, and SiGe. Their results show that germanium atoms have an increased ability (relative to silicon) to displace a substitutional boron atom to an interstitial site. However, the migration energy of a boron interstitial is increased in the presence of Ge, and the boron atom is trapped in the vicinity of the Ge atom on an interstitial

Table 3. Samples Containing Multiple Growth Orientationsa

sample Si0.6Ge0.4 Si Si Si

number percentage of distinct boron orientations growth orientations 0 0.5 1 2

3 2 3 2

[100], mixed, [100] [100], [311] [100], [311], mixed [100], mixed

comment laser-initiated laser-initiated

a Original SW potential models. ‘Mixed’ designates the orientation depicted in Figure 23. Growth orientations listed in the order in which they were present along the length (z-direction) of the sample. The comment ‘laser-initiated’ signifies that the liquid layer in the sample was initially formed through the action of a simulated ‘laser’.

Figure 12. Mixed growth orientation sometimes observed after reorientation. The picture has been rotated to show a cross-section of the crystal in the direction of its growth (Si shown in black, Ge shown in gray).

site. Experimentally, it has been shown that introducing Ge into the Si lattice increases the activation of the sample with respect to boron,38,39 i.e., it increases the number of boron atoms on lattice sites. This result is in contrast to our result, as well as that of Wang and co-workers, which may be due to the short time scales accessible to molecular dynamics. If the interstitial boron atoms take longer than a few nanoseconds to relax to a lattice site, it cannot be observed here. Table 2 also shows an anomalous decrease of the percentage of B atoms occupying interstitial sites from 20 to 0% in B-doped Ge as the doping percentage increases from 0.5% to 1%. We believe this is due to the slower crystallization velocity of the 1% doped sample (3.6 m/s vs 6 m/s for the 0.5% B-doped sample) allowing sufficient time for the boron atoms to locate on lattice sites as the crystallization front passes. Here, an atom is considered to occupy an interstitial site if it is located more than 1 Å from a lattice site. Implications of Halting and Restarting the Crystal Front For several cases listed in Table 3, the crystal front was observed to halt for several picoseconds and then restart. Investigation of the samples showed that, after restarting, the orientation of subsequent crystal growth had changed. In some cases, this change in crystal orientation occurred twice, leaving three different orientations of crystal. The propensity of the system to alter its crystal orientation was observed for most of the systems we considered, including boron-doped Si, certain SiGe alloys, as well as the pure components, Si and Ge, as described in ref 22. The new growth direction was either [311] or the growing crystal was a mixed structure which had both [100] and hexagonal motifs as shown in Figure 12. This reorientation phenomenon may be viewed as the spontaneous formation of polycrystalline material or may arise as an artifact of the simulation. We have observed that the reoriented growth aligns itself to reduce lattice mismatch resulting from the use

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5637 Table 4. Size and Composition of Clusters Formed during Explosive Crystallization of a 1% B-Doped Si0.91Ge0.09 Samplea cluster size (# atoms)

# Si atoms

# Ge atoms

#B atoms

cross-section (Å2)

comment

61 60

28 28

23 21

10 11

154 117

-Laser- initiated

a

‘Laser-initiated’ defined as in Table 3.

Figure 13. Boron profile of a 1% boron-doped Si0.91Ge0.09 sample along the direction of crystallization. Dashed line: initial boron distribution; solid line: distribution after 6 ns. Note boron accumulation at a location 240 Å along the z-axis after the front passes.

of periodic boundary conditions, that is, the lattice at the bottom of the box lines up with the lattice at the top of the box. Periodic boundary conditions undoubtedly affect the crystallization process in these cases; however, studying this phenomenon in detail will involve very large-scale simulation samples beyond our scope. Boron Redistribution during Explosive Crystallization As noted above, in two of the four simulation runs for 1% boron-doped Si0.91Ge0.09 samples, accumulation of Si, Ge, and B into clusters was found to occur. In these instances, the crystal/ liquid front stopped advancing due to an increase of boron atoms in the liquid layer near, defined to be within eight atomic layers of, the crystal/liquid interface. The increase in boron atoms was not initially due to diffusion but to inhomogeneities in the initial ‘doping’ of the sample. The clustering of boron atoms only occurred when the advancing crystal front propagated through a section characterized by a low concentration of boron atoms followed by a region which was boron-rich, relative to the previous slice or two. Figure 13 shows the evolution of the doping profile with time for one of the cases exhibiting the accumulation of boron atoms. There was (by chance) a rather nonuniform initial distribution of boron atoms; there were only two boron atoms in a 44 Å region 138-182 Å from the crystal edge of the box, followed by a sharp increase to 16 boron atoms over the next equivalent distance of 44 Å (182-226 Å from the crystal edge of the sample). When the crystallization front reached the increased concentration of boron atoms at ∼182 Å, the front halted to wait for some of the boron atoms to diffuse away from the crystallization front. While the crystallization front is halted, the extended period of time for which the stagnant liquid layer exists allowed both Ge and B atoms in the liquid layer to form a cluster with the Si atoms at 231-253 Å from the crystal edge of the sample which is at least quasistable and which does not dissolve over the time scale of the simulations performed here (∼7 ns). When the crystallization front reaches this atomic cluster, the crystallization is effectively halted by the presence of the cluster. It is important to note that, in these two cases, the crystal section does not start to recede as is seen in some of the failed conditions mentioned above. Formation of the cluster prevents further growth but does

Figure 14. Atomic representation of a 60-atom Si0.47Ge0.35B0.18 cluster formed during the explosive crystallization processing of a 1% boron-doped Si0.91Ge0.09 sample modeled using original SW potentials. Blue: Si. Yellow: Ge. Red: B.

not cause the crystallization process to become unstable. Clustering effectively reduces the average concentration in the bulk allowing growth in the remaining regions. The two clusters observed in these cases were very similar, both in size and atomic composition. The clusters involved 61 and 60 atoms, respectively, and both clusters have the same structure and composition, as shown in Table 4. Although hard to perceive from the 2D representation of the 3D structure shown in Figure 14, the structure is a diamond lattice made up of alternating boron and germanium atoms on substitutional lattice sites with all of the silicon atoms occupying tetrahedral interstitial positions. The clusters are roughly football shaped, resulting in a different cross-sectional area parallel to the propagating front. The area of the cluster parallel to the plane of the moving interface covers between one-fourth to one-third of the interface, 117 and 154 Å2, compared to the cross-section of the box of 484 Å2. It might be possible that, with a larger sample cross-section, crystallization would be able to occur around the cluster and simply incorporate it into the crystallized material; this was not explored. Last, the long lifetime of the clusters observed here was also seen in simulations of borondoped SW-modeled Si during laser thermal annealing by Wang et al.40 They showed that B-B clusters modeled using a SW potential exhibited a much longer lifetime than those interacting with a tight-binding model (>150 ps for the SW model vs 5 ps for tight-binding). Thus the clusters observed here may have artificially long lifetimes in comparison to more realistically modeled systems or to experimental situations. Conclusions Simulations with two different SW parameter sets indicate that explosive crystallization should occur over the entire compositional range of Si1-xGex alloys at velocities high enough to prevent Ge segregation. In contrast, limited experimental results show quenching of explosive crystallization at xSi g 10% and xGe g 35%. Simulation results using modified SW models

5638

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

showed qualitative agreement with the experimental results for Si-rich systems: Over a range of 0-25% germanium, the simulations predicted a nonlinear decrease in velocity that roughly matched the slope and extent of the rapid drop in interface velocity. The nonlinear dependence of velocity on composition results from slowing of the interface due to the presence of germanium atoms. In contrast, both parameter sets give qualitatively incorrect results for Ge-rich samples. Increasing the amount of silicon in a germanium-rich sample has been shown by Yu’s experiments to sharply decrease the interface velocity, whereas the computational results presented here predict a very gradual increase in velocity from 0 to 50% Si. Overall, these results suggest that the parametrization and/or the potential model may not be appropriate to model the binary mixture of silicon and germanium for alloys containing more than 25% germanium. Alternatively, other mechanisms may limit explosive crystallization that cannot be captured in the finite size of molecular dynamics simulations. It would be helpful for new experiments to be carried out, taking advantage of advances in technique since Yu’s experiments a decade ago. Boron-doped Si1-xGex alloys showed that the limit of incorporation of boron during explosive crystallization is strongly dependent on the amount of germanium present: the more Ge present, the less boron can be incorporated. Successful explosive crystallization runs were performed for boron doping as high as 4% for pure silicon, 1% for pure germanium and a Si0.91Ge0.9 alloy, and 0.5% for a Si0.78Ge0.22 alloy. While undoped samples possessed consistent liquid layer widths and velocities across the composition range, boron-doped samples showed a widening liquid layer and lowering of the velocity as the amount of boron doping was increased. This was especially evident in samples containing germanium. The decrease in velocity was due to slowing of the growth front to incorporate boron atoms into the lattice; the widening of the liquid layer resulted from destabilization due to the presence of the boron atoms. Acknowledgment The authors thank Sarah Meeker (NSF-REU program) for preliminary calculations and the National Science Foundation for financial support of this work through award no. 9980100. The NSF-sponsored Cornell Center for Materials Research provided some of the necessary computational resources. Literature Cited (1) The International Technology Roadmap for Semiconductors; International Sematech: Austin, TX, 2004. (2) Agarwal, A.; Fiory, A. T.; Gossman, H.-L. Ultra-shallow junctions by ion implantation and rapid thermal annealing: spike-anneals, ramp rate effects. Mater. Res. Soc. Symp. Proc. 1999, 568, 19. (3) Fiory, A. T. Recent developments in rapid thermal processing. J. Electron. Mater. 2002, 31, 981. (4) Cullis, A. G.; Webber, H. C.; Chew, N. G. Correlation of the structure and electrical-properties of ion-implanted and laser-annealed silicon. Appl. Phys. Lett. 1980, 36, 547. (5) Wood, R. F.; Kirkpatrick, J. R.; Giles, G. E. Macroscopic theory of pulsed-laser annealing. II. dopant diffusion and segregation. Phys. ReV. B 1981, 23, 5555. (6) Narayan, J.; Holland, O. W.; White, C. W.; Young, R. T. Excimer laser annealing of ion-implanted silicon. J. Appl. Phys. 1984, 55, 1125. (7) Narayan, J.; Holland, O. W.; Christie, W. H.; Wortman, J. J. Rapid thermal and pulsed laser annealing of boron fluoride-implanted silicon. J. Appl. Phys. 1985, 57, 2709. (8) Rajendran, K.; Schoenmaker, W. Measurement and simulation of boron diffusion in strained Si1-xGex epitaxial layers with a linearly graded germanium profile. Solid-State Electron. 2001, 45, 1879. (9) Zangenberg, N. R.; Fage-Pedersen, J.; Hansen, J. L.; NylandstedLarsen, A. Boron diffusion in strained and relaxed Si1-xGex. Defect Diffus. Forum 2001, 194, 703.

(10) Thompson, M. O.; Galvin, G. J.; Mayer, J. W.; Peercy, P. S.; Poate, J. M.; Jacobson, D. C.; Cullis, A. G.; Chew, N. G. Melting temperature and explosive crystallization of amorphous-silicon during pulsed laser irradiation. Phys. ReV. Lett. 1984, 52, 2360. (11) Donovan, E. P.; Spaepen, F.; Turnbull, D.; Poate, J. M.; Jacobson, D. C. Heat of crystallization and melting-point of amorphous-silicon. Appl. Phys. Lett. 1983, 42, 698. (12) Spaepen F.; Turnbull, D. In Laser-Solid Interactions and Laser Processing-1978; Ferris, S. D., Leamy, H. J., Poate, J. M., Eds.; American Institute of Physics Conference Proceedings No. 50, New York, 1979; p 73; Bagley, B. G.; Chen, H. S. ibid. p 97. (13) Lowndes, D. H.; Jellison, G. E., Jr.; Pennycook, S. J.; Withrow, S. P.; Mashburn, D. N. Direct measurements of the velocity and thickness of explosively propagating buried molten layers in amorphous-silicon. Appl. Phys. Lett. 1986, 48, 1389. (14) Stolk, P. A.; Polman, A.; W. C. Sinke, W. C. Experimental test of kinetic theories for heterogeneous freezing in silicon. Phys. ReV. B. 1993, 47, 5. (15) Murakami, K.; Eryu, O.; Takita, K.; Masuda, K. Explosive crystallization starting from an amorphous-silicon surface region during long-pulse laser irradiation. Phys. ReV. Lett. 1987, 59, 2203. (16) Thompson, M. O.; Mayer, J. W.; Cullis, A. G.; Webber, H. C.; Chew, N. G.; Poate, J. M.; Jacobson, D. C. Silicon melt, regrowth, and amorphization velocities during pulsed laser irradiation. Phys. ReV. Lett. 1983, 50, 896. (17) Chojnacka, A. P.; Thompson, M. O. In Growth, EVolution and Properties of Surfaces, Thin Films and Self-Organized Structures; Moss, S. C., Poker, D. B., Ila, D., Eds.; Mater. Res. Soc. Symp. Proc. 648, Warrendale, PA, 2001; p P11.12.1-8. (18) Bostanjoglo, O.; Endruschat, E. Kinetics of laser-induced crystallization of amorphous-germanium films. Phys. Status Solidi A 1985, 91, 17. (19) Bensahel, D.; Auvert, G. In Laser-Solid Interact. Transient Therm. Process. Mater.; Marayan, J., Brown, W. L., Lemons, R. A., Eds.; Mater. Res. Soc. Symp. Proc. 13, New York, 1983; p 165-176. (20) Chojnacka, A. P. Explosive crystallization of germanium films: kinetics and morphologies. Ph.D. Thesis, Cornell University, 2002. (21) Albenze E. J.; Clancy, P. Interface response functions for amorphous and crystalline Si and the implications for explosive crystallization. Mol. Simul. 2005, 31, 11. Albenze, E. J.; Matejik, L. A.; Fynan, N. F.; Clancy, P. In Amorphous and Nanocrystalline Silicon-Based Films; Abelson, J. R., Ganguly, G., Matsumura, H., Robertson, J., Schiff, E. A., Eds.; MRS Symposia Proceedings No. 762 Materials Research Society, Pittsburgh, PA, 2003; p A16.4.1-6. (22) Albenze, E. J.; Thompson, M. O.; Clancy, P. Atomistic computer simulation of explosive crystallization in pure silicon and germanium. Phys. ReV. B. 2004, 70, 094110. (23) Yu, Q. Computer simulation and experimental investigation of crystal growth in group IV semiconductor alloys. Ph.D. Thesis, Cornell University, 1995. (24) Yu, Q.; Thompson, M. O.; Clancy, P. Solidification kinetics in SiGe alloys. Phys. ReV. B. 1996, 53, 8386. (25) Choudhary, D.; Clancy, P. Characterizing the nature of virtual amorphous silicon. J. Chem. Phys. 2005, 122, 174509. (26) Brown, D.; Clarke, J. H. R. A comparison of constant energy, constant temperature and constant pressure ensembles in molecular-dynamics simulations of atomic liquids. Mol. Phys. 1984, 51, 1243. (27) Stillinger, F. H.; Weber, T. A. Computer simulation of local order in condensed phases of silicon. Phys. ReV. B. 1985, 31, 5262. (28) Rasband, P. B.; Roberts, B. W.; Clancy, P. Tight-binding studies of the tendency for boron to cluster in c-Si. I. Development of an improved boron-boron model. J. Appl. Phys. 1998, 84, 2471. (29) Noya, J. C.; Herrero, C. P.; Ramirez, R. Microscopic structure and reorientation kinetics of B-H complexes in silicon. Phys. ReV. B. 1997, 56, 15139. (30) Wang, L.; Clancy, P. Molecular dynamics simulations of boron diffusion in SiGe. J. Appl. Phys. 2004, 96, 1939. (31) Uttormark, M. J.; Thompson, M. O.; Clancy, P. Kinetics of crystal dissolution for a Stillinger-Weber model of silicon. Phys. ReV. B. 1993, 47, 15717. (32) Goodwin, L.; Skinner, A. J.; Pettifor, D. G. Generating transferable tight-binding parameters: application to silicon, Europhys. Lett. 1989, 9, 701. (33) Kwon, I.; Biswas, R.; Wang, C. Z.; Ho, K. M.; Soukoulis, C. M. Transferable tight-binding models for silicon. Phys. ReV. B. 1994, 49, 7242. (34) Sawada, S. New tight-binding model of silicon for theoretical studies of surfaces. Vacuum 1990, 41, 612.

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5639 (35) Bowler, D. R. A simple, effective tight-binding parameterization for Si-Ge interactions on Si(001). J. Phys.: Condens. Matter 2002, 14, 4527. (36) Yu, Q.; Clancy, P. Molecular-dynamics simulation of crystal-growth in (37) Si1-xGex/Si(100) heterostructures. J. Cryst. Growth. 1995, 149, 45. (38) Brunco, D. P.; Thompson, M. O.; Hoglund, D. E.; Aziz, M. J.; Gossman, H.-J. Germanium partitioning in silicon during rapid solidification. J. Appl. Phys. 1995, 78, 1575. (39) Srinivasan, G.; Bain, M. R.; Bhattacharyya, S.; Baine, P.; Armstrong, B. M.; Gamble, H. S.; McNeill, D. W. Tungsten silicide contacts to polycrystalline silicon and silicon-germanium alloys. Mater. Sci. Eng. B. 2004, 114-15, 223.

(40) Dong, L.; Yue, R. F.; Yan, W.; Huang, W. T.; Liu, L. T. Electrical properties of boron and phosphorus doped polycrystalline silicon germanium films. Int. J. Nonlinear Sci. Num. Sim. 2002, 3, 677. (41) Wang, L.; Clancy, P.; Thompson, M. O.; Murthy, C. S. Thermodynamic and kinetic studies of laser thermal processing of heavily borondoped amorphous silicon using molecular dynamics. J. Appl. Phys. 2002, 92, 2412.

ReceiVed for reView December 6, 2005 ReVised manuscript receiVed January 30, 2006 Accepted February 1, 2006 IE051361W