Enhancement of Ion Diffusion by Targeted Phonon Excitation

Ion diffusion is important in a variety of applications, yet fundamental understanding of the diffusive process in solids is still missing, especially considering the interaction of lattice vibrations (phonons) and the mobile species. In this work, we introduce two formalisms that determine the individual contributions of normal modes of vibration (phonons) to the diffusion of ions through a solid, based on (i) Nudged Elastic Band (NEB) calculations and (ii) molecular dynamics (MD) simulations. The results for a model ion conductor of $\rm{Ge}$-substituted $\rm{Li_3PO_4}$ ($\rm{Li_{3.042}Ge_{0.042}P_{0.958}O_4}$) revealed that more than 87% of the $\rm{Li^+}$ ion diffusion in the lattice originated from a subset of less than 10% of the vibrational modes with frequencies between 8 and 20 THz. By deliberately exciting a small targeted subset of these contributing modes (less than 1%) to a higher temperature and still keeping the lattice at low temperature, we observed an increase in diffusivity by several orders of magnitude, consistent with what would be observed if the entire material (i.e., all modes) were excited to the same high temperature. This observation suggests that an entire material need not be heated to elevated temperatures to increase diffusivity, but instead only the modes that contribute to diffusion, or more generally a reaction/transition pathway, need to be excited to elevated temperatures. This new understanding identifies new avenues for increasing diffusivity by engineering the vibrations in a material, and/or increasing diffusivity by external stimuli/excitation of phonons (e.g., via photons or other interactions) without necessarily changing the compound chemistry.


Introduction
Solid state mass and ion diffusion is central to many applications ranging from batteries 1 and fuel cells 2 to sensors 3 and filters. 4 In many of these applications, the performance is limited by the slow diffusion of mass or ions. While some families of solid ion conductors, such as silver iodide and Li10GeP2S12 exhibit high ion conductivities of ~1 S/cm 5 and ~0.01 S/cm 6,7 at room temperature respectively, achieving high ion conductivities for divalent ions including oxides is challenging. For the conductivity of oxygen ions exhibited by yttria stabilized zirconia 8,9 to reach 0.01 S/cm, an elevated temperature of 1100 K is needed. Higher oxygen diffusivity at a lower operating temperature could enable a lower system cost, longer life and greater proliferation. 10 Great advances have been made to increase the ion transport by cation and anion substitution to alter the charge carrier density and the diffusion pathways to facilitate low energy jumps. The classical approaches of chemical modification such as cation and anion substitution have led to tunable conductivity over several orders of magnitude in several families of super-ionic conductors, such as LISICON 6,11 and NASICON. 12,13 Thus, regardless of the mobile ion species, understanding the key factors that dictate ion conductivity is needed to devise new strategies to further increase the conductivity beyond tuning of chemical compositions.
One possible avenue to improve diffusivity without modifying the chemistry stems from considering the local ion movement in the structure. At a given temperature, ions thermally vibrate in their specific crystallographic lattice site until sufficient thermal energy is available for an ion jump. The coupled ion vibration and ion movement that is paramount for ion diffusion leads to the idea that lattice vibrations (phonons) may play a role in the ion jump. Such an idea can further be intuited by considering the role of temperature as a spectrum of scalar contributions, not as a single scalar value. Based on this more detailed view, in solids and rigid molecules, temperature is composed of a summation of individual contributions by modes of vibration. By mode of vibration, it is meant that a group of vibrating atoms can be described as a summation of collective vibrations, each with a specific frequency, often termed eigen modes, normal modes or phonons. According to the traditional description for solid state diffusion based on the Arrhenius relation 14 diffusivity increases exponentially with increasing temperature, where B k and T denote the Boltzmann's constant and the bulk temperature. The term "bulk" temperature here denotes the temperature that would be sensed if the material were in thermal equilibrium. While temperature is generally understood as a proxy for the level of excitation of the phonons or normal modes of vibration within a solid material, it should be noted that each mode has its own time varying amplitude and individual temperature. It is possible for a material to experience some non-equilibrium in the individual mode temperatures, whereby certain modes are highly excited to an effective temperature above all others (e.g., joule heating of optical modes in a transistor 15,16 ). In such a situation, a subset of modes can behave as though they are at a higher effective temperature than the bulk temperature.
In the classical limit, for an individual mode labeled n , the modal temperature n T is obtained from , where n Q and n  are the modal displacement coordinate 17 (which can be interpreted as the mode amplitude) and frequency of vibration of mode n (see the ESI for a derivation of this relation and its quantum analog). The idea of modal temperatures have also been introduced in previous theoretical studies, 18 where for the first time the two-temperature model was extended to study the weak coupling between certain groups of phonons and its effect on thermal transport in solid materials. Our focus in this study is the distinction between the bulk and modal temperatures, and how it influences solid-state ion transport. As the "bulk" temperature increases, the amplitudes of all the modes in the system increase, which facilitates ions to hop over activation barriers, thereby enhancing diffusion. There are practical ways, however, to increase the amplitude and temperature of a small subset of normal modes by exciting them with light either directly through photon-phonon coupling [19][20][21][22] or indirectly through electron-photon coupling [23][24][25] or THz electric fields 26 and keeping them in non-equilibrium with respect to the rest of the modes, without increasing the "bulk" temperature. This idea of selective vibrational mode excitation has been utilized in previous studies, as a novel tool, to modify the electronic properties of the structure through electron-phonon coupling for plasmonics, 27 photovoltaics 28,29 and charge transport applications. 30 Selective modal excitations have also been utilized to induce vibrational and structural changes in the system for phase-change applications. 31,32 Inspired by recent experimental observation of solid-state ion diffusion induced by external THz-electric field illumination, 26 in this study, we aim to examine if only a small subset of modes in the lattice are responsible for the hopping of ions, which would allow the increase of the diffusivity by exciting these subset of modes without increasing the bulk temperature. Experimental realization of such an idea would have far-reaching implications for chemical reactions, phase transitions or other related phenomena, because it would suggest that migration along a reaction/transition coordinate could be accelerated by only exciting the most important normal modes, rather than all modes via an increase in the bulk temperature. This approach would mean that a material and its container could remain "cold" but it could exhibit kinetics associated with a much higher temperature, by only making a subset of the most important modes "hot". To assess this opportunity, we use the example of ion diffusion. We (i) identify the mode level contributions to ion diffusion to determine how many modes are responsible for an ion hop in the system, and (ii) quantify to what extent diffusivity can be increased by targeted excitation of a small subset of these contributing modes without increasing the bulk temperature.
It is important to recognize that the role of specific structural modes and phonons in facilitating solid-state diffusion has been reported in other studies. [33][34][35][36][37][38][39][40][41][42][43][44] The frequency of longitudinal acoustic phonons was shown to exhibit a strong correlation with the enthalpy of self-diffusion in body-centered cubic metals. 33 In addition, octahedral rotation has been shown to promote fast oxygen ion conduction in perovskite-related phases. 34 47,48 and molecular dynamics (MD) simulations. [49][50][51][52] These methods have shown that in different materials, certain group of modes could contribute to heat transfer with higher rates, for instance, the high frequency in-plane phonons in graphene 53 and interfacial modes in semiconductor interfaces. [54][55][56] In this study, we aim to identify modes of vibration that facilitate the ion migration by directly using the knowledge of phonons and their respective eigenvectors, based on the approaches that were recently developed in the context of MD-based phonon analysis. [49][50][51][52] Here, we chose Ge-substituted Li3PO4 as our model system, since Li3PO4 is a well-known parent structure, from which numerous fast solid conductors have been developed including LGPS, 57,58 and its ion conductivity can be tuned up to ~10 orders of magnitude by aliovalent and anion substitution. 7 To identify the modal contributions to lithium ion diffusion we introduce two methodologies, one based on nudged elastic band (NEB) calculations and another based on MD simulations. Following the previously developed approaches for modal decomposition of thermal transport, [49][50][51][52] we project the atomic displacement or velocity fields onto the normal modes of vibration along the hopping trajectory obtained from a NEB calculation. The magnitude of projection determines which phonons/normal modes become excited to facilitate that specific ion hop. Our results suggest that only a small subset of phonons is responsible for ion diffusion in Ge-substituted Li3PO4 and that diffusivity can be increased only by increasing the temperature of these modes instead of the bulk temperature.

Modal Decomposition Based on NEB calculations
Ion diffusion in solids typically involves the translation of an ion from one location to another, where the starting and final locations usually correspond to ion configurations that are local minima in the potential energy surface (see Fig. 1). Consequently, each of the equilibrium ionic configurations at the beginning and end of a hop could have a different set of normal modes of vibration. Thus, although we may have one minimum energy pathway (MEP) between equilibrium configurations 1 and 2 (see Fig. 1), because of the different vibrational modes at these two energy minima, the modes that initiate an ion hop from site 1 to 2 may be different from the ones that initiate a hop from site 2 to 1. More explanation on how two different equilibrium configurations in one structure possess different normal modes of vibration is provided in the ESI. To determine which modes of vibration in a local equilibrium contribute to the ion hop along a specific pathway, we combine NEB and lattice dynamics (LD) calculations.
An NEB calculation [59][60][61][62] determines the MEP between two equilibrium configurations. Each NEB calculation results in a desired number of ionic configurations (snapshots of displaced ions) that are distributed along the MEP between the two end equilibrium configurations. According to the LD formalism, 17 any atomic displacement field can be projected onto any complete set of normal modes of vibration to determine how the normal modes superpose to recreate that exact displacement field. The magnitude of projection, as has been recently learned from modal decomposition of thermal transport, [49][50][51][52] shows the degree to which each normal mode is contributing to that displacement field, which can be quantified using the following expression, 17, 51 where i u is the displacement of atom i from its equilibrium position in the configuration at the hopping origin, interpreted as the contribution by mode n to the ion hop along its migration pathway. To perform the projection of the displaced atoms onto the normal modes, one must make a choice as to which configuration along the MEP should be used to perform the projection. Although this is somewhat arbitrary, for simplicity, we selected the configuration that is halfway between the hopping origin and the transition state (see Fig. 1). We avoided projection on points closer to the saddle point based on the practical consideration that the structure would likely exhibit modes with imaginary frequencies near the transition state, which might introduce undesirable complexity, the effect of which will be analyzed in follow on studies. Nevertheless, we examined the projections closer to the beginning stages of the MEP, compared them to the halfway point, and noticed almost no change in the contributing modes (Fig. S1), which suggests the choice of the halfway point is satisfactory, at least for the system under consideration. To calculate the modal contributions for all the possible hops in the structure, we applied the above procedure to (i) all the distinct equilibrium configurations in the system since each of them would have different vibrational modes, and (ii) all the ion hops in each of these equilibrium configurations. To test our formulations, we chose Ge-substituted Li3PO4 as an example structure. The interatomic potentials for NEB calculations and MD simulations were obtained from a well-established form potential that has been successfully tested in a number of previous studies, 64-68 including Li-conducting compounds. 66,68 This form potential is comprised of a long-range Coulomb term, a short-range Morse function, and a repulsive term. The parameters were obtained from the library of potentials developed by Pedone et al., 69 which have been successfully tested in previous MD simulations of silicates and polyanion type materials including Li3PO4. 67 The exact formulas and parameters of which are provided in Table S1. To further check the accuracy of this utilized interatomic potential, we calculated the lattice constants for the perfect γ-Li3PO4 crystal from the isobaric-isothermal MD relaxation at zero pressure and room-temperature and compared them with the existing first-principles 70 and experimental [71][72][73] values (Table S4). The comparison showed that the lattice constants in our simulations have less than 3% difference with other reported values. The pristine Li3PO4 structure was chosen to be a 3×2×1 supercell containing 192 atoms, (γ-Li3PO4, space group Pnma) with lattice dimensions given is Table S4. In the pristine structure, one P was substituted by one Ge and one Li was added in the form of an interstitial to maintain charge unique interstitial sites to add this Li interstitial. However, among these 36 interstitial sites, only 31 resulted in a stable configuration with positive vibrational frequencies ( Fig. S2 shows schematics of these distinct configurations). Each hopping event brings the structure from one of these equilibrium configurations to another.
As is shown in the results section, the dominant hopping mechanism in our Ge-substituted Li3PO4 is the interstitialcy mechanism. [74][75][76][77] To determine the phonon contributions to all the possible interstitialcy hops in the system, we need to account for all the interstitialcy hops in each of the equilibrium configurations. To do so, we first counted the nearest neighbor Li lattice sites to the interstitial ion. Then, for each of these nearest neighbor Li lattice sites (tetrahedral sites), we counted the nearest neighbor interstitial sites (octahedral sites), with the constraint that the nearest neighbor interstitial site be among the 31 stable interstitial sites. Counting for all these combinations, results in 619 distinct hops in our structure that follow the interstitialcy mechanism. The SCLD calculations determined the normal modes of vibration for all the 31 equilibrium configurations and were performed using the General Utility Lattice Program (GULP) 78

Targeted Excitation of Modes to Enhance Diffusivity
After finding the strongest contributing modes to the ion diffusion using the NEB approach, we examined if external excitation of a small subset of these modes in an MD simulation had any measurable impact on the tracer and conductivity ion diffusivities. At each time step during the simulation, the equilibrium configuration was found by checking where the interstitial ion was located (i.e., which octahedral site was occupied), as the distinction between equilibrium configurations was only based on the occupied interstitial site, and as the relative position of non-interstitial ions remained unchanged across different equilibrium configurations. By knowing the equilibrium configuration, we then excited the 5 highest contributing modes among all the contributing modes to all the possible hops belonging to that configuration, while keeping the bulk temperature fixed. After a hop happened, a new equilibrium configuration was detected, and the process continued by exciting the 5 highest contributing modes among all the possible hops associated with the new configuration.
Since ion diffusivity increases with temperature, we ensured that the bulk lattice temperature remained constant at low temperature, so that any change in the diffusivity could only be attributed to the excitation of the targeted modes and not the bulk temperature. To do so, we kept the total kinetic energy of the system constant via a velocity rescaling scheme, whereby the addition of energy to the top 5 modes was complimented by a uniform reduction in the kinetic energy of all other modes in the system. To change the temperature of mode n to a desired temperature d T , we modified the atomic velocities in the system according to the following formula, where i v is the velocity of atom i , and n Q  is the modal velocity coordinate defined by, 17 * , 1 The derivation of Eq. 3 is provided in the ESI. This rescaling procedure induces only a slight perturbation to the system given that the top 5 modes only comprise 0.86% of the modes in the system, and thus the reduction in energy for all other modes is only <1% of their energy (see the ESI for discussions about the needed thermalization power). In this study, the above mentioned rescaling procedure was applied every 5 times steps (every 5 fs).

Modal Decomposition Based on MD simulations
To directly obtain the contributing modes to the ion diffusion from MD simulations, we used the definition of diffusivity based on the fluctuation-dissipation theorem 80,81 and followed a modal decomposition approach similar to that employed in MD-based thermal transport studies. [49][50][51][52] Consistent with prior work on phonon transport, [49][50][51][52] we have termed this formalism Mass Diffusivity Modal Analysis (MDMA). In MDMA, the modal contributions to the atomic velocities are obtained and then used to decompose the following definition for conductivity where c N is the number of hopping particles in the system,  denotes the auto-correlation operator, and c v is the center of mass velocity for the hopping particles, To calculate the individual contribution by mode of vibration n to Replacing i v in Eq. 6 with its modal contributions   Similarly, the modal contributions to ion conductivity   n  can be obtained, the details of which is provided in the ESI.
By substituting for both center of mass velocities in Eq. 5 with their modal contributions from Eq. 7, the above formulations can be extended to capture the modal contributions to diffusivity in more detail, where individual contributions from the correlation/interaction of eigen mode pairs n and ' which provides additional insight into the degree to which each pair of modes interact, i.e. anharmonicity via terms where ' n n  , to facilitate diffusion.
The combination of the MD and NEB based modal decomposition methods presented herein provide a detailed picture for the ion diffusion. Although both methods are based on projecting the trajectory of the hopping ion on the eigenvectors of vibration, the way they sample the phonon contributions to the ion hop is different. In the NEBbased approach the contributions are calculated based on the atomic displacement field, while in the MD-based method, the contributions are obtained based on the atomic velocity field. In addition, MD samples the dominant mechanisms that are present in a natural simulation of the material, while NEB can be used to study specific mechanisms and ion hop along specific migration pathways, regardless of the likelihood that they would occur at a given temperature.

Diffusivity calculations
Total mean squared displacment method (TMSD) was utilized to measure the diffusivity of Li ions in our structure, the details of which is explained in the existing literature 86,87 and in the ESI. First, the structure is relaxed under the isobaric-isothermal ensemble for 1 ns at zero pressure and T=400 K and under the Nose-Hoover canonical ensemble for another 1 ns at T=400 K, then the structure was simulated for another 10 ns under the same ensemble to obtain the needed dispalcement data for diffusivity calcualtions. The timestep was chosen to be 1 fs, and statistical uncertainty was reduced by considering 5 independent ensembles. 88 For calculating the conductivity diffusion coefficients, total conductivity was calculated following the existing literarure, 84 the details of which are also explained in the ESI.

Results and discussions
Modal contributions obtained from NEB calculations of Ge-substituted Li3PO4 The modal contributions to one example hop in Ge-substituted Li3PO4 shows that a small subset of modes (around 10 THz) contributed much greater than the rest of the modes (Fig. 2a). The dominant hopping mechanism in Ge-substituted Li3PO4 structure, is the concerted hop 89 of two Li ions, one of which is the interstitial, as shown in Fig. 3c and compared to two other hopping mechanisms in Figs. 3a and b. The concerted hop mechanism is known as the interstitialcy mechanism, 74-77 which has been reported 70, 90 as the dominant mechanism for lithium ion diffusion in Li3PO4 and is used as the migration mechanism for the NEB-based modal analyses presented in this study. The normalized integrations of the contributions with respect to frequency (the accumulation) are shown as a continuous curve in Fig. 2b for all the 619 hops in the structure. The majority of the contributions to all the hops in the system come from modes between 8 and 20 THz (Fig. 2b). Although this range is centered around the attempt frequency of lithium ion (~14 THz, see Fig. S3 and the ESI for the calculation procedure), the large range of frequencies for the contributing phonons shows that more intricate phonon processes are at play, which differs from the simpler picture described by traditional approaches, 91 wherein only one value of vibrational frequency (attempt frequency) enters the formulation. By averaging over all the modal contributions calculated for all the hops in the system (Figs. S4), we confirmed that on average >87% of the lattice energy during the ion hop in the structure came from <10% of the modes in the system (Fig. S5). In addition, as is explained and shown in Fig.   S11, in each hopping event, there is at least one mode that has a contribution more than 10 times the average contribution by all the modes. Thus, the results in Fig. 2 show that in Ge-substituted Li3PO4, only a small number of modes are responsible for the interstitialcy diffusion mechanism in Fig. 3c.  Due to the larger amplitudes of vibration that localized modes attribute to the interstitial and other neighboring ions than the delocalized modes 92 (Fig. 4), localized modes might be the most strongly contributing modes to each ion hop. However, inspection of the eigenvectors for the top two contributing modes in Fig. 2a (at 9.64 THz and 10.94 THz) revealed that highly contributing modes can be both delocalized (Fig. 4a) or localized (Fig. 4b) in nature. Since the modal analysis formalism is based on the SCLD calculations combined with Eqns. 1 and 2, all types of vibrational modes in the system with no inherent assumptions about their character were identified from which the contributions to ion diffusion even for localized modes were computed. Having access to such a formalism is crucial, because by breaking the symmetry in the system (e.g., by inclusion of vacancies, interstitial and substitutional ions, which is common in highly conductive solid electrolytes 93 ) the conventional picture of propagating phonons in pure crystalline solids breaks down as other non-propagating and localized modes must be included. 54,63,94 Figure S6 shows the spectral distribution of all the localized modes of vibration for all the equilibrium configurations in our Ge-substituted Li3PO4. Diffusivity enhancement by excitation of the highly contributing modes Figure 5 shows that increasing the energy of the top 5 contributing modes among all the possible hops belonging to the equilibrium configuration present at each time step during an MD simulation, by increasing their temperature to 500 K and 700 K, increases the diffusivity by 2 (~10 -10 cm 2 /s) and 4 (~10 -8 cm 2 /s) orders of magnitude respectively, compared to its value (~10 -12 cm 2 /s) in the unperturbed (natural) simulation at a lower temperature of 400 K. In fact, when the lattice is kept at 400 K and only the top 5 most contributing modes are excited to a particular temperature, the diffusivity becomes what it would have been if the entire material were heated to that temperature. Unexpectedly, even if we only excited the single highest contributing mode in each of the equilibrium configurations, almost the same increase in the diffusivity was observed. The small difference between the increase in diffusivity by exciting the 5 highest contributing modes and the one by exciting the highest contributing mode can be explained by noting that more kinetic energy is input to the hopping ion by exciting 5 modes. How these excited modes interact with the rest of the unexcited modes in the system is an interesting topic for a future study, which can be investigated using the two-temperature model that was presented in a recent study. 18 To illustrate that the increase in diffusivity due to excitation only happens when the specific modes that matter most are excited, Fig. 5 shows that exciting 5 random modes in the 8-20 THz frequency range did not lead to significant enhancement in the diffusivity. The novelty of discovering that only a small subset of vibrational modes are responsible for the orders of magnitude increase in diffusivity by modal excitation can be understood by noting the theoretical definition of diffusion coefficient in solid materials, . 91 In this definition, To check if traditional parameters for diffusivity 96 can explain the enhanced diffusion rates caused by modal excitation in Fig. 5, we quantified the changes in attempt frequency (Fig. S3), jump rate (Table S2), radial distribution function (RDF) (Fig. S7), amplitude of vibration (Fig. S8), Haven ratio (the ratio of the tracer diffusion coefficient to the conductivity (total) diffusion coefficient) 84 (Table S3), and the migration barrier (Fig. S10).
However, the very small changes in these parameters were not able to explain the orders of magnitude increase in diffusivity observed by modal excitation (Fig. 5). This also highlights the novelty of the mechanism discovered herein, illustrating the importance of modal temperatures (i.e., modal energies) and the concept of directionality brought into the analysis by including the eigenvectors of vibration in our phonon-based methodologies. Overall, the 2-4 order of magnitude increase in diffusivity by modal excitation (Fig. 5) can be attributed to increasing the energy of the highly contributing modes in the structure, where these modes can push the ion forward along its migration pathway. These results not only show that the identified modes of vibration using our proposed methodologies are correct and are indeed responsible for the ion diffusivity, but also offer a new strategy to increase the diffusivity of known solid-state ionic conductors by using techniques such as direct [19][20][21] and indirect [23][24][25] coupling to phonons without increasing the temperature of the lattice. Figure 5. By exciting the 5 highest contributing modes or even the highest contributing mode to ion diffusion in the Gesubstituted Li3PO4 structure to 500 K or 700 K and keeping the lattice at 400 K, the diffusivity increases by around 2 and 4 orders of magnitude compared to the unperturbed/natural simulation at 400 K. Exciting 5 random modes to 700 K, however, does not lead to a noticeable increase in diffusivity. The reported values are tracer diffusivity, but the same increase was also observed in the conductivity (total) diffusivity 84 (see Table S3). The activation energy obtained from the linear fitting of the Arrhenius equation 86  The modal contributions to the diffusivity for the Ge-substituted Li3PO4 structure, calculated from the MDMA method, are shown in Fig. 6b. The increase in the diffusivity accumulation function in Fig. 6b shows that the modes with frequencies between ~8-20 THz are responsible for ion diffusivity. This observation is in excellent agreement with the modal contributions obtained from the NEB-based approach (Fig. 2b), which is also averaged and shown in in Fig. 6b. The accumulation of the partial density of states (DOS) for lithium ions shown in Fig. 6b, obtained from the partial DOS in Fig. 6a, also shows reasonable agreement with the modal contributions obtained can be observed better in Fig. 6d after omitting the diagonal terms (i.e., artificially setting them to zero). Even after removing the diagonal terms, all the large contributions are still close to the diagonal, which shows that all the anharmonic terms are also between the modes that have relatively similar frequencies. Quantifying anharmonic interactions is particularly important for fast ion conductors, since the lower migration barriers in these ion conductors have been shown to be correlated with softer lattices, 43, 44 which typically allows for higher degrees of anharmonic interactions. 97 In a recent study, using AIMD simulations and Raman polarization-orientation measurements, Brenner et al. 98 attributed the large broadening of vibrational peaks in AgI superionic conductor to the strong anharmonic interaction of Ag + ion and the rigid iodine lattice. Our chosen Ge-substituted Li3PO4 compound is not a superionic conductor, which can be the reason for the low degrees of anharmonic contributions that we calculated for ion conduction from the MDMA method. Overall, our analysis provides a quantitative way to measure the effect of anharmonicity on the solid-state diffusion which has already been recognized in earlier works, but little in the way of quantitative evidence has been provided. 99,100 This is with the exception of a few compounds, such as CuCrSe2 101 and AgCrSe2 102 where recently, inelastic neutron scattering measurements have revealed a softening of acoustic modes which is associated with the superionic transitions in these materials.
Ultimately, more studies are needed to compare the inelastic interactions observed for our Ge-substituted Li3PO4 compound with other ion conductors with more polarizable anions where the anharmonicity are expected to be more important. 82, 103

Conclusions
We have introduced two modal decomposition formalisms that determine the individual normal mode contributions to ion diffusion. Our results for a model structure of Ge-substituted Li3PO4 show that a small group of modes (<10% of the vibrational modes in the system) are responsible for >87% of the diffusion. Of greatest significance is that when we externally excited a targeted subset of the highest contributing modes, the diffusion increased to the value associated with the effective temperature of the targeted modes even if when the lattice is kept at lower temperature of 400 K. This observation is intriguing, since a high diffusivity can be achieved while the bulk temperature of the material remains at a much lower temperature.
While the selected Ge-substituted Li3PO4 composition in this study comprises a rather slow ion conductor, we believe that the presented rigorous theoretical frameworks in this report can be used to study the fundamental concept of phonon-ion interaction, which is important in the field of superionics. [104][105][106][107] Recent reports have so far only linked the idea that phonons contribute to ion transport, [33][34][35][36][37][38][39][40][41][42][43][44] but this is the first study showing that specific mode excitations can be used to influence diffusion in solids, which is influential to the field of lattice dynamics affecting ion transport and will most likely lead to more work on optoionics. [108][109][110] and the following expression for the mode amplitude,

Calculation of vibrational modes at different local equilibrium configurations
The different modes of vibration at different local equilibrium configurations can be understood and modeled mathematically by expanding the potential energy surface via a Taylor expansion 17, 113 about a given equilibrium configuration as, where E represents the potential energy of all the interacting atoms, 0 E is the energy at the minimum, i u , j u , and k u are the displacements of atoms i , j , and k from their equilibrium position at that particular equilibrium configuration, and  ,  and  are linear, harmonic and cubic force constants. In such an expression, the term including the linear force constants   i  diminishes at a stable equilibrium configuration because the atomic forces decay to zero. The second order terms   ij  are then the largest non-zero terms, and most often the higher order terms, which affect the phonon-phonon interactions, can be safely neglected when we are simply interested in determining the mode shapes and frequencies. We can then solve the equations of motion for the entire system, via a super-cell lattice dynamics (SCLD) calculation. The output of such calculation is then a set of eigenvalues that describe the frequencies of each normal mode and a set of corresponding eigenvectors that describe the direction and magnitude of displacement that each atom in the system will experience when it participates in a specific normal mode (i.e. the mode shape). The modes are then specifically a result of the local equilibrium configuration, and in principle will differ for each distinct configuration.

Derivation of the velocity rescaling formula
After the equilibrium configuration is determined, we excite the modes based on the following formulation.
To change the temperature of the normal mode of vibration n to a desired tempearture d T , we first calculate the modal velocity coordinate of that mode using Eq. 5. Using Eq. 6, n Q  is then used to calculate the modal contributions to atomic velocity i v coming from the normal mode of vibration n   , i n v (Eq. 7 in the manuscript).
Knowing that, for a classical system, the temperature of the normal mode n can be calculated from 17 during the MD simulation, we modify the atomic velocities by first subtracting the effect of the modal components , i n v from i v , and then adding the components that result in the desired modal temperature d T for mode n .
which simplifies to,

Modal contributions to ionic conductivity
Ionic conductivity is proportional to diffusivity through the following relation, 82,83  

Diffusivity Calculations
Total mean squared displacement method (TMSD) is utilized to measure the diffusivity of Li atoms in our Gesubstituted Li3PO4. This method is well explained in existing literature, 86,87 In short, following the Einstein relation, TMSD calculates the diffusivity of the hopping particles based on the slope of their mean square displacment (MSD) data over time interval t  as, For three dimensional systems 3 d  , offset D is the linear shift from the origin, and for tracer diffusivity, (11) where N is the number of hopping atoms/ions n the system, Statistical uncertainty, due to insufficient phase space averaging, has been reduced by considering 5 independent ensembles. 88 The same exact simulation parameters were also utilized for modal excitation simulations with the only difference being that during the simulations the temperature of the desired normal modes in the system was changed to the desired temperature, as exaplined previously in the methods section. For calculating the conductivity diffusion coefficients used in the calculaiton of Haven ratios in Table S3, Eq. (10) is used in combination with the following formula, which is the MSD of the center of mass (COM) of Li + carriers, 84 The curves shown on Fig. S9 suggest that the 10 ns simulation time that we had for this structure was suffcient to get convergence in our Dσ simulations for Haven ratio calculations (Table S3).

Calculations of Attempt Frequency, RDF, Amplitude of Vibration, and Jump Rate
Attempt frequency (Fig. S3) is obtained from the weighted average of the partial phonon density of states of Li + ions, which is calculated from power spectrum of atomic velocities , , x y z v  extracted every 5 time steps (5fs), The total power spctrum is defined by total x y z P P P P    . RDF (Fig. S7) (Table S2) is found following the procedure explained by Klerk et al., 96 where the number of jumps are counted in a known period of simulation time. To detect the occurnace of a jump, the LAMMPS 79 code was modified so it follows which interstitial site is occupied at each time step. A jump happens when the occuiped interstitital site changes, which is detected and counted during the simulation. 96 The reasons for choosing empirical interatomic potentials and LAMMPS package for NEB calculations and

MD simulations
As mentioned in the manuscript, empirical interatomic potential with Pedone's parameters 69 are used for NEB calculations and MD simulations using LAMMPS package 79 . The main reasons that the NEB and MD results presented in this manuscript were obtained from classical interatomic potentials and LAMMPS package are: 1. The MD-based modal decomposition approach presented in the manuscript is based on the conductivity diffusion formula that is based on the fluctuations-dissipation theorem. 84,85 Previous studies have noted that compared to the mean-square displacement approach, 86, 87 the fluctuation dissipation approach 85 used in this study requires longer simulation times for convergence, which would be unachievable using AIMD approaches. 116 Finding a modal decomposition method based on direct mean square displacement (MSD) formulation 86 for diffusivity can help with faster convergence of the modal contributions. However, modal decomposition of the direct MSD formulation of diffusivity seems impossible because of the dependence of the atomic displacement on the reference position at the beginning of the simulation.
2. As explained in previous sections, to examine the effect of the selective modal excitation on the diffusivity, the atomic velocities in the system need to be modified, following the modal atomic velocity rescaling formula (Eqns. 6 and 7 in the ESI). Such a modification to the atomic velocities is achieved by modifying the LAMMPS velocity rescaling subroutine. Modifying the LAMMPS code is more straightforawrd than the other ab-initio codes. This when accompanied by the fact that longer simulation times are needed for modal analysis of lower temperatures such as 400K made us to utilize the LAMMPS package based on empirical interatomic potential.
3. In addition, as explained in the manuscript, 619 NEB and 31 supercell lattice dynamics were performed.
For better analysis of the NEB results to capture the modal contributions to the ion hop, for each of the NEB calculations 17 images were used along the MEP which also adds to the computational expense.
Even if we are using a 193-atoms unit cell, the extent of these calculations and simulations makes ab initio approaches almost impossible. However, nothing restricts the introduced methodologies in this manuscript to be applied to first-principles methods Possible approach for the experimental evaluation of the reported observations A possible approach for experimentally evaluating the reported observations in this study is by illuminating the sample by THz photons [19][20][21][22][23][24][25] or THz electric fields. 26 Figure S5. By changing the x-axis for each of the panels in Fig. S4 from frequency to percent number of modes, and averaging over all the modal contributions to the 619 hops in our system, which are shown in Fig. S4, and sorting the contributions based on their magnitude, it can be seen that more than 87% of lattice energy in the system during the hop comes from less than 10% of the modes. A value of IPR ∼ 1 corresponds to a highly localized mode, and a value of IPR ∼ 0 (more accurately 1 a N where a N is the number of atoms in the system) corresponds to a fully extended (i.e., delocalized) mode. 119         Assessing the thermalization power and excitation rate needed for the modal excitation simulations The power needed to keep the lattice temperature at 400K is shown in Table S5. It can be seen that the power needed for thermalization for our 192 atoms system is impractically large. However, this issue is an artifact of the small size of our simulation cell, and the power deposited for modal excitation does not scale with the system size, because the energy is deposited in a local manner. What matters is that the highly contributing modes are excited and this amount of power does not need to be dumped in such a small volume. For example, exciting modes to 700K is key to getting a diffusivity value comparable to the one that is obtained at a 700K bulk temperature, but the high power density is not necessary, because the energy input is local and thus need not scale linearly with the system size.
This issue can be understood better if we consider again the simulations that were done by exciting the highest single contributing mode in the system. When we increase the energy of the highest singly contributing mode, the energy is deposited in a very localized region in space, however as we increase the size of the system the required power density decreases accordingly (see Table S5). This decrease in power density can be justified by noting that if we increase the number of atoms we correspondingly decrease the size of the eigenvectors on the modes.
However, localized modes remain with large eigenvectors in the local regions. In fact, in larger simulation cells, the majority of the highly contributing modes to the local ion hop become the localized modes, because of their large eigenvectors of vibration. Therefore, by virtue of exciting localized modes, the power density reduces as we increase the system size. However, additional simulations at larger system sizes have shown that the contribution of excited localized modes to the atom hopping remains virtually unchanged even in larger simulation cells, as is evidenced by calculating the ion jump (hop) rate. As Table S5 shows if we excite the highest singly contributing mode to the atom hop in the system, the jump rate stays almost constant as we increase the system size, yet the power density decreases with increasing size. Therefore, a high power density need not be necessary to realize this. Table S5 also shows how low this power density would be if such a small amount of power were extrapolated to a macroscopic sized system. The power density, which in essence would manifest itself like heat generation is only on the order of 10 W/m 3 . To put this in perspective, consider the heat generation rate that occurs in 12 guage copper wires commonly used in residential wiring, which typically sees currents ranging from 1-10 Amps (limit ~ 40A). For such wires, which do not heat up appreciably, and are able to easily cool themselves via natural convection and radiation to the surroundings, currents ranging from 1-10 Amps (e.g., a 100 W light bulb vs. a toaster), the corresponding heat generation rates would be ~1.58-158 kW/m3. Thus, one could afford to increase the power density of targeted mode excitation by 2-4 orders of magnitude from the extrapolated estimation based on the power inputs to the small supercells studied herein. Heat capacity of Li3PO4 (J/g-K) 1.25 The modal excitations in the performed simulations in this manuscript were done every 5 time steps (every 5 fs).
If the excitations rate is decreased from every 20 fs to lower values, such as 100 fs, the increase in hopping rate and diffusivity will diminish. In particular, every 5 fs excitation, cannot be sustained with the existing THz laser technologies with their typical short duty cycles and a cycle repetition period of around milliseconds. 120 However, continuous excitations can still be achieved using continuous wave lasers and indirectly coupling their photons with the highest contributing modes in the system, and thus sustaining the excitations.