Development of Continuum-Atomistic Approach for Modeling Metal Irradiation by Heavy Ions

Over the last several decades active research in the field of materials irradiation by high-energy heavy ions has been worked out. The experiments in this area are labor-consuming and expensive. Therefore the improvement of the existing mathematical models and the development of new ones based on the experimental data of interaction of high-energy heavy ions with materials are of interest. Presently, two approaches are used for studying these processes: a thermal spike model and molecular dynamics methods. The combination of these two approaches – the continuous-atomistic model – will give the opportunity to investigate more thoroughly the processes of irradiation of materials by high-energy heavy ions. To solve the equations of the continuous-atomistic model, a software package was developed and the block of molecular dynamics software was tested on the heterogeneous cluster HybriLIT.


Introduction
The continuous-atomistic approach for modeling the interaction of high-energy heavy ions (HEHI) with various materials combines the equation of thermal conductivity of the electron gas of the thermal spike (TS) model [1,2] and the equations of motion of atoms of the molecular dynamics (MD) method [3].This combination, in contrast to the MD method, allows one to consider processes over a wide range of energy of the incident heavy ion to the target (since inelastic energy losses predominate at high energies of the incident ion and the MD method does not work), and unlike the TS model, it allows to obtain more detailed information on the simulated system.A continuous-atomistic model (CAM) [4] was proposed in 2003 to describe the effect of short-pulsed lasers on metal targets.The use of combined methods requires the development of new numerical schemes, algorithms and program packages.Since the combined models consist of equations of two different classes, special attention is to be paid to the consistency of the numerical methods for both classes of equations.Modeling within the framework of the combined approach requires considerable computational resources, and thus it is necessary to develop new parallel algorithms for multiprocessor computer systems.To carry out the simulation, a software package was developed with the possibility of using it on high-performance computing systems [5].e-mail: zafar@jinr.ru The aim of the present work is to simulate the irradiation of metal targets by HEHI in the framework of the CAM, in terms of the parameters of the source function and the electron-phonon interaction coefficient, as well as the development and refinement of the capabilities of the CAM to describe the available experimental data.

Formulation of the problem
The continuous-atomistic model [4], which describes the processes in a metal film under the action of a short-pulse laser, is taken as basis of the CAM.The model described below assumes a continuous equation for the thermal conductivity of the electronic subsystem of the TS model in Cartesian coordinates and equations of motion of the atoms within the MD method: where A detailed description of equation ( 1) is given in [2], and of equations 2 -in [4].
In this paper we consider, as an example, the irradiation of nickel by uranium ions of energy of 700 MeV for different values of the parameter r 0 of the source function A e (x, y, t) and of the electronphonon interaction coefficient g(T e ).The choice of the target and the incident ion is due to the fact that the authors of [2] performed a simulation of the same ion-target combination within the TS model.The use of the two-dimensional CAM is motivated by the fact that the processes are being modeled in a thin (1-10 nm thick) material.

Numerical methods
To solve the equation ( 1), the computational domain is discretized by a uniform grid over which finite-differences are used [6].The Linked-List Cell MD Algorithm [7] is chosen to integrate the equations (2) of the particle motion.An important step is to calculate the temperature field of the irradiated material.To this aim, the computational domain is divided into superimposed cells so that in each cell the number of particles is 100-1000, and the temperature in each cell is calculated.In contradistinction to [4], this approach makes possible the increase of the density of points of the temperature field of the target.
In the computational schemes for solving equations ( 1)-( 2) it is necessary to preserve: 1) the consistency of the steps in time; 2) the temperature relationship at the node points of the finite-difference method and in the cells of the calculating region of the MD method.
To solve the system (1)-( 2), a software package TS-MD was implemented on a multiprocessor system with shared memory within the framework of the OpenMP standard.Results characterising the efficiency of the MD part of the TS-MD software are given in [5].

Results and discussion
To solve the equations of the CAM, a software package was developed the block of MD software of which was tested on the heterogeneous cluster HybriLIT [8].
The simulation of the process of interaction of uranium ions with a nickel target was carried out as follows.A nickel sample of size 35 × 35 × 1.1 nm with 140700 particles was selected as a target.The physical parameters of equation ( 1) for nickel are taken from [2].In the MD method, the embedded atom model (EAM) potential is used as the interparticle potential.The most important and interesting result in the interaction of HIHE with solids is the formation of spatial fractures of the crystal latticetracks [2].We have investigated this process in two cases.In the first case simulation was made for two values of the parameter r 0 (0.5 nm, 1 nm) of the source function and for the electron-phonon interaction coefficient g = g Ni .In Fig. 1 a) r 0 = 0.5, b) r 0 = 1 the results are given respectively at instants t 1 = 100 fs (left) and t 2 = 300 fs (right).The upper figures are top views and the lower figures are views at an angle of 45 degrees.In these figures, the square is the area of destruction of the crystal lattice.It can be seen that at smaller values of the parameter r 0 , the crystal lattice is severely damaged in a small region.For larger values the damage area increases, but the nature of the fractures does not change significantly.These results show the dependence of the model on the parameter r 0 .When this parameter is changed, the density of energy distribution over space changes.The parameter r 0 of the source function for HIHE is estimated to be within 1-2.5 nm [2].In the second stage, we change the electron-phonon interaction coefficient: g 1 = 0.5g Ni , g 2 = 2g Ni .The results obtained are shown in Fig. 2 a), b), respectively, at time instants t 1 = 100 fs and t 2 = 300 fs.The results show that the increase of the electron-phonon interaction coefficient g results in the increase of the energy transfer to the lattice.

Conclusion
Results of modeling the interaction of a flux of uranium ions of energy 700 MeV with a nickel target in terms of the parameters of the source function and the electron-phonon interaction coefficient are obtained using the developed software package of the CAM.The results of modeling, let us infer the following conclusions: 1. Changing the parameter r 0 in the source function strongly affects the results through the energy dependence of the incident ion and the type of ion.This factor becomes significant when the ion penetrates in depth the target since the moving ion loses more than 70-80 % of its initial energy in the first half of the path.
2. The choice of the coefficient of electron-phonon interaction is determined by straightforward comparison with the experimental data, since it corresponds to the amount of transferred energy from electrons to the lattice.

Figure 1 .
Figure 1.The dynamics of the destruction of the crystal lattice for r 0 = 0.5 nm (a) and r 0 = 1 nm (b) at time instants t 1 = 100 fs (left) and t 2 = 300 fs (right).The upper pictures are top views and the lower pictures are views at an angle of 45 degrees.

Figure 2 .
Figure 2. Dynamics of the destruction of the crystal lattice for parameters g = 0.5g Ni (a) and g = 2g Ni (b) at instants t 1 = 100 fs (left) and t 2 = 300 fs (right).The upper pictures are top views and the bottom pictures are views at an angle of 45 degrees.