Simulations of Formation of Nanostructure in Silicon Surface by Single Slow Highly Charged Ion

To understand the mechanisms of surface erosion induced by slow highly-charged ion (SHCI) bombardment, the surface nanostructure formation in Si (111) surface by single Xe44+ ion was studied by using molecular-dynamics (MD) simulations, based on analyzing the multiple electron emission of the substrate. The time evolutions of the temperature, energy, pressure and density of the substrate have been systematically studied. The results show the bombardment of the incident SHCI resulting in an explosive event in the surface. A shock wave propagating at ~104 m/s is formed in the system during the initial 175 fs. After this initial shock, many particles are ejected from the surface since the extreme non-equilibrium of the system. And at t=370 fs, a crater-like nanostructure with diameter of ~40 Å and depth of ~18 Å is formed at the incident site.


Introduction
The surface nano-structures formation induced by slow highly charged ions (SHCI) has been extensively studied on many kinds of materials, from superconductors to insulators and from amorphous to crystalline materials [1][2][3][4]. These studies not only play an important role in understanding the formation mechanism of surface nano-structures, but also in various latest surface science and technology applications [5][6][7]. A large number of experiments show that the surfaces irradiated by SHCI present different morphology for different materials [5,[8][9][10][11][12], e.g. hillock [13,14] on highly oriented pyrolytic graphite (HOPG) surface and crater [15,16] on silicon surface. The type, size and formation of the surface nano-structures depend on the incident SHCI energy as well as on the material properties [17]. In order to clarify this interesting but puzzling mechanism, the surface nanostructure formation induced by collision of single slow Xe 44+ with Si(111) surface is studied in this work by using molecular-dynamics (MD) method. The time evolutions of the temperature, energy, pressure and density of the substrate during the nanostructure formation are studied, respectively.

Modeling and Simulation
In the present simulations, we prepared a Si(111) surface in the ground state, which is 162 Å×162 Å×37 Å in size, and consists of 50,365 atoms distributed over 28 layers. Non-periodic boundary conditions are used for the surface normal direction (Z direction), while periodic boundary conditions are applied to the two directions (X and Y directions) parallel to the surface. To maintain the surface stability, the speeds of the silicon atoms of the bottom 3 layers are set to zero in the simulations. The atoms of the remaining 25 layers are dynamical according to the force exerted from the neighbors. The deepest 3 dynamical layers, i.e., the 4th-6th layer of atoms from bottom to top were coupled to a thermostat at T 0 =300 K to maintain constant temperature. As a SHCI approaches a semiconductor surface, a lot of electrons are transferred resonantly from the surface to the SHCI. Many experiments show the total number of electrons pulled out of the surface is greater than the initial charge of SHCI [18][19][20][21]. If the charge neutrality of the surface domain that loses electrons can't be re-established on time, a positively charged region will be formed near the incident point on the surface, and the strong repulsion between the newly formed positive ions in the zone will overcome the binding force between them and produce a Coulomb explosion. Therefore, the following two issues have been analyzed before present simulations. i. Whether the characteristic charge neutralization time of the surface is much longer than the interaction time, which determines whether the Coulomb explosion model is suitable for the potential energy deposition of SHCI on the silicon surface. ii. How to determine the size of the positively charged region, which determines the intensity of the explosion and the size of the surface nano-structure. Therefore, we have analyzed the characteristic charge neutralization time by a solution of the static Maxwell equations, and found that it is ~1ps for Si, which is much longer than the interaction time ~10fs. And the Coulomb explosion model is applicable for describing the potential energy deposition of SHCI on the silicon surface in present simulations. And we assume that during the interaction between SHCI and the silicon surface, part of its potential energy is deposited in the form of electrostatic energy in a hemisphere with its equator lying on the upper plane of the sample. The atoms initially inside the hemisphere are singly charged. And the following energy balance relationship is solved to estimate the radius R of the hemisphere, which links the SHCI potential energy to the following MD simulations [22,23].
Where E p and E c are the potential energy of SHCI and the electrostatic energy of the uniform charged hemisphere, N q and I c are the number of Si + and the first ionization energy of Si atom, E el and E ph are the energies used to generate secondary electrons and photons during the interaction between SHCI and the surface, n and ε are the atomic density and dielectric constant of the solid, respectively. To solve the above equations, the experimental electron emission and X-ray yields of Refs. [24] are used to estimate the E el and E ph . The equations solving gives the radius R is about ~ 13Å for Xe 44+ and silicon surface. It should be noted that the kinetic energy of SHCI is not considered in present simulation for its little effect on the surface nano-structure formation. The Tersoff potential [25,26] combined with the repulsive ZBL potential [27] is used to describe the interatomic interactions. And the time step Δt=0.0001 fs is used in the present MD simulations. The simulations run sufficiently long time and terminate when there is no further change in the surface topography. Fig 1 shows the morphology changes of the surface induced by Xe 44+ ion impact, four snapshots of the surface are listed for t=0, 50, 100, 370 fs, respectively. As can be seen, the region filled with Si + ions (indicated by red spheres) expands significantly at t=50 fs, and a pit is formed at t=100 fs. At this time, the size of the pit is approximately equal to the initially charged hemisphere region. At t=370 fs, 119 Si + ions and 17 Si atoms are sputtered out from the surface since the extreme non-equilibrium of the system, leaving a crater-like nanostructure of ~40 Å in diameter and ~18Å in depth in the surface.

Results and Discussion
These results are consistent with the experimental values of Refs [2,24]. To study the time evolutions of the temperature, energy, pressure and density of the substrate during the nanostructure formation, we divided the system into six sub-regions as shown in Fig 1 (a), which are labeled by shell 1, shell 2, and so on, respectively. These sub-regions are three dimensional shells,  3 the center of which is at the center of the upper plane of the sample, and the thickness is 2.5 Å. Fig 2  depicts the temperatures of the sub-regions as functions of time. As seen from the figure, the temperature growth rate of each sub-region is ~ 10 18 K/s, and the temperature of each sub-region reached their maximum values one after the other in the initial 175 fs, then gradually return back to ambient temperature. The highest temperature corresponds to 2.8×10 5 K at t=50 fs for shell1. The above results indicate that a shock wave is formed in the substrate during this process. And the wave speed is about 10 4 m/s. The temperature of shell 1 is much higher than those of the other shells, even at t=350 fs. This may be a result of the detachment of the particles in this region from the substrate within 50 fs. At t=250fs, the temperatures of shells 2-6 converge to the same value, the shock wave dissipated rapidly in the substrate, which basically coincides with the formation time of the nano-pit on the surface.       Fig 5. shows the pressures of sub-regions as functions of time. At t=0 fs, the pressures of shells 1 and 2 are up to 10 5 atmospheres due to the large amounts of repulsive energy stored in the hemispherical charged region. However, the highest pressures quickly drop to ~5×10 4 atmospheres at t=150fs. After t=250 fs, the pressures in all shells seem to converge to the same value. At t=700 fs, the pressures of all shells are still very high, about 10 3 atmospheres, even though they are much less than their initial values. The pressures of shells 5 and 6 undergo a similar evolution, that is to say, the pressures increase first and then decrease, and hit their maximum values at 100 fs and 200 fs, respectively. Plot is for shells 4-6.
At t=0fs, the atomic number density of the substrate is uniform as expected, which is about 0.051 Å -3 .
With the occurrence of the Coulomb explosion [28,29] and the generation of shock waves, the surface undergoes a high temperature and high pressure process. During this process, some atoms are sputtered from the surface, and the atomic number density of each sub-region inside the surface changes. Fig 6. shows the atomic number densities of sub-regions as functions of time. Contrary to the formation of a high density region near shells 4 and 5, the densities of shells 1-3 decrease rapidly in the initial 50 fs. The highest density of shell 4 is ~ 0.064 Å -3 , which is more than 20% of the initial density and occurs at t=20 fs. At t=20~40 fs, the densities of shells 4-5 decrease in different extent. This non-equilibrium process indicates the hole formation due to Coulomb explosion, and accord with the time evolution of the surface morphology in Fig 1. In addition, different from the shells 1-5, the density of shell 6 has no obvious change throughout the simulation.

Conclusions
In this paper, based on analyzing the multiple electron emission of the substrate, we performed MD simulations study on the surface nanostructure formation resulting from the interaction of Xe 44+ ion with silicon surface. The Tersoff potential combined with the repulsive ZBL potential is used. The results show that a crater-like nanostructure with diameter of ~40 Å and depth of ~18 Å is formed on the surface. The time evolutions of the spatial distributions of temperature, energy, pressure and density of the substrate are investigated, respectively. The results show the surface undergoes a high temperature and high pressure process during the formation of the surface nanostructure. A shock wave at a speed of 10 4 m/s propagates in the substrate is formed in the initial 175 fs. The local temperature and pressure of the substrate are up to 2.8×10 5 K and 10 5 atmospheres, respectively. Many particles are ejected from the surface since the extreme non-equilibrium of the system during the nanostructure formation.