Excitation and Relaxation of Swift Heavy Ion Irradiated Dielectrics (cid:3)

The structural modiﬁcations of insulators after swift heavy ion irradiation are studied by applying a combination of the Monte-Carlo method (MC), used to describe the excitation and relaxation induced by the penetrating ion in the electronic subsystem, with the two temperature model (TTM) describing the temporal evolution of the lattice temperature heated by the hot electronic system. This MC-TTM combination demonstrates that secondary ionizations play a very important role for the track formation process and is capable to compute important material parameters which are often unaccessible by experiments. It is found that the secondary electron generation leads to an additional energy source in the heat diﬀusion equation related to energy storage in holes. [DOI: 10.1380/ejssnt.2010.278]


I. INTRODUCTION
When an ion passes through an insulator, it loses its energy as a result of the interaction with the solid. This energy loss is due to two main dissipation channels: elastic collisions with the target atoms (nuclear stopping) and the excitation of the target electrons (electronic stopping). For ions with masses of more than 20 proton masses and kinetic energies in the MeV regime the energy loss due to electronic excitation dominates by about two orders of magnitude [1,2]. Such swift heavy ions (SHI) are known to induce nanoscaled material modifications in insulators, which can be observed using an atomic force microscope for instance [3].
This phenomenon of track creation is often explained in terms of a two temperature model (TTM) commonly called the inelastic thermal spike model [4]. Within this model both the electronic as well as the phononic system are described via two heat diffusion equations coupled by an exchange parameter. While the primary energy loss of the ion is used as a source of energy for the electronic system, the phononic system is heated indirectly by electronphonon coupling. This leads to a local heating of the lattice which may result in a molten area. The thermal spike * This paper was presented at 10th International Conference on Atomically Controlled Surfaces, Interfaces and Nanostructures (ACSIN-10), Granada Conference Centre, Spain, 21-25 September, 2009. † Corresponding author: orkhan.osmani@uni-due.de model has been successfully applied for reproducing experimentally observed track radii and damage thresholds [3][4][5][6].
However, in order to explain the experimentally observed effects of the SHI irradiation, the parameters used in the TTM approach, like the electronic diffusivity and the electron-phonon coupling parameter, need to be handled with care. These parameters are often not known. The difficulty is due to the variable electron density in excited dielectrics, thus these parameters depend strongly on transient excitation dynamics. A reliable theory considering complete excitation dynamics is needed to compute these quantities. Additionally, the TTM is based upon thermodynamical equations, which are not applicable because the system is not in equilibrium after the ion impact. The electron-electron collision time is in the femtosecond regime, and the TTM may be applied after about 100 fs [7] to describe the energy exchange between electrons and phonons which occurs on the timescale of a few picoseconds [8,9]. The difference between these timescales requires special attention when choosing the appropriate theoretical model.
In order to explain the experimental observations, previous applications of the TTM met another problem: it appeared necessary to assume an increased electronphonon coupling rate or a decreased electronic diffusivity to confine the energy in the time-and space domain to achieve the melting temperature in the track core [4,5]. However, the physical origin of this energy confinement has not been further studied. In order to follow the dynamics of SHI irradiated insulators we have performed The Surface Science Society of Japan (http://www.sssj.org/ejssnt) classical Monte-Carlo simulations (MC) which are capable to describe the excitation and transport of the electrons even at highly non-thermal conditions. The output of these MC calculations is then used as initial conditions for a TTM run to compute the lattice heating. This model reveals the physical origin of the energy confinement in the track core as an energy storage in holes and followed by the release of the energy via Auger processes and electron-hole recombinations.

II. MODEL
The Monte-Carlo simulation is one of the classical ways to describe multi-event phenomena, especially for systems of particles not in thermal equilibrium [10][11][12]. It has often been used to describe the processes of excitation and relaxation of electrons irradiated with SHI or laser pulses [13][14][15][16][17]. In the framework of the model every traced particle undergoes a collision process with a certain probability which is randomly chosen. Each such process is simulated event by event. The probability of a particular event depends on its cross-section relative to the total cross-section of possible interactions.
We apply the same Monte-Carlo approach as described in detail in [15][16][17]. The general idea of our developed algorithm is as follows: first, we calculate the free path length between collisions of the ion with the target atoms; then multiple ionization occurs for each collision which results in the creation of the first generation of free electrons. Second, we determine the path lengths between all possible subsequent collisions for all free electrons with the other bound electrons in different energetic states; then the bound electron corresponding to the shortest possible path length is ionized and the transferred energy is calculated. Consequently the mean free path and time for each electron, i.e. for electrons with different energies are also obtained. Each ionization of a bound electron leads to the creation of a hole. The possibility of hole decay by Auger recombination and thus further creation of secondary electrons is also included in the MC simulation, which are repeated many times and finally averaged. The importance of electron-hole pairs has been pointed out already by other authors [30,31]. We will see later that the hole decay indeed plays a crucial role for the energy confinement mentioned above.
For simplicity, we assume perpendicular incidence and neglect any nuclear stopping. Therefore we use a cylindrical geometry with periodic boundary conditions. The projectile is considered as a point-like particle with an equilibrium charge state according to the Barkas formula [2,14,15]. SHI penetration through homogeneous media is described with a Poisson law for the free path between a series of collisions. For solids this mean free path can be chosen equal to the mean interatomic distance [11][12][13][14]. During the collision of SHI with an atom, the electrons of the atom are placed randomly around the nucleus according to their energy levels. The impact parameter between the projectile and the atomic nucleus is chosen randomly within the interatomic distance [12]. According to this impact parameter, the energy transferred to every electron is calculated [15]. The ionized electrons are treated as independent particles, i.e. the energy transfer and an-gles of emission of electrons are uncorrelated [11][12][13][14]. If this energy exceeds the ionization potential of the respective electron, this electron is ionized, otherwise, no energy transfer takes place and no free electron is created. The scattering angle is explicitly determined by the transferred energy, and the polar angle is randomly uniformly distributed within the interval [0, 2π).
For the subsequent electronic propagation, in contrast to SHI, elastic collisions with atoms must be included. Therefore, for electrons there are two channels of energy losses: elastic and inelastic ones. The free path of an electron is defined by the shortest path realized according to Poisson's laws for every possible collision: elastic collisions with any kind of target atoms, and inelastic collisions with bound electrons at different energy levels of every type of atoms.
Mean free paths of elastic collisions are calculated with Mott's cross-section of electronic scattering [11], which depends on the energy of the electron and the atomic number of the colliding atom. The scattering angle is defined by at transferred energy and the polar angle is chosen from the interval [0, 2π). For low energy electrons there is also another possibility for energy loss related to scattering on phonons [10,11]. The electron-phonon scattering process is later incorporated within the TTM.
The mean free paths of the inelastic collisions are calculated using Gryzinski's ionization cross-section [14,19], which depends only on the energy of the respective electron and the ionization potential of the bound electron. The secondary generated electrons can then travel and scatter in the same manner [15][16][17].
The key point of this work is how Auger transitions are taken into account which occur after the ion has ionized the target atoms. To calculate the temporal decay of holes via Auger-recombination the Poisson law is applied [17]. Each atomic shell of each type of atom has a characteristic time [20]. After the calculation of the decay time, the electron gaining the excess energy is chosen. Its final energy is given by the difference between the energy released by the filling of the hole in a deeper shell and its own ionization potential. Its corresponding momentum is chosen uniformly within the solid angle.
For solid targets interatomic Auger-processes, so-called Knotek-Feibelman processes [21][22][23] have to be taken into account as well. In these processes the hole is filled by an electron coming from a neighboring atom. Multiple ionization occurs frequently during SHI irradiation leading to a local depletion of electrons which enhances the occurrence of Knotek-Feibelmann transitions. The characteristic time of these processes are assumed to be equal to usual Auger recombination times. Detailed analysis of the kinetics of relaxational processes of holes due to Auger and Knotek-Feibelman decays is out of the scope of the present work and can be found in [18].
Following every electron and then averaging over the ensemble we determine the spatial and temporal distributions of these particles and their energies. Once the electrons show diffusive (hence thermal) instead of ballistic behavior, we consider the MC part to be finished. Then we use the obtained energy distributions as initial conditions for TTM calculation.
Within this model the temporal and spatial evolution of the electronic and lattice temperature, T e and T p , respectively, are calculated using the TTM given as: Here, D e and D p denote the difussivity of electrons and lattice, respectively, g is the coupling parameter between both subsystems and S(⃗ r, t) is the space-and-time dependend energy source heating the electronic subsystem. For the numerical description the crystal is discretized and equations (1) and (2) are solved fully three dimensional using the finite difference method. During the simulation the system is in contact with an infinite large heat bath at room temperature (assuming the experiment is performed at room temperature) serving as boundary conditions. However electrons with higher energy than the work function of the solid can escape, removing energy from the system. To take that into account we apply open boundary conditions for the electronic system at the surface allowing energetic electrons to be emitted into the vacuum. This set of heat diffusion equations is coupled via the electron-phonon coupling parameter g, which determines the amount of energy transferred from the electrons to the phonons per unit time and volume. It is evident that this parameter plays a key role, as it governs the heating of the lattice. However, this parameter is difficult to come by, as only limited experimental data is available. Also the electronic diffusivity D e poses a problem. It determines the amount of energy that is transported away from the point of release (i.g. the position were the electron was ionized). This transport will effectively limit the amount of energy that is locally available in order to heat the lattice above the melting point. Experimental data on this property is lacking. Although one can extract this number from resistivity measurements for instance, such kind of measurements often do not reflect the transient electronic excitation during the irradiation.
The electronic system is heated while the passing ion loses its energy. This can be caluclated using the SRIM code [27]. The transferred energy enters eq. (1) in the source term S and in the present notation is formulated in terms of a temperature increase. It is therefore necessary to convert the excitation energy into an electronic temperature. This can be done by assuming that the electrons behave like a diluted free electron gas and using that the internal energy of the electrons is related to the temperature by E = γT 2 /2 where γT is the specific heat capacity of the free electron gas: n e k B T F with n e the free electron gas density, m e is the electron mass, k B the Boltzmann constant and T F the Fermi temperature.
The electronic diffusivity is a function of the temperature (see results and discussion section) which is a function of the coordinate: D e (T e (⃗ r)). Therefore the gradient of the heat flux occurs in Eq. (1). In contrast, the lattice diffusivity can be regarded as constant, thus the corresponding term in Eq.
(2) appears in a simplified form. The lattice diffusivity is given as D p = C p /K p ; we assume a linear temperature dependence of the heat conductivity K p [28] and of the specific heat capacity C p . These parameters are taken from [29] while all other parameters entering eqs. (1) and (2) are obtained as an output of the MC part of the model, see below.
As mentioned above, it is important to take into account the energy redistribution related to holes. The holes may decay via the Auger recombination. In this process, the recombination energy is released to an electron. Thus, these recombinations act as an additional source of energy for the electronic system, which is governed by a certain characteristic recombination time. Assuming an exponential decay of the holes we can write this source term as: Here E h (⃗ r) denotes the initial excess energy of the hole at position ⃗ r and τ is the characteristic decay time. Therefore the effective source-term in eq. (1) heating the electronic subsystem can be written as S(⃗ r, t) = S h (⃗ r, t) + S SHI (⃗ r, t) where S SHI (⃗ r, t) is the energy introduced to the electronic system by the primary ion.

III. RESULTS AND DISCUSSION
As a model system we choose the irradiation of solid SiO 2 by a Ca 19+ ions with a total energy of 11.4 MeV/u which corresponds to a recent experiment [32]. By means of Monte-Carlo simulation, within this present work, it was found that only after approximately 1 ps after the ion impact the electronic behavior can be treated as a diffusive, thus the electronic system is assumed to be thermalized (to be published). Therefore a description based upon thermal equations is only valid at times around 1ps and is considered as the starting time for the TTM.
In our simulation we find that the total energy of holes amounts to ≈ 82% of the energy transferred by the SHI. According to eq. (3), these holes act as an additional source of energy for the electronic system. Due to the fact that a major fraction of the transferred energy is released with a delay the heating of the lattice is prolonged to picosecond timescales. We also observed that the radial dependence of the electronic distribution plays a key role in the description of SHI induced track creation. High energetic electrons are escaping from the central track region creating slower electrons in their wake. As a consequence the electronic temperature is spatial dependent, and the time of equilibration for electronic and lattice temperatures will differ for different radii.
The benefit of our combined approach is that parameters needed for the TTM can be calculated within the MC approach. Here we will demonstrate this for the electron phonon coupling parameter g and for the electronic diffusivity D e . As was mentioned before the mean free path λ and the velocity v(E) of the electrons with the energy E are obtained within the MC part of this calculation. Using the well known expression D = vλ/3 we obtain an energy dependent electronic diffusivity D e (E) for SiO 2 . This energy dependence can be transformed into a temperature dependence as was stated in the previous section.
The electronic diffusivity D e (E) shown in Fig. 1 exhibits a linear energy dependence for small energies. This behavior is also observed in [5] where only electronelectron scattering was taken into account. However, within the MC approach we also consider the scattering of electrons with the target atoms. By including this process we find that the diffusivity for electronic energies around 20 eV is at its minimum while for larger energies D e (E) is increasing again. Electrons with energies around 20 eV are located close to the track core while less energetic electrons are found further away, hence electrons close to the track core will have a lower diffusivity compared to those created further away.
We will now evaluate the electron phonon coupling parameter g. The coupling parameter describes the amount of energy transferred from the electrons to the lattice per unit time and volume. From the MC calculation we know the individual time between collisions of electrons and target atoms. The transferred energy between electrons and atoms is calculated within the binary collisions approximation. Because particle trajectories are calculated, the volume in which this energy is transferred is the volume of one atom. Thus we obtain for SiO 2 an electron phonon coupling constant g ≈ 10 18 J·(s m 3 K) −1 . This is comparable to coupling parameters that were obtained experimentally for metals [26,34,35] and is reasonable also for laser-excited dielectrics [8].
The energy transfer from the electrons to the atoms is shown in Fig. 2. The scattering of the data at higher electronic energies is due to lack of statistics at these energies.
With the electron diffusivity D e (T e ) and the electronphonon coupling parameter g all crucial material parameters entering the TTM are thus extracted out of the MC part. The TTM calculations are carried out until the electronic and lattice temperatures are equilibrated. Within this period of time the lattice temperature, if sufficient energy is provided by the electrons, reaches the melting temperature and the material can be assumed to be locally molten. For our model system the size of this modified region is about 4 nm and is in good agreement with the experimental findings.
Additionally, the introduced source term S h related to the holes combined with the temperature-depended electronic diffusivity D e (T e ) confines the electronic energy at the place of excitation. It should be noted here that this is not a confinement in the sense of a stationary electron gas as was presumed in previous works [1,5], but rather a repeated energy emission into the electronic system by Auger recombination of holes at certain positions in the crystal. The reproducibility of the experimental track radius and the successful computation of necessary parameters for the TTM, which have been formerly treated as fitting parameters proof the versatility of the developed MC-TTM model for the discription of swift heavy ion induced processes in irradiated solids.

IV. CONCLUSION
We developed a combination of the Monte-Carlo method with a two temperature model (MC-TTM) to describe the track creation processes in dielectric targets after swift heavy ion irradiation. Within the MC part of the calculations all necessary (and experimentally inaccessible) material parameters are computed so that no fitting is needed within the TTM. Especially the electron phonon coupling parameter is extracted which is an important material dependent parameter. The results of the combined MC-TTM model show that the energy storage in holes and secondary electron creation by the first generation of ionized electrons as well as by Auger decay of holes has to be taken into account. The Auger decay of holes plays a key role as an energy source of the electronic system, acting as a delay and thus confining the energy at the place of its initial excitation, i.e. in the vicinity of the track core. It can be accounted for by means of an additional source term in the heat diffusion equation. The obtained structural modification radius of the material calculated with MC-TTM is in very good agreement with experimental data.
Because universal cross sections enter the MC calculations, this kind of computation can be carried out for any given ion-target combination and thus the new method presented here will prove helpful in predicting ion track properties for any other dielectric material as well.