An algorithm is offered for finding the range within which cumulative particle size distribution functions can be located in consistency with experimental turbidimetric data at a number of wavelengths. It is based on linear programming and minimization techniques. Several tests were performed. The lower righthand branch of the corridor was found to locate near the initial distribution function. © 2014 Optical Society of America OCIS codes: (100.3190) Inverse problems; (290.4020) Mie theory; (290.7050) Turbid media. http://dx.doi.org/10.1364/AO.53.000301
1. Introduction
The problem of retrieving particle size distribution (PSD) functions from optical measurements occupies an important place in colloidal optics [1–3]. Big and small particles scatter light differently. This problem is, in essence, ill-posed and solved by various inversion methods. The Encyclopedia of Analytical Chemistry [3] classifies all inversion techniques into three groups: • Approximating the PSD with a known distributional form. • No assumptions on the shape of PSD. • The method of moments. Another classification [4] also numbers three groups, though slightly differently: • Analytical inversion methods based on the strict solution of Fredholm’s first-kind integral equation with an approximate scattering kernel. • Independent model algorithms using no prior information on the size distribution function, which, therefore, is assumed as a step function, and the
1559-128X/14/020301-05$15.00/0 © 2014 Optical Society of America
heights of steps are calculated by solving a discrete set of linear equations, often with regularization. • Methods taking advantage of some prior information on the size distribution function; in particular, its analytical form with few unknown parameters is preset, the values of these parameters to be estimated by optimization algorithms. All these methods deal with the differential form of PSD as the most useful and suitable for immediate integration to calculate various properties of a system. Every such function obtained from experimental data has a certain accuracy, which, however, cannot be presented in a usual “plus–minus” manner because of the principal property of any differential distribution: its integral over the corresponding range must be equal to unity, and no variations are possible. At the same time, different methods applied to the same experimental data will produce more or less different PSD functions, which together will compose some set. However, to the best of our knowledge, the question of what set of PSD functions would be consistent with some given optical data has been never posed. This is the very set from which different particular algorithms (e.g., the maximum entropy method) select different members (functions) at processing of the “parent” optical data. This question seems important because when few wavelengths are 10 January 2014 / Vol. 53, No. 2 / APPLIED OPTICS
301
selected, transfer from a PSD function to the corresponding turbidity or scattering spectrum is accompanied by information losses, aggravated by the nonmonotonic behavior of Mie’s scattering efficiency factor [5]. The cumulative form of PSD provides more freedom in this respect, having no restrictions on the area under curve. It allows one to outline a corridor enveloping some area with its content consistent with input optical data, with or without experimental errors. The present paper offers an algorithm to process turbidity spectra, which provides limiting cases of the cumulative PSD function (a corridor) and thus envelops the said set. Simultaneously, the maximum values of average particle radius of this or that order (moments of the differential PSD function) can be found.
Consider the first ratio. As α1 λ1 α0 λ0, then, introducing a wavelength ratio κ1
λ0 α1 ; λ1 α0
we obtain for the monodisperse system T1
Kκ1 α0 : Kα0
For the polydisperse system, this looks as R 2 r K2πμ1 r∕λ1 f rdr T1 R 2 r K2πμ1 r∕λ0 f rdr or
2. Theoretical
Experiment provides the value of turbidity τ, which is the coefficient of proportionality in the equation ln
I 0 λ τλl; I l λ
Z and
τ N2π
T 1 K2πμ1 r∕λ0 − K2πμ1 r∕λ1 r2 f rdr 0:
Resubordinating, for the sake of uniformity, the function f r to the argument α0 2πμ1 ∕λ0 r, we obtain
where I 0 λ and I l λ are the intensity of incident and transmitted (optical path l) light, respectively, with a wavelength λ in vacuum. In the absence of optical absorption and multiple (secondary) scattering, the turbidity of monodisperse and polydisperse systems with spherical particles, according to Mie’s theory [5], is τ N 2 πr2 Kα; m
Z
r2 Kα; mf rdr;
Z
T 1 Kα0 − Kκ 1 α0 α20 f α0 dα0 0:
Similar expressions are derived for T 2 …T L−1. Prior to obtaining the limiting distribution curves, it would be useful to find the maximum value of the average radius of particles (here, the relative radius α¯ 0 ): Z
respectively, where N 2 is the numerical concentration of particles, r is their radius, α 2πrμ1 ∕λ is the relative (dimensionless) radius, μ1 is the refractive index of the dispersion medium, K is Mie’s scattering efficiency factor represented by a sophisticated function (the sum of an endless series), m is the relative refraction index, and f r is the numerical differential PSD function with radius as its argument. Formally, the integration is done, here and further, from zero to infinity; actually, the upper limit is restricted by big, quickly sedimenting particles. For brevity’s sake, we will omit the known value m in the arguments of K. Turbidity is measured at several (L) wavelengths, numbering to start from zero (e.g., for a Russian PEC-56 photoelectrocolorimeter, λ0 400 nm). To separate N 2 and r, let us take advantage of the turbidity ratio method [6,7] with λ0 as the base wavelength; (L–1) ratios thus result: T1 302
τλ1 ; τλ0
T2
τλ2 τλ0
T L−1
τλL−1 : τλ0
APPLIED OPTICS / Vol. 53, No. 2 / 10 January 2014
α¯ 0
α0 f α0 dα0 ;
(1)
which is easily recalculable: r¯ λ0 ∕2πμ1 α¯ 0 . This will help to estimate the range that the mentioned functions should be found within. Besides, Kourti et al. [8] mention cases in which an infinite number of distributions could adequately explain their optical data, but all these solutions corresponded to the same weight-average radius. Hence, the range of such radii may also be of interest. The task can be stated as follows: such a function f α0 is to be found that maximizes functional (1). But, as this functional contains no derivatives of the sought function, no variational task (which would suggest itself) can be formulated. Therefore, let us first consider a particular case of a discrete polydisperse system with the aim of future generalization to a continuous f α0 . Such a system comprises k fractions with particles with their relative radii αi and numerical fractions wi, i 1…k. The integration reduces to simple summation, and the task is formulated as the following set:
8P > wi 1; > > P > > > T 1 Kαi − Kκ1 αi α2i wi 0; > > > P < T Kα − Kκ α α2 w 0; 2 i 2 i i i > … > > > P > > T L−1 Kαi − Kκ L−1 αi α2i wi 0; > > > : P α w α¯ → max : i i 0
2
Note that, by physical sense, all wi ≥ 0. Set (2) is easily recognizable as the problem of linear programming (see, e.g., [9]) with respect to the variables wi (they are usually designated as xi in the theory of this method) with a given function Kα and a given value set αi . Linear programming has been used to solve inverse problems [10,11]. Our algorithm is mainly based on Platz and Gordon’s ideas [10]. It is known from the theory of this mathematical discipline that any basic feasible solution, including that providing the maximum, can contain, at most, as many positive (nonzero) components as the number of equations in the set (L). Hence, the average radius at a given value set of T of a discrete polydisperse system is maximum for a system with the number L of fractions or less. Any other moment (α2 , α3 , etc.) will also reach maximum for some system with L or fewer fractions. Now let us increase the number k of points farther to infinity within the same range in order to eliminate the restriction of discreteness and to pass to a continuous distribution function f α0 . This operation, however, will not cancel the above conclusion of the number of nonzero values wi, and in the most general case (the multidimensional problem of linear programming with the space dimensionality tending to infinity), the sought function shall be the sum of L or fewer delta functions. We note in passing that an increased number of points would enable one to achieve a higher maximum (or a deeper minimum). The only difference is in the fact that, in the case of a discrete polydisperse system, the set of αi is given and few members from it are to be selected, while in the case of a continuous function f α0 , these values are to be found within the full range of α0. Therefore, the problem reduces to double maximization, namely, an external one performing by the welldeveloped techniques of linear programming (e.g., the simplex method [9]), and an internal one being the maximization of a function of L variables (the radius values αi ), whose algorithms are also well known [12]. When programming, we did as follows (see also [10]). The maximum number k of fractions was taken divisible by L, with a multiplicity factor of 5 or 7. At the first iteration, αi’s were uniformly spread over a physically reasonable range [0, αmax . The value of αmax was taken small enough (∼5), and then it was increased only if necessary—when the procedure of linear programming failed or produced a nonzero value at the last point (i.e., the optimum might lie
to the right of it). Multivalued solutions described in [13] were thus avoided. At every subsequent iteration, the points with nonzero values wi were “surrounded” with additional “empty” points at the expense of those where wi had been zero at the preceding iteration. Such “surrounding” (an analog of the method of “golden section” or dichotomy [12]) was done more and more densely, and, upon achieving a certain pregiven value of spacing between points, the iteration process finished. Having estimated the maximum average radius and, on its basis, a reasonable range for seeking distribution functions, let us now turn to finding these functions. This could be easily accomplished for cumulative PSD functions designated as Fα0 : in set (2) the expression to be maximized is replaced by
X
U − α0 − αi wi Fα0 → max;
(3)
where U − is a step function [14, 21.9–1], which, in full accordance with the definition of the cumulative PSD function, leaves only those points in the sum, for which αi ≤ α0 , where α0 is the independent variable of the function Fα0 , which changes with a certain step on the range chosen with due account of the maximum α¯ 0 . The suggested algorithm with the usage of linear programming is also applicable when turbidity spectra are processed by Heller’s method [15], which uses taking the logarithm of turbidity for N 2 and r separation; the wavelength exponent n is introduced: n−
∂ ln τ ∂ ln Kα; m ; ∂ ln λ ∂ ln α
which is easily estimated from the double logarithmic ln τ versus ln λ plot and further is taken as the response function. For a polydisperse system, we have [16] R n¯
nαα2 Kαf rdr R 2 ; α Kαf rdr
where nα is the wavelength exponent for the monodisperse case; when integrating with respect to radius, the relation of α with r is taken into account (at the wavelength at which n¯ has been calculated). Hence Z
n¯ − nαα2 Kαf αdα 0:
And in set (2), all the equations with T involved are replaced by one equation: X
n¯ − nαi α2i Kαi wi 0:
However, this reduction of the number of equations affects the precision, especially if n is estimated by numerical or graphical differentiation. 10 January 2014 / Vol. 53, No. 2 / APPLIED OPTICS
303
If the volume fraction Φ of the disperse phase is to be estimated (e.g., in spectroturbidimetric titration of polymer solutions), 4 Φ N2 π 3
Z
r3 f rdr;
then the specific turbidity g [17] can be used: R 2 α Kαf rdr τ 3πμ1 · R 3 2λ Φ α f rdr R 2 α Kαf αdα 3πμ1 : · R 3 2λ α f αdα
gα
The resulting problem of linear-fractional programming reduces to the usual one: 8 P 3 > αi wi 1; >
> : P α2 Kαi wi → max : i
One will have no troubles in rewriting this set to involve T ’s instead of n. Turbidity spectra were simulated with the usage of Rosin–Rammler’s functions [18]:
But the almost zero minimum value of the average radius unambiguously determines the corresponding limiting Fα—starting from the origin of coordinates, it goes upwards along the ordinate axis almost to the F ≈ 1 point, and then rightward, almost horizontally, even slightly rising at those radius values that are necessary to keep the preset T values. In view of a rather wide corridor for F (see below), the said behavior detail of the left-hand upper bound of this function can be neglected and it can be thought of as Γ shaped. Therefore, the task reduces to finding the lower right-hand bound of the corridor of F [the sign “minus” before the sum in Eq. (3)]. Figure 1 shows the source cumulative PSD curve and the calculated right-hand lower limiting PSD curve for the set of wavelengths of the PEC-56 photoelectrocolorimeter: 400, 440, 490, 550, and 585 nm. Though being similar by appearance, the curves have different shapes. Hence, the obtained bound should not be parameterized like a usual PSD. One should bear in mind that every point of the continuous curve of the right-hand lower branch corresponds to its own discrete polydisperse system; hence, this curve cannot be considered as the cumulative PSD curve of any one polydisperse system. Besides, any curve drawn within the corridor according
k r ; Fr 1 − exp − R k k r k−1 r exp − : f r R R R The integration to calculate turbidity was carried out by the five-point Gauss–Laguerre quadratures [14]. Mie’s efficiency factor Kα was calculated by the algorithm from [19]. 3. Results and Discussion
Our preliminary calculations revealed that the minimum average radius was vanishing. This can be explained as follows. Let us take a model system, calculate the turbidity ratios T for it, and start to add particles with a negligibly small radius (e.g., the Raleigh ones). This action will render almost no effect on the values of T, since in the sums of the equation
Fig. 1. Cumulative PSD curves: source (thin) and calculated right-hand lower limiting ones (thick) for the set of wavelengths 400, 440, 490, 550, and 585 nm. Simulation: Rosin–Rammler’s functions, k 10, R 450 nm.
P 2 r K2πμ1 ri ∕λ1 wi T P i2 ri K2πμ1 ri ∕λ0 wi the corresponding r2i K2πμ1 ri ∕λ → 0, but the average radii will decrease to as low values as desired, almost to zero, because the weights wi of the “old” particles will decrease. We note in passing that such a situation is not provided for at the usage of common analytical forms of PSD functions; e.g., in unimodal forms, equally low amounts of very small and very big particles appear. 304
APPLIED OPTICS / Vol. 53, No. 2 / 10 January 2014
Fig. 2. Cumulative PSD curves: source (thin) and calculated right-hand lower limiting ones (thick) for the set of wavelengths 400, 440, 490, 550, and 585 nm. Simulation: a monodisperse system, R 450 nm.
to common restrictions (starting from zero, ending with unity, nondecrease) may not turn out to be the cumulative PSD of some system because different points on such a curve may well correspond to different systems (i.e., systems with different PSDs). If a monodisperse system was modeled, the convergence was fair as well (Fig. 2, 450 nm).
5. 6.
7.
4. Conclusion
In this paper, an algorithm is offered to estimate the right-hand lower bound of cumulative PSD functions consistent with some given turbidity data. The branches of the corridor thus serve as envelopes to the full set of cumulative PSDs consistent with experimental data. Sometimes (see [3]) turbidimetry cannot provide the full PSD but only an average of the PSD. In such cases our algorithm will be helpful enough. We have not considered the effect of experimental errors; they, if included, will cause T’s to vary and expand the calculated range rightward. The algorithm developed can be useful for selecting the number of wavelengths and their values in a certain application. However, even taking an infinite number of wavelengths will not shrink the bounds to the source curve.
8.
9. 10. 11. 12. 13. 14. 15.
References 1. S. Twomey, Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements (Dover, 1977). 2. S. Y. Shchyogolev, “Inverse problems of spectroturbidimetry of biological disperse systems: an overview,” J. Biomed. Opt. 4, 490–503 (1999). 3. T. Kourti, “Turbidimetry in particle size analysis,” in Encyclopedia of Analytical Chemistry: Instrumentation and Applications, R. A. Meyers, ed. (Wiley, 2000), pp. 5549–5580. 4. H. Tang and S.-Y. Ding, “Inversion of particle size distribution with unknown refractive index from specific turbidimetry,” in 2010 International Conference on Computer,
16. 17. 18. 19.
Mechatronics, Control and Electronic Engineering, Vol. 3 (IEEE, 2010), pp. 116–119. G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. 330, 377–445 (1908). B. Sedláček, “Light scattering. XVI. Particle size determination by the turbidity ratio method,” Collection of Czechoslovak Chemical Communications (1929-1939, 1947-) 32, 1374–1389 (1967). B. Sedláček and K. Zimmerman, “Integral and differential turbidity ratio method applied to the analysis of changes in the size and number of scatterers,” Polymer Bulletin 7, 531–538 (1982). T. Kourti, J. F. MacGregor, and A. E. Hamielec, “Turbidimetric techniques. Capability to provide the full particle size distribution,” in Particle Size Distribution II, T. Provder, ed. (ACS, 1991), pp. 2–19. G. B. Dantzig and M. N. Thapa, Linear Programming, 1: Introduction (Springer, 1997). O. Platz and R. G. Gordon, “Rigorous bounds for timedependent correlation functions,” Phys. Rev. Lett. 30, 264–267 (1973). R. McGraw, “Sparse aerosol models beyond the quadrature method of moments,” AIP Conf. Proc. 1527, 651–654 (2013). B. D. Bunday, Basic Optimization Methods (Edward Arnold, 1984). R. L. Zollars, “Turbidimetric method for on-line determination of latex particle number and particle size distribution,” J. Colloid Interface Sci. 74, 163–172 (1980). G. A. Korn and T. M. Korn, Mathematical Handbook for Scientists and Engineers (Dover, 1961). W. Heller, H. L. Bhatnagar, and M. Nakagaki, “Theoretical investigations on the light scattering of spheres. XIII. The “wavelength exponent” of differential turbidity spectra,” J. Chem. Phys. 36, 1163–1170 (1962). D. H. Melik and H. S. Fogler, “Turbidimetric determination of particle size distributions of colloidal systems,” J. Colloid Interface Sci. 92, 161–180 (1983). W. Heller and W. J. Pangonis, “Theoretical investigations on the light scattering of colloidal spheres. I. The specific turbidity,” J. Chem. Phys. 26, 498–506 (1957). P. Rosin and E. Rammler, “The laws governing the fineness of powdered coal,” J. Inst. Fuel 7, 29–36 (1933). C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley-Interscience, 1983).
10 January 2014 / Vol. 53, No. 2 / APPLIED OPTICS
305
