A model and simulation of lattice vibrations in a superabundant vacancy phase of palladium–deuterium

A one dimensional Bravais lattice model is applied to a superabundant vacancy (SAV) delta δ phase (Pd3VacD4—octahedral), in the palladium–deuterium system. SolidWorks is used to simulate the motion of atoms and ions in the lattice. These two approaches give identical results for the vibrations of the deuterons indicating that large vibrations of deuterons are possible when the microstructure is a mixture of beta deuteride and small volume percent delta SAV phase. These conditions result from the unique geometry and crystallography of δ phase. According to both the model and simulation, as the size of δ phase increases, opportunity for high amplitude vibrations of deuterons increases. Increasing temperature should have a similar effect.


Introduction
Superabundant vacancy phases (SAV) offer unique crystallography because the high levels of vacancies (∼25%) are ordered . In the ordered Pd 3 VacD 4 SAV phase, deuterium (D) occupies octahedral interstitial sites of the palladium (Pd) face centered cubic structure (fcc) Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. as positive deuterons (D + ), and vacancies (Vac − ) occupying all unit cell corners with some negative charge. This phase, called delta (δ) phase, is located on the Pd-D phase diagram [22,23] at a nominal D/Pd ratio of 1.33. The unique feature is orthogonal empty channels, along contiguous unit cell edges, 100 directions, occupied only by regularly spaced deuterons with spacing equal to the unit cell lattice parameter. If the deuterons migrate to tetrahedral interstitial sites [24][25][26], the SAV phase is δ [22,23], with empty channels or vacancy tubes with apparent enhanced electronic conduction [23]. In either case, diameters of these tubes vary periodically when traversing 100 directions along contiguous unit cells. Its diameter ranges from a minimum of 0.414 of the Pd atom at the edge midpoint of the fcc unit cell (octahedral position), to a maximum value of the diameter of the Pd atom at corners of the unit cell [23]. Ignoring hydrogen, invisible to x-rays, and with ordered Vac − , the unit cell of SAV phases is simple cubic (sc).
The properties of the SAV phases are not well known except: (1) the crystallography from x-ray diffraction, (2) unit cell dimensional behavior (contraction upon formation from the beta (β) deuteride or hydride phase), and (3) thermal desorption spectral behavior . Three reasons for lack of characterization are: difficulty achieving proper activity levels for SAV phases with electrolysis, the kinetics for their formation in the bulk is slow resulting in very small volume fractions, and their discovery was fairly recent.
The deuterons within the tubes of δ phase can be regarded (and modeled) as a case of a one dimensional Bravais lattice of ions which has been described in solid state physics texts [27][28][29]. This paper develops such a model, examines longitudinal lattice vibrations of deuterons along edges of the unit cell, and uses a commercially available solid modeling computer-aided design software package (SolidWorks published by Dassault Systèmes with settings listed in appendix A) to simulate these vibrations (frequencies, amplitudes and velocities) within the tubes.

Analysis and model development
A one dimensional row of ions in a monatomic Bravais lattice, after Dean [30] and Ashcraft and Mermin [27], is shown in figure 1(a). Each ion of mass m is tied to its neighbor with ideal Hookean massless springs with only nearest neighbor forces considered. All springs and ions are connected, so the extreme end springs are viewed as being connected back to the beginning of the string by a loop with a large number of ions with the end spring (on the right) connected to the beginning ion (on left). Alternatively, connectivity is realized by a massless, perfectly rigid bar (figure 1(a)), capable of motion, assuring any motion on either extreme end is replicated at the opposite extreme end, as in the loop version.
With these end conditions, which assure connectivity from ion to ion only, Newton's second law of motion is applied longitudinally to each ion in the string and yields a dispersion relation of lattice waves or phonons [27,28,30] relating circular frequency ω (rad s −1 ) for normal modes of vibration to the wave vector q and the ion lattice constant a as: Here ω o is (k/m) 1/2 , the square root of the ratio of the Hookean spring constant k to the mass m of the ion and is considered the fundamental circular frequency (rad s −1 ). This dispersion relation assumes a large number of ions and is suitable for monolithic microstructures, but is untenable in this model where phase boundaries pose different string end conditions, and SAV size is often limited from kinetics and processing conditions. To account for phase boundaries and finite number of ions, end conditions are altered: the extreme right and extreme left springs are fixed to rigid constraints (built-in and motionless) as shown in figure 1(b). Then the exact dispersion relation from Kittle [28] and Torre [29] becomes: where N is the total finite integral number of ions in the string. Each ion is numbered as m, so m is integral, 1 m N, and called the mode number of the vibration. This relation is exact for all m and N, even if N is small, e.g. 3 to 10, and has the conditions of figure 1(b) with end conditions fixed and excited with a forcing function on one end (see discussion section for condition of forced vibration at both ends). This is useful since the total number of ions, N 3 , in isolated particles of δ phase, and their volume (N·a) 3 , have been shown to be a very small (e.g. 0.03%) volume fraction of the bulk [22,23]. These small particles of δ phase, distributed within the bulk, are consistent with nucleation and growth mechanisms. Equation (2) is evaluated for each mode m and any N to yield a ratio of normal mode ω to fundamental circular frequency, ω/ω o . Thus the string has other normal modes, also called natural circular frequencies ω [or frequencies f = (1/2π)·ω, in Hz], in addition to the fundamental circular frequency ω o [or fundamental frequency Figure 2 shows equation (2) solutions of natural frequencies f for selected m modes and N ions. For any N, there are N natural frequencies, one for each m, consistent with string length (N·a). If the string is forced to vibrate longitudinally along its length at one of these natural frequencies, it will resonate with amplitude increasing to large values.
In the case of deuterons in the δ phase, the fundamental frequency (1/2π)·(k/m) 1/2 is its thermal vibration frequency, determined from experimental measurements. Table B1 in the appendix B lists fundamental thermal vibrational frequencies f o of isotopes of hydrogen in palladium under various conditions (phases) from the literature as well as the thermal vibration frequency of Pd. These do not account for the string nature, as in this model, but only account The uniqueness of δ phase with ordered Vac − at corners of unit cells supports this correlation to the one dimensional Bravais lattice of ions, and its description by equation (2). Density function perturbation theory (DFT) [12,[15][16][17][18][19][20][21] shows there is a binding energy between positively charged deuterons and Vac − , indicating negative charge associated with Vac − . Pd atoms   (2) one end is fixed by the Pd atom (as in figure 1(b)) and the other end is forced to vibrate by the Pd atom at that end. from β phase fix the end springs because they are 52.8 times more massive and their thermal vibration frequency is lower, but can also exert longitudinal forces on the deuteron string from thermal motion. As Pd atoms vibrate, the deuteron spring compresses and stretches producing Hookean forces on the string equal to displacement (vibration amplitude) times the spring force constant k of the deuteron's spring. In this model, the spring constant for deuteron springs is determined from f o = (1/2π)·(k/m) 1 Of all the normal modes of vibration (natural frequencies) indicated from equation (2)

Simulation results with SolidWorks
Since figure 2 (equation (2)), shows many matches of natural frequencies to Pd at 5.70 THz, a commercial software package called SolidWorks was utilized to verify if these predictions, with large values of amplitude of vibration of deuterons (resonances), could be simulated. The geometry of a ball (mass) and spring model of figure 4 was implemented in the new feature of SolidWorks, motion analysis simulation. Models with 1, 3, 7, 10 deuterons, and other values of N in the string, were executed with results shown in figure 5 through 8 for 1, 3, 10 and 7 deuterons, respectively. These figures show deuteron vibration resulting from affectedly induced frequencies of vibration of the Pd atom at the end of the string. Figure 5 represents Pd vibrating at the fundamental frequency of the deuterons: it can only have a resonance if the Pd were to vibrate in the vicinity of 14.59 THz which is not the real value of vibration of Pd (5.70 THz). This simulation with one deuteron was performed as a control simulation, testing effectiveness of SolidWorks to replicate a known natural frequency. Figure 6 shows three resonances of deuterons at which are also not near that of Pd, but are consistent with the model (no dots in figure 2 at N = 3 are close to 5.70 THz). These resonances cannot be triggered by Pd at 5.70 THz. Figure 7 for a string of 10 deuterons (1000 unit cells of δ) shows 10 resonances (1 m 10) none of which perfectly coincide with the thermal vibration of Pd at 5.7 THz: there is an amplitude minimum between that for m = 1 (at 4.153 THz) and m = 2 (at 8.221 THz). However these resonance peaks, due to their widths, almost capture 5.70 THz. As N gets larger the peaks (from equation (2)) near 5.70 THz more nearly overlap the value of 5.70 THz (referred to as 'near matches' in the previous section). Figure 8 shows a resonance peak centered on Pd vibration frequency of 5.70 THz with half max width on both sides of this frequency. This resonance peak coincides exactly with the model prediction (figure 2 with N = 7 and m = 1), demonstrating Pd atoms do induce resonances at exact matches (and by extension of peak widths, at near matches).
In addition to amplitude of vibration, the velocity of the deuterons were measured in this study. Figure 9(a) shows the velocity of deuteron 7 with respect to deuteron 6. A similar result was obtained for deuteron 3 with respect to deuteron 2. This relative velocity shows the deuterons are effectively vibrating in an optical mode (nearly 180 degrees out of phase) when their relative velocities are near max or min in figure 9(a) (they are separating at positive values and approaching at negative values). They are in acoustical mode when their relative velocities are near zero. Over time they switch back and forth from optical to acoustical mode. From figure 9(a), it can be seen that they approach one another at a velocity of 8.75 × 10 5 m s −1 , however it is possible the relative velocity can be as high as twice the absolute velocity (neighboring deuterons completely out of phase), and figure 9(b) would indicate this to be as high as 1.16 × 10 6 m s −1 or slightly higher since the peak velocity appears to be still increasing at 40 picoseconds.

Discussion
This model has several simplifying assumptions that suggest that a more sophisticated approach could reveal additional details. For example, DFT would elucidate the shielding ability of the charge of the Vac − for deuterons approaching one another on opposite sides of the Vac − and its effect on the spring constant. DFT could also examine how effective the deuterons on mutually perpendicular edges of the unit cell are shielded by both Pd atoms at the center of the face of the unit cell, 1/2, 1/2, 0, of the (001) plane and by the negative charge of the corner Vac − . Note that the deuteron at the octahedral site at the body center of the unit cell, 1/2, 1/2, 1/2, is adequately shielded from the deuterons on the octahedral sites on unit cell edges, being completely blocked by Pd atoms in nearest neighbor positions. DFT would handle the dynamics and effectiveness of the negative charge of the Vac − to shield the interaction of orthogonal strings along mutually perpendicular edges (i.e. [100] and [010], or [100] and [001] or [010] and [001] directions) of each unit cell (i.e. allowing the motion and resonance of mutually perpendicular stings to be independent of one another, or indicate specific interaction). In depth work is underway to address these, but insights gained by this model and simulation may motivate and steer future work on this phenomenon.
The possibility of buckling of the deuteron string in delta phase was compared to buckling in a model of deuterons in a crack [22,23]. Here, delta phase geometry restrains buckling because (a) positive deuterons move inside a symmetric tube (figure 4) having an inside surface with negative charge as a result of the electrons from the surrounding Pd atoms, and this negative charge supports alignment by pushing them inward toward the center of the tube, (b) deuterons gain considerable velocity, as indicated in figure 9 (b) as they approach the vacancy site at the corner lattice position, and this momentum discourages transverse motion especially where most needed, near the vacancy, and (c) the deuteron's attraction to the negative vacancy site deters transverse motion by pulling them inward toward the center of the vacancy. Future DFT calculations may also illuminate these effects.
In this model and in the simulations one end of the string of deuterons is considered fixed while the other end undergoes vibration from the Pd atom. These conditions which fit equation (2) are preserved when both end Pd atoms vibrate if the string has double length (a·[2·N + 1]), wherein the central ion is immobile (from symmetry). Simulation runs in Solid-Works with this geometry of double length confirm the central deuteron is indeed motionless even without constraint due to symmetric stimulation (loading). Thus the double length condition for N = 7 becomes two collinear and contiguous stings, each with 7 deuterons, connected by a stationary deuteron in the middle, with the right set of deuterons loaded on its right end and the left set loaded on its left end, thus each set satisfies equation (2) and the associated SolidWorks simulations reported here.
Resonance theory indicates that resonances should occur at very small amplitudes of driver vibrations when dampening is zero, as in this model. Work is underway to address effects of dampening. In this model the force exerted by thermal vibration of Pd on the ends of δ phase (at β phase border) is about 3 × 10 −10 N, determined from product of spring constant (28.

Conclusions
This model and simulation suggest SAV δ phase crystallography offers a unique geometry for vibrations of deuterons. The results suggest that as the size of δ phase (N·a δ ) 3 increases, the opportunity for high amplitude vibrations of deuterons increases because of near matches (overlapping resonance peak widths). As the temperature of the Pd host increases, the amplitude of Pd vibration increases, making the possibility of resonance more likely due to a stronger driving force. SolidWorks simulations confirm predictions in the one dimensional Bravais lattice model, validating cases with, and cases that lack, resonance from thermal vibration of Pd atoms.

Appendix A. SolidWorks reference information
SolidWorks version 2018 was used with scale factors of: distance = 10 7 , mass = 2.990 430 × 10 26 , force = 10 10 , spring constant = 10 3 , frequency = 1.090 84 × 10 −12 . Motion analysis simulation settings were: accuracy = 0.000 000 01, frames/s = 300, resolution = 50%, and computational times = 5 to 40 s with personal computer real run times lasting about 2 hours for each individual simulation condition.    [15]. This same ratio (1.57) was used to set δ D /β D at 1.57, since the f o of δ H /δ D is assumed to be better represented by f o of β H /β D than by f o of α H /α D , because both β and δ have high x values). Thus the first two entries in table B1 are from values from reference [15] and the use of this ratio (1.57) from the observations of the present work (including all the data from table B1 taken as a whole).