Direct observation of ion-induced self-organization and ripple propagation processes in atomistic simulations

ABSTRACT Patterns on sand generated by blowing winds are one of the most commonly seen phenomena driven by such a self-organization process, as has been observed at the nanoscale after ion irradiation. The origins of this effect have been under debate for decades. Now, a new methodology allows to simulate directly the ripple formation by high-fluence ion-irradiation. Since this approach does not pre-assume a mechanism to trigger self-organization, it can provide answers to the origin of the ripple formation mechanism. The surface atom displacement and a pile-up effect are the driving force of ripple formation, analogously to the macroscopic one. IMPACT STATEMENT The presented model allows to follow the ripple formation and propagation in different steps, at the atomic level, for the first time under low irradiation energies. GRAPHICAL ABSTRACT


Introduction
Self-organization effects are observed in many fields of science, and have the fascinating feature that similar effects can sometimes be observed in systems differing over many orders of magnitude in length and time scales. One such self-organization effect is that of ripple formation on surfaces. This effect is most commonly seen when caused by a moderate wind blowing on water or sand, or water flowing on sandy river bottoms [1]. Surprisingly, ripple formation has also been seen on the surfaces of solids at the nanoscale after prolonged energetic ion irradiation [2][3][4][5][6][7]. A number of ripple formation models based on the numerical solution of differential equations have been developed to explain the phenomenon. The theory proposed by Bradley and Harper [8], later extended by many other groups, explained the rippling effect as the modification of surface curvature due to sputtering under inclined incidence ( dominant mechanism, then the similarity between ionand wind-induced ripples is only coincidental). After the surface morphology has changed, the sputtering rate may also change. A few studies have taken into account this effect as a non-linear extension term related to the erosion rate [9][10][11][12][13][14]. However, it was recently shown the erosive component alone is not able to capture the complex phenomenon [15,16]. Norris et al. [17,18] proposed a mathematical model predicting patterns based on the crater function formalism, which allows also describing atom redistribution as an additional term contributing to ripple formation. Nevertheless, none of the analytical theories is able to capture the full process of selforganization or provide an intuitive understanding of its driving mechanisms. Molecular dynamics (MD) computer simulations of sequential ion impacts can, in principle, describe all the processes taking place concurrently under the ion irradiation condition, and give invaluable insights on the complexity of their interplay. MD simulations have indeed been used to provide important and clarifying insights on built-up stress [19,20] and its relation to ripple formation [21,22]. However, due to the very high computational costs of MD simulations of tens of thousands of keV ion impacts, direct formation and propagation of ripples have until now not been reproduced by this technique. Refs. [23,24] even directly state that MD simulations are not suited for studying the surface rippling due to time and size limitations. However, in recent experiments [15,[25][26][27], it was demonstrated that the 30 eV Ar ion irradiation caused the ripple formation on the surface of amorphous Si as well. In the current work, we gain new insights on the process of ripple formation due to a unique combination of a recently developed speedup scheme for MD simulation [28] to simulate tens of thousands of ion impacts on the same surface, with consideration of the new experimental condition of very low ion energy-induced ripples. The application of this new MD scheme enables us to observe directly the ripple formation and propagation by atomic self-organization at low energy irradiation, and to show that under negligible sputtering, the accumulation of displacement is the driving force of the ripple propagation similar to that observed experimentally.

Materials and methods
The simulations performed in this study were carried out using PARCAS MD code [29,30]. The relaxation process performed to obtain the a-Si structure used in this work is explained in detail in Appendix. More detail can be found in Refs. [15,31]. The size of the flat-surface cell used as a starting point 16.56 × 16.56 × 5.15 nm 3 ; it contained 73,584 atoms.
The consecutive 50,000 impacts by 30 eV Ar ions were effectively simulated at random points of the surface using a fixed azimuthal angle parallel to the x-axis. This was achieved by randomly shifting the cell in xand y-directions over the periodic boundaries prior to each impact. The impact point was always located in the center of the cell far from the border to reduce undesired high energy interactions over the periodic boundaries. The irradiation angle (θ) is 70 • off-normal. This temperature was restored to 300 K in the system after every nine impacts (1 ps long) using the Berendsen thermostat [32] during 30 ps.

Results
Results on the evolution are shown in snapshots in Figure 1, taken after 9000, 27,000 and 50,000 Ar + impacts, respectively. The continuous evolution is also illustrated in a movie, see Supplementary-Movie-1.mp4. Here the atoms are colored according to their zcoordinate (perpendicular to the surface). As we observe in this figure, the prolonged ion irradiation leads initially to random surface roughening. The second frame shows the emerging pattern that becomes fully clear in the last snapshot. The atoms on the ripple crests (the highest z-coordinates) clearly experienced the largest displacements, as shown in Figure 2(a). The further increase of the fluence promotes the propagation of the formed ripples.
Hence, a detailed analysis of surface evolution with the ion fluence may shed the light on atomic-level mechanisms driving the self-organization process.
First, we analyze the role of erosive mechanism in ripple formation. We observe a linear increase in the number of sputtered atoms with the fluence. We estimate the average sputtering yield S y = 0.03. However, when we compare this value to the one obtained in the previous study [15] for the single impact calculations (S y = 0.035 ± 0.013), we see that the increased curvature of the increasingly roughened surface does not affect significantly the sputtering yield and it remains within the error bars. We measure the role of erosion by analyzing the amount of missing atoms in the trough of the ripple, excluding the volume of its crest. The empty volume below the original surface can be found by using the OVITO surface mesh tool [33] similarly to Ref. [16]. The estimated difference in solid volume between the initial and the final configurations results in approximately 4500 atoms. The total number of sputtered atoms during the entire simulation did not exceed 1500. Assuming that the crest is generated by the re-deposited atoms, about 3000 atoms more would need to participate in the sputtering process, which may have ended in re-deposition. Hence, we conclude that the erosion can be safely excluded from the consideration as a prime mechanism of ripple formation at low incident energy. Note that at higher incident energies the sputtering is stronger and may play a more substantial role in pattern formation.
Next, we analyze the total displacement of atoms in a cascade that is found as δ = where w = {x, y, z} are the atom coordinates and N displaced is the number of atoms displaced in a cascade. In Figure 2(b), we draw the evolution of all three components of the total displacement, δ x , δ y (fluctuates around zero due to symmetry considerations) and δ z . In Figure 2(a) and it is noticed that most of the displacement accumulation takes place at the surface. δ z grows slowly in the positive z-direction (outwards from the surface). In the cascades, the atoms are displaced in all directions, however, the displacements towards the open surface are generally larger due to stress relaxation effects at the surface. Moreover, some atoms receive the momentum in the direction out, however, the transfered energy is not sufficient to overcome the surface barrier and sputter. Instead the atoms climb above the initial surface layer adding to the positive component of the δ z .
In Figure 2(b), we observe that the largest component of the total displacement is δ x (parallel to the ion-beam projection). The strong momentum transfer in the direction of the ion beam (positive x-direction) causes strong growth of δ + x with the fluence. Since δ − On the other hand, one might expect only small or no change of δ − x with fluence, since the negative displacements are significantly shorter and have little driving force for accumulation.
Initially, δ − x (Figure 2(b), the right y-axis) decreases linearly fast during the first few hundreds of ion impacts. This linear decrease reflects the accumulation of the total displacement in non-overlapping cascades (see Figure B2 in Appendix). The linear behavior continues until the shape of the surface changes. The ions which hit on the front, induce mainly stronger forward displacements. The ions impacting on the back slopes partially compensate for this effect, as they mainly scatter from the surface (grazing incidence), inducing backward atomic displacements. This is why δ − x still decreases but with a slightly reduced rate.
The picture changes when the cascades begin overlapping with previously affected regions. The total negative displacement does not decrease almost at all due to compensation by the δ + x from the following cascades (δ − x can even occasionally increase). The same process slows down the increase of δ + x , inducing non-linear growth of this component. This continues until about 27,000-30,000 impacts, when both surfaces became equally rough. Notice that, eventually, on average the crest-trough distance is about 1.5 nm, similar to the one measured in Ref. [15].
To verify the model, we performed the simulations for the 85 • incidence as well. The results showed the expected pattern of the ripples in the perpendicular mode (details to be published elsewhere).

Discussion
The analysis of the induced total displacement shows that indeed this parameter reflects the overall behavior of the atomic system.
To understand the different steps of the ripple formation, we focused only on the atoms that eventually form a ripple, following their history from the original positions on the initially flat surface. First, we saw that the motion of these atoms can be described as a biased random walk. In other words, we observe a diffusion-like process of atom migration on the surface. However, thermal diffusion at room temperature is well beyond the time span of MD simulations employed in this work, hence it cannot be used to explain the observed results. Nonetheless, the momentum received by the atoms in cascades stimulates migration of atoms more intensively on the surface than under it. Moreover, the strong momentum component in the x-direction induces a clear bias in this direction. This random motion of atoms creates initial random roughening on the surface, see Figure 3(a). In the following cascades, redistribution of the atoms within the forming random mounds triggers the displacement of the mound as a whole. In some cases, the formed mounds can be destroyed, however, the larger the rough feature grows, the less susceptible it becomes to destruction. While moving on the surface as a cluster, the mounds merge together forming a single ridge impeding the further motion of separate parts, see Figure 3  We also notice that the forward momentum transferred by the ion to the surface atoms depends on the effective incidence angle θ . The smaller the θ, the smaller the forward momentum and the shorter the displacement is induced by ion. However, on the back slope of a mound or a ridge, the ion hits at a grazing incidence which does not transfer to the surface atoms the momentum forward, but rather the momentum inwards stabilizing the back slope. This observation explains the pile-up effect and why the formed mounds and ridges survive and do not flatten down by the impacting ions (this process can be seen in Supplementary-Movie-1.mp4 and Supplementary-Movie-2.mp4). In Figure 4, we summarize the mechanism of ripple formation omitting for clarity the motion of separate mounds prior a ridge formation.
We note here that this process of the radiation-induced biased random walk motion does not require relaxation of internal stresses. However, we analyzed the stress buildup in our simulation cell and found no correlation with the ripple formation. The stress builds up under the surface in the direction of ion beam (σ xx ≈ 2 kbar), but it changes the sign in the direction perpendicular to the surface (σ zz ≈ −0.5 kbar) (see Appendix).
The Carter and Vishnyakov model [34], explained the ripple formation by atom relocation as well but fails to capture the propagation of the ripples, since it does not consider the movement of whole mounds initially formed by piling-up of displaced surface atoms. Our simulation follows all stages of the process: from atom relocation to piling-up and merging of mounds into larger structures eventually forming a continuous ripple, which continues propagating, coarsening and slowing down as well.
As a final point of discussion, we note that this explanation of ripple formations shows that, despite the almost ten orders of magnitude different length scale, there is indeed an analogy of the ion-induced ripple formation mechanism with that of macroscopic ripples on water and sand. Although there are many mechanisms that may cause the macroscopic ripples, the base factors causing these include redeposition of material due to an external driving force [1], in clear analogy with the diffusion-like flow and pile-up effects shown here to occur for an ioninduced driving force. The pile-up effect at the first stages of the simulation can be seen as a product of saltation and reptation processes [35], by analogy of the atoms with sand grains, for instance. However, from this point, we observe a merge of these initial mounds in larger structures due to a reptation-like biased movement, explained in Figure 4.
In summary, using MD simulations we observed directly the pattern formation at low energy (30 eV) irradiation by Ar + ions on a-Si. The observed selforganization was mainly triggered by the radiationdriven random walk biased in the direction of the projection of the ion beam (x-direction) at 70 • off-normal incidence. In the present simulations with low incident energy, we observe a strong effect of atomic redistribution due to collision cascades on ripple formation, but weaker dependence on sputtering or stress buildup. The slowing down in the growth dynamics is explained by increasingly overlapping cascades, where the previous forward displacements are compensated by the backward displacements of the newly coming cascades. Detailed analysis of the origin of the atoms comprising the final ridge revealed that self-organization is fully explained by pile-up and irradiation-enhanced migration of surface atoms.

Disclosure statement
No potential conflict of interest was reported by the authors.