Understanding Muon Diffusion in Perovskite Oxides below Room Temperature Based on Harmonic Transition State Theory

The positive muon spin rotation and relaxation ($\mu^{+}$SR) method has been widely used to investigate the microscopic properties of solids.amato97 ; blundell04 ; cox09 ; yaounc10 ; higemoto16 ; ito20 This technique involves implanting spin-polarized positive muons ($\mu^{+}$) into lattice interstices as microscopic magnetic probes to observe local magnetic fields arising from surrounding electrons and nuclei. This allows for the study of the electronic state of materials from a microscopic point of view, in a manner similar to the nuclear magnetic resonance technique. In most cases, $\mu^{+}$SR is used as a silent probe; the implanted $\mu^{+}$ is assumed to stay at fixed sites and not to alter the electronic properties of the host material. However, in reality, the implanted $\mu^{+}$ may change local electronic structures and atomic arrangements to some extent due to its positive charge.feyerherm95 ; tashma97 ; gygax03 ; ito09 ; foronda15 In addition, $\mu^{+}$ is relatively mobile in crystalline lattices because of its light mass, which is approximately one ninth that of proton. The diffusion of $\mu^{+}$ can cause drastic changes in $\mu^{+}$SR spectra, such as the narrowing of linesstorchak98 and the averaging of spin precession frequencies,alexandrowicz99 ; ito10 even when host materials do not exhibit any anomalous behavior. Therefore, a deep understanding of both $\mu^{+}$ diffusion and the perturbation induced by $\mu^{+}$ on its surroundings is always necessary for researchers using $\mu^{+}$SR spectroscopy, regardless of whether they are interested in such phenomena or not.

The diffusion of $\mu^{+}$ has been extensively studied in elemental metals and its quantum nature at cryogenic temperatures has been unambiguously established.storchak98 For instance, anomalous temperature dependencies of an interstitial $\mu^{+}$ hopping rate in copper were successfully described by a microscopic quantum tunneling model that takes into account interactions between $\mu^{+}$ and the lattice, and conduction electrons.kadono89 ; luke91 ; kondo84-1 ; kondo84-2 ; yamada84 In contrast, $\mu^{+}$ implanted into solid oxides is usually regarded as immobile below room temperature because it is thought to be tightly bound to an oxygen atom in the host lattice to form a muonic analogue of the hydroxy group. Many researchers have accepted this assumption and applied $\mu^{+}$SR spectroscopy to investigations of magnetism and ion diffusion in oxides, where anomalies appearing in $\mu^{+}$SR spectra have simply been interpreted as manifestations of the material’s nature.hord10 ; koda19 ; sugiyama13 On the other hand, $local$ $\mu^{+}$ hopping has undoubtedly been observed below room temperature in transition-metal oxide antiferromagnets with the corundum structure.dehn20 ; dehn21 Moreover, $global$ $\mu^{+}$ diffusion with an activation energy of $\sim$0.1 eV has been reported in some chemically-substituted perovskite oxides, such as Sc-doped SrZrO₃hempelmann98 and BaTiO_3-xH_x oxyhydrideito17 , at cryogenic temperatures. While the quantum effect associated with the light mass is inferred to be important in $\mu^{+}$ diffusion, the reason for the small activation energy despite the formation of the strong O$\mu$ bond has not yet been quantitatively understood.

In this article, we report $\mu^{+}$SR and density functional theory (DFT) calculations studies on $\mu^{+}$ diffusion in insulating KTaO₃ with a cubic perovskite structure to gain a comprehensive understanding of the $\mu^{+}$ diffusion in the perovskite oxide lattice. KTaO₃ is an ideal system for investigating the $\mu^{+}$ diffusion since it has dense nuclear spins and preserves the cubic perovskite structure even at cryogenic temperatures. These features facilitate experimental detection of the $\mu^{+}$ diffusion and its modeling in DFT calculations. DFT has recently been used extensively to provide valuable information for identifying $\mu^{+}$ stopping sites and for estimating their stabilities and charge states as hydrogen-like defects in crystalline lattices from a theoretical viewpoint.foronda15 ; koda19 ; moller13 ; moller13-2 ; vieira16 ; dehn20 ; dehn21 ; ito22 However, the application of DFT to support $\mu^{+}$SR data analysis has mostly been limited to a classical level, despite the importance of quantum effects arising from the light mass of $\mu^{+}$. A recent DFT study by Onuorah et al. highlights the importance of zero-point energies (ZPEs) in finding stable $\mu^{+}$ sites in metals,onuorah19 but their impact on $\mu^{+}$ diffusion has not been discussed in detail. In this study, we handle ZPEs for $\mu^{+}$ in KTaO₃ within the frameworks of double Born-Oppenheimer and harmonic approximations,onuorah19 and show how the effective activation barrier for interstitial $\mu^{+}$ diffusion is lowered due to the difference in ZPEs for $\mu^{+}$ at positions in the O$\mu$ bonding equilibrium state and a bond-breaking transition state based on the harmonic transition state theory.sholl09 ; sholl14

This paper is organized as follows. After a brief description of the experiments in Sec. II, we report our experimental results in Sec. III together with the results of data analysis using the dynamical Gaussian Kubo-Toyabe theory and a two-state model. In Sec. IV and V, we describe computational details and results of DFT calculations performed for providing a semi-quantitative description of long-range interstitial $\mu^{+}$ diffusion in KTaO₃ below room temperature. In Sec. VI, we discuss the implications of these results for interstitial $\mu^{+}$ diffusion in oxides with a specific focus on the impact of ZPEs on the effective barrier height for $\mu^{+}$ hopping. Finally, we summarize our conclusions in Sec. VII.

High-purity (4N) KTaO₃ single crystals of 10$\times$10$\times$0.5 mm³ were obtained from MTI Corporation, USA. The single crystals were grown by the top-seed flux method and cut along the cubic (001) plane. Time-differential $\mu^{+}$SR experiments were performed at Japan Proton Accelerator Research Complex (J-PARC), Japan, and Paul Scherrer Institut (PSI), Switzerland, using a spin-polarized surface muon beam in a longitudinal field configuration with its initial polarization direction nominally parallel to the beam axis. The general purpose $\mu^{+}$SR spectrometers at D1 and S1 areas of J-PARC with a conventional ⁴He flow cryostat and the Low Temperature Facility $\mu^{+}$SR instrument with a ³He-⁴He dilution refrigerator at PSI were used for the temperature ranges of 3-300 K and 0.1-5 K, respectively. The KTaO₃ single crystals were mounted on a silver sample holder with the (001) plane perpendicular to the beam axis. The asymmetry of $\mu^{+}$ beta decay, $A$, was monitored by the “Forward” and “Backward” positron counters as a function of the elapsed time $t$ from the instant of muon implantation.

$A(t)$ can be decomposed into partial asymmetries from muons stopped in the sample, $A_{s}(t)$, and the silver sample holder, $A_{BG}$, the latter of which is substantially a time-independent constant. The magnitude of $A_{BG}$ was determined using a reference sample (holmium) with similar geometry to the KTaO₃ sample and it was subtracted from $A(t)$. The remaining $A_{s}(t)$ was divided by $A_{s}(t=0)$ to obtain the time-evolution of spin polarization $P(t)$ for muons that stopped in the sample. In this normalized form, $\mu^{+}$SR data sets recorded with different spectrometers can be directly compared.

The muons implanted into KTaO₃ depolarize mainly due to dipole-dipole interactions with surrounding ¹⁸¹Ta, ³⁹K, and ⁴¹K nuclei having natural abundances of 99.988%, 93.258%, and 6.730%, respectively.berglund11 $P(t)$ in such a dense nuclear spin system can be approximately described by averaging spin precession signals for muons that feel random nuclear dipolar fields with an isotropic Gaussian distribution (local field approximation).hayano79 When muons are immobile in the $\mu^{+}$SR time window of up to $\sim 20~{}\mu$s, $P(t)$ in a zero applied field (ZF) is expected to have a Gaussian-like shape, which reflects the Gaussian distribution of the local fields. As the hopping rate of muons increases, the shape of $P(t)$ in ZF should gradually change from the Gaussian-like curve to an exponential-like curve with a decrease in a muon spin relaxation rate (motional narrowing). This behavior can be approximately described with the dynamical Gaussian Kubo-Toyabe (DKT) theory based on a strong collision approximation.hayano79 In a faster diffusion regime, the effects of muon trapping at dilute defects and/or detrapping from them may appear in $P(t)$, which can be approximately expressed using a two-state (2S) model.hempelmann98 ; ito17 ; lord01 ; borghini78

$\mu^{+}$SR spectra in ZF exhibit a complicated temperature dependence below room temperature, with the spectral shape changing between Gaussian-like and exponential-like functions (Fig.1(a)). To obtain an overview of the temperature variation, the spectra were tentatively fitted to a stretched exponential function, $p_{0}\exp(-(\Lambda t)^{\beta})$. The initial polarization $p_{0}$ is approximately unity over the entire temperature range (Fig.1(b)), indicating that almost all the muons implanted into the sample are in diamagnetic environments. The exponent $\beta$, which was allowed to vary between 1 and 2, gradually decreases to 1 as the temperature increases from 1 to 10 K (Fig. 1(c)). This is accompanied by a decrease in the muon spin relaxation rate $\Lambda$ (Fig. 1(d)). This behavior strongly suggests that the hopping motion of $\mu^{+}$, which is detectable within the $\mu^{+}$SR time window, is activated even below 10 K. Further information on the $\mu^{+}$ state in this temperature regime can be obtained through curve fitting analysis using a DKT function,hayano79 which is described in Sec. III.1.

The exponential-like $P(t)$ at around 10-75 K changes to a Gaussian-like curve toward 135 K with an increase in $\Lambda$. Similar behavior has been reported in Sc-doped SrZrO₃,hempelmann98 BaTiO_3-xH_x oxyhydride,ito17 ReO₃,lord01 and impurity-doped Nb,borghini78 ascribed to an increasing rate of $\mu^{+}$ trapping at a dilute defect with increasing temperature. Therefore, this behavior in KTaO₃ can also be regarded as a signature of $\mu^{+}$ trapping at a dilute defect (defect A) that occurs via long-range interstitial $\mu^{+}$ diffusion. The exponential-like shape is recovered at around 200 K with a local minimum in $\Lambda$, and then the Gaussian-like shape is restored with increasing temperature toward 300 K. The former and latter behavior can be attributed to an increasing rate of $\mu^{+}$ jumps between defect A through the interstitial diffusion path, and that of $\mu^{+}$ trapping at another more dilute defect (defect B), respectively. All these processes above 75 K can be approximately expressed using the 2S model, which is described in Sec. III.2.

All calculations were performed within the DFT framework using the Quantum ESPRESSO package.giannozzi09 ; giannozzi17 The generalized gradient approximation using Perdew-Burke-Ernzerhof (PBE) exchange correlation functional was adopted. Projector augmented-wave potentials with H($1s$), K(3$s$, 3$p$, 4$s$), Ta(5$s$, 5$p$, 5$d$, 6$s$) and O(2$s$, 2$p$) valence states were used. Electron wave functions were expanded in plane waves with a cutoff of 70 Ry for kinetic energy and 600 Ry for charge density. The KTaO₃ structure was described using a $3\times 3\times 3$ supercell of the 5-atom cubic primitive cell. A H atom was placed at an interstitial position in the host supercell to model a dilute muon defect in the KTaO₃ lattice. While a 2$\times$2$\times$2 supercell was used for similar calculations in Ref. sholl14, ; gomez07, , we more carefully employed the 3$\times$3$\times$3 supercell to exclude the possibility of artificially stabilizing antiferro-distortions, which can occur in the 2$\times$2$\times$2 supercell. Considering the amphoteric nature of the H defect, calculations were performed for defective supercells with $q=+1$, 0, and -1, where $q$ is the total charge of electrons and ions in the defective supercell. For the $q=+1$ or -1 cases, a uniform jellium background with an opposite charge was used to maintain charge neutrality.

Structure optimization. Brillouin zone sampling with a Monkhorst-Pack k-point mesh of $2\times 2\times 2$ and Gaussian smearing with a broadening of 0.01 Ry were applied for structure optimization calculations. Atomic positions were relaxed while the lattice constant for the supercell was kept fixed at the calculated value of 12.024Å  for stoichiometric cubic KTaO₃, which was determined via full structure optimization. This value is consistent with the experimental lattice constant of 3.988Å  for the 5-atom primitive cell.samara73 In all cases, atomic forces were converged to within $7.7\times 10^{-3}$ eV/Å.

Defect formation energy. The self-consistent field calculation for the optimized structure was refined using the optimized tetrahedron methodkawamura14 with a finer k-point mesh of $5\times 5\times 5$ to avoid artifacts associated with the Gaussian smearing. The classical formation energy of the interstitial H defect with a charge $q$, $E_{f}({\rm H}^{q})$, was calculated in accordance with a standard procedure described in Ref. iwazaki10, ; iwazaki14, :

where $E_{t}({\rm H}^{q})$ is the total energy of the defective supercell with a charge $q$, $E_{t}({\rm host}^{0})$ is the total energy of a neutral host supercell, $\mu_{{\rm H}_{2}}$ is the chemical potential of an ${\rm H}_{2}$ molecule, $\mu_{\rm VBM}$ is the chemical potential of the electron reservoir at the valence-band maximum, and $E_{\rm F}$ is the Fermi energy with reference to $\mu_{\rm VBM}$. $E_{corr}({\rm H}^{q})$ is a small potential-alignment term for charged defective supercells,persson05 which was determined from the difference in electrostatic potentials at an interstitial position about 10Å away from the H defect.

Classical minimum energy paths. For the most stable H configuration determined by the above calculations, classical minimum energy paths (MEPs) between equivalent nearest neighbor H sites were searched using the nudged elastic band methodsholl09 with seven images. Brillouin zone sampling with a Monkhorst-Pack k-point mesh of $2\times 2\times 2$ and Gaussian smearing with a broadening of 0.01 Ry were applied. Atomic positions were relaxed with keeping the lattice constant fixed at the optimized value for stoichiometric KTaO₃ through iterations for finding the MEPs. A climbing image algorithmhenkelman00 was used to identify saddle points in the MEPs, which represent transition states (TSs) of corresponding H transfer events, and precisely compute classical barrier heights $E_{0}$ at the saddle points. In all cases, the norm of the force orthogonal to the path was converged to within 0.1 eV/Å.

ZPE-corrected $\mu^{+}$ jump rates. The interstitial diffusion of H (pseudo) isotopes in the classical over-barrier hopping regime can be approximately explained using the harmonic transition state theory, which takes into account isotopic effects arising from zero-point vibrations within the harmonic approximation.sholl14 We adopted the double Born-Oppenheimer approximation for computing muon vibrational frequencies, which requires, in addition to the usual Born-Oppenheimer approximation being satisfied, that the vibrational motion of muons can be decoupled from that of other nuclei.onuorah19 This is justified because the muon’s mass is much lighter than the mass of any other nucleus. In the double Born-Oppenheimer and harmonic approximations, ZPEs for muons in the initial (and final) equilibrium state (ES) and TS are expressed as $\sum_{i=1}^{3}\frac{1}{2}\hbar\omega_{i}^{\rm ES}$ and $\sum_{i=1}^{2}\frac{1}{2}\hbar\omega_{i}^{\rm TS}$, respectively, where $\omega_{i}^{\rm ES}$ and $\omega_{i}^{\rm TS}$ are harmonic vibrational angular frequencies for muons in the ES and TS. It should be noted that the sum for the TS runs over only two modes because the other mode is associated with a displacement along the MEP and therefore has an imaginary frequency at the saddle point. Consequently, a ZPE-corrected barrier height for the MEP with the classical barrier height $E_{0}$ is given as,

According to the harmonic transition state theory, a ZPE-corrected muon jump rate $k_{\rm ZPC}$ can be calculated via the following equation,

where $f(x)=\sinh(x)/x$.sholl14 At high temperatures, $f(\hbar\omega/2k_{B}T)\rightarrow 1$, so that eq. 3 reduces to the simple Arrhenius law with the activation energy of $E_{0}$. At low temperatures,

where the preexponential factor $k_{B}T/h$ is in the order of 10¹² Hz in the temperature range 75-300 K, which is close to the Debye frequency of KTaO₃. To calculate $k_{\rm ZPC}$ via eq. 3, $\omega_{i}$s for the ES and TSs in all possible MEPs were computed using a density functional perturbation theory module in the Quantum ESPRESSO package within the framework of the double Born-Oppenheimer and harmonic approximations. In the density functional perturbation theory calculations, the muon was treated as a H isotope having the mass of 0.113$m_{p}$, where $m_{p}$ is the proton’s mass.

KTaO₃ exhibits quantum paraelectric behavior at low temperatures, where the cubic structure and a virtual ferroelectric structure are nearly degenerate and the cubic symmetry is maintained via zero-point quantum fluctuations over these configurations. The atomic displacement in the virtual ferroelectric phase from corresponding atomic positions in the cubic phase was estimated from DFT calculations to be at most 0.026Å.esswein22 This is one order of magnitude smaller than the local displacement of muon-neighboring atoms from our calculations, which reflects the chemical bonding state with muons. Therefore, the dynamical ferroelectric displacement associated with the quantum paraelectricity is likely to have only minor effects on muon kinetics, which involves the formation and breaking of the chemical bonds.

Since the potentials experienced by oxygen-bound protons in perovskite oxides generally reveal significant anharmonic terms along some directions,iwazaki10 it is not obvious whether the harmonic method can serve as a good approximation for describing ZPEs of muons in KTaO₃. In particular, the adiabatic potentials corresponding to the O$\mu$ stretching vibrations (the $\omega_{1}$ modes for ES and TS(R) in Table 3) are highly asymmetric and anharmonic, and the deviation from the quadratic curve for the harmonic approximation should be significant except at the bottom of the potential wells.

A similar situation has been reported for vibrational modes of muons at around equilibrium positions in some fluorides.moller13-2 In that paper, a comparison was made between harmonic ZPEs and those obtained by solving the Schrödinger equation for a muon in one-dimensional adiabatic potentials calculated along each vibration direction. Despite significant anharmonicity identified in some vibrational modes, there was surprisingly good agreement between the two sets of ZPEs.

We adopted the same strategy for KTaO₃ and solved the Schrödinger equation for a muon in adiabatic potentials corresponding to the seven muon vibrational modes in Table 3 with real frequencies using a finite-differences method. As a result, we obtained ZPEs of 0.9632, 0.5488 and 0.8621 eV for ES, TS(T), and TS(R), respectively, and the bottleneck barrier height of 0.1173 eV. The good agreement with corresponding harmonic values in Table 3 strongly suggests the validity of the harmonic approximation for describing the vibrational states of muons in perovskite oxides.

The ZPE correction based on the harmonic approximation successfully accounted, at least semi-quantitatively, for the long-range interstitial $\mu^{+}$ diffusion in KTaO₃ with the activation energy of $\sim 0.1$ eV. This relatively low-cost approach is probably also applicable to explaining that in Sc-doped SrZrO₃ and BaTiO_3-xH_x, which is characterized by similar activation energies.hempelmann98 ; ito17 Furthermore, we can gain general insight into $\mu^{+}$ states in oxides from the calculation results for KTaO₃ as described below.

The long-range interstitial $\mu^{+}$ diffusion in oxides can usually be expressed as the combination of the transfer and reorientation-type hopping processes as in the case of KTaO₃. The transfer-type process is characterized by the temporary breaking of the strong O$\mu$ bond in the TS(T), which results in a marked difference between ZPEs for the ES and TS(T). This can not only cause a significant reduction in the effective barrier height, but also change the ES for muons in extreme cases; from a viewpoint of ZPEs, muons prefer to stay in large interstitial or vacancy sites rather than forming the strong O$\mu$ bond at the conventional proton site adjacent to oxygen. This means that muons and hydrogen may be found in different stable sites due to the ZPE effect.

The barrier height for the transfer-type path between the conventional proton sites can vary widely depending on host crystal structures. When the proton sites are dense in the host lattice, as in the KTaO₃ perovskite, the barrier height can be relatively low in comparison with dilute cases; when the proton sites are far apart from each other, a higher energy for completely breaking the O$\mu$ bond ($\sim$5 eV, estimated from the average OH bond energy in H₂O) is expected to be required for transferring a muon from one site to another. Therefore, whether global muon diffusion becomes important in muon spin relaxation below room temperature should strongly depend on the crystal structure of host oxides.

On the other hand, the energy required for transferring a muon through the reorientation-type path is considered to be relatively independent of the crystal structure of host oxides because this process does not involve the temporary breaking of the O$\mu$ bond. Thus, provided that the crystal symmetry allows for such reorientation, the local motion of muons around oxygen may be characterized by a small barrier height of 0.1$\sim$0.2 eV. Since $\mu^{+}$SR is a local probe, even such local motion can cause significant changes in $\mu^{+}$SR spectra.dehn20 ; dehn21 Therefore, researchers who use $\mu^{+}$SR to study oxides need always keep in mind the possibility that muons can move even below room temperature when examining $\mu^{+}$SR data, and then carefully discuss physical properties of the oxides, such as magnetism and ionic conductivity.

It is worth mentioning that the above discussion on global and local motion of muons in oxides may also be applied to materials that involve elements with high electronegativity, such as C, N, S, Se, and halogens, which tend to create strong chemical bonds with implanted muons.brewer86 ; nishiyama03 ; lancaster07 It should also be noted that quantum tunneling of muons is beyond the scope of the DFT and harmonic transition state theory framework and higher-level calculations are necessary for correctly handling this phenomenon, such as using the ab initio path-integral molecular dynamics method.kimizuka18

The harmonic ZPEs for muons can be converted to those for protons under the double Born-Oppenheimer approximation by multiplying by the square root of the muon-proton mass ratio ($\sqrt{m_{\mu}/m_{p}}\sim$ 1/3). From these scaled values, the bottleneck process for long-range diffusion of protons is inferred to remain the transfer type, with an activation barrier of about 0.3 eV. This value is considerably larger than that for muons, suggesting that the over-barrier-type diffusion of hydrogen (pseudo) isotopes in KTaO₃ exhibits a normal kinetic isotope effect. However, the justification for applying the double Born-Oppenheimer approximation to protons remains to be established, and more careful treatment may be necessary.

In summary, we carried out $\mu^{+}$SR and DFT studies on $\mu^{+}$ diffusion in the cubic perovskite KTaO₃ to gain a comprehensive understanding of the $\mu^{+}$ diffusion in the perovskite oxide lattice. In the $\mu^{+}$SR experiment, we observed motional narrowing behavior, the evidence of interstitial $\mu^{+}$ diffusion, at cryogenic temperatures. The $\mu^{+}$ diffusion observed below room temperature was classified into two types; one is possibly the quantum-tunneling type, which is dominant below 75 K, and the other is the classical over-barrier type, which shows a major contribution at higher temperatures. The latter was indirectly identified through interactions between muons and defects, characterized by a small activation energy of $\sim 0.1$ eV. We successfully accounted, at least semi-quantitatively, for the small activation energy associated with the long-range interstitial $\mu^{+}$ diffusion using the DFT and harmonic transition state theory framework. The significant reduction in the effective barrier height for the transfer-type $\mu^{+}$ hopping path due to the ZPE correction is the key to understanding the small activation energy for KTaO₃ and other perovskite oxides.

We also obtained important insight into $\mu^{+}$ states in other oxides from the calculation results for KTaO₃. The effective barrier height for the transfer-type path can vary widely depending on the crystal structure of host oxides, whereas that for the reorientation-type path is expected to be relatively structure-independent and in the range of approximately 0.1 to 0.2 eV if the crystal symmetry allows for such reorientation. The low barrier height suggests that significant changes in $\mu^{+}$SR spectra can arise even below room temperature due to the local $\mu^{+}$ motion regardless of whether the global motion is activated or not. This may prompt a revision of previous interpretations of $\mu^{+}$SR anomalies in some oxides, which have been claimed to be directly associated with the intrinsic properties of the materials.