## Abstract

This paper investigates several highly accurate algorithms which can be used to calculate the morphology in a wide range of thin film process simulations, and which require minimum computational effort. Three different algorithms are considered, namely the kinetic energy corrector (KEC) algorithm, the thermal control layer marching (TLM) algorithm, and the thermal control layer marching algorithm with an incorporated KEC function (TLMC). A common characteristic of these algorithms is that they all address the recovery of the impact incident energy within the free reaction layer. However, they differ in their treatment of the thermal control layer. The TLM and TLMC algorithms consider this layer to be moveable, whereas the KEC algorithm regards it as being fixed. The advantage of employing a moveable thermal control layer is that the computational effort required to carry out simulation is reduced since the atoms lying below this layer are excluded. The relative accuracy and efficiency of the proposed algorithms are evaluated by considering their use in the simulation of the trench-filling problem associated with the damascene process. The results of the present investigation indicate that the TLM algorithm has the ability to provide an accurate morphology calculation for low and medium energy incident atoms. However, for higher incident energy impacts, the TLMC algorithm is found to be a more appropriate choice because the incorporated energy corrector function is required to remove the higher energy accumulation which occurs within the deposited atoms. Furthermore, for all three algorithms, it is noted that a suitable specification of the free reaction layer thickness is essential in determining the accuracy and efficiency of the simulation. Finally, this paper discusses the relationship between the energy absorption rate and the thickness of the free reaction layer, and presents the optimal free reaction layer thickness for different incident energy intensities.

## 1. Introduction

Molecular dynamics (MD) simulation has been widely employed in the study of thin film growth mechanisms for various deposition processes. The purpose of these studies has been to determine how the growth mechanisms influence the form of the final film morphology. MD simulation is an appropriate technique in such investigations since it has the ability to model the dynamic behaviour at a microscopic level. Therefore, it provides results which, given the nano-scale of the components involved, would be impossible to achieve experimentally. However, an inherent disadvantage of the MD approach is that the time step used in the simulation must be less than 10^{−13} s (Hockni & Istwood 1981), i.e. the cycle time of atomic vibration. Consequently, the deposition rates adopted in MD simulations are generally higher than those employed in practice, and this has sometimes caused the accuracy of the results to be called into question. Typically, the deposition rate employed in the physical vapour deposition (PVD) process is 0.05 μm min^{−1} (film thickness per min). However, the deposition rates specified in the current MD simulation, as well as in deposition simulations performed by other researchers, are of the order 10^{7}–10^{8} μm min^{−1}.

If realistic deposition rates were adopted in the simulation, the resulting computational task would be complex and would take a very long time to execute. Since the impact energy imparted by the incident atoms is only slowly transferred away from the impact zone, the deposited atoms take a long time to return to their equilibrium state. Consequently, the majority of the atoms would remain in a temporarily increased state of random thermal vibration throughout most of the simulation. A simple remedy can be used to compress the temporal domain of the simulation, i.e. since collision-induced mobility takes place over a short duration, i.e. short-term behaviour, the energy transfer which occurs after this time can be accelerated such that it takes place in a compressed time period. This allows a rapid return of the atoms to their equilibrium state. Adoption of this approach preserves the qualitative accuracy of the simulation results. Thus, a more accurate film growth morphology can be obtained by increasing the time interval between two successively deposited atoms. This allows the impact energy of the first incident atom to be transferred away from the impact zone completely before the second atom impacts the surface, and this accurately reflects the long-term recovery of the atoms to their equilibrium state. To determine a suitable deposition rate, most previous studies progressively reduced the rate in order to locate the point at which there was no further change in the final morphology for any further reduction in the deposition rate (Kelchner & DePristo 1997). This is usually referred to as the independent deposition rate morphology. Although this empirical method is convenient since it avoids the need to specify the deposition rate precisely, it suffers the drawback of requiring a greater computational effort to obtain an accurate result. Moreover, it is invalid in certain deposition situations, e.g. when high incident energy atoms are deposited at a long distance from the thermal control layer. The current paper provides a more rigorous approach to solving this problem. It introduces a kinetic energy corrector (KEC) algorithm, which checks that the residual energy of the deposited atoms is lower than a specified energy tolerance value before allowing deposition of a new incident atom to proceed. However, although it is demonstrated that the KEC algorithm improves the accuracy of the morphology calculation, it also requires an excessive computational time.

Finally, it is important to compare the methods considered in the present study with the MD simulation of long-term thin film growth addressed previously by Voter (1997). Voter considered the acceleration of the deposition process simulation from a more rigorous physical and mathematical perspective when conducting the simulation over an extended deposition time (20 μs). This situation cannot be modelled using the conventional MD time stepping approach due to the limitations of present computing ability. Voter's study also involved considerations relating to the long-term behaviours of atomic movement, including the thermally induced atomic rearrangement process, diffusion, down-stepping, etc. in addition to the short-term behaviour of kinetic energy induced atomic movement. The approaches considered in the present study, which operate from a basis of physical intuition to accelerate the deposition simulation, neglect the long-term behaviour of atomic movement. However, from an engineering perspective, the current methods nevertheless provide valuable solutions for large-scale morphology experimental comparison since the morphology changes caused by long-term atomic movement are sufficiently small that they can be ignored.

A variety of physical vapour deposition processes have been developed for the fabrication of thin films, including sputter deposition (Ju *et al*. 2002), ionized cluster beam deposition (Hwang *et al*. 2002*a*), and ionized physical deposition (Hwang *et al*. 2002*b*,*c*). In addition to thin film fabrication on a flat substrate, some processes have been developed specifically to grow thin films on corrugated substrates, e.g. the trenches or vias used in the damascene process (Su *et al*. 2002), and the V-shape substrate employed in giant magneto resistance (GMR) applications (Weng *et al*. 2002). In the case of films grown on corrugated surfaces, the objective of the damascene process is to deposit a void- or defect-free film in a specific geometry. However, the aim of the V-shaped substrate deposition process is to fabricate Co/Cu films with enhanced physical properties capable of fully satisfying the device's functional requirements for resistance and MR ratio. Successful operation of the device is related directly to the quality of the final film morphology. Therefore, an accurate prediction of the final morphology is essential if the necessary control information is to be provided to the practical manufacturing process.

The primary objective of the current investigation is to use the concept of a moveable thermal control layer as the basis for developing algorithms which are able to provide an accurate calculation of the final morphology within an MD simulation. Section 2*a* of this paper provides a brief overview of the standard MD simulation model used in the simulation of trench filling with a constant deposition rate. Section 2*b* describes how the use of the KEC algorithm improves the accuracy of the simulated film growth morphology. Section 2*c* introduces a thermal control layer marching algorithm, which eliminates unnecessary atom trajectory computations, and therefore reduces the computational time required. Section 2*d* presents a thermal control layer marching algorithm which incorporates the KEC function. Section 2*e* provides a quantitative evaluation of the relative error of each algorithm in calculating the morphology. Section 3 presents a detailed discussion of the simulation results. A comprehensive analysis of the performance of the three algorithms for various incident energy intensities is provided, and the influence of the free reaction layer thickness on the accuracy of the results and the computational time for each of the different algorithms is discussed. Section 3 also considers an appropriate choice of free reaction layer thickness for the TLMC algorithm under high incident energy conditions. Finally, §4 provides a brief summary of the principal conclusions of the present investigation.

## 2. Molecular dynamics simulation methodology

### (a) MD simulation model

Figure 1*a*,*b* present schematics of the trench used in the damascene process, and the MD simulation model, respectively. The MD simulation involves two models, namely the trench model and the deposition model. The trench model consists of a barrier layer of Ti atoms with an underlying thermal control layer. The outermost layer of the trench model is fixed in order to prevent the trench atoms from shifting. It is noted that the trench model is a three-dimensional model with a thickness of 1.5 nm, i.e. in the *y* direction. The barrier layer of Ti atoms has an hcp arrangement and a (0001) surface. Finally, approximately 1500 Cu atoms are required to fill the trench. The deposition model is mainly concerned with simulating the spatial and angular distribution of the incident Cu atoms. It uses the many body, tight-binding second-moment approximation (TB-SMA) potential model to simulate the interatomic forces, which exist within the simulated system, i.e. the Ti–Ti, Cu–Cu and Cu–Ti interactions. This potential model is given by(2.1)where *ξ* is an effective hopping integral, *r*_{ij} is the distance between atoms *i* and *j*, and *r*_{o} is the first-neighbour distance. The parameters *A*, *p*, *q* and *ξ* which are used in the simulation to describe the Cu–Cu and Ti–Ti interactions are listed in table 1, and are taken from Cleri & Rosato (1993). Rather than resorting to rigorous mathematical manipulation or theoretical analysis to determine the parameters of interatomic potential for Cu–Ti, an alternative approach is merely to take the average value of the two pure interatomic potential parameters, i.e. the relevant values of *A*, *p*, *q* and *ξ* for Cu–Cu and Ti–Ti. This procedure has been used in previous literature to obtain the interaction parameters of Sb–Si (Ladeveze *et al*. 1998) and Cu–Ru, Cu–Ir, Cu–W, etc. (Iwasaki 2000). It might be argued that the inaccuracies arising as a result of adopting this simplified approach should be reflected within the simulation. However, Cu–Ti interaction is restricted to the initial stages of the deposition process when Cu atoms are first deposited on the Ti trench. As deposition continues, the distance between the point of impact of the incident Cu atoms and the Ti trench exceeds the truncation distance of Ti. Therefore, Iwasaki (2000) states that the effects of employing the simplified approach are negligible and may be ignored. The atom trajectories are calculated in accordance with Newton's equation of motion by manipulating the potential model given in equation (2.1) to obtain the force term of the motion equation, i.e.(2.2)where(2.3)

Other than the additional treatment required to adjust the velocity of the atoms to reflect the desired equilibrium temperature, this trajectory calculation is also valid for the atoms which lie within the thermal control layer. This treatment is usually referred to as the scaling technique (Haile 1992; Rapaport 1997; Frenkel & Smit 1996). The scaled velocity is given by(2.4)where *v*_{i} is the original velocity obtained from equation (2.2), *E*_{k}^{d} is the desired kinetic energy of all of the thermal control layer atoms and is determined by the specified equilibrium temperature, and *E*_{k}^{p} is the original kinetic energy of the thermal control layer atoms based upon their velocities given in equation (2.2). In addition, the velocity of the thermal control layer is updated according to equation (2.4) at each time step.

The MD deposition simulation model also describes the characteristics of the incident atoms. These characteristics depend upon the incident energy of the atoms and their spatial distribution. In the traditional MD deposition simulation, the spatial distribution of the incident atoms is generated by a random distribution function, and the incident atoms are deposited at a constant deposition rate. The incident angle of the deposited atoms is generated by a random function and falls within a range of −5° to +5° (where the vertical *z* axis represents the zero angle reference). A cutoff angle of −5° to +5° is also applied.

### (b) Kinetic energy corrector algorithm

The adoption of a constant deposition rate in the traditional simulation may result in an accumulation of the incident energy of the impacting atoms within the deposited atoms. Depending upon the severity of this energy accumulation, this may introduce errors into the simulation results. These errors tend to become more significant as the deposited film thickness increases, particularly in the case of the deposition of high incident energy atoms. This paper introduces a KEC algorithm into the simulation computations in order to address errors of this type. The function of this algorithm is to determine whether the kinetic energy of each individual atom in the system is lower than a given kinetic energy tolerance, whose value is chosen such that the energy accumulation within the overall simulation system is minimized. Only if this condition is met, is a new incident atom introduced into the system. The KEC is expressed as(2.5)where *m* is the time step, *E*_{j}^{i} is the kinetic energy of the *j* th atom at the *i* th time instant, *n* is the total number of deposited atoms at a specific time instant, and *E*_{p} is the specified tolerance for the kinetic energy, which is determined by multiplying the kinetic energy at the equilibrium state by a scaling factor, *f*(≥1). For example, in the present simulation the equilibrium temperature is taken to be 300 K, which corresponds to an equilibrium kinetic energy, *E*_{s}, of 0.026 eV atom^{−1}. It is very difficult to attain a state of complete equilibrium within an acceptable simulation time duration, and so the value of the scaling factor is always specified to be greater than or equal to 1 since this reduces the time required for computation. Note that in the current investigation, several values of *E*_{p} are used, namely 0.0475, 0.0593, 0.0711, and 0.0830 eV atom^{−1}. The equivalent temperatures are 548, 684, 820, and 958 K, respectively. It is noted that although the temperatures are quite high, the kinetic energies corresponding to these equivalent temperatures are smaller than that associated with the incident energy, e.g. 10 eV (equivalent temperature, 1.15×10^{5} K). For the *E*_{p} values used in the present calculations, even for the highest value used, the average temperature is equal to approximately 300 K for statistic by Gaussian distribution.

From equation (2.5), it can be seen that the energy corrector function checks that the kinetic energy of each atom, *j*, averaged by the *m* time steps, satisfies the required criterion. In the calculations, the KEC is checked once for every 5 times of the atomic vibration cycle time. Each cycle time is divided into discrete 100 time steps. Hence, the KEC equation (2.5) is employed in a total of 500 time steps, i.e. *m*=500, and a temporal average is then obtained. Note that the energy of each individual atom is checked. This is necessary in order to take into account the high localized energy accumulation which tends to occur around the incident atom impact areas.

### (c) Thermal control layer marching algorithm

The thermal control layer, as depicted in figure 1*b*, in the traditional simulation deposition model (and in the KEC algorithm presented above) is stationary. This is a straightforward treatment of the thermal control layer and is consistent with the approach adopted in almost all of the simulation models presented in the literature. However, a disadvantage of this approach is that it increases the computational time required to simulate the complete outward transfer of the impact energy through the thermal control layer. This is particularly the case when the incident atoms impact at a location far removed from the thermal control layer. The thermal control layer marching (TLM) algorithm has been developed to address this problem. The TLM algorithm employs a moveable thermal control layer whose location shifts upwards during the deposition process such that it always remains at a constant distance from the surface. The adoption of a moveable thermal control layer is shown to be a convenient and effective means of removing the impact energy of incident atoms, thus preventing an accumulation of energy within the deposited atoms.

The first step of the TLM algorithm determines the locations of the surface atoms by comparing the coordination number of individual atoms with that of atoms within the inner trench. In the case of an atom lying on the surface:(2.6)where *C*_{j} is the coordination number of atom *j*, *C*_{inner} is the coordination number of an atom in the inner trench, and *n* is the total number of previously deposited atoms.

Once the locations of the surface atoms have been determined, the lower boundary of the free reaction layer can be located. The thickness of the free reaction layer, *r*_{free}, is specified by the user. The upper boundary of the free reaction layer is given by the position of the surface atoms, while its lower boundary is located at a distance *r*_{free} below the surface, i.e.(2.7)where *r*_{free}^{s} is the thickness of the surface atoms layer, *r*_{free}^{i} is the thickness of the inner free reaction layer; in which ‘inner’ means lying beneath the layer of surface atoms. The surface atoms layer and the inner free reaction layer constitute a free reaction layer, as described by the left side of equation (2.7). *r*_{free}^{ev} is the minimum thickness of the free reaction layer which gives optimal simulation results for a specified incident energy.

The thermal control layer lies below the free reaction layer, and so once the location of the lower boundary of the free reaction layer has been established, the lower limit of the thermal control layer can also be located, i.e. at a distance of *r*_{thermal} below the last layer of atoms within the free reaction layer. Finally, a fixed layer is located beneath the thermal control layer, which prevents any further onward energy transfer.

Once the location of the total reaction layer, which includes the surface layer, the free reaction layer, the thermal control layer and the fixed layer, has been determined, the position of the non-reaction layer can be identified. This layer is composed of the deposited atoms which lie beneath the fixed layer. As the deposition process continues, the total reaction layer progressively moves upwards, and the number of atoms within the non-reaction layer increases accordingly. The treatment of these atoms is different from those in the total reaction layer, i.e. their trajectories (given by equation (2.2) or (2.4)) are no longer computed. They are excluded for the simple physical reason that the impact energy of subsequent atoms deposited on the surface is too weak to greatly influence their motion. Excluding these atoms results in a considerable saving of computational effort. A further advantage is that the computational time varies linearly with the total number of deposited atoms, rather than exponentially, as is the case for simulations which use a fixed control layer. Finally, the moveable thermal control layer is updated at a time period corresponding to the time required by the subsequently deposited atoms to form at least one complete layer over the configuration existing at the previous update. From the previous description of the moveable structure, it is clear that if an atom moves outside the thermal control layer, it will move into the free reaction layer, where it can move freely and interact with other atoms.

### (d) Thermal control layer marching with kinetic energy corrector algorithm

This section introduces the TLMC algorithm, which is a thermal control layer marching algorithm with an incorporated KEC function. Basically, it represents the integration of the KEC and the TLM algorithms. However, the scope of the KEC algorithm is slightly different in this implementation since it only corrects the kinetic energy of atoms lying within the free reaction layer and in the thermal control layer, whereas it will be recalled that in §2*b*, the algorithm considered all of the deposited atoms. As stated previously, high incident energy impacts tend to result in an accumulation of energy within the free reaction layer. Incorporating the KEC algorithm into the TLM algorithm allows this energy accumulation to be removed from the free reaction layer and transferred into the thermal control layer, and therefore the accuracy of the morphology calculation is improved.

### (e) Quantitative evaluation of morphology difference

Figure 2 presents a comparison between the current morphology (as calculated by the TLM or TLMC algorithm) and a ‘reference’ morphology, which is obtained from the KEC algorithm. The difference between any two filling morphologies is established by counting the number of Cu atoms which lie in the shaded region bounded by the two outer extremities of the individual morphologies. The relative error of the TLM and the TLMC algorithms is then determined by dividing the number of Cu atoms within this shaded area by the total number of atoms in the current morphology.

## 3. Result and discussion

### (a) Traditional treatment of a constant deposition rate

#### (i) Kinetic energy variation

Figure 3 shows the average kinetic energy history of the deposited Cu atoms for an incident energy of 5 eV atom^{−1}, and a constant deposition rate of 1.25 atoms ps^{−1}. The figure shows quite clearly that the average kinetic energy of the deposited atoms increases monotonically over time. In the case of the traditional deposition process, which uses a constant deposition rate, the impact energy of an incident atom will result in a sudden increase in the kinetic energy of the deposited atoms, corresponding to the peaks visible in figure 3. The impact energy is then transferred into the thermal control layer via a series of inter-atomic collisions. Depending upon the intensity of the incident energy, the time required to remove sufficient incident energy to allow the impacted atoms to return to their thermal equilibrium state may be long. If the deposition rate is too rapid, there is insufficient time for the incident energy to be fully transferred into the thermal control layer before the next incident atom impacts on the surface arrives, and so an energy accumulation occurs within the free reaction layer.

### (b) Kinetic energy corrector algorithm

Figure 4 explores the relationship between the relative error in morphology and the number of deposited atoms for three different kinetic energy corrector tolerances. The figure indicates the error relative to a morphology calculated using a reference kinetic energy tolerance of 0.0475 eV atom^{−1}. Note that this value is higher than the ‘true’ kinetic energy tolerance of 0.026 eV atom^{−1}, which corresponds to the energy of the atoms in their thermal equilibrium states. Thus, this result can just be represented as an approximate measure of the error, based on comparsion to the 0.0475 eV atom^{−1}. The higher tolerance value is chosen in order to avoid the excessive computational time and complexity which would otherwise result.

Figure 4 plots the relative error in morphology against the total number of deposited atoms for kinetic energy corrector tolerances of 0.0593, 0.0711 and 0.0830 eV atom^{−1}. It can be seen that the relative error resulting from an energy tolerance of 0.0830 eV atom^{−1} is slightly higher than the error associated with tolerances of either 0.0593 or 0.0711 eV atom^{−1}. It is also noted that the use of the latter two tolerance values produces very similar relative error results. Since the magnitude of the relative error is small in these two cases, the accuracy of the morphology calculation generated using these two tolerances may be taken to be almost as precise as that obtained using the smaller tolerance value of 0.0475 eV atom^{−1}. In order to avoid the excessive computational effort associated with a smaller kinetic energy tolerance value, the following discussions relating to the KEC and TLMC algorithms adopt a kinetic energy tolerance of 0.0711 eV atom^{−1} unless stated otherwise.

### (c) Thermal control layer marching algorithm

#### (i) Low incident energy and medium incident energy

Figure 5 indicates the location of the thermal control layer at three instants within a simulation of the deposition process using the TLM algorithm, and an incident energy level of 5 eV. The figure shows that the deposited atoms gradually accumulate along the trench wings on either side of the trench opening. In a traditional simulation, which adopts a fixed thermal control layer, this accumulation makes it increasingly difficult to remove the incident energy in a tolerable time duration because the distance between the point of impact of the incident atoms and the control layer increases as more atoms are deposited. However, as may be seen in figure 5, the TLM algorithm overcomes this limitation by shifting the location of the control layer upwards as deposition of new incident atoms continues. This approach ensures an effective absorption of the impact energy, and prevents an accumulation of the incident energy within the deposited atoms.

Figure 6 shows the relative error of the morphology, as calculated by the TLM algorithm and by the traditional deposition simulation, compared to a reference morphology obtained from the KEC algorithm based upon an energy tolerance of 0.0711 eV atom^{−1}. In the simulation, a three-layer thickness of the free reaction layer is used. The results presented for the low incident energy of 5 eV atom^{−1} and the medium incident energy of 10 eV atom^{−1} are based upon deposition rates of 1.25 and 0.9 atoms ps^{−1}, respectively, for both types of simulation. Observation of the figure indicates that the TLM algorithm provides a more accurate calculation of the morphology than the traditional treatment.

For an incident energy of 5 eV atom^{−1}, it can be seen that the relative error of the TLM algorithm decreases continuously as deposition continues. By contrast, the relative error of the traditional treatment first falls, and then rises as more atoms are deposited. This increase can be attributed to the accumulation of energy within the deposited atoms. A similar trend in the morphology error is noted for the traditional simulation process using a medium incident energy of 10 eV. However, in this case, it is noted that the magnitude of the relative error is higher, and that an abrupt increase in the error occurs at an elapsed time of 3000 ps. It may be seen in figure 7*a* that atoms accumulate along the trench wings during the traditional simulation. As has been discussed previously, this results in an accumulation of energy. The atoms in these regions are less able to resist the impact of subsequent incident atoms, and accordingly, the atom clusters tend to collapse when impacted by subsequent incident atoms. As evidenced by the abrupt increase in relative error shown in figure 6, this greatly increases the inaccuracy of the morphology calculation. A comparison of the elapsed time at which this abrupt increase takes place, indicates that the atom clusters collapse earlier for a higher incident energy.

Prior to an elapsed time of 1500 ps, the magnitude of the relative error of the TLM algorithm with an incident energy of 10 eV atom^{−1} is approximately 1% higher than when the lower incident energy of 5 eV atoms^{−1} is adopted. However, as deposition continues, the discrepancy between the two sets of results decreases continuously, and the results shows that the final relative error of the TLM algorithm is approximately 2% for both low and medium incident energy atoms.

Figures 7 and 8 present a comparison between the calculated morphologies resulting from three different treatments (traditional, KEC, and TLM) at incident energies of 5 and 10 eV atoms^{−1}, respectively. Note that the total number of deposited atoms is the same in each figure. It is immediately obvious that a touching of the two overhanging clusters only occurs within the traditional simulation (i.e. figures 7*a* and 8*a*). Furthermore, although not presented in the current paper, the transient morphology evolution shows that the deposited atoms above the trench opening exhibit a softer material behaviour in the case of the traditional treatment, and that previously deposited atoms readily ‘collapse’ when impacted by subsequent incident atoms. Finally, it is noted that there is little difference between the KEC algorithm results and those generated by the TLM algorithm at either level of incident energy.

#### (ii) High incident energy

In order to investigate the influence of the free reaction layer thickness on the accuracy of the morphology calculation at a high incident energy of 20 eV atom^{−1}, figure 9 shows the relationship between the relative morphology error and the number of deposited atoms for the TLM and TLMC algorithms with different free reaction layer thicknesses (i.e. 2-, 3-, 4- and 5-layer thicknesses). The figure also plots the results for the traditional simulation for comparison purposes. The results for the traditional treatment and the TLM algorithm are both based upon a deposition rate of 0.67 atom ps^{−1} while those of the TLMC algorithm are based upon a smaller average deposition rate of 0.2–0.3 atom ps^{−1}.

Considering the case of the TLM algorithm, figure 9 shows that the morphology calculation becomes less accurate as the thickness of the free reaction layer increases. However, even with the minimum free reaction layer thickness of 2 atoms, the figure shows that for a high incident energy of 20 eV atom^{−1}, the TLM algorithm is unable to match the 2% accuracy yielded for low and medium incident energy atoms. From figure 9, it can be seen that a five-layer free reaction thickness results in the least precise morphology calculation. This may be explained in the same way as for the case of the traditional simulation, i.e. energy tends to accumulate within the deposited atoms as the distance between the surface atoms and the thermal control layer increases, and this results in a collapse of the overhanging atom clusters.

In order to address the problem of accumulated energy within the free reaction layer, the kinetic energy corrector function is incorporated into the TLM algorithm. The resulting algorithm, referred to as the TLMC algorithm, prevents deposition of an additional atom until the incident energy imparted by the previous atom has been fully transferred outwards through the thermal control layer. In this way, the accumulation of energy within the deposited atoms is reduced, and the accuracy of the morphology calculation is improved. Referring to figure 9, it is noted that the accuracy of the morphology calculation using the TLMC algorithm improves as the thickness of the free reaction layer increases. Although this trend is the exact reverse of that noted for the TLM algorithm, it is not unexpected since the greater layer thickness used in the TLMC algorithm is similar to that used in the accurate KEC algorithm, in which the free reaction layer thickness increases continuously. Observation of figure 9 shows that incorporating the kinetic energy function into the TLM algorithm is successful in reducing the relative error of the morphology calculation to less than the minimum value of 2% achieved by the TLM algorithm for low and medium incident energies.

Close inspection of figure 9 shows that the TLMC algorithm is less accurate than the TLM algorithm when a free reaction layer thickness of 2 or 3 is adopted. In other words, it would appear that there is no benefit in incorporating the kinetic energy check function within the TLMC algorithm when the free reaction layer is thin. In order to investigate this issue further, figures 10 and 11 plot the decay of the accumulated energy, which is calculated by summing the kinetic energy of each individual atom and then subtracting the total equilibrium energy of these atoms at 300 K, for the TLM algorithm and the TLMC algorithm with different free reaction layer thicknesses.

Figure 10 plots the decrease in the accumulation energy over the first three cycles of the TLM algorithm for four different free reaction layers. The results are plotted against the energy transfer achieved over the first cycle of the TLMC algorithm for comparison purposes. The latter algorithm adopts a seven-layer thickness since this thickness is expected to give results which are as accurate as those provided by the KEC algorithm. It can be seen that the thicker free reaction layers are less efficient in removing the impact incident energy. Furthermore, a comparison of the results suggests that the thicker layers become increasingly inefficient compared to the TLMC algorithm with a seven-layer thickness as the number of cycles increases. Regarding the TLM algorithm, the figure shows that the degree of energy accumulation is less when the layer thickness is thinner. This observation is in full agreement with the results presented in figure 9, which showed that the TLM algorithm yields more accurate results for thinner free reaction layer thicknesses.

Figure 11 plots the incident energy transfer against the elapsed time from impact for the TLMC algorithm with different free reaction layer thicknesses. The time required for atoms within the free reaction layer to return to their equilibrium state after impact depends upon the thickness of the layer. Although the results shown in figure 11 were all generated using a computational time sufficiently long to allow atoms within even the thickest reacted layer to reach their equilibrium state, for reasons of convenience, the figure truncates the decay plots, and the time scale refers only to a free reaction layer of two atom thickness. It can be seen that the free reaction layer with a layer thickness of two deviates most significantly from the solid line representing a seven-layer thickness, which was used for reference purposes in figure 10. Inspection of the slope of the various decay curves indicates that the two-layer thickness results in the greatest energy absorption rate. However, the results presented in figure 9 suggest that the effect of this high absorption rate is to introduce large errors into the morphology calculation. In order to demonstrate that incident energy absorption is strongly dependent upon the free reaction layer thickness, figure 11 also plots the kinetic energy variation for the TLMC algorithm with a 13-layer free reaction layer thickness. An excessive thickness results in a slower deposition rate, and requires a longer computational time since it takes longer for the kinetic energy of the deposited atoms to fall below a given kinetic energy tolerance.

In figure 9, it was shown that the TLMC algorithm with a five-layer thickness reacted layer thickness provides the most accurate morphology calculation. Therefore, referring back to figure 11, it would appear that the slope of the energy decay curve of the free reaction layer with a five-layer thickness represents an optimal absorption rate. Absorption rates in excess of this value will result in inaccurate morphology calculations, while less rapid absorption will yield more accurate results. To illustrate this point, figure 12 presents the reduction in the total number of atoms whose energy exceeds 0.0711 eV atom^{−1} over elapsed time from impact for the TLMC algorithm with various free reaction layer thicknesses. Note that a reference energy value of 0.0711 eV atom^{−1} is deliberately adopted in order to distinguish the number of atoms which are activated to a higher level of energy than the energy tolerance of 0.0593 eV atom^{−1} used in the corrector algorithm. Observation of the figure shows that the decay curves tend to approach each other as the number of layers increases. It is also noted that the number of activated atoms for a two-layer thickness is far smaller than in the case of the thicker layers in the immediate post-impact stage. This observation can be attributed to the high absorption rate of the thin free reaction layer, which was discussed previously. Figure 12 suggests that specification of a suitable layer thickness within the TLMC algorithm should consider the total number of atoms which will become activated as a result, and should endeavour to activate as many atoms as are generated by a highly accurate algorithm such as the KEC algorithm. Finally, it should be noted that although the use of a thicker layer does not greatly increase the number of atoms with a higher activated energy, the number of atoms within the inner regions of the deposited atoms which have a smaller activated energy does increase. However, the intensity of their activated energy is insufficient to cause migration from their equilibrium position, and therefore the accuracy of the morphology calculation is unaffected.

From the preceding discussions, it is clear that two factors are important in determining the accuracy of the morphology calculation. In the case of the TLM algorithm, the energy accumulation in the free reaction layer is key, while in the case of the TLMC algorithm, the absorption rate of the thermal control layer is the overriding concern. The TLM algorithm does not include a kinetic energy corrector function to allow the kinetic energy to be transferred prior to the arrival of a subsequent incident atom, and so the resulting energy accumulation determines the accuracy of the results. In the TLMC algorithm, the kinetic energy is transferred outwards until the energy of the deposited atoms reaches the specified tolerance value before a further incident atom is introduced into the system. Therefore, the absorption rate governs the accuracy of the morphology calculation. It has been shown that the TLM algorithm is most accurate when using a two-layer thickness of the free reaction layer since the degree of energy accumulation is minimum in this case. Conversely, use of a two-layer thickness in the TLMC algorithm provides the least accurate set of results due to the excessive energy absorption rate.

Figure 13 compares the morphology calculation for four different algorithms (i.e. traditional, KEC, TLM and TLMC) at a high incident energy of 20 eV atom^{−1}. The total number of deposited atoms is the same in all four cases. Note that a five-layer free reaction layer thickness is used in the TLM and TLMC algorithms. A comparison of figure 13*a*,*b* shows that in the case of the traditional simulation, the accumulated energy within the atoms on the trench wings causes these clusters to collapse upon impact from subsequent incident atoms, and to flow towards the trench cavity. This results in the trench cavity being overfilled in comparison to the results predicted by the KEC algorithm. With a five-layer thickness of the free reaction layer, the results for the TLM algorithm shown in figure 13*c* again show that the trench is overfilled. Finally, a comparison of figure 13*b*,*d* show the morphology calculations to be very similar when using the KEC algorithm or the TLMC algorithm with a five-layer free reaction layer thickness.

### (d) Computational time comparison for different algorithms

Figure 14 investigates the efficiency of the various simulation approaches by considering the relationship between computational time and the number of atoms deposited with an incident energy of 20 eV atom^{−1} under a variety of different conditions. These conditions include different thickness of the free reaction layer for the TLM and TLMC algorithms, and different kinetic energy tolerances for the KEC algorithm. The computational time for the traditional simulation approach is also included for comparison purposes, and is seen to lie between the KEC algorithm and the TLM algorithm. Observation shows that the TLM algorithm requires the shortest computational time, while the KEC algorithm is far more demanding in its computational requirements. Furthermore, the computation time required for the TLMC algorithm is always greater than the time required for the TLM algorithm, regardless of the free reaction layer thickness. This is to be expected since the TLMC algorithm includes the additional kinetic energy corrector function. Additionally, it can be seen that the computational times consumed by the KEC and TLM algorithms increase as either the given kinetic energy tolerance is reduced, or as the thickness of the free reaction layer increases, respectively. In the case of the KEC algorithm, this observation is explained by the fact that it takes longer for the kinetic energy of the activated deposited atoms to fall below the specified kinetic energy tolerance when the value of this tolerance is lower. Regarding the TLM algorithm, the computational time increases with the thickness of the free reaction layer since the number of atoms which must be considered within the computation is higher.

Considering the relationship between the computational time and the total number of atoms deposited, it is observed that the time required for the KEC algorithm is an exponential function of the number of deposited atoms, whereas the relationship for the TLM algorithm is linear. The linear increase in the computational time is a particular characteristic of the TLM algorithm. In this algorithm, the total number of atom trajectory calculations remains almost the same at any time instant because the thermal control layer morphology advances progressively as deposition continues. Conversely, in the case of the traditional treatment and the KEC algorithm, the thermal control layer is fixed, and so the number of atom trajectory calculations required increases continuously as more atoms are introduced into the system.

Although the TLM algorithm is convenient from the point of view of its computational efficiency, the results have shown that it is difficult to achieve a relative morphology error lower than 2%. Conversely, the KEC algorithm, which proves to be highly accurate in its calculation of morphology results, is inefficient in its use of computing resources. In fact, analysis of the results presented in figure 14 indicates that the KEC algorithm requires approximately four times as long as the traditional method to complete processing.

Comparing the relative efficiencies of the traditional simulation process and the TLMC algorithm reveals an interesting result. Figure 14 shows that the computational time for the TLMC algorithm is greater than for the traditional method when the thickness of the free reaction layer is either 2, 3 or 5 atoms, but is less than the traditional simulation run time when the thickness of the layer is 4 atoms.

Figure 14 shows that the longest computational time for the TLMC algorithm occurs with a two-layer free reaction layer thickness, and that an increase in the layer thickness to three atoms reduces the computational time. In other words, the computational time reduces as the layer thickness increases. However, the results for a free reaction layer thickness of four and five atoms show that the computational time required for the five-layer thickness is actually longer than for the four-layer thickness. This unanticipated behaviour is associated with the degree of recovery within the thermal control layer. When the free reaction layer is very thin, the high energy of the incident atoms may cause the deposited atoms to penetrate through the free reaction layer and to activate the atoms in the thermal control layer into a state of violent motion. In the part of the free reaction layer comprised of deposited atoms, the excited atoms are able to reach their equilibrium positions ‘freely’ and quickly, and are governed by the equation of motion given in equation (2.2). However, in the thermal control layer, the excited atoms are subjected to the additional scaling law, and so the elapsed time of the atoms' recovery to their equilibrium positions is longer. This is due to the usage of a global kinetic energy scaling factor to scale the velocity of a local atom, which is a violent motion state. This gradual scaling procedure makes the recovery of the atom to its equilibrium position very slow. A detailed investigation of the kinetic energy of the thermal control layer atoms finds that there exists an abrupt increase in the atoms' energy when two- or three-layer thickness free reaction layers are used in the TLMC algorithm. Therefore, the computational time is longer when these layers are used in the TLMC simulation, as shown in figure 14. The fact that the computational time is greater when the TLMC algorithm adopts a five-layer free reaction thickness rather than a four-layer thickness is in accordance with the observations made for the TLM algorithm, where the computational time was seen to increase with increasing layer thickness. It is noteworthy that the computational time required for a four-layer thickness of the free reaction layer in the TLMC algorithm is shorter than that of the traditional treatment. In other words, the TLMC algorithm is able to satisfy the twin requirements of accuracy and efficiency.

Although the results are not presented within this current article, a comparison of the computational time for the traditional method and the TLM algorithm at incident energies of 5 and 10 eV atom^{−1} revealed results similar to those presented in figure 14 for an incident energy of 20 eV atom^{−1}.

### (e) Optimal free reaction layer thickness

Referring back to the discussions of figure 12, it can be seen that the TLMC algorithm yields an accurate morphology calculation when the thickness of the free reaction layer is such that the number of activated atoms is close to the number of atoms represented by the thick solid line (i.e. the TLMC algorithm with a thirteen-layer free reaction layer thickness). This search principle may be repeated for different values of incident energy in order to determine the approximate thickness of the free reaction layer which results in the majority of atoms being activated by the incident impact energy. This allows the relationship between the optimal free reaction layer thickness and the incident energy intensity to be determined, as shown in figure 15. It can be seen that the optimal thickness of the free reaction layer increases as the incident energy increases. Furthermore, the curve fitting result shows that the relationship between the optimal thickness of the free reaction layer and the incident energy is very nearly linear.

## 4. Conclusion

This paper has presented the use of thermal control layer marching algorithms in the morphology calculations of a MD simulation of the damascene process. Two marching algorithms have been considered, namely the basic TLM algorithm, and the TLMC algorithm, which is essentially an integration of the TLM algorithm and the incorporated kinetic energy corrector function. The current results, relating to energy levels of 5-, 10-, and 20 eV atom^{−1}, indicate that the TLM algorithm provides accurate results when the incident energy of the deposited atoms is low. In this case, the error introduced by energy accumulation is very small and can be neglected. However, for higher incident energies, the addition of a kinetic energy corrector function is required to address the energy accumulation which arises within the free reaction layer. Although it has been noted that the use of a thinner free reaction layer in the TLM algorithm can effectively, and quickly, absorb the incident energy, it has also been observed that the relative morphology error cannot be reduced below 2%. The present results seem to indicate that only the TLMC algorithm has the ability to meet the twin requirements of computational efficiency and morphology accuracy.

If these two requirements are to be satisfied, the correct choice of free reaction layer thickness is crucial since this determines the energy absorption rate of the thermal control layer. If the thickness is too small, the energy absorption rate is too high and this introduces errors into the morphology calculation. Conversely, if the layer thickness is too great, then the time to generate an accurate result is excessive. The optimal thickness of the free reaction layer is defined as the thickness for which the total number of activated atoms is close to the number obtained when using a very thick free reaction layer in the TLMC algorithm, as shown in figure 12, or when using the KEC algorithm.

Although trench filling of the damascene process has been used as the example in this present investigation, the proposed algorithms can also be applied to the simulation of thin film deposition simulation on a flat substrate. It is worth mentioning that the advantage of the constant computational task associated with the TLM algorithm renders this particular algorithm suitable for use in the calculation of morphologies for large-scale deposition problems. However, it is clear that the inaccuracy of this algorithm for deposition under high incident energy conditions must be addressed as a matter of immediate priority.

## Acknowledgments

The authors gratefully acknowledge the support provided to this research by the National Science Council, Republic of China, under grant nos. NSC90-2212-E006-161.

## Footnotes

As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs.

- Received January 12, 2004.
- Accepted May 12, 2005.

- © 2005 The Royal Society