Linear vs. Non-linear Behaviour in Ion Irradiation Nanostructuring of Nickel and Silicon Surfaces

Spontaneous nanometre-scale quasi-periodic ripple-like structures are formed at the surface of polycrystalline Ni films and Si(111) single crystal wafers by irradiation with a broad Ar+ ion beam at room temperature and studied with Atomic Force Microscopy as a function of fluence. The development of these structures can be reproduced by numerical solution of a continuum equation describing the evolution of surface morphology under ion irradiation, using realistic coefficients derived from material properties. In particular, we demonstrate that differences observed in pattern formation on the two surfaces under the conditions studied, such as wavelength stability and exponential growth of interface width for the Ni surfaces compared with wavelength coarsening and interface width saturation on Si(111), can be understood in terms of a cross-over between linear and non-linear behaviours.


Introduction
Irradiation of a solid surface by a beam of ions with sufficiently high kinetic energy can lead to the removal of atoms from its surface, a process known as ion beam sputtering.Under certain conditions, it has been found that ion irradiation can lead to the formation of self-organized regular patterns such as nano-ripples (one-dimensional quasi-periodic nanostructures), nano-dots, and nano-holes (zerodimensional nanostructures) [1][2][3].Such self-organized structures can be used for a range of applications, including templates for the growth of metal nanowire and nanodot arrays [4,5], a bioactive surfaces [6], patterns for evaluation of scanned probe microscopy tips [7], luminescent emitters [8], and highly optically absorbing/anti-reflective surfaces [9,10].Producing nanostructures by this approach has several advantages in comparison with focused ion beam patterning [11], electron beam lithography and optical lithography [12] due to its low cost and flexibility.In particular, ion-beam induced nanopatterning involves a single step and only requires a broad ion beam, with the result that alignment and focusing are not significant issues.
Unfortunately, there are some drawbacks associated with this process, such as a lack of predictability as to which conditions are appropriate for the formation of a desired structure of selected pitch/wavelength, as our understanding of the mechanism of self-organisation is far from complete.An in-depth grasp of how to optimise primary experimental parameters such as the ion energy, incidence angle, dose, flux, sample temperature, sample type (i.e.amorphous, polycrystalline or crystalline) and sample rotation during irradiation is a key requirement for future application.Such knowledge requires continued development of the theoretical basis for the formation process of these induced self-organized surfaces.
Significant theoretical efforts have been expended to explain ion bombardment-induced periodic nanostructures on various substrates.The first successful theoretical approach for the understanding of induced self-organized ripple formation was formulated by Bradley and Harper [13] (known as the BH continuum model).The linear BH partial differential equation describes the development of selforganised ripple structures as resulting from the balance between two competing mechanisms: sputterinduced surface roughening, which can be understood as producing a negative effective surface tension, and smoothing due to an isotropic thermal surface relaxation originating from thermally-induced diffusion.Although this linear continuum theory was successful in the prediction of many features of ripple formation, it has several shortcomings: (1) it does not predict the formation of ripples at low temperature, where thermal diffusion is negligible (e.g., on most surfaces at room temperature); (2) it does not correctly foretell the surface evolution at long time scales, in particular the experimentally observed power-law coarsening of ripple wavelength; (3) it does not predict the saturation of the interface width, w, at high ion fluence; (4) it does not predict the observed rotation of ripple orientation by 90° at high ion dose and temperature; (5) and it cannot explain the formation of nano-dots (quantum dots) and nano-holes [7,8,[14][15][16].
To address these failings several theoretical models have been developed which extend the BH model [17][18][19][20].In particular, Makeev, Cuerno and Barabasi (MCB) [20] derived a continuum equation to describe the self-organisation process by considering the microscopic elementary processes occurring during sputter-induced surface morphology evolution which led to inclusion of higher-order derivatives of the local surface height.This model successfully explained the formation of quasi-periodic structures at temperatures below which the thermal surface relaxation mechanism can play a sufficient role to counterbalance roughening through introducing a non-thermal effective ion-induced diffusion [21].In this model, the temperature-independent fourth derivative terms in the description of the rate of evolution of the surface topography (i.e., local surface height) play the same role of surface smoothing as conventional diffusion but without involving any mass transport.The continuum equation developed by Makeev and co-authors addressed many of the limitations of the BH continuum equation such as successfully describing the rotation of ripple orientation at high incident beam angle and high fluence, ripple saturation, and ripple formation at room temperature.
Although considered for the case of nanodot formation [22], to the best of our knowledge there is currently no direct comparative work between the results of experiment and theory regarding ion irradiation-induced ripple formation presented in the literature, which we address in this work.The experimentally-observed evolution of the surface topography of polycrystalline Ni thin films and Si(111) surfaces under inert gas ion irradiation as a function of ion fluence are compared with the MCB model, which we solve numerically.In particular, the coefficients for the various terms in the MCB equation are derived either from experimentally determined parameters or the results of simulations of the sputtering process using the Stopping and Range of Ions in Matter (SRIM) packages [23].Calculations of interface width and wavelength show the same general trends as the experimental data.Moreover, we find that, over the range of irradiation conditions studied the experimental results, particularly ripple wavelength and interface width, can be modelled in purely linear terms for the case of Ni but require the inclusion of non-linear terms if the behaviour of the Si surface is to be correctly described.These differences can be explained by the values of the characteristic "cross-over time",y , between linear and nonlinear regimes which we find numerically for both materials from the Makeev model.It would be expected that the higher surface erosion rate and lower ion penetration depth in Ni compared with Si would favour non-linearity in the former case, however the difference in incident ion beam angle between experiments on the two surfaces is found to dominate, indicating the delicate balance of factors governing ion-induced pattern formation and the need to carefully consider all experimental variables when comparing measurements to theoretical predictions in these systems.

Methods
Samples consisting of virgin Si(111) single crystal wafers and 100 nm thick polycrystalline Ni films deposited by thermal evaporation on Si(111) substrates were cut to 1×1 cm 2 and introduced into a homebuilt ultra-high vacuum (UHV) system (base pressure 1×10 -7 mbar) at Durham University where they were irradiated with a defocussed Ar + beam produced by a cold cathode ion source (Vacuum Generators AG21) such that the ion beam uniformly illuminated the surface of the samples.Ion bombardment was undertaken at ion energy, Eion, of 5.0 keV, ion flux, f, of 2.1×10 13 ions cm -2 s -1 and an incident angle,  of 75° with respect to the surface normal for the Ni surfaces.For Si(111) the irradiation conditions were Eion = 3.8 keV, f = 1.6×10 13 ions cm -2 s -1 and  = 60°.In both cases, a range of ion fluences was studied by varying the total irradiation time at fixed fluence.The chamber pressure typically increased to 1×10 -5 mbar during ion bombardment and the current incident upon the sample was measured directly using a Keithley 445 electrometer.
After ion bombardment samples were removed from the UHV system and imaged, in Durham, by Atomic Force Microscopy (AFM) under ambient atmosphere using a Veeco NanoScope IIIa scanned probe microscope in tapping mode with sharp tapping etched silicon probes (TESP).Care was taken to image all samples at a range of scan angles, scan speeds and feedback loop conditions to confirm that all images of surface structures observed reflected the topography of the surfaces rather than experimental artefacts.Multiple different areas on each samples were imaged to ensure unformity.After imaging, analysis of the surface topography of both Ni and Si(111) samples was undertaken using Gwyddion [24] from which a number of key parameters were extracted, such as interface width and average surface height.
Surface structural evolution was explored using a special case of the MCB model.The linear erosion rate expected for ion sputtering of a planar surface, -, along with first and third-order terms which contribute only to the motion of the quasi-periodic ripples parallel to the surface, but do not otherwise influence surface morphology [20], are neglected.The neglect of these terms is equivalent to a coordinate transformation into a frame of reference moving at the same velocity as the eroding surface along the (global) surface normal and at the velocity of propagation of the bombardment-induced surface features in the plane perpendicular to the normal and leads to equation 1, which is similar to the noiseless anisotropic Kuramoto-Sivashinsky equation [ where vx, and vy are effective surface tension coefficients in the x and y directions respectively.Dxx and Dyy are the orthogonal effective ion-induced surface diffusion coefficients and Dxy is the diagonal effective ion-induced surface diffusion coefficient.λx, λy are lateral coefficients in x and y directions respectively.Values for these parameters were determined for the specific cases of Ni and Si using results determined from SRIM [23] simulations and the defintions of Makeev, Cuerno and Barabasi [20] to determine sputter yield, average ion penetration depth, longitudinal and transverse ion straggles.
For negative values of νx, vy and positive values of Dxx, Dyy periodic ripples will form with wavevector oriented along the x or y direction.The characteristic wavelength, ℓ of these ripples is given by: where kmax is the fastest growing wavevector.Numerical integration was performed using Exponential Time Differencing (ETD) with Runge-Kutta accuracy of order 4 (ETDRK4).Discretization of equation 1 was performed using the Spectral Fourier Method.Since equation 1 is categorized as a stiff semi-linear partial differential equation, we choose a small temporal increment Δt = 0.01 to get a stable numerical solution and impose periodic boundary conditions h(x, y, t)= h(x + Lx, y, t) = h(x, y + Ly, t), where the domain size is Lx×Ly=180×180 nm 2 .The initial surface was selected to be a (Gaussian) randomly rough surface with a local surface height uniformly distributed between [-0.04, 0.04] nm with a zero average.The computational codes were written in Python and run on the Hamilton supercomputer at Durham University.

Results and Discussion
AFM topographs illustrating the evolution of the surfaces of the polycrystalline Ni thin film and Si(111) surface with ion fluence are presented in Figure 1 and Figure 2. The direction of the incident ion beam, projected onto the sample surface, is illustrated by an arrow in each case.For both surfaces the formation of one-dimensional quasi-periodic structures -or ripples -is clearly apparent.Two possible ripple orientations can be predicted theoretically [13,20] known as "parallel mode" and "perpendicular mode" ripples, depending on the orientation of the ripple wavevector with respect to the projection of the ion beam direction onto the surface.Defining the direction of the wave-vector of the ripple structures, k max , at right angles to the average direction of the "wave-fronts" it can be seen that perpendicular-mode ripples form in both cases.
As fluence increases the wavelength, ℓ, of the ripples formed on the Ni film surface initially remains constant, as shown in Figures 1(a-c) and Figure 3(a), followed by a sudden and abrupt increase in value, Figure 1(d) which is accompanied by a reduction in interface width, Figure 3(c), which otherwise increases exponentially.Examination of the thin film samples post-irradiation indicated that at fluences beyond approximately 5×10 16 ions cm -2 the 100 nm thick Ni film was completely eroded and the transition in wavelength and interface width reflects the formation of periodic structures in the Si(111) substrate supporting the Ni film.In contrast to the constant wavelength of the self-organised ripples formed on polycrystalline Ni, ℓ is seen to change significantly during irradiation of the Si(111) surface, coarsening as a function ion fluence, Figure 2(a-d) and Figure 3(b).The variation of wavelength with ion fluence can be fitted with a power-law dependence ℓ ~  n (red curve, Figure 3(b)) with a coarsening exponent n ≈ 0.2.Likewise, the evolution of the interface width on the bare Si(111) surface (Figure 3(d)) is very different to the exponential increase observed for polycrystalline Ni and is observed to saturate at large fluence.A somewhat similar saturation of interface width by Ziberi and co-workers [14] for irradiation of Si(100) surfaces.The similarity between the behaviours of the (100) and (111) surfaces is to be expected due to amorphization under ion bombardment, in comparison with the facet-dependent structural evolution seen on some metal surfaces [1].Note that we required a rather higher range of fluence to observe structures on the Si(111) surface in comparison with the polycrystalline Ni films and the interface width remains small.The behaviour of the surface topography of the polycrystalline Ni film under ion irradiation, at least until the Si(111) substrate is reached, is consistent with the linear BH model which predicts a constant ripple wavelength and exponentially increasing interface width [13].However, the ripple coarsening observed for irradiation of the bare Si(111) surface is similar to that observed for medium energy ion irradiation of Si(100) in the work of Datta and Chini [26] and requires the inclusion of non-linear terms in the description of surface evolution under irradiation, as does the saturation of interface width shown in Figure 3(d).To explain the differences in the evolution of wavelength and interface width between the two materials, and in particular to understand the differences in the relative importance of non-linear terms between the two systems the time-evolution of the local surface height, h(x,y,t) was by numerical integration of equation 1 using parameters derived for both the experimentally studied systems, Ni and Si, as described above.The sputtering yields for Ni and Si surfaces as a function of incident ion beam energy and angle, obtained by carrying out SRIM simulations [23], were used to produce erosion rate factors, F, which are corrected from those presented by Makeev et al. [20] to accurately describe the variation of sputter yield at high angles found from the SRIM simulations.The full set of parameters used for equation 1, derived for each of the materials, is shown in Table 1.
Table 1.Calculated coefficients used in the numerical solution of the MCB model.

Ni Si
Target density, ρN (atoms/cm 3 ) 9.1×10 22   Numerical solutions for surface topography as a function of iteration (equivalent to increasing fluence) for the Ni surface are presented in Figure 4 and those for Si(111) in Figure 5.For the Ni surface, the ripple wavelength, ℓ, remains constant even after approximately 10,000 iterations, the increasing number of computational cycles up to this point leading only to more clearly-defined and ordered structures.Figure 6(a) shows, for the Ni surface the evolution of the ripple wavelength with computational iterations for selected iteration times.The wavelength of the ripple features is initially constant then, after a certain computational time begins to coarsen.This behaviour can be understood in terms of a regime of linear behaviour below a "cross-over time", y ~1 × 10 6 iterations (or equivalently, for fixed flux, a cross-over fluence, c), in which the contribution of non-linear terms in equation (1) are negligible, followed by a regime dominated by non-linear terms for t > y [27].However, in our experiment the latter case can only be reached at fluences beyond that at which we completely erode the Ni film, Figure 1 and Figure 3(a).
The behaviour of the interface width, w, and the average interface height, ℎ ̅ (), is shown in Figure 6(c).ℎ ̅ () is defined as: where ℎ(, , ) is the local interface height at (x,y) and time t and ℎ  () is the (decreased) height which would be expected due to sputter erosion of a uniform planar surface under the same irradiation conditions.Below y the interface width, w, is seen to increase exponentially and ℎ ̅ () remains zero, before w saturates and ℎ ̅ () decreases rapidly at larger numbers of iterations.The decrease in ℎ ̅ () can be understood in terms of an enhancement of erosion rate in comparison with a planar surface when non-linear terms in the MCB equation become dominant and the surface coarsens.Such variations in sputter yield have been observed theoretically by Makeev and Barabasi for both rough surfaces [28] and those with ion irradiation-induced periodic structures [29].A similar computational analysis of the evolution of the Si(111) surface, Figures 6(b) and 6(d), shows qualitatively similar behaviour of ℓ, w and ℎ ̅ (), with the key difference that the cross-over time between the dominance of linear and non-linear terms occurs earlier at   ~6 × 10 5 iterations.The earlier crossover time for the Si(111) surface coupled with the higher range of fluences employed experimentally (1.7×10 17 ions cm -2 to 4.6×10 17 ions cm -2 for Si(111) compared with 1.1×10 16 ions cm -2 to 3.5×10 16 ions cm -2 for polycrystalline Ni) explains the significant differences in the behaviour of interface width and wavelength observed.
where   ,   , a and F are defined in Table 1 and Z is a function of a, ion beam incidence angle, , and   ,   which are the ratios of the longitudinal and lateral ion straggles to the average ion range, respectively.From the examination of the material-dependent terms a and F it might be initially thought that the cross-over time for Si should be greater than that of Ni due to the larger average ion penetration depth and smaller erosion rate factor (Table 1).This would certainly be the case if the experiments were performed with identical incident ion beam angles.However, the difference of 15º between the two sets of measurements and the resulting change in Z is enough to reverse the order at which the cross-over time is reached, indicating the delicate dependence of ion beam-induced pattern formation on the conditions of irradiation as well as the intrinsic properties of the surfaces and the ions incident upon them.
A comparison of the experimentally determined values of ripple wavelength for the polycrystalline Ni surface, Figure 3(a) with the linear region of the computationally-determined wavelength for the same surface, Figure 6(a) shows that the latter is systematically smaller by a factor of ~5.Likewise, a comparison of the experimental interface width obtained on the silicon surface at saturation, Figure 3(d), with the computationally determined value, Figure 6(d) shows that the latter is a factor of ~50 smaller.Given the explicit use of materials parameters to determine the coefficients of equation (1) the disagreement between the calculated evolution of the silicon and nickel surfaces under ion bombardment and those observed experimentally point to either inaccuracies in determining the ion straggle when using TRIM, as suggested by Facsko and co-workers [22], or additional processes which have not been explicitly incorporated into the description modelled by equation (1).We can expect that our neglect of thermally induced diffusion at room temperature is unlikely to impact significantly in describing surface evolution at room temperature.However, the current model neglects the influence of viscous surface relaxation and flow [20,30,31] and ion-induced mass transport [32] which, through changing the balance between competing processes of roughening and smoothing which lead to ripple formation, would influence the theoretically determined values for ripple wavelength and interface width.

Conclusions
Experimental observations of the formation of quasi-periodic patterns (ripples) due to Ar + ion irradiation of the surfaces of two different materials, polycrystalline Ni and Si(111), show significant differences.In particular, ripples formed on polycrystalline Ni display a wavelength which is constant with fluence and an interface width (amplitude) which increases exponentially, whereas the structures formed on Si(111) coarsen with increase irradiation and show amplitude saturation.
Numerical simulations performed by solving the continuum equation for surface evolution under ion irradiation derived by Makeev, Cuerno and Barabasi (MCB) using realistic coefficients derived from the physical properties of the two materials and the results of calculated ion irradiation parameters show a good agreement with experimental trends.The simulations indicate that, under the conditions studied, the differences between experimentally observed pattern formation on polycrystalline Ni and Si(111) surfaces can be understood in terms of differing values of the cross-over time from behaviour dominated by linear to that controlled by non-linear terms in the MCB equation [20,27].Moreover, our results indicate the subtle interplay between materials properties, ion beam and incident beam angle in determining the fluence (or time) evolution of ripple structures under ion bombardment.
Numerical modelling using the MCB equation with coefficients derived for specific physical systems and experimental validation as reported in this work paves the way for quantitative predictions of the evolution of ripple-like surface structures produced by ion irradiation.However, it is clear from the differences between experimentally and theoretically determined values of wavelength and interface width that further work is required to fully capture all contributions to morphological evolution in a quantitative continuum model.Fully understanding the conditions for the formation of specific surface structures with selected amplitude and wavelength should ultimately enable ion bombardment to be used as an effective tool for the application-specific tailoring of surfaces.

Figure 3 .
Figure 3. Evolution of ripple wavelength for (a) polycrystalline Ni and (b) Si(111) surfaces and interface width, w, for (c) polycrystalline Ni and (d) Si(111) as a function of ion fluence.Note the sudden transition in both wavelength and interface width in (a) and (c) at the point the Ni film is completely eroded and the Si substrate reached.The ripple wavelength on Si(111) shows a power law dependence ℓ ~  n , with n ≈ 0.2, as shown by the red line in (b).In the early stages of ripple formation on polycrystalline Ni the increase in interface width increases exponentially, as shown by the fitted line in (c).The line in (d), indicating saturation of the interface width on Si(111) is a guide to the eye.Irradiation conditions were: ion energy E ion = 5.0 keV, ion flux f = 2.1×10 13 ions cm -2 s -1 , incident angle θ = 75º and room temperature for Ni; ion energy Eion = 3.8 keV, ion flux f = 1.6×10 13 ions cm -2 s -1 , ion angle θ = 60º and room temperature for Si(111).

Figure 4 .
Figure 4. Simulated surface evolution during ion bombardment of Ni, under the same conditions as those employed experimentally, shown at selected computational iterations (cycles): (a) 5.0×10 4 , (b) 2.0×10 5 , (c) 1.0×10 6 , and (d) 1.9×10 6 iterations.The projection of the ion beam onto the target surface is indicated by the arrow in (a) and is the same in all panels.

Figure 5 .
Figure 5. Simulated surface evolution during ion bombardment of Si, under the same conditions as those employed experimentally, shown at selected computational iterations (cycles): (a) 5.0×10 4 , (b) 2.0×10 5 , (c) 6.0×10 5 , and (d) 1.9×10 6 iterations.The projection of the ion beam onto the target surface is indicated by the arrow in (a) and is the same in all panels.

Figure 6 .
Figure 6.Results from simulations of surface evolution: Ripple wavelength as a function of computational time (iterations) for (a) Ni and (b) Si; Variation of mean surface height, ℎ ̅ (), and interface width, w, for (c) Ni and (d) Si. 25]: