Emergence of non-linear effects in nanocluster collision cascades in amorphous silicon

Cluster ion beams create considerably more damage in silicon and other substrates and eject more material than single ions that deposit at the same kinetic energy on the substrate. The mechanisms that causes the non-linear growth of damage and sputtering are interesting from the point of view of both basic materials research and industrial applications. Using classical molecular dynamics, we analyse the dynamics of collision cascades that are induced in amorphous silicon by small noble gas nanoclusters. We show that the sputtering and other non-linear effects emerge due to the high-energy density induced in a relatively small region in the substrate during the cluster stopping phase and because of the timing of consequent events that dissipate the energy over a larger volume of the substrate.


Introduction
Atomic cluster beams are nowadays used in many important applications including surface processing and secondary ion mass spectrometry [1,2]. Therefore, it is important to understand in detail the atomic level mechanisms of cluster stopping in the target material and how the impacts of clusters permanently change the surface and bulk structure. Among the various phenomena induced by nanocluster bombardment of solid materials, the ejection of atoms during the bombardment, which is usually called sputtering, is one of the most interesting and most studied phenomena due to its applications and also because the sputtering yields can be measured.
Non-linearity refers to the observations that various quantities like sputtering yields and dimensions of the damaged region do not scale linearly with the number of atoms N (nuclearity) in the impacting nanocluster, when the impact velocity is constant. If the constituents of a nanocluster interacted independently with the target substrate at the same velocity, the resulting sputtering yield, for example, would be the sum of the yields induced by the constituents. The challenge for cluster impact physics is to understand the mechanisms that lead to the non-linear behaviour.
Non-linear increases of damage and sputtering are observed in many experiments. The main results are summarized in a recent review [1], so we mention here only some examples of the results. It is now clear that Au clusters induce non-linear damage in Si [3], although in the first experiments these effects were not seen [4]. Au cluster ions produce 50% more damage in crystalline Si than monomer bombardment [3]. In addition, it has been experimentally verified that the number of displaced Si target atoms in a C cluster impact does not increase linearly with the nuclearity at 6 keV cluster −1 and when N = 1-10 [5]. It has also been found experimentally, that both the range and damage produced by B N cluster ion implantation increases with cluster nuclearity N [6]. The strength of non-linear effects also depends on the ratio between the masses of projectile and target atoms [7], which indicates that the effects are related to dynamics of binary collisions.
The aim of this study is to examine the emergence of non-linear effects in nanoclusterinduced collision cascades in amorphous Si using classical molecular dynamics (MD). A particularly interesting energy and nuclearity region is where the change from linear to nonlinear collision cascades occurs. According to our MD simulations, this occurs at small cluster sizes (N < 20) and at energies around 1 keV atom −1 . We have chosen amorphous Si as the model system, because it can be reasonably well modelled in MD and collision cascades appear in their simplest form in a random mono-atomic substrate. For example, channelling of the cluster constituents affects cluster stopping in fcc metals [8], which makes the analysis of the cascade formation more complicated than in amorphous Si.
Already in 1998, Ihara et al [9] investigated the cluster implantation mechanism in Si(100) with MD. They found that the cascade temperature increases with cluster nuclearity. Because of the large computer time required by the simulations they only did a qualitative analysis of collision cascades based on single simulation runs. They also observed shock waves and the amorphization of crater walls, which is typical for cluster impacts on Si. These phenomena have been observed and analysed in many other simulation studies since then, for example in [7], [10]- [16]. High-energy densities induced by the cluster impacts are thought to be the reason for the non-linear behaviour. However, no systematic study of the cluster energy deposition mechanisms has been published. In this paper, we present a detailed study of the knock-on atom trajectories, which reveals the mechanisms that produce the non-linearity. 3 The rest of this paper is organized as follows. The simulation methods and data analysis are described in section 2. The results of the simulations are discussed in section 3. In section 4, we introduce a continuum model for cascade development and show that its predictions agree qualitatively with the simulations. Finally in section 5, we compare the results and conclusions with the experimental and theoretical results obtained elsewhere.

Methods
The simulated amorphous Si substrate was a built by copying a 3 × 3 × 3 nm cube of amorphous Si that was optimized with the algorithm of Wooten, Winer and Weaire (WWW) [17,18], and then relaxed with MD to have an amorphous structure, where 97% of the Si atoms in the bulk are fourfold coordinated. Our test simulations using target structures annealed from liquid to amorphous state using only MD show that the structure is not relaxed as well as using the WWW algorithm, and voids are left inside the material. In the simulations of unoptimized Si systems, cluster impacts often induce structural changes that are not present when an optimized WWW structure is used as a target. This can prevent the precise calculation of collision cascade dimensions.
For each combination of definite energy and cluster size, 27 simulations were run to get statistically significant results. The initial positions and orientations of the clusters were varied randomly between the simulations. The size of the rectangular substrate was 20 × 20 × 20 nm. The simulation time was 30-40 ps depending on the cluster energy, which is long enough to reach the phase when the sputtering is over. Berendsen temperature control was used to cool the sides and the bottom of the simulation cell to 300 K [19]. The thickness of the cooled region was 1.0 nm. This technique to prevent shock waves reflecting back and disturbing the cascade is discussed in [20] and other methods used in the simulations are described in detail in [21]- [23].
Positions, energies and velocities of atoms were saved every femtosecond during the first 500 fs, and the knock-on atoms were identified from these data. An atom was labelled as a primary knock-on atom (PKA), if a cluster atom was closer than 3 Å at the moment when the kinetic energy of the atom exceeded 0.1 eV for the first time. Similarly, an atom was identified as a secondary knock-on atom (SKA), if there was a PKA but no cluster atoms in the neighbourhood. If neither cluster atoms nor PKAs were within the 3 Å radius from the atom, and the atom gained higher energy than 0.1 eV, it was considered a tertiary knock on atom (TKA). Therefore, the tertiary collision process covers several collision steps in this study. The rather low energy threshold 0.1 eV ensured that all knock-on atoms were detected. The 3 Å neighbourhood radius is considerably longer than the repulsive interaction distance between Ar and Si atoms. This value was chosen because some knock-on atoms are induced in many-body collisions where they gain energy indirectly from the colliding primary atom. If a shorter neighbourhood radius were used, a considerably smaller number of PKAs would be detected. Due to the many-body collisions, it is not possible to detect precisely which atoms are SKAs and which are primary atoms that gain their energy in primary many-body collisions. The impossibility of distinguishing the collision steps after the secondary collisions is the reason for the decision to bundle all later collisions in the same tertiary collision category. Immediately after collisions, some knock-on atoms have only positive potential energy, because they are pushed near other atoms and are not yet moving. Therefore, the kinetic energies used in the analysis are the maximum kinetic energies that the atoms have during their displacement. 4 The interatomic potential used in the simulations can considerably affect the results of cluster bombardment [24]. We chose the environment-dependent interatomic potential (EDIP) [25], because it provides single ion sputtering yields which agree well with the experimental yields at the energies used in this study, and because in cluster impact simulations it produces craters of average sizes and forms compared to the other interatomic potentials in the comparison [24]. In addition, the melting point of EDIP Si is only 10% below the experimental value, which is better than values given by some commonly used DFT methods [26], and which ensures that the description of displacement cascade development is realistic. In addition, it describes the crystalline phase, amorphous phase and point defects very well. However, the average coordination of the liquid EDIP Si is 4.5 while it is 6.5 experimentally [26]. This may be a problem in high energy (E > 10 eV atom −1 ) cluster impact simulations, where the behaviour of the liquid flows are important. The EDIP potential has a fairly similar functional form for the two-body interaction as the Stillinger-Weber potential [27,28], but it is modified according to the local coordination of the atoms [25]. This flexibility makes it a good choice for impact simulations where both solid and liquid phases as well as the transitions between them must be described. In addition to the EDIP potential, a short-range repulsive potential [29] was smoothly joined to the EDIP potential and used to prevent high energy Si atoms moving too close to each other. Electronic stopping was applied as a non-local frictional force to all atoms having a kinetic energy larger than 10 eV [30,31]. According to our tests, this limit does not affect the results presented in this study, although it has been shown recently that electronic stopping is an important effect also at low particle velocities [32].
The Ar clusters were prepared using a Lennard-Jones potential [33]. Binding energies of Ar-Ar and Ar-Si are weak compared with that of Si-Si, so only repulsive potential [29] was applied to Ar atoms during the simulations. In spite of the lack of attraction between the cluster atoms, the prepared clusters completely maintained their spatial configurations before they arrived at the surface, because the repulsive potential mainly comes into play at short ranges. Figure 1 shows how the sputtering yield per cluster Y increases non-linearly with the nuclearity N , when the cluster size is larger than seven atoms and the impact energy is greater than 500 eV atom −1 . At lower energies and smaller cluster sizes, Y is almost linearly related to N . The deviations from linear behaviour are due to the random variations in the simulations. In general, the yields continue their non-linear growth also when the energy per atom E increases beyond the highest energy used in this study. At very high impact energies, yields start to decrease, because high energy clusters penetrates so deep into the substrate that the main collision cascade does not reach the surface any more [8,24]. Therefore, we conclude that the change from linear to non-linear sputtering regime occurs at the energies and nuclearities shown in figure 1, and the sputtering mechanisms must change in this region.

Results
Before drawing any conclusion from MD simulations, one should verify that the results are in reasonable agreement with experiments. The EDIP potential used in the simulations gives rather good agreement with the experimental sputtering yields for single Ar ions [24]. However, that comparison was made using Si(111) substrate instead of the a-Si substrate used in this study. For 1 keV Ar ions and Si(111) substrate, the sputtering Y = 0.75 [24], while in this study the a-Si substrate gives a lower yield Y = 0.3. Rubio et al [34] have also compared simulated  1 keV Ar ion yields in Si(100) and a-Si. They report values Y = 0.39 for a-Si and Y = 0.46 for Si(100) using the Stillinger-Weber potential. In general, the simulated yields for a-Si are lower than the simulated and experimental yields for c-Si. In conclusion, the magnitude of the sputtering yields (figure 1) is consistent with other MD simulations and in modest agreement with experiment. In addition to sputtering, cluster impacts induce craters on the surface, if the impact energy is high enough [24]. One indication of cratering is a high number of substrate atoms ejected above the original surface plane. The right frame of figure 1 shows that the ejection increases non-linearly with N , and is proportional to N 2 in the non-linear sputtering regime (N > 7, E > 500 eV). Clear and regular-shaped craters are observed in the visualizations of the 2 keV Ar 16 simulations, whereas the craters become more irregular and smaller with decreasing N . Thus, the change from linear to non-linear sputtering is related to displacement of a large number of atoms, which also leads to crater formation. In other words, the sputtering yield is almost linearly proportional to N until the energy density in the collision cascade is high enough to induce a coherent flow of atoms to vacuum. The flow begins at lower energies and nuclearities than non-linear sputtering, indicating that non-linear sputtering does not begin until the flow is intensive enough to eject the atoms to vacuum. At lower-energy densities, only individual highenergy atoms are able to leave the surface. Next, we analyse the collision cascade development which leads to this transition from surface evaporation to the coherent flow of atoms. The flow of atoms is observed also in other substrates [35,36].  figure 2, the initial position distributions of the displaced atoms are smaller than the corresponding TKA distributions. Thus, the atoms in the outer shell of the collision cascade receive kinetic energy, but the original atomic structure is not significantly modified. Figure 3 shows that the large clusters lose a larger portion of their energy at the surface than the small ones, but otherwise the stopping of the cluster atoms is very similar at various nuclearities. However, a small clearing-the-way effect [37,38] increases the projected ranges of the cluster atoms. This is seen in figure 3, where the distribution of final positions of the cluster atoms just after the stopping phase is shown to become wider and deeper with increasing N . The  Note that the depth scale is logarithmic in (a) in order to show the differences of energy deposition at the surface. The statistical variation makes it impossible to define the lateral width precisely. The positions are detected just after the stopping phase is over. Some cluster atoms can sputter after that among the flow of atoms to vacuum. enlargement of the distribution is around 10% between N = 1 and 16, which is not enough to explain the non-linearities alone. The explanation to this effect can be seen in the visualizations of the simulated events. A frontier of PKAs clears the way for the cluster atoms, thus the  atoms of large clusters penetrate deeper into the substrate. This clearing-the-way effect is a consequence of the classical collision dynamics of the particles and is not based on the coherent electronic effects described in the theories of swift cluster stopping [39].
The kinetic energies of PKAs, SKAs and TKAs increase with N ( figure 4). When the cluster atoms move in the material close to each other, they occasionally collide simultaneously with the same Si atoms, which on average gain more energy than in the case when only twobody collisions occur. Naturally, the SKAs and TKAs also gain more energy with the primary atoms. Figure 4 also shows that there are less low energy PKAs and SKAs at high nuclearities. All Si atoms in the vicinity of the track of a large cluster become PKAs at the layer just beneath the surface where the cluster atoms are not yet separated from each other. We can say that the primary and secondary knock-on processes saturate near the surface, and after that, the number of these knock-on atoms cannot increase any more with nuclearity.
Clearly, the increase of energy lengthens the ranges of knock-on atoms. In addition, the ranges become longer because the relation between the knock-on atom energy and its range changes with nuclearity. The following scaling law between the kinetic energy of a PKA E kin and its projected range R is fitted to the simulated data The projected range is the distance between the initial position and position at 500 fs of the knock-on atoms. At 500 ps, the PKAs have lost their initial energy in collisions, and after that they may move among the coherent flow of atoms. Thus, their positions at 500 fs gives a reasonable approximation of their ranges due to their initial kinetic energy. The parameters q and p are plotted as a function of N in figure 5. The exponential dependence on the energy reaches its constant value at N = 4. However, the parameter b increases with both energy and nuclearity. A coherent motion of knock-on atoms outwards from the cluster track can explain the result. The energy density in the saturated primary knock-on region is very high, which induces the outward pressure that lengthens the ranges. In the light of these results, it is clear that cluster collision cascades cannot be described neither with binary collision models nor with linearly overlapping cascades induced by individual cluster atoms. However, the lengthening of  (1)), that approximates the scaling between kinetic energy of PKAs and their projected ranges. Right: an example fit of equation (1) to PKA range values. For 1000 eV Ar 7 , p = 0.30 ± 0.01 and q = 0.38 ± 0.10 nm.
the primary and secondary ranges cannot alone explain the non-linear growth of the number of displaced atoms.
The number of PKAs, SKAs and TKAs per cluster atom decrease with N ( figure 6). The explanation is the saturation discussed earlier in this section. There is always a limited number of substrate atoms within the interaction distance from the cluster atoms. When the energy increases, the cluster penetrates deeper and the number of potential PKAs can increase slightly, as shown in the top left frame of figure 6. The same mechanisms affect the numbers of the SKAs and explain their decrease with N . In the tertiary collision process, the collision cascade expands radially and the number of available atoms does not limit the growth. The tertiary process continues until the foremost tertiary atoms do not have enough energy to knock more substrate atoms from their locations. Therefore, the number of TKAs is linearly proportional to N and E, as shown in the bottom right frame of figure 6.
The number of displaced atoms increases with N indicating that the volume of the destroyed region inside the collision cascade becomes larger (figure 6). The increase is clearly non-linear when N < 7, and after that it is almost linear. The scaling law for the number of displaced atoms is where r = 0.40 ± 0.02 and s = 0.23 ± 0.02 are parameters obtained by fitting the equation to the data shown in figure 6. At low nuclearities, only a small volume inside the knock-on region is destroyed, but at high nuclearities, the volume of the destroyed structure grows with the volume of the knock-on region and linearly with N . Because the sputtering yield grows non-linearly, the volume of the destroyed region inside the collision cascade alone cannot explain its non-linearity.
When N 4, the diameters and depths of the PKA, SKA and displaced atom distributions follow the scaling law  sputtering yield and the number of atoms ejected above the original surface plane. However, figure 7 shows that the linear scaling with N is probably not valid at higher nuclearities. The non-linear regime of these diameters and depths starts at around N = 7, whereas the area of the displaced atoms scales non-linearly with N even at lower values than N = 7. Thus, the coherent ejection of particles is triggered by some additional mechanism. Experimental velocity distributions of sputtered material and sizes of ejected clusters indicate that the flow of atoms is the dominant sputtering mechanism at high energies [40]- [42]. Figure 8 shows that the total energy deposited to the TKAs in the 3 nm surface layer forms a pattern where the maximum is located at a certain distance from the origin. Naturally, the maximum energy increases non-linearly with the distance from the origin, because the area of the tertiary knock-on circle increases as the distance squared. On the other hand, the vertical velocity maximum is located closer to the origin than the energy maximum, and its maximum is already reached at N = 7 (right frame of figure 8). When N > 7, there is an outer ring of TKAs, which expands laterally, and an inner ring, whose area does not increase with N . This separation of vertically moving hub and laterally moving ring emerges at the same N as the nonlinear sputtering (figure 1). The sputtering is proportional to the intersection area between the displacement cascade and the original surface plane. After the formation of vertically expanding ring of TKAs, the intersection area grows as N 2 , as shown earlier. Therefore, we conclude that the formation of the expanding outer ring explains the sputtering behaviour shown at high energies in figure 1. However, the sputtering is not directly induced by the collision process but  by the pressure of the Si gas in the displacement cascade. Therefore, the average velocity of the sputtered particles does not increase with N , because the cluster energy is evenly deposited in the knock-on process. This lateral expansion of the displacement cascade and the pressureinduced flow of atoms from the cascade to the vacuum becomes the most important sputtering mechanism at higher energies, as in fcc metals [8,35].
Ranges of displaced atoms increase non-linearly with nuclearity ( figure 9). At 1000 eV, the ranges are less than 1 nm, if N = 1-2, and less than 2 nm, if N = 1-7. This supports the conclusion that the velocities of displaced atoms increase with N . Together with the previous observation that the vertical velocities increase with N at the surface, we conclude that the nonlinear increase of sputtering yield is a consequence of the shift of the velocity distribution of displaced atoms towards high velocities.
The cluster and primary knock-on stopping occurs mostly during the first 150 fs, whereas the tertiary knock-on process takes several hundred femtoseconds, as shown in figure 11. Figure 10 visualizes the process where the knock-on atoms included at later times are mostly located at the boundaries of the collision cascade. The tertiary knock-on process takes over 1 ps, because it covers all other collisions except the primary and secondary collisions. Eventually, it ceases down to ordinary heat transport in the outer boundary of the collision cascade, when the colliding atoms no longer have enough energy to displace other atoms from their positions in the atomic structure. Figure 11 shows that the timing of the tertiary process does not change at all with energy and changes only slightly with nuclearity. At high nuclearities the maximum collision frequency occurs later than at lower nuclearities, but the difference is very small. Therefore, we can conclude that the increasing energy density inside the collision cascade really induces the non-linear sputtering. The flow to vacuum is induced before the energy diffuses to the inner parts of the substrate. Figure 12 describes the qualitative differences between collision of Ne, Ar and Xe clusters with the same kinetic energy. Because of the differences in masses of cluster atoms, the maximum possible scattering angle varies so that it is largest with Ne and smallest with Xe. Therefore, the trajectory bundle shapes are different. Ne atoms scatter at large angles forming an almost spherical bundle. On the other hand, the Xe track bundle is very conical in shape and the primary energy deposition occurs in a small volume compared to the Ne and Ar systems. The displacement cascade of the Ne impacts is dispersed, which leads to a smaller ejection yield. The Ar displacement cascade has largest volume, but the Xe cascade has the largest energy density, thus the sputtering yields are about the same at these energies. In all cases, the course of cascade development is qualitatively similar regardless of the noble gas used in the cluster beam.

Model
In order to verify that the non-linear behaviour mainly arises from the increasing energy density inside the destroyed region and from the timing of energy dissipation process, we have constructed a simple continuum model for collision cascade development, which we introduce and discuss in this section. It covers the phases from cluster impact on the surface to the start of the flow of atoms to vacuum. The model is not intended to be a precise theoretical description of the phenomenon. Its purpose is to demonstrate the emergence of non-linearity from some simple assumptions. It describes the geometrical shape of the energy distribution inside the substrate in three discrete phases. Many of the minor effects discussed in the previous section, like the clearing-the-way effect, are not included in the model. The model is based on the same physical assumptions as the models introduced for high-energy cluster impacts on fcc metals [8,35]. The energy is deposited from the cluster atoms in the primary collisions in a very limited region inside the substrate. This induces a high-energy displacement cascade, that cannot be considered as a linear combination of cascades of individual cluster atoms, but rather a single, combined cascade. In the third phase, the energy is deposited over a larger volume. This hemispherical liquid or gaseous volume releases its contents to vacuum leaving a crater on the surface. The assumptions are the same as in the MEDF model [43] discussed in the next section and developed independently of our model.
The model includes three discrete energy deposition phases, which is a rough approximation because the primary, secondary and tertiary processes overlap in time. In the first phase, the energy is deposited from the cluster atoms to the PKAs. The energy E 1 and momentum p 1 distributions depend on the cluster energy per atom E, nuclearity N , incident angle and ratio of masses of the cluster and substrate atoms m = M 1 /M 2 . The energy and momentum densities at a point rv in substrate are N , m, ), p 1 (r) = g 1 (E, N , m, ). Notice that, in this study, the tertiary process includes all other collision events except the primary and secondary collisions; thus TKAs are induced even after 1 ps. The middle frame shows a comparison of TKA induction processes when the nuclearity changes and the energy per atom is constant. In the right frame, the processes are compared when the energy varies.
In the second phase, the energy is deposited from the PKAs to SKAs, which produces the second distributions representing the initial positions, energies and momentums of the SKAs The third phase is a combination of all remaining collisions in the collision sequences. The result is the energy and momentum distributions just before the flow of atoms begins The number of atoms to be ejected above the original surface plane after the third phase is approximated by the volume that has higher energy density than a certain threshold value. In reality, the phases overlap in time ( figure 11).
The model is difficult to solve analytically, if the functions f i and g i are reasonably realistic. Therefore, we have solved the model calculating the distributions separately one after another in a three-dimensional rectangular grid.
In the first phase, the cluster energy is uniformly distributed in the region shown in the upper left frame of figure 14. Its boundary is defined by the following function: where a i are parameters, z is the depth from the original surface plane and x is the lateral distance. The parameters are chosen to provide the shape of the track bundles shown in figure 2, because the primary energy deposition occurs around the cluster atom tracks.
In the second phase, the range of the PKAs induced in the point r is approximated to be within the spherical zone between 2a 6 √ E 1 /3 and a 6 √ E 1 around the point. However, the momentum gained from the cluster atoms is on average perpendicular to the direction defined from the point of impact to the point r, thus the PKAs move only in the directions which are within the angle a 7 radians from that direction. This is a rough approximation, because the simulated asymmetric range distribution is more complicated (figure 4). The energy E 1 (r) is distributed uniformly into this part of the zone. The angle parameter a 7 can be used to fit the model to other projectile types which have different masses. It is also possible to model oblique incident impacts using asymmetric boundaries instead of equation (7). Integration over all knock-on atom source points gives the second energy distributions (E 2 ). As shown in the middle top frame in figure 14, the energy distribution is now moved downwards on average.  Notice that the simple distribution function employed here does not imply that the only energy deposition mechanism is the binary collision mechanism. The function only describes where the energy is deposited from each point.
In the third phase, the energy is again deposited in a spherical zone around the point of its origin. All scattering angles are now equally probable. The outer range is a 7 √ E 2 , where a 7 > a 6 because this deposition phase represents several sequential collisions. The parameters are again chosen to give distributions whose sizes agree with the distributions in the MD simulations ( figure 2). The resulting energy distribution E 3 is shown in the rightmost upper frame in figure 14. The frames in the lower frame show how the distribution grows with N . The fact that in the second and third phases the energy in each cell is distributed into the surrounding cells is just a statistical description for the various atomistic processes, and does not imply that binary collisions are the only mechanism.
The model is simulated using various N values. The volume where the energy density is above an arbitrary threshold values is plotted in figure 15 as a function of N . The shape of these curves are qualitatively similar to the curves describing the number of knock-on atoms, displaced atoms and ejected atoms in figures 1 and 6. Thus, the model gives qualitatively similar energy density scaling to the atomistic simulations. Because the energy density is shown to be the reason for non-linearities, the model demonstrates that the non-linearity emerges mostly because of high initial energy density and the timing of the energy dissipation process. The threshold energies plotted in figure 15 do not have any precise physical reason, but the energy of sputtered atoms is often less than 3 eV.
The model is not fitted to the simulated data to get quantitative agreement because the aim is only to demonstrate the effect of high initial energy density. However, the model in its present simple form gives almost correct diameters for the secondary and tertiary distributions, if the parameters of equation (7) are fitted to the simulated cluster track bundle shape. A more precise model is needed to make quantitative predictions, and it is questionable whether this kind of model is needed at all, because atomistic simulations can be used instead.

Discussion
The most important reason for the emergence of non-linearity is that the energy density in the inner parts of the displacement cascade increases with nuclearity. This happens because the cluster atoms stop in a relatively small volume and the size of the volume does not increase linearly with nuclearity ( figure 2). Instead of linear combination of collision cascades, a single, combined collision cascade is created. The consequent collision process dissipates the energy. The duration of the dissipation does not increase linearly with the nuclearity, which maintains the high energy density inside the displacement cascade and atoms start to escape to vacuum.
In addition to this simple picture, several other phenomena occur and affect the displacement cascade development. The number of PKAs decreases with the cluster nuclearity, but their energy increases, because they gain energy from several cluster atoms. Thus, their ranges become longer. The PKAs clear the way slightly increasing the range of cluster atoms, and this coherent motion also increases their own range. Similarly, the number of the SKAs decreases with the nuclearity and their average range increases slightly.
The timing of the energy dissipation processes is an important factor. The primary and secondary collision processes occur quite fast, and a distribution of high-energy atoms is formed in the substrate. The energy is deposited further in the tertiary knock-on process, which in this study includes all collisions other than the primary and secondary collisions. It lasts longer than the primary and secondary processes and the number of TKAs increases linearly with the nuclearity. Meanwhile, the original structure of the inner parts of the collision cascade is destroyed and the atoms start to flow to vacuum due to the pressure inside the collision cascade. The diameter of the destroyed structure grows faster than linearly with the nuclearity because of the increasing energy density. This gives rise to the non-linear sputtering.
It is clear that the interactions between atoms in the displacement cascade are so frequent that the process cannot be described as a binary collision cascade. It is also demonstrated that the displacement cascade develops as one entity and cannot be described as a linear combination of independently developing cascades induced by single cluster atoms.
Shao et al [44] have recently constructed a theoretical model for the overlap of collision cascades of individual cluster atoms. The model agrees with the experimental damage values measured in Si under carbon cluster bombardment. The results of our simulations support the model, but provide a more detailed atomic level description of the cascade development.
The Sigmund-Clausen model for thermal sputtering [45] applies to systems where material evaporates from the surface of the collision cascade. In this model, the cluster track is cylindrical and the sputtering increases more rapidly than E 2 d , where E d is the energy deposited per unit track length. The evaporation rate in the model is approximated by assuming that the collision cascade is a container of an ideal gas and the temperature in the intersection area with the substrate surface develops according to the equation of heat conduction. The Sigmund-Clausen model does not describe the total destruction of the substrate and the coherent flow of atoms, thus it cannot be directly applied in this study. In spite of this, the simulations and the Sigmund-Clausen model both give a similar result, that the material ejected above the original surface level and the atoms sputtered increase faster than linearly, when N > 7 and E > 500 eV. This is shown in figure 16, where the simulated sputtering yield and number of atoms ejected above the original surface plane are divided by the factor E 2 d = (E N ). At higher energies than used here (E > 10 keV atom −1 ), the largest flow of substrate atoms occurs at the boundaries of the crater, thus the evaporation approximation is not valid at these energies [35].
In our earlier study, we investigated the non-linear sputtering mechanism in Au(111) when it is bombarded with 20-280 keV atom −1 Au clusters (N = 5, 13) [8]. The conclusion was that the diameter of the displaced region is linearly proportional to N , and because the sputtering yield in that case is approximately proportional to the intersection area between the displaced region and the Au(111) surface, the yield is proportional to N 2 . This explains the Y ∝ N 2 relationship observed experimentally [46]. Because of channelling, the cluster atoms move in an Au target more independently than the Ar cluster atoms in amorphous Si. The knock-on atom stopping is also more complicated in Au because the atoms can channel in the lattice and the ranges can vary considerably. There are also differences in the mechanisms of how the target structure is destroyed. Therefore, the scaling laws for cascade and crater size, as well as for sputtering yields, vary depending on the target substrate. However, the common features are the stopping of cluster atoms in a relatively small region and the emergence of a single high-energy region and displacement cascade, which induces the non-linear scaling.
The simulations, our model and the other models of non-linear sputtering discussed here give a very coherent description of the emergence of non-linear sputtering. However, a fundamental question regarding all simulation studies is, are the results in agreement with reality. Because no direct experimental measurements of cascade development are available, the comparison with experimental results is indirect. The EDIP potential combined with the repulsive short-range potential used in the simulations reproduces experimental sputtering yields very well in the case of single Ar ion bombardment of Si(111) [24], thus it describes both the stopping phase and the cascade development phase reasonably well. Using larger cluster sizes in this study, it was found that the number of displaced atoms calculated in MD simulations generally are in good agreement with the experimental results [47]. Craters induced by Ar clusters on the crystalline Si surface in MD simulations have a similar form and size to the real craters [48]. However, Mazzarolo et al [49] have found that atom displacement energies in crystalline Si given by the commonly used Si interatomic potentials and a tight binding method are significantly different. Because the first phases of collision cascade development depend mostly on the repulsive Ar-Si and Si-Si potentials and the difference between the threshold displacement energies given by the EDIP potential is reasonably close to the values of the tight binding model, we conclude that the displacement cascade development would be very similar, if the potential gave exactly the same threshold energies as the tight binding method. All these results together make us confident about the overall reliability of the picture of collision cascade development presented in this paper.
In spite of the fact that the systems simulated in different research groups are not directly comparable, there are common findings that support the picture presented here. Aoki et al [13] have found in their MD simulations that the maximum penetration depth of Ar cluster atoms in Si(100) crystal does not depend on nuclearity and is proportional to the cube root of kinetic energy of the cluster atoms. The depth of SKA distribution calculated in this study depends on cluster atom penetration depth and ranges of the PKAs. This depth is also proportional to the cube root of the cluster atom energy and depends only weakly on nuclearity, as shown in table 1. Aoki et al [14] have also shown that when the total energy is constant, the number of displacements increases with increasing cluster size, but starts to decrease after the cluster size is large enough, for example larger than 10 000 at 50 keV cluster −1 . When the cluster size is constant, the number of displaced atoms per incident atom increases almost linearly with energy per atom. Medvedeva et al [7] have compared 1.5 keV atom −1 Au 2 and Al 4 impacts on Si(100) using MD. Al dimer atoms (27 amu) quickly disintegrate on the initial part of the trajectory while Au dimer atoms (197 amu) stay together, because the maximum scattering angle for an Au atom hitting a Si atom is only 8 • . Therefore, the energy deposited from an Au dimer is on average localized within a smaller region compared to the energy deposited from an Al dimer, and the cascades induced by the dimer atoms overlap producing a very compact cascade of energetic atoms. The enhancement factor for the sputtering yield is 7.0 for Au 2 but only 1.7 for Al 2 . The qualitative comparison shown in figure 12 is in agreement with these findings. It is found in simulations that the phosphorous dimer implant profile is almost indistinguishable from the atomic implant profile suggesting that the lattice-mediated vicinage effects are insignificant [15]. Our results are in agreement with that conclusion because the effects depending on N are found to occur at higher values than N = 2.
Andersen et al [50] have measured vertical ranges and range stragglings of small Au clusters in amorphous silicon. They found that the average range of cluster atoms is the same as the range of single ions at the same energy/atom. Also some broadening of the range distributions was found for the clusters. Three things should be considered when our simulations are compared with these results. Firstly, the mass of Au atoms is larger than Ar atoms, thus the trajectory bundle of Au atoms is even more compact and narrow than the Xe trajectory bundle in figure 12. As seen in figure 12, the range straggling decreases when the mass of cluster atoms increases. Secondly, the measured ranges are sums of the ranges in the stopping phase and the diffusion ranges during the cascade expansion phase. At 10.0-44.3 keV atom −1 , which is the energy in the experiments, the craters are deeper than the simulated craters in this work. Therefore, it is not probable that the cluster constituents will diffuse long distances towards the surface, because the cascade is expanding mainly laterally or 'downwards', and the atoms are ejected to vacuum mainly from the near surface parts of the cascade. Thirdly, at energies higher than 10 keV atom −1 , the cluster constituents travel in the substrate so fast that the small clearing-the-way effect reported in this work becomes negligible. Because of these three facts, we can conclude that the vertical ranges measured by Andersen et al [50] represent the stopping phase ranges of the cluster atoms and the clearing-the-way effect is probably very weak at the energies they used. On this basis, their main conclusion, that the ranges of the cluster atoms are the same as the ranges of single ions, is in agreement with our result shown in figure 3(a). This verifies that atoms of small clusters stop like single atoms, although they move near each other and often collide with the same substrate atoms. A similar behaviour is observed in simulations of Au cluster impacts on Au(111) [8].
Russo and Garrison [43] have developed the mesoscale energy deposition footprint (MEDF) model, which describes the relative yield and ejection volume size and shape differences between various cluster/substrate combinations. The key idea of the model is that the energy distribution formed in the substrate during the first phases of cluster impact (footprint) determines the position and shape of the displacements cascade and therefore also the sputtering yield. The model is in agreement with the findings reported here and our recent studies of cluster stopping in Au, which are summarized in the droplet model [8,35]. The MEDF and droplet models predict sputtering yields based on the mesoscale geometry of the cluster collision phenomenon. The geometry seems to be very universal among various substrates and cluster species. The cluster track bundles shown in figure 2 roughly correspond to the footprints in the MEDF model.
Finally, we make a remark about the validity of the scaling laws. The experimental value for the damage layer thickness in 60 eV atom −1 Ar 500 impacts is 20 nm [47]. If we use the scaling law (3) to predict the corresponding depth of the displacement region at the same energy and nuclearity, the result is 71 nm, which is clearly too much. This shows the scaling of the cascade dimensions with N is different for N 16. Nakayama et al [47] do not provide experimental values for the cluster sizes used in this study, thus direct comparison with our results is not possible. After simulating an ideal van der Waals solid, Anders et al concluded that the sputtering yield is linear in the incident energy above an energy threshold that depends slowly on cluster size. However, at lower energies, the yield is nonlinear. Our results support this conclusion, that both linear and nonlinear behaviours of sputtering yield and other quantities can be observed in the same impactor/substrate system depending on impactor energy and size.

Conclusions
Atomic level processes leading to the emergence of non-linear damage and sputtering when the cluster size increases are investigated using classical MD. The findings are summarized in a collision cascade development model, that is solved numerically in its most simple form. The results are compared to other theoretical sputtering models and results of simulations. The validity of simulations is discussed.
The results presented in this paper together with the results of other research groups give a very coherent picture of the rather complicated chain of events that lead to the emergence of a flow of particles from the substrate to vacuum after the impact of an energetic atomic cluster. Briefly, the non-linearities emerge due to the high-energy density induced in a relatively small region in the substrate in the cluster stopping phase, and because of the timing of consequent events that dissipate the energy over a larger volume of the substrate.
In the non-linear regime, the sputtering yield is often proportional to the cluster nuclearity squared N 2 . This dependence arises from the geometry of the non-linear cascade. The diameter of the hemispherical displacement cascade grows linearly with the number of atoms N in the cluster. Thus, the intersection of the cascade and the original surface plane grows proportional to N 2 . However, our current results show that the scaling of quantities with cluster atom energy and cluster nuclearity varies depending on the substrate and energy regime.