Correlated terahertz phonon–ion interactions control ion conduction in a solid electrolyte

Abstract

Ionic conduction in solids that exceeds 1 mS cm⁻¹ is predicted to involve coupled phonon–ion interactions in the crystal lattice. Here, we use theory and experiment to measure the possible contribution of coupled phonon–ion hopping modes which enhance Li⁺ migration in Li_0.5La_0.5TiO₃ (LLTO). The ab initio calculations predict that the targeted excitation of individual TiO₆ rocking modes greatly increases the Li⁺ jump rate as compared to the excitation of vibrational modes associated with heating. Experimentally, coherently driving TiO₆ rocking modes via terahertz (THz) illumination leads to a ten-fold decrease in the differential impedance compared to the excitation of acoustic and optical phonons. Additionally, we differentiate the ultrafast responses of LLTO due to ultrafast heating and THz-range vibrations using laser-driven spectroscopy (LUIS), finding a unique long-lived response for the THz-range excitation. These findings provide new insights into coupled ion migration mechanisms, indicating the important role of THz-range coupled phonon–ion hopping modes in enabling fast ion conduction at room temperature.

---

New concepts

Solid-state electrolytes (SSEs) with superionic conductivity (>1 mS cm⁻¹) are required to improve the safety and cycle life of state-of-the-art battery technology. However, the question of how the lattice dynamically interacts with mobile ions has remained unclear for decades. The extent to which lattice-ion coupling allows or hinders transport is an open question and, to this point, has not been directly probed experimentally. In this work, we develop an experimental method to measure phonon–ion coupling strengths for different phonon modes by exciting specific modes and measuring the corresponding change in ion conduction. Comparison to theory allows the degree of correlated phonon–ion motion to be determined. We find that targeted excitation of TiO₆ rocking modes reduces resistance by an order of magnitude compared to excitation of acoustic phonons in Li_0.5La_0.5TiO₃. The results indicate that structures which enable cooperative phonon–ion processes and have a large population of coupled phonon–ion hopping modes in the low-energy, THz regime have the potential to enhance ion conduction significantly. The developed methodology can probe dynamical process in any system in which ions conduct, such as interfacial battery processes, polymers, fuel cells, and artificial or natural membranes.

Introduction

Solid-state ion conduction has applications in Li-ion batteries,¹ biological membranes,² solid oxide fuel cells,³ and more. Ongoing research in this field has led to the continual improvement of superionic conductors that can enable energy dense solid-state electrochemical systems while maintaining fast conductive properties akin to liquids.^4–8 Realizing an all-solid-state battery with a superionic solid-state electrolyte and a Li metal anode can meet the necessary energy and power requirements for developing long-range electric vehicles, and improve safety.⁹

The fundamental mechanism of ion conduction (σ_ion) relies on the Arrhenius form equation as shown in eqn (1), where σ₀ is the Arrhenius pre-factor, T is temperature, E_a is the activation energy, and k_B is the Boltzmann constant.¹⁰

Researchers have sought out different strategies to increase σ_ion by lowering the E_a of the ion migration pathway by targeting crystal structures with corner sharing frameworks⁶ and connected migration channels,^9,11 employing highly polarizable polyanions in the sub-lattice,^12,13 introducing disorder to cation sites,^6,8 and leveraging polyanion–motion-mediated ion transport.^14–16 These principles have led to the discovery of several superionic conductors through aliovalent/isovalent substitution strategies¹⁷ and a general search for ideal structural families.^8,13,14,16 Li₁₀GeP₂S₁₂,^4,6,7 Li₇P₃S₁₁,^6,9,14,15 Li-argyrodites,^6,8,12 and Li-oxyhalides^18–20 have very high conductivity values, but the exact mechanisms supporting superionic conductivity remain debated.^20,21 In particular, the coupling of phonon modes with ion hopping has been suggested to play a significant role in the ionic conduction mechanism.^13,22–26

To understand phonon–ion coupling, we can expand the pre-factor σ₀ in eqn (1). As shown in eqn (2), σ_ion depends on the entropy of migration (ΔS_m) and enthalpy of migration (ΔH_m).¹³

The value of ΔS_m is determined by collective phonon–ion interactions²⁷ which are predicted to correlate with and enhance ionic conductivity.^{13–17,22–25,27,28} The value of ΔH_m is experimentally measured as E_a (determined by eqn (1)).

Several experimental efforts have attempted to explore the phonon–ion couplings predicted to relate to ΔS_m. Nuclear magnetic resonance (NMR) techniques have been used to probe phonon–ion coupling in Li-argyrodites^22,29 and other conductive Li-sulfides, -nitrides, and -oxides.³⁰ Neutron scattering techniques have been applied to Na-sulfides and -selenides,^16,23,31,32 as well as Na-orthophosphates and -sulphates,^16,33 specifically by characterizing the reorientation of anions as quasi-elastic structure factors.³³

Another approach to measuring how phonon–ion coupling enhances Li⁺ hopping is to drive relevant modes with light and measure the resulting change in conductivity. The terahertz (THz) region, in particular, can selectively target polyanion modes theorized to mediate ion transport.^28,34–37 Specifically, modes that involve (PS₄³⁻),^14,23,26 (SiS₄⁴⁻),¹⁴ (PO₄³⁻),³⁶ (YBr₆³⁻),³⁸ and (ErCl₆³⁻)³⁸ polyanion vibrations or rotations occur with vibrational frequencies at or below 10 THz^28,36—similar to the picosecond timescales of ion hopping.¹⁶ The same modes are theorized to be responsible for coupled hopping in LISICON-like conductors (Li₁₀SiP₂S₁₂,³⁹ Li₁₀GeP₂S₁₂,^4,28 Li₃PO₄,^25,36 Li₇P₃S₁₁^14,28), Li-argyrodites (Li₆PS₅Cl, Li₆PS₅Br),^12,22 and superionic Li-halides (Li₃YBr₆, Li₃ErCl₆).^28,38 Similarly, octahedral rotations are also reported to promote fast ionic transport in several perovskite-related phases,^24,40,41 making phonon-driven ion migration a promising approach for unveiling ion-coupled mechanisms that may lead to superionic conductivity.

While driving acoustic and optical modes to change electronic and magnetic properties is common in condensed matter physics,^42–44 very few studies have studied the effects of directly driving phonon–ion couplings to promote ion transport.^34,37 In one report, Poletayev et al. employed a THz pump to induce a Na⁺, K⁺, and Ag⁺ ion hop in β-Al₂O₃ and then measured the change in birefringence using an 800 nm probe. The THz Kerr Effect measured the anisotropy caused by the THz excitation, which they calculated to lead to a 20× increase in hop anisotropy compared to isotropic heating.³⁷ In addition, Gordiz et al. used nudged-elastic band (NEB) calculations and molecular dynamic (MD) simulations to predict that over 87% of the lattice energy required for the ion hop in Ge-substituted Li₃PO₄ came from <10% of the modes in the system.³⁶ By computationally exciting highly contributing modes, they enhanced the ionic conductivity by orders of magnitude without increasing the bulk temperature.

In this paper, we use ab initio calculations, MD simulations, and experiments that combine laser-driving and electrochemical impedance spectroscopy (EIS) to determine the role of octahedral rocking TiO₆ phonon modes which couple with the hopping of Li⁺ in LLTO. The bottleneck of Li⁺ migration in LLTO (Fig. S1) is comprised of 4 oxygen atoms formed by 4 corner sharing TiO₆ octahedra^40,41 that are reported to “rock” in place at low THz frequencies. These “rocking modes” are depicted in Fig. 1. To determine how well the TiO₆ modes couple with the hopping of Li⁺ ions, we first excite these modes with THz irradiation from 0.7 to 4.5 THz, or 67–428 µm (general schematic shown in Fig. 1). We compare the driving of targeted phonon modes predicted to couple with Li⁺ conduction to laser-induced heating, which is performed via the indirect excitation of acoustic phonon modes with NIR light (800 nm or 0.8 µm). Acoustic phonons can mediate ion conduction via increased random thermal vibrations.⁴⁵ Lastly, we excite relevant optical phonon modes within the MIR light range. The La–O mode (11.6 µm) is predicted to weaken if Li⁺ interactions with the TiO₆ octahedra are strong,⁴¹ so we test whether driving the La–O vibrations with MIR light (9.4–14.16 µm) can weaken Li⁺–TiO₆-coupling and modulate Li⁺ transport.

Fig. 1 Schematic of the laser-driven impedance method with LLTO. An excitation source (λ) excites the sample between two in-plane Au electrodes while EIS measurements are collected. Depending on the wavelength of the excitation source, coupled phonon-ion rocking modes, optical phonons, and acoustic phonons can be excited to influence the Li⁺ migration in LLTO.

We simultaneously use EIS to probe the change in impedance due to NIR–THz excitation, correlating these excitations with changes in the ion transport behavior. We quantify the role of the differing phonon couplings that mediate ion conduction by measuring shifts in impedance at each excitation wavelength and normalizing the measured enhancement by the power density of the incident laser beam. Additionally, we utilize laser-driven ultrafast impedance spectroscopy (LUIS)^46,47 to differentiate the ultrafast responses of the NIR and THz excitations. We separately quantify the impacts of laser-induced heating on the sample, finding that long-lasting thermalization effects do not account for the differences in the measured enhancement, indicating that the targeted photoinduced lattice dynamics predominantly affect ion migration. We note that the aim of this study is not to significantly reduce the sample's resistance using the laser. Instead, we use a differential measurement to compare the impact of phonon excitations with that of laser-induced heating, providing an initial proof-of-concept demonstration that vibrational modes play a significant role in enabling fast ion conduction.

In qualitative agreement with the ab initio calculations, our measurements indicate that selectively exciting the rocking modes in the THz range which enable phonon–ion coupled Li⁺ hopping leads to an order of magnitude decrease in the bulk resistance, or R_bulk, as compared to exciting acoustic or optical modes. Additionally, our LUIS measurements demonstrate a unique long-lived response time in the case of the THz excitation. Reminiscent of other metastable phases induced by light in inorganic–organic hybrid perovskites,^48–51 the enhancement that we measure is persistent and reversible, even though the recorded changes are averaged over the 250 Hz repetition rate for the THz source and 1 kHz for the generated NIR-MIR light. Our results provide initial evidence that laser driving phonon–ion coupled hopping modes can facilitate fast ion transport in a solid, suggesting that structures which have these modes in the THz range may make for fast ion conductors with ionic conductivities greater than 1 mS cm⁻¹.

Results

Using ab initio calculations, the contributions of phonon–ion coupled hopping modes to the energy needed by Li⁺ to reach its transition state (herein referred to as ion hopping) in LLTO are calculated based on the magnitude of the projection of the atomic displacement field obtained for the transition state on the eigenvectors of vibration for each specific phonon, as explained in detail in the Methods. The structural parameters for the calculations are obtained from quantitative Rietveld refinements of experimental synchrotron diffraction data (Fig. S2). To efficiently sample different possible Li⁺ hops and relevant contributing vibrational modes, 22 unique Li⁺ hops based on two different hopping mechanisms (single and concerted) across three degrees of La–Li disorder in the LLTO crystal are chosen and analyzed (Fig. S3).

Fig. 2a shows the normalized cumulative contribution of different phonon modes to the Li⁺ hopping process, here defined as the energy needed for a Li⁺ to reach its transition state. Fig. 2b shows the individual breakdown of the total phonon mode contributions that sum to Fig. 2a. Within the experimentally excited region of <4.5 THz, marked by the dashed line in Fig. 2a, there are numerous phonon modes which couple strongly with the hopping of Li⁺ ions. These modes are further separated based on their projection on different atomic species in Fig. S4. The modes below 4.5 THz contribute to 17.5% of all sampled hops, corroborating that low frequency phonon modes strongly couple with ion transport as studied in other theoretical reports.^15,28

Fig. 2 Calculated phonon contributions to Li⁺ hopping and representative LLTO phonon modes. (a) The cumulative contributions of phonon modes to the normalized energy required for ion hopping. The black dashed line marks the experimental limit for THz generation (4.5 THz). The cumulative effects of the highest 5% and 10% contributing modes, as well as the octahedral rocking modes are specifically shown. (b) The individualized contributions of phonon modes to the normalized energy required for ion hopping where ρ_Oh represents octahedral rocking modes. The data is normalized over all 22 Li⁺ hop investigations. (c)–(e) The eigenvectors of three vibrational modes show how the lattice displacements project along with the Li⁺ hop, under (c), a thermal mode (2.66 THz), (d) the highest contributing non-rocking mode (4.29 THz), and (e), the highest contributing rocking mode (3.46 THz). The red rectangle highlights the hopping Li⁺ which is being investigated. The enhancement provided by each of these modes is outlined in Table 1 and the hop to which these modes contribute is depicted in Fig. S5.1.

The significance of the phonon contribution is further supported by the increase in ion hopping rate determined using MD simulations (see Methods and Table S1). The effects of targeted excitation of a thermal mode, the highest contributing non-rocking mode, and the highest contributing rocking mode are depicted in Table 1, showing that a three order of magnitude increase in the jump rate is found through targeted excitation of a phonon–ion coupled rocking mode relative to exciting a random phonon mode. The modes in Table 1 are further examined in Fig. 2c–e where the random mode at 2.66 THz does not appear to align with the hopping path of the Li⁺ ion (Fig. 2c), which is highlighted by the red box, while the highly contributing modes, both rocking and non-rocking, appear to have crystal lattice vibrations aligned with the Li⁺ jump path (Fig. 2d and e), again highlighted in red. These figures show how increased alignment between the hopping Li⁺ and the excited phonons lead to an enhanced ion hopping rate and increased phonon–ion coupling. The overall hopping path that these modes couple with is depicted in Fig. S5.1. Similar hopping pathways are further examined in Table S1 and Fig. S5.2 and S5.3 indicating similar levels of enhancement due to the targeted excitations of coupled phonon–ion hopping modes which have a high degree of alignment between the phonon and the hopping ion. Due to the low frequency nature of modes below 4.5 THz, the inclusion of quantum effects changes their contribution by less than 1% (see Methods and Fig. S6), making the calculations relevant to the room temperature experiments.

Table 1 MD-simulated enhancement in the ion hopping rate of one of the 22 investigated hops by targeted excitation of selected modes which have frequencies in the experimentally accessible frequency range. The visualization of this hop, which is a single hop in the disordered LLTO structure, is provided in Fig. S5.1. Details on the calculation are provided in the Methods

		400 K (natural MD simulation)	700 K (natural MD simulation)	400 K (with mode excited to 700 K)
Jump rate (#jumps/ps)	Thermal mode (2.66 THz)	0.0011	1.45	0.0027
	Highest contributing non-rocking-mode (4.29 THz)			1.27
	Highest contributing rocking-mode (3.46 THz)			1.39

---

To experimentally test the role of phonons in coupling with Li⁺ conduction in LLTO, we first synthesize the material in accordance with literature⁵² and structurally characterize it using synchrotron X-ray diffraction, as described under Methods.

We then irradiate the LLTO using THz, NIR, and MIR light using an open electrochemical cell described in the Methods and described previously in literature.⁴⁷ We adopt an in-plane electrode geometry, as depicted in Fig. 1. This sample geometry allows for the irradiation of the sample but has irregular (i.e. nonplanar) surfaces with no accurately defined thickness, meaning the ionic conductivity of the laser-excited measurements cannot be calculated, requiring the use of R_bulk to study the photoinduced changes in ion migration. Accordingly, we quantify the laser-induced enhancement on ion conduction through the decrease in R_bulk, which we determine by fitting the Nyquist plots (obtained with EIS) to a R1 + R_bulk/Q2 circuit, shown in both plots of Fig. 3. The high frequency data points, which represent the bulk feature, between 32 MHz and 803 kHz, are partially affected by electronic induction which arises due to the wiring of the experimental setup not physical phenomena that we aim to probe in our samples. Consequently, we include only approximately 5 points in our fit. We are still confident, however, that we are probing the bulk region as EIS measurements on LLTO have long-reported low-impedance bulk and high-impedance grain features,⁵³ which we have replicated in Fig. S7.

Fig. 3 Nyquist plots of LLTO under varied light intensities. The zoomed-in Nyquist plots of the bulk ion migration feature measured over a total frequency range of 32 MHz to 1 Hz with a sinus amplitude of 50–100 mV to minimize electronic noise for (a), THz light excitation powers varied between 0 and 1.6 mW and (b), 800 nm light excitation varied between 0 and 25 mW. Because the absolute changes in impedance are difficult to discern, insets are shown for data points measured at 803 kHz, also marked by an arrow on the corresponding Nyquist plot. The average value of each data point under each light intensity condition across three samples and three trials each is plotted.

We first drive the THz rocking modes predicted to couple with Li⁺ hopping (see Fig. S8 and S9) with 0.24–1.6 mW of THz power focused to a 250 µm beam diameter with a 250 Hz repetition rate (Methods and Fig. S10). The bandwidth and field strength of the THz radiation are determined by the generation process, described in more detail in the Methods. When the 0.7–4.5 THz modes are driven, a linear decrease in the Z′ on the order of hundreds of Ωs is measured, as seen in Fig. 3a, scaling with increased power (Table S2). We next measure the enhancement induced by heating using NIR light. We excite the sample with 5–25 mW of 800 nm pulsed light focused into a 250 µm beam diameter at a 1 kHz repetition rate to populate acoustic phonon baths associated with random thermal vibrations, or heat³² (see Methods). The Nyquist plots for LLTO under 800 nm laser excitation are shown in Fig. 3b, with the replicable enhancements tabulated in Table S3. We additionally rule out any enhancement from structural changes in the material during and after excitation. We perform a damage test, as outlined in the Methods, and find that the Raman spectra acquired at the laser spot before and after excitation are identical (see Fig. S11).

Using the fitted Nyquist plots shown in Fig. 3, we compare the percent change in R_bulk across different powers and excitation wavelengths (Fig. S12). We find that incoherently heating the acoustic phonon bath with the 800 nm light induces a 0.10% decrease in R_bulk per mW of power. Comparatively, a 0.78% decrease in R_bulk per mW is measured when exciting phonon–ion coupled hopping modes with THz light, in line with the predictions from ab initio computational modeling which predict a significant enhancement associated with the excitement of single phonon modes. Fig. S12 not only demonstrates the greater enhancement induced by the THz-range coupled phonon–ion modes, but also rules out any nonlinear effects, given the linear correlation seen between the excitation enhancement and power for both wavelengths.

In Fig. 4, we normalize the enhancement by the average power density of the beam used to excite the sample. Normalization by average power density is utilized because it allows for comparison across different excitation wavelengths which have varied spot sizes and repetition rates, as mentioned previously and outlined in the Methods. With this normalization, we find a 4.8 × 10⁻⁴% change in R_bulk (mW cm⁻²)⁻¹ for the broadband THz excitation of phonon–ion coupled hopping modes, while laser-induced heating gives a calculated 5 × 10⁻⁵% change in R_bulk (mW cm⁻²)⁻¹. We also attempt to excite optical phonon modes with lower energy MIR 1100 and 1500 nm light, finding no detectable change (Fig. 4a). The results are consistent with our theoretical calculations where the optical modes contribute significantly less to ion hopping above 10 THz (∼30 µm) shown in Fig. 2b, which possibly justifies why no detectable changes are observed.

Fig. 4 Enhancement in bulk ion migration due to illumination across near-IR (NIR), mid-IR (MIR), and THz light. (a) Enhancement ratios (cyan squares) reported in %ΔR_bulk (mW cm⁻²)⁻¹ × 10⁻⁴, calculated from exciting LLTO between NIR – THz light. The width of the horizontal orange bar represents the spectral bandwidth of the excitation pulse, with the frequency range probed in the THz regime also noted. The NIR enhancement corresponds to incoherent heating of the acoustic phonon bath. The MIR excitation corresponds to coherent driving of optical phonon modes, which was inconclusive for the same excitation density as the THz modes. The THz light coherently drives highly contributing phonon–ion coupled hopping modes, showing the largest relative enhancement. (b)–(d) Absorption spectra of LLTO (magenta) in the NIR, MIR, and THz regions obtained using UV-vis spectroscopy, FTIR, and THz transmission, respectively. The NIR absorption plot was measured on a Shimadzu SolidSpec-37000iDUV spectrophotometer, the MIR was measured on the Bruker Vertex 80, and the THz absorption was measured on a custom THz time-domain spectroscopy set up which is presented in more detail in the Methods. The use of different spectrometers necessitates the differing units on each plot.

The small magnitude of the normalized enhancements in Fig. 4a satisfies our aim of utilizing the differential nature of the measurement to compare the impact of different excitations without damaging the sample (Fig. S11) and while avoiding any nonlinear effects (Fig. S12). Still, we acknowledge that the absorption depths of the 800 nm and THz excitations (7 and 6 µm, respectively) are roughly 1% of the total sample thickness and thus may diminish the measured enhancement. We chose this thickness due to the brittle nature of LLTO causing thinner pellets of this size to crack during measurements. Still, we aim to lessen this spatial mismatch by using an in-plane geometry, which is expected to primarily apply the EIS probe near the sample surface between the sputtered electrodes, where the laser most strongly interacts with the sample (Fig. S13). A graphical depiction of the laser-driven EIS experiment can be found in Fig. S14. The entire sample holder is depicted in detail in previous work.⁴⁶

Next, we spectroscopically explore what types of modes explain the differential enhancement measured in Fig. 4a using UV-VIS, FTIR, THz transmission, and Raman measurements. In Fig. 4b, one can see there is low absorption in the NIR used for incoherent thermal laser heating. This indicates why higher powers of the incident laser are needed in measurements (Fig. S12). In Fig. 4c, one can see the vibrations of Ti–O and La–O vibrational modes.⁴¹ Evidently, these modes do not affect Li⁺ hopping within the experimental error bars. Finally, in Fig. 4d, we see high absorption coefficients, on the order of 1000 cm⁻¹ for THz light excitation. Given that there are many modes in this range, as indicated by Fig. 2a and b, we cannot discern individual modes via THz absorption measurements. However, in Fig. 5, we identify Raman-active modes that exist within the THz excitation energy range. The first of these is a shoulder peak <100 cm⁻¹, marked with an asterisk, which is hypothesized to be related to octahedral tilts and small oxygen shifts,⁵⁴ very similar to the key modes which we believe couple strongly with Li⁺ conduction. Additionally, the E_g symmetry mode is associated with Ti vibrations, which affect the TiO₆ polyhedra,⁴¹ similar to the computationally identified TiO₆ rocking modes. While not all Raman-active phonon modes can be directly excited by THz irradiation, it has been proposed that THz excitation of IR-active modes can decay into available Raman modes via coupled pathways, which we hypothesize to be relevant in this study.⁵⁵ We further discuss the spectroscopic evidence for highly contributing modes in LLTO, particularly the THz absorption measurements in Fig. 4d, at the beginning of the Discussion.

Fig. 5 Raman spectrum of LLTO. The Raman spectrum of LLTO, demonstrating available phonon modes which can be accessed in the excitation range, from 0.7–4.5 THz, as highlighted in the figure.

Finally, we utilize LUIS to differentiate the response of LLTO when exciting THz-range phonon–ion coupled hopping modes, as compared to 800 nm NIR light, which is meant to represent the effects of laser-induced heating. In LUIS, as described in previous work,^46,47 a single frequency AC field with similar amplitude to the laser-driven EIS experiments, is used to initiate ion hopping on site-to-site timescales of tens of picoseconds. Then, a laser excitation is used to induce a change in the crystal lattice, in this case specifical vibrational modes with THz light and laser-induced heating with 800 nm light. The corresponding impedance response is then measured with a high-speed oscilloscope. This technique is described in further detail in the Methods, and is depicted in Fig. S15. We provide a summary comparing the laser excitation and AC probing field parameters for the laser-driven EIS and LUIS measurements in Table S4.

For these experiments, we utilize a 19 GHz, −5 dBm (178 mV sinus amplitude) alternating current (AC) field to maximize the signal-to-noise ratio. The probing depth of the 19 GHz AC field in LLTO is 1.15 cm,⁵⁶ allowing for the full 450 µm-thick sample to be probed before reaching the metal short, which reflects the microwave signal. For these measurements, a thinner sample can be used as compared to the laser-driven EIS experiments. For the LUIS experiments, the pin used to transmit the AC signal is 350 µm in diameter, fully sampling the 6–7 µm laser absorption depth for the THz and 800 nm excitations, allowing for through-plane measurements, which differs from the in-plane laser-driven EIS measurements. The signal is further processed using averaging and amplitude demodulation, as seen in Fig. S16a, b and S17a, b.

We fit the signal responses to a Gaussian response function convolved with an exponential decay, detailed in the Methods. The width of the Gaussian is reported as the rise time, limited by the sampling rate of the oscilloscope, which is 7.8 ps. Further subsampling is possible when recording sub-nanosecond responses, leading to the differing sample rates seen in Fig. 6a and b and rise times reported in Table S5. We also report the exponential decays with a 95% confidence interval and attribute these responses to dynamics induced by each laser excitation. The deconvolution of these fits is seen in Fig. S16c and S17c.

Fig. 6 Laser-driven ultrafast impedance response of LLTO. The LUIS response of LLTO due to (a) broadband THz excitation, with a wavelength range of 54.5–599 µm (0.5–5.5 THz) at a field strength of 744 kV cm⁻¹ and (b) an 800 nm 8 mJ per pulse excitation. The red and blue signals represent the averaged amplitude demodulation response due to a 19 GHz, −5 dBm reflected voltage signal with a root-mean-squared envelope function applied to the response signal. The inset in (b) highlights the sub-nanosecond response of the 800 nm excitation. Both signals were fit to a Gaussian convolved with an exponential decay, described in detail in the text, and denoted by the solid black lines. The dashed black lines at 0 mV demonstrates the baseline.

We find, in Fig. 6a, that the THz excitation induces in a 1.03 ± 0.03 ns decay time, which we hypothesize is due to increased ion hopping due to the excitation of strongly absorbing phonon–ion THz-range rocking modes (Fig. 4d), as well as potential decay pathways associated with the phonon–ion coupled hopping modes that are theorized to enhance ion hopping (Fig. 2a and b). As seen in Fig. 6b, the 800 nm excitation induces a much sharper response with a decay time of 0.15 ± 0.01 ns, which we attribute to the excitation of an acoustic phonon bath associated with heating.⁴⁷ We summarize the fitted rise and decay times in Table S5. In short, the differing ultrafast responses of the two excitations demonstrate that the THz excitation induces a unique, long-lived response that is distinct from that due to laser-induced heating and is responsible for the outsized impact of phonon–ion coupled hopping modes in the THz energy range, as detailed throughout this work.

Discussion

The relative decrease in R_bulk with THz light excitation of rocking modes is almost ten-fold that seen with 800 nm excitation, which is predicted to incoherently heat the sample. The significant difference in enhancement emphasizes the potential for theorized phonon–ion coupled hopping modes to promote ion migration at higher rates than uncorrelated thermal, acoustic phonon-induced hops.

At the macroscopic scales probed by the several hundred micron diameter spot size of the THz absorption data shown in Fig. 4d, the absorption coefficient of LLTO above 2–3 THz exceeds 1000 cm⁻¹, values in excess of those seen in high dipole moment liquids such as water and acetonitrile.⁵⁷ The multi-octave bandwidth of the THz pulse enables both resonant electric dipole-allowed transitions and sum frequency excitation via the polarizability operator, that is Raman active phonon modes.⁵⁸ The absolute values and smooth nature of the cross sections suggest that the transition dipoles are large, and that local variability in the Li/La ratio impacts the phonon mode frequencies, leading to a broad distribution of modes in this energy range, which cannot be cleanly decoupled in the THz absorption spectra (Fig. 4d). This is further corroborated by Raman measurements of Li_3xLa_2/3−xTiO₃ (0.04 ≤ x ≤ 0.15), where the local distribution of Li⁺ and La³⁺ alters the spectroscopic signature of modes associated with the TiO₆ octahedra.⁴¹ We hypothesize that the absorption coefficients at the level of Fig. 4d require significant coupling of octahedral rocking and ion displacement, likely with potential energy surfaces that are quite flat for small displacements. This suggests the coupling of relatively higher frequency phonons with Li⁺ modes at relatively lower frequencies, as proposed by Xu et al.²⁸ Furthermore, recent work indicates that in LLTO, modes such as the key TiO₆ rocking modes we investigate here can distort Li sites, and the corresponding distortion can lead to a mechanism termed “ion-lattice-ion concerted diffusion”.⁵⁹ The activation of modes that give rise to this concerted mechanism upon THz excitation may be responsible for the long decay time found in the LUIS experiments and the strong relative enhancement observed in laser-driven EIS. Determining the fraction of THz modes that couple with actual ion hopping/conductivity requires careful comparison to the theory of the hopping mechanism and action spectroscopy to assess the role of THz excitation on ion hopping quantitively.

Having performed these measurements, we propose that the enhancement due to THz laser excitation can be attributed to phonon–ion couplings which enhance Li⁺ hopping. However, the generation of electronic carriers and heating effects must also be considered. Any electronic contribution toward the total measured impedance, including the excitation of electronic carriers due to THz fields, are unlikely in LLTO, because (1) there is no change to the measured open circuit voltage during the EIS experiments, (2) the THz radiation is <0.025 eV⁶⁰ which lies sufficiently below the 3.3 eV band gap of LLTO,⁴⁷ and (3) no nonlinear effects are observed (Fig. S12). Furthermore, we have previously demonstrated that, even in the case of a band gap excitation, electronic contributions to changes in the conductivity are negligible.⁴⁷ Therefore, the THz fields are unlikely to generate electronic carriers or have sufficient energies to promote enhanced electronic conductivity in LLTO. In the presented work, the measured changes in impedance should be attributed to changes in ion hopping. Additionally, no dependence on excitation near the contacts on the LLTO was found at the metal-material junction, ruling out the generation of carriers on the metal contacts.

Deconvoluting thermal enhancements of ion hopping from those associated with exciting targeted lattice-driven dynamics is non-trivial. A complete Nyquist plot that captures the bulk semi-circular feature requires tens of seconds for each measurement, preventing direct probing of the “peak” temperature change induced by the laser which occurs in less than 100 ns (Fig. S18 and S19). On the milliseconds-to-seconds timescales of the EIS measurements, we expect both excitations to have contributions from the decay of the initial phonon excitation into heat. However, if heating is the only mechanism by which the laser enhances ion hopping for both the 800 nm and THz light excitations, we would expect these excitations to cause similar levels of enhancement when normalized, given they have similar penetration depths and laser-sample interaction volumes. Instead, the THz light leads to an order of magnitude larger enhancement compared to the 800 nm light when normalized by power density (Fig. 4a). Additionally, the timescale of the ion hopping response due to the THz excitation is an order of magnitude longer, as measured by LUIS (Fig. 6). Such a discrepancy requires further investigation into the differences in laser-indued heating and longer-term photo-modulation effects between both excitation energies.

To more accurately compare the enhancements in ion migration between the laser-induced heating and THz-range phonon–ion coupled hopping modes and separate contributions due to heating and photo-modulated ion hopping at milliseconds-or-longer timescales, we adopt a finite-difference time domain (FDTD) method, detailed in the Methods. The method aims to capture resistor heating, while accounting for wavelength-dependent experimental parameters, outputting a single metric of heating normalized by the average laser power (Fig. S20). We can then separate the effects of heating from the different excited phonon baths at each wavelength.

We first identify the absorption depths of both 800 nm and THz light. Then, using the model, we produce cross-sectional heating maps of the sample, which show differences in the temperature gradient between the THz and 800 nm light (Fig. S21 and S22). Given the challenges in reconciling the region probed electrochemically with that heated by the laser, we aim to quantitatively determine an upper bound on the heating relative to the milliseconds or longer timescales probed by EIS. We note that ultrafast laser pulses will induce peak temperatures well above the metric we define, but this heating is unlikely to be captured by the EIS measurement. As such, we instead capture the average maximum temperature in the sample at all time points during and after the final pulse, once the sample has reached a baseline thermal steady state, which occurs within 100 ms (Fig. S23 and S24). Using this metric, we find that 800 nm light induces a temperature change of 0.03 K mW⁻¹. In contrast, THz light induces a change of 0.02 K mW⁻¹. We confirm the accuracy of the method by calibrating it against measurements with an IR thermal gun for the 800 nm excitation (Methods and Fig. S25), and as described in previous work.⁴⁷

In Table 2, we use the thermal enhancement factors in K mW⁻¹ to differentiate the computed thermal effects from the measured total enhancement, the slope in Fig. S12. To ascertain the ion migration enhancements associated with resistor heating, we utilize eqn (3).

Table 2 Deconvolution of thermal-induced ion enhancement from the total measured enhancement of ion migration

Excitation source	Calculated thermal effects (%ΔR/mW)	Measured total mnhancement (%ΔR/mW)
800 nm	0.07 ± 0.01	0.10 ± 0.01
0.7–4.5 THz	0.04 ± 0.01	0.78 ± 0.04

---

We measure the decrease in sample resistance per degree Kelvin, the slope in Fig. S26, and multiply this factor by the computed wavelength-dependent heating factors, discussed above, thus making up the right side of eqn (3). The value provided by this equation is a quantitative measure of the upper limit for resistor-like thermal-induced heating effects on timescales relevant to EIS, in units of %ΔR/mW, as depicted in the second column of Table 2. Finally, we compare these computed thermal effects to the slopes of Fig. S12, the third column in Table 2.

From Table 2, we can determine the fraction of total enhancement attributable to resistor-like heating relative to the total measured enhancement, finding that 60–80% of the enhancement at 800 nm is associated with thermal heating. We note that the approximately 30% discrepancy likely indicates that laser and resistor heating affect ion hopping differently. However, the model still allows us to attribute most of the 800 nm excitation enhancement to purely thermal effects. Given this, we believe that this wavelength measures the enhancement associated with populating the acoustic phonon bath, or heating the sample.⁴⁷ For the THz excitation, we can again find the fraction of the total enhancement which is associated with resistor-like heating effects. We find that approximately 5% of the THz-induced enhancement is attributable to heating. In comparison, the remaining 95% should be attributed to additional photo-modulation effects associated with the population of phonon–ion coupled hopping modes, in qualitative agreement with the computational work presented in this paper. By separating the thermal contributions from the total measured enhancement, we demonstrate that we are indeed probing the photo-modulation of ion hopping rather than simply the effects of laser-induced heating.

We further aim to account for grain boundary scattering by varying the thermal conductivity of the sample in the FDTD model. Grain boundaries can suppress phonon scattering and thus reduce the thermal conductivity of the sample.⁶¹ We find that altering thermal conductivity does not meaningfully change the wavelength-dependent enhancement, as shown in Tables S6 and S7. This model does not explicitly account for wavelength-dependent phonon scattering and instead generalizes these effects in the remainder of the total enhancement, beyond the calculated resistor-like heating contribution. This remainder not only includes any directly excited dynamics in the crystal lattice, but also the associated decay pathways, as seen in Fig. 6a and b. Based on the previously presented computational results, we hypothesize that the modes seen in Fig. 5 and the highly absorbing nature of the phonons found in Fig. 4d are part of the THz-specific decay pathway and affect the TiO₆ octahedra in ways that are theorized to couple with ion hopping.

In addition to ruling out electronic carriers and heat, we also aim to address discrepancies between the theoretical and experimental data. The results in Table 1, which demonstrate a three order of magnitude increase in ionic conduction, represent limiting cases where energy can be deposited into three well-defined comparative cases: the excitation of the full phonon bath, individual highly contributing modes, and non-contributing modes. In comparison, the THz excitation excites all phonon modes which can absorb THz light which does not match the computationally excited cases. Rather, we must use THz absorption and Raman spectroscopy to gain insight into what our light excites. As shown in Fig. 4d, we observe strong but smooth absorption across the THz spectrum, preventing the identification of individual modal components via deconvolution. This smooth absorption is due to local variability in Li/La ratios throughout LLTO. However, when investigating what modes exist in the low-energy THz regime, we can consider the Raman spectrum (Fig. 5), where we identify modes that may couple strongly to hopping Li⁺-ions and represent plausible decay pathways.⁵⁵ The optical excitation of these modes leads to a factor of ten increase in the differential enhancement as compared to laser-induced heating, which quantitatively differs from Table 1 by two orders of magnitude. We believe that this disagreement can mainly be attributed to the differences in the computational and experimental excitations discussed above.

Still, we consider alternative explanations below. One can rule out the mismatch in EIS probe and laser absorption depth because both the 800 nm and THz light have similar absorption depths, so the relative enhancement ratios should be maintained, given that the laser pumps the same sample volumes. We also note the mismatch between the timescales of the computational and experimental results. While the computed enhancement is reported in hops/ps, the 803 kHz EIS probe used in Fig. 4a probes processes occurring on microsecond timescales, meaning that we may not perfectly probe bulk properties, but instead may be probing local structural effects. Impurities inherent in LLTO⁶² and domain boundaries⁶³ may be captured by EIS but are not captured computationally, even when considering both ordered and disordered structures (Fig. S3). As discussed, we attribute quantitative discrepancies to differences between computational and optical excitations, as well as to differences between computational and experimental probes, while emphasizing the qualitative agreement across the various portions of this work. This, when combined with the long-lived response measured in the LUIS experiments, paints an intriguing picture regarding the important role of THz-range vibrations on enhancing ionic conduction.

Conclusion

In conclusion, we compare the relative role of acoustic, optical, and THz phonon modes in ion migration through laser-frequency-selective perturbation of EIS and LUIS measurements. By directly driving highly contributing phonon–ion coupled hopping modes, which are most likely rocking modes, in the 0.7–4.5 THz region, a ten-fold enhancement in ion migration is measured relative to incoherent heating. Using an FDTD model, we attribute most of this enhancement to laser-driven photo-modulation effects and find that these excitations lead to notably different responses on ultrafast timescales. The anomalous contribution of rocking modes to ion hopping in the <4.5 THz region is qualitatively consistent with ab initio calculations, validating the use of laser-driven methods to probe the mechanism of phonon–ion-coupled interactions.

Furthermore, the findings herein can aid in the design of future solid-state ion conductors by emphasizing the design of structures with a large population of phonon modes that couple to ion hopping and are inherently active at room temperature (<6 THz). More specifically, our study emphasizes the role of complex phonon–ion interactions, indicating that, at least in perovskite structures, engineering structures with highly active THz rocking modes may increase ionic conductivity by coupling lattice dynamics with the hopping ion. Additionally, THz-range phonon–ion coupling has been found to be relevant in polymer electrolytes, potentially allowing for new design rules which emphasize key vibrational modes which can enhance ion hopping rates.^64,65 Similarly, in biological systems, THz light may enable the modulation and control of ion channels.⁶⁶ This work provides an initial foray into the direct experimental measurement of how THz-range vibrational modes affect ion hopping. We emphasize this exploratory motivation, which still offers intriguing results, demonstrating the outsized role of low-energy phonon–ion coupled hopping modes in enhancing ionic conduction in a solid-state ion conductor.

Methods

General DFT setup

Li⁺ ion hop simulations in this study were conducted under the spin-polarized density functional theory (DFT) approximation using the Vienna ab initio simulation package (VASP)⁶⁷ within the projector augmented-wave approach. Perdew–Burke–Ernzerhof (PBE)⁶⁸ generalized-gradient approximation (GGA) functionals, electronic KPOINTs equal to 4 × 4 × 4, an electronic energy cutoff of 500 eV, and Gaussian smearing with sigma of 0.05 eV were employed for the calculations. The valence electrons for lanthanum, lithium, titanium, and oxygen atoms were generated in the 5s² 5p⁶ 5d¹ 6s², 1s² 2s¹, 3p⁶ 3d² 4s², and 2s² 2p⁴ configurations, respectively.

Details for NEB calculations

The minimum energy pathway for Li⁺ ion hops in the lattice were calculated using the nudged elastic band (NEB) methodology following the formulations by Henkelman et al.⁶⁹ to ensure the correct determination of the saddle point as implemented in VASP. Three middle images were considered for the NEB calculations. Structures at the end of the hops were relaxed to lower the forces and variations in energy below 10⁻⁴ eV Å⁻¹ and 10⁻⁵ eV, respectively. Obtaining these end-point relaxed structures are critical in a successful NEB calculation to obtain meaningful migration barriers (MB). Because of the challenges in the convergence of DFT relaxation simulations, which is mainly caused by the multitude of relaxation sites available in the lattice of highly Li⁺ conductive compounds,⁷⁰ different optimization algorithms were tested for the convergence of both local relaxed structures (corresponding to the beginning and end of the hops) and NEB calculations (reaction coordinates). The structure relaxation and NEB calculation were initially tested for convergence using conjugate-gradient energy-based optimization method (IBRION = 2). If no convergence was achieved using energy-based optimizers, force-based optimizers were then employed. Two force-based optimizers⁷¹ were identified to be effective in helping the convergence: (i) QM = Quick-Min (IOPT = 3, IBRION = 3, POTIM = 0), and (ii) FIRE = Fast Inertial Relaxation Engine (IOPT = 7, IBRION = 3, POTIM = 0). We reported the relaxed structure, the migration barrier (MB), the reaction coordinates using whichever method that converged. Afterwards, we utilized the reaction coordinates for the phonon modal analysis of the respective ion hop.

Choosing the initial structure for computations

The same experimental structure identified from synchrotron X-ray diffraction of La_0.5Li_0.5TiO₃ in this study is used for the ab initio NEB and phononic calculations. For ease of calculations and better visualizations, the experimentally determined R3̄c structure (a = b = 5.4711 Å, c = 13.4040 Å, α = β = 90°, γ = 120°) is transformed to its pseudo-cubic (perovskite-looking) counterpart⁷² (a = b = c = 7.7378 Å, α = β = γ = 90°) through the rotation/transformation of the lattice based on the eqn (4) where, α = 0°, γ = −45°, and β = −35° correspond to yaw, roll, and pitch angles, respectively. We considered 40 atoms in this pseudo-cubic unit cell, and we used it for the rest of our simulations. After the transformation, only internal atomic coordinates were relaxed during the calculations.

Maximizing the sampling of Li⁺ hop in the LLTO lattice through different choices of lattice disorder, Li⁺ hopping mechanisms, and hopping pathways

The experimentally determined coordinates for La and Li contain partial occupancies and need to be substituted with atoms in order to perform atomistic simulations. It is important to maintain the exact experimental stoichiometry (La_0.5Li_0.5TiO₃) during this process. Depending on which La and Li partial occupancy sites are selected for atom placement in the simulation, locally different structures of LLTO can be formed. Experimental⁷³ and computational^74–76 studies have shown that the various structures of LLTO can be classified into three categories based on the degree of La|Li disorder in the lattice (Li|La orderings): (i) fully ordered structure, (ii) partially ordered structure, and (iii) fully disordered structure. However, due to the limited size of the supercell used in our simulations (40 atoms), it is difficult to distinguish between the partially ordered and fully disordered structures. Therefore, to maximize our sampling of Li⁺ ion hops in the LLTO lattice within the constraints of our supercell size, we included three different La|Li orderings in our analysis of Li⁺ ion hops, as illustrated in Fig. S3. These chosen structures included one structure with ordered Li|La layers, as well as two structures with disordered Li|La layers.

Since our calculations do not involve any hopping for La atoms, the position of La atoms is the primary factor that distinguishes these three structures with different La|Li orderings. Therefore, these structures can be regarded as parent structures, under which we have defined several substructures by shuffling Li⁺ ions among unoccupied sites in the LLTO lattice. Although our simulation cell sizes may not be large enough to fully capture the exact differences between these parent structures with different La|Li orderings, we believe that incorporating these different ordered and disordered parent structures, along with their substructures and the different Li⁺ hopping mechanisms included in these substructures, has increased our sampling of different local chemical environments in the LLTO lattice in our calculations. The schematics of the parent structures with different Li|La ordering layers, the number of Li⁺ hops calculated for different substructures under these parent structures obtained by shuffling Li⁺ ions in the lattice or along different hopping pathways, and the different considered hopping mechanisms are illustrated in Fig. S3.

In addition to the different degrees of disorder in the LLTO lattice, we also considered two different hopping mechanisms: (i) single and (ii) concerted ion hops. For the single ion hop mechanism, the hop occurs through a vacancy. However, there are no vacancy sites in the stoichiometry considered in this study (i.e., the nominal number of vacancies in the stoichiometric La_2/3−x/3Li_x□_1/3−2x/3TiO₃ for x = 1/2 is zero). To allow the occurrence of the vacancy-mediated mechanism, a La cage must be doubly occupied by two Li⁺ ions. This results in an unoccupied La cage that was previously occupied, which is now available for other Li⁺ ions to hop into. Thus, the choice of the doubly occupied La cage and the resulting hopping pathway will be the determining factor for performing our NEB calculations. We combined different (i) orderings of Li|La, (ii) choices of which La cage to doubly occupy, and (iii) hopping mechanisms to choose 22 different Li⁺ hopping events to maximize our sampling efficiency of the LLTO lattice. The distribution of the different hopping mechanisms (single and concerted) across the chosen structures with different degrees of disorder is illustrated in Fig. S3.

DFPT setup for modal calculations

To obtain the vibrational modes (phonons) for each of the structures considered in this study, we first calculated the Hessian matrix using density–functional–perturbation theory (DFPT) with the VASP software and the Phonopy package⁷⁷ to capture the matrix. Next, we calculated the phonon frequencies and eigenvectors using a custom lattice-dynamics Python code.⁷⁸ It was important to obtain these values for the exact supercell containing 40 atoms, since the NEB calculations were also performed on cells of this size.

Phonon modal analysis of ion hop

To identify the contribution of different phonons to the Li⁺ hop in the lattice, we used the methodology recently proposed based on combining the lattice dynamics and NEB calculations.³⁶ Following this method, after the minimum energy pathway for the Li⁺ hop is identified using the NEB calculation, we project the displacement field obtained from the transition state (saddle point) on the eigenvectors of vibration belonging to the structure at the beginning of the hop. The magnitude of projection shows the degree to which each normal mode is contributing to that displacement field, which can be quantified using the following expression,^79,80where u_i is the displacement of atom i from its equilibrium position in the configuration at the hopping origin; e_i,n is the eigenvector for mode n, assigning the direction and displacement magnitude of atom i obtained from the lattice dynamics calculation^79,81,82 for the structure at the hopping origin; * denotes the complex conjugate operator; and m_i is the mass of atom i. Q_n is the modal displacement coordinate, the square of which is proportional to the mode potential energy E_n according to the following equation,⁸⁰The total energy of the system E is equal to the summation over all modal energy values ; here, E_n is the contribution by mode n to the potential energy of the displaced lattice during the ion migration, which can be interpreted as the contribution by mode n to the ion hop along its migration pathway. In addition, ω_n is the frequency of vibration of mode n. E_n can be divided by E to calculate the normalized contribution of that phonon to the Li⁺ hop in that specific ion hop event. Moreover, by subsequently normalizing E_n through division by the total number of ion hopping events considered in the study, the resultant normalized E_n can be regarded as the normalized contribution of mode n to the Li⁺ hopping events in the LLTO compound under investigation. The code used to identify the modal contributions of Li⁺ hop in the LLTO lattice, following the outlined methodology, was written in Python and is available as a GitHub repository.⁷⁸

Algorithm for identifying the rocking modes

In this study, we identified the octahedral modes of vibration by quantifying the degree of circulation the eigenvectors of vibration impose on the octahedron units (TiO₆) in the lattice. To this end, the eigenvectors of vibration for the 6 oxygen atoms belonging to the TiO₆ octahedron unit were utilized to quantify the circulation of the unit according to the definition of circulation (Γ), e.g., , where Γⁱ_n is the circulation belonging to the octahedron unit i (in our case we have 8 octahedron units in our simulation cells, thus i = {1, 2, 3,…,8}) coming from mode of vibration n. In addition, is the path vector along which the circulation is being evaluated, and v⃑ is the eigenvector of vibration belonging to the specific atoms along the direction. In practicality, to simplify the evaluation of the circulation formula above, we considered the perfect representation of a perovskite structure (e.g., with untitled octahedron units) and evaluated Γⁱ_n along the circulation pathways normal to x, y, and z directions, which results in Γ^i,α_n {α = x, y, z}. By averaging the obtained values of Γ^i,α_n over α, a representative circulation value is found for octahedron unit number i (Γⁱ_n). By averaging again over all the octahedron units in the simulation cell, i, a representative value of circulation is obtained for all the octahedron units in the simulation cell (Γ_n). If this value is larger than a specified criterion, the mode of vibration number n is identified as an octahedral mode of vibration. Example MATLAB code block below demonstrates this calculation:

% mondenum = 3*Natoms % number of modes of vibration in the system

for n = 1: modenum

 vortTi = zeros(8, 3);

 % TiO₆ unit 1 (atomids listed below)

 ix1 = 23;

 ix2 = 26;

 iy1 = 21;

 iy2 = 27;

 iz1 = 34;

 iz2 = 36;

 vortTi(1,1) = (ev((iy2-1)*3 + 3, n) − ev((iy1-1)*3 + 3, n)) − (ev((iz2-1)*3 + 2, n) − ev((iz1-1)*3 + 2, n));

 vortTi(1,2) = (ev((iz2-1)*3 + 1, n) − ev((iz1-1)*3 + 1, n)) − (ev((ix2-1)*3 + 3, n) − ev((ix1-3)*3 + 3, n));

 vortTi(1,3) = (ev((ix2-1)*3 + 2, n) − ev((ix1-1)*3 + 2, n)) − (ev((iy2-1)*3 + 1, n) − ev((iy1-1)*3 + 1, n));

⁝

 for ii = 1 : 8

  vortTi_mag(ii) = (vortTi(ii,1)^2 + vortTi(ii,2)^2 + vortTi(ii,3)^2)^.5;

 end

 vortTi_ave(n) = mean(vortTi_mag);

 vortTi_mag_allmodes(n,:) = vortTi_mag;

end

tagOctRot(1:modenum, i) = 0;

for n = 1: modenum

 vorticity_criteria = 0.25;

 if vortTi_ave(n)^2 > vorticity_criteria

  tagOctRot(n) = 1;

 end

end

Calculating the imparted energy on different atomic species by different modes of vibration

To calculate the energy that each mode of vibration imposes on the hopping Li⁺, we leverage the following identity equation existing from lattice dynamics theory,⁸⁰, where v⃑ⁱ_n is the eigenvector of vibration belonging to mode n, imposing vibrations on atom i. The vibrations of atoms in the crystal are governed by the eigenvectors of vibration, and the imparted vibrations on an atom will add up to give an energy equal to k_BT to that atom. In this view, the v⃑ⁱ_n·v⃑ⁱ_n term in the summation above can be considered to be proportional to k_BT. This methodology can be extended to the eigenvectors of vibration on the hopping Li⁺, where they can be projected along the hopping direction of the Li⁺ to give the energy imparted on Li⁺ along the hopping pathway, e.g., . Similarly, the term can be considered to be proportional to k_BT.

Quantifying the change in O-4 bottleneck area due to phonon modes in the structure

To evaluate the impact of each phonon on the O-4 bottleneck area during ion hopping, we applied a displacement to the lattice along the eigenvectors associated with that phonon. This displacement can be considered as an induced excitation of the lattice. The extent of this excitation can be calculated for each vibrational mode using the following equations.

The potential energy for an oscillator can be written as arising from the sum of harmonic and anharmonic contributions,⁸⁰

〈H_pot〉 = 〈H_harmonic〉 + 〈H_anharmonic〉 (7)

where 〈…〉 represents the ensemble average.⁸³ The harmonic portion has been shown to dominate the potential energy.⁸⁴ Thus, we will also leverage the harmonic potential energy to evaluate the mode amplitude that is needed to calculate the degree of displacement field in the structure. Based on the equipartition theorem, the average harmonic energy of a classical oscillator is equal to⁸⁵where k_B is the Boltzmann constant and T is the temperature of the system. The potential energy of a harmonic oscillator can also be calculated from the normal mode amplitude analysis viaeqn (9).⁸⁰where Q_n is the modal displacement coordinate (eqn (5)), and ω_n is the frequency of the eigen mode. It should be noted that calculating the harmonic energy of an eigen mode using the knowledge of force constant matrix and the respective atomic displacements for an eigen mode is equal to the approach based on the normal mode amplitude analysis (eqn (10)). The modal displacement coordinate can be explicitly obtained from eqn (5), which is rewritten here.⁸⁰where u_i is the displacement of atom i from its equilibrium position in the configuration at the hopping origin; e_i,n is the eigenvector for mode n, assigning the direction and displacement magnitude of atom i obtained from the lattice dynamics calculation^79,81,82 for the structure at the hopping origin; * denotes the complex conjugate operator; and m_i is the mass of atom i. By combining eqn (8) and (9) we then have,

Q_n²ω² = k_BT (11)

To determine the distribution of the harmonic energy for an eigen mode amongst the atoms in the system, one can envision a state of the system whereby all of the atoms in the system are displaced from equilibrium in the direction of their respective eigen vectors from mode n (i.e., e_i,n). In this view, the attributed displacement to an atom would be equal to,

u_i = α_ne_i,n (12)

where α_n is a scaling factor that associates a certain degree of displacement with the mode's amplitude at a given temperature. The exact value of the scaling factor can then be calculated from the combination of eqn (5), (11) and (12). Replacing for the atomic displacement in eqn (5) with the definition in eqn (12), we would have,which by the substitution of , yields a simpler form,Using eqn (11) we then have and by incorporating eqn (14), we can calculate the scaling factor as,By using the obtained α_n in eqn (15), the displacement field is obtained. The O-4 bottleneck area will then be calculated for both the equilibrium (no displacement) lattice and the lattice for which the atomic positions have been displaced according to eqn (15) to assess the capability of each phonon in increasing the O-4 bottleneck area against the ion hop. It should be noted that the four oxygen atoms forming the bottleneck area are not necessarily on the same plane in space, making it impossible to calculate the planar area connecting them. Therefore, four smaller triangular areas inside the bottleneck area are defined (as shown in Fig. S1), and their areas are calculated and averaged to approximate the bottleneck area in this situation, as demonstrated in Fig. S1. The percent change between the O-4 bottleneck area in the equilibrium lattice and the one in the displaced lattice is then calculated and reported as the change in the ability of that specific phonon in altering the bottleneck area against the hop of Li⁺ ion.

Quantum correcting the calculated modal contributions to the ion hop

To include the quantum effects in our reported modal contributions to the Li⁺ hop in LLTO, we multiply our obtained modal contribution by mode n (E_n, obtained from eqn (6)) with frequency ω by a correction coefficient that is based on the Bose–Einstein distribution function at temperature T (f_Q(ω, T)) based on the following formula to obtain the quantum corrected contribution of that mode n to the hop in the lattice (E^qc_n),⁷⁹

E^qc_n(ω) = f_Q(ω, T)E_n(ω) (16)

The function f_Q represents the ratio of quantum to classical specific heat for mode n. It restricts the contributions of the high frequency modes at low temperatures and modulates the modal contributions obtained from the NEB-based ion hop simulations. The quantum expression of volumetric specific heat, based on Bose–Einstein statistics, is given byand the classical volumetric specific heat is given by , where V is the volume of the supercell under study. Thus, the quantum heat capacity correction factor which is the ratio of C_q and C_c is

Enhancing ion hopping rate by targeted excitation of phonons in MD simulations

To assess if any increase in the hopping rate of Li⁺ in LLTO lattice can happen by exciting the identified contributing phonons to the ion hop, we chose three hops from the 22 investigated hops in this study. We then excited the highly contributing phonons in their respective structures. We performed the following three tests: (i) exciting the top contributing (non-rocking) mode overall, (ii) exciting the top contributing rocking mode, and (iii) exciting a random mode.

To perform the molecular dynamics simulations, we utilized classical molecular dynamics simulations based on the Morse and Buckingham form potentials as proposed by Cormack and coworkers^86,87 and successfully utilized by Chen and Du⁸⁸ to investigate the diffusion of Li⁺ diffusion in the LLTO lattice. The simulations were performed using the same DFT simulation cells using the large-scale atomic/molecular massively parallel simulator (LAMMPS)⁸⁹ package. The modes that were excited were the same modes that were detected from the DFT calculations. The goal of these experiments was to only assess if any extra energy (excitation) given to these highly contributing modes will result in any enhancement in the hopping rate of the ion in the lattice.

First, the structure was relaxed under the Nosé–Hoover canonical ensemble for 100 ps at T = 400 K; then, the structure was simulated until a Li⁺ hops in the structure. Once the hop gets detected through continuous examination of equilibrium MD positions during the MD simulation,³⁶ the simulation gets interrupted and restarted with a new initial condition for atomic velocities, which results in exploring a new region of the phase space. This process gets repeated until a total of 50 ns has been simulated. Then, the total number of detected ion hops will be compared based on the group of phonons that were excited during the simulations. The time step was chosen to be 1 fs.

Because ion diffusivity increases with temperature, we ensured that the bulk lattice temperature remained constant at a low temperature (T = 400 K), so that any change in the diffusivity could be attributed only to the excitation of the targeted modes and not to the bulk temperature. To do so, we kept the total kinetic energy of the system constant via a velocity-rescaling scheme, in which, the addition of energy to the intended 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 T_d, we modified the atomic velocities in the system according to the following formula:

where v_i is the velocity of atom i, and Q̇_n is the modal velocity coordinate of mode n defined by,⁸²

In this study, the above-mentioned rescaling procedure was applied every five time steps (every 5 fs). Although a velocity-perturbation scheme has been used in this study, other methods for modal excitation, such as atomic position perturbation⁹⁰ or simultaneous perturbation of atomic positions and velocities,⁹¹ have also been employed in previous studies, particularly those with the goal of investigating the phonon–phonon interactions in MD simulations. However, the chosen method for modal excitation should not change the reported observations in this study because modal excitation (increasing the energy of a mode) can be achieved either by increasing the potential (position perturbation) or the kinetic (velocity perturbation) energy of the mode; the choice of which should not matter.

The results for the excitations are shown in Table S1. By targeted excitation of the detected contributing phonons in the LLTO structure (including the contributing rocking modes), the hopping rate of the Li⁺ increases to the same order magnitude of increasing the hopping rate by increasing the temperature, but this time, by keeping the lattice at the base (T = 400 K) temperature.

Development of finite-difference time-domain methods to deconvolute thermally induced impacts on ionic conduction

The heat accumulation and subsequent change in temperature of the samples was simulated using a finite-difference time-domain method which has previously been developed to reproduce the results of a three-dimensional two-temperature model.^92,93 This method treats the laser pulse as instantaneous given that the timesteps in the simulation are significantly longer than the length of the pulse. For this paper, 100 ns timesteps were taken for a total simulation time of 1 s. A subsection of the sample with a 1.5 mm radius and 0.3 mm depth was simulated. The total simulation time was chosen because it allowed the sample to reach a clear steady state in its baseline temperature. Additionally, the timesteps and total volume chosen allowed heat to effectively diffuse through the sample, while also balancing the large computational cost of simulating for long periods of time.

For this model, the thermal conductivity and heat capacity of LLTO were estimated to be similar to other oxide solid-state electrolytes.^94,95 Given the measured porosity of our synthesized LLTO (∼0.264), the thermal conductivity was then adjusted using previous relations developed for other oxide materials.⁹⁶ While these properties fluctuate greatly with temperature, given that most laser-induced heating effects from ultrafast lasers dissipate in far less than 1 ms, the time between the 800 nm laser pulses, these variations are likely to be very short-lived and have minimal impact on the final temperature values.

The penetration depth of the 800 nm light was estimated to be 7 µm using the following equation:⁹⁷

where c is the speed of light, κ is the imaginary part of the refractive index, and ω is the angular frequency of light. For this calculation, the imaginary part of the refractive index was taken from a previous study on the optical properties of LLTO,⁹⁸ while the absorption depth of THz light was found to be 6 µm via broadband THz absorption spectroscopy on LLTO.

The simulations produce a micron-resolution heating map, as seen in Fig. S21 and S22, showing the spatial variation of temperature both on the surface and into the depth of the sample. It is difficult to ascertain the exact thermal effects of the laser since there is a fast rise and decay of the lattice temperature caused by each pulse, pictured for the final pulse in Fig. S18 and S19. However, as seen in Fig. S25, the peak temperature in the sample was found to be between 0.4 and 0.8 for a 20 mw 800 nm beam, and ∼0.5 from simulation. It was not possible to take a similar measurement for the THz light due to instrument constraints. As such, using this value as a calibration, we take the average of the sample's peak temperature, during and after the final pulse, as an approximate temperature in the region we are probing, allowing us to separate thermal effects due to laser illumination. As seen in Fig. S23 and S24, by the time 1 s has passed, a baseline steady state has been reached, making this average temperature a reasonable approximation for the temperature in the sample. This characteristic temperature increase is found to scale linearly with the average laser power as seen in Fig. S20. The slopes of the lines in Fig. S20 are utilized to find the calculated thermal effects presented in Table 2.

Synthesis

Li_0.5La_0.5TiO₃ (LLTO) was synthesized according to literature.⁹⁹ A stoichiometric amount of La₂O₃, Li₂CO₃, and TiO₂ were mixed in an agate mortar and pressed into pellets under 100 MPa of pressure. The pellets were placed on a bed of sacrificial powder and calcined at 800 °C for 4 h then 1200 °C for 12 h at a ramp rate of 1 °C min⁻¹. The following material was routinely characterized using a Rigaku X-ray diffractometer with CuKα radiation and scanned from 10 to 70 2θ at a scan rate of 0.4° per second. The following x-ray diffraction pattern was fitted using a Rietveld refinement with the GSAS II software. High resolution synchrotron powder diffraction data was collected using beamline 11-BM at the Advanced Photon Source (APS), Argonne National Laboratory.

Sample preparation

The resulting powder was pressed into a pellet with a diameter of 9 mm and a thickness of 0.7 mm under 2 5 minute cycles of 2 tons of pressure, yielding a 75–76% pellet density. The pellet was subsequently annealed at 1100 °C for 6 h at a ramp rate of 2 °C min⁻¹ over a bed of sacrificial powder. A 1.6 mm thick strip of Au was sputtered onto one side of the pellet with a 1 mm gap in the center for the light source to excite the LLTO to eliminate possible effects on the impedance from illuminating the Au. Cu and stainless steel contacts were used to contact the sputtered Au to the 1260A Solartron impedance analyzer for EIS measurements. An open cell set-up was employed to allow THz irradiation to excite the sample without compromising average power.

The conductivity was then measured by an AC impedance method over a frequency range of 32 MHz to 1 Hz with an applied 50–100 mV sinusoidal amplitude with and without the excitation source. We confirm the linearity of this 100 mV sinus amplitude by comparing results for R_bulk at 20, 50, and 100 mV in a swagelok cell with Au-sputtered electrodes via EIS measurements between 3 MHz and 1 Hz. The results are seen in Fig. S27 and Table S8. A copper mesh Faraday cage was custom-made with a 1.4 mm copper wire spacing to reduce the noise in EIS measurements due to electromagnetic interference from the environment.

Heating cell set up

A heating cell was custom made to heat the sample between 298–333 K. The cell temperature was controlled using a TC-48-20 OEM temperature controller, 12 V power supply, and corresponding TC-48-20 OEM software and thermocouple. The heating cell was placed inside a Faraday cage for all experiments to reduce noise from electromagnetic interference. This set up was used to acquire the data in Fig. S26 and is depicted in previous work.⁴⁶

Raman characterization

Raman measurements were performed on a Horiba Instruments PLUS Raman spectrometer with a 532 nm laser focused on uncompressed LLTO powder on a glass microscope slide. The signal was acquired by averaging 200 acquisitions lasting 1 s each with a 50 µm slit and 500 µm hole and 10% average power.

Damage testing

Pressed pellets of LLTO were characterized using Raman, as described above. Then, we focused 25 mW of the 800 nm laser at a spot on the sample for 30 minutes on and 10 minutes off for 6 rounds. We then characterized this spot using Raman after the 800 nm excitation to ascertain if there was any damage. We note that the peaks in Fig. S11 are much sharper than in Fig. 5 because pressed pellets increase the intensity of otherwise hard-to-resolve peaks in the Raman spectra.¹⁰⁰

Terahertz Kerr effect (TKE) setup

The full THz set-up is pictured in Fig. S28. A 1 kHz, 38 fs, 13 mJ pulse centered at 800 nm is produced via regenerative amplification in a Ti:Sapphire laser (Coherent Legend Elite Duo). The output beam was split by a 25 : 75 beam splitter, with the more intense ∼10 mJ beam used to pump an optical parametric amplifier (OPA, TOPAS, Light Conversion Inc.), generating a 2.1 W, 1440 nm output signal beam. The signal beam travelled through a mechanical delay stage and a chopper set at 250 Hz to reduce the average power to 1.1 W to not burn the DAST (4-N,N-dimethylamino-4′-N′-methyl-stilbazolium tosylate) crystal, which was purchased from Swiss Terahertz. The µW output of the DAST crystal, was focused through three OAPs into a ∼250 µm beam spot, providing a peak field intensity in the 100 kV cm⁻¹ range. The chamber was purged under inert gas (N₂) to minimize water absorption and to increase the field strength of the THz source. The THz field was characterized using electro-optical detection in 30 µm GaP.¹⁰¹ The time-domain trace from this generation set-up is found in Fig. S10. We note that there is a water signal in the time-domain trace due to the ∼10% relative humidity inside of the purge box. The 10 dB bandwidth of the generated light is 0.7–4.5 THz,¹⁰² as noted in Fig. 4 and 5 in the main manuscript.

For the LUIS measurements, the same beam path to Fig. S28 was used, however, the chopper was set to 500 Hz, and the nonlinear organic crystal used was 4-N,N-dimethylamino-4′-N′-methyl-stilbazolium 2,4,6-trimethylbenzenesulfonate (DSTMS). The DSTMS crystal was pumped with 820 mW of 1300 nm light (410 mW after chopping) generated by an OPA, and the signal was charcterized using electro-optic sampling, with 100 µm GaP. The time trace and frequency character of the signal are seen in Fig. S29. The peak field strength of the THz field for these measurements was 744 kV cm⁻¹.

IR power to THz field strength

Terahertz pulse energies were measured using a calibrated pyroelectric THz joulemeter (Genctec EO SDX 1211). The subsequent THz field strength values were used to calculate an average power for comparison to the non-resonant heating.

THz average power calculation

Avg. power (W) = J × Rep. rate (Hz) (23)

Using a pyroelectric THz joulemeter, a measured voltage was used to calculate the pulse energy of the THz field using eqn (22) where 614 000 V J⁻¹ is the conversion factor according to the manufacturer specifications. The average power of the laser was then calculated using eqn (23) where J is the pulse energy and the repetition rate of the pulses after passing the chopper is 250 Hz. The conversion between THz field strength and laser power (in mW) is seen in Table S9.

Near-IR (NIR) to Mid-IR (MIR) light set-up

The NIR and MIR light sources were generated using a regenerative Ti:Sapphire laser amplifier operating at 1 kHz. The 800 nm output was sent through an optical parametric amplifier for the NIR light and a difference frequency generation unit for the MIR light.

THz absorption spectrometer

The Terahertz (THz) absorption spectrum shown in Fig. 5d was acquired using a home-built transmission-geometry THz time domain spectrometer (THz-TDS) using optical light (Coherent Astrella 800 nm, 35 fs pulse width, 1 kHz rep-rate Ti:Sapphire regenerative amplifier) to illuminate a spintronic THz emitter (TeraSpinTech T-Spin2) to generate p-polarized THz radiation, which is passed through a wire grid polarizer to assure polarization purity. Generated THz light is subsequently collimated and refocused down onto the sample position using a train of off-axis parabolic mirrors (OAPs). In the sample position, we mount the sample to a translation stage to alternate the acquisition of reference (dry nitrogen) spectra and sample spectra by moving the sample in and out of the sample position to account for any systematic drift through the data-acquisition process. After passing through the sample position, THz light is once again collimated and focused onto a GaP THz detector crystal (Swiss THz, 100 µm thick active cut optically contacted onto 1 mm thick inactive cut) where it is profiled by electro-optic sampling. In this process, the arrival time of a linearly polarized optical probe (800 nm, 35 fs) is scanned using a mechanical delay stage. As the arrival time of the two pulses is scanned, the polarization of the optical probe is rotated due to the transient birefringence induced in the detector crystal by the THz pulse. This polarization shift is then detected in a pair of balanced photodiodes.¹⁰³ The combination of the THz emitter and detector detailed here results in a calculated field strength of ∼15 kV cm⁻¹, and a 10 dB bandwidth of 0.5–7.3 THz.¹⁰² We present the absorption spectrum in the range of 0.7–4.5 THz in Fig. 5d to match the bandwidth of the higher field strength radiation that was used to drive the rocking modes in LLTO. Finally, the THz-TDS spectrometer is fully enclosed in a dry nitrogen purge box and is monitored with a hygrometer to assure that all spectra are taken in a dry environment.

Interpretation of THz-TDS

Upon capturing the time trace of the THz pulses for our reference (dry nitrogen) and sample (pressed LLTO:PTFE pellet), we take the fast Fourier transform (FFT) to get the complex frequency spectra E_ref(ω) and E_s(ω), respectively. By dividing the sample spectrum by the reference, we arrive at the transfer function H(ω) defined below:¹⁰⁴where l is the sample thickness, c is the speed of light, and n_s(ω) and k_s(ω) are the real and imaginary refractive index of the sample, respectively. By taking the argument and logarithm of the transfer function, we arrive at the following two equations:From these, we then extract the optical constants and related absorption coefficient directly from the transfer function as shown below:¹⁰⁴

n_s(ω) = 1 − c∠H(ω) (27)

By comparing the pure absorbance spectra of ball-milled PTFE with that of the dilute LLTO in PTFE mixed pellets and assuming linear contributions to absorption by concentration, it was possible to extract the pure LLTO absorbance spectra.

THz-TDS sample preparation

To prepare the pellet samples for THz-transmission data, LLTO powder was loaded into a ball-mill grinder alongside the THz transparent polymer PTFE (Teflon) in a known concentration and allowed to mix for 20 minutes. PTFE was chosen as the dilutant polymer due to its low porosity in pressed pellets, high tensile strength, and low propensity to capping defects at high pressure.¹⁰⁵ Following the mixing process, the PTFE and LLTO grains shared an average sub-wavelength, size of ∼1 µm, measured with scanning electron microscopy. Following mixing in the ball-mill, pellets were pressed for 1 minute under 1 ton of pressure. To account for scattering effects, multiple pellets of varying concentrations of LLTO in PTFE, ranging from 0% to 7.5%, were analyzed.

LUIS

In LUIS, a 19 GHz, −5 dBm AC carrier is utilized to initiate ion hopping under a fast-oscillating electric field. Then, an ultrafast laser, as discussed in the TKE Setup and NIR-MIR light sections of the Methods, is used to excite the sample, inducing a change in ion hopping. Due to the impedance mismatch at the LLTO sample, there is a modulated reflection of the AC carrier signal due to the photoexcitation. The GHz signal is carried via coaxial cables and routed to the sample and the oscilloscope, using a hybrid coupler. Averaging, amplitude demodulation, and a root-mean-square (RMS) envelope are used to extract the signal (Fig. S16a, b and S17a, b). This modulated response is captured using a real-time oscilloscope, which samples at a rate of 128 gigasamples per second (every 7.8 ps), and, at shorter timescales, applies further subsampling. This setup is depicted in Fig. S15 and the sample holder is detailed in previous work.⁴⁶ We fit the response to a Gaussian convolved with an exponential decay, where the width of the Gaussian represents our overall response time and the exponential decay characterizes the laser-induced dynamics that we measure. This is outlined below:where C is the fitted amplitude, t₀ is the time at which the laser excites the sample, σ is the system response and width of the Gaussian, and τ_d represents the characteristic decay time of the laser-induced dynamics.

Author contributions

S. K. C. and K. A. S. supervised the project. S. K. C. and K. A. S. conceived the work. S. K. C., K. A. S., and K. H. P. designed the experiments. K. H. P. developed and performed the laser-driven EIS experiments. K. H. P., A. K. L., and N. A. S. synthesized and characterized the material. J. M. M. and H. L. generated the THz radiation, characterized the THz spectrum, and supported the integration of THz with the impedance technique. N. A. S. and A. K. L, performed the THz absorption and ultrafast impedance measurements. K. G. and D. V. performed the atomistic simulations. N. A. S. performed the heating model simulations. K. H. P., K. G., and N. A. S., wrote the manuscript. S. K. C., K. A. S., G. B., A. K. L., D. V., A. H., and Y. S. H. edited the manuscript.

Conflicts of interest

The authors declare no competing interests.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. The SI containts details on the X-ray diffraction, DFT and MD simulations, laser-driven EIS experiments, THz light generation, the LUIS experiments, and the FDTD heating model. See DOI: https://doi.org/10.1039/d5mh01990g.

The code is available under the MIT license and can be found at https://github.com/kgordiz/modalNEB.

Acknowledgements

KHP acknowledges financial support from the National Science Foundation, Air Force Office of Science & Research under award number FA9550-21-1-0022, and David & Lucile Packard Foundation. AKL acknowledges support by the California Institute of Technology through the Beckman-Gray Graduate Fellowship and the Natural Sciences and Engineering Research Council of Canada, funding reference number 599279. NAS acknowledges financial support by the National Defense Science and Engineering Graduate Fellowship. NAS and AKL acknowledge additional personnel support from the U.S. Army Research Office (ARO), grant number W911NF-23-1-0001. We also acknowledge support for the equipment used in this work by the Air Force Office of Science Research (AFOSR), DURIP grant number FA9550-23-1-0197. FT-IR data was collected at the Laser Resource Center in the Beckman Institute of the California Institute of Technology with the assistance of Dr Jay Winkler. Solid-state UV-VIS data was collected at the Earle M. Jorgenson Laboratory of the California Institute of Technology with the assistance of Dr Weilai Yu. We thank Prof. David Hsieh and Dr Omar Mehio for assistance in using the difference frequency generation unit for IR light. We thank Ricardo Zarazua at the Chemistry and Chemical Engineering Machine Shop for machining the custom heating cell in this work. We thank Dr Zachery W. B. Iton for assistance with the characterization of LLTO. We thank Jadon Bienz for assistance with analysis of the XRD data including multi-phase Rietveld refinements of LLTO. We thank Jax Dallas for generating the THz light for the THz absorption and ultrafast impedance measurements. Use of the Advanced Photon Source at Argonne National Laboratory was supported by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH1135. The finite-difference time-domain heating model computations were conducted in the Resnick High Performance Computing Center at the Resnick Sustainability Institute at the California Institute of Technology.

References

1. S. Song, N. Hu and L. Lu, Solid electrolytes for solid-state Li/Na–metal batteries: inorganic, composite and polymeric materials, Chem. Commun., 2022, 58, 12035–12045.
2. T. Mei, H. Zhang and K. Xiao, Bioinspired Artificial Ion Pumps, ACS Nano, 2022, 16, 13323–13338.
3. S. Sun, et al., A focused review on structures and ionic conduction mechanisms in inorganic solid-state proton and hydride anion conductors, Mater. Adv., 2023, 4, 389–407.
4. A. Bhandari and J. Bhattacharya, Origin of Fast Ion Conduction in Li₁₀GeP₂S₁₂, a Superionic Conductor, J. Phys. Chem. C, 2016, 120, 29002–29010.
5. Y.-S. Choi and J.-C. Lee, Electronic and mechanistic origins of the superionic conductivity of sulfide-based solid electrolytes, J. Power Sources, 2019, 415, 189–196.
6. K. Jun, et al., Lithium superionic conductors with corner-sharing frameworks, Nat. Mater., 2022, 21, 924–931.
7. N. Kamaya, et al., A lithium superionic conductor, Nat. Mater., 2011, 10, 682–686.
8. M. A. Kraft, et al., Inducing High Ionic Conductivity in the Lithium Superionic Argyrodites Li_6+xP_1–xGe_xS₅I for All-Solid-State Batteries, J. Am. Chem. Soc., 2018, 140, 16330–16339.
9. A. Manthiram, X. Yu and S. Wang, Lithium battery chemistries enabled by solid-state electrolytes, Nat. Rev. Mater., 2017, 2, 1–16.
10. J. B. Goodenough, Review Lecture – Fast ionic conduction in solid, Proc. R. Soc. London, Ser. A, 1997, 393, 215–234.
11. D. A. Weber, et al., Structural Insights and 3D Diffusion Pathways within the Lithium Superionic Conductor Li₁₀GeP₂S₁₂, Chem. Mater., 2016, 28, 5905–5915.
12. M. A. Kraft, et al., Influence of Lattice Polarizability on the Ionic Conductivity in the Lithium Superionic Argyrodites Li₆PS₅X (X = Cl, Br, I), J. Am. Chem. Soc., 2017, 139, 10909–10918.
13. S. Muy, R. Schlem, Y. Shao-Horn and W. G. Zeier, Phonon–Ion Interactions: Designing Ion Mobility Based on Lattice Dynamics, Adv. Energy Mater., 2021, 11, 2002787.
14. Z. Zhang, et al., Targeting Superionic Conductivity by Turning on Anion Rotation at Room Temperature in Fast Ion Conductors, Matter, 2020, 2, 1667–1684.
15. J. G. Smith and D. J. Siegel, Low-temperature paddlewheel effect in glassy solid electrolytes, Nat. Commun., 2020, 11, 1483.
16. Z. Zhang and L. F. Nazar, Exploiting the paddle-wheel mechanism for the design of fast ion conductors, Nat. Rev. Mater., 2022, 7, 389–405.
17. J. C. Bachman, et al., Inorganic Solid-State Electrolytes for Lithium Batteries: Mechanisms and Properties Governing Ion Conduction, Chem. Rev., 2016, 116, 140–162.
18. Y. Tanaka, et al., New Oxyhalide Solid Electrolytes with High Lithium Ionic Conductivity >10 mS cm⁻¹ for All-Solid-State Batteries, Angew. Chem., 2023, 135, e202217581.
19. B. Singh, et al., Critical Role of Framework Flexibility and Disorder in Driving High Ionic Conductivity in LiNbOCl₄, J. Am. Chem. Soc., 2024, 146, 17158–17169.
20. K. Jun, G. Wei, X. Yang, Y. Chen and G. Ceder, Exploring the soft cradle effect and ionic transport mechanisms in the LiMXCl₄ superionic conductor family, Matter, 2025, 102001,  DOI:10.1016/j.matt.2025.102001.
21. K. Jun, B. Lee, R. L. Kam and G. Ceder, The nonexistence of a paddlewheel effect in superionic conductors, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2316493121.
22. I. Hanghofer, B. Gadermaier and H. M. R. Wilkening, Fast Rotational Dynamics in Argyrodite-Type Li₆PS₅X (X: Cl, Br, I) as Seen by ³¹P Nuclear Magnetic Relaxation—On Cation–Anion Coupled Transport in Thiophosphates, Chem. Mater., 2019, 31, 4591–4597.
23. T. Krauskopf, et al., Comparing the Descriptors for Investigating the Influence of Lattice Dynamics on Ionic Transport Using the Superionic Conductor Na₃PS_4–xSe_x, J. Am. Chem. Soc., 2018, 140, 14464–14473.
24. X. Li and N. A. Benedek, Enhancement of Ionic Transport in Complex Oxides through Soft Lattice Modes and Epitaxial Strain, Chem. Mater., 2015, 27, 2647–2652.
25. S. Muy, et al., Tuning mobility and stability of lithium ion conductors based on lattice dynamics, Energy Environ. Sci., 2018, 11, 850–859.
26. J. Ding, et al., Liquid-like dynamics in a solid-state lithium electrolyte, Nat. Phys., 2025, 21(1), 118–125.
27. G. H. Vineyard, Frequency factors and isotope effects in solid state rate processes, J. Phys. Chem. Solids, 1957, 3, 121–127.
28. Z. Xu, X. Chen, H. Zhu and X. Li, Anharmonic Cation–Anion Coupling Dynamics Assisted Lithium-Ion Diffusion in Sulfide Solid Electrolytes, Adv. Mater., 2022, 34(49), 2207411.
29. M. Brinek, C. Hiebl, K. Hogrefe, I. Hanghofer and H. M. R. Wilkening, Structural Disorder in Li₆PS₅I Speeds ⁷Li Nuclear Spin Recovery and Slows Down ³¹P Relaxation–Implications for Translational and Rotational Jumps as Seen by Nuclear Magnetic Resonance, J. Phys. Chem. C, 2020, 124, 22934–22940.
30. M. Wilkening and P. Heitjans, From Micro to Macro: Access to Long-Range Li⁺ Diffusion Parameters in Solids via Microscopic ^6,7Li Spin-Alignment Echo NMR Spectroscopy, ChemPhysChem, 2012, 13, 53–65.
31. M. K. Gupta, et al., Fast Na diffusion and anharmonic phonon dynamics in superionic Na₃PS₄, Energy Environ. Sci., 2021, 14, 6554–6563.
32. M. K. Gupta, et al., Investigation of Low-Energy Lattice Dynamics and Their Role in Superionic Na Diffusion and Ultralow Thermal Conductivity of Na₃PSe₄ as a Solid-State Electrolyte, Chem. Mater., 2024, 36, 11377–11392.
33. D. Wilmer, H. Feldmann, R. E. Lechner and J. Combet, Sodium Ion Conduction in Plastic Phases: Dynamic Coupling of Cations and Anions in the Picosecond Range, J. Mater. Res., 2005, 20, 1973–1978.
34. T. Morimoto, et al., Microscopic ion migration in solid electrolytes revealed by terahertz time-domain spectroscopy, Nat. Commun., 2019, 10, 2662.
35. A. X. Gray, et al., Ultrafast terahertz field control of electronic and structural interactions in vanadium dioxide, Phys. Rev. B, 2018, 98, 045104.
36. K. Gordiz, S. Muy, W. G. Zeier, Y. Shao-Horn and A. Henry, Enhancement of ion diffusion by targeted phonon excitation, Cell Rep. Phys. Sci., 2021, 2(5), 100431.
37. A. D. Poletayev, et al., The persistence of memory in ionic conduction probed by nonlinear optics, Nature, 2024, 625, 691–696.
38. R. Schlem, A. Banik, S. Ohno, E. Suard and W. G. Zeier, Insights into the Lithium Sub-structure of Superionic Conductors Li₃YCl₆ and Li₃YBr₆, Chem. Mater., 2021, 33, 327–337.
39. B. Zhang, L. Yang, L.-W. Wang and F. Pan, Cooperative transport enabling fast Li-ion diffusion in Thio-LISICON Li₁₀SiP₂S₁₂ solid electrolyte, Nano Energy, 2019, 62, 844–852.
40. S. Stramare, V. Thangadurai and W. Weppner, Lithium Lanthanum Titanates: A Review, Chem. Mater., 2003, 15, 3974–3990.
41. X. Guo, P. S. Maram and A. Navrotsky, A correlation between formation enthalpy and ionic conductivity in perovskite-structured Li_3xLa_0.67−xTiO₃ solid lithium ion conductors, J. Mater. Chem. A, 2017, 5, 12951–12957.
42. M. Först, et al., Nonlinear phononics as an ultrafast route to lattice control, Nat. Phys., 2011, 7, 854–856.
43. T. F. Nova, et al., An effective magnetic field from optically driven phonons, Nat. Phys., 2017, 13, 132–136.
44. W. Hu, et al., Optically enhanced coherent transport in YBa₂Cu₃O_6.5 by ultrafast redistribution of interlayer coupling, Nat. Mater., 2014, 13, 705–711.
45. M. J. Rice and W. L. Roth, Ionic transport in super ionic conductors: a theoretical model, J. Solid State Chem., 1972, 4, 294–310.
46. K. H. Pham, A. K. Lin, N. A. Spear and S. K. Cushing, Laser-driven ultrafast impedance spectroscopy for measuring complex ion hopping processes, Rev. Sci. Instrum., 2024, 95, 073004.
47. K. H. Pham, et al., The Dynamical Role of Optical Phonons and Sublattice Screening in a Solid-State Ion Conductor, J. Am. Chem. Soc., 2025, 147, 26456–26467.
48. T. F. Nova, A. S. Disa, M. Fechner and A. Cavalleri, Metastable ferroelectricity in optically strained SrTiO₃, Science, 2019, 364, 1075–1079.
49. Y. G. Frank, et al., Snapshots of a light-induced metastable hidden phase driven by the collapse of charge order, Sci. Adv, 2022, 8(29), eapb9076.
50. Z. Chen, G. Brocks, S. Tao and P. A. Bobbert, Unified theory for light-induced halide segregation in mixed halide perovskites, Nat. Commun., 2021, 12, 2687.
51. S. Draguta, et al., Rationalizing the light-induced phase separation of mixed halide organic–inorganic perovskites, Nat. Commun., 2017, 8, 200.
52. J. L. Fourquet, H. Duroy and M. P. Crosnier-Lopez, Structural and Microstructural Studies of the Series La_2/3−xLi_3x□_1/3−2xTiO₃, J. Solid State Chem., 1996, 127, 283–294.
53. Y. Inaguma, et al., High ionic conductivity in lithium lanthanum titanate, Solid State Commun., 1993, 86, 689–693.
54. M. L. Sanjuán, M. A. Laguna, A. G. Belous and O. I. V’yunov, On the Local Structure and Lithium Dynamics of La_0.5(Li,Na)_0.5TiO₃ Ionic Conductors. A Raman Study, Chem. Mater., 2005, 17, 5862–5866.
55. D. Nicoletti and A. Cavalleri, Nonlinear light–matter interaction at terahertz frequencies, Adv. Opt. Photonics, 2016, 8, 401.
56. A. Bojtor, et al., Millisecond-Scale Charge-Carrier Recombination Dynamics in the CsPbBr₃ Perovskite, Adv. Energy Sustainable Res., 2024, 5, 2400043.
57. D. S. Venables and C. A. Schmuttenmaer, Far-infrared spectra and associated dynamics in acetonitrile–water mixtures measured with femtosecond THz pulse spectroscopy, J. Chem. Phys., 1998, 108, 4935–4944.
58. G. Mead, H.-W. Lin, I.-B. Magdău, T. F. Miller and G. A. Blake, Sum-Frequency Signals in 2D-Terahertz-Terahertz-Raman Spectroscopy, J. Phys. Chem. B, 2020, 124, 8904–8908.
59. B. Gao, R. Jalem and Y. Tateyama, Ion diffusion driven by dynamic lattice deformations in perovskite solid electrolytes, J. Mater. Chem. A, 2025, 13, 30370–30381.
60. B. I. Greene, P. N. Saeta, D. R. Dykaar, S. Schmitt-Rink and S. L. Chuang, Far-infrared light generation at semiconductor surfaces and its spectroscopic applications, IEEE J. Quantum Electron., 1992, 28, 2302–2312.
61. A. Sood, et al., Direct Visualization of Thermal Conductivity Suppression Due to Enhanced Phonon Scattering Near Individual Grain Boundaries, Nano Lett., 2018, 18, 3466–3472.
62. A. Jonderian, M. Ting and E. McCalla, Metastability in Li–La–Ti–O Perovskite Materials and Its Impact on Ionic Conductivity, Chem. Mater., 2021, 33, 4792–4804.
63. S. Kobayashi, et al., Lithium Lanthanum Titanate Single Crystals: Dependence of Lithium-Ion Conductivity on Crystal Domain Orientation, Nano Lett., 2022, 22, 5516–5522.
64. J. Weidelt, et al., Fundamental Picture of the Conduction Mechanism in Solid-State Polymer Electrolytes Revealed by Terahertz Spectroscopy, J. Phys. Chem. C, 2024, 128, 6868–6876.
65. H. A. Hafez Eid, et al., Terahertz conductivity of polymer electrolytes, in Terahertz Photonics III, ed. M. Jarrahi, D. Turchinovich and S. Preu, SPIE, Strasbourg, France, 2024, 24 DOI:10.1117/12.3016564.
66. Z. Song, L. Xue, Q. Ouyang and C. Song, Impact of a Terahertz electromagnetic field on the ion permeation of potassium and sodium channels, Commun. Chem., 2025, 8, 101.
67. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
68. J. P. Perdew, M. Ernzerhof and K. Burke, Rationale for mixing exact exchange with density functional approximations, J. Chem. Phys., 1996, 105, 9982–9985.
69. G. Henkelman, B. P. Uberuaga and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys., 2000, 113, 9901–9904.
70. X. He, Y. Zhu and Y. Mo, Origin of fast ion diffusion in super-ionic conductors, Nat. Commun., 2017, 8, 15893.
71. D. Sheppard, R. Terrell and G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys., 2008, 128, 134106.
72. E. Gil-González, A. Perejón, P. E. Sánchez-Jiménez, J. M. Criado and L. A. Pérez-Maqueda, Thermoanalytical Characterization Techniques for Multiferroic Materials, in Handbook of Thermal Analysis and Calorimetry, ed. S. Vyazovkin, N. Koga and C. Schick, Elsevier Science B.V., 2018, ch. 16, vol. 6, pp. 643–683.
73. O. Bohnke, H. Duroy, J.-L. Fourquet, S. Ronchetti and D. Mazza, In search of the cubic phase of the Li⁺ ion-conducting perovskite La_2/3−xLi_3xTiO₃: structure and properties of quenched and in situ heated samples, Solid State Ionics, 2002, 149, 217–226.
74. M. Catti, Short-range order and Li⁺ ion diffusion mechanisms in Li₅La₉□₂(TiO₃)₁₆ (LLTO), Solid State Ionics, 2011, 183, 1–6.
75. M. Catti, First-Principles Modeling of Lithium Ordering in the LLTO (Li_xLa_2/3−x/3TiO₃) Superionic Conductor, Chem. Mater., 2007, 19, 3963–3972.
76. M. Catti, Ion Mobility Pathways of the Li⁺ Conductor Li_0.125La_0.625TiO₃ by Ab Initio Simulations, J. Phys. Chem. C, 2008, 112, 11068–11074.
77. A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater., 2015, 108, 1–5.
78. K. Gordiz, kgordiz/modalNEB, 2025.
79. H. R. Seyf, K. Gordiz, F. DeAngelis and A. Henry, Using Green-Kubo modal analysis (GKMA) and interface conductance modal analysis (ICMA) to study phonon transport with molecular dynamics, J. Appl. Phys., 2019, 125, 081101.
80. M. T. Dove, Introduction to Lattice Dynamics, Cambridge University Press, 1993.
81. K. Gordiz and A. Henry, Phonon transport at interfaces: Determining the correct modes of vibration, J. Appl. Phys., 2016, 119, 015101.
82. H. R. Seyf and A. Henry, A method for distinguishing between propagons, diffusions, and locons, J. Appl. Phys., 2016, 120, 025101.
83. K. Gordiz, D. J. Singh and A. Henry, Ensemble averaging vs. time averaging in molecular dynamics simulations of thermal conductivity, J. Appl. Phys., 2015, 117, 045104.
84. K. Gordiz and A. Henry, Phonon Transport at Crystalline Si/Ge Interfaces: The Role of Interfacial Modes of Vibration, Sci. Rep., 2016, 6, 23139.
85. W. Greiner, L. Neise and H. Stöcker, Thermodynamics and Statistical Mechanics, Springer Science & Business Media, 2012.
86. J. Du and A. N. Cormack, The medium range structure of sodium silicate glasses: a molecular dynamics simulation, J. Non-Cryst. Solids, 2004, 349, 66–79.
87. A. N. Cormack, J. Du and T. R. Zeitler, Alkali ion migration mechanisms in silicate glasses probed by molecular dynamics simulations, Phys. Chem. Chem. Phys., 2002, 4, 3193–3197.
88. C. Chen and J. Du, Lithium Ion Diffusion Mechanism in Lithium Lanthanum Titanate Solid-State Electrolytes from Atomistic Simulations, J. Am. Ceram. Soc., 2015, 98, 534–542.
89. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
90. A. J. H. McGaughey and M. Kaviany, Observation and description of phonon interactions in molecular dynamics simulations, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 184305.
91. A. Raj and J. Eapen, Deducing Phonon Scattering from Normal Mode Excitations, Sci. Rep., 2019, 9, 7982.
92. L. L. Taylor, R. E. Scott and J. Qiao, Integrating two-temperature and classical heat accumulation models to predict femtosecond laser processing of silicon, Opt. Mater. Express, 2018, 8, 648.
93. L. L. Taylor, J. Qiao and J. Qiao, Optimization of femtosecond laser processing of silicon via numerical modeling, Opt. Mater. Express, 2016, 6, 2745.
94. C.-W. Wu, X. Ren, W.-X. Zhou, G. Xie and G. Zhang, Thermal stability and thermal conductivity of solid electrolytes, APL Mater., 2022, 10, 040902.
95. J. Neises, et al., Study of thermal material properties for Ta- and Al-substituted Li₇La₃Zr₂O₁₂ (LLZO) solid-state electrolyte in dependency of temperature and grain size, J. Mater. Chem. A, 2022, 10, 12177–12186.
96. S. K. Rhee, Porosity—Thermal conductivity correlations for ceramic materials, Mater. Sci. Eng., 1975, 20, 89–93.
97. A. Zong, et al., Spin-mediated shear oscillators in a van der Waals antiferromagnet, Nature, 2023, 620, 988–993.
98. A. Chouiekh, et al., Experimental and DFT analysis of structural, optical, and electrical properties of Li_3xLa_2/3−xTiO₃ (3x = 0.1, 0.3 and 0.5) solid electrolyte, Ceram. Int., 2023, 49, 25920–25934.
99. Y. Inaguma, L. Chen, M. Itoh and T. Nakamura, Candidate compounds with perovskite structure for high lithium ionic conductivity, Solid State Ionics, 1994, 70–71, 196–202.
100. J. R. Ferraro, Obtaining solid Raman spectra by pelletizing techniques, Spectrochim. Acta, 1964, 20, 901–907.
101. G. Gallot, J. Zhang, R. W. McGowan, T.-I. Jeon and D. Grischkowsky, Measurements of the THz absorption and dispersion of ZnTe and their relevance to the electro-optic detection of THz radiation, Appl. Phys. Lett., 1999, 74, 3450–3452.
102. D. S. Sitnikov, et al., Estimation of THz field strength by an electro-optic sampling technique using arbitrary long gating pulses, Laser Phys. Lett., 2019, 16, 115302.
103. J. B. Baxter and G. W. Guglietta, Terahertz Spectroscopy, Anal. Chem., 2011, 83, 4342–4368.
104. W. Withayachumnankul and M. Naftaly, Fundamentals of Measurement in Terahertz Time-Domain Spectroscopy, J. Infrared Milli Terahz Waves, 2014, 35, 610–637.
105. K. N. Murphy, M. Naftaly, A. Nordon and D. Markl, Polymer Pellet Fabrication for Accurate THz-TDS Measurements, Appl. Sci., 2022, 12, 3475.

---

Footnote

† These authors contributed equally to this work.

---