- Email: [email protected]
On the numerical evaluation of the dielectric relaxation time distribution function from permittivity data
363 Advances in Molecular Relaxation Processes, 5 (1973) 363-373 0 Elsevier Scientific Publishing Company, Amsterdam - Printed
in The Netherlands
ON THE NUMERICAL EVALUATION OF THE DIELECTRIC RELAXATION TIME DISTRIBUTION FUNCTION FROM PERMPTTIVITY DATA
K. GIESE III. Physikalisches
Institut
der Universitiit,
34 GGttingen
(W. Germany)
CONTENTS
I. Introduction ................ II. The determination of g(t) from e(f)
...... of eqn. (5) ......
A. The Fourier inversion B. Computation of the discrete Fourier C. The resolving capacity ........... References ...................
transforms
. . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.
. . . .
363 365 365 367 370 373
1. INTRODUCTION
In the linear response theory a description of the macroscopic properties of materials is provided by some characteristic functions, each completely determines the behaviour of the material. Depending upon of experiment, performed either in the frequency or time domain, the ments are discussed in terms of the complex permittivity or in terms of
dielectric of which the kind measurethe pulse
or step response function. If a time-dependent electric field E(t) acts upon the material, the course in time of the dielectric displacement o(t) is given as a convolution of the pulse response function s(t) and the field E(t) according to D(r) = so&(t) * E(t)
(1)
where s0 is the permittivity of vacuum. Integration of the pulse response function yields the step response function. The Fourier transform of eqn. (1) may be written as D(f)
= W(f)W)
(2)
where the distinction between time functions and spectra is expressed only by the different arguments, time t or frequency f, whereas equal symbols are used to emphasize the connection between corresponding functions. The Fourier transform of the pulse response function s(t) is the complex relative permittivity s(f)
= s’(f) -j&“(f)
= /yIc(r)
exp ( - j2$r)
dr
(3)
364
K.
the real and imaginary
parts of which are the Hilbert
transforms
GIESE
of one another,
so that, in principle, the knowledge of one of the functions I, c’(f) and s”(f) allows the computation of the others. In practice, however, these functions are frequently obtained from experiments within only a finite range of time or frequency at a finite number of measuring points. The evaluation of the integral transformations then gives rise to errors by truncation and aliasing’, for which in many cases limits can be given2. Until recently only frequency domain measurements of the complex permittivity were applicable to the investigation of fast dielectric relaxation processes as given by the orientational motion of small polar molecules. Due to improvements in instrumentation, time domain techniques are available now in a range corresponding to an upper frequency limit3-’ of about 15 GHz. The evaluation of the measuring data obtained by these fast time domain spectroscopy methods, however, requires mathematical operations which have to be performed in the frequency domain. Even here, therefore, the information about the dielectric material is primarily preserved in terms of the complex permittivity. A direct measurement of the pulse or step response function is restricted to relatively slow polarization processes. In the investigation of the relaxation behaviour of dielectric materials the function of greatest physical relevance is the relaxation time distribution function g(z), which has to be determined from the measurable time or frequency functions. Apart from an additional contribution due to conductivity the pulse response function s(t) depends upon g(z) by [E(Z) 5 0 for t < 0] c(t) = c,a(t)+(&~-I:I)SY(exP(-t/T)g(T)d~)/?
(4)
0
and the complex
c(f) =
relative
G, +
(G-4
permittivity
by
sm(g(T)dT)/(l 0
+GWT)
(5)
with
s
omg(r)dr = 1
and g(z) 2 0 in the case of pure relaxation behaviour. The evaluation of g(r) from E(t) or s(f) requires the inversion of the integral transforms eqn. (4) or eqn. (5). The inversion of the former equation may be accomplished by expanding E(t) in a sequence of functions which themselves are Laplace transforms of known functions. The technique presented by Twomey6 yields an expansion of the relaxation time distribution function in terms of Laguerre functions. Frequently the continuous distribution function g(r) is replaced by a number of &functions. If the pulse response function is numerically given by 2N equally
DIELECTRIC
spaced
RELAXATION
data points,
TIME
the application
positions
of N &functions.
presented
by Prony’
DISTRIBUTION
of Prony’s
A modern
365
FUNCTION
method’
formulation
yields the weights
of this method,
and
originally
in 1795, has recently been given by the use of Z-transform and Jen” discussed in detail the accuracy of the time function
theoryg. Strehlow analysis in the case of a superposition of some exponentials, one of which has a considerably smaller amplitude than the others. Similar inversion techniques have been proposed for the determination of the distribution function from the permittivity data by Lassier and Brotl’ and by Linden”. Even here the relaxation time distribution function is approached by a sequence of &functions. As can be seen from the work of Lindon, however, the resolving capacity of the method seems to be unsatisfactory for practical applications in the case of fairly narrow relaxation time distributions. In 1941 Fuoss and Kirkwood’ 3 presented a unique solution for the evaluation of the continuous distribution function from the I” curve. They were especially interested in the problem of computing g(r) from an analytical function which was found empirically to represent the measured loss-number curve. The method of fitting empirical functions to experimental permittivity data, performed by proceeding either from a complex plane representation of the data or from the individual s’(f) and a”(f) curves, is surely the most powerful method if s(f) is known at only a few frequencies. It becomes unsatisfactory if both the measuring range and the measuring accuracy are increased. Frequently the difficulty arises of finding suitable formulae which, on the one hand, describe the experimental complex I data and, on the other hand, have a real and imaginary part which are Hilbert transforms of each other14. An increase in measuring accuracy together with an increase in measuring range and sampling rate, to be obtained by the new time domain spectroscopy methods at even higher frequencies, means, however, that numerical methods become applicable to the integral transformation proposed by Fuoss and Kirkwood. It is the purpose of the present paper to discuss in detail the possibility of a numerical evaluation of the integral transforms obtained by Fourier inversion of both the real and the imaginary parts of eqn. (5). A first application of this method was given by the work of E. Marchal and J. Marchal15, who obtained results showing a significantly smaller resolving capacity than is obtained by other numerical methods. As will be seen below, this is caused mainly by errors due to aliasing and truncation which in most cases can be reduced to a reasonable degree.
II.
THE DETERMINATION
A. The Fourier
Equation
OF g(7)
FROM E(f)
inversion of eqn. (5)
(5) can be written
in a simpler
notation
by the use of logarithmic
366
K.
scales for both the relaxation s = In (r/r,)
time z and the frequencyf.
and
If one defines
r = ln (LX!!)
with 27&r,, = 1, both scales are logarithmic time scales. The frequency be chosen to correspond to the maximum of the s”(f) curve. With h(s) = rg(r)
(6) f0 may
(7)
the real and imaginary &‘(I-)-&,
GIESE
parts of eqn. (5) are given by
= $(l +h(r)
* tanh r)
%--Em
(8)
and d’(r)
___~
= *h(r)
* sech r
Es--Em
respectively. permittivity
A symmetric notation for both the real and imaginary parts of the is obtained by proceeding from the derivative of the dispersion curve
d d(r) -------_ dr
4 h(r)
* sech2 r
E,-E,
Starting from eqn. (9) Fuoss and Kirkwood evaluated h(r) by deconvolving the loss-number function d’(r) from the function sech r. lf we use, in a short-hand notation, the symbols S{f(u)} and s-‘{f(r)} for the Fourier transform and the inverse Fourier transform, respectively, with the corresponding Fourier variables II and r, their result may be written as (Es-8,)
h(r) = (2/z)s{cosh
rc2u S- ‘{d’(r)}>
(11)
or more simply as (Ed-
E,)
h(r) = (2/n) Re {d’(r +jx/2)}
(12)
Here d’(r+jn/2) is the analytic continuation of the function d’(r) to complex values of the argument, the real part of which yields the relaxation time distribution function h(r). Similarly one obtains, proceeding from eqns. (10) and (8)
(13) and (s,--E,)h(r)
= (2/7c) Im {c’(r+jn/2)}
(14)
With d(r) or d’(r) given as analytical functions, as obtained from experimental values by curve fitting, eqn. (14) or (12) yields straightforwardly the relaxation time distribution function h(r). A numerical evaluation has to be based on eqn. (13) or (1 I), possibly by the use of the Fast Fourier Transform algorithm16.
DIELECTRIC
RELAXATION
B. Computation
TIME
DISTRIBUTION
of the discrete Fourier
FUNCTION
367
transforms
For the computation of the discrete Fourier transform to be practicable we will assume that the permittivity functions e’(r) and E”(r) are obtained from measurements as sequences of data points equally spaced by an amount 6r = In formulation (fi+l/fi). The e’(r) d a t a may be processed by the use of Samulon’s of the sampling theorem17. The relation
holds for small values of U. The evaluation of the discrete Fourier transforms SF-‘{de’(r)/dr} and 3-l (E”(r)} has to be restricted to u values corresponding to lu16r 6: 1 in order to keep the errors due to aliasing small. For practical applications the admissible range may be extended up to the limit given by lu16r = 0.1. In general E’(r) and E”(r) are determined only within a finite measuring range Ar = In (f,,Jfmin). Th e introduction of truncated sums yields the second source of error, which for larger u values influences the computation of amount. As the E”(r) curve is ;F- ’ {ds’(r)/dr) and 3-l (E”(r)} by a considerable broader than the de’(r)/dr curve, one might expect that the truncation error is of greater importance if one proceeds from the former curve. With u0 given by the condition that for IuI 5 u0 the computation of the discrete Fourier transforms of the truncated permittivity functions yield spectra nearly unaffected by the aliasing and truncation errors, it is convenient to introduce a window W(U) with w(u) = 0 for ju\ 2 u0 and to write eqns. (11) and (13) as (s,---E,)h(r)
* ;F{w(u)}
= (Z/rc)~F(w(u) cash rr2u 3-‘{&‘(r)})
(1%
Here u,, determines the useful range of definition of the spectra and, as follows from eqns. (16) the resolving capacity of the numerical method for the evaluation of the relaxation time distribution function h(r). In general the value of u0 to be inserted into the calculation depends in a complicated manner upon the measuring range Ar, the relative shift of this range against the permittivity curves, the sample width 6r, and upon the shape of the curves, which means upon h(r) itself. As will be shown below, in simple cases a quantitative estimation of u0 is possible. A usual window w(u) is the rectangular window according to
1’ w1(u) = \O the spectrum
for IuJ _I for IuI >
of which is given
uo uo by
K. GIESE
368
Compared to other window functions, the spectrum of the rectangular window has a rather narrow main lobe. The slow decrease of the side lobes is a disadvantage. To get an idea of the quantity u0 it seems to be useful to consider the case of a simple Debye relaxation process. If its relaxation time is equal to zO, eqn. (6), one has
S-‘(h(s)}
(17a)
3 1 = 1. Es--Em
(17b)
The result of a discrete Fourier transform computation of the right-hand sides of eqns. (17a) and (17b), the former being evaluated by the use of eqn. (15), is shown in Figs. 1 and 2. Outside a range of width Ar = 8, symmetrically situated with respect to the frequency so = 1/27cz, on a logarithmic frequency scale, the functions de’(r)/dr and E”(Y) are assumed to be zero. The width Ar = 8 corresponds to a frequency range off max/f m,n = 3 x 103. Three different sampling rates are used:
6r = 0.8 (fi+Jf,
= 2.25),
6r = 0.4 (fi+i/fi
= 1.49) and the limiting
___,,,r.a., y... ‘.” a
-~
:
Fig. 1. Evaluation of the discrete Fourier transform of the truncated function ds’(r)/dr = (sech2r)/2 (Debye relaxation process). The function is truncated outside the interval jr1 5 4 and the sample width is 6r = 0.8 (a), 6r = 0.4 (b). Curve (c) is calculated for the case of an infinite sampling rate.
DIELECTRIC
RELAXATION
TIME
DISTRIBUTION
FUNCTION
369
Fig. 2. Evaluation of the discrete Fourier transform of the truncated function E”(I) = (sech r)/2 (Debye relaxation process). The function is truncated outside the interval Irl 5 4 and the sample width is 6r = 0.8 (a), Br = 0.4 (b). Curve (c) is calculated for the case of an infinite sampling rate.
rate Sr + 0 (fi+l/fl --t1). The choice of a measuring range situated symmetrically with regard to the permittivity curves yields spectra which are symmetric functions in the variable u so that only the values for u 2 0 are plotted in Figs. 1 and 2. While the exact evaluation of the right-hand side of eqn. (17a) yields the value 1 for any U, the discrete Fourier transform of the truncated function de’(r)/dr shows a strong dependence upon U. The course of curve (a) in Fig. 1 is mainly determined by errors due to aliasing. If the sampling rate is increased by a factor of two or more (curves (b) and (c)), the truncation error becomes dominating, resulting in an exponentially increasing oscillation of period 2/Ar. In the range JUl 5 2.40w 1, however, the deviation from the exact value 1 remains small. Proceeding from the d'(r) curve, eqn. (17b), the result is influenced more strongly by the truncation process. Within the considered measuring range of width Ar = 8 the function d'(r) decreases to only about 2 ‘A of its maximum value so that even for u = 0 a visible deviation from the value 1 is obtained. The oscillation amplitude of curves (a), (b) and (c) in Fig. 2 remains within a 10 % limit for I2415 u0 z 0.35, which is less than half the range as obtained by proceeding from the d(r) curve.
370
K.
GIESE
C. The resolving capacity
To demonstrate the resolving capacity of the method, we have evaluated eqn. (16b) by considering the following relaxation processes: (i) the superposition of two Debye relaxation processes, h(r) = (d(r-r,)+d(r+r,))/2, and (ii) a relaxation process corresponding to a distribution function h(r) = 1/2r, for Irl 5 rz and h(r) = 0 for Irl > r2 as introduced by Frohlich”. With rl = 0.5 and r2 = 0.902 both distribution functions yield the same value for the maximum of the loss number function s”(O) = 0.443, and the d(r) curves, calculated from the two distribution functions, deviate from each other only by a few percent. By the use of a rectangular window wr(u) with u0 = 1 the evaluation of eqn. (16b) yields the curves for the relaxation time distribution function h(r) as plotted in Figs. 3 and 4. The functions de’(r)/dr were assumed to be zero outside a range given by Irl >= 4.6 (fmax/fmin = 104) and the sample width was 6r = 0.1. With the choice of u0 = 1 errors due to aliasing and truncation are sufficiently suppressed. In addition to the calculated curves the figures show the distribution functions on which the computation was based. As can be seen from Fig. 3, the by the two lines of distance two Debye relaxation processes, characterized resolved from each other. Unfortunately the 2ri (rk = 2.72), are distinctly side lobes are large and decrease to zero only slowly. Smaller oscillations are obtained in the case of the continuous relaxation time distribution function.
b
hlr)
\An
0
2
5
-r Fig. 3. Approximation of the relaxation time distribution function as given by the superposition of two Debye relaxation processes. The curve is calculated by the use of a rectangular window wI(u) with u,, = 1.
DIELECTRIC
-5
RELAXATION
-2.5
TIME
0
~lsTRtBuTl0N
371
FUNCTION
5
2.5
-,
Fig. 4. Approximation of the Friihlich relaxation time distribution culated by the use of a rectangular window w1 (u) with uO = 1.
function.
The curve
is cal-
Figure 4 shows the rectangular distribution function as approached by the calculated function. Though the resolving capacity obtained by the use of the rectangular window is rather good, the large side lobes characterizing the spectrum of this window are disadvantageous. By the use of other windows the side lobes can be reduced. However, parallel to this desirable effect, the width of the main lobe increases. The spectrum of the “hanning” window” (l+cos(7r~/u~))/2
w2(u)=
() i
for 1~1 5 uO for (uJ > u.
shows side lobes of which the heights remain
within
a 3 “/, limit as compared
to
the height of the main lobe. Its main lobe has twice the width of the main lobe obtained for the rectangular window. Figure 5 shows the result of an evaluation of eqn. (16b) by the use of the hanning window We. The computation was based on the same dispersion curve as considered in the case of Fig. 3. In view of the smaller resolving capacity to be obtained by the use of the hanning window both the measuring range Ar and the sampling rate l/& were increased by a factor of two. Proceeding from the loss number function s”(r), the application of eqn. (16a) together with a corresponding smaller value of u,, reduces the resolving capacity by a considerable amount. The loss number function therefore seems to be of less value for determining the relaxation time distribution function, as the resolving capacity obtained above by proceeding from the dispersion curve would require the loss number curve to be measured within an unrealistically large frequency range. Now in many cases of dielectric relaxation the permittivity curves show an asymp-
372
K. GIESE
1
h(r) 0.5 T
-
1
-2.5
-5
0
5
25
-r
Fig. 5. Approximation of the relaxation time distribution of two Debye relaxation processes. The curve is calculated with u,, = 2.
function as given by the superposition by the use of a hanning window wz(u)
totic behaviour which allows an extrapolation of the curves. The frequently used empirical formulae of Cole and Cole 2o, Davidson and Cole” and Havriliak and Negami 22, for example, correspond to a frequency dependence of the loss number function which in a double logarithmic plot allows a linear extrapolation if the virtual measuring range has a sufficient extension. The possibility of an extrapolation of the measured curves, however, means that the error due to truncation can be reduced to any small amount. A quantitative description may be given for relaxation processes corresponding to relaxation time distribution functions h(s) which are zero outside a range given by sr 5 s S s2. Then the low-frequency asymptotic behaviour is described by e(r) -00
E
r
l-e-2r
J SI
G--E*
and the high-frequency
E(r)--Em _ e2r %--Em
r m
s2
=
e2”h(s)ds-je-’
asymptotic
r s*
JSl
ee2”h(s)ds-je’
J s,
behaviour
r
e”h(s) ds
(18)
by
s2
JS,
e-“h(s)ds
(19)
If the s”(r) curve is assumed to be zero outside an interval Ar = r2 -rl, large enough for the asymptotic formulae as given by eqns. (18) and (19) to be applicable, the evaluation of the inverse Fourier transform causes an error as given by
DIELECTRIC
RELAXATION
TIME
DISTRIBUTION
FUNCTION
313
with E(u) given by er~(1 + j2nu) E(u) = 1 +j2rcu
s2 s
e-“h(s)ds s1
e-rz(l
-j2nu)
s*
e”h(s) ds
+ l-j2rru
If the width of the relaxation time distribution is integrals in eqn. (21) are not very different from can be used to estimate the dependence of the Similar relations hold if the calculation is based on
s s,
(21)
not too large, the values of the unity. Equations (20) and (21) truncation error on rl and r2. the truncated function de’(r)/dr.
REFERENCES
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
A. M. NICOLSON, IEEE Trans. Instrum. Meas., IM-17 (1968) 395. F. R. SCHWARZL AND L. C. E. STRUIK, Advan. Mol. Relaxation Processes, 1 (1968) 201. H. FELLNER-FELDEGG,J. Phys. Chem., 73 (1969) 616. H. FELLNER-FELDEGG,J. Phys. Chem., 76 (1972) 2116. A. SUGGETT, in M. DAVIES (Ed.), Dielectric and Related Molecular Processes, Vol. 1, The Chemical Society, London, 1972. S. TWOMEY, J. Franklin Inst., 275 (1963) 121. R. A. BUCKINGHAM, Numerical Methods, Sir Isaac Pitman and Sons, London, 1962. R. PRONY, J. Ecole Polytech., 1 (1795) 24. C. S. BURRUS AND T. W. PARKS, IEEE Trans. Audio Electra Acoustics, AU-18 (1970) 137. H. STREHLOWAND J. JEN, On the accuracy of chemical relaxation measurements, Report of the Max-Planck-Institut fiir Physikalische Chemie, Gattingen, 1970. B. LASSIER AND C. BROT, in Magnetic resonance and relaxation, Proc. Colfoq. AMPERE, 14 (1967) 693. P. LINDON, Proc. IEEE, 58 (1970) 1389. R. M. Fuoss AND J. G. KIRKWOOD, J. Amer. Chem. Sot., 63 (1941) 385. K. GIESE, U. KAATZE AND R. POTTEL,J. Phys. Chem., 74 (1970) 3718. E. MARCHAL, AND J. MARCHAL, Arch. Sci. (Geneoa), 13 (1960) 82. J. W. COOLEY AND J. W. TUKEY, Math. Compm., 19 (1965) 297. H. A. SAMULON, Proc. IRE, 39 (1951) 175. H. FRBHLICH, Theory of Dielectrics, Clarendon Press, Oxford, 1958. R. B. BLACKMAN AND J. W. TUKEY, The Measurement of Powder Spectra, Dover Publications, New York, 1959. K. S. COLE AND R. H. COLE, J. Chem. Phys., 10 (1942) 98. D. W. DAVIDSON AND R. H. COLE, J. Chem. Phys., 19 (1951) 1484. S. HAVRILIAK AND S. NEGAMI, J. Polym. Sci., Part C, 14 (1966) 99.
in The Netherlands
ON THE NUMERICAL EVALUATION OF THE DIELECTRIC RELAXATION TIME DISTRIBUTION FUNCTION FROM PERMPTTIVITY DATA
K. GIESE III. Physikalisches
Institut
der Universitiit,
34 GGttingen
(W. Germany)
CONTENTS
I. Introduction ................ II. The determination of g(t) from e(f)
...... of eqn. (5) ......
A. The Fourier inversion B. Computation of the discrete Fourier C. The resolving capacity ........... References ...................
transforms
. . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.
. . . .
363 365 365 367 370 373
1. INTRODUCTION
In the linear response theory a description of the macroscopic properties of materials is provided by some characteristic functions, each completely determines the behaviour of the material. Depending upon of experiment, performed either in the frequency or time domain, the ments are discussed in terms of the complex permittivity or in terms of
dielectric of which the kind measurethe pulse
or step response function. If a time-dependent electric field E(t) acts upon the material, the course in time of the dielectric displacement o(t) is given as a convolution of the pulse response function s(t) and the field E(t) according to D(r) = so&(t) * E(t)
(1)
where s0 is the permittivity of vacuum. Integration of the pulse response function yields the step response function. The Fourier transform of eqn. (1) may be written as D(f)
= W(f)W)
(2)
where the distinction between time functions and spectra is expressed only by the different arguments, time t or frequency f, whereas equal symbols are used to emphasize the connection between corresponding functions. The Fourier transform of the pulse response function s(t) is the complex relative permittivity s(f)
= s’(f) -j&“(f)
= /yIc(r)
exp ( - j2$r)
dr
(3)
364
K.
the real and imaginary
parts of which are the Hilbert
transforms
GIESE
of one another,
so that, in principle, the knowledge of one of the functions I, c’(f) and s”(f) allows the computation of the others. In practice, however, these functions are frequently obtained from experiments within only a finite range of time or frequency at a finite number of measuring points. The evaluation of the integral transformations then gives rise to errors by truncation and aliasing’, for which in many cases limits can be given2. Until recently only frequency domain measurements of the complex permittivity were applicable to the investigation of fast dielectric relaxation processes as given by the orientational motion of small polar molecules. Due to improvements in instrumentation, time domain techniques are available now in a range corresponding to an upper frequency limit3-’ of about 15 GHz. The evaluation of the measuring data obtained by these fast time domain spectroscopy methods, however, requires mathematical operations which have to be performed in the frequency domain. Even here, therefore, the information about the dielectric material is primarily preserved in terms of the complex permittivity. A direct measurement of the pulse or step response function is restricted to relatively slow polarization processes. In the investigation of the relaxation behaviour of dielectric materials the function of greatest physical relevance is the relaxation time distribution function g(z), which has to be determined from the measurable time or frequency functions. Apart from an additional contribution due to conductivity the pulse response function s(t) depends upon g(z) by [E(Z) 5 0 for t < 0] c(t) = c,a(t)+(&~-I:I)SY(exP(-t/T)g(T)d~)/?
(4)
0
and the complex
c(f) =
relative
G, +
(G-4
permittivity
by
sm(g(T)dT)/(l 0
+GWT)
(5)
with
s
omg(r)dr = 1
and g(z) 2 0 in the case of pure relaxation behaviour. The evaluation of g(r) from E(t) or s(f) requires the inversion of the integral transforms eqn. (4) or eqn. (5). The inversion of the former equation may be accomplished by expanding E(t) in a sequence of functions which themselves are Laplace transforms of known functions. The technique presented by Twomey6 yields an expansion of the relaxation time distribution function in terms of Laguerre functions. Frequently the continuous distribution function g(r) is replaced by a number of &functions. If the pulse response function is numerically given by 2N equally
DIELECTRIC
spaced
RELAXATION
data points,
TIME
the application
positions
of N &functions.
presented
by Prony’
DISTRIBUTION
of Prony’s
A modern
365
FUNCTION
method’
formulation
yields the weights
of this method,
and
originally
in 1795, has recently been given by the use of Z-transform and Jen” discussed in detail the accuracy of the time function
theoryg. Strehlow analysis in the case of a superposition of some exponentials, one of which has a considerably smaller amplitude than the others. Similar inversion techniques have been proposed for the determination of the distribution function from the permittivity data by Lassier and Brotl’ and by Linden”. Even here the relaxation time distribution function is approached by a sequence of &functions. As can be seen from the work of Lindon, however, the resolving capacity of the method seems to be unsatisfactory for practical applications in the case of fairly narrow relaxation time distributions. In 1941 Fuoss and Kirkwood’ 3 presented a unique solution for the evaluation of the continuous distribution function from the I” curve. They were especially interested in the problem of computing g(r) from an analytical function which was found empirically to represent the measured loss-number curve. The method of fitting empirical functions to experimental permittivity data, performed by proceeding either from a complex plane representation of the data or from the individual s’(f) and a”(f) curves, is surely the most powerful method if s(f) is known at only a few frequencies. It becomes unsatisfactory if both the measuring range and the measuring accuracy are increased. Frequently the difficulty arises of finding suitable formulae which, on the one hand, describe the experimental complex I data and, on the other hand, have a real and imaginary part which are Hilbert transforms of each other14. An increase in measuring accuracy together with an increase in measuring range and sampling rate, to be obtained by the new time domain spectroscopy methods at even higher frequencies, means, however, that numerical methods become applicable to the integral transformation proposed by Fuoss and Kirkwood. It is the purpose of the present paper to discuss in detail the possibility of a numerical evaluation of the integral transforms obtained by Fourier inversion of both the real and the imaginary parts of eqn. (5). A first application of this method was given by the work of E. Marchal and J. Marchal15, who obtained results showing a significantly smaller resolving capacity than is obtained by other numerical methods. As will be seen below, this is caused mainly by errors due to aliasing and truncation which in most cases can be reduced to a reasonable degree.
II.
THE DETERMINATION
A. The Fourier
Equation
OF g(7)
FROM E(f)
inversion of eqn. (5)
(5) can be written
in a simpler
notation
by the use of logarithmic
366
K.
scales for both the relaxation s = In (r/r,)
time z and the frequencyf.
and
If one defines
r = ln (LX!!)
with 27&r,, = 1, both scales are logarithmic time scales. The frequency be chosen to correspond to the maximum of the s”(f) curve. With h(s) = rg(r)
(6) f0 may
(7)
the real and imaginary &‘(I-)-&,
GIESE
parts of eqn. (5) are given by
= $(l +h(r)
* tanh r)
%--Em
(8)
and d’(r)
___~
= *h(r)
* sech r
Es--Em
respectively. permittivity
A symmetric notation for both the real and imaginary parts of the is obtained by proceeding from the derivative of the dispersion curve
d d(r) -------_ dr
4 h(r)
* sech2 r
E,-E,
Starting from eqn. (9) Fuoss and Kirkwood evaluated h(r) by deconvolving the loss-number function d’(r) from the function sech r. lf we use, in a short-hand notation, the symbols S{f(u)} and s-‘{f(r)} for the Fourier transform and the inverse Fourier transform, respectively, with the corresponding Fourier variables II and r, their result may be written as (Es-8,)
h(r) = (2/z)s{cosh
rc2u S- ‘{d’(r)}>
(11)
or more simply as (Ed-
E,)
h(r) = (2/n) Re {d’(r +jx/2)}
(12)
Here d’(r+jn/2) is the analytic continuation of the function d’(r) to complex values of the argument, the real part of which yields the relaxation time distribution function h(r). Similarly one obtains, proceeding from eqns. (10) and (8)
(13) and (s,--E,)h(r)
= (2/7c) Im {c’(r+jn/2)}
(14)
With d(r) or d’(r) given as analytical functions, as obtained from experimental values by curve fitting, eqn. (14) or (12) yields straightforwardly the relaxation time distribution function h(r). A numerical evaluation has to be based on eqn. (13) or (1 I), possibly by the use of the Fast Fourier Transform algorithm16.
DIELECTRIC
RELAXATION
B. Computation
TIME
DISTRIBUTION
of the discrete Fourier
FUNCTION
367
transforms
For the computation of the discrete Fourier transform to be practicable we will assume that the permittivity functions e’(r) and E”(r) are obtained from measurements as sequences of data points equally spaced by an amount 6r = In formulation (fi+l/fi). The e’(r) d a t a may be processed by the use of Samulon’s of the sampling theorem17. The relation
holds for small values of U. The evaluation of the discrete Fourier transforms SF-‘{de’(r)/dr} and 3-l (E”(r)} has to be restricted to u values corresponding to lu16r 6: 1 in order to keep the errors due to aliasing small. For practical applications the admissible range may be extended up to the limit given by lu16r = 0.1. In general E’(r) and E”(r) are determined only within a finite measuring range Ar = In (f,,Jfmin). Th e introduction of truncated sums yields the second source of error, which for larger u values influences the computation of amount. As the E”(r) curve is ;F- ’ {ds’(r)/dr) and 3-l (E”(r)} by a considerable broader than the de’(r)/dr curve, one might expect that the truncation error is of greater importance if one proceeds from the former curve. With u0 given by the condition that for IuI 5 u0 the computation of the discrete Fourier transforms of the truncated permittivity functions yield spectra nearly unaffected by the aliasing and truncation errors, it is convenient to introduce a window W(U) with w(u) = 0 for ju\ 2 u0 and to write eqns. (11) and (13) as (s,---E,)h(r)
* ;F{w(u)}
= (Z/rc)~F(w(u) cash rr2u 3-‘{&‘(r)})
(1%
Here u,, determines the useful range of definition of the spectra and, as follows from eqns. (16) the resolving capacity of the numerical method for the evaluation of the relaxation time distribution function h(r). In general the value of u0 to be inserted into the calculation depends in a complicated manner upon the measuring range Ar, the relative shift of this range against the permittivity curves, the sample width 6r, and upon the shape of the curves, which means upon h(r) itself. As will be shown below, in simple cases a quantitative estimation of u0 is possible. A usual window w(u) is the rectangular window according to
1’ w1(u) = \O the spectrum
for IuJ _I for IuI >
of which is given
uo uo by
K. GIESE
368
Compared to other window functions, the spectrum of the rectangular window has a rather narrow main lobe. The slow decrease of the side lobes is a disadvantage. To get an idea of the quantity u0 it seems to be useful to consider the case of a simple Debye relaxation process. If its relaxation time is equal to zO, eqn. (6), one has
S-‘(h(s)}
(17a)
3 1 = 1. Es--Em
(17b)
The result of a discrete Fourier transform computation of the right-hand sides of eqns. (17a) and (17b), the former being evaluated by the use of eqn. (15), is shown in Figs. 1 and 2. Outside a range of width Ar = 8, symmetrically situated with respect to the frequency so = 1/27cz, on a logarithmic frequency scale, the functions de’(r)/dr and E”(Y) are assumed to be zero. The width Ar = 8 corresponds to a frequency range off max/f m,n = 3 x 103. Three different sampling rates are used:
6r = 0.8 (fi+Jf,
= 2.25),
6r = 0.4 (fi+i/fi
= 1.49) and the limiting
___,,,r.a., y... ‘.” a
-~
:
Fig. 1. Evaluation of the discrete Fourier transform of the truncated function ds’(r)/dr = (sech2r)/2 (Debye relaxation process). The function is truncated outside the interval jr1 5 4 and the sample width is 6r = 0.8 (a), 6r = 0.4 (b). Curve (c) is calculated for the case of an infinite sampling rate.
DIELECTRIC
RELAXATION
TIME
DISTRIBUTION
FUNCTION
369
Fig. 2. Evaluation of the discrete Fourier transform of the truncated function E”(I) = (sech r)/2 (Debye relaxation process). The function is truncated outside the interval Irl 5 4 and the sample width is 6r = 0.8 (a), Br = 0.4 (b). Curve (c) is calculated for the case of an infinite sampling rate.
rate Sr + 0 (fi+l/fl --t1). The choice of a measuring range situated symmetrically with regard to the permittivity curves yields spectra which are symmetric functions in the variable u so that only the values for u 2 0 are plotted in Figs. 1 and 2. While the exact evaluation of the right-hand side of eqn. (17a) yields the value 1 for any U, the discrete Fourier transform of the truncated function de’(r)/dr shows a strong dependence upon U. The course of curve (a) in Fig. 1 is mainly determined by errors due to aliasing. If the sampling rate is increased by a factor of two or more (curves (b) and (c)), the truncation error becomes dominating, resulting in an exponentially increasing oscillation of period 2/Ar. In the range JUl 5 2.40w 1, however, the deviation from the exact value 1 remains small. Proceeding from the d'(r) curve, eqn. (17b), the result is influenced more strongly by the truncation process. Within the considered measuring range of width Ar = 8 the function d'(r) decreases to only about 2 ‘A of its maximum value so that even for u = 0 a visible deviation from the value 1 is obtained. The oscillation amplitude of curves (a), (b) and (c) in Fig. 2 remains within a 10 % limit for I2415 u0 z 0.35, which is less than half the range as obtained by proceeding from the d(r) curve.
370
K.
GIESE
C. The resolving capacity
To demonstrate the resolving capacity of the method, we have evaluated eqn. (16b) by considering the following relaxation processes: (i) the superposition of two Debye relaxation processes, h(r) = (d(r-r,)+d(r+r,))/2, and (ii) a relaxation process corresponding to a distribution function h(r) = 1/2r, for Irl 5 rz and h(r) = 0 for Irl > r2 as introduced by Frohlich”. With rl = 0.5 and r2 = 0.902 both distribution functions yield the same value for the maximum of the loss number function s”(O) = 0.443, and the d(r) curves, calculated from the two distribution functions, deviate from each other only by a few percent. By the use of a rectangular window wr(u) with u0 = 1 the evaluation of eqn. (16b) yields the curves for the relaxation time distribution function h(r) as plotted in Figs. 3 and 4. The functions de’(r)/dr were assumed to be zero outside a range given by Irl >= 4.6 (fmax/fmin = 104) and the sample width was 6r = 0.1. With the choice of u0 = 1 errors due to aliasing and truncation are sufficiently suppressed. In addition to the calculated curves the figures show the distribution functions on which the computation was based. As can be seen from Fig. 3, the by the two lines of distance two Debye relaxation processes, characterized resolved from each other. Unfortunately the 2ri (rk = 2.72), are distinctly side lobes are large and decrease to zero only slowly. Smaller oscillations are obtained in the case of the continuous relaxation time distribution function.
b
hlr)
\An
0
2
5
-r Fig. 3. Approximation of the relaxation time distribution function as given by the superposition of two Debye relaxation processes. The curve is calculated by the use of a rectangular window wI(u) with u,, = 1.
DIELECTRIC
-5
RELAXATION
-2.5
TIME
0
~lsTRtBuTl0N
371
FUNCTION
5
2.5
-,
Fig. 4. Approximation of the Friihlich relaxation time distribution culated by the use of a rectangular window w1 (u) with uO = 1.
function.
The curve
is cal-
Figure 4 shows the rectangular distribution function as approached by the calculated function. Though the resolving capacity obtained by the use of the rectangular window is rather good, the large side lobes characterizing the spectrum of this window are disadvantageous. By the use of other windows the side lobes can be reduced. However, parallel to this desirable effect, the width of the main lobe increases. The spectrum of the “hanning” window” (l+cos(7r~/u~))/2
w2(u)=
() i
for 1~1 5 uO for (uJ > u.
shows side lobes of which the heights remain
within
a 3 “/, limit as compared
to
the height of the main lobe. Its main lobe has twice the width of the main lobe obtained for the rectangular window. Figure 5 shows the result of an evaluation of eqn. (16b) by the use of the hanning window We. The computation was based on the same dispersion curve as considered in the case of Fig. 3. In view of the smaller resolving capacity to be obtained by the use of the hanning window both the measuring range Ar and the sampling rate l/& were increased by a factor of two. Proceeding from the loss number function s”(r), the application of eqn. (16a) together with a corresponding smaller value of u,, reduces the resolving capacity by a considerable amount. The loss number function therefore seems to be of less value for determining the relaxation time distribution function, as the resolving capacity obtained above by proceeding from the dispersion curve would require the loss number curve to be measured within an unrealistically large frequency range. Now in many cases of dielectric relaxation the permittivity curves show an asymp-
372
K. GIESE
1
h(r) 0.5 T
-
1
-2.5
-5
0
5
25
-r
Fig. 5. Approximation of the relaxation time distribution of two Debye relaxation processes. The curve is calculated with u,, = 2.
function as given by the superposition by the use of a hanning window wz(u)
totic behaviour which allows an extrapolation of the curves. The frequently used empirical formulae of Cole and Cole 2o, Davidson and Cole” and Havriliak and Negami 22, for example, correspond to a frequency dependence of the loss number function which in a double logarithmic plot allows a linear extrapolation if the virtual measuring range has a sufficient extension. The possibility of an extrapolation of the measured curves, however, means that the error due to truncation can be reduced to any small amount. A quantitative description may be given for relaxation processes corresponding to relaxation time distribution functions h(s) which are zero outside a range given by sr 5 s S s2. Then the low-frequency asymptotic behaviour is described by e(r) -00
E
r
l-e-2r
J SI
G--E*
and the high-frequency
E(r)--Em _ e2r %--Em
r m
s2
=
e2”h(s)ds-je-’
asymptotic
r s*
JSl
ee2”h(s)ds-je’
J s,
behaviour
r
e”h(s) ds
(18)
by
s2
JS,
e-“h(s)ds
(19)
If the s”(r) curve is assumed to be zero outside an interval Ar = r2 -rl, large enough for the asymptotic formulae as given by eqns. (18) and (19) to be applicable, the evaluation of the inverse Fourier transform causes an error as given by
DIELECTRIC
RELAXATION
TIME
DISTRIBUTION
FUNCTION
313
with E(u) given by er~(1 + j2nu) E(u) = 1 +j2rcu
s2 s
e-“h(s)ds s1
e-rz(l
-j2nu)
s*
e”h(s) ds
+ l-j2rru
If the width of the relaxation time distribution is integrals in eqn. (21) are not very different from can be used to estimate the dependence of the Similar relations hold if the calculation is based on
s s,
(21)
not too large, the values of the unity. Equations (20) and (21) truncation error on rl and r2. the truncated function de’(r)/dr.
REFERENCES
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
A. M. NICOLSON, IEEE Trans. Instrum. Meas., IM-17 (1968) 395. F. R. SCHWARZL AND L. C. E. STRUIK, Advan. Mol. Relaxation Processes, 1 (1968) 201. H. FELLNER-FELDEGG,J. Phys. Chem., 73 (1969) 616. H. FELLNER-FELDEGG,J. Phys. Chem., 76 (1972) 2116. A. SUGGETT, in M. DAVIES (Ed.), Dielectric and Related Molecular Processes, Vol. 1, The Chemical Society, London, 1972. S. TWOMEY, J. Franklin Inst., 275 (1963) 121. R. A. BUCKINGHAM, Numerical Methods, Sir Isaac Pitman and Sons, London, 1962. R. PRONY, J. Ecole Polytech., 1 (1795) 24. C. S. BURRUS AND T. W. PARKS, IEEE Trans. Audio Electra Acoustics, AU-18 (1970) 137. H. STREHLOWAND J. JEN, On the accuracy of chemical relaxation measurements, Report of the Max-Planck-Institut fiir Physikalische Chemie, Gattingen, 1970. B. LASSIER AND C. BROT, in Magnetic resonance and relaxation, Proc. Colfoq. AMPERE, 14 (1967) 693. P. LINDON, Proc. IEEE, 58 (1970) 1389. R. M. Fuoss AND J. G. KIRKWOOD, J. Amer. Chem. Sot., 63 (1941) 385. K. GIESE, U. KAATZE AND R. POTTEL,J. Phys. Chem., 74 (1970) 3718. E. MARCHAL, AND J. MARCHAL, Arch. Sci. (Geneoa), 13 (1960) 82. J. W. COOLEY AND J. W. TUKEY, Math. Compm., 19 (1965) 297. H. A. SAMULON, Proc. IRE, 39 (1951) 175. H. FRBHLICH, Theory of Dielectrics, Clarendon Press, Oxford, 1958. R. B. BLACKMAN AND J. W. TUKEY, The Measurement of Powder Spectra, Dover Publications, New York, 1959. K. S. COLE AND R. H. COLE, J. Chem. Phys., 10 (1942) 98. D. W. DAVIDSON AND R. H. COLE, J. Chem. Phys., 19 (1951) 1484. S. HAVRILIAK AND S. NEGAMI, J. Polym. Sci., Part C, 14 (1966) 99.