Langmuir 1996, 12, 6691-6700
6691
Experimental Study of Noise on Photon Autocorrelation Functions Daniel Harrison† and Michael R. Fisch*,‡ Departments of Physics, John Carroll University, University Heights, Ohio 44118, and Case Western Reserve University, Cleveland, Ohio 44106 Received June 25, 1996X A theoretical model (Scha¨tzel, K. Quantum Opt. 1990, 2, 287) for the covariance of real quasielastic light-scattering correlation functions has been experimentally verified. The variance calculated from this model was found to be comparable to the variance obtained from a set of experimental data runs while requiring only one or a few runs. The covariance was also verified by direct comparison to the averaged covariance calculated directly from the set of data. The effects of pinhole size upon the covariance were also studied. The nondiagonal covariance was approximately independent of pinhole diameter, while the variance was inversely proportional to the pinhole diameter. Thus, relative noise reduction was verified, and criteria for a suitable compromise between the pinhole diameter and signal-to-background noise ratio was established. This criteria suggests that the optimal pinhole diameter is such that six-ten coherence areas of light are collected. Multiple scattering solutions, with transmissions g40%, were found to have a negligible effect upon this model. The effects of correlation on the error in the fitted parameters were investigated and found to be significant in some cases. Finally, a Mathematica implementation of the covariance into least-squares fitting to account for correlated noise was implemented. The effects of correlated noise on both the fit and the errors in the fit are analyzed. The results are compared to several standard fitting routines.
I. Introduction Photon correlation spectroscopy, also known as dynamic light scattering, provides a very powerful technique for determining the time correlation functions of processes occurring in solutions which scatter light.1 There has been continuing interest in extracting physically interesting information, such as diffusion coefficient distributions, radii or size distributions, and the like, from these correlation coefficients. This process is greatly complicated because the inversion of the experimental correlation function to obtain the distribution of scatterers is mathematically ill-conditioned.2 Many different techniques and algorithms have been developed to perform such analysis. Most of these algorithms either estimate the noise as arising solely from the photon noise or they ignore noise altogether. This is almost always incorrect, as in most experiments which use modern correlators, there are two types of noise: intensity and photon. Furthermore, due to the discrete nature of modern autocorrelation techniques, intensity noise, which is strongly correlated, may dominate the photon noise.3,4 The effects of such correlated noise on data reduction may be rather severe and may result in parameters being unknowingly over-fitted. Lastly, one of the most common assumptions made in data analysis algorithms of all sorts is that the noise is random; intensity noise violates this assumption. A theory which incorporates both types of noise on such autocorrelation functions has been developed by Scha¨tzel.5 †
Case Western Reserve University. ‡ John Carroll University. X Abstract published in Advance ACS Abstracts, December 15, 1996. (1) See, for example: Schmitz, K. S. Dynamic Light Scattering by Macromolecules; Academic: New York, 1990; Pecora, R., Ed. Dynamic Light Scattering; Plenum Press: New York, 1985. (2) Pike, E. R.; et al. In Measurement of Suspended Particles by QuasiElastic Light Scattering; Dahneke, . E., Ed.;Wiley: New York, 1983; pp 107-128. Also see: Provencher, S. W. Comput. Phys. Commun. 1982, 27, 213. (3) Peters, R. In Dynamic Light Scattering; Brown,W., Ed.; Oxford University Press: Oxford, 1993; Chapter 3, pp 149-176. (4) Scha¨tzel, K. In Dynamic Light Scattering; Brown, W., Ed.; Oxford University Press: Oxford, 1993; Chapter 2, pp 76-148.
S0743-7463(96)00631-2 CCC: $12.00
It is reviewed and extended by Scha¨tzel4 and by Peters3 in separate chapters of a recent text. The first purpose of this paper is to show that, by using more powerful PCtype computers than previously used and new computational tools such as Mathematica, one may rather quickly obtain the full covariance matrix of photon autocorrelation functions in a format which is easily transportable across platforms. A second purpose is to study the effects of varying pinhole size on the covariant matrix and on fitted parameters in samples containing both a single population and a bimodal population of polystyrene spheres. Lastly, we will discuss scattering from monodisperse samples which exhibit multiple scattering and present and discuss the effects of multiple scattering on both the covariance matrix and the resulting analysis of the correlation functions. II. The Variance and Covariance Matrices of Photon Autocorrelation Functions The experimental results of a typical homodyne lightscattering spectroscopy experiment produce the secondorder correlation function of the spatially integrated intensity, which has the following form:
G2(τ) ) 〈I〉2(1 + βg2(τ))
(1)
where
∫0∞ exp(-Γτ) p(Γ) dΓ]-2 + noise
g2(τ) ) [
In these expressions 〈I〉2 is the average scattered intensity squared, β is the intercept of the correlation function at τ ) 0, g2(τ) is the experimental normalized second-order autocorrelation function (henceforth referred to as “the correlation function”), Γ is the decay time, and p(Γ) is the decay time distribution. The goal of data analysis is to obtain p(Γ) from the experimentally known g2(τ). Clearly, an understanding of the magnitude and sources of the (5) Scha¨tzel, K. Quantum Opt. 1990, 2, 87.
© 1996 American Chemical Society
6692 Langmuir, Vol. 12, No. 26, 1996
Harrison and Fisch
various contributions to the noise is necessary to obtain a meaningful inversion of g2(τ) to obtain p(Γ). Peters3 notes that there are two ways to determine the variance of the signal in such autocorrelation functions. The first is to make a large number of independent measurements of the correlation function and then calculate the actual variance at each sample time τ. Obtaining such data is rather time consuming; however, the data analysis is extremely fast. In the present study we calculated the variance of each channel using this technique and found that the calculated variances exhibited sufficiently large variations from a smooth curve (even when fifty 30 s runs were averaged). After comparing the first technique with the second technique of calculating the variance from a model, we decided that the first technique, even if there were no correlated noise, is not as useful as the calculations based on noise models. The second technique uses a model which incorporates both photon and intensity noise and the experimentally known g2(τ) to calculate the (diagonal) variance matrix. Assuming a serial digital correlation with the independent channels labeled with the index k, such that τ ) k∆t, and by using the notation of Peters, this matrix may be written as
∑i |χi|4 + ∑i |χi|2|χi+2k|2 + 2β|χk|2∑|χi2| + 2βχk2∑χiχi+2k + 2β∑χiχ2i+kχi+2k + i i i 4β|χk|4∑|χi|2 - 8β|χk|2χk∑χiχi+k + 2β-1n-1[|χo|2 + i i
var(g2(k)) ) M-1β2{
|χ2k|2 + 2β|χk|2 + 2βχk2χ2k - 2β|χk|4] + β-2n-2[1 + β|χk|2]} (2) In this expression M is the total number of samples collected at each channel, n is the count rate per sample time ∆t (when a multiple sample time correlator, also known as a “muti-τ correlator”, is used this will vary in different subcorrelators), k is the channel number, and β is the experimentally determined intercept. The terms labeled by the various χ’s are assumed to be real for diffusing particles and are defined as follows: exp 2 χk ) xgexp 2 (k), provided g2 (k) > 0, |χk| ) -2 gexp 2 (k), and |χi-j)0| ) (∆t)
∆t |χ(t)|2(∆t - |t|) dt ∫-∆t
(3)
in these three expressions we have assumed that the noise is small enough that the ideal correlation function χk is very well approximated by the experimentally determined correlation function gexp 2 (k). We have also included the “triangular averaging” that will occur in our estimate of the zero time correlation function. All of our data were obtained using a commercial correlator that allows configuration in both a serial and multi-τ configuration. When data is analyzed from multi-τ runs, the above form must be modified to include the variation of sample time and count rate in the different octaves of the correlation function. This will be discussed in greater detail later. The expression for the covariance matrix is somewhat more complex; using the same notation as above this may be expressed as
∑i |χi|2|χi+k-l|2 + ∑i |χi|2|χi+k+l|2 - 4β|χk|2 χl∑i χiχi+l + 2βχkχl × (∑χiχi+k-l + ∑χiχi+k+l) + 2β∑χiχi+kχi+lχi+l+k + i i i
Cov(g2(k)g2(l)) ) M-1β2{
2β-1n-1[|χk-l|2 + |χk+l|2 + 2βχkχk-lχl + 2βχkχk+lχl -
∑i |χi|2 - 4β|χl|2χk ∑i χiχi+k +
2β|χk|2|χl|2] + 4β|χk|2|χl|2
β-2n-2δkl[1 + β|χk|2]} (4) where Cov means covariance and the elements of the matrix are labeled by k and l. Peters3 discusses the importance of calculating χm in the above equation using the above forms for χm that include triangular averaging. The above expressions are valid for a serial correlator, a correlator in which there are n channels all with the same constant sample time of ∆t. In practice, this clock time is chosen so that the characteristic decay time of the autocorrelation function τc is ,∆t. A typical procedure is to pick ∆t so that the correlation function falls to a predetermined fraction of its initial value, τ ) 0, at the last channel. When the correlation function contains more than one characteristic decay times, this is not the optimal arrangement of the correlator, because a sample time which is appropriate for a long decay time will be too long for a short decay time, and vice versa. This potential problem has been overcome by introducing mutiple-tau techniques, in which the correlator is configured as a series of serial correlators with different channel times and with overlapping channels. It is common to use eight channels in each “subcorrelator” and to have the sample time in each subcorrelator be twice that of the preceding correlator. As discussed by Peters3 and others,4,5 this leads to a minimal amount of triangular averaging distortion and a reduction in photon noise by a factor of 21/2 for each block of eight channels. Calculating the variance and covariance for a multiple-τ correlation function is much more complicated than calculating those discussed for a serial, or single sample time, correlator. The variation of both the mean count rate and the number of samples with subcorrelator is challenging because of both notational problems and the fact that much of the needed information is not transparently available as an output in our Malvern System 4700 correlator. There are minimal triangular averaging effects on the variance, even when using the muti-tau mode. The last term in eq 2, in which a χ0 ≡ χn-p)0 term appears, is the only one which is affected, and the calculation of this term is numerically straightforward. It is important to realize that this term approaches 1 for ∆t much less than τc and is less than 1 for ∆t ≈ τc. The most significant changes in the covariance matrix are in the off-diagonal elements. The largest contribution is in the terms in which a combination of subscripts, i, j, k, and l sum to zero. The details of this averaging are available in Peters’ chapter in Brown’s text.3 The implementation scheme for the Malvern System 4700 correlator, used in the present experiments, is described in the Appendix. III. Results and Analysis Apparatus and Experimental Details. The data described in this paper were obtained using a Malvern System 4700 correlator and goniometer and a Lexel Ar+ laser operating at a wavelength of 5145 Å and a constant power of approximately 0.15 W. The intensity was adjusted downward using a half-wave plate and Glan-
Photon Autocorrelation Functions
Langmuir, Vol. 12, No. 26, 1996 6693
the run, g2i(k) is the kth channel of the ith run of the correlation function, and avg(g2i(k)) is the average value of the kth channel of the correlation function
avg(g2i(k)) )
1
N
∑g2i(k)
N i)1
where N is the number of runs, i is the index number of
The much smoother curve is obtained from the same data by using the noise model for the variance, eq 2. These two calculations of the variance have the same general functional form. However, the variance predicted by the model appears to be approximately two times smaller than that determined directly. Checks of the computer code and studies of the noise from the laser and from other sources indicate that these are not the problem. Nevertheless, the very near agreement gives credence to the correctness of both the model and the computer code used to calculate it. The fact that there is extensive covariance between data obtained at different decay times may be seen in several ways. The first is to graph the quantity Dj(k) for any run j vs channel number k. A typical graph of this type, Dj(k) vs channel number, is shown in Figure 2b. The important observation is that the deviation of this and the jth correlation function from the mean correlation function is correlated. That is, when the value of Dj(k) in one channel is above the mean, the data in adjacent channels are also very likely to be above the mean and vice versa. A different run, from the same sample and experiment, shown in Figure 2a, shows a variation that is more like what one might expect; the variation appears random and of zero mean. A very convenient way to visualize this correlation in the deviation is to construct a lag plot.10 Such plots are shown in Figures 2c and 2d for the data shown in Figures 2a and 2b, respectively. If the Dj(k) were totally random, such a plot would be a random distribution about the origin. A linear fit to such data would obtain an intercept, slope, and coefficient of determination10 R2 very close to zero. The coefficent of determination R2 is a measure of the dependence of one variable to another: a value of one implies complete dependence, and a value of zero indicates no dependence.11 The data shown in Figure 2a and lag plotted as Figure 2c obtain a coefficient of determination of 0.002. The data of Figures 2b and 2d obtain a coefficient of determination of 0.21. A second technique, which is very similar to the one used to calculate the variance, is used to experimentally determine the covariance matrix as the straightforward generalization of eq 5. Such calculations are neither computationally difficult or very time consuming. However, they result in an extremely noisy covariance matrix. The third technique is to use a noise model and calculate the covariant matrix using eq 4. This technique yields a covariant matrix that is sufficiently smooth to be used in later calculations. A matrix obtained for 110 nm polystyrene spheres is shown in Figure 3. The experimental details are included in the figure caption. These experiments were repeated with the correlator operating as eight 16-channel subcorrelators. The code was modified as discussed in the Appendix, and the data were analyzed to determine their variance and covariance. The results of a calculation of the variance using both eqs 2 and 5 are shown in Figure 4. Notice once more the ability of the model to accurately predict the experimentally determined variance, only with significantly less
(6) Haller, H. R.; Destor, C.; Cannell, D. S. Rev. Sci. Instrum. 1983, 54, 973. (7) Koppel, D. E. J. Chem. Phys. 1973, 57, 4814. (8) McWhirter, J. C.; Pike, E. R. J. Phys. A: Math. Gen. 1982, 11, 1729. (9) Provencher, S. W. Comput. Phys. Commun. 1982, 27, 213.
(10) Bates, D. M.; Watts, D. M. Nonlinear Regression Analysis and Its Applications; Wiley: New York, 1988. (11) Baird, D. C. Experimentation, 2nd ed.’ Prentice Hall: Englewood Cliffs, NJ, 1988. (12) Schmitz, K. S. Macroions in Solution and Colloidal Suspension; VCH: New York, 1993; p 342.
Figure 1. Variance vs channel number for a serial correlator. The sample consists of 110 µm polyspheres in water; the pinhole diameter is 500 µm.
Thompson polarizer as described by Haller et al.6 The data were collected at a constant scattering angle of 90° and temperature of 25 °C. Typically fifty 30 s runs were obtained for each sample. Whole runs of data were rejected on the basis of two criteria. The first criterion required the calculated and measured far-point (or base line of the correlation function) be identical within 0.5%. The second was a visual inspection of a surface plot of the correlation functions vs run number and channel k. Any run which lay substantially above or below the rest was discarded. Typically 47 or 48 runs were analyzed, except in very dilute samples where the number might fall to 20-30. Similar data were obtained at other scattering angles; these studies were in the context of other experiments and will be reported elsewhere. The samples consisted of polystyrene spheres (Interfacial Dynamics Corp.) that were diluted in water, obtained from a Millipore Milli-Q water purification system, to the desired concentration. Except in those experiments in which multiple scattering was explicitly desired, the concentration of spheres was adjusted so that multiple scattering effects were negligible. Variance and Covariance. The effects of correlated noise and the existence of a nondiagonal covariant matrix were not generally included in standard analysis packages such as cumulant,7 multiexponential,8 and CONTIN.9 These standard analyses assume that the noise associated with the signal at each channel is uncorrelated with the noise at any other channel. In the following it will be demonstrated that the existence and magnitude of such correlated noise, while directly discernible from studies of the variance of correlation functions, are more easily and precisely obtained from the noise models presented above. The variance of a typical normalized correlation function, obtained from a serial run with 128 channels, is graphed vs channel number, which is proportional to time, for a monodisperse solution of 110 nm polystyrene spheres in water in Figure 1. This figure consists of points which are the variances obtained experimentally and a smooth curve which is based on the model discussed above. The data, which appear to represent a noisy decaying exponential, are determined from the variance obtained by averaging fifty 30 s runs and then calculating the variance as
var(k) )
1
N
∑(Di(k))2, with Di(k) )
N - 1 i)1
g2i(k) - avg(g2i(k)) (5)
6694 Langmuir, Vol. 12, No. 26, 1996
Harrison and Fisch
Figure 2. (a and b) Deviation D(k) between an individual run and the average run vs channel number k. (c and d) Lag plots of the deviations shown in a and b, respectively. The experimental conditions are the same as those in Figure 1.
Figure 3. Covariance for the same experimental conditions as those in Figure 1.
Figure 4. Variance vs channel number for a multi-τ correlator. The sample is 110 µm polyspheres in water; the pinhole diameter is 400 µm.
random variation from channel to channel than the direct calculation. Figure 5 is a graph of the Dj(k) quantities and their corresponding lag plots for two different correlation functions obtained from the same set of sample
runs in the parallel mode. The first, shown in Figure 5a and as a lag plot in Figure 5c, is fairly typical and shows the same type of correlated difference from the mean between different channels that was seen in Figure 2. For the particular run shown, the coefficient of determination is 0.48. The runs in this set obtained coefficients of determination between 0.01 and 0.48, with an average value of 0.21. The run which obtained a coefficient of determination of 0.08 is shown in Figures 5b and 5d. The difference between the run shown in Figure 5a and analyzed in Figure 5c and the runs of Figures 5b and 5d is obvious and clearly indicates that even at the longer decay times associated with the parallel correlators the correlation in the noise does not decrease substantially. The logarithm of the covariance for this measurement is shown in Figure 6. Notice once more that there are substantial off-diagonal elements in this matrix. In summary, there is significant correlation between the noise in different channels. This noise is clearly obtainable in both parallel and serial correlators and leads
Photon Autocorrelation Functions
Langmuir, Vol. 12, No. 26, 1996 6695
Figure 5. (a and b) Deviation D(k) between an individual run and the average run vs channel number k. (c and d) Lag plots of the deviations shown in a and b, respectively. The experimental conditions are the same as those in Figure 4.
Figure 6. Log covariance vs channel number for a multi-τ correlator. The experimental conditions are the same as those in Figure 4.
to a covariance matrix in which there may be significant off-diagonal elements. Pinhole Effects. Scha¨tzel5 predicts that, contrary to common wisdom,13 larger apertures collecting approximately 10 coherence areas will lead to a decrease in higher order relative moments and hence reduce the estimator variance for the autocorrelation function, improve the signal-to-noise ratio, and result in a more diagonal covariance matrix. To better quantify this in the present instrument, the covariant matrix was determined for a range of pinhole and hence coherence area sizes. In these measurements, 120 nm polystyrene spheres in water were used, and to keep dead time effects and the like constant, the mean count rate was kept constant to better than 10%. To visualize the effect of pinhole size on the covariance matrix, cross-diagonal slices, from back left to (13) Maret, G.; Wolf, P. E. Z. Phys. B: Condens. Matter 1987, 65, 409.
front right, of the matrix were made. These slices, for pinhole diameters between 50 and 500 µm, are shown in Figure 7 for the serial correlator. Note that the nondiagonal part of these matrices was essentially unchanged by pinhole diameter. The major effect of increasing the pinhole size is to increase the magnitude of the diagonal (variance) terms with respect to the off-diagonal terms which indicate the covariance of the signal. For a pinhole 500 µm in diameter, the diagonal terms are found to be approximately six times as large as the largest off-diagonal term. This suggests, as does the analysis below, that larger pinholes leading to a value of β of approximately 0.15 are preferable in most analysis schemes. Moreover, a larger pinhole size is more ideal for static light-scattering experiments and allows the incorporation of lower power lasers into a given experimental apparatus. The results of this series of experiments are summarized in Figure 8, where the ratio of a typical diagonal term to the largest off-diagonal term is plotted vs β. The region 0.15 e β e
6696 Langmuir, Vol. 12, No. 26, 1996
Harrison and Fisch
Figure 7. Cross sectional cuts of covariant matrices vs pinhole diameter. The particle size is 110 nm.
Figure 8. Variance divided by the maximum covariant term as a function of β, the signal-to-background ratio.
0.3 shows a good compromise between a reasonable value of β and large values of diagonal to off-diagonal terms. Multiple Scattering Effects. In recent years there has been an increase in interest of trying to understand the effects of multiple scattering on the inversion of lightscattering data.13-17 The effects of multiple scattering were investigated by studying polystyrene sphere/water samples that had normalized light transmissions in the range from approximately 40 to 95%, at λ ) 5145 Å; such solutions look somewhat cloudy when a light is viewed through them. Multiple scattering in these sample was verified by measuring the normalized transmission of light and comparing that to the transmission calculated by assuming single scattering described by the Mie theory.18 Multiple scattering was presumed to exist when the measured transmission was less than that calculated using the known number of particles per cm3 and the scattering cross section of the particles.19 As the number of spheres per unit volume became smaller, the samples transmitted more light, and the calculated and measured values approached one another as expected. Samples which transmit g95% of the incident light were essentially single scattering. The effects of multiple scattering on the covariance matrix are shown in Figure 9, which has calculated covariance matrices from data obtained at a scattering angle of 90° and a pinhole diameter of 300 µm for samples which transmit 38, 92, and 95% of the incident (14) Edrei, I.; Kaveh, M. J. Phys. C: Solid State Phys. 1988, 21, L971. (15) Bailey, A. E.; Cannell, D. S. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1994, 50, 4853. (16) Dhont, J. K. G. Physica 1983, 120A, 238. (17) Vrij, A.; Jansen, J.; Dhont, J. K. G.; Pathmamanoharan, C.; Kops-Werkhoven, M. M.; Fijnaut, H. M. Faraday Discuss. Chem. Soc. 1983, 76, 19. (18) Bohren, C. F.; Huffman, D. R. Absorption and Scattering of Light by Small Particles; Wiley: New York, 1983; pp 83-129. (19) van de Hulst, H. C. Light Scattering by Small Particles; Dover: New York, 1981; pp 33 ff, 338 ff.
Figure 9. Variance vs channel number for a 110 µm polyspheres in water for several concentrations of spheres (‚‚‚, 95% of the incident light transmitted; s, 90% transmission; - -, 38% transmission).
light. There is remarkably little variation in these matrices. The effects on the dynamics of the particles are much greater.13,17 As an aside, the mean scattered light intensity was measured as a function of scattering angle and the resulting data fit the predictions of the Mie theory. Even in the most concentrated solution the difference between the fit and the data were quite good. This suggests that even in highly scattering solutions the size of a monodisperse population of particles may be determined by static light scattering. IV. Data Fitting Overview. One of the goals of the present work is to find the appropriate weights to use to invert the measured correlation function and hence to find the most appropriate fitted parameters and estimates of their errors. The goal of least-squares fitting is to, in a robust manner, find the global minimum of M
χ2 )
∑ k)1
M
mod 2 [wk(gexp 2 (k) - g2 (k)) ] )
[wk(∆k)2] ∑ k)1
(6)
is the where M is the number of data points used, gmod 2 model function for the correlation function (cumulants, biexponential, etc.), and wk is the weight of the kth data point. For noncorrelated noise wk ) 1/σ2k, and ∆k ≡ mod gexp 2 (k) - g2 (k). Furthermore, the model function contains parameters that are to be determined by the fit. This formula may be rewritten in matrix form as
χ2 ) ∆Tk W∆k
(7)
where the weights wk have been replaced by a diagonal weight matrix, W ) wkδjk, where δjk is the Kronecker δ, and where the ∆k now forms a 1 × M vector. This suggests
Photon Autocorrelation Functions
Langmuir, Vol. 12, No. 26, 1996 6697
the following generalization to include a covariant matrix with off-diagonal terms such as used in the present analysis:
χ2 ) ∆Tk Cov-1∆k ) ∆Tk W∆k
(8)
where Cov is the covariant matrix calculated above and Cov-1 ) W is the weight matrix.20 For strongly diagonal covariant matrices the corresponding inverse matrix (weight matrix) is also strongly diagonal. A program to implement this minimization on Mathematica was written and used to obtain the results discussed below. Analysis of Mondisperse Populations of Particles. To obtain a better understanding of the effects of nondiagonal weights on the performance of the present fitting routine, the following analysis was undertaken. Lightscattering data, the average of 50 serial runs obtained using various pinhole diameters obtained from a dilute 120 nm polystyrene sphere/water sample, were analyzed using the full (calculated) covariant matrix, the calculated diagonal variance matrix, and a diagonal variance matrix in which all the elements have the same value, the average of the diagonal variance matrix. Because the polystyrene spheres are rather monodisperse, the stated variance on the bottle was 8.6%, and the model used for g2, gmod 2 , was chosen to be the three-term nonlogarithmic version of the cumulants7 expansion (see eq 9). In this expression, K1, 2 1 3 gmod 2 (τ) ) A exp(-2K1τ + K2τ - /3K3τ )
Figure 10. First cumulant K1 vs pinhole diameter for fits using only the variance and the covariance. The data have been dithered for better visibility, and error bars representing the uncertainties obtained in the fitted parameters are shown. An error of 1% is shown on the right of the figure (2, variance matrix; [, covariant).
(9)
K2, and K3 are the cumulants in terms of the average diffusion coefficient, and we have self-consistently (as in our definition of gexp 2 ) had the model approach zero as τ f ∞. The results of such an experiment are interesting. First, there is a very small systematic variation in the parameters A and K1 with the type of weight matrix used. The parameter A is systematically smaller when the full covariant matrix is used by between 0.1 and 0.2%. However, there is little difference between the diagonal and average diagonal covariant matrix schemes. This matrix type variation should be compared to a typical uncertainty in A which is approximately 0.002 (A is typically very close to 1.000) regardless of the type of weighting. Thus this is a very minimal effect. The typical uncertainty in the first cumulant depends on the weighting scheme. The first cumulant K1 obtained from fits to the full covariant matrix is systematically smaller than that obtained by the other two covariant matrix schemes. The value of K1 also depends upon whether K3 is zero or not. In those fits in which K3 is zero (by either forcing it to zero or through the fitting procedure), K1 is roughly 1% larger than in those cases in which K3 is nonzero regardless of the weighting scheme. The uncertainty in K1 depends on the pinhole size and is always between 1.5 and 2 times as small when the diagonal covariant matrix schemes are used as that obtained using the full covariant matrix. The uncertainty in the first cumulant, as determined by the curvature matrix, is between 0.5 and 1%; furthermore, the uncertainty of the first cumulant is rather bimodal. The larger pinhole sizes, 300 and 400 µm, obtain nearly identical errors which, while dependent on the weighting scheme, are essentially independent of pinhole size. The smaller pinholes (200, 100, 50, and 35 µm) also have uncertainties that depend on the weighting scheme. However, in all cases the uncertainties are roughly one-half of those found using (20) Lupton, R. Statistics in Theory and Practice; Princeton University Press: Princeton, 1993; pp 81-97.
Figure 11. Second cumulant K2 vs pinhole diameter for fits using only the variance and the covariance. The data have been dithered for better visibility, and error bars representing the uncertainties obtained in the fitted parameters are shown. An error of 10% is shown on the right of the figure (2, variance matrix, [, covariant).
the larger pinholes. These results are shown graphically in Figure 10. The situation for the second and third cumulants is much more dramatic. The values of both the second and third cumulants depend rather strongly on both the weighting scheme and the pinhole size. This can be a rather serious effect as the polydispersity parameter P often defined through P ≡ K2/K21 depends upon the pinhole size.. Thus the apparent polydispersity of a sample depends on pinhole sizesa very unsatisfactory situation. First, we will consider the situation in which K3 is forced to equal zero. In this case, K2, the second cumulant, is a very strong function of pinhole size, decreasing from approximately 375 000 to 112 000 s-2 as the pinhole diameter is decreased from 400 to 35 µm. This corresponds to the polydispersity parameter decreasing from 0.08 to 0.03. Another significant observation is that the uncertainties in the parameters are also a strong function of the weighting scheme. The uncertainty in K2 is approximately 15-33% for the full covariant matrix while it ranges from only 7 to 9% for the diagonally weighted fits. These trends are shown in Figure 11. The fits and the resulting uncertainties also change when K3 is included in the model. First, and not at all surprisingly, the uncertainties in K1 and K2 both increased dramatically. More significantly, K2 decreased by approximately an order of magnitude. From a practical point of view, the best fit value of K3 is zero for diagonal weighting at the smaller pinhole sizes. For the 300 and 400 µm pinholes, the values of K1-K3 are all very close (0.05% in K1, 2% in K2, and less than 3% in K3) regardless of the weighting scheme. When observed as a whole this analysis strongly suggests the following. First, the optimal pinhole size, in
6698 Langmuir, Vol. 12, No. 26, 1996
Harrison and Fisch
the present apparatus is roughly 300-400 µm. This corresponds to an intercept in the correlation function β of roughly 0.1-0.15. This verifies Scha¨tzel’s proposition. Second, for these larger pinhole sizes the variance dominates the covariance sufficiently so that the weighting schemes which utilize only the diagonal weights produce optimal fit values that are indistinguishable (for nearly monodisperse populations) from those produced using the whole covariant matrix. Third, the situation with regard to the uncertainties of the fitted parameters is more problematic. The covariant matrix produces larger uncertainties, by 1.5-2 times greater than those of the diagonal covariant matrix schemes. Fourth, in the case of fairly mondisperse polystyrene spheres in water the polydispersity determined at larger pinhole sizes appears to be a more correct estimate of the width of the particle distributions. Analysis of Bimodal Populations of Particles. In many situations the population of the scatterers is unknown. For this reason we applied a modification of this analysis scheme to polystyrene sphere/water samples which contain two different, nearly monodisperse, populations of particles. The present experiments used mixtures of 73 and 221 nm spheres, mixtures of 87 and 507 nm , and mixtures of 120 and 541 nm particles. These size ratios were studied because the former is near the minimum that we hoped to separate, while the latter two should be easily distinguishable as two populations. The measured hydrodynamic diameters and polydispersity parameters are summarized in Table 1. The number of particles per unit volume was adjusted so that nearly identical mean scattered intensities were individually obtained, from each component of a mixture, at a scattering angle of 90°. Once more, 50 runs were made at each pinhole size. Typically 30 runs were accepted. Because the characteristic relaxation time of the two populations differs by at least a factor of three, the parallel correlator mode was utilized. The data from the sample with a population ratio of approximately 3.1:1 were analyzed in several ways. First, the data were analyzed using the Malvern version of CONTIN and the Malvern exponential sampling, the model free analysis scheme, with the weighing schemes inherent to these commercial packages. The CONTIN analysis almost always predicted a broad monodisperse population of particles. On occasion a bimodal population was the optimal fit. However, in those cases, one population was always very small, roughly 25-30 nm, and represented less than 6% of the scattered intensity. The other population of particles was very polydisperse and had a mean radius near that obtained by a cumulant fit. The exponential sampling scheme gave either a bimodal population, in which the smaller sized population was predicted to be too large and the larger sized population was predicted to be too small, or a unimodal peak similar in mean value and width to that obtained from all cumulant fits and CONTIN fitting. This value was roughly 125 nm, the geometric mean of the two mean particle sizes. The exponential sampling scheme gave relative intensities that were closer to the know intensities than those obtained by CONTIN analysis, differing by only a factor of 2. The present analysis scheme was used in the following modified manner. The data was first analyzed using the cumulants scheme discussed above. The data was then fit to the model correlation function of eq 10. where the “intensity constants” Ai, the inverse decay
Table 1. Measured Particle Sizes and the Corresponding Polydispersities of Polystyrene Latex Spheres particle diam (nm)
polydispersity param
particle diam (nm)
polydispersity param
72.6 87.0 119.8
0.013 0.007 0.038
221.4 503.0 543.0
0.002 0.032 0.020
the number of terms N was a variable between 2 and 5. The constant c was included to allow for the fact that the data did not always decay totally to zero, and its inclusion is suggested in the work by Stock and Ray.21 Using both the variance and the full covariance matrix, the fits almost always converged to two populations with mean inverse decay times very close to those predicted from the known mean particle sizes. A further criteria for a good fit was to exam a lag plot of the residuals: any large correlation in the residuals, as indicated by a large slope or intercept, could be clearly identified and suggested that the “best” fit had not yet been obtained. This criterion was applied to all fits in the present analysis. These results are summarized in Table 2. Table 2 has two parts. Part A summarizes the mean size of the two populations, 〈Dsmall〉 and 〈Dlarge〉, where D is the diameter, as determined from eq 10 for the various number of exponentials N. In all cases the best fit is obtained from the N ) 2 sum of exponentials fit. Note also that the mean sizes obtained from the fit of the two population sample are very close to the values measured independently, 73 and 221 nm, respectively. In the one case where the N ) 3 fit is as good, it is clear that the population of larger sized particles was split into two populations whose geometric mean is 222 nm. Thus the present analysis scheme can identify a bimodal population, when presented with the possibility of multimodal populations. An equally important test of a fitting routine is its ability to correctly distinguish between a broad unimodal distribution and bimodal distribution with a small size ratio. In particular, it is of interest to determine the minimum size ratio of the two populations of a bimodal distribution that a given analysis can separate. To determine this ratio in the present analysis scheme, we fit the same data to the cumulants form, eq 9 with K3 ) 0 and with K3 * 0, and the following form which allows for a bimodal population, where each population has a width: 2 gmod 2 (t) ) [A1 exp(-b1t + K21t ) +
A2 exp(-b2t + K22t2)]2 + c (11) This expression, a straightforward extension of eq 9, has seven adjustable parameters. The fit to this form obtained mean sizes of roughly 71 and 220 nm and polydispersity parameters of approximately 0.01 and 0.03, respectively. Once more the mean sizes are in very good agreement with the individual size measurements; however, the polydispersity parameter is somewhat smaller than that determined on the samples individually. The second part of Table 2, part B, compares the χ2 and F test22 results of the cumulant fit eq 9, the sum of exponential fits, eq 9, and the two exponentials with width, eq 11, schemes. The cumulant fit with K3 ) 0 is so poor that it is not included in the present analysis. The conclusions to be obtained from this table are (a) the twoexponential fit is almost always the best fit, even better than the two exponential with width fit, and (b) the F test
N
gmod 2
Ai exp(-bix))2 + c ∑ i)1
)(
(10)
times bj and the constant c were all free parameters and
(21) Stock, R. S.; Ray, W. H. J. Polymer Sci., Part B: Polym. Phys. 1985, 23, 1393. (22) Bevington, P. R.; Robinson, D. K. Data Reduction and Error Analysis for the Physical Sciences, 2nd ed.; McGraw-Hill: New York, 1992.
Photon Autocorrelation Functions
Langmuir, Vol. 12, No. 26, 1996 6699
Table 2 Measured Particle Sizes and the Corresponding Polydispersities of Polystyrene Latex Spheres A. Predicted Mean Diameters of Sample Containing Both 72 and 220 nm Monodisperse Populations of Particles as a Function of Number of Exponentials and Pinhole Size N)2
N)3
N)4
N)5
pinhole diam (µm)
type of weight
〈Dsmall〉 (nm)
〈Dlarge〉 (nm)
〈Dsmall〉 (nm)
〈Dlarge〉 (nm)
〈Dsmall〉 (nm)
〈Dlarge〉 (nm)
〈Dsmall〉 (nm)
〈Dlarge〉 (nm)
500
variance covariance variance covariance variance covariance variance covariance variance covariance
73.8 73.8 70.5 70.2 72.8 70.8 71.8 72.2 69.4 67.5
251 250 222 220 235 224 226 228.9 219 212.9
71 75.2 62.5, 91.5 65.7 73 63.8, 120 69.3 63.8, 87.4 69.6 67.4
159, 323 161 240 177.3 236.5 97442 163, 303 240.7 219.6 214.9
66, 113 62.4, 103 65.2, 125 66.8 66.8, 122 75 58.6, 95.2 66.7, 146 70.8 63.1, 143
307 299 275 185 290 210 171, 417 5772.8 170.6 52897
67.5, 126 49.4, 99 66.8 49.5, 112 67.6, 135 60.2, 124 74.5, 134 68.8, 138 69 51.5, 118
322 224 142, 280 10755 300.7 56103 262.8 275 238 398.8
400 300 150 50
B. Comparison of χ2 and F Values Using the Covariant Matrix and Both Cummulants and a Sum of Exponentials Fits pinhole diam (µm)
cummulant fitb
500 400 300 150 50
37.3 25.0 33.4 19.2 26.4
2 exponentials χ2 F 38.1 24.0 30.1 18.7 24.1
0.465 0.41 0.46 0.46 0.35
3 exponentials χ2 F 9 × 10-11 0.46 0 0.492 0.41
191 23.6 211 18.7 24.2
4 exponentials χ2 F
5 exponentials χ2 F
2 exponentialc χ2 F
37.7 23.2 37.8 33.9 69.6
40.4 25.7 40.3 18.6 31.9
38.0 23.7 32.2 18.6 29.6
0.38 0.48 0.23 0.006 0
0.24 0.31 0.13 0.44 0.12
0.42 0.44 0.49 0.49 0.28
a The preferred fit, based on χ2 and the F value of a given sum of exponentials, as compared to the cummulant fit is in bold type. b This is the fit to a model in the form of eq 8 in the text. c This is the fit to a model in the form of eq 11 in the text. The fits with N ) 2 had very nearly identical weights for the large and small size populations ((3%). In the other fits the intensity ratio was very nearly 1:1 when there were two populations.
Table 3. Comparison of CONTIN, Malvern Multiexponential, and Cumulants, and 2-Exponential Fitting Schemes for Three Two-Particle Mixtures with Different Size Ratios. The Ratio 3.1:1 Corresponds to a 72 and 220 nm Mixture, 4.5:1 Corresponds to a 110 and 530 nm Mixture, and 5.8:1 Corresponds to a 87 and 507 nm mixture. The Best Fits Are in Bold Type CONTIN
multiexponential
cumulants
2 exponential
ratio
〈D〉 (nm)
width (nm)
〈D〉 (nm)
width (nm)
〈D〉 (nm)
width (nm)
〈D〉 (nm)
P
K3 (s-3)
〈D〉 (nm)
P
〈D〉 (nm)
P
3.1:1 4.5:1 5.8:1
140 253 102.6
68 146.6 16
117 135 90
75 40 68
0 406 196
0 46 30
110.1 192.9 112.0
0.297 0.414 0.33
2.2 × 108 1.2 × 108 4.1 × 108
72.9 116.7 84.3
0.006 0 0
223.5 495.4 465.6
0.018 0.057 0.044
results indicate that even though the sample is known to contain two populations of particles, one cannot statistically determine between a fit to a broad unimodal population and a bimodal population at this size ratio. While we cannot statistically determine this difference, it is important to realize that, because we know the nature of the sample, we can, and did, obtain the correct sizes. The F test results suggest that this size ratio is somewhat smaller than that which can be separated using the present analysis scheme. The final test that was performed was to examine other particle size ratios to determine the minimum size ratio that, based on meaningful statistical criteria, could be separated. Data for mixtures of 110 and 530 nm and 87 and 507 nm spheres were examined using a the same analysis procedures as described above. The results of these tests are shown in Table 3. It is significant to note that the current scheme separates two nearly monodisperse populations with size ratios of 4.5:1 and 5.8:1, while the other analysis schemes were unable to achieve such separation. Once more the best fits are shown in bold. The results of the F test for these two size ratios indicate only a 10% probability that the cummulant fit is better than the two exponential fit.
the signal to background ratio β of the correlation function is roughly between 0.1 and 0.2. At smaller β the signal is too small, while at larger β there is significant covariance. Second, for β in this range the covariant matrix is sufficiently diagonal that there is little difference in the determined hydrodynamic diameter of the particles between the complete covariance and the variance. However, there are sufficient differences in the uncertainties and higher order cummulant terms to more than justify the use of using the covariance matrix. Third, the covariant matrix is changed very little by multiple scattering effects. Fourth, using the covariant matrix one may clearly separate two nearly monodisperse populations of particles that differ in size ratio by approximately 4:1. Fifth, if it is necessary to use a small pinhole, the full covariant matrix, which may be calculated from a single run, should be used in data reduction algorithms.
V. Conclusions
Acknowledgment. We thank Carol Mertz of Argonne National Laboratories for her assistance with the program CONTIN and Professors D. Keith Robinson and Charles Rosenblatt of Case Western Reserve University for their helpful conversations. This work was done under National Science Foundation Advanced Liquid Crystal and Optical Materials Grant DMR8920147 and National Science Foundation Grant DMR9321924.
The following conclusions may be drawn from the present study. First, the optimal pinhole size for quasielastic light scattering experiments should be chosen so that
Appendix To correctly construct a covariant matrix, information, in addition to the correlation function such as the total
6700 Langmuir, Vol. 12, No. 26, 1996
Harrison and Fisch
Figure 12. (a) Covariance vs channel number for a multi-τ correlator without smoothing. The smoothed matrix is shown in Figure 6.
number of samples, the far points, and the signal-tobackground noise ratio, is needed. This information was extracted from the header information of Malvern’s data files, which had to be converted from a binary format to an ASCII format. Obtaining the information in a usable form was somewhat more difficult as Malvern’s software “normalizes” the data after collecting it, and the method of normalization is poorly documented. The ASCII file was then imported to Mathematica where runs were eliminated as previously stated. The data were then parsed into a set containing the above mentioned additional information and correlation function, which was then normalized and finally averaged. Any points that were less than zero in the normalized correlation function were set to zero.3-5 To ensure that the summations in eq 4 were able to be carried out to the number of points in the correlation function, as in the sums where there are combinations of i, k, and l, the normalized correlation function was extrapolated out to three times the number of points used. This was done by fitting the last two octaves to a decaying exponential of the form g2 ) a(e-bt). These points, as well as the extrapolated ones, are close to zero. Hence, this had a negligible effect, less than 1%, on the covariant matrix. The extrapolated normalized correlation function was then used in eq 4 to calculate the covariant matrix. The preceding straightforward calculation must be modified in order to analyze the multi-τ runs. First, triangular averaging effects, as discussed by Scha¨tzel,4,5 must be implemented. The term ∆t |χ(t)|2(∆t - |t|) dt ∫-∆t
|χ0| ) (∆t)-2
(A1)
cannot be estimated as 1 for long lag times,4,5 so it was calculated directly from the data using numerical integration. Because of the symmetry of the correlation function around t ) 0, and to accommodate the fact that χ is always positive, it was necessary to split the integral eq A1 into four integrals. Scha¨tzel states that two other terms,
∑i |χi|4 ) 2(∆t)-2∫2∆t|χ(t)|4 dt + 2(|χl|2)2 + (|χ0|2)2 t
and ∆t |χ(t)|2(∆t - |t - ∆t|) dt ∫-∆t
|χl|2 ) (∆t)-2
(A2)
are affected by triangular averaging, but as we observed less than 1% difference between the integral and the summation evaluations in our samples, we left them as sums. In multi-τ correlation functions the time per channel depends on the particular octave of data. The average count rate must be changed to reflect this time dilation. The average count rate in the jth octave was j 0, where j is the octave number, n j0 changed to n j j ) 2j-1 n is the mean count rate in the first octave, and 2j-1 is the time dilation factor. The total number of samples also changes in the multi-τ mode, depending on the octave. The correct total number of sample in the jth octave is Mj ) jM0, where M0 is defined as total number of samples in the first octave. Introduction of these quantities gives the covariant matrix a steplike look, shown in Figure 12. Since the weight matrix is the inverse of the covariant matrix, this causes some unusual artifacts in the weight matrix To minimize the impact of these steps upon inversion, the covariant matrix was smoothed by a moving average over the rows and columns as shown in Figure 6. Our fitting procedure used eq 8, whereas almost all standard fitting routines use eq 6. A nonlinear fitting routine using the Levenberg-Marquardt algorithm, within the Mathematica package nonlinearFit, was modified to incorporate a matrix of weights, rather than a single weight per point. This was accomplished rather simply by replacing the original χ2 function with the Mathematica equivalent expression of eq 8. The present program is easily modified for other correlators. Modifications allowed data from a Brookhaven BI 9000I correlator to be also analyzed, obtaining results similar to those obtained in the present experiment. The present fitting routine has also been modified to allow CONTIN-type fits with up to 13 adjustable intensities and six adjustable exponentials and amplitudes. Lastly, the present analysis scheme allows noise on a correlation function to be modeled on the basis of a single run. We are making the modified nonlinear fit package available at the Wolfram WWW site. LA9606310
