Frenkel pairs cause elastic softening in zirconia: theory and experiments

Recent results from molecular dynamics simulations have shown that significant concentrations of Frenkel pairs can be introduced by the proliferation of phonons lying at the edge of the Brillouin zone and when above the Debye temperature. Following the work of Granato (2014 Eur. Phys. J. B 87 18) we extend those calculations to the influence of Frenkels on the elastic modulus. Significant softening is predicted which is confirmed by in situ measurements of the elastic modulus during flash. Frenkel pairs have been proposed to play a central role in the flash phenomena.


Introduction
In 2018 the authors have published results from molecular dynamics (MD) simulations [1] that predicted the generation of Frenkel pairs by the proliferation of phonons with wavelengths lying at the edge of the Brillouin zone when the temperature rose above the Debye temperature. Results from this first work, on aluminum single crystals, are reproduced in figure 1. (Similar results were obtained in titania, for both cations and anions, that were published soon thereafter [2].) In this figure we note that: (i) injection of phonons at sufficiently high rate and at temperature below the Debye temperature θ D raises the temperature but only until it crosses the Debye threshold, where it levels off and is accompanied by the formation of Frenkel pairs in increasing concentrations, (ii) if the starting temperature is above Debye-but not too close to the melting temperature-the temperature falls and then hugs Debye again producing Frenkels. The Debye asymptote of the temperature suggests that the injected energy is consumed by endothermic formation of Frenkel pairs.
In the nineties Granato [3] was intrigued by the local strain field of Frenkels in metals and argued that they increased the entropy thereby lowering the melting temperature. He considered thermodynamic equilibrium for the formation of Frenkels, which having a high energy of formation, would develop close to the melting temperature. Thus, the lowering of the melting point would have been slight, which could not be confirmed by experiment.
However, the analysis of the strain field of Frenkels and their influence on the elastic modulus emerged as a standalone issue. Granato predicted a softening of the shear modulus with increasing concentration of interstitials, c I , which he described by [4] G (c I ) /G 0 = exp (−βc I ) . (1) Thus the influence of the Frenkels on elastic behavior emerges as an area of research on its own. It is relevant to Frenkels produced not only in radiation damage but also by another athermal mechanism, as hypothesized during flash experiments, where it was demonstrated that, in three instances, flash occurred only if the specimens were above the Debye temperature [5,6], as predicted in figure 1.
The hypothesis of Frenkels in flash experiments was also supported by comparing the measurements of excess lattice expansion (measured in situ with x-ray at Argonne National Laboratory) with predictions from ab initio calculations [7]. More recently measurements of an energy deficit between the measured temperature and the temperature predicted from a black body radiation model, gave quantitative agreement between the extent of defect generation and the acceleration in the sintering rate [8].
In this article, the above theoretical and experimental results are extended to the influence of Frenkels on the elastic modulus. We begin with MD calculations to predict the elastic modulus with Frenkels and compare subsequently with in situ measurements of the modulus in specimens held in a steady state of flash (so called stage III). The experiments were carried out by measuring the change in the free vibration frequency of cantilever beam specimens. Baseline data were obtained with furnace heating. They were compared with measurements of the modulus during in-flash experiments. The results from MD simulations gave good agreement with experiment.

The method
We did MD simulations of ZrO 2 in order to answer the following questions: can the presence of Frenkel defects account for elastic softening in zirconia, and if so, how does the extent of softening depend on their concentration?
The simulations were done with the program LAMMPS [9]. As interaction potential we took where q i and q j are the effective charges of the ions i and j, and r ij is their distance. The parameters A ij ,C ij and ρ ij for Zirconia were determined in [10]. For the two non-Coulombic terms in equation (2) we chose a short-range cut-off distance of 10 Å. The Coulomb interactions were evaluated by means of the Ewald sum method [11] with a precision of 10 −5 . We constructed initial configurations with Frenkel defects by the following procedure: following [10], we started from a perfect cubic lattice of ZrO 2 with periodic boundary conditions. It consisted of 10 × 10 × 10 unit cells with a lattice constant of a = 5.082 Å, each containing 4 Zr-atoms and 8 O-atoms, so that the total number of molecules was N = 4000. The zirconium and the oxygen atoms sit in alternating planes separated by the vector 0.25a(0, 0, 1). In the zirconium planes the atoms form a square lattice with relative coordinates of nearest neighbors 0.5a (1, 1, 0) and 0.5a (−1, 1, 0). The oxygen atoms form a square lattice with relative coordinates 0.5a(1, 0, 0) and 0.5a (0, 1, 0). Then we selected n Zr zirconium and n O = 2n Zr oxygen atoms at random and displaced each of them by the vector r I = 3.25a(1, 1, 0) from their perfect lattice positions to preliminary off-lattice ones. This configuration was then allowed to relax anisotropically at fixed temperature T = 1400 • C and zero pressure [12] for 0.3 ns allowing the displaced atoms to settle into preferred interstitial positions. Choosing different random seeds we produced six initial configurations for n Zr = 3, and four for n Zr ∈ {15, 25, 45, 60} this way.
During the relaxation the initial cubic structure transformed spontaneously into a tetragonal structure, as also reported in [10], regardless of the concentration c I of interstitials. Both, the vacancies and the interstitials are mobile during the relaxation, leading to clustering and partial recombination. The surviving numbers of Zr-and O-interstitials, n true Zr and n true O , were measured for each initial configuration by means of the Wigner-Seitz defect analysis [13] (after a quench to T = 0 at zero pressure). The molar concentration of interstitials was then calculated by assuming that no crystal defects are lost during the quench. Deviations from n true Zr = 2n true O occurred, but they were not systematic.
For each relaxed initial configuration (not quenched) the six independent single-crystal elastic constants C 11 , C 12 , C 13 , C 33 , C 44 , and C 66 (Voigt notation [14]) were calculated. These results are reported in supplementary information (https://stacks.iop.org/NJP/23/053013/mmedia). In order to compare with the experimental results for polycrystalline ZrO 2 we used the Voigt-Reuss-Hill approximation [15] to calculate the shear modulus, G, and the bulk modulus, B [16], where 2 13 , and, M = C 11 + C 12 + 2C 33 − 4C 13 . The Young modulus E was calculated from bulk and shear modulus: We obtained the elastic constants from the slope of the stress(σ)-strain(ε)-diagrams according to Hook's law in Voigt notation [14], The initial configurations were deformed in an MD-simulation without a thermostat and a barostat. For the dilatational strains ε 1 , ε 2 or ε 3 we increased the total x-, y-respectively z-length of the simulation box by Δ = 0.5 · 10 −3 Å in each MD-timestep t n and rescaled the x-, y-respectively z-coordinates of all atoms accordingly. For ε 4 or ε 6 we sheared the yz-respectively xy-face of the simulation box by Δ in each MD-timestep t n and imposed corresponding displacements of the atomic coordinates, y (t n+1 ) = y (t n ) + Δ L z z (t n ) leaving the x-and z-coordinates unchanged (for ε 4 ); and x (t n+1 ) = x (t n ) + Δ L y y (t n ) leaving the y-and z-coordinates unchanged (for ε 6 ). For each timestep the stress components were measured as the sum of a kinetic and a virial contribution as implemented in LAMMPS.
During the elastic deformation we monitored the temperature (derived from kinetic energy as usual). Since the deformation happens adiabatically, the temperature decreases slightly with increasing strain, as shown in figure 2 (left). The decrease is less than 3% for maximal strain.

Results from MD simulations
A specific example of how we calculated one of the elastic constants, C 11 , is shown in figure 2 (right). It shows the stress-strain curves for a perfect crystal (red line) and a crystal with an interstitial-concentration of c I = 0.03475 (blue line). The elastic constant C 11 is the slope.
In figure 3, the shear modulus G, bulk modulus B and Young modulus E are plotted as functions of the concentration of interstitials c I . The error bars are calculated from the errors in the linear fits to the stress-strain curves. Data with the same value of n Zr are marked by the same color. Their spread indicates the variability of the true interstitial concentration (3) after relaxation of the initial configuration, arising from different degrees of recombination.
The simulation results show that the zirconium oxide crystal softens in the presence of interstitials. The data can be well fitted by the functions  which in contrast to equation (1) takes the moduli for amorphous zirconia, G am , B am , and E am into account for large interstitial concentrations c I . The fit parameters are (7) calculated with equation (4) as a function of the molar concentration of interstitials c I . The result shows that the softening of the shear modulus increases with the concentration of interstitials. The full lines are the fits, equation (6), to all data.
A linear approximation of (6), G (c I ) ≈ G 0 (1 − βc I ) leads to the shear susceptibility constant β = 51, while the linear fit to the data up to c I < 0.006 gives β = 33.79, which shows how sensitive this parameter is to the way the fitting is done. The absolute values for the defect-free case, G 0 = 137 GPa, and B 0 = 251 GPa are larger than the experimental values. This is probably due to the choice of the potential, equation (2), an experience in line with previous studies [17].

Experiments
The experiments measured the change in the vibration frequency of a cantilever beam, with and without the electric field applied to the mid-portion along the length of the specimen. The change in the frequency was related to the change in the modulus with modal finite element analysis, which is described in detail in the following section.
The flash behavior of ceramics is orchestrated by placing the specimen within a furnace held at a specific temperature and applying an electric field. After an incubation time, there is an abrupt rise in conductivity signaling the onset of the flash. This is called stage I. The rise in current is controlled by switching the power supply from voltage to current control. This transition has been called stage II. Thereafter the specimen can be held nearly indefinitely under current control, which has been known as stage III [18]. While the original experiments were designed for flash sintering of powder pressed samples, the more fundamental studies of the flash mechanism, as in the present work, have been carried out with dense polycrystals. Flash experiments with single crystals have provided further insights into the science of the flash phenomenon [19].
It is of note that flash in stage III can be maintained with the current supplied directly to the specimen from a power supply operated under current control, without any auxiliary heating. At this point the furnace can be switched off. The in situ experiments described below were carried out with furnace-off, which rendered considerable simplicity to the measurement of the specimen temperature.
The objective of the experiments was to measure the change in the elastic modulus under the influence of flash, relative to the intrinsic modulus of zirconia. Therefore, the experimental protocol consisted of two steps. First the intrinsic modulus was measured with furnace heating-these results compared reasonably well to literature data. Next the specimen was activated with flash and the modulus measured as a function of temperature. These two measurements were compared at same temperature, the first obtained with the furnace and the second with heating from the current deployed to sustain the flash.

The method
In situ elastic modulus was measured from the change in the frequency of (free) vibration of a cantilever beam specimen. The change in the frequency was related to the change in modulus with the help of finite element modal analysis, as explained later.
Displacement-time signals for the vibrations were collected by a Polytech Portable Digital Vibrometer (PDV-100). It uses a helium neon laser of wavelength 633 nm with a frequency measurement ranging from 0.5 to 20 000 Hz. It was kept at a distance of about 20 inches from the sample as proposed by Polytech in order to get optimal measurements. An impact model-hammer (PCB 086C02) from PCB Piezotronics was used to give an impulse to the sample so that it would vibrate in a direction parallel to the laser beam. A Polytech Data Acquisition System (VIB-E-220) was used to acquire data from both the Vibrometer and the model-hammer. VibSoft 5.3 and VibSoft 5.3 Desktop (software) along with keys (VIBSOFT-20 and VIBSOFT-desktop) were needed to access the data on a computer.
The beam was vertically aligned. It was held securely at the bottom. The top, the free end of the beam, was ground to a flat surface. A reflective tape was attached to this surface. The laser beam was focused on this flat surface so that the reflection could be read by the vibrometer. In this way the back-and-forth vibrations of the free end of the beam parallel to the laser beam could be specifically measured.
The experimental set-up is shown in figure 4. A rod of 3YSZ, 1.83 mm in diameter and 137.6 mm long, obtained from CoorsTek, Inc., Golden, CO, served as the cantilever. The beam was fixed at the bottom with the top end being free to vibrate. The fixed end was secured within a short piece of alumina tube with Omega Bond OB-600 (Omega Engineering, Norwalk, CT) high temperature cement. The sample was placed on a support structure for 24 h to allow the cement to dry after which it is was cured at 82 • C for 4 h and then further at 104 • C for another 4 h. The beam was placed within a home-built platinum furnace with a total height of 56 mm.
The vibration frequency of the beam at ambient temperature was checked against the theoretically expected value from the Euler-Bernoulli beam equation where ω n (rad sec −1 ), E (GPa), I (m 4 ), M (kg), and (m) are the natural frequency (first harmonic), the elastic modulus, the area moment of inertia, and the mass and the length of the sample respectively. Here, I   [20]. Experimental values of the frequency were determined by fast Fourier transform (FFT) with a program written in MATLAB. The first fundamental frequency of vibration was chosen. The change in frequency was translated into the change in modulus of the hot section of the beam by modal finite element analysis. The method is illustrated in figure 5. The graph in 5(a) on the left predicts the change in frequency as a function of the change in the modulus of the mid-section of the beam, based upon modal finite element analysis. The graph 5(b) shows the experimental data for the change in the frequency when the mid-section is heated to a certain temperature. Cutting across the two figures at the same frequency yielded the modulus as a function of temperature. The same approach yielded the modulus when the specimen was heated with a flash current.
The baseline temperature dependence of the modulus was measured by heating the furnace from ambient to 1200 • C. A comparison of these data with literature values is shown in figure 5(c). The difference between the current data and literature values [20] lies within 15 GPa. This difference is much smaller than the change in the modulus induced by flash, as shown later. The furnace, which was held at 750 • C to initiate the flash, was switched off upon reaching stage III, and allowed to cool down before measuring the vibration frequency. The current density was now changed in steps to change the specimen temperature. At each step the vibration frequency was measured.

Measurements of the elastic modulus with a flash current
To induce flash two thin platinum wires were attached to the mid-section of the beam creating an effective gage length for the flash current which was controlled by the power supply. Three sets of data for flash lengths of approximately 10 mm, 20 mm and 30 mm, were obtained. The wires were wrapped around the thin rod with a touch of platinum paste to promote electrical contact. The effect of the platinum wires on the vibration was determined by measuring the frequency with and without the platinum wires. At room temperature the frequency with the wires was lower by 0.5 Hz than without the wires; this drop corresponds to a difference in effective modulus of 2 GPa, which is insignificant in comparison to the softening induced by flash.
Flash was triggered by heating the furnace to 750 • C and applying a field of 100 V cm −1 . The steady state of the flash under current control (stage III) was quickly established. At this point the furnace was turned off. The specimen was now heated with the flash current supplied directly to the specimen from the power supply. Experiments were carried out for four levels of current density 120, 100, 80, and 60 mA mm −2 . The voltage expressed across the gage length, was measured. The product of the measured field and the injected current density gave the power density. The data were analyzed with the following nomenclature and units for the electrical parameters: P W is the power density in units of mW mm −3 , φ is the field in V cm −1 , and J the current density in mA mm −2 . Therefore P W = (φ/10) × J mW mm −3 . The power density is used to estimate the temperature of the specimen with a black body radiation model [21].
The field, the current density and the power density for the case of 20 mm gage length are shown in figure 6. Note that the furnace was turned off after ∼2000 s and allowed to cool down to ambient temperature before measuring the vibration frequencies. The temperature of the specimen was reduced in steps by changing the current flowing through the specimen.
The spike in the power density seen during the transition from stage II to stage III is characteristic of flash experiments [18]. The downward spike is seen when the current is dropped down to 100, 80 and 60 mA mm −2 ; but it rises nearly immediately to a steady state value. Note that the field generated across the specimen under current control does not change with the current density. The reason for this effect, which has been observed in other flash experiments in our laboratory, remains unclear. Experiments similar to those in figure 6 were carried out for gage lengths of 10 mm and 30 mm. Those results can be found in [22].
The temperature of the in-flash specimen was estimated by two methods: with a black body radiation model (BBR) [21], and then by direct measurement of the temperature with a pyrometer. The BBR model has been shown to be consistent with synchrotron measurements with a platinum standard [23]. The model invokes the concept that the black body radiation loss should be equal to the electrical power dissipation, which leads to the following expression, Here, T K is the specimen temperature and T 0 K is the ambient temperature in K. In the present analysis T 0 K = 300 K (since the furnace had been turned off). P W is the power density calculated from the prescribed gage-section, V is its volume and A the surface area of the gage-section. The Stefan Boltzmann constant is σ = 5.687 × 10 −8 W m 2 K −4 , and ε is the emissivity, assumed to be 0.9 for zirconia [24].
The error in T K for uncertainty δP W is obtained by taking logarithm both sides in equation (2), assuming that T K T 0 K , and differentiating Equation (10) was applied to estimate the error in the temperature from uncertainty in the power density. Since the power density is equal to the watts expended divided by the volume of the flashed section of the cantilever beam, an uncertainty in the gage length can lead to an error in P W . The electrical contact was made by wrapping platinum wire around the cylindrical sample by hand and applying some platinum paste to establish good contact. The error in the placement of the wires would have been most likely for the shortest gage length. Indeed this was apparent in the measurement of the power density: it was 700 mW mm −3 for 10 mm but 560-600 mW mm −3 for 20 and 30 mm gage lengths, leading the an uncertainty in the power density of about 15%. Substituting in equation (10) it would have led to an error of ∼4% in the estimate of the temperature by the BBR model.
The specimen temperature was also measured with a pyrometer (Micro-Epsilon, Model 1MH, Focus: CF4; Raleigh, NC). The emissivity at the pyrometer was set to 0.9. The pyrometer measurements, which are explained in detail in [22], are compared with the BBR estimates in figure 7.
The in situ modulus of the specimens was obtained by the method explained in figures 5(a) and (b). First the vibration frequency was predicted for different values of the elastic modulus in the mid-section of the cantilever beam. Then the vibration frequency was measured as a function of the specimen temperature. In this way, the elastic modulus was linked to the temperature, via the frequency.
Twelve sets of experiments were carried out, corresponding to three gage lengths and four current densities for each of them. These results are reported in figure 8. They show a significant softening of the modulus under flash current. Furthermore, the rate of decrease of modulus steepens with temperature whereas the intrinsic modulus tends to flatten out at higher temperature.
While the data for 20 mm and 30 mm are consistent, those for 10 mm fall significantly lower. As described earlier there was some uncertainty in the power density for this specimen but that is accounted for in the error bar for the temperature. Our view is that this difference arises from the behavior of the platinum-zirconia interface. As the gage length becomes narrower the relative influence of the interface Figure 8. The softening of the elastic modulus in stage III of flash. Note that the rate of change of modulus with temperature in flash is steeper than the intrinsic behavior which appears to be flattening. The data for the 10 mm specimen appears to give lower values than the 20 and 30 mm specimens, which is discussed in the text.  (which does not depend on the gage length) becomes more important. These measurements may suggest that the softening in the region close to the interface is stronger than within the gage section.

Measurements of damping
We have measured the damping of free vibrations. It was analyzed in two ways: (i) by the logarithmic decrement defined as δ = (1/n) n(x 0 /x n ), where x 0 and x n are the amplitudes separated by n number of oscillations, and (ii) by measuring the full width half maximum (FWHM) of the frequency peak in the Fourier transform space; the logarithmic decrement is related to the FWHM, Δf, and the peak frequency, f, by δ = 1.814(Δf/f ) [24]. The damping coefficient, ζ used in the description of the damped simple harmonic motion is then related to the loss tangent by ζ = δ/ √ 4π 2 + δ 2 . The damping in internal friction is usually expressed as a Q factor, which is related to the damping coefficient: The damping profiles, with flash and without flash are shown in figure 9. Both sets of data were obtained at approximately the same temperature of 1200 • C. Note that the amplitude decays more quickly under flash, and that the width of the Fourier peak is broader. Expansive data for different temperatures and gage lengths are given in table 1.

Discussion and conclusion
This paper presents a couple of new facts, which give further insight into the flash phenomenon in yttria stabilized zirconia. We were interested to study if steady state of flash (so called Stage III) softens the elastic modulus. The experiments measured the change in modulus as a function of temperature, first with furnace heating, and then with the flash current without furnace heating. These results are summarized in figure 8. The data show that the intrinsic change in modulus up to a temperature of 1300 • C were in good agreement with published results [20]. Next the experiments were repeated with the sample held in a steady state of flash (stage III) supplied with constant current. The sample temperature was measured with a pyrometer; and then also estimated from a black body radiation model. These two measurements are compared in figure 7; they differ by less than 10%.
The elastic modulus in the state of flash was measured at temperatures above 1000 • C, this high temperature being needed to sustain the flash. The modulus was significantly lower than the measurements of the intrinsic modulus without the flash. At 1100 • C the modulus under flash is about 13% lower, at 1200 • C about 28% lower, and for 1400 • C, with slight extrapolation it is lower by 33% in comparison to the intrinsic modulus.
In order to answer the question, whether the softening of the modulus could be explained by Frenkel defects, we calculated the elastic constants of tetragonal zirconia at 1400 • C with molar interstitial concentrations varying between 0 and 0.045. The numerical values were combined to obtain approximate values of the elastic moduli of polycrystalline zirconia. The dependence on the interstitial concentration could be fitted to an exponential decay for defect concentrations of <0.01. At higher concentrations the softening approaches asymptotic values between 20% and 30% of the moduli for a crystal without Frenkel defects. These predictions should be valid for other temperatures, as well, and are in qualitative agreement with the experimental values described above.
In conclusion, we can interpret the elastic softening observed in flashed zirconia as an indication of a concentration of Frenkel defects that increases with the current density (and hence the sample temperature) supplied to maintain stage III of flash. According to the simulation, the interstitial concentration at 1400 • C should be of the order of 0.01. Such high molar concentrations of interstitials cannot be expected to be thermally excited at 1400 • C (the applied electric fields were too weak-they corresponded to just a few milli-eV; while the formation energy of interstitials is several eV). Therefore, our results point towards a far-from-equilibrium mechanism for the generation of Frenkels, as discussed in the results from MD simulations reported in figure 1. The MD simulations invoke the concept that the current preferentially excites short wavelength lattice vibrations creating a state that is far from equilibrium. This non-equilibrium proliferation of phonons was recently shown to lead to Frenkel concentrations in titania [2] that are of the same order as the ones reported in this paper.
Finally, it is notable in figure 8 that the 10 mm specimen consistently yielded a lower elastic modulus than specimens with much longer gage lengths. We attribute this difference to interface behavior at the electrodes. As the specimen length decreases interfaces play an increasingly dominant role in the overall behavior. Thus, we anticipate that there is greater softening in the modulus near the interfaces. The defect generation near interfaces in flash experiments has been of interest recently in explaining the presence of some blackening at the cathode [25].
The measurements of damping we present are also likely to be of interest in understanding the movement of defects in the crystals. It is possible that internal friction experiments can be designed where the applied frequency can be resonated with the intrinsic jump frequency of the defects [26].
training us in its use. Finally RR is grateful to Professor R W Balluffi for directing him to the work of Professor Andy Granato at the University of Illinois on the relationship between interstitialcy and the elastic properties.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).