Simulations of silver-doped germanium-selenide glasses and their response to radiation

Chalcogenide glasses doped with silver have many applications including their use as a novel radiation sensor. In this paper, we undertake the first atomistic simulation of radiation damage and healing in silver-doped Germanium-selenide glass. We jointly employ empirical potentials and ab initio methods to create and characterize new structural models and to show that they are in accord with many experimental observations. Next, we simulate a thermal spike and track the evolution of the radiation damage and its eventual healing by application of a simulated annealing process. The silver network is strongly affected by the rearrangements, and its connectivity (and thus contribution to the electrical conductivity) change rapidly in time. The electronic structure of the material after annealing is essentially identical to that of the initial structure.


Background
Chalcogenide materials are among the most flexible and useful in current technology. Certain GeSbTe alloys are the basis of phase change memory technology [1] (now a credible alternative to conventional FLASH memory) and DVDs [2]. Amorphous Se is the active element for digital x-ray radiography [3], and metal-doped chalcogenide glasses are among the best known solid electrolytes or 'fast ion conductors' [4] and form the basis for another quite promising class of FLASH memory devices, 'conducting bridge' memory. The basic science of the material is just as interesting as the other phenomena such as the optomechanical effect [5] and photomelting [6]. Recently, a new application has emerged: the use of chalcogenide glasses for the detection or sensing of radiation (a dosimeter) [7,8]. The electrical conductivity is found to be well-correlated to radiation dose [9]. With annealing, the damage is readily reversed so that the device may be reused. This important discovery is presently empirically understood, suggesting the need for theoretical research both to understand the basic process and to aid in optimizing the materials for future device application.
In this paper, we undertake the first simulation to understand the atomistics of the response of chalcogenide glasses to highly energetic events. Like many other challenging material problems, we find it helpful to use multiple methods, in this case both empirical potentials, and ab initio techniques. We also have taken advantage of the contributions of others, such as the use of an appropriate 'heat bath' to handle the excess thermal energy after the thermal spike [10]. We detail the disordering process from a knock-on event to the subsequent recovery process. We show that the spike is indeed reversible upon annealing and discuss the electron states near the Fermi level -those responsible for the changes in (electronic) conduction after radiation exposure. The picture that emerges is that the electronic transport is significantly determined by the connectivity of the Ag subnetwork. Thus, if clusters percolate through the entire system, we have a network of nanowires that provide a low resistivity. As these nanowires break, form, or otherwise change, the carrier transport changes accordingly.
This paper is organized as follows. First, we discuss the formation of a suitable model of the material. Next, we track the damage wrought by a thermal spike http://www.nanoscalereslett.com/content/9/1/594 and its healing. The electronic properties are discussed briefly and conclusions are provided at the end of the paper.

Model formation
Simulations have been widely used to characterize amorphous materials and satisfactory ab initio models of GeSeAg glasses for various compositions have already been reported [11,12]. Molecular dynamics (MD) simulation is a natural approach to simulate high-energy processes because MD offers detailed trajectories of the atoms as the system evolves after the radiation induced event. The pitfall of MD is that it is only as good as the force field used and it is as computationally costly as it is detailed. Furthermore, large models are needed to realistically simulate radiation events in any material and this greatly increases the computational cost.
Many simulations of radiation damage have been presented, with varying details and system sizes ranging from 446 to 2.3 million atoms; the following is a highly incomplete list [10,[13][14][15]. Early simulations applied many approximations like the binary collision approximation (BCA) [16], linear interactions [16,17], and others to reduce the computational demand. Clever algorithms and parallel machines have enabled full simulations on large models [15]. We used the potential of Iyetomi et al. [18] to model the interatomic attractions. This potential is simple in its form containing a Coulomb interaction term, charge-dipole interaction term, and a short range repulsion term, and yet commendable in its ability to predict wide range of properties of this material.
The primary model of a silver rich glass (with stoichiometry close to that used for a detector) was fabricated as follows. We used a cubic supercell containing 5,184 atoms with periodic boundary condition to represent bulk Ge 3 Se 9 Ag 4 (the cell had 984, 2,888, and 1,312 atoms of Ge, Se, and Ag, respectively). The model described in this work is obtained by using the melt-quenching method [19]. Starting from a randomly placed collection of atoms, we performed 10 5 steps of MD with constant NVE. Then the atoms were given random velocities corresponding to a macroscopic temperature of 5,000 K and were allowed to equilibrate for another 10 5 steps. The system was then cooled to 1,200 K over 3.8 × 10 5 steps and equilibrated at 1,200 K for another 10 5 steps. The system was then cooled to 300 K over 0.75 × 10 5 steps and equilibrated at 300 K for 10 5 steps. Finally, the system was relaxed using a conjugate gradient algorithm. The MD simulations described in this work were performed using the classical molecular dynamics simulation package LAMMPS [20]. A time step of 1 fs was used throughout, except when variable time steps were required.
This model faithfully reproduces many features of the material (see also reference [18]). The total radial distribution function (RDF) given by this model is in reasonable agreement with the experimental RDF [21] including a two-peak first neighbor feature in g(r) (see Figure 1). Our model also reproduces the experimental nearest neighbor distances of Ge-Se, Se-Se, and Se-Ag correlations [21]. It overestimates the Ag-Ag correlation distance, but this can be understood considering the broad peak of Ag-Ag pair correlation function (see Figure 2). The first peak of the partial RDF of Se-Ag displays a companion alongside the first peak in total RDF. The Se-Ag interatomic potential has a shallow minimum near the Se-Ag correlation distance, whereas all other interactions of Ag-atoms are repulsive. As a result, Ag atoms are very mobile in the network and Se atoms see high coordination with respect to Ag. The exceptional mobility of Ag atoms with respect to the host atoms is a widely reported phenomenon [11,18,22], and the basis of many applications, and even the accelerated crystallization of phase-change memory materials [23]. The mean-squared displacement (MSD) of atoms calculated for our model also predicts the high mobility of silver at all temperatures below the melting point of the material. Figure 3 shows the MSD for Ag atoms at different temperatures. The diffusion coefficients and conductivity calculated using the total MSD of the system compare favorably with the corresponding experimental and ab initio values reported for most similar system as ours (see Table 1).

Damage simulation using thermal spike
We carry out a thermal spike simulation using the model above. The damage inflicted on a material by high-energy  radiation starts with a sudden transfer of kinetic energy from the incoming particle to an atom or a group of atoms that happen to suffer a collision with the incoming particle. For incident particles of energy in the range of MeV, the first interaction with the atoms on the target is entirely ballistic and the detailed role of interatomic potential between the impinging projectile and the target can be neglected. So, following Rubia et al. [14], we modeled the onset of radiation damage by igniting a thermal spike at the center of the supercell. We defined a sphere of radius 2.5 Å located at the center of the supercell, i.e., at d/2, d/2, and d/2, where d is the size of the cubic supercell to receive the thermal spike. For the particular configuration we modeled, this sphere contained two atoms. These two atoms were given an initial velocity consistent with 1 MeV of energy, and the rest of the atoms were assigned a veloc-ity distribution associated with a temperature of 300 K. These conditions are thus intended to mimic a damage event at 300 K with the center two atoms representing the primary knock-on atoms (PKAs). The system was then allowed to evolve. For an isotropic material like GeSeAg and energy of PKA being as high as MeV, the direction of the initial velocity of the PKA should have no observable effect on the damage production.
In view of the large velocity imparted to the atoms and the unknown behavior of the empirical potential under extreme conditions, we used a variable time step as in [25]. To avoid the diverging cascade of damage from bouncing back from the boundary, we used a damped outer layer of 0.5 Å thickness. The velocities of the atoms falling in this boundary region were rescaled at every dynamical step [26]. The schematic diagram of the simulation setup is shown in Figure 4.

Results and discussion
Immediately after the detonation of a thermal spike at the center, the hot atoms' trajectories resemble a projectile, and they pass through the network transferring huge amount of momentum to the atoms along their trajectory. These recipients of momentum get knocked from the network and, in turn, start moving like a projectile, thus , where the velocity rescaling was applied (the outer boundary), and where the normal molecular dynamics was performed. http://www.nanoscalereslett.com/content/9/1/594 creating a cascade. At this stage, the velocity of the PKA is highly non-canonical as indicated by the observation that temperature of the system excluding the PKA lags behind the temperature of the whole system for a brief period of time (see Figure 5). An animation tracking the positions of atoms and the average temperature of the system reveals rapid local melting. The average temperature of the system remains above the melting point of GeSe 3 Ag (approximately 1150 K, [21]) for 11.9 ps. To visualize radiation damage inflicted on the network, we present snapshots of a slice of the supercell containing central 6 Å of the simulation box at different times after the detonation (Figure 6). The damage is most conspicuous at 2.25 ps, and the image clearly shows the formation of voids and internal surfaces. At 10 ps, the system can be seen returning to its initial structure as also indicated by first peak of total RDF gaining height around that time.
With the onset of local melting, the system loses its short-range order. The temporal evolution of the RDF in Figure 7 shows an interesting recovery. The short-range order can be seen evolving continuously with the first peak gaining height and the first minimum continuously deepening. The evolution of second peak follows. The hump in the first peak originating from Se-Ag coordination, however, does not reappear as late as 50 ps.
To investigate the effect of radiation events on the phase separation of Ag atoms, we performed a cluster analysis of Ag atoms over the entire evolution of damage. We defined a cluster as a group of all atoms lying within a cutoff distance from at least one other atom belonging to the same cluster. We chose a cutoff distance of 3.7 Å (the Ag-Ag correlation distance is 3.55 Å) to define the cluster. Defined in this way, the initial configuration has one large cluster of about 750 atoms and other numerous smaller clusters. At just 2 ps into the damage evolution, the clustering of Ag is lost and the cluster distribution at this point essentially resembles a random configuration. The clusters, however, begin to grow as the network rearranges and recovers its initial connectivity (see Figure 8). The size of the biggest cluster in the network hits its minimum between 1.6 and 2.6 ps (see Figure 5). This is an interesting observation in that it lags behind the time the system has the highest temperature. We observe that the electrical conductivity presumably depends sensitively on the connectivity of the Ag subnetwork, the number, and the structure of the Ag filaments.
We performed ab initio calculation of the electronic densities of states (EDOS) of a 648-atom model prepared using the same empirical potential and the damaged snapshots of this model. The size of this model is a compromise between being large enough for damage production and being small enough for an ab initio calculation. Our calculations and reference [18] have confirmed that the 648-atom model is statistically similar to  the 5184-atom model we discussed above. EDOS calculation is done using plane-wave basis code VASP [27][28][29]. Plane waves of up to 350 eV and PAW potentials were used [30,31].
EDOS and inverse participation ratio (IPR) (a measure of the spatial localization of the Kohn-Sham electronic states) of our 648-atom model is plotted in Figure 9. A comparison with EDOS from an ab initio model [11] reveals that our model lacks an energy gap though there are fewer slightly localized states in the expected gap region. Tracing the structural origin of these states, we note these are p-orbitals localized around Se atoms that are under coordinated with respect to Ge atoms (either one-fold or no Ge coordination). Our calculation on a smaller 100-atom model of [11], prepared using first principles method, also shows that gap states are contributed by Se atoms that are under-coordinated with respect to Ge atoms.
We also analyzed the evolution of the electronic structure as the system passes through damage and restoration (see Figure 10). The damaged structures lead to gap states arising from under-coordinated Se atoms. It is interesting to note that the electronic structure recovers very quickly compared to the network itself. The RDF evolution of Figure 7 suggests that although the first peak and first minimum begin to take qualitative shape very early, the network attains similar structural order as the starting configuration as late as 25 ps. The electronic states however show a remarkable healing and reversibility even at 10 ps after the incidence of radiation. This naturally points Figure 9 Density of electronic states and Localization. The electronic density of states and the IPR, which gauges spatial localization of states of 648-atom models prepared using empirical potential of reference [18]. Fermi level is at 2.01 eV. http://www.nanoscalereslett.com/content/9/1/594

Figure 10
The evolution and reversibility of electronic structure. Electronic density of states of six instantaneous configurations of the model at different times with respect to damage event. The model structure and damaged structures were produced using empirical potential, and the electronic structure was calculated for these structures using first principle methods. Note the high degree of reversibility (comparing t = 0 and t = 50 ps).
to an interesting physics underlying the damage recovery that some features of network recover early whereas some others (eg., the Se-Ag correlation) are not recovered even after some equilibration of the system. Prasai and Drabold [11] have pointed that Se atom bonded with Ag atom contributes to widening the gap. Our work predicts that Se-Ag correlation is lost in the network as a result of radiation-induced damage. The details of the physics of these two separate observations are not fully understood yet.
We have not yet modeled the conductivity of the system. This is complicated by the existence of both ionic and electronic conductivities for some parts of the simulation. Realistic calculations are currently being formulated [32].

Conclusions
This paper is the first word on the atomistic processes underlying the fascinating experiments and device of Mitkova and coworkers [7,9]. To fully realize the potential of our approach, many issues such as cell size, composition, details of the modeling of thermal spikes (and subsequent relaxation), and material composition must be explored. Observables like the transport and optical properties should be extracted at representative moments in the simulation. Nevertheless, this work reveals key aspects seen in the experiments including a remarkable reversibility upon annealing. We show that judicious use of the empirical potential of Iyetomi et al. [18] leads to a credible model of the dynamical processes and correctly reproduces many aspects of the material. Ab initio methods are an important tool to augment this work and to understand its limitations and the electronic and optical properties.