Distinguishing Ion Dynamics from Muon diffusion in Muon Spin Relaxation

We propose a model to describe the fluctuations in the internal magnetic field due to ion dynamics observed in the muon spin relaxation ($\mu$SR) by an Edwards-Anderson type autocorrelation function that separates the quasi-static and dynamic components of the correlation by a parameter $Q$ (where $0\le Q\le1$). Our Monte Carlo simulations for this model showed that the time evolution of muon spin polarization deviates significantly from the Kubo-Toyabe (KT) function. To further validate the model, the results of simulations were compared with the $\mu$SR spectra observed in a hybrid organic-inorganic perovskite FAPbI$_3$ [with FA referring to HC(NH$_2)_2$], where local field fluctuations associated with the rotational motion of FA molecules and quasi-static fields from the PbI$_3$ lattice are presumed to coexist. The least-squares curve fitting showed reasonable agreement with the model with $Q=0.947(3)$, and the fluctuation frequency of the dynamical component was obtained. This result opens the door to the possibility of experimentally distinguishing fluctuations due to dynamics of ions around muons from those due to self-diffusion of muons. Meanwhile, it suggests the need to carefully consider the spin relaxation function when applying $\mu$SR to the issue of ion dynamics.


I. INTRODUCTION
Muon spin rotation (µSR) is an experimental method to probe magnetic fields in matter, where the spin-polarized muons (µ + ) stopped in the sample can directly provide information on the magnitude of the local field at the interstitial site(s) via the frequency of their Larmor precession [1].The muon gyromagnetic ratio γ µ (= 2π × 135.539MHz/T) is 3.18 times higher than that of protons (= 2π × 42.577 MHz/T, which is the highest among all the nuclei of stable elements); and thus, it is the most sensitive probe of the internal magnetic field upon their implantation to materials.Muons are provided as a 100%-spin-polarized beam by proton accelerator facilities, which enables µSR measurements in zero magnetic field.In addition, µSR has various advantages: the muon implantation energy is high enough (≥ 4 MeV, corresponding to the stopping range ≥ 0.1 g/cm 2 ) to be surface-independent (bulksensitive), and the implanted muons (volume concentration ≃ 10 5 cm −3 for the high-flux beams at J-PARC MLF) decay to positrons and neutrinos with an average lifetime of 2.198 µs; therefore, they do not accumulate in the sample unlike other ion beam irradiation methods.Most notably, the spatial distribution of the emitted high-energy decay positrons has a large asymmetry (1/3 when averaged over the positron energy) with respect to the muon spin polarization, and the time evolution of spin polarization can be observed by measuring the timedependent asymmetry (which is called the "µSR time spectrum").
Recently, attempts have been made to observe the diffusive motion of ions in battery materials and ionic conductors by µSR measurements under zero or longitudinal magnetic field (ZF/LF-µSR) [2][3][4][5][6][7][8][9][10][11], where the time-dependent fluctuations of a weak random local field H(t) (∼10 −4 T) due to the magnetic dipole moments of cation nuclei are monitored via the spin relaxation of muons implanted into the materials of interest.In such studies, it is often an issue whether the observed fluctuations are caused by cation diffusion or muon self-diffusion, but it was believed that the distinction could not be made solely from µSR spectra.Therefore, the attribution of the origin of the fluctuations has been based on information from other experimental methods and, more recently, on inferences from ab initio density functional theory (DFT) calculations.
The muon spin relaxation induced by H(t) under ZF/LF is considered well approximated by the Kubo-Toyabe (KT) function [12], and it has been routinely used in the analysis of µSR time spectra by curve fitting based on the leastsquares method.In deriving the KT function, the effect of the fluctuations is incorporated using a strong-collision model (an approximation of the Markovian process), where H(t) is assumed to have no correlation before and after a change.As illustrated in Fig. 1(a), this approximation is considered to hold well in the case of the jump diffusion of muons, because the relative configuration of all nuclear magnetic moments is swapped simultaneously in a single jump.On the other hand, in the case of ion diffusion around muons, it is naturally expected to be a rare event that all ions move simultaneously (especially for slow diffusion), resulting in a strong autocorrelation in H(t) before and after its change [see Fig. 1(b)].This is even more so when multiple kinds of nuclei contribute to the internal magnetic field (e.g., 7 Li and 59 Co in Li x CoO 2 ) while only some of them exhibit diffusion.Nevertheless, the KT function has traditionally been used in the study of ion diffusion.
A similar situation has been observed in very recent µSR studies of hybrid organic-inorganic perovskites (HOIPs) [13,14], where the thermally activated local rotational motion of organic cation molecules contributes to the fluctuation of H(t) while the internal fields from the PbI 3 lattice remain quasi-static [see Fig. 1(c)].Moreover, it has been suggested that the observed µSR time spectra may not always be satisfactorily reproduced by the KT function upon the onset of molecular motion.These circumstances led us to examine the corresponding longitudinal spin relaxation function for such cases.
In this paper, we propose a model in which the fluctuations of H(t) are described by an Edwards-Anderson-type autocorrelation function with the parameter Q used to divide the quasi-static and dynamical components in the correlation.Our Monte Carlo simulations on the time evolution of muon spin polarization for this model indicate that the relaxation lineshape could be significantly different from that of the KT function.This difference will open the door to the possibility of experimentally distinguishing between muon self-diffusion and cation dynamics based solely on µSR measurements.We show that the origin of H(t) fluctuations can indeed be inferred by comparing relaxation functions based on this model with high-precision µSR time spectra observed in formamidinium lead iodide [FAPbI 3 , where FA denotes HC(NH 2 ) 2 ] [14].Furthermore, on the basis of this result, we will discuss the problems encountered in the preceding studies of ion diffusion by µSR.

II. MODEL A. Relaxation function for the static internal field
Since our goal is to investigate how different autocorrelations of H(t) can change in the spin relaxation function, we will introduce the KT function in a simple classical way in the following.When muons are initially polarized along the z-direction and subjected to a homogeneous magnetic field H = (H x , H y , H z ), the time evolution of muon polarization is described by where H = |H|.Such a single-frequency rotation is observed when muons are subjected to a uniform external magnetic field (not parallel to z).In the case of random local fields exerted from nuclear magnetic moments, H varies with the muon position r µ so that the µSR spectrum corresponds to a random sampling of H(r µ ) = (H x (r µ ), H y (r µ ), H z (r µ )), which is described by a density distribution: The time evolution of spin polarization projected to the z-axis is then given by ) where α = x, y, and z.Unless n(H) is a delta function, g z (t) accompanies relaxation due to the loss of phase coherence in the spin precession.
Specifically, provided that the static distribution of H(r µ ) is approximated by an isotropic Gaussian distribution, the corresponding relaxation function is obtained analytically by executing the integral in Eq. ( 4) to yield the static KT function [12], where the term 1/3 originates from the spatial average of cos 2 θ in Eq. ( 2), corresponding to the probability that H(r µ ) is parallel to z.The linewidth ∆ is determined by the second moment of the dipole field from the nuclear magnetic moments of the neighboring atoms (including noncationic ions), where is the dipole tensor for the ith nuclear magnetic moment γ I I i situated at r i from the muon site.Since the magnetic dipolar field decays proportionally to 1/r 3 i , the magnitude of ∆ is approximately determined by the nearest neighboring (nn) nuclear magnetic moments around a muon.The Gaussian distribution for n(H) works reasonably well when the number N of nn nuclear magnetic moments satisfies the condition N ≥ 4. On the other hand, it is a relatively poor approximation for N ≤ 3, and it is necessary to calculate g z (t) by treating the muon-nucleus cluster as a few-quantum spin system using the density matrix method [15,16].With recent advances in the computational environment, such calculations have been performed for nuclear magnetic moments located farther from a muon, and good agreement of calculation results with experimental results has been obtained [17].

B. Effect of fluctuating internal fields
In general, the effect of fluctuations in linear response theory is treated by a Gaussian-Markov process.In this model, the autocorrelation of fluctuations is given by where ... denotes the thermal average over the canonical ensemble.(For a more rigorous treatment, see, for example, Ref. 18.) The strong-collision model (also called the random Markov process or random phase approximation) is a simplified version of this model, where it is assumed that the density distribution n α (H) is constant at any time t and that H(t) changes with an average rate ν = 1/τ without correlation before and after a change in H.The resulting relaxation function is a good approximation of a Gaussian process except in the νt ≪ 1 region [19].In this case, the relaxation function G z (t), which incorporates the effect of fluctuations on the static relaxation function g z (t), is derived by solving the following integral equation (a linear Volterra equation of the second kind): This integral equation can be solved analytically by the Laplace transform of g z (t) and G z (t) [12] or direct numerical integration [20].The dynamical KT function G KT z (t; ∆, ν) obtained in this way successfully reproduced the spin relaxation due to muon jump diffusion in a nonmagnetic metal [21][22][23].This indicates that fluctuations in the internal magnetic field due to muon self-diffusion can be well described by the strong-collision model with Eq. ( 9).
On the other hand, fluctuations in H due to ion diffusion do not necessarily occur in its entire amplitude; let us consider the situation in Fig. 1(b) where part of H is due to contributions from cations or anions that remain stationary while other ions exhibit local motion.A similar situation is expected in the case of FAPbI 3 shown in Fig. 1(c), where cation molecules are in random rotational motion relative to the static PbI 3 lattice.In these cases, it is clear that the autocorrelation function differs from Eq. ( 9), and it may be given by where Q is the parameter that describes the fractional amplitude exhibiting fluctuations with the rate ν i , and ν µ is the fluctuation rate due to muon diffusion.The equation with ν µ = 0 corresponds to a model introduced to describe the slowing down of magnetic fluctuations (with the rate ν i ) in spin glasses, where 1 − Q is interpreted as an order parameter that describes "ordering in time" [24,25].This model was also applied to deriving relaxation functions for analyzing the µSR spectra observed in dilute-alloy spin glasses [26].However, the corresponding G z (t) has only been discussed in approximate forms obtained as analytical expressions when ν µ ≪ ∆ ≪ ν i [13, 26], and an extension to the general case is still awaited.Therefore, in the following, we investigate the behavior of G z (t) corresponding to arbitrary Q and ν i by Monte Carlo simulations.

C. Monte Carlo simulations
We performed numerical simulations of the spin relaxation function for Q > 0 and ν µ = 0, which we call G ID z (t; ∆, ν i , Q), where the ion dynamics is the sole origin of fluctuations with the amplitude Q∆ 2 and the frequency ν i : Figure 2 schematically shows the situation in Eq. ( 12) under these conditions.In our Monte Carlo simulation, the local field that each muon experiences was approximately expressed as a vector sum of static and dynamical fields.These fields were randomly generated according to the Gaussian probability distribution n α (H) with the widths of (1 − Q)∆/γ µ and √ Q∆/γ µ .Note that the assignment of 1 − Q and Q is opposite to that in Refs.24 and 25, where Q represents the static component; this is so defined such that the case Q = 0 should be reduced to the conventional KT function (with ν µ > 0 for muon diffusion).The fluctuation of the dynamical component was implemented in accordance with the strong-collision model by regenerating the dynamical field at the average rate of ν i .The time evolution of a muon spin, initially oriented to the z-direction, in the local field was calculated and its z-component was averaged for 10 8 muons to obtain G ID z (t) for ZF.For finite LFs, the time evolution was calculated in the field obtained by the vector sum of the local field and the LF along the z-direction.
As shown in Figs.3(a)-3(c), in the case of 0 < Q < 1, the 1 3 term once decays as ν i increases, then begins to recover with the apparent linewidth decreasing to √ 1 − Q∆.Provided that the µSR spectrum showing such temperature dependence is analyzed with the usual KT function [Fig.3(d)], the apparent fluctuation rate would increase and then decrease.Thus, it is likely to lead to an incorrect interpretation that the cause of the fluctuations is different from a simple thermal excitation process.An important clue to avoiding this is to determine whether the linewidth changes at both ends of the apparent increase and decrease in ν derived from the KT function.In Figs.3(b) and 3(c), it can be seen that as ν i increases, the Q∆ 2 contribution disappears owing to motional narrowing, and the relaxation rate decreases (converging to √ 1 − Q∆) while maintaining the Gaussian linshape.Such behavior is actually observed in HOIP [13,14], and the change in the overall lineshape is found to correspond approximately to the case of Q ≈ 0.95 (see below).Note here that the linewidth behaves nonlinearly with respect to Q.As shown in Fig. 2(b), it changes steeply around Q → 0 for The LF dependence of G ID z (t; ∆, ν i , Q) is indispensable in- formation for a more accurate determination of fluctuation frequency in µSR measurements.Let us compare the field dependence of the time spectrum for Q = 0.5 and ∆/ν i = 0.5 shown in Fig. 3(e) with that of the conventional KT function in detail.First, the KT function calculated for ∆/ν i = 1 2 √ 0.5 is shown in Fig. 3(f); this is because the linewidth for the fluctuations in the former case is √ Q∆ rather than ∆.The oscillations around ∆ •t ≃ 2 [originating from the static component corresponding to (1 − Q)] seen in Fig. 3(e) are absent, indicating a qualitatively different behavior.Furthermore, assuming practical situations of curve fitting with the KT function, the results with freely varying ν i to reproduce the behavior in Fig. 3(e) are shown in Figs.3(g) and 3(h).Figure 3(g) reproduces the low-field spectrum well, whereas the fit is poor on the high-field side; Fig. 3(h) shows the opposite trend.From these comparisons, we can conclude that the curve fitting with the KT function does not accurately reproduce the field dependence of the spectrum when Q < 1.Now, let us discuss the behavior of the relaxation function for some extreme cases, for which we define ∆ µ = (1 − Q)∆ and ∆ i = √ Q∆.When the ion dynamics is fast enough to satisfy ν i ≫ ∆, the spin relaxation for that part is Then, if the muon is quasi-static, the corresponding total relaxation function is indicating that the curve fitting using the conventional KT function will yield a reduced linewidth (1 − Q)∆.
Conversely, in the limit of slow ion dynamics (ν i ≪ ∆), the change in relaxation function is independent of linewidth, and the exponential relaxation of the 1  3 component dominates: which is valid for Q ≈ 1.This suggests that it is difficult to distinguish the ion dynamics from muon diffusion in the slow exponential damping.As demonstrated in the case of FAPbI 3 below, the difference from the conventional KT function becomes notable in the intermediate region of ν i ∆, so the spectra in such a temperature region are important in data analysis to distinguish between ion dynamics and muon selfdiffusion.

D. Comparison with experimental results
One of the fields of µSR application that would greatly benefit from the possible distinction between local ion dynamics and muon diffusion would be in the study of ion diffusion.From the viewpoint of the present model, however, we are still unable to find convincing examples of µSR data that exhibit the characteristic behavior of ion diffusion revealed in our simulations.
For instance, an earlier µSR study on Li 0.73 CoO 2 reported that Li ions exhibit diffusive motion whereas µ + remains bound to an O atom and immobile at low temperatures [2].Above T ≈ 300 K, the linewidth obtained by KT function analysis (≡ ∆ KT ) shows a significant decrease with increasing T , which is attributed to the onset of muon diffusion in the literature.Since such a narrowing behavior is also expected to occur in the limit of large ν i (∆/ν i → 0) in our model [see Figs.3(a)-3(c)], let us examine the consistency of our model with the reported µSR result across the entire T range from 10 to 400 K.
From the value of ∆ KT extrapolated to T → 0 K and that to ∼400 K [2], the total linewidth ∆ and its static part √ 1 − Q∆ in the limit of ν i → ∞ for our model are respectively estimated to be 0.26 and 0.07 MHz, yielding Q = 1 − (0.07/0.26) 2 = 0.93.However, this value disagrees with that obtained by another evaluation, Q = 1 − (0.23/0.43) 2 = 0.71 derived from ∆ KT = 0.23 and 0.43 MHz respectively calculated for x = 0 and 0.75 in Li x CoO 2 at the presumed muon site adjacent to oxygen [2,27,28].Moreover, in our model with Q = 0.93, ZF-µSR time spectra should retain Gaussian-like features for any ∆/ν i value, as exemplified in Fig. 4. On the other hand, the spectrum observed at 275 K [27,28] exhibits a exponential lineshape similar to that for ∆/ν i < 0.5 in Fig. 3(d) (Q = 1, where ∆ = ∆ KT ), and the lineshape is suggested to become more exponential above 275 K as inferred from the decrease in ∆/ν i to ∼0.08 at 400 K. Thus, it is difficult to explain the behavior of the ZF-µSR spectra in Li 0.73 CoO 2 over the entire T range in a unified manner with our immobile-µ model, and it may be necessary to consider the effect of muon diffusion at lower temperatures, as has been discussed recently [29].
On the low Q side, our simulation indicates that the change in lineshape is strongly suppressed; as seen in Fig. 3(a) (Q = 0.2), the time spectral change is so small that it appears to correspond to quasi-static KT functions regardless of ∆/ν i .This still seems inconsistent with the significant time spectral change (and associated ν) in Li x CoO 2 .Besides this, such small deviations from the KT function for a small Q are also expected from other causes, e.g., a smaller number of nearestneighbor nuclear magnetic moments around a muon, making it more difficult to identify the true cause of the spectral change.In contrast, as discussed below, FAPbI 3 turns out to be an excellent example for testing our model because it suggests a large Q.It has been inferred from a previous µSR study that muons in FAPbI 3 are exposed to random local fields originating from the nuclear magnetic moments of both FA molecules and the PbI 3 lattice [14].Among these fields, the magnetic field from the FA molecule (mainly due to the proton magnetic moment) is expected to exhibit fluctuations due to the rotational motion induced by thermal excitation.The analysis of the µSR time spectra in the tetragonal phase of FAPbI 3 using the KT function showed that the linewidth ∆ appears to decrease monotonically from ∼50 to 140 K, whereas the fluctuation frequency ν exhibits a weak peak at around 100 K.Moreover, the global curve fitting including the spectra under various LFs are relatively poor in the 50-120 K range where ∆ changes significantly.These observations point to the situation that the actual relaxation function deviates from the KT function, as demonstrated in Fig. 3. Therefore, we focused on the temperature variation of the ZF-µSR spectrum and investigated whether the spectral changes be reproduced by using the appropriate Q, ∆, and ν i when the data were analyzed with the new model.Specifically, we prepared a multidimensional numerical table of G ID z (t; ∆, ν i , Q) for the parameters ∆ • t, Q, and ν i calcu- lated by Monte Carlo simulation [30], and we performed the TABLE I. Nuclear dipolar linewidth calculated for the muon sites inferred from DFT calculations, where ∆ is the total linewidth and ∆ PbI 3 is that without the contribution of FA molecules [14]."Fixed FA" and "VCA" correspond to those shown in Figs.1(c) and 1(d), respectively, where the FA molecular orientation was set to be common in the calculations of ∆.
least-squares curve fitting of the time-dependent asymmetry using the equation where A 0 is the initial asymmetry depending on the instrument and A b is the background signal from muons that missed the sample.In addition, a global fitting for multiple spectra was performed, where A 0 , ∆, and Q were common for the spectra at all temperatures.As shown in Fig. 4(a), all the spectra were well reproduced (reduced net chi-square = 1.33), and the values of the two common parameters were obtained as Q = 0.947(3) and ∆ = 0.191(1) MHz.
The theoretical evaluation of ∆ in the tetragonal phase of FAPbI 3 has a problem in that it is necessary to take into account the multiple orientations of FA molecules in the PbI 3 lattice, which makes the ab initio DFT calculations difficult in the Tetra-LT phase (< 140 K).Therefore, we discuss the results of DFT calculations for the Tet-HT phase (140-280 K), where the contribution of FA molecules to ∆ is quenched by the motional effect.More specifically, the hydrogen defect calculations (for mimicking muons) were conducted in two cases, one with the FA molecules fixed to one of the eight equivalent orientations and the other with all possible orientations taken into account using the virtual crystal approximation (VCA).(See Ref. 14 for the details of the DFT calculations.)The predicted muon sites in each case are shown in Figs.1(c) and 1(d), and the corresponding ∆ values are summarized in Table I.Although the agreement between the calculation and experimental results is not necessarily good for these absolute values, the result of VCA calculation is in reasonable agreement with experimental result for √ 1 − Q, which indicates the relative change in ∆ with and without the contribution of FA molecules.
More importantly, as shown in Fig. 4(b), we successfully demonstrated using µSR that ν i exhibits a tempera- ture dependence consistent with the thermally activated process in the temperature range where the rotational motion of FA molecules is strongly suggested from other experiments; this was not obtained directly by analysis using the KT function [14].Curve fitting using the Arrhenius function, ν i = ν 0 exp(−E a /k B T ) + ν c , yields excellent agreement with ν 0 = 3.6(6) × 10 8 s −1 , E a /k B = 630(18) K, and ν c = 1.62(9) × 10 5 s −1 [see the dashed line in Fig. 4(b)].A b is small (≃ 0.015 ≪ A 0 ) and nearly independent of temperature as expected, supporting the present model.

III. DISCUSSION
We have shown by Monte Carlo simulations that when the fluctuations of H(t) due to ion dynamics correspond to 0 < Q < 1 in Eq. ( 12), we encounter a ∆/ν i region where analysis using the KT function provides unsatisfactory results.This immediately points to the possibility that a careful analysis of the time spectra in the temperature region showing such ∆/ν i may enable us to determine the origin of the fluctuations.We then compared our simulation data with previous µSR experimental data for the Li ion battery material Li x CoO 2 and the hybrid perovskite FAPbI 3 .As a result, it was demonstrated in the latter that the lineshapes of the time spectra are coherently reproduced by the simulations and that a reasonable temperature dependence of the cation molecular fluctuation frequency is obtained, whereas it was suggested in the former that the experimental results exhibit a qualitatively different behavior from the simulations.
Historically, the possibility of applying µSR to the study of ion dynamics was first discussed for the diffusion of Li + ions in Li x CoO 2 [2], and then the application of µSR to other transition metal oxides containing Li has been attempted [3][4][5][6][7][8][9][10][11].In such studies, the µSR time spectra reflecting H(t) originating from the 7 Li and other nuclear magnetic moments (e.g., 59 Co in Li x CoO 2 ) were analyzed using the dynamical KT function, and it was commonly assumed that the obtained fluctuation rate ν directly corresponds to the jumping frequency of Li + ions.
However, a very recent paper reporting the results of an "operando" µSR experiment in the charging and discharging of Li x CoO 2 incorporated in a half battery suggests an alternative interpretation that the fluctuation rate ν (derived from the analysis using the KT function) reflects the contributions of both Li + and muon jump diffusion, and that ν divided by the square root of the Li-muon mass ratio ( m Li /m µ ≃ 7.9) is the jump frequency of Li [29].In this regard, it is worth mentioning the established fact that muons diffuse in crystalline solids via the thermally activated tunneling process at temperatures higher than about a fraction of the Debye temperature Θ D [31,32].Although the above interpretation is based on different muon sites inferred from simplistic DFT calculations [29], we need to consider the general possibility that muon diffusion via the thermally activated tunneling process prevails over ion diffusion in the relevant temperature range.
Another important factor in terms of identifying the origin of fluctuation is the temperature dependence of ∆ derived from the analysis using the KT function.Although the decrease in apparent ∆ with increasing temperature may indicate a contribution of ion dynamics to spin relaxation, muon diffusion is also expected in the temperature region of interest.Our simulation results show that when the muon diffusion is dominant (ν µ ≫ ν i ), the spin relaxation can be well repro- duced using the KT function.Therefore, when the spectra are well reproduced using the KT function even though ∆ ex-hibits a change (as seen in the case of Li x CoO 2 ), it is more likely that the cause of the change in ∆ is attributed to the issue of muon site(s) rather than that of ion dynamics; note that ∆ can vary with the muon position r µ in inhomogeneous materials.For example, muons may prefer to occupy atomic vacancies or voids to reduce zero-point energy, which has a significant impact on the behavior of light particles [33].Provided that ∆ = ∆(r µ ) is small for such defect sites, the mean ∆ will decrease as the probability of a muon reaching these sites increases with temperature in a diffusion-limited process.Such behavior has been actually observed, for example, in InGaZnO [34] and VO 2 [35].
In examining previously reported µSR studies of Li ion diffusion, there are many cases in which ∆ decreases with increasing temperature, but thus far there has been little discussion on such behavior from the aforementioned viewpoint.If the observed spin relaxation is well reproduced by the KT function, it is incompatible with the ion diffusion for which the autocorrelation function is assumed to be given by Eq. ( 12) with Q < 1.To determine whether the observed fluctuations are due to ion diffusion, it will be necessary to examine the reproducibility of µSR spectra using the model presented here with Q < 1.
Finally, as a more general possibility, note that Q may not necessarily be constant over the temperature range of interest.For example, there could be a situation where the diffusion motion of two kinds of ions (other than µ + ) occurs in stages as the temperature increases because of different activation energies.Developing a model that takes such a possibility into account is a future issue.

IV. CONCLUSION
We have performed Monte Carlo simulations to investigate how the residual correlations in the fluctuations described by the Edwards-Anderson type parameter Q affect the lineshape of the relaxation function observed in the µSR experiment.As a result, we found a significant deviation from the KT function for 0 < Q < 1.We then compared the simulated relaxation function, G ID z (t; ∆, ν i , Q), with the µSR spectra observed in FAPbI 3 in which internal field fluctuations associated with the rotational motion of FA molecules and the quasi-static internal field from the PbI 3 lattice are presumed to coexist.The curve fitting shows that the model can provide an appropriate description of the temperature-dependent lineshape with reasonable ∆, ν i and Q, which is in line with cation molecular motion.This result paves the way for the possibility of experimentally distinguishing fluctuations due to ion dynamics around muons from muon self-diffusion, which has been considered difficult.It also calls into question the attribution of internal field fluctuations entirely to ion diffusion for the µSR time spectra well reproduced by the KT function, suggesting the need for a major reexamination of the results of previous analyses.

FIG. 1 .
FIG. 1. (Color online) Schematic of time-dependent fluctuation of local field H(t) at muon sites for (a) muon self-diffusion and (b) ion diffusion, where ions are located at corners on a square lattice.The difference in the dipolar field from each ion (with the nuclear magnetic moment orientation shown by an arrow) is shown in different colors (or by brightness in gray scale).Although the combination of dipolar fields at muon sites changes randomly with muon diffusion in (a), the contributions from two ions at k and l sites remain unchanged, wheras those at i and j sites migrate (b).(c), (d) Structures of formamidinium lead iodate [FAPbI 3 , with FA denoting HC(NH 2 ) 2 ] in the tetragonal phase, and associated muon sites (areas indicated by arrows) suggested from first-principles DFT calculations (c) for the FA molecules fixed to one of the eight equivalent orientations and (d) for multiple FA orientations using VCA (see Sect.II D).Internal field fluctuations associated with the rotational motion of FA molecules and the static internal field from the PbI 3 lattice are presumed to coexist.

FIG. 2 .
FIG. 2. (Color online) (a) Schematic of local fields consisting of static and fluctuating components, where Q is the parameter representing the fractional amplitude of the local fields fluctuating at a frequency ν i .(b) Q dependence of the linewidths of the dynami- cal ( √ Q∆) and static ( √ 1 − Q∆) components corresponding to (a) [∆ = 1, see Eq. (12) in the main text].Dotted arrows correspond to the case for FAPbI 3 (Q = 0.947).

IDFIG. 3 .
FIG. 3. (Color online) Some examples of the muon spin relaxation function G ID z (t; ∆, ν i , Q) incorporating the effect of ion dynam- ics simulated using the parameter Q, where Q∆ 2 is the amplitude of fluctuation.(a)-(d) under zero field with different Q and ν i values, and (e)-(h) under various longitudinal fields B LF with Q = 0.5 and ∆/ν i = 0.5.Note that Q = 1 (and 0) corresponds to the conventional KT function.
approximated to

FIG. 4 .
FIG. 4. (Color online) (a) Zero-field µSR spectra observed at low temperatures in FAPbI 3 fitted with a relaxation function calculated by Monte Carlo simulation (global fitting assuming that the initial asymmetry, Q, and ∆ are common for all temperatures).All the spectra are reproduced by using Q = 0.947(3) and ∆ = 0.191(1) MHz.(b) Arrhenius plot of the internal field fluctuation frequency ν i at- tributed to the FA molecular motion obtained in this analysis.The dashed line is the fitting by the sum of the thermally activated function and the constant.Inset: temperature dependence of the constant background (see text for details).