Morphology change of the silicon surface induced by Ar$^+$ ion beam sputtering

Two-level modeling for nanoscale pattern formation on silicon target by Ar$^+$ ion sputtering is presented. Phase diagram illustrating possible nanosize surface patterns is discussed. Scaling characteristics for the structure wavelength dependence versus incoming ion energy are defined. Growth and roughness exponents in different domains of the phase diagram are obtained.


Introduction
It is well known that low and medium energy ion sputtering may induce a fabrication of periodic nanoscale structures on an irradiated surface [1]. Depending on the sputtered substrate characteristics and sputtering conditions, different types of nanoscale structures such as ripples, nanoholes and nanodots can grow on a target during ion beam sputtering [2][3][4][5][6]. These patterns have been found on both amorphous and crystalline materials including insulators, semiconductors and metals (see reference [7] and citations therein). Main theoretical models describing ripple formation are based on the results of famous works by Bradley and Harper [4], Kardar et al. [8], Wolf and Villian [9], and Kuramoto et al. [10]. The main control parameters in these models reduced to surface tensions, tilt-dependent erosion rates and diffusion constants are determined by sputtered substrate characteristics and sputtering conditions (see for example [7]).
Among theoretical investigations there are a lot of experimental data manifesting a large class of patterns formed due to the self-organization process. It was experimentally shown that the main properties of pattern formation processes depend on ion-beam parameters such as ion flux, energy of deposition, angle of incidence, and temperature of the substrate (target). Therefore, to study the ion beam sputtering processes theoretically one needs to determine the mentioned parameters of the model according to the physical conditions related to concrete materials.
One of the most frequently used materials for ion beam sputtering is silicon because it is the mainstream material in modern microelectronic industry and it is readily available with high purity and quality. Nanostructuring of silicon has received much attention due to its potential application in developing the Si light sources [11]. Various techniques such as acid etching, ion implantation, reactive evaporation, chemical vapor deposition and molecular beam epitaxy have been used in developing the Si nanomaterials (porous Si, Si nanocrystal-doped dielectrics and Si quantum dots) (see reference [11] and citation therein).
In this paper we study the properties of the formation of nanoscale patterns on a silicon target sputtered by Ar + ions. To this end, we use a two-level scheme, based on Monte-Carlo simulations and the modified Bradley-Harper theory. In the first approach we compute the ion energy dependent penetration depth, widths of the ion energy distribution and sputtering yield. Next, we exploit these characteristics as input data for the continuum approach describing the evolution of the surface height field. We define the domains of values for the angle of incidence and ion energy where different nanoscale structures can be formed. The dynamics of nanoscale pattern formation is discussed. We will show that at fixed values for the incidence angle, one has two scaling exponents for wavelength related to small and large values for ion energy. In addition, we obtain roughness and growth exponents.
The work is organized in the following manner. In section 2 we present the theoretical model in the framework of the modified Bradley-Harper approach. In section 3 using Monte-Carlo modeling we compute the main characteristics incorporated into the continuum theory. The related phase diagram, the dynamics of the formation of nanoscale structures, the ion energy dependent wavelength and scaling properties of the surface morphology are discussed in section 4. We conclude in the last section.

Theoretical model
Let us consider a d-dimensional substrate and denote with r the d-dimensional vector locating a point on it. The surface is described at each time t by the height z = h(r, t). If we assume that the surface morphology changes during ion sputtering, then we can use the model for the surface growth proposed by Bradley and Harper [4] and further developed by Cuerno and Barabasi [2]. We consider the system where the direction of the ion beam lies in x − z plane at an angle θ to the normal of the uneroded surface. Following the standard approach one assumes that an averaged energy deposited at the surface (let say at point O), due to the ion arriving at the point P in the solid, follows the Gaussian distribution [4] ε denotes the kinetic energy of the arriving ion, σ and µ are the widths of distribution in directions parallel and perpendicular to the incoming beam. Parameters σ and µ depend on the target material and can vary with physical properties of the target and incident energy. The erosion velocity at the surface point O is described by the formula v = p R drΦ(r)E(r), where integration is provided over the range of the energy distribution of all ions; here Φ(r) are the corrections for the local slope dependence of the uniform flux J. The material constant p is defined as: where U 0 and C 0 are the surface binding energy and the constant proportional to the square of the effective radius of the interatomic interaction potential, respectively [12]. The general expression for the local flux for surfaces with non-zero local curvature is [13]: Hence, the dynamics of the surface height is defined by the relation ∂ t h ≃ −v θ − ∇ x h, ∇ 2 x h, ∇ 2 y h and is given by the equation ∂ t h ≃ −v(θ) 1 + (∇h) 2 , where 0 < θ < π/2 [2-4, 8, 14]. The linear term expansion yields Here v 0 is the surface erosion velocity; γ = γ(θ) is a constant that describes the slope depending erosion; ν α = ν α (θ) is the effective surface tension generated by erosion process in α direction.
If one assumes that the surface current is driven by differences in chemical potential µ, then the evolution equation for the field h should take into account the term −∇ · j s in the right hand side, where j s = K∇(∇ 2 h) is the surface current; K > 0 is the temperature dependent surface diffusion constant. If the surface diffusion is thermally activated, then we have K = D s κρ/n 2 T , where D s = D 0 exp (−E a /T ) is the surface self-diffusivity (E a is the activation energy for surface diffusion), κ is the surface free energy, ρ is the areal density of diffusing atoms, n is the number of atoms per unit volume in the amorphous solid. This term in the dynamical equation for h is relevant in high temperature limit which will be studied below.
Assuming that the surface varies smoothly, we neglect spatial derivatives of the height h of third and higher orders in the slope expansion. Taking into account nonlinear terms in the slope expansion of the surface height dynamics, we arrive at the equation for the quantity h ′ = h + v 0 t of the form [2,4] where we drop the prime for convenience. Here we introduce the uncorrelated white Gaussian noise ξ with zero mean mimicking the randomness resulting from the stochastic nature of the ion arrival 23602-2 to the surface. In equation (1) the effective surface tensions ν x and ν y generated by the IBS, the tilt-dependent erosion rates λ x and λ y are defined through the incident angle θ, penetration depth of incident ion a, distribution widths σ, µ and the sputtering yield Y 0 as follows [2,4]: Here we have used the following notations: Let us perform the stability analysis for a system with additive fluctuations. To this end, we average the Langevin equation (1) over noise and obtain Considering the stability of the smooth surface characterized by h = 0, we can rewrite the linearized evolution equation in the standard form: with notationsν It is easy to see that equation (7) admits a solution of the form h = A exp[i(k x x+k y y−ωt)+χt]. Indeed, substituting it into equation (7) and separating real and imaginary parts we get As far as F 0 , f , a are positive values, hence one has ν y < 0, whereas ν x can change its sign. Therefore, the Bradley-Harper model does not provide for stable smooth surface. Hence, we can conclude that if ν x > 0, then ripples (wave patterns) appear in x-direction. On the contrary, when ν x < 0, equiaxed structures (nanodots/nanoholes) can be formed on an eroded surface. In addition, the sign of the product λ x · λ y can play a crucial role in ripple formation processes [15].
For the noiseless nonlinear model (1) it was shown that as the sets ν α and λ α are the functions of the angle of incidence θ ∈ [0, π/2] there are three domains in the phase diagram (a σ , θ) where ν x and λ x change their signs, separately [2]. This results in the formation of ripples in different directions x or y varying a σ or θ.
One needs to note that in the Bradley-Harper approach describing the processes of ripple formation on amorphous substrates the penetration depth can be approximated as a(ε) ∼ ε leading 23602-3 to the power-law asymptotics for the wavelength of the ripples Λ versus ion energy as follows: [4]. In the next section, performing calculations for the sputtering of the silicon target by Ar + ions we shall show that there are deviations from these asymptotics due to power-law dependence of a, σ, µ versus ion energy ε. Moreover, it will be shown that the sputtering yield depends on both incident angle θ end ion energy ε in a power-law form. To define a, σ, µ and the sputtering yield, we shall use Monte-Carlo approach.

Monte-Carlo modelling
To study the evolution of the silicon surface morphology during Ar + ion beam sputtering one needs to know the energy of ions and target characteristics such as: penetration depth of the Ar + ions into the silicon target a; widths of the distribution in parallel and perpendicular directions of the incoming beam (σ and µ) and sputtering yield Y 0 . Moreover, to simulate the target morphology evolution one should define temperature T , uniform flux J, atomic density of the target N , surface binding energy U 0 and the effective radius of interatomic interaction potential. Using data from reference [5] at T = 550 C one has K = C ′ 8.49 × 10 3 nm 4 /s, where addimer concentration is C ′ = 0.04 atoms/site (≈ 4% coverage), or, C ′ = 0.07 atoms/nm 2 . To define the material constant p we shall use N ≃ 50 atoms/nm 3 , U 0 = 4.73 eV. For the effective radius of interatomic interaction potential, we put the length of the main diagonal of silicon primitive cell with lattice parameter 0.5437 nm. To compute the time evolution of the silicon surface morphology, we should calculate the effective surface tensions ν x and ν y generated by the IBS, and the rates λ x and λ y [see equation (2)]. Following relations (3) and (4) dependent on ion energy, parameters a, σ and µ, as far as Y 0 = Y 0 (θ, ε), can be computed. From experimental point of view, the control parameters at IBS are the energy of ion beam, off-normal incidence angle and ion flux. In all our calculations we put J = 20 ions/(nm 2 s). Thereafter we vary the ion energy in the interval 100 eV÷10 keV and In further study, we use the well-known program codes (TRIM and SRIM) to calculate the stopping range of ions in matter and transport range of ions in matter. Description of algorithms and the basic principles for Monte-Carlo calculation of both the transport range of ions in matter and the stopping range of ions in matter can be found in [16]; TRIM and SRIM codes can be found on the web-site www.srim.org.
Values for parameters a, σ and µ for silicon target sputtered by Ar + ions were obtained with the help of SRIM code (a program for calculating the stopping range of ions in matter). Results for relative penetration depths a σ ≡ a/σ and a µ ≡ a/µ versus ion energy are shown in figure 1 (dependencies a(ε), σ(ε) and µ(ε) are shown in the insert). In figure 1 it is seen that longitudinal and transverse widths σ and µ, respectively, as far as a σ and a µ satisfy the following relations, are as follows: σ < µ and a σ < a µ . In reference [17] it was shown that the rotated ripple structures formed when λ x λ y < 0 with rotation angle ϕ = tan −1 −λ x /λ y can be observed at small incidence angles θ when a µ < a σ (a σ = 1) and at intermediate and large θ when a µ > a σ . Hence, one can expect the appearance of rotated ripple structures in our system. It is principally important that dependencies of penetration depth a and longitudinal and transverse widths σ and µ, respectively, versus ion energy ε deviate from the linear law predicted by Bradley and Harper [4]. For a silicon target sputtered by Ar + ions we have obtained a power-law approximation of the form: φ(ε) = A + Bε C , where φ = {a, σ, µ}, constants A, B and C are fitting parameters. So, we can expect that the wavelength dependence Λ(ε) ∼ ε −δ can be characterized by the exponent δ = 1/2. In our further continuum approach we shall use the obtained power-law asymptotics for a, a σ and a µ from Monte-Carlo simulations.
To compute the dependence of the sputtering yield versus ion energy and angle of incidence we use Monte-Carlo approach realized in TRIM code (program for the calculation of transport range of ions in matter). The results of calculations for sputtering yield versus incident angle at fixed ion energy and sputtering yield versus ion energy at fixed incident angle are shown in figures 2 (a) and 2 (b), respectively. In figure 2 it is seen that sputtering yield depends on both ion energy and incidence angle in accordance with a power law as follows Y 0 (ψ) = A ′ + B ′ ψ C ′ , where ψ = {θ, ε}, constants A ′ , B ′ and C ′ are fitting parameters. Therefore, all parameters (ν x , ν y , λ x , λ y ) required to monitor the time evolution of silicon surface morphology during IBS are well defined.

Phase diagram and typical patterns
Firstly, let us compute a phase diagram ε(θ) defining domains for different surface patterns of silicon sputtered by Ar + ions. To this end, we shall monitor a sign change of surface tension ν x and tilt-dependent erosion rates λ x and λ y (as it was mentioned above ν y is always less than 0). The corresponding phase diagram indicating possible patterns is shown in figure 3. We need to stress that in the related interval for both the ion energy and the incidence angle except ν y < 0, one has λ y < 0. From figure 3 it follows that plane (θ, ε) is divided by three curves into five domains A, B, C, D and E. If one crosses the dash-dot curve, then quantity ν x changes it sign. Therefore, in the linear regime at small incidence angles θ (domain A), instability of the silicon surface occurs in both x and y directions due to ν y < 0 and ν x < 0. In the domain E (at large θ) in the linear regime, patterns are stable in x-direction due to ν x > 0. At large times (nonlinear regime) the 23602-5 surface morphology is governed by nonlinear parameters λ x and λ y . Solid curve in figure 3 divides domains characterized by λ x < 0 and λ x > 0. Therefore, between solid and dash-dot lines only λ x is positive (domains C and D), whereas in the domain E both ν x and λ x are positive. Dash curve corresponds to the condition ν x = ν y . Hence, before the dash curve (domains A and C) when ν x < ν y , vertical elongated surface structures should be formed, whereas after the dash curve (domains B, D and E), the corresponding structures should be of a horizontal elongated type. To illustrate typical structures in each domain in figure 3 we numerically solve equation (1) on quadratic lattice L × L of the linear size L = 256 with periodic boundary conditions. Spatial derivatives of the second and fourth orders were computed according to the standard finitedifference scheme; the nonlinear term (∇h) 2 was computed according to the scheme proposed in references [18,19]. We have used Gaussian initial conditions taking h(r, t = 0) = 0 and (δh) 2 = 0.1; the integration time step is ∆t = 0.005 and the space step is ℓ = 1.
Typical surface patterns in domains (A-E) are shown in figure 3. It is seen that on the left hand side of the solid curve when ν y < 0 and ν x < 0, pattern type of holes is realized (see snapshots A and B). It follows that patterns realized at high energy ions are characterized by small size (see snapshot A), whereas at small ε one has large-scale patterns (see snapshot B) 1 . Moreover, orientation of holes in points A and B is different. It is defined by a minimal value of both ν x and ν y . Structures, shown by snapshots C and D (ripples) are characterized by positive value of parameter λ x , which defines nonlinear effects in x-direction. An orientation of the corresponding ripples is defined by a minimal value of both ν x and ν y as in the previous case. Hence, as far as ν x < ν y from the left of the dashed curve in figure 3, the related patterns in domains A and C are elongated in y direction. On the contrary, in snapshots B, D, and E there are horizontal elongated structures. Structures in snapshots A, B, C and D are characterized by instabilities in both x and y directions due to ν x < 0 and ν y < 0. In the domain, indicated by point E due to ν x > 0, structures are stable in x direction.
The obtained phase diagram is in good correspondence with the results of experimental studies of the dynamics of the surface Si(001) sputtered by Ar + ions [20], where according to the experimentally obtained a phase diagram in the plane "ion energy -angle of incidence" it was shown that if the angle of incidence or ion energy varies, then orientation of ripples can be changed.
It is important that in the considered interval of incidence angle, the obtained phase diagram in figure 3 is topologically similar to the experimental one. However, nonlinear KS equation (1) with parameters defined by equations (2)-(4) does not presume a stable smooth surface because ν y is a negative quantity. To prove that holes and ripples are stable in time, let us consider the dynamics of the surface morphology change. We analyze two representative kinds of patterns shown in figure 3 as snapshots A and E and compute the number of islands for each pattern in time. To this end, we have cut the surface h(x, y) at an average height level h and calculated the relative number of islands N/N max at fixed times, where N max is a maximal value of islands. In our computation scheme we used the following definition for the island: all points on the surface with h < h belonging to one manifold having a closed boundary, form an island. The corresponding boundary of the island was obtained according to the percolation model formalism. Results for relative number of islands were averaged over 20 independent runs. Typical evolution of the number of islands is shown in figure 4 at ε = 2 keV for θ = 50 • and θ = 63 • . It is seen that the relative number of islands grows at small time interval that corresponds to processes of the formation of islands. At intermediate times, the relative number of islands decreases which means a realization of coalescence processes. It is important that in the process of ripple formation the coalescence regime is well pronounced (see empty circles). On the contrary, for the process of nanohole formation (filled circles), this such regime is only weakly observed. At large times one has a stationary behavior of the relative number of islands. Hence, processes of ripple and nanohole formation are stationary ones: at large time intervals the averaged number of islands does not change in time. Snapshots of the silicon surface morphology for θ = 50 • and θ = 63 • at t = 0, 40, 100 and 400 seconds are shown in figure 4 in the top and in the bottom of the figure, respectively. 23602-7

Wavelength dependence on the ion energy
Next, let us study the wavelength dependence on the incident ion energy and on the angle of incidence. As it was shown earlier in the Bradley-Harper theory, a relation between parameters ν x and ν y determines the orientation of surface patterns. The wavelength of selected patterns in the corresponding direction is defined as follows: Λ x,y = 2π 2K/|ν min x,y |, where ν min x,y = min(ν x , ν y ). One needs to note that following the phase diagram shown in figure 3, a variation in the ion energy at fixed angles of incidence causes a change in the orientation of structures: at small ε one has structures elongated in y direction, whereas at large ε, structures are horizontally elongated. Corresponding dependencies of the wavelength versus ion energy at fixed values for incidence angle are shown in figure 5. It is seen that the wavelength decreases with the ion energy growth according to a power law and varies in the interval from 100 nm to 1 µm. This result is in good correspondence with experimental data for sputtering of the silicon target by Ar + ions [5]. It is principally important that as far as the penetration depth depends on the ion energy in a nonlinear manner (see the insert in figure 1) one can expect a deviation from the Bradley-Harper wavelength asymptote Λ ∼ ε −1/2 . In figure 5 it is seen that at small and large incidence angles (see dot and dash lines, respectively) one has linear dependencies in log-log plot characterized by the corresponding unique slope. However, at intermediate values for θ related to the dash curve in figure 3, the dependence Λ(ε) has a kink. This kink means a change in the orientation of patterns. In such a case one has two slopes at small energies, i.e., before kink, one has selected the patterns characterized by Λ x , whereas at large energies the patterns are defined by Λ y (see solid line and asymptotics in figure 5). Therefore, at small ion energies the patterns are oriented in y direction, whereas at large ion energies they are oriented in x direction. Hence, for the wavelength dependence on the ion energy one can write Λ ∼ ε −δ where the scaling exponent δ is defined as a slope of the dependence Λ(ε) in double logarithmic plot before and after the kink. The dependence of the scaling exponent versus incidence angle is shown as an insert in figure 5. One can see that for the described interval for the angle of incidence δ > 1/2. Moreover, at small and large θ the exponent δ does not essentially change it values, whereas in the interval for θ when Λ(ε) has a kink, the exponent δ varies from 0.65 toward 1.05. We should note that the obtained picture is realized when the incoming ion flux J and temperature T are constants. In the opposite case, variation in J and T leads to the known asymptotes: where E a is an activation energy.

Scaling properties of patterns
Finally, let us study the scaling properties of the surface patterns, computing growth and roughness exponents. To this end, we analyze a height-height correlation function C h (r, t) = [h(r+ r ′ , t) − h(r ′ , t)] 2 . In the framework of dynamic scaling hypothesis following references [21,22], one arrives at scaling relations C h (t) ∝ t 2β , C h (r) ∝ r 2α , allowing one to define the growth exponent β and the roughness exponent α.
In reference [23] it was shown that there is a set of exponents {β} describing the universal behavior of the correlation function at early stages of the system evolution. At late times where a true scaling regime is observed there is a unique value for β. The roughness exponent α takes similar values at different time windows and can be considered as a constant depending on the system parameters only. From practical viewpoint, the analysis of the surface growth is urgent at large time intervals where the true scaling regime is observed and there is no essential difference in values β at different time windows. It is known that anisotropic surfaces studied in this paper may exhibit a more complex dynamic scaling behaviour than isotropic ones because anisotropy of the surface is reflected in lateral correlations of the surface roughness [24,25]. In reference [6] it was proposed to use local roughness scales α x , α y in the directions normal and parallel to the projection direction of the ion beam. Values for growth and roughness exponents together with surface tensions ν x and ν y at ε = 2 keV and fixed values for incidence angle are presented in table 1. It is seen that when ν x < ν y , a relation α x > α y is realized due to orientation of the structures in y direction. On the contrary, if ν x > ν y holds, then one has α x < α y . Hence, making an analysis of the obtained scaling exponents, one can conclude that if structures are oriented in y direction, then roughness is larger in x-direction and vice versa (compare patterns in snapshots A, C, D and E in figure 3 with exponents in table 1). The obtained results for growth and roughness exponents are in good correspondence with experimental studies of the silicon target sputtered by Ar + ions (see [6,26]).

Conclusions
Two-level modeling for nanoscale pattern formation on silicon target induced by Ar + ion sputtering has been reported. We have used Monte-Carlo simulations and a continuum approach based on the Bradley-Harper theory. It was shown that for the described system, the dependencies of the averaged penetration depth of the incident ion and the corresponding distribution widths of the deposited energy in directions parallel and perpendicular to the incoming beam versus ion energy are of the power-law form. Varying the incoming ion energy and ion incidence angle, we have defined the sputtering yield with the help of Monte-Carlo simulations. The obtained results have been used in the modified Bradley-Harper theory within the framework of two-scale modeling scheme.
We have computed a phase diagram for control parameters: i.e., incidence angle and ion energy that defines possible patterns on silicon target sputtered by Ar + ions. It was shown that at small incidence angles, nanohole patterns are realized, whereas at large incidence angles, pattern type of ripples is observed. Analyzing the morphology change of silicon surface we have shown that during the system evolution, the number of nanoholes/ripples becomes constant, indicating stability of the obtained structures in time.

23602-9
We have found that there are deviations from the Bradley-Harper asymptotics for the wavelength dependence on the ion energy. Moreover, when the orientation of patterns changes, a kink is realized in such asymptotics. The exponent of such power-law asymptotics depends on the angle of incidence. At fixed values for incidence angle one has two scaling exponents related to small and large values for the ion energy according to a change in the orientation of structures. While studying the scaling characteristics of the height-height correlation function, the growth exponent together with longitudinal and transverse roughness exponents are obtained for different values of incidence angle at a fixed ion energy. It was shown that relations between roughness exponents are defined through relations between corresponding effective surface tensions.
The results obtained in a two-scale modeling scheme are in good correspondence with the known theoretical and experimental data for sputtering of silicon target by Ar + ions in the considered interval of values for incidence angle of ions, the incoming ion energy, temperature and ion flux [5,6,20,26,27].