MD-Predicted Phase diagrams for Pattern Formation due to Ion Irradiation

MD-Predicted Phase diagrams for Pattern Formation due to Ion Irradiation



Energetic particle irradiation of solids can cause surface ultra-smoothening [1], self-organized nanoscale pattern formation [2], or degradation of the structural integrity of nuclear reactor components [3]. Periodic patterns including high-aspect ratio quantum dots [4], with occasional long-range order [5] and characteristic spacing as small as 7 nm [6], have stimulated interest in this method as a means of sub-lithographic nanofabrication [7]. Despite intensive research there is little fundamental understanding of the mechanisms governing the selection of smooth or patterned surfaces, and precisely which physical effects cause observed transitions between different regimes [8, 9] has remained a matter of speculation [10]. Here we report the first prediction of the mechanism governing the transition from corrugated surfaces to flatness, using only parameter-free molecular dynamics simulations of single-ion impact induced crater formation as input into a multi-scale analysis, and showing good agreement with experiment. Our results overturn the paradigm attributing these phenomena to the removal of target atoms via sputter erosion. Instead, the mechanism dominating both stability and instability is shown to be the impact-induced redistribution of target atoms that are not sputtered away, with erosive effects being essentially irrelevant. The predictions are relevant in the context of tungsten plasma-facing fusion reactor walls which, despite a sputter erosion rate that is essentially zero, develop, under some conditions, a mysterious nanoscale topography leading to surface degradation. Our results suggest that degradation processes originating in impact-induced target atom redistribution effects may be important, and hence that an extremely low sputter erosion rate is an insufficient design criterion for morphologically stable solid surfaces under energetic particle irradiation.


Southern Methodist University, Dallas, TX 75205

2 Department of Physics and Helsinki Institute of Physics, P.O. Box 43, FIN-0014 University of Helsinki, Helsinki, Finland

Harvard School of Engineering and Applied Sciences, Cambridge, MA 02138

At irradiation energies between , the irradiation process is dominated by the nuclear collision cascade caused by ion impact [11, 12]. Displaced atoms that reach the surface with enough kinetic energy to leave are permanently sputtered away; all other displaced atoms come to rest within the solid or on the surface after phonon emission times of seconds. These processes contribute prompt erosive [13, 14] and prompt redistributive [15, 16, 10] components of morphology evolution, respectively, and are collectively denoted . For most materials other than elemental metals, the damage resulting from these collisions quickly ( seconds) leads to the amorphization of a thin layer of target material. Over much longer time scales ( seconds), mass transport by kinetic relaxation processes causes a gradual relaxational effect [13, 17].2 Hence, to the prompt term we add a phenomenological term for the gradual relaxation regime denoted , assuming a mechanism of ion-enhanced viscous flow, which is expected to predominate in irradiated amorphous materials near room temperature [17]. The prompt and gradual contributions to the rate of motion of the surface in the normal direction superpose:


The prompt regime may be characterized using molecular dynamics (MD) simulations [16, 19] or experimental methods [20]. Given data from many impact events, we may obtain the “crater function” describing the average change in local surface height at a point resulting from a single-ion impact at , with incidence angle .3 We then upscale the crater function into a continuum partial differential equation for the surface evolution using a multi-scale framework. The theoretical formalism for this process is described elsewhere [21]; here we provide a brief summary of the important points for the linear case. Given the crater function and the flux distribution , we write the prompt contribution to surface evolution as a flux-weighted integral of the crater function [22, 10]:


A well-known observation in the field is that scale of the craters is much smaller than the scale of the resulting pattern and of the flux distribution. To exploit this fact, we employ a formal multiple-scale analysis, based on a small parameter related to the ratio of impact scales to pattern scales. This formalism allows ready separation of the spatial dependence of the crater function (fast) from that of the surface shape (slow), and leads eventually to an upscaled description of the prompt regime of the form


Here the represent surface divergences, and the (increasing-order) tensors are simply the angle-dependent moments of the crater function . This compact formulation is interesting for two reasons. First, the moments are readily obtainable directly via MD simulation, and they converge with far fewer trials than do descriptions of the entire crater function. Second, while atomistic methods have been used in the past to obtain the amplitude of a single term in a PDE obtained via phenomenological modeling [23, 16, 24, 25], we believe this is the first derivation of an entire PDE from molecular dynamics results.

Equation (3) describes the prompt regime . To fully capture the surface dynamics, we add to this a relaxation mechanism associated with ion-enhanced viscous flow [17]. Together, and completely determine surface morphology evolution via Equation (1). From this (nonlinear) equation, pattern-forming predictions are then obtained by examining stability of the the linearized equation as a function of the laboratory incidence angle . The derivation follows that in [13], and in an appropriate frame of reference one finds that the magnitude of infinitesimal perturbations away from a flat surface evolve, to leading order in , according to the PDE


where the angle-dependent coefficients


are determined from the first moments obtained via MD, and the constant coefficient is estimated from independent experiments. The structure of Equation (4) indicates that linear stability is determined strictly by the signs of the calculated coefficients : for values of where either of these coefficients is negative, linearly unstable modes exist and we expect patterns, whereas for values of where they are both positive we expect flat, stable surfaces.

Existing uses of MD crater data for investigations of surface pattern-forming are entirely numerical in nature [26], and could be viewed as a scheme for numerically integrating Equation 2. In contrast, our analytical upscaling of Equation 2 illustrates exactly which qualities of the crater – namely, its moments – play the dominant role in surface evolution. Furthermore, this analytical form can be linearized, allowing predictions of stability boundaries, and changes to those boundaries as crater shape is varied. A crucial component of our approach is that the crater function – and hence the moments – contains the contributions of both erosion and mass redistribution. Whereas these effects have traditionally been treated separately via unrelated phenomenological models, viewing the crater function as fundamental integrates erosion and redistribution into a unified description, allowing both processes to be treated identically and readily separated and compared. Indeed, this approach has permitted us to confirm for the first time conjectures [22, 10, 19] that the stability of irradiated surfaces could be dominated by redistributive effects, with erosion – long assumed to be the source of roughening – being essentially irrelevant.

We obtain angle-dependent moments from a series of MD simulations, in an environment consisting of an amorphous, 20x20x10 nm Si target consisting of 219,488 atoms created using the Wooten/Winer/Weaire (WWW) method [27], and then annealed with the EDIP potential [28]. This gives an optimized amorphous structure where most of the Si atoms have coordination number 4. The target was then bombarded with Ar at 100 eV and 250 eV. During bombardment, the interaction between Si atoms was again described using the EDIP potential, which gives a good agreement between simulated and experimental sputtering yields [29], whereas the Ar-Si interaction was a potential calculated for the Ar-Si dimer [30]. The kinetic energy was gradually removed during the simulations from the 1 nm borders of the substrate to prevent it re-entering the impact area via the periodic boundary conditions used in the simulations. The ambient temperature in the simulations was 0 K. The simulation arrangements and their suitability for cluster and ion bombardment simulations are discussed in more detail in the supplement and in Refs. [31, 32, 33, 29].

For each energy, two hundred impacts were simulated at each incidence angle in 5-degree increments between 0 and 90 degrees, yielding moments as summarized in Figure 1. From the initial and final atomic positions, moments were obtained using the method described in [21]. Briefly, by assuming that densities in the amorphous layer attain a steady state (i.e., that defect distributions immediately project to the target surface), we obtain erosive moments by assigning a height loss at each location proportional to the number of sputtered atoms originating from that location, while redistributive moments were obtained by assigning height losses at initial atomic positions and height gains at final atomic positions. For use within our analytical framework, we also fit the moments to Fourier series constrained by symmetry conditions and by the observation that all moments tend to zero at . For both energies, the redistributive first moments are much larger in magnitude than – and have the opposite sign of – the erosive moments. The implication is that redistributive effects completely dominate erosive effects, except possibly at the highest (grazing) angles where all moments tend to zero.

To corroborate this finding, we calculate in Fig. 2 the coefficients in equation (5) for the 250 eV moments, and compare the pattern wavelengths they predict to experimental observations in the same environmental conditions [34] (clean linear experimental data at 100 eV are not currently available). The agreement is generally good at intermediate angles, but several discrepancies between theory and experiment should be addressed. First, the small quantitative difference between predicted and observed bifurcation angles – which depends only on the shape of – could readily arise from the approximate nature of the classical potential, on which our simulations are based (unlike the bifurcation angle, precise wavelength values depend on the value of , which could only be estimated). Second, the measured moments do not predict a transition to perpendicular modes at the highest angles; this could be due to our neglect of explicit curvature-dependence in the crater function, but additional physical effects such as shadowing and surface channeling – not addressed here – are known to be important at grazing angles. Third, we find no prediction via MD of the experimentally-observed perpendicular-mode ripples at low angles[9]; indeed, our redistribution-dominated continuum PDE is maximally stable at low angles, and equations of the general form 4 are anyway unable to generate the constant-wavelength or “Type I” bifurcation that is observed [10]. These observations, together with the experimental observation that low-angle ripples develop over much longer timescales than their high-angle counterparts [34], suggest that low-angle perpendicular-mode ripples are not due to crater-function effects at all. It has already been observed [10] that the low-angle Type I bifurcation is consistent with any of several non-local physical mechanisms such as gradual stress buildup and relaxation [10], or non-local damping [35]. None of these effects would be captured in the (prompt, local) crater function, and this observation motivates future studies incorporating such physics.

Despite the limitations of our approach, when one considers the lack of any free parameters in the theory, the agreement for the diverging-wavelength or “Type II” bifurcation [36] near 50 degrees is remarkably good. The agreement remains good even when the erosive coefficients are omitted, and the similar shapes of the redistributive moments at 100 and 250 eV is consistent with the reported [34] energy-insensitivity of the stability boundary. The most striking aspect of this result is its logical conclusion that erosion effects are essentially irrelevant for determining the patterns: according to Fig. 2 the contributions of redistributive effects to the coefficients, which determine stability and patterns, are about an order of magnitude greater and opposite in sign. This conclusion overturns the erosion-based paradigm that has dominated the field for two decades [13] and we suggest its replacement with a redistribution-based paradigm. An important direction for future research is identifying the energy range over which this conclusion holds. Although it is conceivable that erosive effects might become non-negligible at higher ion energies, or for dense metals where heat spikes enhance sputtering yields, it remains to be determined whether erosion is actually important for stability or pattern formation in any physical experiment to date. Preliminary analysis of simulations at 1 keV are consistent with our conclusions from the results at 100 eV and 250 eV, and experimental results at 1 keV lead to the same conclusion [37].

Non-erosive ion impact-induced atom redistribution at surfaces has heretofore not been considered in the design of plasma-facing fusion reactor wall materials, where low sputter yield has been an important design criterion in the selection of tungsten for stable surfaces that must be exposed to large plasma particle fluxes for extended periods. Because the average helium ion energy is only and the threshold energy for sputter removal of tungsten is , this material has been considered impervious to the effects of erosion. Within the erosion-based paradigm of pattern formation, the nanoscopic surface morphology that evolves on tungsten surfaces under some such conditions has therefore appeared mysterious [3]. However, because the negligibility of erosive effects does not prevent redistributive effects from causing pattern-forming instabilities, as we have shown here through a crater-function analysis, atom redistributive effects may be important contributors to the origin of these mysterious morphologies, and we present in the Supplement a re-analysis of data gathered in [38] that supports this idea. If this conjecture turns out to be correct, then ultimately crater function engineering considerations may provide a more refined materials design criterion than simply a low sputter yield.


S.A.N. and M.P.B. were supported by the National Science Foundation through the Division of Mathematical Sciences, M.P.B. was additionally supported through the Harvard MRSEC and the Kavli Insitute for Bionano Science and Technology at Harvard University. M.J.A. and C.S.M. were supported by Department of Energy grant DE-FG02-06 ER46335. The work of J.S., L.B., M.B., D.F., and K.N. was performed within the Finnish Centre of Excellence in Computational Molecular Science (CMS), financed by The Academy of Finland and the University of Helsinki; grants of computer time from the Center for Scientific Computing in Espoo, Finland, are gratefully acknowledged. We also thank M.J. Baldwin, N. Kalyanasundaram, and H.T. Johnson for helpful discussions.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to M.J.A. (email:

Figure 1: Fitted angle-dependent moments of the crater function as determined from molecular dynamics, for both 100 eV and 250 eV. (a) Zeroeth erosive moment (sputter yield times atomic volume). (b,c) First erosive (b) and redistributive (c) moments . Each of the latter contain components in both the (downbeam) and (crossbeam) directions, with the latter expected to be zero from symmetry arguments [21]. At both energies, redistribution dominates erosion.
Figure 2: Comparison between predicted and measured wavelength for 250 eV Ar Si. Dimensional coefficients were calculated with a flux of for comparison with experiment [34]. (a,b) Coefficients of and in the linearized evolution equation 4 using the experimental flux. These coefficients are dominated by redistributive effects. (c) Comparison of predicted ripple wavelengths with average experimentally-observed wavelengths. Circles / squares indicate experimental patterns with wavevector parallel / perpendicular to the ion beam, and the vertical dashed black lines indicate experimental phase boundaries. On top of this, the solid line indicates our predicted wavelengths (which are all parallel-mode), and the dashed blue line indicates the wavelengths predicted if erosion were neglected entirely.

Appendix A Simulation Methods

In the simulations, the actual impact induced structural deformations can be detected only in a well relaxed silicon structure. A stress existing within the material as a result of incomplete relaxation of the structure before the simulation can induce displacements of silicon atoms upon impact. For example, our test simulations showed that a structure created by annealing silicon in MD is not dense enough to model the real amorphous silicon. The surface of such a structure collapses upon impact. It is possible to relax the internal stresses by bombarding the surface with silicon atoms before the actual simulation. However, the relaxation will be not uniform. Therefore, we have used the Wooten/Winer/Weaire (WWW) method [27].

The WWW method is a computer algorithm that generates realistic random network models of amorphous silicon. In this method, the structure is described with positions of N atoms and a list of bonds between the atoms. After a random switch of two bonds, the structure is relaxed using an interatomic potential [39, 40]. In connection to the structural relaxation, the Berendsen pressure control algorithm is used to relax the diagonal components of the pressure tensor to zero [41]. The algorithm is computationally demanding and therefore it is possible to fully optimize only relatively small structures of a few thousand atoms. Therefore, the optimized block must be copied and the copies joined to create an optimized silicon structure large enough for impact simulations. Our tests showed that the best approach is to partially optimize a rather large (10x10x10 nm) building block instead of building the structure of very small fully optimized blocks. The latter approach can induce artificial internal shear stresses in spite of the optimization (note that pressure control at periodic boundaries an an orthogonal cell does not necessarily remove the nondiagonal, shear components of the stress tensor). The optimization of a 10x10x10 nm amorphous silicon structure was achieved using a parallelized implementation of the WWW algorithm.

After the WWW optimization phase, the structure and the surface were relaxed in MD using the EDIP potential, before it was used in the actual impact simulations. The structure used in the simulations was built of four identical optimized blocks. The density was , which indicates that the structure is dense and not likely to collapse upon impact. Two-thirds of the silicon atoms had four neighbours, the others had five neighbours, which is also a sign of a dense structure. In a perfect amorphous network, all silicon atoms should have four neighbors. However, the test with structures built of smaller blocks that were better optimized than the 10x10x10 nm structure showed that the behavior of the moments as a function of the impact angle are very similar as those reported in the main article.

To simulate ion irradiation of bulk silicon samples, we have used periodic boundary conditions with slab boundary conditions (free boundary in the incoming ion direction, and periodic boundaries in the other two lateral directions). The bottom 1 nm layer of atoms in the simulation cell were held fixed, and temperature scaling was also applied to the atoms in a 1 nm thick layer above it. A border region of 1 nm thickness in the lateral directions was cooled during the simulations. The role of the cooling zones is to prevent shock waves and phonons to re-enter the impact area via the periodic boundaries. The size of the simulation box (20x20x10 nm) was chosen to be large enough to contain not only the area containing the crater but also the area of small deformations which reach 5-7 nm from the impact point.

In addition to the good amorphous silicon structure and cooling of the boundaries, the quality of the repulsive Ar-Si potential affects the outcome of impact simulations. The kinetic energy of the impacting Ar atom is deposited first in relatively strong Ar-Si collisions where the atoms come close each other. The repulsive Ar-Si potential used in this study is calculated using density-functional-theory calculations utilizing a numerical basis sets. This method is shown [30] to give a more accurate repulsive potential than the standard universal ZBL potential [42]. The same method was used to create Si-Si short-range repulsive potential which was smoothly joined to the EDIP potential. This ensures that the possible collisions between high-energy primary knock-on silicon atoms are correctly modelled.

The same initial structure was used for all simulations. However, to reduce as much as possible the effect of the initial surface structure on the measured crater statistics, the following steps were taken. First, the impact point was varied randomly over the entire surface, with (periodic) cooling regions dynamically re-assigned for each simulation so as to be maximally distant from the impact point. Second, the azimuthal angle was varied randomly, and only the atoms within 7.5 nm of the impact point were included in the moment calculations (i.e., cooling zones and corners of the square box were ignored). Last, the global average – or “background” – displacement (over all simulations at a given impact angle) was subtracted from each individual displacement before analysis. Because we varied the impact point and azimuthal angle randomly, this global average ought to be zero for a perfect target over enough trials, and subtracting it helps to remove target-specific effects from our measurements. With this setup, 200 simulations at each impact angle were performed, which was sufficient to isolate the background displacement, and to cancel the effect of thermal vibrations which were present in the structure after impact.

A final challenge in analyzing the data arose due to a combination of the amorphous nature of the target with the periodic boundary conditions. On an ideal, very large MD target, the effect of the ion impact would only be felt within a finite distance from the impact point. However, to allow the gathering of data within reasonable time, a limited box size must be used, and the periodic box described above – with cooling zones – appears to be the most physically plausible way of accomplishing this. However, a periodic box always means that one is truly simulating an infinite number of simultaneous side-by-side impacts, and with a small enough domain, these impacts can generate enough co-ordinated momentum transfer to shear the entire target, especially for an amorphous target. Indeed, within our target, the average downbeam displacement of atoms was consistently a linear function of distance from the target’s rigid floor, except near the surface, where larger displacements were concentrated.

To combat the problem of ”global shear,” we measured the average downbeam displacement within different annular slices of target parallel to the surface, using inner/outer radii of 2nm/9nm (i.e., away from both the impact and the cooling boundaries). We then fit the bottom half of the resulting depth-dependent profile to a line using least squares, and finally subtracted the extrapolation to the surface of this fit from the overall displacement field, as illustrated in Figure 3. The results localize the displacements to within a few nanometers of the surface, which is consistent with measured amorphous layer thicknesses of 3 nm for Si irradiated by Ar at 250 eV [34]. In the future, we will explore the response of a larger, hybrid target consisting of a 3 nm layer of amorphous Si atop a crystalline (but not rigid) base. However, we believe our existing measurements are within the accuracy level of the other estimates in the paper.

Appendix B Moment Capture and Fitting

Here we describe in more detail how we obtain moments, how we fit them, and how we obtain the final linearized evolution equation (4). For each simulated ion impact with our initially flat MD target, we define a co-ordinate system centered at the impact point with pointing normal to the surface, pointed in the direction of the projected ion path, and perpendicular to both and so as to complete a right-handed co-ordinate system. Hence, the ion is always arriving from the negative -direction. Hence, the crater function describes the height change associated with an impact at the origin, of an ion with indicence angle of from the vertical, coming from the negative -direction.

After impact, we extract the moments


using the method described in [21]. Here, it is important to note that the first moments contain both erosive and redistributive components (because the zeroeth moment describes mass loss, and because redistribution is mass-conserving, has no redistributive component).

Simulation sets were performed at 5-degree increments, and the average of the resulting moments were fit to Fourier functions of the form


These fittings reflect the observation that all moments tend to zero at , and also their symmetries about (because a positive theta indicates an ion beam coming from the negative -direction, moments that are odd/even in should also be odd/even in ). As seen in the main text, this method does not produce perfect fits, but the use of simple Fourier modes eliminates potential model-bias, while the restriction to a small number of terms excludes inter-angle noise from the fitted curve.

For our data, the moment fits were given by:


The moments are zero to within sampling error, as expected from symmetry considerations.

Appendix C Analysis: from Moments to Coefficients

For the general reduction of moments to (nonlinear) PDE terms , we refer to the framework derived in [21]. However, in the linear case discussed here, it is sufficient to consider the linearization of the first-order term obtained by combining Equations (1) and (3) from the main text:


(in the linearization, the zeroeth- and second-moment terms do not contribute to stability). Now, indicates a surface divergence, and indeed this calculation is most naturally performed in a local co-ordinate system associated with the surface normal and projected beam direction. In particular, both the flux and the moments depend on the local incidence angle , while the vector is observed to always point in the direction of the projected ion beam.

Following [21], surface velocities at an arbitrary point will be calculated in a local co-ordinate system centered at , where corresponds to the surface normal, corresponds to the downbeam direction associated with the projected ion beam, and . In this system the surface is described locally by the equation , the ion flux is , and the first moment is , where as measured from MD. Now, as described in [21], it is sufficient for the purposes of calculating one surface divergence to approximate via


where , , and describe the surface curvature at . All other variable quantities can then be approximated in the vicinity of to first order in via:


When we now take the surface divergence and evaluate at (i.e., at ), we obtain in the local co-ordinate system




While linear in the local co-ordinates, Equation (12) is in general nonlinear in the lab fram. However, for stability studies we need only the linearization of (12) in the lab frame, which in dimensional form is simply


because, to linear order,


To expression (14) for the prompt regime we add the linearization of the gradual regime associated with ion enhanced viscous flow [17], which is a lubrication approximation with the form


Adding the prompt and gradual regimes, we obtain the evolution equation (4) in the main text.

For completeness, we conclude with the functional forms of the coefficients associatd with our fittings, which are obtained from equations (13):


Appendix D Estimation of Viscous Flow coefficient

The materials parameter appearing in Equation (4) of the main text is defined [43] as


where is the surface free energy, is the thickness of the thin amorphous layer that is being stimulated by the ion irradiation, and is the layer’s viscosity. We assume the surface free energy of amorphous silicon under ion irradiation to be equal to its value in the absence of irradiation; the value of measured via molecular dynamics simulations by Vauth and Mayr [44] happens to be numerically equal to that measured experimentally for single-crystal Si(001) [45]. For the amorphous layer thickness, we directly measured via cross-sectional transmission electron microscopy on samples irradiated at normal incidence and 30 degrees from normal. Finally, we estimate the viscosity of the top amorphous Si layer during irradiation to be , as shown below.

The reciprocal of viscosity is the fluidity , which is generally understood to scale with the flux , and can be expressed in the form


where is the radiation-induced fluidity, and is the average number of displacements per atom per second. Using molecular dynamics simulations, Vauth and Mayr [44] report at an energy and temperature – we use this value for our comparison with experiment, with the caveats discussed below. The average number of displacements per atom per second is given by


where is the atomic volume of silicon, is the experimental flux in the plane perpendicular to the ion beam, is the amorphous layer thickness, and is the number of recoils generated per ion impact. To estimate , we use the Kinchin-Pease [46] model for the gross number of Frenkel pairs per incident ion, obtaining , where is the ion beam energy and is the displacement threshold energy of Si [47]. Taking all of these quantities, we obtain a value of the viscosity of


As discussed above, the value for listed in (21) is associated with Vauth and Mayr’s value of at an energy and temperature . In contrast, our experiments were carried out at , and the irradiated sample is observed to reach temperatures of approximately . Hence, there is some uncertainty in our value of , which translates to uncertainty regarding the vertical position of the theoretical curve in Figure 2 of the main text. For the temperature difference, Vauth and Mayr observe to increase weakly with increasing substrate temperature, which would shift the curve upward; however, the shift would likely be less than a factor of two. As for the energy difference, in a study of CuTi, another amorphous material, Mayr et al [48] found to either increase or decrease with recoil energy depending on the details of the simulations; hence the theoretical curve in Fig. 2 of the main text could shift either upward or downward for MD simulations at 250 eV, with a potential magnitude of perhaps a factor of two.

Appendix E Preliminary Supportive Data

We believe the phenomena reported in this paper are applicable to a wide variety of systems. However, the main data sets leading to the discovery of this result both describe low-energy irradiation of amorphous silicon. Therefore, we provide here some preliminary results associated with higher-energy argon irradiation of silicon, and also for the case of helium-irradiated tungsten which was referred to in our conclusion.

For silicon irradiated by argon at 1 keV, a much larger target must be used, with dimensions of 40 x 40 x 10 nm, which makes simulations expensive. However, we have performed 30 simulations at 5-degree increments between 30 and 65 degrees, which spans the experimental transition, and the results so far are consistent with those at 100 and 250 eV: the redistributive contributions to the first moment dominate the erosive contributions.

For tungsten irradiated by helium at 100 eV, average displacements are so small that we were not able to isolate a clear signal against background noise after 200 simulations. However, using existing data from earlier work [38], we were able to calculate the cumulative displacement field at a single angle after 4,000 simulations. In all of these simulations, not a single tungsten atom was sputtered, yet a downbeam bias in the displacement field is clearly observed in Figure 3. Hence a crater-function approach may enable better engineering of surface stability under conditions where the sputter yield is zero but impact-induced target atom displacements still occur.

Figure 3: Illustration of the removal of shear at 60 degrees for irradiation at 250 eV. Green dots are the original measurements, and the blue line represents a linear fit to the bottom half of the dots, extrapolated to the surface. The red line is the result of subtracting this extrapolation from the original data, which localizes the displacements to the vicinity of the surface. All data are associated with a 2nm/9nm annulus that masks the impact zone and the cooling boundary zone.
Figure 4: Cumulative displacement field for 100 eV He W after 4,000 impacts. Impinging helium ion comes from upper right, at an angle of 25 degrees from normal. The sticks combine the initial and final position of atoms that were initially in a half a unit cell thick cross section through the simulation cell. The displacements are analyzed for data initially simulated in Ref. [38]. The blue ends of the sticks indicate the initial and the red ends the final positions of the atoms. The surface is located at the top. No W atoms were sputtered, but a clear downbeam bias in the displacement field is seen.


  1. footnotetext: Submitted Sept. 15, 2010
  2. Note that our focus on amorphous materials precludes potentially confounding effects due to singular crystallographic energetics and kinetics, such as the Villain instability [18], in which surface diffusion is destabilizing and may be responsible for pattern formation on faceted single-crystal metal surfaces at certain temperatures [2].
  3. In principle, also contains an explicit dependence on the initial surface shape (as opposed to the implicit dependence that arises via the angle-dependence). This effect is described in [21], but is very difficult to capture with molecular dynamics, so we leave its analysis for future work.


  1. Yamada, I., Matsuo, J., Toyoda, N. & Kirkpatrick, A. Materials processing by gas cluster ion beams. Materials Science & Engineering R-Reports 34, 231 (2001).
  2. Chan, W. & Chason, E. Making waves: kinetic processes controlling surface evolution during low energy ion sputtering. Journal of Applied Physics 101, 121301 (2007).
  3. Baldwin, M. & Doerner, R. Helium induced nanoscopic morphology on tungsten under fusion relevant plasma conditions. Nuclear Fusion 48, 035001 (2008).
  4. Facsko, S. et al. Formation of ordered nanoscale semiconductor dots by ion sputtering. Science 285, 1551–1553 (1999).
  5. Ziberi, B., Frost, F., Höche, T. & Rauschenbach, B. Ripple pattern formation on silicon surfaces by low-energy ion-beam erosion: Experiment and theory. Physical Review B 72, 235310 (2005).
  6. Wei, Q. et al. Ordered nanocrystals on argon ion sputtered polymer film. Chemical Physics Letters 452, 124 (2008).
  7. Cuenat, A. et al. Lateral templating for guided self-organization of sputter morphologies. Advanced Materials 17, 2845 (2005).
  8. Ziberi, B., Frost, F., Tartz, M., Neumann, H. & Rauschenbach, B. Ripple rotation, pattern transitions, and long range ordered dots on silicon by ion beam erosion. Applied Physics Letters 92, 063102 (2008).
  9. Madi, C. et al. Multiple bifurcation types and the linear dynamics of ion sputtered surfaces. Physical Review Letters 101, 246102 (2008).
  10. Davidovitch, B., Aziz, M. & Brenner, M. On the stabilization of ion sputtered surfaces. Physical Review B 76, 205420 (2007).
  11. Sigmund, P. Theory of sputtering. I. Sputtering yield of amorphous and polycrystalline targets. Physical Review 184, 383–416 (1969).
  12. Sigmund, P. A mechanism of surface micro-roughening by ion bombardment. Journal of Materials Science 8, 1545–1553 (1973).
  13. Bradley, R. & Harper, J. Theory of ripple topography induced by ion bombardment. Journal of Vacuum Science and Technology 6, 2390–2395 (1988).
  14. Makeev, M., Cuerno, R. & Barabási, A.-L. Morphology of ion-sputtered surfaces. Nuclear Instruments and Methods in Physics Research B 197, 185–227 (2002).
  15. Carter, G. & Vishnyakov, V. Roughening and ripple instabilities on ion-bombarded si. Physical Review B 54, 17647–17653 (1996).
  16. Moseler, M., Gumbsch, P., Casiraghi, C., Ferrari, A. & Robertson, J. The ultrasmoothness of diamond-like carbon surfaces. Science 309, 1545–1548 (2005).
  17. Umbach, C., Headrick, R. & Chang, K.-C. Spontaneous nanoscale corrugation of ion-eroded : The role of ion-irradiation-enhanced viscous flow. Physical Reveiw Letters 87, 246104 (2001).
  18. Villain, J. Continuum models of crystal growth from atomic beams with and without desorption. Journal de Physique I (France) 1, 19–42 (1991).
  19. Kalyanasundaram, N., Ghazisaeidi, M., Freund, J. B. & Johnson, H. T. Single impact crater functions for ion bombardment of silicon. Applied Physics Letters 92, 131909 (2008).
  20. Costantini, G., de Mongeot, F. B., Boragno, C. & Valbusa, U. Is ion sputtering always a “negative homoepitaxial deposition”? Physical Review Letters 86, 838 (2001).
  21. Norris, S., Brenner, M. & Aziz, M. From crater functions to partial differential equations: a new approach to ion bombardment induced nonequilibrium pattern formation. Journal of Physics: Condensed Matter 21, 224017 (2009).
  22. Aziz, M. Nanoscale morphology control using ion beams. Matematisk-fysiske Meddelelser 52, 187 (2006).
  23. Enrique, R. & Bellon, P. Compositional patterning in systems driven by competing dynamics of different length scale. Physical Review Letters 84, 2885–2888 (2000).
  24. Zhou, H., Zhou, L., Özaydin, G., Ludwig Jr., K. & Headrick, R. Mechanisms of pattern formation and smoothing induced by ion-beam erosion. Physical Review B 78, 165404 (2008). Note: this reference contains an error in the calculation of . As seen in the supplemental material, using a framework involving local co-ordinates and separation of large and small spatial scales, one obtains for this value , and not the stated value of .
  25. Headrick, R. & Zhou, H. Ripple formation and smoothening on insulating surfaces. JOURNAL OF PHYSICS: CONDENSED MATTER 21, 224005 (2009).
  26. Kalyanasundaram, N., Freund, J. & Johnson, H. A multiscale crater function model for ion-induced pattern formation in silicon. Journal of Physics: Condensed Matter 21, 224018 (2009).
  27. Wooten, F., Winer, K. & Weaire, D. Computer generation of structural models of amorphous Si and Ge. Physical Review Letters 54, 1392 (1985).
  28. Bazant, M. Z., Kaxiras, E. & Justo, J. F. Environment-dependent interatomic potential for bulk silicon. Phyiscal Review B 56, 8542 (1997).
  29. Samela, J., Nordlund, K., Keinonen, J. & Popok, V. Comparison of silicon potentials for cluster bombardment simulations. Nuclear Instruments and Methods in Physics Research B 255, 253 (2007).
  30. K. Nordlund, N. R. & Sundholm, D. Repulsive interatomic potentials calculated using hartree-fock and density-functional theory methods. Nucl. Instr. Meth. Phys. Res. B 132, 1–8 (1997).
  31. Nordlund, K. et al. Defect production in collision cascades in elemental semiconductors and fcc metals. Phyiscal Review B 57, 7556 (1998).
  32. Ghaly, M., Nordlund, K. & Averback, R. S. Molecular dynamics investigations of surface damage produced by kiloelectronvolt self-bombardment of solids. Philosophical Magazine A 79, 795 (1999).
  33. Samela, J., Kotakoski, J., Nordlund, K. & Keinonen, J. A quantitative and comparative study of sputtering yields in Au. Nuclear Instruments and Methods in Physics Research B 239, 331 (2005).
  34. Madi, C., George, H. & Aziz, M. Linear stability and instability patterns in ion-sputtered silicon. Journal of Physics: Condensed Matter 21, 224010 (2009).
  35. Facsko, S., Bobek, T., Stahl, A. & Kurz, H. Dissipative continuum model for self-organized pattern formation during ion-beam erosion. Physical Review B 69, 153412 (2004).
  36. Cross, M. & Hohenberg, P. Pattern formation outside of equilibrium. Reviews of modern physics 65, 851–1123 (1993).
  37. Madi, C., Anzenberg, E., Ludwig Jr., K., & Aziz, M. (2010). (unpublished).
  38. Henriksson, K., Nordlund, K. & Keinonen, J. Molecular dynamics simulations of helium cluster formation in tungsten. NUCLEAR INSTRUMENTS & METHODS IN PHYSICS RESEARCH SECTION B 244, 377 (2006).
  39. Mousseau, N. & Barkema, G. T. Binary continuous random networks. J. Phys.: Condens. Matter 16, S5183–S5190 (2004).
  40. von Alfthan, S., Kuronen, A. & Kaski, K. Realistic models of amorphous silica: A comparative study of different potentials. PHYSICAL REVIEW B 68, 073203 (2003).
  41. Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A. & Haak, J. R. Molecular dynamics with coupling to external bath. J. Chem. Phys. 81, 3684 (1984).
  42. Ziegler, J., Biersack, J. & Littmark, U. The Stopping and Range of Ions in Matter (Pergamon Press, New York, 1985).
  43. Orchard, S. On surface levelling in viscous liquids and gels. Appl. Sci. Res. 11A, 451 (1962).
  44. Vauth, S. & Mayr, S. Relevance of surface viscous flow, surface diffusion, and ballistic effects in kev ion smoothing of amorphous surfaces. Physical Review B 75, 224107 (2007).
  45. Eaglesham, D. J., White, A. E., Feldman, L. C., Moriya, N. & Jacobson, D. C. Equilibrium shape of si. Physical Review Letters 70, 1643 (1993).
  46. Gnaser, H. Low-energy ion irradiation of solid surfaces, vol. 146 of Springer tracts in modern physics (Springer, 1999).
  47. Wallner, G. et al. Defect production rates in metals by reactor neutron irradiation at 4.6 K. Journal of Nuclear Materials 152, 146 (1988).
  48. Mayr, S. G., Ashkenazy, Y., Albe, K. & Averback, R. S. Mechanisms of radiation-induced viscous flow: Role of point defects. Physical Review Letters 90, 055505 (2003).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description