• View
• Edit
• History
• Print

KTDyn

Index

The derivation of the dynamical relaxation by Kubo

We first recall how to calculate zero field relaxation in the strong collision model, according to Kubo and Hayano. Then we show how to implement it in matlab.

The essence of this phenomenological model is that the spin dynamics is modelled by random jumps: at each jump the muon picks up a new random field value from the characteristic distribution (for the sake of the argument we consider here the Gaussian case). In the limit of infinite mean time {$\tau$} between jumps the relaxation function is just the static Kubo-Toyabe distribution {$g(t)$}. The dynamic relaxation may be written as a sum over events with {$0,1,2,\cdots,n,\cdots$} jumps, each described by a characteristic relaxation {$g_n(t)$}

{$(1)\qquad \qquad G(t)=\sum_{n=0}^\infty g_n(t)$}

where the function with no jumps is just again the static Kubo-Toyabe relaxation times the Poissonian probability for no events, {$g(t)e^{-t/\tau}$}. Likewise the relaxation function with one jump is the average over all possible times {$x<t$} for that single jump, which may be written as the evolution with no jumps until time {$x$}, followed again by the evolution with no jumps from time {$x$} on:

{$\qquad\qquad g_1(t)=\int_0^t \frac {dx} {\tau} e^{-(t-x)/\tau}g(t-x)e^{-x/\tau}g(x) $}

where, after the jump, the relaxation restarts from the residual polarization, {$g(x)$}, with the same functional form {$g(t-x)$}. Supposing {$\tau\ll t$} we can take the upper limit of integration to infinity.

{$\qquad\qquad g_1(t)=\int_0^\infty \frac {dx} {\tau} e^{-(t-x)/\tau}g(t-x)e^{-x/\tau}g(x) $}

Taking the product of the two exponentials {$g_1(t)$} reduces to

{$ \qquad\qquad g_1(t)= \frac 1 \tau e^{-t/\tau} \int_0^\infty dx g(t-x)g(x)$}

where we recognize the convolution of two {$g$}. The Laplace transform {$ f_1(s)=\int_0^\infty dt g_1(t) e^{-st}$} may be written therefore, thanks to the convolution theorem, as

{$ \qquad\qquad f_1(s)= \frac 1 \tau \int_0^\infty dt \left[ \int_0^\infty dx g(t-x) g(x) \right] e^{-t(s+\frac 1 \tau)} \,=\,\tau \frac {(f(s+\frac 1 \tau))^2}{\tau^2} $}

in terms of the Laplace transform of {$g(t)$}, {$f(s)=\int_0^\infty g(t) e^{-st} dt$}. Working ot the general terms one recognizes that it may be written in a similar way

{$ \qquad\qquad f_n(s) = \tau \frac {(f(s+\frac 1 \tau))^n}{\tau^n} $}

so that the whole series 1 may be rewritten as a geometrical series of ratio {$\nu f(s+\nu)$}, where {$\nu=\frac 1 \tau$}, starting from {$n=1$}.

Hence it is easy to obtain the Laplace transform {$F(s)$} of the dynamical relaxation function {$G$}

{$ \qquad \qquad F(s) = \tau\sum_{n=1}^\infty (\nu f(s+\nu))^n = \frac {f(s+\nu)}{1-\nu f(s+\nu)}$}

top

How to implement the scheme numerically

Finally this algorithm may be directly implemented nuumerically in the following way. Let us consider a complex frequency {$s+i\omega$}. The Laplace transform {$f(s,\omega)$} of the function {$g(t)$} is then:

{$f(s,\omega)=\int_0^\infty g(t) e^{-st}\,e^{-i\omega t} dt = \int_{-\infty}^\infty g(t) \theta(t) e^{-st}\,e^{-i\omega t} dt $}

where {$\theta(t)$} is the step function. Since {$g(t)$} behaves as a constant at long times this integral converges for any value of the real frequency {$s>0$}, i.e. for any physical value of interest of the real frequency. It is also easy to recognize the integral as the Fourier transform of {$g(t)\theta(t) e^{-st}$} at frequency {$\omega$}, which may therefore be obtained numerically by the FFT algorithm. The quantity {$f(s+\nu)$} is likewise calculated as the FFT of {$g(t)\theta(t) e^{-(s+\nu)t}$} and one can directly obtain

{$F(s)=\frac{f(s+\nu)}{1-\nu f(s+\nu)} \qquad\qquad (1) $}

which is regular as long as {$f(s+\nu)\ne \nu^{-1}$}. The dynamic KT function {$G(t)$}, multiplied by {$e^{-st}$}, may be then obtained by the inverse fast Fourier algorithm as the IFFT of the array {$F$}, or else directly, performing the calculation for {$s=0$}. This reflects the fact that for any Laplace transform, {$f(s,\omega)$} the inverse Laplace transform is given by

{$ g(t)e^{-st}=\int_{-\infty}^{\infty} f(s,\omega)e^{i\omega t} d\omega$}

where {$s$} must be chosen large enough to let the integral converge (for {$G(t)$} this amounts to the same criterion as for the Laplace transform, i.e. {$s+\nu>0$}, provided {$f(s+\nu)\ne \nu^{-1}$}).

In performing the finite sums, however, as it is implied by the fast Fourier algorithm, G. Allodi noted that one must carefully consider what replaces each quantity in Eq. 1. In first approximation, when {$F, f$} replace the FFT sums, the term {$\nu$} must be multiplied by {$dt$} (originating from the replacement of the integral, {$f(s,\nu)$}, by the finite FFT sums). However this term must correctly represents the probability of one scattering event taking place in time {$dt$} while the latter, originally an infinitesimal quantity, becomes finite. Therefore it must be replaced by the cumulative probability density, {$1-e^{-\nu dt}$}, to yield

{$ F_k = \frac {f_k(\nu)}{1-(1-e^{-\nu dt}) f_k(\nu)}$}

Actually this procedure produces the correct behavior at large scattering frequency {$\nu$} (checked againts the approximate Abragam formula, {$G(t)=\exp(-2\frac{\Delta^2}{\nu^2}(e^{-\nu t} -1 +\nu t))$}), when rescaled to have unit initial polarization. For the algorithm itself to produce a function starting at t=0 from unity, it requires larger number of points {$N$} and finer resolution {$dt$} as {$\nu$} increases above 1.

Index