Generation of Frenkel defects above the Debye temperature by proliferation of phonons near the Brillouin zone edge

A novel, non-radiative mechanism is reported by which Frenkel pairs of vacancies and interstitials are generated in molar concentrations far above thermal equilibrium. This mechanism is demonstrated in molecular dynamics (MD) simulations of an aluminum single crystal with a free surface. They suggest that three conditions must be fulfilled: (i) lattice vibrations near the Brillouin zone edge are being excited, (ii) these vibrations proliferate at a sufficiently high rate, and (iii) the sample temperature is above the Debye temperature (but significantly below the melting point). The simulations employed an EAM potential for Al. We attempt to draw a confluence between our MD simulations and recent experiments on flash sintering of aluminum. The simulation results are also consistent with flash experiments on polycrystals and single crystals of zirconium and titanium oxides where the Debye temperature was discovered to be the lower limit for the onset of the flash.


Introduction
Point defects, like vacancies or interstitials, play an important role in crystals. Their entropy imparts thermodynamic stability with Arrhenius temperature dependence related to the energy of formation. Usually, vacancies have a lower energy barrier and their concentrations exist in meaningful numbers at high homologous temperatures, enabling solid-state diffusion, which is the cornerstone of several high temperature phenomena in materials science. Interstitials however have a higher energy of formation, and in most instances, their concentration is too low to influence solid-state diffusion. A higher order defect configuration is the Frenkel pair where the interstitial and the vacancy are created simultaneously in a closed system by the displacement of an atom into an adjacent interstitial site. This event is energetically difficult, and the entropy it generates is lower than for single point defects; thus, the equilibrium concentrations of Frenkel pairs are scanty for most metals [1,2]. Frenkels are usually created by forced events such as radiation damage [3,4].
In this paper, we will present a new mechanism that generates lattice defects by the proliferation of short wavelength phonons at a temperature between the Debye and the melting temperature. In contrast to the thermodynamically equilibrated defects described above, this mechanism rests on a non-equilibrium process that is related to the excitation rate of the phonons. It has implications in processes such as sintering and high temperature deformation and fracture that rely on defects to induce solid-state diffusion. We show, through molecular dynamics (MD) simulations, that such defects, akin to Frenkel pairs, can be generated at concentrations that are far above equilibrium at relatively low temperatures without melting the crystal.
In the following sections the MD calculations are presented which lead to three necessary conditions for the dynamic generation of these 'Phonon Frenkels'. They are (i) that the specimen temperature is above the Debye temperature, (ii) the phonons proliferate at a high rate, and (iii) that the phonon wavenumbers are at the edge of the Brillouin zone. The mechanism has important implications in physical phenomena such as sintering and deformation in that they can be induced at relatively low temperatures. Recent work on flash sintering of aluminum [5] supports this conclusion.

Simulation model
We choose single crystal aluminum for the MD simulations since the results could be relevant to recent work where high heating rates have been shown to induce ultrafast sintering of aluminum powders [5]. The hypothesis is that defect generation by short wavelength phonons, above the Debye temperature could be the underlying reason.
MD simulates the movement of atoms; thus, the structure of the simulated crystal is known at all times. In this way atomic mechanisms, e.g. the generation of crystal defects, can be probed. Such information can then be used to understand the fundamental processes for the generation of crystal defects. To our knowledge, the work presented here is the first attempt to study the effects of proliferation of phonons of specific wavelengths in MD simulations.
The MD simulations are carried out with the well established program LAMMPS [6]. To simulate an aluminum crystal, we choose a well-tested version of the EAM potential [7,8] for Al crystals [9]. The initial configuration of our system is shown in figure 1. The Al atoms are placed on a perfect face centered cubic (fcc) lattice with a lattice constant for aluminum a 4.05 = Å [10]. The total system size is a a a 20 17 20´. Altogether, a block of a a a 20 2.5 20´atoms are set to be immobile. We choose periodic boundary conditions (PBC) along the aand c-axis, leaving a free surface on the 010 ( )-plane. The free surface was chosen in order to allow the necessary volume changes for the generation of Frenkel-interstitials. Note that, alternatively, such necessary volume changes could also have been achieved by implementing full PBCs and using a barostat. We expect, that this alternative setup would not change the basic results discussed in the next chapter. But a benefit of our chosen initial configuration is, that it allows us to investigate the generation of surface defects such as adatom-hole pairs along with the bulk defects.
The integration time-step is set to t 0.65 fs D = . This time-step was chosen in order to simulate the oscillations of phonon modes, without encumbering undue error.
The goal of our MD study is to investigate how non-equilibrium proliferation of lattice vibrations, characterized by the wave vector q and mode number j, modifies the crystal structure. This is achieved by carrying out the simulations in the following steps: at first, the initial system is equilibrated in the canonical ensemble at a certain temperature T 0 . The temperature is controlled by means of the Nosé-Hoover chain thermostat [11]. The simulations were carried out at values of T 0 , below and above the Debye temperature in order to determine its role in defect generation. At this stage the phonon distribution in the crystal corresponds to the temperature T 0 which is established in t 0.2 ns = . After this initialization of the system in thermal equilibrium, we drive the system out of equilibrium and explore how additional phonon modes of different wavelengths, excited at different rates, influence defect generation. These driven MD simulations are carried out by switching the thermostat off and extending the velocity-Verlet algorithm [12] to include the incremental velocity representing the additional phonons. This velocity, specified by V D , is calculated by the standard lattice dynamics theory [13][14][15], as described below. We begin by describing our system as a Bravais lattice where the initial position of the mth primitive unit cell is R m0 . A primitive unit cell contains r basis atoms which are labeled as μ=1, K , r. In our specific case study, there is one Al basis atom, thus r=1. By this description, we can find (in the harmonic approach) an analytical solution for a displacement u m j m ( ) of an atom μ in unit cell m as a function of the parameters, that characterize a certain lattice vibration, i.e. the wave vector q and mode number j r 1, , 3 = ¼ . As it is known from standard lattice dynamics theory [13][14][15], such displacements describe the wavelike solution for the lattice vibration: Here, Λ and j are the amplitude and phase of the wave, M μ and e j m ( ) denote the mass and polarization vector of atom μ and ω j =2π ν j is the vibration frequency.
This expression of a displacement is utilized to develop a specific lattice vibration as follows: for chosen values of j and q, we calculate the real part of the time derivative from (1) Then, we add the incremental velocity given by to all atoms at the end of each MD-step. Each oscillation period, given by q 1 j n ( ( )), is discretized into time steps t n =t 0 +nΔt, where n n t 1, , 1 j in order to maintain continuous oscillations. The excitation rate is then defined by where NΔt is the waiting time between two excitations (i.e., N is the number of MD-steps between two excitations). The phases of these excitations are uncorrelated since we choose a random value for Here, q D( ) is the dynamical matrix, which is determined numerically (for aluminum) by means of MD simulations as described in [16]. We are now able to calculate all quantities in (1) to simulate the excitation of a lattice vibration. The dispersion of the vibration frequencies q j n ( ) at T 300 K = is shown in figure 2. These results compare well with the experimental data in [17]. They serve to give the frequency of the vibrations which allows us to study whether or not lattice vibrations at the Brillouin zone edge have an influence on defect generation.

Results and discussion
Our MD simulations show that lattice defects are generated if three necessary conditions are satisfied: (i) the specimen temperature T is above the Debye temperature θ D , which is 428 K D q = for Al [14], (ii) the phonons n ( ) at T 300 K = with the method described in [16] (solid black line). Some dispersion curves are degenerate (as indicated j=1, 2). The red data points indicate experimental data taken from [17]. The green and blue circles indicate the modes that are being excited in this work.
proliferate at a high rate, q, and (iii) the phonon wavelengths, q, are near the edge of the Brillouin zone. We will establish these three necessary conditions in sections 3.2-3.4 by performing a parameter study of T 0 , q and q.
But first, in section 3.1, we illustrate the way we detect and visualize lattice defects in simulation snapshots. This is done for a specific case study where the parameters are chosen such that lattice defects can be generated. This method is applied in all other sections. We find, that the generation of lattice defects can quantitatively be measured by kinetic energy versus time plots, where a time of t 0 ns = corresponds to the time after the system had been equilibrated at T 0 .
For the simulations described in the following sections, the chosen values of T 0 , q and q are summarized in table 1. In addition to these variable parameters, we keep j=1 and M 0.1 L = m Å fixed in all of our simulations. For the transversal acoustic mode, j=1, the experimental and simulated dispersion curves agree very well (see figure 2). The amplitude Λ is chosen such that it is far below Lindemann's melting criterion [18].
The temperatures T 0 are chosen to be below and slightly above the Debye temperature θ D . The rates are chosen such, as to show a difference in the generation of lattice defects for high and low excitation rates. The 3.1. The generation of lattice defects shown in a specific case study The parameters (see table 1) for this specific case study are chosen such that the generation of lattice defects can be observed. In figure 3, the temporal evolution of the kinetic energy is shown as a blue curve. Since we are adding kinetic energy to the system by exciting lattice vibrations, the kinetic energy increases rapidly; but interestingly, this happens only until a time t 0.09 ns D » , although the excitation-and therefore the adding of energy to the system-is not interrupted. After t D it is observed, that the kinetic energy increase stagnates. As it is shown by the red curve of figure 3, the energy increase ceases when the generation of lattice defects, more precisely, Frenkel-interstitials begins. At first, we will illustrate, how we find these lattice defects. Then, we will discuss the correlation between the energy leveling off and the generation of lattice defects.  have been calculated numerically with the method described in [16].  After every 5000 MD-time-steps we save a snapshot of the current atom configuration. In figure 4(a), such a snapshot can be seen. The configurations are then quenched down to T 0 K = . This is done by a minimization of the system's potential energy. Subsequently, we perform a common neighbor analysis (CNA) [19] with the method implemented in OVITO [20] to characterize the atom structure. The CNA of the minimized system can be seen in figure 4(b). In a perfect crystal without defects, all atoms-except the ones on the surface-would be characterized as fcc atoms. But in our case, we find that about 0.6% of the atoms are characterized as non-fcc atoms. In figure 4(c) all perfect fcc atoms were removed from the snapshot. Here it can be seen clearly, that our system has vacancies and interstitials. By means of the Wigner-Seitz defect analysis implemented in OVITO [20], we have counted 12 Frenkel-interstitials (i.e. 12 vacancies and 12 interstitials). In addition to the bulk defects, we also observe the generation of adatoms at the free surface. Their generation is either due to migration of bulk defects to the surface or due to direct formation of vacancies at the surface. The latter case is observed more often, since the formation energy for such adatom-hole pairs is about an order of magnitude smaller than for Frenkel pairs in the bulk.
Note that the number of crystal defects is far from equilibrium. In this specific case study, the crystal defects are generated at a temperature of about T 430 K » . By extrapolating experimental data taken from [21], we determine the ratio of vacant sites to the total number of atoms in Al to be n N 1. =´-, which is about 5 orders of magnitudes higher than the concentration of vacant sites in an equilibrium crystal. The stagnation of energy increase in figure 3 can be explained as follows: the energy injected when exciting lattice vibrations can be dissipated as heat or, alternatively, as formation energy of lattice defects. In the former case, the kinetic energy increases. After t 0.09 ns D » , crystal defects are generated, which explains that the kinetic energy levels off.

The Debye temperature as a lower bound for the defect generation
In this section, we will show for Al, that crystal defects in concentrations far beyond equilibrium are only generated above the Debye temperature 428 K D q = [14]. This does not mean that defects cannot be generated at lower temperatures at all, but this is significantly less likely. In our driven MD simulations defects are generated by non-equilibrium energy fluctuations assisted by the sample temperature. We will argue now, that in this scenario it is not by accident, that the system's temperature has to be above the Debye temperature to generate a non-equilibrium concentration of crystal defects. To show this, we have chosen values of T 0 as given in table 1).
In figure 5, the results of our simulations for different starting temperatures T 0 are shown in a kinetic energy versus time plot. At first, we will discuss the four curves where T 0 <θ D was chosen. As it was already discussed in the previous section, the kinetic energy for these simulations increases; but interestingly, this, again, happens only until a time t 0.09 ns D » . After t D it is observed, that the kinetic energy does not increase anymore, but reaches a plateau for all different values of T 0 , indicating the generation of lattice defects (as it has been described in the previous section). More interestingly, the temperature of this energy-plateau, and thus the generation of lattice defects, corresponds to the value of the Debye temperature 428 K D q = for Al [14]. The role of the Debye temperature can be explained as follows: the Debye temperature characterizes a threshold above which all possible lattice vibration modes in an equilibrated crystal are thermally excited. Below (a) Simulation snapshot at t 0.13 ns » . This snapshot is quenched down to T 0 K = , so that the atoms are placed in a local potential minimum. Subsequently, we perform a CNA as shown in (b). If we delete all perfect fcc atoms, vacancies and interstitials can clearly be seen in (c). By a Wigner-Seitz defect analysis, we count 12 Frenkel-interstitials for this simulation snapshot. The visualization and analysis has been realized with OVITO [20]. θ D , there are lattice vibration modes-for example those near the Brillouin zone edge-which are only exponentially weakly occupied in an equilibrated crystal. In the simulations discussed here, non-equilibrium lattice vibrations with a wave vector of q 0.98, 0, 0.98 a = p ( ) are generated. We propose that, when T 0 <θ D , the crystal lattice is too cold to keep these edge modes for long enough; in other words, below θ D , the excited edge modes decay, and the injected energy is just dissipated into heat. As soon as the Debye temperature is reached, the lattice is hot enough that the additional non-equilibrium excitations of Brillouin zone edge vibrations lead to the generation of crystal defects.
We will investigate this assumption in more detail: instead of T 0 <θ D , we choose T 430 K 0 D q = > and leave the excitation rate at q 0.42 THz = . In this case, the excitation of edge modes does not heat up the lattice anymore, but primarily generates lattice defects. The kinetic energy curve essentially is a plateau during the whole simulation. This is different if T 0 is even higher. For a starting temperature of T 600 K 0 = , we find that eventually, the kinetic energy decreases to a value corresponding to the Debye temperature, although the system is additionally driven by exciting short wavelength lattice vibrations. Our explanation of this at first sight paradoxical phenomenon is the following: the homologous temperature of the short wavelength lattice vibrations decreases due to defect generation below the temperature T 0 . As a consequence, energy will be redistributed from the hot lattice modes of longer wavelengths to the cool ones near the Brillouin zone edge. This leads to a cooling of the lattice as a whole until the Debye temperature is reached. If, however, the initial temperature is T 700 K 0 = , then the crystal melts. We conclude that in this case, the defect generation is too slow to create a negative temperature gradient between the long and short wavelength lattice vibrations.

The influence of different excitation rates q
In figure 6, we have plotted the effect of different rates (see table 1) on the generation of lattice defects. As highlighted, depending on q, we have observed two events: the generation of lattice defects or the melting of the crystal without the generation of lattice defects. The five curves, where the generation of lattice defects has been observed (q 0.34 THz, , 0.42 THz show, that the rate has an interesting effect: the lower we choose q, the later the generation of lattice defects starts. This can be seen by the plateau of the kinetic energy, which occurs later and at higher temperatures the lower we choose q. The trend of this rate-dependency leads to the question, what the necessary conditions are, to actually generate lattice defects. For example, does decreasing or increasing the rate even further, also lead to generation of lattice defects eventually? The answer is also given in figure 6. We have performed simulations with a very high (q 2.2 THz = ) and a low rate (q 0.22 THz = ) at T 430 K 0 = and q 0.98, 0, 0.98 a = p ( ) . In both cases, we observed qualitatively the same: the kinetic energy continuously increases until the lattice is so hot, that it starts to melt. During the heating of the lattice, no Frenkel-interstitials were found. This leads to the following conclusion: if the rate is too low, the crystal has enough time to equilibrate between two excitations. The phonon distribution is not dominated by the excited mode anymore, but due to phonon-phonon scattering, all vibration modes are approximately occupied as in an equilibrium state. The excitation of a lattice vibration at a low rate thus only leads to heating the lattice until the crystal melts. Yet, if the excitation rate is too high, the lattice does was chosen. We found that the beginning of the plateau is highly correlated with the generation of lattice defects, as highlighted in the plot. The corresponding temperature for the generation of the lattice defects is about at the Debye temperature 428 K D q = for Al [14]. For T 700 K 0 = , the crystal melts due to the driving of the lattice. not have enough time to generate a crystal defect. Therefore, we also see a melting of the system for very high rates.

The influence of a lattice vibration mode close to the Γ-point
In addition to changing q, we have also investigated, how the choice of q affects the generation of lattice defects. For this purpose, we have changed the wave vector to lie close to the Γ-point, and therefore close to the middle of the Brillouin zone (see table 1). We found that, independent of q, the crystal only heats up and that no Frenkelinterstitials are generated during the heating. Representatively for all different rates, this phenomenon is shown by the red curve in figure 7 for q 0.38 THz = and T 430 K 0 = . This curve indicates, that the system just heats up until it eventually melts. For comparison, the blue dashed line shows a simulation for the same rate q 0.38 THz = , where the wave vector was chosen to lie close to the Brillouin zone edge, q 0.98, 0, 0.98 a = p ( ) . We conclude, that it is necessary, to excite a lattice vibration mode near the Brillouin zone edge in order to generate lattice defects.

Conclusions
By means of MD simulations of single-crystalline Al we showed, that driving a lattice out of equilibrium by exciting certain phonons more than others can lead to the generation of interstitials and vacancies before the crystal melts. The number of lattice defects was found to be far beyond the equilibrium concentration; we found concentrations about 5 orders of magnitude higher than in thermal equilibrium. Three conditions are met, when crystal defects are generated: the temperature of the system is above the Debye temperature, the lattice and different values of q. As highlighted, depending on q, two events can be observed: the generation of lattice defects or the melting of the crystal without the generation of lattice defects. vibration mode lies close to the Brillouin zone edge, and the excitation rate is high enough, so that a nonequilibrium phonon distribution exists.
Our findings are relevant to recent results on flash experiments (for an overview on flash sintering, see e.g. [22][23][24]): in [25] Yadav and Raj hypothesized, that the start of the flash event is related to the onset of 'nonlinear lattice vibrations', similar to the proliferated Brillouin zone edge modes we have excited in our MD simulations. Moreover they found, that the Debye temperature is a lower bound for the onset of the flash for poly-and single crystals of zirconium and titanium. It is assumed [26], that nonlinear lattice vibrations soften the shear modulus of the crystal, leading to the generation of Frenkel-interstitials. While this issue remains controversial, the present work suggests the significance of Frenkel pair generation in flash sintering [27][28][29]. Therefore, our results may be relevant to understand the mechanisms for recent flash sintering experiments on Al alloys, performed by McWilliams et al [5].
The present work analyses a mechanism, how driven lattice vibrations spend their energy on the generation of point defects. The driving mechanism itself, however, was not addressed. In order to substantiate the relation to flash experiments, it will be vital to explain the driving mechanism in terms of the electrical field applied to the sample. In this context an important hint is given in a recent paper by Meldonado et al [30]. The ab initio calculation shows that electrons scatter particularly strongly from phonons at the Brillouin zone edge. This suggests that a high current density selectively excites short wavelength phonons. This happens on femtosecond timescale while the phonons have a lifetime of several picoseconds [30] such that they can proliferate, as needed for the mechanism explained in the present work. This idea remains to be worked out in the future. It should be noticed, however, that such a scenario is not in conflict with a recent calculation for the insulator HfO 2 that concluded that field-enhanced generation of point defects was negligible [31].