Motional Narrowing under Markovian and Non-Markovian Hopping Transitions in Inhomogeneous Broadened Absorption Line Shape

TL;DR

This paper investigates motional narrowing effects in inhomogeneously broadened spectra, deriving exact solutions for spin correlation functions considering Markovian and non-Markovian transitions, and explaining experimental ESR/EPR line width minima.

Contribution

It provides an exact analytical model for motional narrowing incorporating static inhomogeneous broadening and non-Markovian effects, linking theory with experimental observations.

Findings

Exact solution for spin correlation functions considering both Markovian and non-Markovian processes.

Demonstration that line width minima occur as a function of inverse temperature under thermal activation.

Showing that Gaussian line shape narrowing does not necessarily produce Lorentzian line shapes, especially with heavy-tailed waiting time distributions.

Abstract

Inspired by recent experiments showing a minimum of electron paramagnetic resonance (ESR/EPR) line width as a function of inverse temperature, we studied the motional narrowing effect by considering a combined model of carrier transitions and static dispersion of the angular frequency giving rise to an inhomogeneous broadening in the spectrum. The dispersion of the angular frequency results from the distribution of the local field. The transition between the sites under inhomogeneous static local field induces adiabatic relaxation of the spin. We also considered the on-site inherent (nonadiabatic) relaxation of the spin. We obtained the exact solution of the spin correlation function by explicitly considering transitions between two sites for both Markovian and non-Markovian transition processes. The absorption line shape is expressed in terms of the Voigt function, which is a convolution of a Gaussian function and a Lorentzian function. Using the known properties of the Voigt function, we discuss the correlation between the change in the full-width at half-maximum and the change in line shape, both of which are induced by motional narrowing. By assuming thermal activation processes for both the hopping transition and the on-site inherent relaxation, we show that the minimum of the width appears as a function of inverse temperature as observed experimentally in organic materials.Contrary to the general belief, we also show that the narrowing of the Gaussian line shape under a local random field did not necessarily lead to a Lorentzian line shape in particular under the presence of heavy tail property in the waiting time distribution of hopping transitions.

Equations

$$\frac{d}{dt}{M}_a(t)=i{\omega }_a{M}_a(t)-{\gamma }_s{M}_a(t)-\frac{\gamma }{2}\left(M_a(t)-M_b(t)\right)$$

$$\frac{d}{dt}{M}_b(t)=i{\omega }_b{M}_b(t)-{\gamma }_s{M}_b(t)-\frac{\gamma }{2}\left(M_b(t)-M_a(t)\right),$$

$$⟨M(t)/M(0)⟩=\varphi (t)$$

$$I(\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }dt\varphi (t)\exp (-i\omega t).$$

$$\varphi (t)=\exp \left(-\frac{1}{4}{\Delta }_0^2{t}^2-2{\gamma }_{s}t\right)\sum _{m=0}^{\infty }\frac{(-1{)}^m}{m!}{\left(\frac{{\Delta }_{0}t}{2}\right)}^{2m}\mathcal{M}\left(m;2m+1;-\gamma t\right).$$

$${\varphi }_s(t)\approx \exp \left(-\frac{1}{2}{\Delta }_0^2{t}^2-2{\gamma }_{s}t\right),$$

$${\varphi }_f(t)\approx \exp \left(-\frac{1}{4}{\Delta }_0^2{t}^2-2{\gamma }_{s}t\right),$$

$${\varphi }_f(t)\approx \exp \left(-\frac{1}{4}{\Delta }_0^2{t}^2-2{\gamma }_{s}t\right)\left[1-\frac{1}{2}{\left(\frac{{\Delta }_0}{\gamma }\right)}^2\left(\gamma t-1+e^{-\gamma t}\right)\right],$$

$$1-\frac{1}{2}{\left(\frac{{\Delta }_0}{\gamma }\right)}^2\left(\gamma t-1+e^{-\gamma t}\right)\approx 1-{\left(\frac{1}{2}{\Delta }_{0}t\right)}^2,$$

$$\varphi (t)\approx \exp \left[-\frac{1}{4}{\Delta }_0^2{t}^2-2{\gamma }_{s}t-\frac{1}{2}{\left(\frac{{\Delta }_0}{\gamma }\right)}^2\left(\gamma t-1+e^{-\gamma t}\right)\right],$$

$$I(\omega )\approx V(\omega ;{\Delta }_0,2{\gamma }_s),$$

$$I(\omega )\approx V(\omega ;{\Delta }_0/\sqrt{2},2{\gamma }_s)$$

$$V(x;\sigma ,r)=\text{mbox}Re\left[\exp \left(-\frac{(x-ir{)}^2}{2{\sigma }^2}\right)\text{mbox}erfc\left(\frac{(ix+r)}{\sqrt{2}\sigma }\right)\right]/\sqrt{2\pi {\sigma }^2}.$$

$$V(x;\sigma ,r)={\int }_{-\infty }^{\infty }dy{P}_G(y;\sigma )P_L(x-y;r),$$

$$P_G(x;\sigma )$$

$$P_L(x;r)$$

$$\Omega =\frac{\omega }{{\Delta }_0},{\Gamma }_s=\frac{2{\gamma }_s}{{\Delta }_0}\text{mbox}and\Gamma =\frac{\gamma }{{\Delta }_0}.$$

$$\varphi (\tau )\approx \sum _{j=0}^{\infty }\frac{1}{j!}{\left(-\frac{1}{2{\Gamma }^2}\right)}^j\exp \left[-\frac{1}{4}{\tau }^2-\left({\Gamma }_s+\frac{1}{2\Gamma }+\Gamma j\right)\tau +\frac{1}{2{\Gamma }^2}\right].$$

$$I(\Omega )=\sum _{j=0}^{\infty }\frac{1}{j!}{\left(-\frac{1}{2{\Gamma }^2}\right)}^j\exp \left(\frac{1}{2{\Gamma }^2}\right)V\left(\Omega ;\frac{1}{\sqrt{2}},{\Gamma }_s+\frac{1}{2\Gamma }+\Gamma j\right).$$

$$f_{Vw}$$

$$f_{Vn}$$

$$f_V$$

$$\frac{\partial }{\partial x}V(x;\sigma ,r)$$

$$=\frac{-xV(x;\sigma ,r)+rL(x;\sigma ,r)}{{\sigma }^2},$$

$$L(x;\sigma ,r)=\text{mbox}Im\left[\exp \left(-\frac{(x-ir{)}^2}{2{\sigma }^2}\right)\text{mbox}erfc\left(\frac{(ix+r)}{\sqrt{2}\sigma }\right)\right]/\sqrt{2\pi {\sigma }^2}.$$

$$I^{\prime }(\Omega )$$

$$\frac{-\Omega V\left(\Omega ;\frac{1}{\sqrt{2}},{\Gamma }_s+\frac{1}{2\Gamma }+\Gamma j\right)+({\Gamma }_s+\frac{1}{2\Gamma }+\Gamma j)L\left(\Omega ;\frac{1}{\sqrt{2}},{\Gamma }_s+\frac{1}{2\Gamma }+\Gamma j\right)}{{\sigma }^2}.$$

$$I^{\prime }(\omega )$$

$$\psi (t)={\int }_0^{\infty }dEg(E)\gamma (E)\exp \left(-\gamma (E)t\right)\sim \frac{\alpha \Gamma \left(\alpha +1\right)}{{\gamma }_{\mathrm{r}}^{\alpha }{t}^{\alpha +1}},$$

$$\hat{\psi }(s)\sim 1-\left[\pi \alpha /\sin (\pi \alpha )\right]{\left(s/{\gamma }_{\mathrm{r}}\right)}^{\alpha }\text{mbox}fors/{\gamma }_{\mathrm{r}}<1.$$

$$\varphi (t)=\exp \left(-\frac{1}{4}{\Delta }_0^2{t}^2-2{\gamma }_{s}t\right)\left[1+\sum _{m=1}^{\infty }\frac{(-1{)}^m}{m!}{\left(\frac{{\Delta }_{0}t}{2}\right)}^{2m}\frac{(2m)!}{(m-1)!}{E}_{\alpha ,2m+1-\alpha (m-1)}^{(m-1)}\left(-{\gamma }_{\alpha }{t}^{\alpha }\right)\right].$$

$${\varphi }_s(t)\approx \exp \left(-\frac{1}{2}{\Delta }_0^2{t}^2-2{\gamma }_{s}t\right).$$

Full text

I Introduction

Electron paramagnetic resonance (ESR/EPR) is a powerful tool in the study of charge transport in organic electronic materials because ESR selectively detects unpaired conducting electrons and is insensitive to bounded paired electrons. The ESR technique has been applied to reveal molecular orientation with respect to the substrate, carrier doping effects, delocalized carrier states, and carrier transport properties.Namatame et al. (2017); Marumoto et al. (2005, 2006); Tanaka et al. (2006); Marumoto et al. (2011); Matsui et al. (2010, 2012, 2008) Recently, motional narrowing was observed in organic materials such as pentacene and rubrene by varying temperature; Matsui et al. (2008, 2012, 2010); Marumoto et al. (2011); Tanaka et al. (2006) the width of the line shape shows a minimum as a function inverse temperature. Matsui et al. (2008, 2012, 2010) In some instances, this motional narrowing was associated with an absorption spectral change from a Gaussian absorption line shape to a Lorentzian line shape but, interestingly, it was sometimes observed without an apparent change in the line shape.Matsui et al. (2008, 2012, 2010); Marumoto et al. (2011); Tanaka et al. (2006) The motional narrowing has been investigated theoretically, mainly in the time domain, including the analytical work of Kubo and Anderson;Kubo (1954); Kubo and Tomita (1954); W. Anderson (1954); Kubo (1969) the asymptotic time dependence of the spin correlation function has been investigated but the complete kinetics has received less attention. In ESR, motional narrowing is measured using the spectrum of the first derivative of the absorption line shape. Although motional narrowing has been known for decades,Slichter (1990); Bloembergen et al. (1948) its effect on the absorption line shape has not been fully determined compared with the magnetic relaxation in the time domain.

Qualitative features of motional narrowing can be understood by analyzing the Kubo–Anderson model.Kubo (1954); Kubo and Tomita (1954); W. Anderson (1954); Kubo (1969) In this model, the influence of the local field fluctuation is assumed to involve Gaussian and Markovian random processes.Kubo (1954); Kubo and Tomita (1954); W. Anderson (1954); Kubo (1969) The fluctuating field is thought to reflect the change in the magnetic field imposed on the spin by transport in an inhomogeneous medium. In the limit of slow decay times of the local field fluctuations, a Gaussian decay of the spin correlation was obtained, whereas in the limit of fast decay, an exponential decay of the spin correlation was obtained. In this stochastic description, the change in the decay of the spin correlation function corresponds to motional narrowing of the absorption line shape. The Gaussian decay and the exponential decay lead to a Gaussian and a Lorentzian absorption line shape, respectively. Qualitatively different line shapes are expected. Nevertheless, we find the more reserved statement ”narrowing process of the Gaussian spectrum in the presence of Markovian motion” in the original paper by Kubo;Kubo (1954) the nature of the narrowed spectrum was not specified.

It has been recognized that the fluctuating field does not obey Gaussian and Markovian random processes when field fluctuations originate from the movement of a magnetic carrier in inhomogeneous magnetic-field environments.Czech and Kehr (1984, 1986); Mitra and Le Doussal (1991) In these studies, an exponential decay is obtained as the long-time asymptotic decay of the spin correlation regardless of the dimensionality of the kinematic space. Although the asymptotic decay was exponential as predicted assuming both Gaussian and Markovian random processes, the relaxation rate constant differed from the results obtained under this assumption in particular for low dimensions. These studies pointed out the importance of taking into account the actual movement of spin carriers in inhomogeneous local-field environments. Moreover, absorption line shapes have received less attention in theoretical work.

In this paper, we take into account explicitly the motion of the spin carrier under random local fields. For simplicity, we consider two-site transition processes. The local field at each site obeys a Gaussian distribution. We obtain the exact analytical expression for the spin correlation function for both Markovian and non-Markovian transition processes. Using the exact expression for Markovian transition processes, the spin correlation function is expressed by the function similar to that given by the familiar expression of the Kubo–Anderson model.Kubo and Tomita (1954); Kubo (1954) The slight difference manifests itself when the transition rate between sites is large. Therefore, the narrowed absorption line shape from a fast transition between sites may differ from the Lorentzian line shape. The absorption line shape is expressed in terms of the Voigt function, which is a convolution of Gaussian and Lorentzian functions. Using the known properties of Voigt function,Olivero and Longbothum (1977) we discuss the correlation between the change in the full-width at half-maximum (FWHM) and the change in line shape, both of which are induced by motional narrowing. To confirm the above-mentioned results, we also present the spectrum of the first derivative, which is relevant to ESR measurements. By using the above mentioned results, we show that the minimum of the width appears as a function of inverse temperature as observed experimentally in organic materials. We also show that the narrowing of the Gaussian line shape does not necessarily lead to a Lorentzian line shape under a random local field.

II Theory

Introduced by Kitahara,W. Kehr and Kitahara (1987); W. Kehr et al. (1987); Kitahara (1994) the model comprises a spin executing transitions between a pair of statistically equivalent sites labeled by a and b in Fig. 1. We consider a large number of the pairs and study relaxation of spins by transitions between two sites with local field at each site. We study the ensemble average of spins when a large number of sites are exposed to different local field environments. The local field at each site obeys a Gaussian distribution (Figure 1).

Adiabatic relaxation occurs via spin transitions under an inhomogeneous local magnetic field; the time correlation of spins results from recursive nature of spin transitions even though local fields are spatially uncorrelated, i.e. spin correlations are generated by returning to the originally occupied site where the time independent and random magnetic field is assigned. The ensemble average of the magnetic moments of the carriers precesses under different local fields and decays by Markovian transitions. In addition, we consider on-site inherent (nonadiabatic) relaxation characterized by the decay rate constant denoted by $\gamma_{s}$. Within linear response theory, the width of the absorption lines is characterized by the transverse relaxation time.Kubo et al. (1991) Therefore, we consider transverse relaxation as in the Kubo–Anderson model,Kubo (1954, 1969); Anderson and Weiss (1953); W. Anderson (1954) where the external magnetic field is applied in the z-direction in the Cartesian coordinate system that is erected. For the spin labeled by the site $i=a\mbox{ or }b$, we introduce $M_{i}=M_{i,x}+iM_{i,y}$ in this coordinate system. We denote by $\omega_{i}$ the angular frequency of the spin rotating around a fixed magnetic field in the z-direction. We assume that a Gaussian distribution of the magnetic fields at all sites; that is, $\omega_{a}$ and $\omega_{b}$ are statistically independent and obey Gaussian distributions. The mean of $\omega_{i}$ is zero and the variance is denoted by $\Delta_{0}^{2}$. The transition rate between the sites is denoted by $\gamma$. When a static magnetic field is applied along the z-axis, the magnetization $M_{i}^{(z)}(t)$ undergoes a precession around the field direction with Larmor angular frequency $\omega_{0}$. Because we are interested in the dispersion of the magnetization around the coherent motion with the Larmor angular frequency, we study the magnetization given by $M_{i}(t)=M_{i}^{(z)}(t)\exp(i\omega_{0}t)$. The kinetic equation for the spins is written

[TABLE]

where we consider the initial conditions at equilibrium given by $M_{a}(0)=M_{b}(0)=M_{0}/2$.

We calculate the spin correlation averaged over the distribution of the local Larmor angular frequency and the initial condition

[TABLE]

and the line-shape function using the Wiener–Khinchin theorem,

[TABLE]

The term $\langle\cdots\rangle$ in Eq. (3) can be expressed as $1/(2\pi\Delta_{0}^{2})\int_{-\infty}^{\infty}d\,\omega_{a}\int_{-\infty}^{\infty}d\,\omega_{b}\exp\left[-(\omega_{a}^{2}+\omega_{b}^{2})/(2\Delta_{0}^{2})\right]\cdots$. The power spectrum gives the absorbance of the microwave radiation multiplied by $\omega_{0}$ if $\omega_{0}\gg\omega_{i}$ is satisfied for the majority of the distribution of $\omega_{i}$.Kubo and Tomita (1954) The ESR absorption is essentially proportional to the power spectrum of the autocorrelation function of $M(t)$.

As shown in the Appendix A, the spin correlation function becomes

[TABLE]

For a slow transition rate expressed by $\gamma t\ll 1$, taking the limit simplifies Eq. (5) to

[TABLE]

where $\sum_{m=0}^{\infty}(-1)^{m}\left(\Delta_{0}t/2\right)^{2m}/m!=\exp\left(-\Delta_{0}^{2}t^{2}/4\right)$ is used. In contrast, for a fast transition rate expressed by $\gamma t\gg 1$, taking the limit reduces Eq. (5) to

[TABLE]

where the asymptotic expansion given by ${\cal M}\left(m;2m+1;-\gamma t\right)\approx\Gamma(2m+1)/[m!(\gamma t)^{m}]$ is used. This equation indicates a slowing down of the spin relaxation in the fast-transition-rate limit because the factor given by $(1/2)\Delta_{0}^{2}$ in Eq. (6) changes to $(1/4)\Delta_{0}^{2}$ in Eq. (7).Kitahara (1994); W. Kehr and Kitahara (1987); W. Kehr et al. (1987) Correspondingly, the absorption spectral line is narrowed.Kitahara (1994); W. Kehr and Kitahara (1987); W. Kehr et al. (1987) The result of Eq. (7) can be interpreted by averaging $\exp(i\omega_{m}t)$ using the distribution of $\omega_{m}$ given by Eq. (A.7). Here, $\omega_{m}$ denotes the mean angular frequency between the sites. The result implies that the spin precesses at the mean angular frequency between the sites in the fast-transition-rate limit. The limit was also found for quantum tunneling.W. Kehr and Kitahara (1987); W. Kehr et al. (1987) By further calculating an ensemble average with respect to the angular frequency, the relaxation becomes a Gaussian function. In the fast-transition-rate limit, the spin precesses at the mean angular frequency between the sites due to recursive nature of transitions. Recursive nature of transitions might still remain for spatially extended systems, in particular in low-dimensional systems; the result could depend on the lattice structure. We also reserve the possibility that the spin relaxation function given by Eq. (7) could result from confinement effect; the recursive feature of transitions could be enhanced by confinement in a finite volume.

A correction to Eq. (7) may be calculated from Eq. (5), the result being

[TABLE]

where the summation in Eq. (5) is taken up to $m=1$; ${\cal M}\left(0;1;x\right)=1$ and ${\cal M}\left(1;3;-x\right)=2[x-1+\exp(-x)]/x^{2}$ are introduced. Using the approximation given by

[TABLE]

we expand Eq. (8) and find $1-\frac{1}{2}\Delta_{0}^{2}t^{2}$. This expression coincides with the short time expansion of Eq. (6). On the basis of this observation, Eq. (8) is modified to

[TABLE]

which, given in terms of $\gamma$, makes the expression applicable to longer time durations regardless of the limit. The result of Eq. (10) is close to the numerically exact result [Figure 2(a)]. In the absence of the $-(1/4)\Delta_{0}^{2}t^{2}$–term and $\gamma_{s}$, the right-hand side of Eq. (10) is a well-known function describing the narrowing of the Gaussian profile through increasing $\gamma$.Kubo and Tomita (1954); Kubo (1954); W. Anderson (1954) In the absence of the $-(1/4)\Delta_{0}^{2}t^{2}$ term and $\gamma_{s}t$, the functional form was regarded as the basis for understanding motional narrowing of line shapes. Here, Eq. (10) including the $-(1/4)\Delta_{0}^{2}t^{2}$ term was obtained in an approximation by considering a combined model of carrier transitions and a static dispersion of the angular frequency that gave rise to inhomogeneous broadening in the spectrum. The exact result for the model is given by Eq. (5), which along with Eq. (10), constitute in part the main results of this paper.

III Spectrum

The spectrum is obtained numerically by substituting Eq. (5) into Eq. (4); we also obtained analytical expressions. We have

[TABLE]

in the limit of small value of $\gamma$, and

[TABLE]

in the limit of large value of $\gamma$. In both limits, the spectrum is expressed by the Voigt function,Olver et al. (2010)

[TABLE]

Denoting the complex conjugate of a function $f$ by $f^{*}$, we have $Re(fg)=Re(f^{*}g^{*})$. Therefore, Eq. (13) can be rewritten as $V(x;\sigma,r)=\mbox{Re}[\exp(-z^{2})\mbox{erfc}(-iz)]/\sqrt{2\pi\sigma^{2}}$, where $z=(x+ir)/(\sqrt{2}\,\sigma)$. $w(z)=\exp(-z^{2})\mbox{erfc}(-iz)$ is called the Faddeeva function.Olver et al. (2010) As a convolution of a Gaussian profile $P_{G}(x)$ and a Lorentzian profile $P_{L}(x)$, the Voigt profile is

[TABLE]

where the Gaussian and Lorentzian profiles are given by

[TABLE]

respectively. The Voigt function describes a combined broadening of the Gaussian and Lorentzian line shapes.

For the numerical calculation, we introduce dimensionless parameters using $\Delta_{0}$ as

[TABLE]

We systematically vary these parameters to study the absorption line shape. Dimensionless parameters are also used to obtain an approximate expression valid for a wide range of $\gamma$. We first express Eq. (10) using $\tau=\Delta_{0}t$ as

[TABLE]

By substituting Eq. (18) into Eq. (4), we obtain the approximate expression

[TABLE]

The power spectra were expressed in terms of the weighted sum of the Voigt distribution, where the weight is expressed as a power of $-(1/2\Gamma^{2})$. As shall be shown later [Figure 2(b)], the approximate expression reproduces the numerically exact result by summing the first four terms as long as the condition $\Gamma>1/\sqrt{2}$ is satisfied. Equation (19) and the simplified expressions given by Eqs. (11) and (12) contribute to the main results of this paper.

In Fig.  2(a), we compare the exact result obtained from Eq. (5) and the approximate result obtained from Eq. (10). The results are shown for $\Gamma_{s}=0.01$; for a fixed value of $\Gamma$, the approximation is better as $\Gamma_{s}$ increases. We also present the result for a short-time expansion given by Eq. (8). We find that the degree of the approximation improves by modifying Eq. (8) as in Eq. (10). In Fig.  2(b), we compare the numerically exact spectrum obtained using Eq. (5) with Eq. (4) and the approximate one obtained using Eq. (19) by summing the first four terms. The results are shown again for $\Gamma_{s}=0.01$; the approximation is better as $\Gamma_{s}$ increases for the same value of $\Gamma$. When $\Gamma$ is varied, the convergence of the sum is obtained for $\Gamma>1/\sqrt{2}$ regardless of the values of $\Gamma_{s}$.

IV Results of FWHM

We study the FWHM and the equivalent width as (or the integral breadth) indicated in Fig.  3. The equivalent width can be found by drawing a rectangle with a height equal to that of the absorption line shape and the area equal to that under the absorption line; the width of the rectangle thus specified is called the equivalent width. When the area is normalized, the equivalent width is equal to the inverse of the height. The kurtosis parameter (the inverse of the form factor) of the spectrum, defined by the equivalent width divided by the FWHM, Langford (1978) reflects the spectral shape variation by motional narrowing. The kurtosis parameter value is larger when the kurtosis of the peak is larger. The kurtosis parameter of the Lorentzian spectrum is $\pi/2\approx 1.57$ and that of the Gaussian spectrum is $\sqrt{\pi/\log 2}/2\approx 1.06$. The kurtosis parameter value of the Lorentzian spectrum is larger than that of the Gaussian spectrum.

Using the empirical equation given by $r+\sqrt{r^{2}+8\ln(2)\sigma^{2}}$ to express the FWHM of the Voigt function $V(x;\sigma,r)$,Olivero and Longbothum (1977) the FWHM of the Voigt function given in Eq. (11) for small $\gamma$ is written

[TABLE]

From Eq. (12), the FWHM of the Voigt function given in Eq. (12) for large $\gamma$ is narrower than the above result and is expressed as

[TABLE]

In Fig.  4, we show the FWHM normalized by $\Delta_{0}$ as a function of the dimensionless transition rate. The exact numerical results were calculated by substituting Eq. (10) into Eq. (4) and performing a numerical integration. Judging from Fig.  4, we note that the fraction of narrowing effect is larger when $\Gamma_{s}$ is smaller. When $I(\Omega)$ is dominated by the first term given by $j=0$ of Eq. (19), the FWHM can be approximated by

[TABLE]

The (red) thin lines indicate the results of the approximate solution given by Eq. (22). We find that the approximation is satisfactory for $\Gamma>\mbox{Max}(\Gamma_{s},1)$. The (red) circles represent the FWHM calculated numerically using Eq. (19) by summing the first four terms. As long as $\Gamma>1/\sqrt{2}$ is satisfied, the red circles approximate the exact numerical results. Because the integration of the spectrum over the angular frequency is normalized by $1$, the equivalent width can be obtained from $1/I(0)$.

We define the kurtosis parameter by the equivalent width divided by the FWHM. Langford (1978) In Fig. 5, the kurtosis parameter is shown against $\Gamma$ for various values of $\Gamma_{s}$. The variation of the kurtosis parameter values against $\Gamma$ is larger as the value of $\Gamma_{s}$ is smaller. For $\Gamma_{s}=10$, we find a monotonic increase in the kurtosis parameter as $\Gamma$ increases; the line shape approaches a Lorentzian line shape through motional narrowing. For $\Gamma_{s}=0.1$, we find the maximum value of the kurtosis parameter as a function of $\Gamma$. The appearance of the maximum and the relatively small kurtosis parameter value at the maximum indicate that the narrowing is not the qualitative change in the spectrum from the Gaussian function to the Lorentzian function but the Gaussian-like shape is retained during the narrowing process.

In Fig.  6(a), we present the temperature variation of the spectral width. We assume that the on-site inherent (nonadiabatic) relaxation is due to thermal activation with $\Gamma_{s}=c_{s}\exp[-E_{as}/(k_{\rm B}T)]$,Orbach and Blume (1962) where $E_{as}$ denotes the activation energy, $k_{\rm B}$ the Boltzmann constant, $T$ the temperature, and $c_{s}$ a parameter. For a given hopping transition rate, we assume thermal detrapping from a trap given by $\Gamma=c\exp[-E_{a}/(k_{\rm B}T)]$Kehr and Wichmann (1996) with $E_{a}$ the activation energy and $c$ a parameter. For another hopping transition rate, we assume the Marcus equation described by $\Gamma=c\sqrt{T_{m}/T}\exp[-\lambda/(4k_{\rm B}T)]$,Marcus (1956, 1964); Levich and Dogonadze (1956); Hol (1959); Holstein (1959) with $\lambda$ the reorganization energy and $T_{m}$ a parameter. In regard to the temperature dependence, the difference between these two transition rates turns out to be small (see Fig.  6). As the temperature is increased from a low temperature limit, the minimum of the width appears as a result of motional narrowing. In the high temperature regime, the increase in width obtained by increasing the temperature originates from the temperature dependence of the on-site relaxation rate constant denoted by $\Gamma_{s}$. The similar minimum of the FWHM as a function of inverse temperature is observed experimentally in organic materials.Matsui et al. (2008, 2012, 2010) In Fig.  6(b), we show the corresponding values of the kurtosis parameter as a function of inverse temperature. Judging from these values, the line shape is closer to a Gaussian profile than a Lorentzian profile over the entire temperature range. In the above analysis, the on-site inherent (nonadiabatic) relaxation rate was kept small within the temperature variation and this is the reason that the near-Gaussian profile is retained. If the on-site inherent (nonadiabatic) relaxation rate is large, as characterized by $\Gamma_{s}\gg 1$, the narrowing of a near-Lorentzian line shape can be obtained judging from Figs. 4-5. For such parameters, the profile is closer to being Lorentzian but the variation of the FWHM becomes smaller, as indicated in Fig.  4.

V First derivative of the power spectrum

When the recorded ESR spectra are the first derivatives of the absorption spectra, it is more convenient to calculate the first derivative of the spectrum to compare directly with the measured spectrum. This first derivative is calculated using the derivative of the Voigt function given byOlver et al. (2010)

[TABLE]

where the imaginary Voigt distribution function is given by,Olver et al. (2010)

[TABLE]

Using Eq. (19), the first derivative of the power spectrum is obtained

[TABLE]

In the limit of small and large values of $\gamma$, the first derivative of the power spectrum is given by

[TABLE]

respectively, where we have used Eqs. (11) and (12). In Fig.  7, we plotted the first derivative of the power spectrum when the Gaussian spectrum is partially narrowed. To avoid the effect of on-site nonadiabatic relaxation rate constant given by $\Gamma_{s}$, we set a small value for $\Gamma_{s}$. The value of $\Gamma$ is chosen so that partial narrowing of the Gaussian spectrum occurs (see Fig.  4). We find that the first derivative of the power spectrum calculated from Eq. (26) is close to the first derivative of a Gaussian profile compared with that of a Lorentzian profile. The result is consistent with the low kurtosis parameter value given by $1.11$.

VI Generalization to disordered systems

So far, we have considered Markovian transitions. Non-Markovian transitions have been considered for the Kubo-Anderson model using waiting time distribution functions. Jung et al. (2002); Loring and Mukamel (1985); Šanda and Mukamel (2005) Here, we take into account the averaging with respect to local field at each site in addition to Non-Markovian transitions.

In the Markovian model, we have considered the case that the transitions between the sites are characterized by a rate constant associated with a single activation energy. When the activation energy is distributed, the waiting time distribution for the transitions becomes non-exponential. We consider the waiting time distribution of spin at site a for the transition to site b denoted by $\psi(t)$. For simplicity, we consider the same waiting time distribution of spin at site b. Algebraic asymptotic tails could result from the distribution of the activation energy. Pfister and Scher (1978); Jakobs and Kehr (1993) When the activation energy distribution function is the exponential function with the energy depth given by $k_{\rm B}T_{\rm c}$, Pfister and Scher (1978) $g(E)=\exp\left[-E/(k_{\rm B}T_{\rm c})\right]/(k_{\rm B}T_{\rm c})$, and the activated release rate is assumed for transition between the sites, $\gamma(E)=\gamma_{\rm r}\exp\left[-E/(k_{\rm B}T)\right]$, the waiting time distribution is expressed as Pfister and Scher (1978)

[TABLE]

where $\alpha\equiv T/T_{\rm c}$ and $\Gamma(z)$ is the Gamma function. Olver et al. (2010) The waiting time distribution shows the power law decay with the absolute value of the exponent given by $\alpha+1$. Owing to the activation energy distribution, the waiting time distribution exhibits heavy tail property. By decreasing the temperature, the absolute value of the exponent becomes lower. The effect of heavy tail property on the diffusion has been studied by the continuous time random walk model using the waiting time distribution function for transition to a neighboring site. Montroll and Weiss (1965) According to the continuous time random walk model, the exponent $\alpha$ characterizes the subdiffusive growth of the mean square displacements given by $t^{\alpha}$, where $\alpha=1$ for the normal diffusion and $\alpha<1$ for the dispersive diffusion. In the small $s$ limit, the Laplace transform of the waiting time distribution function is given by Pfister and Scher (1978)

[TABLE]

The asymptotic expression will be used below. In this section, we consider the waiting time distribution function having an algebraic asymptotic tail and being normalized; $\int_{0}^{\infty}d\,t\psi(t)=1$.

As before, we calculate the spin correlation averaged over the distribution of the local Larmor angular frequency and the initial condition, where the even site occupation probability is assumed. As shown in the Appendix B, the spin correlation function becomes

[TABLE]

This is the exact solution and is generalization of Eq. (5) when the waiting time distribution has heavy tail property due to, e.g., the distribution of the activation energy. Equation (31) is the main result of this section.

For a slow transition rate expressed by $\gamma_{\alpha}t^{\alpha}\ll 1$, Eq. (31) can be simplified by introducing the lowest order series expansion given by $E_{\alpha,3}\left(-\gamma_{\alpha}t^{\alpha}\right)\approx 1/2$ as

[TABLE]

In contrast, for a fast transition rate expressed by $\gamma_{\alpha}t^{\alpha}\gg 1$, taking the limit reduces Eq. (31) to

[TABLE]

As explained below Eq. (7), this equation indicates a slowing down of the spin relaxation and narrowing of the absorption spectral line. It should be noted that the condition $\gamma_{\alpha}t^{\alpha}\gg 1$ could be hardly attained when the value of $\alpha$ is lowered; the results suggest that the motional narrowing can be hardly observed when the hopping motion is sluggish in disordered media.

A correction to Eq. (33) may be calculated from Eq. (31), the result being

[TABLE]

where the summation in Eq. (31) is taken up to $m=1$. As before Eq. (34) can be modified to

[TABLE]

By setting $\alpha=1$ and using $E_{1,3}\left(-x\right)=[x-1+\exp(-x)]/x^{2}$, Eq. (35) reduces to Eq. (10), where $\gamma=\gamma_{1}$. By applying the asymptotic expansion expressed by

[TABLE]

the asymptotic expression of Eq. (35) can be given by

[TABLE]

The last term in Eq. (37) indicates the exponential decay only when $\alpha=1$. For $\alpha<1$, non-exponential decay can be obtained owing to heavy tail property in the waiting time distribution. The exponent approaches to $2$ as $\alpha$ decreases; the exponential decay for $\alpha=1$ changes to the Gaussian-like decay as $\alpha$ decreases. Therefore, the corresponding power spectrum will deviate from Lorentzian form when $\alpha<1$.

In Eq. (35) [See also Eqs. (A.4) and (A.8)], $\langle\phi_{r}(t)\rangle_{\delta}$ is approximated by

[TABLE]

It can be rewritten as

[TABLE]

where $\psi(t)$ corresponds to the correlation function when the fluctuating angular frequency is modeled as the Gaussian-Markovian fluctuation. $\psi(t)$ is given by

[TABLE]

where the Mittag-Leffler function $E_{\alpha}(-\gamma_{\alpha}t^{\alpha})=E_{\alpha,1}(-\gamma_{\alpha}t^{\alpha})$ is defined through the inverse Laplace transformation as $E_{\alpha}(-\gamma_{\alpha}t^{\alpha})={\cal L}^{-1}\left[1/\left(s+\gamma_{\alpha}s^{1-\alpha}\right)\right]$. Initial decay obeys the stretched exponential form $\psi(t)=\exp\left[-\gamma_{\alpha}t^{\alpha}/\Gamma(1+\alpha)\right]$, which suggests that the Mittag-Leffler function is the natural generalization of the exponential correlation function. Indeed, $E_{\alpha}(-\gamma_{\alpha}t^{\alpha})$ represents a solution of the fractional relaxation equation given by $(\partial/\partial t)\psi(t)=-\gamma_{\alpha}\mbox{ }_{0}D_{t}^{1-\alpha}\psi(t)$, Metzler and Klafter (2000) where the Riemann-Liouville fractional derivative is defined through the convolution to an arbitrary function $f(t)$ as $\mbox{ }_{0}D_{t}^{1-\alpha}f(t)=[1/\Gamma(\alpha)](\partial/\partial t)\int_{0}^{t}d\,t_{1}f(t_{1})/\left(t-t_{1}\right)^{1-\alpha}$. The conventional exponential relaxation function is obtained for $\alpha=1$. On the other hand, the asymptotic expansion given by $\psi(t)=\left(\Gamma(1-\alpha)\gamma_{\alpha}t^{\alpha}\right)^{-1}$ leads to $\langle\phi_{r}(t)\rangle_{\delta}\approx\exp\left[-(\Delta_{0}^{2}t^{2-\alpha})/\left\{2\Gamma(3-\alpha)\gamma_{\alpha}\right\}\right]$, which shows intrinsic non-exponential decay when $\alpha<1$.

The conclusion of this section is as follows. In disordered media, power-law waiting time distribution is obtained when the activation barrier for hopping is exponentially distributed. The initial decay of the spin correlation function was Gaussian function of time and was not influenced by the presence of heavy tail property in the waiting time distribution of hopping transitions. The asymptotic decay was influenced by the presence of heavy tail property in the waiting time distribution; by increasing heavy tail property in the waiting time distribution function, the exponential-like decay component influenced by the hopping frequency approaches to the Gaussian-like decay. The results indicate that the motional narrowing to Lorentzian absorption line shape could be hardly observed in disordered media, though the narrowing of Gaussian absorption line shape while retaining the Gaussian-like line shape could be possible. Moreover, the on-set time of slowing down of spin relaxation increases with decreasing the absolute value of the exponent of the power-law waiting time distribution; motional narrowing phenomena could be hardly observed when $\alpha\ll 1$.

VII Conclusion

We studied motional narrowing of spins using the Kubo–Anderson–Kitahara model. In the model, the transition between two sites and the distribution of local field on each site were explicitly taken into account. The local field distribution of spins was assumed to be Gaussian. The transition between the sites having an inhomogeneous local field induces an adiabatic relaxation of the spin. We also considered the on-site inherent (nonadiabatic) relaxation of the spin.

We obtained the exact solution of the spin correlation function. [Equation (5)] We also derived an approximate expression given by Eq. (10); the spin correlation function was expressed by a well-known function describing the narrowing of a Gaussian profile with a multiplicative correction term arising from the transition between sites under inhomogeneous local field environments.Kubo and Tomita (1954); Kubo (1954) Because the spin correlation function was expressed in the form similar to the well-known function derived by assuming a Gaussian–Markovian fluctuation of the field, the narrowing process caused by transitions between sites having an inhomogeneous local field can be understood based on Gaussian–Markovian field fluctuations. More importantly, some deviations from this simplified view are also clearly recognized. The additional term $\exp\left[-(1/4)\Delta_{0}^{2}t^{2}\right]$ can be regarded as an ensemble average of the mean angular frequency between the sites. The effect of carrier transitions between the sites having distributed angular frequency gives rise to inhomogeneous broadening in the ESR spectrum. When the distribution of the on-site local field is Gaussian, the relaxation is given by a Gaussian function. As we discussed below Eq. (7), recursive nature of transitions manifest itself in the spin correlation function in the fast-transition-rate limit. For extended systems, recursive nature of transitions should depend on the dimensionality and the lattice structure. Unfortunately, the exact analytical treatment is exceedingly difficult for spatially extended systems.

The motional narrowing has been studied using the absorption spectrum rather than the spin correlation function. We derived the approximate expression of the absorption spectrum in terms of the Voigt function. The Voigt function, as a Gaussian–Lorentzian convolution, is suited to describing the combined effect of Gaussian line-shape and the Lorentzian line-shape broadening. Using an empirical equation describing FWHM of the Voigt function,Olivero and Longbothum (1977) we examined the motional narrowing for a wide range of values of the on-site inherent (nonadiabatic) relaxation rate constant. We also analyzed the kurtosis parameter defined using the equivalent width and the FWHM. The value of the kurtosis parameter increases as the Gaussian line shape changes to a Lorentzian line shape.

When the on-site inherent (nonadiabatic) relaxation rate is small, the Gaussian line shape changes to a narrowed Gaussian line shape by increasing the transition rate. The result is confirmed by studying the first derivative of the absorption spectrum. The kurtosis parameter shows a maximum deviation from the Gaussian line shape in the course of narrowing.

When the on-site inherent (nonadiabatic) relaxation rate is large, the kurtosis parameter increases monotonically by increasing the transition rate; the line shape approaches a Lorentzian line shape through motional narrowing.

The fraction of the FWHM variation plotted against its maximum value decreases by increasing the on-site inherent (nonadiabatic) relaxation rate. Therefore, the narrowing effect is larger by decreasing the on-site inherent (nonadiabatic) relaxation rate. By assuming activation type for both adiabatic and nonadiabatic relaxation rates, the FWHM shows a minimum as a function of inverse temperature; the result is consistent with the experimental results.Matsui et al. (2008, 2012, 2010)

We also considered non-Markovian transitions using power-law waiting time distribution. In disordered media, power-law waiting time distribution results from the exponential distribution of activation barrier for hopping. The initial decay of the spin correlation function was not influenced by the presence of heavy tail property in the waiting time distribution. The asymptotic decay component influenced by the hopping frequency approaches to Gaussian-like decay by increasing heavy tail property in the waiting time distribution function. Moreover, the on-set time of slowing down of spin relaxation increases as the hopping motion is sluggish in disordered media.

For simplicity, we only considered two-site transitions. As pointed out previously, there may be an extra subtle effect for extended systems and the results would depend on the dimensionality.Czech and Kehr (1984, 1986); Mitra and Le Doussal (1991) When the mobility of a spin carrier is high, such effects may be important. We also assumed the Gaussian random local field to be in the direction of the external field. Although we considered a simplified model, our results captured an additional effect often ignored in assuming Gaussian–Markovian field fluctuations by explicitly considering adiabatic transitions in locally inhomogeneous field environments.

The motional narrowing effect has been observed in the optical resonance absorption line shapes as well as the magnetic resonance absorption line shapes. Sumi (1977); Blumen and Silbey (1978); Jackson and Silbey (1981); Cao and Silbey (2009); Dutta et al. (2017) In optical line shapes, the quantum coherent effects on line shapes give additional difficulties to analyze the motional narrowing effect. The role of quantum coherence in the narrowing of optical line shape under both the site energy fluctuation and the off-diagonal coupling between sites has been treated recently. Dutta et al. (2017) As a future problem, it is necessary to take quantum coherence into account to study the narrowing effect on the optical line shapes.

Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number 15K05406. We would like to thank Professor K. Kitahara for informing us Ref. W. Kehr et al., 1987.

Appendix A. Derivation of Eq. (5)

We define $M(t)=M_{a}(t)+M_{b}(t)$ and $q(t)=M_{a}(t)-M_{b}(t)$. The equations of motion for $M(t)$ and $q(t)$ are obtained as

[TABLE]

where the initial condition is given by $M(0)=M_{0}$ and $q(0)=0$; $\omega_{m}$ and $\delta$ are defined by

[TABLE]

Using the Laplace transformation, defined by $\hat{f}(s)=\int_{0}^{\infty}dt\exp(-st)f(t)$ for a function $f(t)$, and eliminating $\hat{q}(s)$ by taking the Laplace transformation of Eqs. (A.1)–(A.2), we obtain $\hat{\phi}(s)=\hat{\phi}_{r}(s-i\omega_{m}+2\gamma_{s})$, where the Laplace transform of $\phi_{r}(t)$ is given by

[TABLE]

By the inverse Laplace transformation, we obtain

[TABLE]

The term $\langle\cdots\rangle$ in Eq. (3) can be calculated by first averaging with respect to $\omega_{m}$ and subsequently averaging with respect to $\delta$; the averaging procedure can be expressed as $\langle\langle\cdots\rangle_{\omega_{m}}\rangle_{\delta}$, where each average is given by

[TABLE]

The averaging of $\phi(t)$ with respect to $\omega_{m}$ is obtained as

[TABLE]

The averaging with respect to $\delta$ can be calculated after the Laplace transformation. $\langle\hat{\phi}_{r}(s)\rangle_{\delta}$ is obtained as

[TABLE]

for $s>0$ and $\gamma>0$. Introducing the asymptotic expansion,

[TABLE]

and denoting the inverse Laplace transformation by ${\cal L}^{-1}\left[\cdots\right]$, we then obtain

[TABLE]

where ${\cal M}\left(a;b;z\right)$ is the Kummer confluent hypergeometric function. Olver et al. (2010) Finally, we obtain Eq. (5).

Appendix B. Derivation of Eq. (31)

Spin precession and relaxation during the waiting time denoted by $\Delta t$ can be expressed as $\exp\left[(i\omega_{b}-2\gamma_{s})\Delta t\right]$ for the spin at site b. The occupation probability of spins at each site at time $t$ can be obtained by the time convolution of the time distribution of spin just arriving at the site and the remaining probability of the arrived spin. First, we note that the spin state $m_{a}(t)$ for the spin just arriving at site a at time t obeys the relation given by,

[TABLE]

Then, we calculate the remaining probability from the waiting time distribution as $\varphi(t)=1-\int_{0}^{t}d\,t_{1}\psi(t_{1})=\int_{t}^{\infty}d\,t_{1}\psi(t_{1})$; the remaining probability is the fraction of spin escaping from transition to site b until time $t$. The spin state at site a at time t is obtained as $M_{a}(t)=\int_{0}^{t}d\,t_{1}\exp\left[(i\omega_{a}-2\gamma_{s})\left(t-t_{1}\right)\right]\varphi(t-t_{1})m_{a}(t_{1})$, where the spin precession and relaxation during the remaining period are taken into account. We also define the similar properties for the site b.

By applying Laplace transformation to Eq. (A.14) and subtracting $\hat{\psi}\left(s+i\omega_{a}+2\gamma_{s}\right)\hat{m}_{a}(s)$, we obtain,

[TABLE]

Using $\hat{M}_{i}(s)=\left[1-\hat{\psi}\left(s+i\omega_{i}+2\gamma_{s}\right)\right]\hat{m}_{i}(s)/\left(s+i\omega_{i}+2\gamma_{s}\right)$ for $i=a,b$, Eq. (A.14) can be rewritten as

[TABLE]

By introducing Eq. (30), the inverse Laplace transform of Eq. (A.16) is obtained as

[TABLE]

where the generalized hopping frequency is defined by Eq. (A.19). By exchanging a and b we also obtain the similar equation for $M_{b}(t)$. Equation (A.19) is the generalization of Eq. (1) when the waiting time distribution for transitions has heavy tail property. We assume the even site occupation probability of spins in equilibrium. The equilibrium initial condition can be expressed as $M_{a}(0)=M_{b}(0)=M_{0}/2$. As before, we introduce $M(t)=M_{a}(t)+M_{b}(t)$, $q(t)=M_{a}(t)-M_{b}(t)$, $\omega_{m}$ and $\delta$ defined by Eqs. (A.3)-(A.4). We calculate the spin correlation averaged over the distribution of the local Larmor angular frequency and the initial condition.

Using the Laplace transformation, we obtain $\hat{\phi}(s)=\hat{\phi}_{r}(s-i\omega_{m}+2\gamma_{s})$ and

[TABLE]

where the generalized hopping frequency is defined by

[TABLE]

By the inverse Laplace transformation, we obtain

[TABLE]

The averaging of $\phi(t)$ with respect to $\omega_{m}$ is obtained as

[TABLE]

The averaging with respect to $\delta$ can be calculated after the Laplace transformation applying to $\phi_{r}(t)$. $\langle\hat{\phi}_{r}(s)\rangle_{\delta}$ is obtained from Eq. (A.18) as

[TABLE]

for $s>0$ and $\gamma_{r}>0$. By introduce the asymptotic expansion given by Eq. (A.11), we obtain

[TABLE]

We note the useful relation given by Podlubny (1999)

[TABLE]

where $E_{a,b}(x)$ denotes the generalized Mittag–Leffler function and $E_{a,b}^{(n)}(x)$ is its derivatives of order $n$. The generalized Mittag–Leffler function is defined by $E_{a,b}(x)=\sum_{n=0}^{\infty}x^{n}/\Gamma(an+b)$. The inverse Laplace transform of Eq. (A.23) can be performed by using Eq. (A.24). The result is given by Eq. (31).