Multiple scattering calculations for technology II

Multiple Scattering Calculations for Technology II P. S. M U D G E T T AND L. W I L L A R D R I C H A R D S

National Research Corporation, (A Subsidiary of Cabot Corporation), 70 Memorial Drive, Cambridge, Massachusetts 02142 Received November 1, 1971; accepted January 31, 1972 Equations are developed for calculating the reflectance and transmittance of media which scatter and absorb light and are bounded by parallel planes. A comparison between the differential equations of the many-flux and the two-flux (Kubelka-Munk) calculation methods shows that the two-flux scattering and absorption coefficients are usually different for the forward and reverse flux, and gives values for these coefficients which make the two-flux equations give the same results as the many-flux equations. This comparison can be made in the limit of an infinite number of manyflux channels, and it provides a straightforwaid theoretical derivation of the two-flux scattering and absorption coefficients. The comparison is extended to obtain accurate coefficients for the four-flux, six-flux, etc., calculation methods. It is shown that the six-flux calculation is particularly useful when the internal reflection at the boundaries of the scattering medium depends strongly on the angular distribution of the radiation. INTRODUCTION Quite elaborate calculations are required to obtain accurate results for the intensity, polarization, and angular distribution of the radiation in a medium which absorbs and scatters light. However, many problems in technology require only part of this information and put a great premium on calculation methods which are simple. Therefore approximate methods, such as the two-flux calculation, are widely used. The purposes of this paper are to: (a) explore the properties of the two-flux calculation further to see if its accuracy can be improved without spoiling its simplicity; and (b) devise other calculation methods which are nearly as simple but much more reliable. Providing calculation methods with a range of compromises between simplicity and accuracy makes it possible to do only as much work as necessary for the particular problem at hand. Progress is made in this paper in developing an understanding of the two-flux calcuCTEyright (~) 1972 b y Academic Press, Inc.

lation and in determining the information which must be known to do accurate calculations with it. As a corollary, these results give information about the things which can be learned about a scattering medium from a two-fux calculation. It turns out that the two-flux calculation can be used quite reliably when there is little or no absorption, but that the more complete calculation methods described here should be used if good results are desired for absorbing media. TWO-FLUX CALCULATIONS The two-flux calculation was first introduced by Schuster (1), and since then many papers have been written to describe elaborations on the original equations and ways in which their accuracy can be tested (2). The most widely known relations derived from these equations are due to Kubelka and Munk. The two-flux calculation procedure is restricted to a homogenous scattering medium bounded by parallel planes. In addition, the medium should extend for distances Journal of Colloid and Interface Science, Vo[. 39, No. 3, J u n e 1972

551

552

MUDGETT AND RICHARDS

large compared to its thickness and be uniformly illuminated so that M1 of the calculated quantities are independent of the coordinates parallel to the boundary planes. The simplicity of the calculation procedure arises from the fact that only the light fluxes in the forward and reverse directions, F / a n d F,,, are determined. These quantities give the rate of energy transfer across a unit area parallel to the boundary planes. The differential equations on which this calculation method is based use the coefficient K to describe the weakening of each flux due to absorption, and the coefficient S to describe the scattering of light from the direction of one flux into the direction of the other. The equations are dFffdx = -- (K --~ S ) F / + SFr -dF,/dx

= SF s -

(K Jr S ) F , ,

[11

and their solution is easily obtained and fit to any of a variety of boundary conditions (3). We have not seen a statement that K and S must be different for the forward and reverse fluxes, but so many papers have been written o n tw0-flux calculations that such a statement undoubtedly exists. It is shown here t h a t K and S are almost always different for the forward and reverse fluxes in situations of practical interest. Therefore, we shall write dFffdx = - (KI + Ss)F/ A- S~F~

[2]

- d E ~ / d x = SsFs -- (K~ Jr S~)F~.

Scattering Coe~cients

The usefulness of these equations depends on the ability to obtain good scattering and absorption coefficients, so this part of the problem occupies the major portion of the remainder of the paper. Absorption Coe~cients

If a collimated beam of light of flux F passes through a layer of infinitesimal thickness dx and absorption coefficient ]c, the change in flux due to absorption is dF = - k F dx/I cos 0 I,

where 0 is the angle between the beam of light and the normal to the layer, and dx/l cos 01 is the pathlength of the beam in the layer. If light is incident over a range of angles, it is possible to divide the flux into many collimated beams, use the above equation for each beam, and sum the results to obtain the total amount of radiation absorbed. As is customary in the two-flux calculation, the absorption is described by the coefficient K, which is defined in such a way that the change in the diffuse flux F due to absorption is given by - - K F dx. It is apparent that the value of K depends on the angular distribution of the light, and that it becomes larger as more of the light travels in directions far from the normal. The relation between the angular distribution of the flux and K has been derived by Kubelka (4), and it is shown that K = 21~ when the radiation is perfectly diffuse. This ease is not often encountered in practice. For example, if there is a net transport of flux in the forward direction, then it is usual for the forward flux to be more intense near the normal and the reverse flux to be less intense near the normal. As a consequence, it is usually found in this ease that Ks < 2k < K~. If the direction of the net transport is reversed, then these inequalities are reversed. Formulas for the numerieal evaluation of Ks and K, are given later in the paper in Eqs. [35] and [36].

[3]

Journal of Colloid and Interface Science, VoI. 39, No. 3, J u n e 1972

The scattering coefficients are more difficult to deal with because it is necessary to determine where the light goes as well as how much light is scattered. The calculation of the amount of light scattered into all directions is parallel to the calculation of the amount of absorption. Thus the weakening of a collimated beam of flux F on passing through an infinitesimal layer is dF = - - s F dx/] cos 01,

[4]

where the scattering coefficient s describes

MULTIPLE

SCATTEP~ING

the extinction due to scattering and the other symbols are as in Eq. [3]. For the moment it is possible to simplify the problem of determining where the radiation goes by considering only Rayleigh scattering. Then for any ray, half of the light t h a t is scattered goes in the forward direction and half in the reverse direction. In this case the derivation of the two-flux scattering coeffieients Ss and & parallels the derivation of Ks and K~, except that only half of the light scattered from one flux goes into the range of angles included in the other flux. Thus for Rayleigh scattering Ss/s = ~Ks/£

and

S~/s = ~/~K~/£

[5]

for any angular distribution of the flux. For diffuse flux, this would mean that Ss ~ S~ 8.

Attempts to carry out two-flux calculations in cases where the radiation is very nearly diffuse clearly show that the S which must be used in Eqs. [1] to obtain good results is about 0.75 s (5). This apparent discrepancy has been the source of some confusion, and arises from the attempt to use one scattering coefficient S for the forward and reverse fluxes as in Eqs. [11 rather than different values for Ss and & as in gqs. [2]. One way of illustrating the source of this discrepancy is to imagine a situation in which everything is kgown about a multiple scattering problem except the values of Ss and & . Then Eqs. [2] form two equations with two unknowns. The attempt to solve them for Ss and & leads to only one constraint: & F s - S~F~ = -(~lFs/dx ) -- K s F s

[6J

= -- (dF~/dx) + K~F~.

This shows that for any value of Ss, there is a value for S~ which makes the two-flux equations give the correct values for F s and F~. However only one pair of coefficients will give the correct values to the amount of flux scattered from the forward direction into the reverse S j F s , and the amount scattered from the reverse into the forward direc-

IN

TECHNOLOGY

553

tion SrF~. If these "correct" coefficients are known, then any other pair of scattering coefficients Sf and & which satisfy the equation $sFs - &Fr = SsFs - SrFr

[7]

will give the proper values for the forward and reverse flux throughout the scattering medium. In particular, if it is desired to replace Ss and Sr by one scattering coefficient S, the value which should be used is S = (SsF s - S r F ~ ) / ( F s - Fr).

[8]

The ideas just introduced can be illustrated by some numerical values obtained from a many-flux calculation for Rayleigh scattering using 26 channels. In this example k / s = 1.0 X 10 5, the medium has a very great optical thickness, and the surface which is illuminated with 100 units of perfectly diffuse flux is smooth so that the Fresnel equations can be used to calculate the reflection at the surface. The calculations were done using surface reflection coefficients appropriate for an index of refraction ratio of 1.54 as reported in Table IV in (5). The results obtained from this many-flux calculation at a distance x from the illuminated surface such that sx = 10 are given in Table I. The value of S / s calculated from these data using Eq. [8] is 0.753. According to these results, the correct values for S s / s and S , / s are quite close to unity. However they do differ from each other slightly, and that difference makes it necessary to use a value for S / s which is quite different from unity. Thus, if only one scattering coefficient is to be used for the forward and reverse fluxes, it seems wise to think of it as a number which makes the TABLE Two-FLux

I

Q U A N T I T I E S FRO~I ~.

26-FLUX F s = 222.26 Fr = 220.65 dFs/dx = -1.21997 dF~/dx = - - 1 . 2 1 1 1 2

C~kLCULXTION KS Kr SS S~

= 1.99819 = 2.00179 = 0.99908 = 1.00088

k k s s

Journal of Colloid and Interface Science, VoL 39, No. 3, June 1972

554

MUDGETT

AND

two-flux calculation give good results rather than as a number which describes the total scattering between the forward and reverse directions. For a scattering coefficient to give good results, it is only necessary that it correctly describe the net scattering of flux between the forward and reverse directions. Figures 1-3 give additional examples of the scattering coefficients Ss and Sr which make the two-flux equations give the same values for F/ and Fr at every place in the medium as are obtained from a 26-flux calculation. Except as noted in Table II, these calculations were done exactly in the same way as those which gave the results in Table I. The coefficients shown in the figures are not constant, so the two-flux differential equations might be solved most conveniently by numerical integration. The figures also show the value of the single coefficient S which can be used in place of Ss and Sr without changing the values calculated for Fs and F,.. Figure 1 gives results from the same calculation as was used to obtain the data in

RICHARDS

1.3

I

I

i

1

i

I

(

I

I

I

I

I

I

I

Sr/S

/

s/~

J

1.2

1.1

I

1.0 i

>--%/s 0.9" -

i

0.8-

S o.6

0.5

I

I 2

o

I

I 4

I

I 6

I

I

I

8

zo

OPTICAL THICKNESS (s + k)x ~ sx

FIG. 2, T w o - f l u x scattering coefficients obtained from a 2 6 - f l u x c a l c u l a t i o n . ST = 1.89 a t s x = 9.99.

OPTICAL THICKNESS (s + k)x 8 12 16

4 I

1,3 i

I

I

I

E

I

I

I

I

20 I

l

1.01-

1.2 Sr/S

1.00

~2

S/s

1,1

sf/s

Y l

f

1.O

0.99

~ o

~

0.8

i

0.7

N

s/~

/

0.9"

! 0.8-

0.6

m

Sf/s

0.7

/ 0.5

I

I 2

I

I

4

1

I

6

I

!

8

I

io

0.6

F

s/s

OPTICAL THICKNESS (s 4-k)x ~ sx

FIG. 1. T w o - f l u x scattering coefficients obtained from a 2 6 - f l u x calculation as described in the text. N o t e that the scale is expanded by a factor of 10 near unity to show the difference between S / and S~ .

Journal of Colloid and Inter'~ace Science, Voh 39, No. 3, June 1972

O.5

t

I

i

I

I

I

I

OPTICAL THICKNESS (SCATTERING ONLY)

sx

2

4

I

6

I

8

lo

FIG. 3. T w o - f l u x scattering coefficients obtained from a 26-flux c a l c u l a t i o n . S r = 1.76 a t ~x = 9.99.

MULTIPLE TABLE

SCATTERING

II

VALUES PERTAINING TO FIGURES 1-3

Figure

k/s

sX

1 2 3

1.0 X 10 -~ 1.0 X 10-5 1.0

~ 10 10

Reflectance Total a t the back reflectance boundary

IN TECHNOLOGY

555

weakened by only 2.6 % on going from the front surface to sx = 5. Therefore, two-flux calculations are most useful when k/s is small. The Differential Equations

-0.0 0.0

0.983 0.616 0.144

Table I. It can be seen that the angular distribution of the forward flux, and hence the value of Sffs, changes over the first few units of optical thickness ; and then becomes constant through the rest of the medium. Figures 2 and 3 show results of calculations for media which have a thickness X such that sX = 10, and which are bounded at the side away from the illumination by a perfect absorber. Therefore, the reverse flux is zero and S becomes equal to Sz at x = X. In Fig. 2, the absorption due to the backing has an appreciable effect on the angular distribution of the flux and hence on Si and S,. It is interesting that the values of S are similar to those in Fig. 1. Figure 3 shows an example of the general result that as h/s increases, the angular distribution of the flux becomes more asymmetric, and the scattering coefficients take on quite different values from those appropriate for diffuse radiation• To help interpret these data, Table II also gives the reflectance calculated for each medium (including the specular reflectance at the front surface). The dependence of S/s and K / h for the forward and reverse fluxes upon k/s makes it difficult to use the two-flux calculations for media with appreciable absorption. This problem is compounded by the fact that the observed reflectance of an absorbing medium is primarily determined near the surface where the angular distribution of the radiation, and hence the scattering and absorption coefficients, are rapidly changing With depth. For example, in the case shown in Fig. 3 only 2.3 % of the radiation penetrates to a depth of sx = 2. By way of contrast, the forward flux in the case shown in Fig. 1 is

It is shown later in the paper that it is simpler to use the one scattering coefficient S rather than the two different ones, Sf and S,. When only one scattering coefficient but two different absorption coefficients are used for the forward and reverse fluxes, Eqs. [2] become dFs/dx = -- (Ks -}- S ) F / + SF, [9] - dFr/dx = SFs -- (Kr + S)Fr. The general solution to these equations can conveniently be expressed in the form

F~ = C+(1 - ~)e x+" + C_(1 + ~)e ~-x I10j F, = C+(1 + ~)e ~+~ + C_(1 - ~)e ~-',

where the eigenvalues h+ and h_ are given by 2h~ = (K~ - Ks) ~ [(Ks + K ~ ) • (K/+

[11]

K , + 4S)] "~

and

= [(Ks + K , ) / ( K f + K, + 4S)] "2.

[12]

The constants C+ and C_ are determined by fitting the general solution to the boundary conditions of the problem being solved. In the special case where the scattering material is infinitely thick, then C+ = 0 and the reflectance R~ is equal to (1 - 3 ) / ( 1 + fl). This relation can be rearranged to the more familiar form 2R~/(1 -- R=) 2 = 2 S / ( K / q- K,).

[13]

It is thought that these expressions will receive only limited use. If the angular distribution of the light is far enough from being diffuse that it is important to take into account the fact that K~ ~ K:., then it is likely that the angular distribution of the radiation is poorly enough known that a calculation with more than two channels Journal of Colloid and Interface Science, Vol. 39, No. 3, June 1972

556

MUDGETT AND I~ICHARDS

would be preferred. The six-flux calculations described later in this paper are successful enough that no calculations have been done to verify this opinion. MANY-FLUX CALCULATIONS

is the angle between channel i and the positive direction normal to the planes bounding the scattering medium, and the az are coefficients which express the phase function (scattering function) as a series of Legendre polynomials ~ = 0 azPz(cos 0). In all cases in this paper it is assumed that the phase function is normalized so that a0 = 1, that an even number of channels is used, and that channels i and n + 1 - i have the same solid angle. These restrictions on the channel divisions introduce a useful symmetry about the forward and reverse directions. The general solution to the many-flux equations gives the flux in all channels at any distance x into the scattering medimn

The many-flux calculation can be regarded as an extension of the two-flux method for solving multiple scattering problems fer uniformly illuminated media that are bounded :by parallel planes. This section of the paper extends the ideas introduced in the preceding discussion to this more general calculation method. Since a rather complete discussion of the many-flux equations has been published by Mudgett and Riehards (5), the description of the equations here is quite brief. However, there is one change in Fi = ~ CjAije xjx i =- 1, 2, " . , n, [19] j=l notation. Here the distance x is measured in units of length, as in the two-flux and four- where the Xj are the eigenvMues and the Ai~. flux equations in (5) and much of the tech- are the corresponding eigenvectors of the nological literature. Therefore the many- matrix "yiS~j, and the constants Cj. are obflux scattering and absorption coefficients t a i n e d by requiring that the solution i~t the are a factor of s larger here than in (5). boundary conditions. The differential equations for the manyThere are numerous ways in which the flux calculation with n channels can be scattering coefficients in the many-flux equawritten as tions can be changed without influencing the calculated results. The procedure chosen "Yi(dF~/dx) = ~ S~jFj in this work can be introduced by combining 3"=1 [14] Eqs. [14] and [17] to obtain i = 1,2, . . . , n , "y~(dFi/dx) =

where ~,~ = 1

for

i <= n/2,

~,~ = --1

for

i>

~

j=l,j#i

S~jF~ - K i F i

[15]

So"

s~i ~ at P~(cos Oi) 4v I cos 03. I l=0 Pl(COS03.)

--

n/2,

for

=

--Ki --

~

m=l,m#i

SmYi ,

which can be written as [16]

"Yi (dFi/dx )

[21]

i#j = -KiF~ +

Sii

~ m=l,m~i

L _

[20]

n

Smi,

[17]

~

j~l,jT~i

(S~jFj-- SjiF~).

The quantity Si3.Fj -- S/~FI gives the net flux transfer from channel j to channel i due to scattering. As long as the value of this Ki = k/I cos Oi I" [18] quantity is preserved for all i and j, both the The s and k are the extinction coefficients net flux transfer between all combinations of for scattering and absorption as in Eqs. [3] channels and also the rate of change of the and [4], ~ is the solid angle in channel i, 0~ flux in each channel is not changed. Thereand

Journal of Colloid and Interface Science, Vol, 39, No. 3, J u n e 1972

M U L T I P L E S C A T T E R I N G IN T E C H N O L O G Y

fore it is appealing to use this as a constraint on the scattering coefficients if they are changed from their correct values. It is apparent that even with this constraint there is still considerable freedom, in the values that may be given to the coefficients. It often simplifies the calculations if the scattering coefficients that are symmetrical to each other with respect to the forward and reverse directions have the same value. Therefore, new scattering coefficients are introduced which have the property that ,S~,. = S.~s

and

[22] S~e = Ssi,

where the subscripts are related as follows: q=n+l--i

[23] r =n+l--j. The~eonstraint that the new coefficients give the correct net flux transfer between channels, and hence the correct calculated flux, leads to the equations SoF.i

-

Sj~F~ = goF,

-

-

$iiFi

-

and

[24] SqrFr

-

SrqFq = 8qrFr

-

-

grqFq.

-

Here $ is used to designate the correct scattering coefficients which are being replaced. The new scattering coefficients are now uniquely specified, and given by the equation &i = [F~ ( ~ F j -

-

SjIF~)

F~(~rF~ -

(FiFq

S~.~,F~)]/ -

-

[9~]

F~Fi)

for i ~ j. In the special c~se where j = n q1 i then q = j, r = i, and Eq. [25] simplifies to -

-

S~j = (&yFy - gj~F~)/ (Fy -- F~).

[26]

It is apparent that Eq. [8] is a special case of Eq. [25].

557

Before these relations can be used, it is necessary to know the correct values for the scattering and absorption coefficients for each channel. Equations for these quantities will now be derived. T H E G E N E R A L D E R I V A T I O N OF T H E SCATTERING AND ABSORPTION

COEFFICIENTS The approach used in the following mathematics is to establish a correspondence between a many-flUX calculation with a large number of channels and one with fewer channels. Tl~is makes it possibleto find the coefficients for the simpler calculation which make it give the same results as the one with the larger number of channels. Once these general results are obtained, it is possible to let the larger number of channels become very large, or even go to infinity. It is also possible to use any convenient number of channels in the smaller calculation. In the extreme case, this permits obtaining the same forward and reverse flux from a two-flux cMculation as from a many-flux calculation with an infinite number of channels, provided the necessary information for the angular distribution of the radiation is available. In this comparison the subscripts i and j are used to designate channels in the larger many-flux calculation, and Greek letters to designate the channels in the few-flux calculation. The notation }-~,~=~is used to indicate that the summation extends over all values of i corresponding to channels which fall within the channel e. For convenience, it is assumed that the channel divisions are arranged so that all few-flux channel boundaries fall on a many-flux channel boundary. The results of this comparison are simple enough that they are first asserted and then justified by a more careful, but somewhat involved, derivation. The total flux in an arbitrary few-flux channel e is obtained by adding the flux in the corresponding manyflux channels, F~ = ~ F ; .

[27]

Journal of Colloid and Interface Science, Vol. 39, No. 3, June 1972

558

MUDGETT AND RI~HARDS

Adding the flux absorbed in each of the many-flux channels contained within channel ~ gives ~.=~ K3F~ dx. If we wish to represent this absorption by the expression K~F, dx, then the coefficient K~ must have the value

K~ = ~ gjFj/~-~ F~.

[29]

The demonstration that these coefficients do indeed cause the two calculations to give the same results for the flux in each few-flux channel is straightforward but lengthy, so it is only outlined here. The rate of change of the flux in an arbitrary channel e can be found by adding the rates of change of flux in the corresponding many=flux channels. Then using Eq. [14] gives

"=

= ~

~

= - Xjx= c,

m=l,m#j

r

-X so, j=~ m=l,m#~

Again the terms in the double sum are to be grouped according to the few-flux channel divisions to prepare the expression for comparison with

--K,F,-

~

S,,F, = S,,F~.

[32]

a=l,a#~

By actually carrying out the comparison it can be shown that Eqs. [28] and [29] give the few-flux coefficients that cause the two calculation methods to give the same results for the flux in the few-flux channels. In order to use these equations it is necessary to have some information about the dependence of the flux on the polar angle 0. Therefore, it is convenient to introduce a general expression for the flux ..

F¢ = I #J J o J (~j),

[33]

where ~ = cos 0 and the function I ( ~ ) can be expanded as a series of Legendre polynomials

j~l

S~ F~

[31]

j=e

[28]

Similarly, the total flux scattered from the few-flux channel e into a different channel is found by adding the scattering for all combinations of many-flux channels to obtain ~--~4~ ~-~i=~ S~Fs dx. If this is to be represented by the expression S~,F, dx then,

= E E

= E~=~ ~=~,j~i E

[30]

M

I(~) = ~

c2,,,(,).

[34]

m=0

i=~ j = l , j # e

Since there is a few-flux channel division at 0 = ~r, the Vi will always have the same sign as ~ . The terms in the second double sum can be grouped according to the few-flux channel divisions to prepare it for comparison with ~--:~=~,~ S~F~, where n is the number of few-flux channels and the sum extends over all few-flux channels except channel e. The first double sum in Eq. [30] can be split into two parts and Eq. [17] used to obtain

E Es, F

If the radiation is completely diffuse, Co is the radiance and all other cm equal zero. With the aid of this expression and Eq. [18], Eq. [28] can be written as

K¢/k

=

[35]

In the limit as the solid angle ~j of each many-flux channel within the few-flux channel e goes to zero, ~, becomes 27r dm and the Sums can be expressed as integrals to obtain K.

I I ( g ) dg -

k i=~ j ~ e , j # i

j=~

Journal of Colloid and Interface Science, Vol. 39, N o . 3, J u n e 1972

"

.I I I

,

d.

[36]

MULTIPLE SCATTEP~ING IN TECHNOLOGY which is the general expression for the fewflux absorption coefficient. The formula derived by Kubelka (4) for the two-flux absorption coefficient is a special ease of this relation. Similar substitutions can be made in Eq. [29] for the scattering coefficient using Eqs. [16] and [33]. Here it is convenient to rearrange the order of summation to obtain

559

The reason for this is that the function [ (u) can be adequately approximated in advance of the calculations. A good example of this appears in the following section, where it is shown that the two-flux scattering coefficient is independent of the angular distribution of the radiation as long as this distribution can be expressed by only the first three terms of Eq. [34].

L

The Two-Flux S It turns out that Eqs. [8], [34], and [38] can be combined to obtain a simple, completely general expression for the two-flux scattering coefficient to be used for both the forward and reverse fluxes. Expressions for S s and S~ are obtained from Eq. [38] by making the following substitutions in the integration limits: For Ss, E is replaced by 0 t o l and v by - l to 0, a n d f o r & , e i s r e placed by - 1 to 0 and ~ by 0 to 1. Then Eq. [8] can be written as

4=e

for

u¢ e

which becomes L

S~. ~-'~=oa' .I I(")P'
S 8

, ~ e

1

Pz(~) d~ /=i

i m=O

P 1

M

d-- I

m~O

in the limit as the solid angle of each manyflux channel goes to zero. It is not necessary to retain the subscripts i and j on ~ because it is possible to completely separate the integrations. Equations [36] and [38] give the values of the scattering and absorption coefficients which cause the few-flux equations to give the same results for the flux in each fewflux channel as are given by the many-flux equations in the limit of an infinite number of channels. Therefore, it is recommended that Eq. [16] be replaced by Eq. [38] (with the appropriate replacement of ~ by i and e by j), and that the K~ be calculated from Eq. [36] (with e replaced by i). Equation [17] remains unchanged. The fact that the function I (~) in Eq. [34] must be known to use these relations is not as much of a problem as might at first appear.

[39]

where use has been made of p 0

1

J_, p,(.)

---

-./o- ~ P,(.) d.

[40]

for 1 # 0

and integrals combined where possible to give an integration over the range from - 1 to 1. Since P1 (~) --- t* and I

L P,,~(u)Pt(I~) dl* = 0 1

-

for 2

m¢ l [41] for m = l,

2m+l

Eq. [39] can be simplified to S 8

3 i 251 m=l

c,,[ao -- am/(2m + 1)] 1

[421

• £ P~(.) d~, which is the general expression for the twoJournal of Colloid and Interface Science, V o l . 39, N o . 3, J u n e 1972

560

MUDGETT AND RICtIARDS

flux S. The remaining integral has the values

THE SIX-FLUX CALCULATION

The preceding derivations have shown that a many-flux calculation can be done with any even number of channels and the = 0 for m = 2, 4, 6 , . . . same results obtained as would be obtained = (1-3.5 .... [43] in the limit of an infinite number of channels, provided sufficient information is available m ( m -~- 1)m! about the angular distribution of the radiafor m = 1 , 3 , 5 , . - . tion. As the number of channels is increased, and a0 = 1. Therefore, the first few terms in the scattering and absorption coefficients for each channel become less sensitive to changes Eq. [42] are in the angular distribution of radiation. In addition, more information about the angular s - 4 ( 1 - ax/3) - ~ 1 1-- distribution of the radiation is generated by the calculation. Here it is shown that a sixc0 ( a s ) 5c7 ( a T ) [44] flux calculation represents an excellent com1- H 1-T5 promise between simplicity and reliability. The equations are simple enough that the eigenvalues can be found analytically, and yet there are enough channels to allow the It is interesting to note that if the angutar use of scattering and absorption coefficients distribution of the intensity can be described which are independent of the angular distriby only the first three Legendre polynomials bution of the radiation in a rather wide so that c~ = 0 for m > 3, Eq. [44] reduces to range of practical problems. In addition, the S / s = 3~ (1 -- ax/3) [451 calculation generates useful information about the angular distribution of the radiaand the two-flux scattering coefficient de- tion. pends only on the second term of the series It has long been known that the twoexpansion of the phase function. flux calculation is not adequate when colFor optically thick media in which k / s << 1, limated illumination is used, and this probthe angular distribution of the flux is de- lem has been solved by the addition of scribed quite well by only the first two terms channels for the collimated flux (6). The of Eq. [34] when the optical depth is greater resulting four-flux calculation is still not than 3. This explains why the above expres- adequate when the scattering media have sion was obtained in the empirical correla- an optical thickness less than a few units tions described in the previous paper (5). and a larger index of refraction than the surWe believe that this is the first satisfacroundings. In this case, radiation in the tory theoretical derivation of an expression medium which approaches the boundaries at for the scattering coefficient to be used in the directions farther from the normal than the two-flux equations. This expression makes it possible to do accurate, absolute calculations critical angle is totally internally reflected, with the two-flux equations in a rather broad whereas only a small portion of the remainrange of practical problems, and clearly de- ing radiation is internally reflected. Therefines the limitations of the two-flux method. fore, variations in the angular distribution The discussion now turns to a compromise of the radiation approaching the boundabetween the simplicity of the two-flux equa- ries from the inside cause large variations in tions and the reliability of a calculation with the fraction of the flux that is internally reflected. Numerical data showing the iraa large number of channels. Journal of Colloidand InterfaceScience,Vol.39, No. 3,June1972

fo x

MULTIPLE SCATTERING IN TECHNOLOGY portance of this effect are given in Table VI of Ref. (5). This problem is solved in the six-flux calculation by placing a channel division at or near the critical angle. Then the calculation generates enough information about the angular distribution of the flux at the boundaries to obtain good values for the internal reflection. The division of the forward and reverse diffuse flux into two channels makes it possible to do calculations of sufficient accuracy for all but the most demanding purposes without a priori information for the angular distribution of radiation within the scattering medium. The calculation is reliable for tc/s values less than 5 and for optieM thicknesses for scattering greater than 2. It is sufficiently difficult that a small computer is required, but the coding of the programs is not difficult, and nine decimal digit arithmetic is adequate. The differential equations and the general solution for the problem are given by Eqs. [14], [15], and [19]. Since the analytical solution of the equations is simplified by making the scattering and absorption coefficients symmetrical with respect to the forward and reverse directions, the adjustment of the coefficients described in Eq. [25] is used here.

The Absorption Coe~cients The absorption coefficients for channels of large solid angle depend on the angular distribution of the radiation in the manner described by Eq. [36]. When the radiation is perfectly diffuse

K~ = 1~/1 m 1,

[46]

where #~ is half way between the values of ~ = cos 0 at the boundaries of channel i. This coefficient is symmetric with respect to the forward and reverse directions. When the radiation is not diffuse, the absorption coefficients derived from Eq. [361 will not in general have the desired symmetry. If the coefficients are required to be symmetrical and the sum of the flux ab-

561

sorbed in channels i and q is correctly described, then

Ki = Kq = (XiFi + XqG)/(F~ + G),

[47]

where ~ and ;Eq are the "COlTeet" coefficients from Eq. [36] which are being replaced, and q = n + 1 - i. If the angular distribution of the radiation is such that the cm in Eq. [34] are nonzero only for m = 0 and for odd values of m, then by using the relations

~ Pl(,) dtt = (-1)~ f P~(~) d~ and

[48]

it can be shown that Eq. [47] gives the same coefficients as Eq. [46]. Therefore Eq. [46J, which is the same as Eq. [18J, was used in the six-flux calculations in this paper. This approximation is slightly better than the one used in Eq. [51] in the following derivation of the scattering coefficients.

The Scattering Coe~cients The calculation is simplified if the solid angles for channels 1 and 6 which contain the collimated flux are small enough that scattering into them can be neglected. Here it is fm'ther assumed that the collimated illumination is normal to the boundaries of the scattering medium. As a result of these restrictions, 812 =

S13 =

$14 =

S15 =

S16 =

$61

[49]

and Eq. [38] becomes L

S~=G~=~

= at

P~(~)d.

[50]

with q = n q- 1 -- i and i ~ 1. The scattering coefficients in Eq. [50] are symmetric with respect to the forward and reverse direction, so the use of Eq. [25] causes no change in them. In the derivation of the remaining scatterJournal of Colloid and Interface Science, Vol. 39, No. 3, June 1972

562

MUDGETT AND RICHAI~DS

ing coefficients it angular distribution scribed by only the [34], so this equation I

(~)

=

is proposed that the of the radiation be defirst two terms of Eq. becomes [51]

CO -}- el ].t.

is the unit matrix. Since the Sij are symmetric with respect to the forward and reverse directions, the odd powers of X do not appear when the characteristic equation is expanded to give

Then Eqs. [25], [33], and [38] give

f

s Ji l u ]du ~ al ~ij

( S ~ , - X2)(X4+ VX2+ W) = 0, [56] uPI(,) du

r l ( , ) d, -

uPl(,) d,

Pl<,) d~

"~

[521

where use has been made of Eqs. [48J and the following relation, which uses the notation of Eqs. [23]:

where V=

-

$2

2

2

22- Saa+Sa4+825

2

+ 2&s&4 - 2&a&2, f ~Pl(u) du fiPl(~) d~ [53]

= --£~'Pl(u) d u f Pl(Adu

[57]

and

w = l sl = 2 { & & (&2&~ + 835823)

for all 1. It should be noted that the scattering coefficients are independent of the values of Co and cl, and hence independent of the angular distribution of the radiation as long as Eq. [51] is satisfied. The integrals involving the Pz in Eq. [52] can be evaluated by using the expressions

&2&3(&s&4 + &3&2) + &s&3 (&a&~ + & ~ & ) -

2 2 2 2 2 2 + &2S33 - Sa4822 - &sSa3

72 Pl(,) du - -- 1 l+1

[gPl(~) - Pl-l(g)]~2~

[54]

2

2

2

2

2

2

2

2

-

832&4 + Sa4&5.

--X6 =

Sii

)M ---- --X4 = [ ( - V + I{

_

122

F~

- V / Vz -

%,/V~ -

4 W ) / 2 ] ~/2

4W)/2] 1/2.

l

+ (t-

1+ 1

[59]

1)(2z + 1)

[55]

} p,(~)

(l + 2 ) ( 2 / + 1) 21 + 1

2

,.-%5&3 + Sa2&a + 835&4

X2 = --X~ = [ ( - - r

=

2

-

The eigenvalues may then be expressed by Xl =

for 1 ~ 0 and

[581

- &&~ (&3&2 + &4&5)}

+

The eigenveetors corresponding to each X~are determined by the equations: 6

1

~Pl-l(~) ,,

~_, "r~SijAj, -

The eigenvalues are the roots of the characteristic equation / 7S - E x I = 0, where E J o u r nal of Colloid and Interface Science, VoL 39, No. 3, J u n e 1972

[60] i - - 1,2, - . . , 6 .

f o r / _-> 2.

The Eigenvalues and Eigenvectors

X~Ai~ = O,

j 21 Although due to a S matrix, Eq. [60]

some simplification is introduced number of zeros appearing in the the determination of the Aj-k from is sufficiently difficult to require

563

M U L T I P L E S C A T T E R I N G IN T E C H N O L O G Y

flux equations to any of a variety of boundary conditions is the same as for the manyflux equations (5). Therefore a discussion of this part of the calculation need not be repeated here.

numerical solution. The work in this paper used the subroutine SIMQ from the IBM 1130 Scientific Subroutine Package (7). Each eigenvector can be multiplied by an arbitrary constant and still satisfy Eq. [60], and it is this freedom which is used later in the calculation to fit the general solution to the boundary conditions. In order to assign definite values to the eigenvectors, we arbitrarily set All, A~G, and A2a for k = 2, 3, 4, and 5 equal to unity. A check for a satisfactory solution is obtained by evaluating the left hand side of Eq. [60] with i = 2. These checks are usually less than 10-5 with k/s in the range from 10-6 to 1, but become much larger with increasing k/s.

RESULTS

The results of calculations are presented here which demonstrate the reliability and usefulness of the six-flux calculation over an extended range of conditions. These comparisons and those previously reported for the four-flux calculation are primarily useful as an aid in selecting a computation method that gives the required accuracy with the least effort. The effect of increasing the number of channels in a many-flux calculation is shown in Table IIi. The exact values of the reflectance for the three phase functions were

The Cj The procedure for finding the values of Ci which fit the general solution of the six-

TABLE III A COMPARISON WITH I~,EFLECTANCES CALCULATED BY ORCHARD Collimated illumination Optical thickness

Exact reflectance

4-Flux

Diffuse illumination

Error X 104a Exact 6-Flux 22-Fluxb 22-Flux" reflectance

4-Flux

Error X 104 6-Flux 22-Fluxb 22-Fluxc

p(cos ¢) = 1.0 + 0.5P~(cos ¢) 1 2 4 6 8 10

0.3403 0.5162 0.6897 0.7729 0.8211 0.8524

--21 --27 --11 --1 3 5

1 2 4 6 8 10

0.2346 0.3990 0.5886 0.6892 0.7504 0.7916

--66 --77 --37 --16 --4 +1

1 2 4 6 8 10

0.1678 0.3163 0.5085 0.6189 0.6891 0.7375

--100 --113 --67 --33 --16 --7

--7 --5 1 3 4 4

--1 --1 1 0 0 0

0 0 0 0 0 0

0.4468 0.6106 0.7541 0.8204 0.8585 0.8833

--182 --101 --41 --22 --14 --10

--63 --33 --14 --8 --5 --3

--5 --2 --1 --1 --1 --1

14 12 9 6 5 4

--35 --51 --24 --14 --8 --7

--6 --4 --1 --1 0 --1

18 18 15 13 1i 9

--101 --66 --34 --20 --14 --10

--7 --5 --2 --1 --1 --1

22 23 22 19 16 14

p(cos ~) = 1 . 0 + 1.0Pi(eos ¢) + 0 . 5 P 2 ( c o s @) --23 --19 --6 --1 +3 +3

--3 --1 --1 0 1 0

1 3 4 4 5 4

0.3580 0.5157 0.6739 0.7541 0.8026 0.8352

--247 --157 --73 --41 --26 --19

p(eos @) = 1 . 0 + 1.5Pi(eos ~) + 0.5P~(eos ~) --35 --32 --15 --5 --1 2

--3 --3 --2 -- 1 0 0

4 7 9 10 9 9

0.3020 0.4490 0.6104 0.6985 0.7541 0.7924

--293 --204 --104 --62 --41 --29

T h e error is t h e m a n y - f l u x m i n u s t h e exact value. b M e t h o d of t h i s paper. ° F r o m tZef. 5. Journal of Colloid and Interface Science, Vol. 39, No. 3, June 1972

564

MUDGETT AND RICEARDS

calculated by Orchard (8). In the exact calculations there is no surface reflection at the boundaries and no absorption. The manyflux calculations were done with k/s = 10.6 to prevent degeneracy in the eigenvalues, but this small amount of absorption has no effect on the significant figures shown in the table. A comparison of the 22-flux calculation with the data of Orchard has been reported by Mudgett and Riehards (5). The calculations are repeated here with the S matrix obtained from Eq. [52]i rather than Eq. [16], and it can be seen that this change reduces the errors by about a factor of three. In all four-flux and six-flux calculations in this paper the channels for the collimated flux have a small enough solid angle that scattering into them is neglected. In the six-flux calculation the boundary between channels 2 and 3 is taken at cos 0 = V'2/2, so that when the illumination is diffuse the incident flux is divided equally between these channels. Since there is no flux in the channels for the collimated radiation when the illumination is diffuse, the four-flux calculation reduces to a two-flux calculation in this ease, and the errors reported in Table I I I for diffuse illumination also describe the

accuracy of the two-flux calculation. The use of the two-flux calculation for collimated illumination leads to even greater errors, and is not recommended. I t can be seen from Table I I I that the errors for the six-flux calculation are about

a factor of three smaller than for the fourflux calculation. This is some improvement, but it does not entirely justify the present discussion of the six-flux method. The motivation for introducing the six-flux equations arises primarily from the consideration of cases in which the index of refraction of the scattering medium is larger than that of its surroundings, so that some of the radiation is totally internally reflected at the boundaries. The following comparisons show the strength of the six-flux calculation in dealing with this problem. When comparing these results with the two-flux and four-flux results reported previously (5), it should be remembered that the surface reflection coefficients there were obtained from a manyflux calculation with a large number of channels. Figures 4 and 5 give a comparison between six-flux and 26-flux calculations. Here the scattering medium has an index of re1,01 T

1.01-

5O

D

x :::) hO0i,

5 -.u~&

1.00 -

$

0.99 -

"~ 0.99X

d ,.b 0.9s-

i

o

0.98 .

i-

n.- O. 97-

0,97

z ~ ~t t-

0.96

0,96

03

z,< 0 , 9 5 n-" 10.94

0,95

COLLIMATED I LLUM INATION I I 10 -6

I I I0 -4

I I 10-z

I I

I

I0 I 2

k/s

O. 9 4

u

DI FFUSE I LLUMINATION I I I0 -6

I I 10~4

I l 10-2

I

I

[

t 102

k/s

FIG. 4. The ratio of t r a n s m i t t a n c e s calculated from the 6-flux and 26-flux equations as a f u n c t i o n of

k/s. The numbers beside each curve give sX, the optical thickness for s c a t t e r i n g only. Journal of Colloid and Interface Science,

VoI. 39, N o . 3, J u n e 1972

MULTIPLE SCATTERING IN TECHNOLOGY

i. Oi

-

1.00-

,~ 0550:

~

I. O0

/ \

I

o_ x I- ~
0,99-

004 ~X

0,98-

Od Ldb_

2

•

0.99

S

0,98

I

0.97

t~ ~) 0.97" [K

0.98

0.96 • COLLIMATED ILLUMINATION BLACK SUBSTRATE 0,95

• ] ] 10-6

]

I

I I I0 "2

I0 -4

I I

I

] ]0 z

1.02.

o_

0.95

t

J DIFFUSE BLACK

I

I0 "6

I

i

I0 -4

ILLUMINATION SUBSTRATE I

I

-

I

I

-4 102

I

I IO z

I

i'%.

1.01

w~

•

1,00.

I

IO "2

1,02

x,0I

Z~l

565

-~

I,O0

~

5

a

U.ld Ld~

0.89'

0.98

/" COLLIMATED ILLUMINATION WHITE SUBSTRATE l I I I I I I I I0 "6 I0 -4 I0 "2 I

O, 9 9

I 102

k/s

0,98

DIFFUSE ILLUMINATION WHITE SUBSTRATE [ I [ I ] I I 10 "6 10"4 I0 -z 1

k/s

FIG. 5. The ratio of diffuse reflectances calculated from the 6-flux and 26-flux equations as a function of k/s. The media are in optical contact w i t h a substrate w i t h a reflectance of zero (black) or a specular reflectance of 0.94 (white). The numbers beside each curve give sX.

fraction of 1.54 and the angular dependence of the single scattering is typical of a rutile titanium dioxide pigment at a wavelength of 546 rim. The phase function and data for the 26-flux channel divisions and surface reflection coefficients are given in Tables I I I and IV of Ref. (5). In the six-flux calculation channels 2 and 3 are divided at the critical angle of 40.5 °, and the Fresnel reflection coefficients are 0.0452, 0.09763, and 1.0 for channels 1, 2, and 3, respectively. When the illumination is diffuse, the fraction 1 - 0.09763 = 0.90237 enters channel 2 and the rest is specularly reflected. The values for the reflectance in Fig. 5 do not include the flux that is reflected from the illuminated surface. The transmittances in Fig. 4 are calculated for free-standing media so that the surface reflections occur both at the front and back boundaries.

Figure 4 shows that the six-flux equations give transmittances which are reliable to within 1% of the value of the transmittance as long as l¢/s is less than 10-2. As the absorption increases, the transmittance becomes small and the relative error in the 6-flux calculations increases. Figure 5 shows that the calculated reflectances are nearly as reliable. In all reflectance and transmittance calculations here the magnitude of the difference between the 6-flux and 26-flux results is less than 0.006. Comparisons with experimental data are somewhat less informative because of a lack of accurate data which cover a range of optical thicknesses and ratios of absorption to scattering. In the past the data of Atkins (9) for a series of plastic discs containing titanium dioxide and a red dye have been found useful, and they are used again here Journal of Colloid and Interface Science, Vol. 39, No. 3, J u n e 1972

566

MUDGETT AND RICHARDS

even though it is believed that the experimental data are less reliable than the results of the six-flux calculation. These data have been discussed by Atkins and Billrneyer (10), and it has been shown that they can be fit to within the experimental accuracy by a 26-flux calculation (11). In this fit, only one parameter was adjusted, the relation between the amount of titania present and the coefficient s which describes the extinction due to scattering. In the present six-flux calculations no additional parameter was adjusted. The boundary between channels 2 and 3 was set at cos 0 = ~ / 5 / 3 , which corresponds to the critical angle for an index of refraction of 1.50. The surface reflection coefficients are 0.040, 0.091786, and 1.0 for channels 1, 2, and 3, respectively. The absorption coefficient of the dye was measured by Atkins, and the scattering properties of the titania were obtained from the Mie calculations and fit mentioned above. These data are given in (5) and (11). The six-flux calculations fit the experimental data about as well as the 26-flux calculations. This is shown by the results in Table IV, where the largest difference beTABLE IV A COMPARISON BETWEEN

MANY-FLux

CALCULATIONS AND EXPERIMENTAL D A T Aa

4-Fluxc 4-Flux~ 6-Flux 26-Flux °

Error of largest magnitudeb

SE

Reflec- Transtanee mittance

Reflec- Transtanee mittance

0.0540 0.0201 0.0166 0.0141

0.0261 0.0071 0.0065 0.0064

--0.0508 --0.0229 --0.0181 0.0213

0.0325 0.0147 0.0110 0.0104

The experimental data are from Ref. (9). b The error is the calculated minus the experimental value. o Constant surface reflection coeffÉcients were used. From (5). The surface reflection coefficients were obtained from the 26-flux calculation. From (11). Journal of Colloid and Interface Science, Vol. 39, No. 3, June 1972

tween experiment and theory is given along with the standard error (which is the rootmean-square difference between experiment and theory). I t was demonstrated in the previous paper that the four-flux calculations could fit the data provided good surface reflection coefficients were known beforehand. Here it is shown that the accuracy of the four-flux calculation suffers if these coefficients are not known, and constant values appropriate for diffuse illumination are used. 1 Comparing these results, which were all obtained from the same absorption coefficient, scattering coefficient, and phase function, clearly demonstrates that a major strength of the six-flux calculation is its ability to determine the fraction of the flux which is internally reflected at the boundaries. SUMMARY All of the multiple scattering calculation methods discussed in this paper fit into the general form of the many-flux equations, which are summarized beginning with Eq. [14]. I t is shown here that these equations give the best results when the scatterhag coefficients are chosen to give the correct net flux transfer due to scattering between channels rather than the correct total flux transfer. B y comparing calculations with differing numbers of channels, and allowing the number of channels in one of the calculations to go to infinity, it is possible to derive scattering coefficients which make the calculation with the smaller number of channels give the same results as would be obtained from the use of an infinite number of channels. Among other things, this comparison provides a theoretical derivation of the twoflux (Kubelka-Munk) scattering and absorption coefficients which makes the connection between them and the scattering and absorbing properties of the medium 1 The internal reflection coefficients at both boundaries are 0.040 for the collimated flux and 0.5964 for the diffuse flux.

MULTIPLE SCATTERING IN TECHNOLOGY definite and unambiguous. Unfortunately, this relation depends on the angular distribution of the radiation in the scattering medium, and this limits the usefulness of the two-flux equations. When k / s << 1 (i.e., the single scattering albedo is nearly unity), the angular distribution of radiation in optically thiek scattering media is such that the two-flux scattering coefficient becomes independent of the angular distribution of the radiation in the medium and the two-flux equations ean be used quite reliably. However, when the extinction due to absorption becomes comparable to the extinction due to scattering, the scattering coefficient depends on the angular distribution of the radiation rather strongly and usually varies with depth in the scattering medium. Under these conditions it is not advisable to use the two-flux calculation if accurate results are desired. The six-flux calculation is singled out for special discussion in this paper because it is believed that it is an excellent solution to a long-standing multiple scattering problem. M a n y scattering and absorbing media, such as paints and plastics, have an index of refraction appreciably larger than their surroundings so that a large fraction of the radiation is internally reflected at the boundaries. The fraction of the radiation internally reflected is quite sensitive to its angular distribution, so this information must be generated by any successful multiple scatter-

567

ing calculation method. The six-flux calculation is complete enough to generate the necessary information, yet simple enough so that the eigenvalues can be found analyticMly and the use of a large and time-consuming computer subroutine to determine them numerieally is avoided. It is believed t h a t these calculation methods are simple and reliable enough that they should find numerous practical applications. REFERENCES 1. SCHUSTER, A., A~trophys. J. 2i, 1 (1905). Reprinted in "Selected Papers on the Transfer of Radiation," D. H. Menzel, Ed., Dover Publications, New York (1966). 2. For a review, KOR~ffM, G., "Reflectance Spectroscopy," Springer-Verlag, New York (1969). 3. HAMAKlrR, H. C., Philips Res. Rep. 2, 55 (1947). 4. KUBELKA, P., J . Opt. Soc. A m e r . 38, 448,

1067 (1948). 5. MUDGETT, P. S., ~'kNDRICHARD.S, L. W., Appl.

Opt. I0, 1485 (1971). (Part 1) 6. SILBERSTEIN, L., Phil. Mat. 4, 1291 (1927). 7. "1130 Scientific Subroutine Package (1130CM-02X) Programmer's Manual," Doc. No. H20-0252-3, IBM Teeh. Publ. Dept. White Plains, NY (1968). 8. ORCHARD,S. E., Astrophys. J. 149,665 (1967). 9. Atkins, J. T., "Absorption and Scattering of Light in Turbid Media," thesis, Univ. Delaware (1965). 10. ATKINS, J. W., AND BILLM~:~na, F. W. JR., Color Eng. 6, 40 (1968). 11. RICRAaDS, L. W., J. Paint Technol. 42, 276 (1970).

Journal of Colloid and Interface Science, Vol. 39, No. 3, June 1972