- Email: [email protected]
A backward Monte Carlo study of the multiple scattering of a polarized laser beam
J. Quanr. Specwosc. Radiar. Tramfeer Vol. 58. No. 2, pp. 171-192. 1997 0 1997Elsevier Science Ltd. All rights reserved Printed in Great Britain
Pergamon
PII: !300224073(!37)00034-4
0022-4073/97 $17.00+ 0.00
A BACKWARD MONTE CARLO STUDY OF THE MULTIPLE SCATTERING OF A POLARIZED LASER BEAM A. AMBIRAJAN?
and D. C. LOOK
Thermal Radiative Transfer Group, Mechanical and Aerospace Engineering and Engineering Mechanics Department, University of Missouri-Rolla, Rolla, MO 65401, U.S.A. (Received 15 November 1996)
Abstract-A backward Monte Carlo estimator is developed to describe the multiple scattering of a polarized, narrow light beam by a plane-parallel medium. The case of a right circularly polarized beam is analyzed in this paper. Results indicate that the diffuse light field is partially polarized even at significant optical radii from the incident light beam. The degree of polarization of the diffuse light field is dependent on the optical thickness of the medium and the size parameter of the scatterers. 0 1997 Elsevier Science Ltd
NOMENCLATURE D = geometric diameter of scatterer e[...] = I,Q,CJ, V = I, = I= I, = I,, = K= L= L= III = M = n, = p,I,pzI,s21,d~l=
P = P, = P, =
P Z r s,~,; x X, Y,Z
= = = = = =
estimator components of the Stokes vector incident intensity Stokes vector Stokes vector for the delta function loading zeroth-order Stokes vector identity matrix optical thickness of slab rotation matrix complex refractive index Stokes matrix refractive index of liquid carrier medium elements of the Mie scattering matrix degree of polarization Degree of linear polarization degree of circular polarization phase matrix scattering matrix optical radius = m optical coordinates also used to denote the size parameter (n,nD/i,) Geometric coordinates
Greek s~~mbols aJ?.y = variables of integration for higher scatter orders
y,,y: = 6(.r) = 0= 0 = u= 4 = $ = w=
limits to p for single scattered radiation dirac delta function polar angle included angle during scattering rotation angle azimuthal angle parameter used in the arguments of I,, single scattered albedo
Subscripts
i = incident quantity t = total quantity Superscripts
’ = first scatter point i = photon number T =
transpose
tTo whom all correspondence
should be addressed. 171
172
A. Ambirajan and D. C. Look 1. INTRODUCTION
The multiple scattering of light within a material is an area of considerable interest. In medical diagnostics, the analysis of scattered laser light from tissue samples is increasingly being used as a diagnostic tool.‘m4The use of lasers in surgery also requires an understanding of the way in which laser light is scattered and absorbed by tissue.4 Lasers are also being used for remote sensing in lidar systems.5 Multiple scattering of light has also been used in process applications to determine particle concentrations.6 An issue that has not received much attention in the literature is the effect of multiple scattering on the polarization of a narrow light beam such as a laser. When a linearly polarized light beam is incident on a scattering medium, the backscattered radiant intensity exhibits a highly anisotropic pattern7-12 (the “bow-tie”). Look* observed that this bow-tie had a cosine variation. When viewed through a crossed polarizer, a four-lobed pattern was also observed. Dogariu and Asakura” showed that the bow-tie was a consequence of the polarizing effect of single scattering, and was thus mainly observable for lower orders of scatter and small size parameters. Dogariu and Asakura” also used the diffusion approximation to evaluate the cross-polarized component of backscattered light from a plane-parallel medium and obtained good results. These analyses thus indicate that polarization effects are propagated to significant distances in a scattering medium. Mueller and Crosbie have presented a formulation for three-dimensional radiative transfer with polarization by using the integral transform method.‘j They also developed certain spatial symmetry relationships for the reflection and transmission matrices. In a subsequent paper, Mueller and Crosbie14 used this integral transform method and the symmetry relationships to obtain the functional form (i.e., of the bow-tie) of the spatially varying reflection and transmission matrices for the multiple scattering of a polarized laser beam from a plane-parallel scattering medium composed of scatterers with a plane of symmetry. Their method is valid for any incident state of polarization and any axisymmetric spatial loading. Their transform technique, though exact, is time consuming to implement. In an interesting paper, Zege and Chaikovskya’5 simplified the vector radiative transport equation by using the properties of the Mie phase matrix and obtained an approximate system of equations whose solution describes the multiple scattering of a polarized narrow light beam. The backward Monte Carlo method has been used in studying the propagation of polarized light in scattering atmospheres.‘6, ” Ambirajan and Look’* suggested an algorithm to deal with the multiple scattering of an unpolarized narrow light beam neglecting polarization effects. The problem to be tackled in this paper is that of a polarized searchlight beam incident on a plane-parallel medium composed of scattering particles. In reasonably dilute suspensions, analysis of the multiple scattering of light can be carried out with radiative transfer theory.19 For higher number densities of scatterers, the analytical wave theory” may be more appropriate. In the analysis that follows, it is assumed that the scattering is independent, no coherence effects are considered, and the light beam is spatially incoherent; hence, the radiative transfer theory is used. In the analysis to be presented in this paper a technique is used that was developed earlierI to analyze the multiple scattering of a narrow light beam including polarization effects. The numerical data presented are for the case of a right circularly polarized searchlight beam. 2. SINGLE
SCATTERING
OF POLARIZED
LIGHT
In this section, some basic theory used in this paper to describe the multiple scattering process will be presented. First, the description of polarized light in terms of Stokes vectors will be outlined. Subsequently, the interaction of polarized light with a scattering medium will be characterized in terms of a matrix. 2.1. Light field In the theory of radiative transfer by light that is polarized, the state of polarization and intensity of a beam is specified by the 4 x 1 Stokes vector (I) in the following form:
I
=
G ii U’ V
(1)
Multiple scattering of a polarized laser beam
173
The elements Q, U and V define the state of polarization of the light beam and I denotes the intensity.*’ Another quantity of interest in radiative transfer of polarized light is the degree of polarization (P) of the light. This is defined as
*+ u* + v*
.JQ
I
and indicates the fraction of the light intensity that is polarized regardless of whether it is linear and/or circular. Similarly, the degree of linear polarization (P,) is defined by
The quantity ,/m represents the magnitude of linear polarization in the light field. However, it does not specify the orientation of the electric vector. Further, the degree of circular polarization (PC) is defined by
Also note that P = JPf
+ Pf.
2.2. Scattering medium In the theory of radiative transfer of polarized light, the scattering properties of a medium are defined by a quantity called the scattering matrix. The 4 x 4 scattering matrix transforms the Stokes vector of an incident light beam to the Stokes vector of a scattered light beam. The scattering matrix can be obtained experimentally,** or derived in a rigorous manner.” For scattering by a medium composed of spherical scatterers, this matrix can be derived rigorously by use of the Maxwell equations with appropriate boundary conditions.?’ The form of this matrix is given by P,lW
P*d@)
P*d@)
P,,W
0 0
sz,(@) &(@)
(5)
where Z(0) is a function of the included angle (0) between the incident and scattered light beams as shown in Fig. 1. All the elements of the scattering matrix, Z(O), are real. The scattering matrix, Z(O), operates on the incident Stokes vector (the I axis is assumed to be parallel to the scattering plane as shown in Fig. 1). Note in Fig. 1 that the shaded region is the scattering plane. Each of the four independent functions in Eq. (5) is a function of the size parameter (s = n,,nDji., where D is the geometric diameter of the scatterer) and the relative complex index of refraction (m = n - ix) of the spherical scatterer. Note that n,, is the refractive index of the medium surrounding the scatterers and is taken to be 1.333 in this study. Each of these four functions can be expressed in terms of a series of Legendre polynomials, and the coefficients of this series have been obtained rigorously. 24A number of standard computer codes are available for the purpose of computing the coefficients of the Legendre series expansions of the above elements.” It should be noted that p,,(O) is normalized as
&
p,,(O) dn = 1
s4n where fi is the solid angle and p,, is the scalar Mie phase function.” In the case of Mie scattering, examples of the four elements of the scattering matrix are shown in Fig. 2 for a variety of size parameters.
174
A. Ambirajan and D. C. Look
Incident beam
Fig. 1. Geometry of the rotation angles necessary in the phase matrix. Note that the shaded plane is the scattering plane. Light beams propagate in the direction r x I.
In multiple scattering problems, the I axes of the Stokes vectors, both incident and scattered, are assumed to be parallel to the meridonal planes (planes of constant 4 as seen in Fig. l), and designated 4, and & respectively. To compute the scattered Stokes vector at a particular location, the coordinate frame of the incident light beam must be rotated by an angle z - cr, (see Fig. 1) such that the 1 axis is parallel to the scattering plane. The scattering matrix then operates on this rotated Stokes vector. A subsequent rotation by - cZ brings the 1 axis parallel to the meridonal plane of the scattered light beam. The above rotations are mathematically equivalent to pre- and post-multiplying the scattering matrix by rotation matrices. Various papers outline this operation ,3,2’hence only the final result will be given here. If the rotation operator is denoted by L(a),for a rotation of U, then the phase matrix is given by (7) In the above equation, note that p is the cosine of 13and that (Tcan be found by using spherical trigonometry. The relevant formulae for (T,0 and L(a) are listed in appendix A. So, if C, represents cos 2a,, C, equals cos 2a,, S, equals sin 20, and S, equals sin 2a,, c, PA(@)
Pl,W w2,hw#4
GP2,W
C2P,,WG
GP2lW
S2P,,@>G
-
-
0
SIP2,W
~2~2dW,
-
C~P~@)S
-
+ C2~2,vvs
-
f32p,,WS
+ C2~2,WG
s232,WG
W2,W
=
0
Wd@)
GM@)
s24,(@) s2,w
’ I (8)
For Rayleigh scattering, the d2, term is zero, which results in the decoupling of the circular polarized component. I4 The phase matrix is the vector analog of the phase function used in the scalar radiative transfer theory. 3. MULTIPLE
SCATTERING
OF POLARIZED
LIGHT
In this section, the transport equation will be described along with a solution in terms of the Successive Order of Scattering Series (SOSS). This series will later be evaluated by the Monte Carlo
Multiple scattering of a polarized laser beam
-Rayleigh
_______.~&43~ . . .-. . . x=0.813 10’ : ;
m=l.l97+Oi
..X=1.5?2 ------x=2.379 ---x=3.264
,
0.5.
175
m=l.l97+Oi.
-Royleigh _**.430 -------x&813
.. ..__
X=1.592
------x=2.379 ----x=3.264
‘*r-----_____-__~4).3,3 rd.592 ------1~2.379 -x=3.264
8
-I m=l.l97+Oi
-Rayleigh -------xdl.430
10
i 6 K
:
4
:
i
P
CL
Fig. 2. Elements of the Mie phase matrix for m = 1.197 + Oi.
technique. The problem considered in this paper is to describe the multiple scattering of a polarized narrow beam of light (such as a laser) incident on a plane-parallel scattering medium. Consider a polarized incident light beam with a spatial distribution given by
I,“C(X,Y,W+#J) = LWPc_YP(P - lP(dJ).
(9)
The term I, indicates the Stokes vector of the incident light beam, and the terms involving the delta functions indicate that the light beam is concentrated at the origin of the coordinate system and is incident normal to &he XY plane. The above stated problem is illustrated in Fig. 3. Mathematically, the backward Monte Carlo technique is equivalent to Monte Carlo integration of the SOSS.‘7~26 In general, the Stokes vector I, at any location and in any direction of propagation in a scattering medium can be represented by an SOSS relationship like
I, = f1. = I, + I , a JQSRT 58 ?-~B
(10)
A. Ambirajan and D. C. Look
176
where the subscript t denotes the total (inclusive of all orders of scatter) Stokes vector at some location which is the sum of the diffuse and collimated components. The diffuse component of the light field is denoted without the subscript “t” and is given by
(11) Further, I, is the Stokes vector for photons which have been scattered n times and is represented by a set of multidimensional integrals. The zeroth-order Stokes vector, I,,, refers to the contribution of the incident polarized searchlight beam to the total Stokes vector and is given by’* = WX)&.Y)&~ - 1)@4) e-’ .
I&,YJ;&)
(12)
The Monte Carlo method is used to evaluate the integrals associated with each order of scatter. Thus the purpose of this section is to present these integrals in a general form. The derivation of these integrals follows a similar procedure used for the unpolarized searchlight beam problem.‘8 In the subsequent sections, a 4 x 4 matrix called the “Stokes matrix” (M) will be evaluated. The Stokes matrix is distinct from the scattering matrix, Z, and the phase matrix, P. The latter two matrices refer to the response of a single scatterer to an incident light beam whereas the Stokes matrix refers to the response of a multiple scattering medium to an incident light beam. The scattering matrix (Z), phase matrix (P) and the Stokes matrix (M) all have the same properties as a Mueller matrix.2’ The scattered beam is thus described by a Stokes matrix defined by the equation: I = MI,.
(13)
Once M is found, Eq. (13) can be used to obtain the scattered Stokes vectors, I. The advantage of this method is that once M is found, one can operate on M with I, (the incident polarization state) to find the resultant Stokes vector for any incident state of polarization. Further, similar to the definitions for the various orders of Stokes vectors (I”), various orders of Stokes matrices (M,) can be defined as M,=iM,=M,+M 0
I,
w
Wl) wi9
S(Y) I
cI L
Scattering layer Fig. 3. Geometry under consideration.
(14)
Multiple scattering of a polarized laser beam
177
Fig. 4. Coordinate system.
Note that,
as in Eq. (13)
I, = ;, I,, = i M,I, = MJ, , n=ll
( 15a)
,I=0
I = f I, = i MJ, = MI,, “=I PI= 1
(15b)
1, = MA, .
(15c)
and
In the subsequent
subsections,
formulae
will be presented
for M, M,,, M,, Mz and M,,, respectively.
3.1. Optical coordinate system In this paper, all coordinates will be specified in optical units. An optical coordinate is defined as the product of the extinction coefficient, be, and the geometric coordinate. The optical coordinate in the X direction is defined to be x, (i.e., x = &I’where fie is constant). The same convention will apply for the Y and Z directions. The rectangular optical coordinate system chosen for this problem is shown in Fig. 4, which depicts a beam of light, originating at the point (x,,y,,z,,), passing through point (xJ,z). The azimuthal angle of the direction of propagation of the beam of light is 4, the polar angle is 19(p is the cosine of this polar angle). The azimuthal and polar angles are specified with respect to a local coordinate system (X/J’,; ‘) with its origin at (XJJ). A useful relationship will be obtained later by noting now that x0 = x + (zO - z) tan (3 cos C#I,
(164
y. = y + (z,, - z) tan 0 sin C$
(16b)
178
A. Ambirajan and D. C. Look
and tan
e = J1 - j.?/p
.
(16~)
Equations such as these will be used to specify the straight lines along which beams of light travel. In general, the radiation field is a function of both the direction of propagation (P,c$) and the location coordinates (x,y,z). Thus the Stokes vector (I) and the Stokes matrix (M) are usually represented as I(x,y,z;p,t$) and M(x,y,z;&), respectively. 3.2. Vector radiative transport equation Because the intensity at any location can be broken into the sum of a diffuse component and a collimated component,‘8 the intensities without subscripts represent diffuse components only. The matrix form of the vector radiative transport equation follows as Eq. (17). Note that the order of the source terms is indicated.
P
$f +M(x,y,z;wb)
=
z P(P,A~,W(X)~(Y)e-’ (17)
The relevant boundary conditions are for ,U> 0 WV,O;P,+)
= 0
(18a)
M(s,y,&,#~)
= 0.
(18b)
and for p < 0
The bold faced 0 of Eqs 18a and 18b is the zero or null matrix and is given by 0
o=
0
0
0000. 0 0 0 0 0 i 0
0
0 0
I
(19)
Thus, as seen in Eq. (14) the total Stokes matrix at any location is the sum of the diffuse and collimated components, CoIlmated
DlfTUX
where the first term on the right-hand side represents the attenuated but unscattered Stokes matrix (M,). The second term on the right-hand side represents the diffuse part of the light field (M), thus
M = f M,,(x,y,wc,~) II =
(21)
I
In this paper, Eq. (17) will be solved for the diffuse Stokes matrix (M). Note that K is the 4 x 4 identity matrix [i.e., K = diag(l,l, l,l)]. 3.3. The first order of scatter The derivation of the expression for the first order of scatter is similar to that for the scalar searchlight problem,” and will be omitted here for brevity. Thus the first-order diffuse Stokes matrix (M,) is given by
Multiple scattering of a polarized laser beam
[ jf$
x exp -
179
x2 + y2 6[+M - arctan(y/x)]
yI < p < y2 ,
(22)
d---l
where YI=
L-z
-
(234
x2 + y* + (L - z)*
and
y2=J-&q.
Wb)
3.4. The second order of scatter
As with the preceding section, Ref. 18 presents details of the derivation. The formula of the second order of scatter for the diffuse Stokes matrix (M,), for p > 0, is
M&,)~,r;/4)
= ( E >‘“-;_>“’
[[arcsin
- arcsm(y,)]
where
z, = &
ln[a(e”~)-
l)+
11.
(25)
Similarly, for the p < 0 case, the double scattered diffuse Stokes matrix is Ze-: -~lL~‘Lm””
l[arcsin(y?)
- arcsin(
where
(27) Note that, in the above equations, j_4= arcsin{ /I[arcsin(y,) - arcsin( x, = x + (z, - 2)
Jl-P2 P
+ arcsin(
cos 4
9
,
(284 (28b)
180
A. Ambirajan and D. C. Look
Yl = y +
YI=
(2, - z)
Jl -CL2sin 4
,
P
L - z,
-
(284
x: + y: + (L - zJ2 ’
Pe)
and 4, = arctanQ,/x,)
.
(280
3.5. Higher orders of scatter
As mentioned previously, the equations for M, (n > 2) can be derived in a similar manner to that done in an earlier paper” and only the resultant form will be presented here. Thus, for p > 0 M,(x,y,z;p&)
=
~(1 - e-““)
x M,-I(xn-I,~n-l,~n-l + w-~~L~)
dLI @-, da,-, . (294
Similarly, for the case of p < 0, M,(x,Y,z;~,~)
=
041 - e-(‘-L”r)
x Mn-,(xn-,a-,,$,-,
+w,-,AJdL,
dyn-, da,-,
W’b)
where PL,-I =2/L, 4”-, X,-I =x+*,_,
Wd
- 1,
= 27v-, JI-Tci P
W’b)
9
(3Oc)
cosl#iJ
and (304
For p < 0, the expression for
rc/,,_ , is $.
, = p ln[a, _ ,( 1 - e -‘i,‘) + e - Z,l’l]
(30e)
and for ~1< 0, the expression for tin _, is $, _ , = p ln[g _ ,( 1 _ e - (; - L)!tl)+ e - CT - L)rLl]
(300
3.6. Arbitrary incident spatial distribution In many practical situations, the incident light beam may not be concentrated at one spatial location as presented in Eq. (9), instead the incident beam may be spread out spatially.‘7~28This situation can be handled in a manner analogous to that used in Ref. 18. Thus all orders of the diffuse Stokes vector field for such a situation (spread out though collimated incident light beam) can be computed by using the theory developed for the searchlight problem. So, the situation previously discussed for the case of a normally incident delta function loading [see Eq. (9)] can be generalized to that of a normally incident beam of arbitrary incident spatial distribution of intensity and polarization. This can be achieved by using the superposition principle.29,30 The incident beam is assumed to have the form 1inc(x~Y~z;P~4)
=
lof(x#P(x)G(Y)G(P
-
1)6(4)
7
(31)
Multiple
scattering
of a polarized
181
laser beam
wheref(x,y) is the incident spatial distribution of energy. Let all of the total (inclusive of all orders) diffuse Stokes vector field due to this spatial distribution be represented by I,(.u,y,z;~,~).The corresponding Stokes matrix for this case is M,(x,y,z;p,t$).If the searchlight loading is visualized as an impulse load, then the total (inclusive of all orders) diffuse Stokes matrix for this case, M(.ug,z;p,c#~),can be considered as the impulse response. Hence, by superposition, the multiply-scattered Stokes matrix for an incident distribution given by Eq. (31) is obtained as f(s’,f)M(x
- X’J - y’,q&)
dx’ dy’ .
(32)
The above equation is analagous to conoo1z+zg2ythe impulse response with the incident load distribution. 4.
MONTE
CARLO
ESTIMATOR
Before outlining the Monte Carlo estimator characteristics of the problem being presented, a brief outline of the basis behind the Monte Cario estimator used in this paper will be given and then the various estimators will be introduced. The goal of this section is to present estimators for the first-, second- and nth-order Stokes matrices, namely M,, M, and M,, respectively. In this section, the symbol (I) N implies the mean value of Z after N simulations. The symbol E[M] is the estimated value of the integral involved in the matrix M. Consider the integral
s
’wfx)g(x) dx ,
Q=
(33)
0
where w(x) is the weighting factor such that w(x) > 0
(34)
and
s
w(x) dx = 1 .
(35)
0
Thus w(x) has the characteristics of a probability density function” while g(x) is an arbitrary function. If x, is a random variable from a distribution w(x,), then an estimator of Q is given by
edPI = &xl)
(36)
and the mean value of the estimator after N events is
(Q>N = (l/W~eJQl. ,=I
Alternatively,
(37)
it is possible to sample xi from a uniform distribution, and define the estimator as
edQ1= dO+fx,) .
(38)
The latter estimator [i.e., Eq. (38)] will, in general, lead to a larger variance in the solution when compared with use of the former estimator [i.e., Eq. (36)]. Thus the selection of Eq. (36) to represent the estimator will result in a smaller value and is therefore termed a variance reduction technique.” Consider the cumulative distribution of w(x)
u = F(x)
=
s
wb)dy
0
Olxll
(39)
182
A. Ambirajan and D. C. Look
and the inverse transformation x=F-‘(U).
(40)
Thus Eq. (33) can be rewritten, by changing variables, as Q= [g[F-‘(U)]dU.
(41)
If U is sampled from a uniform distribution, then this is equivalent to sampling x from the distribution w(x).” Thus the estimator has the form
4Ql = gF’- ‘(UJI9
(42)
where U, is the Chsample from a uniform distribution and is equivalent to the estimater given by Eq. (36). This is the basis for O’Brien’s26 assertion that most variance reduction techniques are essentially equivalent to integral transformations. This is also the reason for the rather lengthy transformations carried out in the derivation of the integral forms for M,, M, and M, in Ref. 18. Now, from Eq. (6) and the fact that pI, > 0, it can be seen that the element p,, of the scattering matrix [Eq. (5)] has the characteristics of a probability density function. It would thus seem reasonable to sample the direction of a scattered photon from this “distribution”. Remember that p,, is the scalar Mie phase function.23 In the stochastic simulation of scalar (unpolarized) photon trajectories, the direction of the photon after each scattering event is obtained by sampling the scalar phase function32 (p,,). The same procedure will be used in the present paper to simulate the trajectory of a polarized photon. In essence, the trajectory of the polarized photon will be identical to that of a scalar photon. However, the weights associated with the trajectory will be different from the scalar case.” Another algorithm for the simulation of polarized photon trajectories, proposed recently, is the method of symmetrized trajectories. 33This method apparently produces reduced variances for problems involving spheroidal shell atmospheres. It is not used here. For the second order of scatter, the estimators are obtained based on Eqns (24) and (26). Thus the estimator for the double scattered Stokes matrix when ,u > 0 is e,[M,(x,y,z;p&)]
= (z
)I ‘~1~~“’
[arcsin(y2)
- arcsin(
where tl and p are sampled from uniform distributions. Note that x,, y,, 4,, y, and yz are obtained from Eqs. (28). Similarly, for ~1< 0 the following estimator is obtained 2e-:_e
1 ILi’L-““’ [arcsin
- arcsin(
exp[ - js
JG] . (44)
Since the direction of the sampled photon is obtained from p,,, Eq. 29a indicates that for p > 0
where Ii/b_, is obtained by selecting a random number from a uniform distribution and x;,_, and JJ;_, are derived from I&-, by means of Eqs. (30). pb _, and #, _, are obtained by sampling the phase function distribution pII by using an algorithm suggested by Barkstrom.-‘? The process of dividing the estimator in Eq. (45) by p,, is a consequence of the fact that the sampling is obtained
183
Multiple scattering of a polarized laser beam
from this (p,,) distribution. estimator for p < 0 is
This process is often called Importance
Sampling.3’ Similarly, the
(46) In Eqns (45) and (46) the angle 0 is given by cos 0 = /yin_, - (1 -/1)“2(1 -p,-,)‘~cos(~
- f#&_,).
(47)
Thus the estimator for higher order scatters is derived recursively from lower order scatters. 5. RESULTS In this section, the theory outlined above will be used to analyze the multiple scattering of a narrow, circularly polarized light beam by a plane-parallel medium composed of spherical dielectric scatterers. The backscattered component of the radiation will be studied. The spherical scatterers were assumed to have a relative refractive index of m = 1.197 + O.OOOi,which corresponds to latex particles suspended in water. A variety of situations were studied as shown in Table 1. For each size parameter, the backscattered Stokes vector was obtained for four optical thicknesses (,C) of the slab. The various elements of the Mie phase matrix are shown plotted in Fig. 2 for the size parameters of interest in this study. Further, all results presented in this paper are for an albedo (w) of one. Results for smaller albedos can easily be obtained since the intensity at each order of scatter is stored.” It should be noted that the numerical technique used computes the entire backscatter reflection matrix (M) of the medium. As pointed out earlier, M operates on the incident Stokes vector (IJ as shown in Eq. (13). If the incident state of polarization is circular, the diffuse radiation field is symetric about the laser beam.“4 In this study, the analysis will focus on the case of a circularly polarized searchlight beam incident on the medium. The Stokes vector of the incident light beam is
(48)
where 1, is the intensity of the incident beam. The aim of this study is to examine the physical nature of the scattering of polarized light and to see how multiple scattering affects the degree of polarization. Knowledge of this fact would aid in determining the validity of using scalar theories in analyzing the multiple scattering of polarized light. The quantities of interest are the elements of the Stokes vector, the degree of linear polarization (P,), and the degree of circular polarization (PC). Another quantity of interest is the mean number of scattering events (I) undergone by a photon. Since Z,(x,y,z;p,~$) is the intensity (note that this represents the first element of the Stokes vector) of light which has undergone n scattering events, Table I. Range of values of the size parameter (x) and optical thickness (L). The number of orders of scatter are the numbers in the body of the table. L Y
0.430 I .592 2.379 3.264
0.1
0.5
I.0
5.0
25 25 25 25
50 50 50 50
50 50 50 50
100 100 100 100
A. Ambirajan
184
and D. C. Look
and X:= ,Z,(x,y,z;&) is the intensity of light which has been scattered at least once (i.e., n 2 l), the probability that a photon gets scattered n times (where n 2 l), pnr is
Thus the mean number of scattering events [I(x,y,z;p,+)] undergone by a photon at a particular location of interest is the estimated value of n, and is evaluated with the equation
In the following sections, the results for multiple scattering of a circularly polarized laser beam will be presented. Comparisons with the scalar theory will also be presented. A FORTRAN code was written to implement the above theory and was executed on a HP-9000 work-station. The random numbers were generated by a multiplicative congruential generator with
(b) x=0.430 -+-~~=I.592 -x=2.379 -x=3.264
IO”
(Osqll (0.241~) (0.36Opj (0.494fi)
-
GO.430 ~I.592 x=2.379 ~~3.264
IO”
IO”
Optical radius, r
(0.065p (0.241~ (0.36Op (0.494p
lo”
IO”
Optical radius, r
(4
(cl ---cGo.430 --a-x=1.592 -x=2.379 ~~3.264
(O.Ot+) (0.241p) (0.3flOp) (0.494@
-t --o-c ---A---
Optical radius, r Fig. 5. Normalized
intensities
for a number
GO.430 x=1.592 x=2.379 ~~3.264
(0.065# (0.241~) (0.36Op) (0.494v)
Optical radius, r of size parameters
and optical
thicknesses.
18.5
Multiple scattering of a polarized laser beam
shuffling, and a multiplier of 16 807.35 For each case considered in this paper, 1 400 000 photons were traced. The computer code used in this study was validated in two ways. The first method of validation was to assume a phase matrix such that only the p,, term was non-zero. Thus the intensity predicted by the code should yield results which are identical to the scalar case for the Z term. Further, the Q, U, and V terms should all be zero. This indeed was the case. Secondly, to check the validity of the polarized part of the code, the data were compared with the polarized one-dimensional code of Evans? This was done by convolving [see Eq. (32)] the polarized searchlight solution with a uniform distribution [f(x,y) = l] over the incident surface for the phase matrices of interest in this study. The results from the polarized searchlight code compared well with the one-dimensional data. 5.1. ZntensitJ
The radial (see Fig. 3) variation of the backscattered polarized intensity appears to have the same form as that observed for the scalar problem. 27.28 This can be seen in Fig. 5, where the normalized intensity increases and then decreases with an increase in the optical radius (r). Some interesting features can be noted. For increasing particle size parameters, the backscattered intensity decreases. This can be attributed to the fact that larger particles have a larger asymmetry factorz3 which implies that a greater portion of the energy of the incident beam is scattered in the forward direction. For larger optical depths [see Fig. 5(d)], the backscattered intensity increases because a larger proportion of the energy of the incident beam is backscattered.
. A v x
40 -
5
IA.5
x=0.430(0.065&L) x=1.592(0.241~) ~2.379 (0.36Op) ~3.264 (0.494~)
CJe1.0
I”=( 1.0,O.l)
z ee F “, ‘:
30.4
1 201, 1
=
A
*
A
A:.
n
: v
I.. IO-
,
‘X ov IO“
%
.
I
,xX y
Y”
x
;
-
v )i
+ 7,
v IO”
IO”
ld
Optical radius, r 50
. A V q
40 :
E %
30
.
..,,....,
. . . . . ...
. . . . .. IA.0
x=Oo.430 (0.065p) x=1.592(0.241~) x=2.379(0.36Q1) x=3.264(0.49+)
.,. 6Fl.O
I,=( l,O,O,I1
.
1
3
F 7 20‘:I A;:
;::
A ’
t’
b t
=
IO-
, w IO”
“XI .w . . - v lo”
I x ; .r...J IO’
7
‘I
_,, ,,
t
,:.
loo
Optical radius. r
Fig. 6. Percentage difference between the vector and scalar radiative transfer solutions for the diffuse light field.
186
A. Ambirajan and D. C. Look
-
-x=0).430 -x=1.592 -x=2.379 -x=3.264
x=0).430(0.065 x=1.592(0.241 x=2.379(0. x=3.264(0.
(0.065&l) (0.2418) (0.3m) (0.494N)
Optical radius, r
r t
0.8
0.6
L=l.a=l.O
1
Io=ll,o,o,l)
j
-c
-c x=0.430(0.065~) --a-.-x=1.592 (0.241~) --c x=2.379 (0.36@) -----b x=3.264(0.499)
-
x=0.430(o.Ofjsp) x=1.592(0.241~) x=2.379(0.36Op) x=3.264(0.494~)
d
0.0: lo”
IO”
IO”
Optical radius, r
IO0
lo”
IO”
I@’
IO”
Optical radius, r
Fig. 7. Degree of linear polarization for a number of size parameters and optical thicknesses.
In Fig. 6 the vector and scalar” estimates of the intensity for two optical thicknesses (L) are compared. Clearly, for smaller size parameters, the difference between the vector and scalar theories is significant whereas it is smaller for larger size parameters. Also, it appears that for very large optical radii the difference between vector and scalar solutions decreases. However, this does not necessarily imply that the diffuse light field at large optical radii or small size parameters is unpolarized. Quite the contrary is observed in the following presentation. Thus even though the vector and scalar solutions approach each other as far as the intensity is concerned, the diffuse light field may still be polarized. In other words, the Q, CJ, and V terms may be significant with respect to the intensity I. This indeed was the case in our numerical simulations. 5.2. Degree of linear polarization In Fig. 7 the degree of linear polarization for the above cases is outlined. The behavior here is very interesting. The degree of linear polarization at small optical radii is the same irrespective of optical depth. This is because at small optical depths, the multiple scattering is primarily second-order (see Fig. 8). When the optical radius is small (i.e., x’ + JJ* is small), the quantities y, and y2 in Eqs 23a and 23b tend to one and zero, respectively, when z = 0 (upper surface) and x2 + y2+0 (small optical radii). Thus the double scattered Stokes matrix for backscatter [Eq. (26)]
187
Multiple scattering of a polarized laser beam
becomes essentially independent of optical depth. However, the variation with respect to particle size is harder to explain, particularly at small radii. One clue to this variation is the sharp variations in the functions comprising the phase matrix for Mie scatterers as shown in Fig. 2. This type of oscilation with respect to particle size for side scatter in a dilute suspension was first observed experimentally by Look and Chen.j4 For an optical depth of 0.1, the degree of linear polarization at large optical radii exhibits a behavior that is monotonic with respect to particle size. The degree of linear polarization decreases with an increase in particle size. The case of x = 3.264, however, shows an anomalous behavior in the sense that the data are highly oscillatory at even small optical radii. This is because the intensity is very low in this region for x = 3.264. Thus when the degree of linear polarization is computed by dividing by the intensity, the errors in this quantity lead to unpredictable results. Similar trends are observed for larger optical thicknesses. However, the degree of linear polarization decreases for a given size parameter with an increase in the optical thickness. The interesting feature in Fig. 7 is due to the fact that, even at reasonably large optical radii, the light field exhibits a high degree of linear polarization. In Fig. 8, the mean number of scattering events (I) for a medium with an optical thickness of 5.0 at an optical radius of 1.0 is greater than 5. Thus even for quite a large number of scattering events (see Fig. 8) the light field retains its polarized character (Fig. 7).
25.0
I
,,
1
111..(
‘....‘,
,....,.,
,?-
M.1, 20.0 -
__ .-...----------
-
uF1.0
.
lo=1 1,O,O.l I
. ----...-. ---
GO.430 (0.065p) x=l.592(0.241~) x=2.379 (0.3fX&, x=3.264 (0.494P)
15.0 -
x=0 430 (O.Oql) n=I.J92(0.241p) x=2.379 (0.360)~) x=3.264 (0.494~)
15.0 -
L
L 10.0 -
. 10.0 -
5.0 -
0.0 lo”
, “...,J
, ‘...*d 10-l
I@*
0.0’
“Lei,J’ Id
“11,,1’
lo”
1V’
Optical radius, r
‘ll(l.f lo”
‘.J
‘l(((ii’ IO”
Optical radius, r
25.0
. . . . ..I
.,.,..,
v...,,
r
Ls5.0, co=l.O _ I,=(l.O.O,l) 20.0
15.0
0.0’ lo”
20.0 GO.430 (0.065~) ._------. x=1.592 (0.241~) -----.x=2.379 (0.3cXy~) .. x=3.264 (0.494~)
’
,....A lo”
.
..‘,,,I
......--. p
.
x=0.430 (0.065p) x=1.592 (0.241~) x=2.379 (0.3w) x=3.264 (0.4%)
IS.0 -
I
lo“
Optical radius, r
,..1,d loo
.
,I
I 10“
. . .,..I
. . . ,,,, I
IO“
loo
Optical radius, r
Fig. 8. Mean number of scattering events (r) for a number of size parameters and optical thicknesses.
188
A. Ambirajan
and D. C. Look
x=0.430(0. -x=1.592 (0. --c ~22.379(0. --b ~3.264 (0.
-
0.5
-----o-
-
x=0.430(0.065p) x=1.592(0.241~) x=2.379(0.3m) ~3.264 (0.494~)
.I..-n
lo”
lo”
IO’
lo” Optical radius, r
lo“
l@’
IO”
Optical radius, r
“OV] xd.430 -x=1.592 ~2.379 ~3.264
lo”
-t -
(O.O+&) (0.241~) (0.3q) (0.494~)
100
10”
Optical radius, I Fig. 9. Degree
of circular
polarization
x=0.430(O.OLqL) x=1.592(0.241~) x=2.379(0.36Op) x=3.264(0.4%)
lo”
10-l Optical radius,
for a number
of size parameters
IO”
I
and optical
thicknesses.
5.3. Degree of circular polarization
In Fig. 9, the degree of circular polarization is shown. In this case PCis identical for all optical depths at small optical radii. The anomalous behavior for a size parameter of s = 3.264 can be explained in a similar manner to the previous paragraph. At large optical radii, the trend is an increase in the degree of circular polarization for larger size parameters. This can be explained as follows: for small size parameters, the dz, term in the scattering matrix [see Eq. (5)] is negligible as shown in Fig. 2. This implies that the circular component of the vector radiative transfer equation is essentially decoupled.14 Thus the circular component (V) of the vector radiative transfer equation behaves as a scalar intensity with an effective albedo of less than one. It thus approaches zero much faster than the actual intensity (I). However, with an increase in the size parameter, the quantity d2, becomes more significant. The physical effect of d?, is the coupling of the circular and linear components of the light field, as pointed out by Zege and Chaikovskya.‘5 Thus, for a larger size parameter, an increase in the conversion of linear light to circular light can be expected. Another interesting feature to be observed is that the sign of the degree of circular polarization goes negative as the optical radius is increased. In other words, the direction of rotation of the electric vector changes at larger optical radii. However, the degree of circular polarization does not seem to decrease appreciably at large optical radii for increased optical depths of the medium. Again, the interesting physical effect is that the light field is polarized at large optical radii.
189
Multiple scattering of a polarized laser beam
5.4. Analysis for various orders of scatter In the above analysis, it is seen that the light field retains its polarized character at quite significant optical radii from the incident light beam. This is a very interesting phenomenon because one would expect that, at large optical radii from the beam, the increased number of scattering events would increasingly depolarize the light field. This does not seem to be the case. To understand this phenomenon further, the Stokes vector will be analyzed with respect to the order of scatter. Plots of some polarization parameters for two optical depths of the medium are shown in Figs 10 and 1I. The data in these figures are plotted against the order of scatter (N). The quantity I, [see Fig. IO(a) and Fig. 1l(a)] is the backscattered intensity for photons that have been scattered N times and the quantity ,/m [see Fig. IO(b) and Fig. 1l(b)] is the backscattered linearly polarized light that has been scattered N times. It can be observed that IN and dm decrease after N = 8 exponentially with respect to the order of scatter. A similar trend was observed for scalar multiple scattering in one-dimensional problems” and for scalar multiple scattering by a searchlight beam.‘” The linearly polarized component (JmV) and the on the same graph in Fig. 10(c) and Fig. 1l(c). The ratio of these two quantities (
0
’
“I’,
”
10
20
s ‘a
30
Order of scatter, N
1.01 L 0.9 :
I ,
Cd) I
I ,
k0.S ~1.0 x=2.379(0.36%) r&O01
0.8
loa
0.6
. .
lo”
. .
0.5 :
7 lo* r
b
. v
. ’
lo“
. .
101
1, w (Q;+U5’” l
.
0
’ 5 Order
.
8 ’ 10
of scatter, N
.
0.3 :
. .:
loa IV0
:
0.4 : . .
0.2 : . . v . . . . ’
.
_) . 8 . 3
.
15
0
. l m. ‘L”“““‘._ s
10
15
Order of scatter, N
Fig. IO. Various polarization parameters plotted as a function of the order of scatter for an optical depth of L = 0.5.
A. Ambirajan and D. C. Look
190
. . .,.
.
I,..
1
.I....,.
.
.
W .
L.&o ~1.0 x=2.379(0.360~)
Order of scatter, N
Order of scatter, N
(cl
1 ” ” I.. ” IA.0 o=l .O x=2.379(0.360~) r=O.Ool
100 . IO”
1.0.
0.9
7
(4 1
,
I
7
1
L=S.O~1.0 x=2.379(0.360~) PO.001
0.8 :
. 10-Z
0.7 y .
lo-’
..
I
. 10” lo” r lo1 r
. . .
.
. .
’ .
v.7.
8 ‘, . (Q;+U$‘”
0.6 ; 0.5 :
. . .
. S . . 1
lo” r
Order of scatter, N
Order of scatter, N
Fig. 11. Various polarization parameters plotted as a function of the order of scatter for an optical depth of L = 5.0.
I,,,) is P,N [see Fig. 10(d) and Fig. 1l(d)], which is the degree of linear polarization of a photon that has been scattered N times. For all the cases considered here, the incident light beam is circularly polarized and the intensities for higher orders converge for all optical radii. This was observed in an earlier paper by the present authors.‘* The same is observed here for the linearly polarized component. The curious feature illustrated in Fig. IO(d) and Fig. 1l(d) is that the degree of linear polarization initially decreases with an increase in the order of scatter and then actually increases after approximately the fifth order for both the optical thicknesses considered. This was observed for the other size parameters also. This would seem to go against intuition, which suggests that a higher number of scattering events should lead to a decrease in the polarization of the light field. It should be noted here that since most of the energy is concentrated at lower orders for smaller optical radii, the net (inclusive of all orders) light field will still be largely determined by the lower orders. The idea behind these plots, however, is to illustrate the fact that higher orders are polarized and would thus influence the degree of polarization of the diffuse light field. For a given order of scatter, the degree of linear polarization decreases with an increase in optical thickness. Another noteworthy point is that, with an increase in the optical radius, the higher orders of scatter gain
I91
Multiple scattering of a polarized laser beam
increasing significance in specifying the diffuse light field as shown in Fig. 10(a) and Fig. 1l(a). The increase in the polarization of higher order photons would tend to cause the the net light field to stay polarized to significant optical radii. The large scatter in the higher order terms is indicative of the decreased accuracy in the evaluation of these terms. 6. CONCLUSIONS
In this paper a backward Monte Carlo estimator was developed for the multiple scattering of a narrow, polarized light beam. It was used to analyze the multiple scattering of a circularly polarized light beam. Results from this study were compared with data in which polarization effects were ignored (scalar theory). Certain other polarization parameters were studied. Data were compared for a variety of size parameters (using Mie theory), optical thicknesses, and optical radii. It should be emphasized that the numerical data presented in this study are for the backscattered component of the diffuse light field. Significant differences were observed between the intensities computed with the vector radiative transfer theory and the scalar theory. Differences between the two theories as great as 25% were observed for the intensity and were greater for smaller optical thicknesses and smaller size parameters. The diffuse light field was polarized even at significant optical radii from the incident light source. The degree of polarization (both linear and circular) was strongly dependent on the size parameter. At smaller optical radii, where double scattering was predominant, the various polarization parameters were essentially independent of the optical thickness of the medium and depended only on the size parameters of the scatterers. The influence of the order of scatter on the various polarization parameters was presented. The intensity and other components of the Stokes vector were observed to decrease exponentially with respect to the order of scatter for the optical thicknesses considered here after N = 8. Further, it was observed that light that has been scattered more than five times is still polarized. Acknow[edgemPnrs-The authors of this paper wish to thank Dr F. K. Evans of the University making his one-dimensional nit.colarado.edu/ - evans).
polarized
radiative
transfer
code
freely
available
on
REFERENCES I. Kienle, A., Lilge, L., Patterson, M. S., Hibst, R., Steiner, R. and Wilson,
the
of Colorado for kindly world wide web (http:‘/
B. C., Appl. Opt., 1996, 35,
2304. 2. Wang, L. and Jacques, S. L., J. Opt. Sot. Am. A, 1993, 10, 1746. 3. Groenhuis, R. A. J., Ferwada, H. A. and Ten Bosch, J. J., Appl. Opt., 1983, 22, 2456. 4. Wilson, B. C. and Adams, G., Med. Phys., 1983, 10, 824.
5. Bissonnette, L. R., Bruscaglioni, P., Ismaelli, A., Zaccanti, G., Cohen, A., Benyahu, Y., Kleiman, M., Egert, S., Flesia, C., Schwendimann, P., Starkov, A. V., Noormohammadian, M., Oppel, I_. G., Winkler. D. M., Zege, E. P., Katsev, 1. L. and Polonsky, I. N., Appl. Phys., 1995, B60, 355. 6. Yamazaki, H., Tojo, K. and Miyanami, K., Powder Technol., 1992, 70, 93. 7. Carswell, A. I. and Pal, S. R., Appl. Opt., 1985, 24, 3464. 8. Look, D. C. Jr. Opt. Eng., 1989, 28, 160. 9. Look, D. C. Jr. and Chen, Y. R., J. Thermophys. Heat Transfer, 1993, 7, 631. 10. Look, D. C. Jr. and Chen, Y. R., Appl. Opt., 1995, 34, 144. II. Dogariu, M. and Asakura, T., J. Opt. (Paris), 1993, 24, 271. 12. Dogariu, M. and Asakura, T., Waves in Random Media, 1994, 4, 429. 13. Mueller, D. W. Jr. and Crosbie, A. L., submitted to J. &ant. Spectrosc. Radiat. Tram/k. 14. Mueller, D. W. Jr. and Crosbie, A. L., submitted to J. Quant. Spectrosc. Radiat. Transfer. 15. Zege, E. P. and Chaikovskaya, L. I., J. Quant. Spectrosc. Radiat. Tramfer, 1996, 55, 19. 16. Collins, D. G., Bllttner, W. G., Wells, M. B. and Horak, H. G., Appl. Opt., 1972, 11, 2684. 17. Marchuk, G. I., Mikhailov, G. A., Nazaraliev, M. A., Darbinjan, R. A., Kargin, B. A. and Elepov, B. S., The Monte Carlo Methods in Atmospheric Optics. Springer-Verlag, Berlin, 1980. 18. Ambirajan, A. and Look, D. C. Jr. J. Quant. Spectrosc. Radiat. Transfer, 1996, 56, 317. 19. Ishimaru, A., Part. Part. Syst. Charact., 1994, 11, 183. 20. Ishimaru, A., Wave Propagation and Scattering in Random Media, Vol. I. Academic Press, New York, 1978. 21. Hovenier, J. W. and van der Mee, C. V. M., Astron. Astrophys., 1983, 128, I. 22. Voss, K. J. and Fry, E. S., Appl. Opt., 1984, 23, 4427. 23. van de Hulst, H. C., Light Scattering by Small Particles. Dover, New York, 1981. 24. Dave, J. V., Appl. Opt., 1970, 9, 1888.
192
A. Ambirajan and D. C. Look
Wiscombe, W. J., Appl. Opt., 1980, 19, 1505. O’Brien, D. M., J. Quant. Spectrosc. Radiat. Transfer; 1992, 48, 41. Crosbie, A. L. and Dougherty, R. L., J. Quant. Spectrosc. Radiat. Transfer, 1982, 27, 149. Crosbie, A. L. and Dougherty, R. L., J. Q uant. Spectrosc. Radiat. Transfer, 1982, 28, 233. Hecht, E., Optics, 2nd edn. Addison-Wesley Publishing Company, New York, 1993. Spanier, J. and Gelbard, E. M., Monte Carlo Principles and Neutron Transport Problems. Addison-Wesley Publishing Co., Reading MA, 1969. 31. Kales, M. H. and Whitlock, P. A., Monte Carlo Methods, Vol. I: Basics. John Wiley & Sons, New York, 1986. 32. Barkstrom, B. R., J. Quant. Spectrosc. Radiat. Transfer, 1992, 53, 23. 33. Voshchinnikov, N. V. and Karjukin, V. V., Astron. Astrophys., 1994, 288, 1. 34. Look, D. C. Jr. & Chen, Y. R., paper presented at 6th AIAA/ASME Joint Thermophysics and Heat 25. 26. 27. 28. 29. 30.
Transfer Conference, Colarado Springs, Colarado,
1994. AIAA Paper No. 94-2093.
Users Manual, Version 1.1, Vol. 3. Houston, January 1989, Ch. 8. 35. IMSL MATH/LIBRARY, 36. Evans, F. K. and Stephens, G. L., J. Quant. Spectrosc. Radiat. Transfer, 1991, 46, 413. 37. van de Hulst, H. C., Multiple Light Scattering. Academic Press, New York, 1980.
APPENDIX
A
In this appendix formulae for the rotation angles used in the radiative transfer of polarized light will be given. Note that 0 is the included angle between the incident and scattered light beams and that cos 0, and cos 19~are represented as p, and p>, respectively. The formula for cos 0 can be obtained by using trigonometry and is given by cos 0 = cos 8, cos 92 - (1 - co? &)‘:‘(I - co? e,)’ 2cos(fjz- 4,)
(AlI
By using spherical trigonometry, the following expressions can be obtained cos u, =
cos
e2 - cos 8, cos 0 sin 8, sin 0 ’
cos u2= cos 8, -
cos
e2 cos 0
sin e2 sin 0
sin u, =
sin e2 sin(@, - &) sin 0
sin (r2=
sin 8, sin(f$, - &) sin 0
’
64’4 643)
Thus we can obtain cos 2a,, cos 20, sin 2a, and sin 2~~ by using the standard trigonometric identities cos 20 = 2 cos’ 0 - I
W)
sin 2a = 2 sin D cos u
647)
and Further, the rotation matrix L(u) is given by I L(o) =
; 0
0
0
0
sin 2u 0 cos 2u _ sin 2u cos 20 0 0 0 1I
W9
Pergamon
PII: !300224073(!37)00034-4
0022-4073/97 $17.00+ 0.00
A BACKWARD MONTE CARLO STUDY OF THE MULTIPLE SCATTERING OF A POLARIZED LASER BEAM A. AMBIRAJAN?
and D. C. LOOK
Thermal Radiative Transfer Group, Mechanical and Aerospace Engineering and Engineering Mechanics Department, University of Missouri-Rolla, Rolla, MO 65401, U.S.A. (Received 15 November 1996)
Abstract-A backward Monte Carlo estimator is developed to describe the multiple scattering of a polarized, narrow light beam by a plane-parallel medium. The case of a right circularly polarized beam is analyzed in this paper. Results indicate that the diffuse light field is partially polarized even at significant optical radii from the incident light beam. The degree of polarization of the diffuse light field is dependent on the optical thickness of the medium and the size parameter of the scatterers. 0 1997 Elsevier Science Ltd
NOMENCLATURE D = geometric diameter of scatterer e[...] = I,Q,CJ, V = I, = I= I, = I,, = K= L= L= III = M = n, = p,I,pzI,s21,d~l=
P = P, = P, =
P Z r s,~,; x X, Y,Z
= = = = = =
estimator components of the Stokes vector incident intensity Stokes vector Stokes vector for the delta function loading zeroth-order Stokes vector identity matrix optical thickness of slab rotation matrix complex refractive index Stokes matrix refractive index of liquid carrier medium elements of the Mie scattering matrix degree of polarization Degree of linear polarization degree of circular polarization phase matrix scattering matrix optical radius = m optical coordinates also used to denote the size parameter (n,nD/i,) Geometric coordinates
Greek s~~mbols aJ?.y = variables of integration for higher scatter orders
y,,y: = 6(.r) = 0= 0 = u= 4 = $ = w=
limits to p for single scattered radiation dirac delta function polar angle included angle during scattering rotation angle azimuthal angle parameter used in the arguments of I,, single scattered albedo
Subscripts
i = incident quantity t = total quantity Superscripts
’ = first scatter point i = photon number T =
transpose
tTo whom all correspondence
should be addressed. 171
172
A. Ambirajan and D. C. Look 1. INTRODUCTION
The multiple scattering of light within a material is an area of considerable interest. In medical diagnostics, the analysis of scattered laser light from tissue samples is increasingly being used as a diagnostic tool.‘m4The use of lasers in surgery also requires an understanding of the way in which laser light is scattered and absorbed by tissue.4 Lasers are also being used for remote sensing in lidar systems.5 Multiple scattering of light has also been used in process applications to determine particle concentrations.6 An issue that has not received much attention in the literature is the effect of multiple scattering on the polarization of a narrow light beam such as a laser. When a linearly polarized light beam is incident on a scattering medium, the backscattered radiant intensity exhibits a highly anisotropic pattern7-12 (the “bow-tie”). Look* observed that this bow-tie had a cosine variation. When viewed through a crossed polarizer, a four-lobed pattern was also observed. Dogariu and Asakura” showed that the bow-tie was a consequence of the polarizing effect of single scattering, and was thus mainly observable for lower orders of scatter and small size parameters. Dogariu and Asakura” also used the diffusion approximation to evaluate the cross-polarized component of backscattered light from a plane-parallel medium and obtained good results. These analyses thus indicate that polarization effects are propagated to significant distances in a scattering medium. Mueller and Crosbie have presented a formulation for three-dimensional radiative transfer with polarization by using the integral transform method.‘j They also developed certain spatial symmetry relationships for the reflection and transmission matrices. In a subsequent paper, Mueller and Crosbie14 used this integral transform method and the symmetry relationships to obtain the functional form (i.e., of the bow-tie) of the spatially varying reflection and transmission matrices for the multiple scattering of a polarized laser beam from a plane-parallel scattering medium composed of scatterers with a plane of symmetry. Their method is valid for any incident state of polarization and any axisymmetric spatial loading. Their transform technique, though exact, is time consuming to implement. In an interesting paper, Zege and Chaikovskya’5 simplified the vector radiative transport equation by using the properties of the Mie phase matrix and obtained an approximate system of equations whose solution describes the multiple scattering of a polarized narrow light beam. The backward Monte Carlo method has been used in studying the propagation of polarized light in scattering atmospheres.‘6, ” Ambirajan and Look’* suggested an algorithm to deal with the multiple scattering of an unpolarized narrow light beam neglecting polarization effects. The problem to be tackled in this paper is that of a polarized searchlight beam incident on a plane-parallel medium composed of scattering particles. In reasonably dilute suspensions, analysis of the multiple scattering of light can be carried out with radiative transfer theory.19 For higher number densities of scatterers, the analytical wave theory” may be more appropriate. In the analysis that follows, it is assumed that the scattering is independent, no coherence effects are considered, and the light beam is spatially incoherent; hence, the radiative transfer theory is used. In the analysis to be presented in this paper a technique is used that was developed earlierI to analyze the multiple scattering of a narrow light beam including polarization effects. The numerical data presented are for the case of a right circularly polarized searchlight beam. 2. SINGLE
SCATTERING
OF POLARIZED
LIGHT
In this section, some basic theory used in this paper to describe the multiple scattering process will be presented. First, the description of polarized light in terms of Stokes vectors will be outlined. Subsequently, the interaction of polarized light with a scattering medium will be characterized in terms of a matrix. 2.1. Light field In the theory of radiative transfer by light that is polarized, the state of polarization and intensity of a beam is specified by the 4 x 1 Stokes vector (I) in the following form:
I
=
G ii U’ V
(1)
Multiple scattering of a polarized laser beam
173
The elements Q, U and V define the state of polarization of the light beam and I denotes the intensity.*’ Another quantity of interest in radiative transfer of polarized light is the degree of polarization (P) of the light. This is defined as
*+ u* + v*
.JQ
I
and indicates the fraction of the light intensity that is polarized regardless of whether it is linear and/or circular. Similarly, the degree of linear polarization (P,) is defined by
The quantity ,/m represents the magnitude of linear polarization in the light field. However, it does not specify the orientation of the electric vector. Further, the degree of circular polarization (PC) is defined by
Also note that P = JPf
+ Pf.
2.2. Scattering medium In the theory of radiative transfer of polarized light, the scattering properties of a medium are defined by a quantity called the scattering matrix. The 4 x 4 scattering matrix transforms the Stokes vector of an incident light beam to the Stokes vector of a scattered light beam. The scattering matrix can be obtained experimentally,** or derived in a rigorous manner.” For scattering by a medium composed of spherical scatterers, this matrix can be derived rigorously by use of the Maxwell equations with appropriate boundary conditions.?’ The form of this matrix is given by P,lW
P*d@)
P*d@)
P,,W
0 0
sz,(@) &(@)
(5)
where Z(0) is a function of the included angle (0) between the incident and scattered light beams as shown in Fig. 1. All the elements of the scattering matrix, Z(O), are real. The scattering matrix, Z(O), operates on the incident Stokes vector (the I axis is assumed to be parallel to the scattering plane as shown in Fig. 1). Note in Fig. 1 that the shaded region is the scattering plane. Each of the four independent functions in Eq. (5) is a function of the size parameter (s = n,,nDji., where D is the geometric diameter of the scatterer) and the relative complex index of refraction (m = n - ix) of the spherical scatterer. Note that n,, is the refractive index of the medium surrounding the scatterers and is taken to be 1.333 in this study. Each of these four functions can be expressed in terms of a series of Legendre polynomials, and the coefficients of this series have been obtained rigorously. 24A number of standard computer codes are available for the purpose of computing the coefficients of the Legendre series expansions of the above elements.” It should be noted that p,,(O) is normalized as
&
p,,(O) dn = 1
s4n where fi is the solid angle and p,, is the scalar Mie phase function.” In the case of Mie scattering, examples of the four elements of the scattering matrix are shown in Fig. 2 for a variety of size parameters.
174
A. Ambirajan and D. C. Look
Incident beam
Fig. 1. Geometry of the rotation angles necessary in the phase matrix. Note that the shaded plane is the scattering plane. Light beams propagate in the direction r x I.
In multiple scattering problems, the I axes of the Stokes vectors, both incident and scattered, are assumed to be parallel to the meridonal planes (planes of constant 4 as seen in Fig. l), and designated 4, and & respectively. To compute the scattered Stokes vector at a particular location, the coordinate frame of the incident light beam must be rotated by an angle z - cr, (see Fig. 1) such that the 1 axis is parallel to the scattering plane. The scattering matrix then operates on this rotated Stokes vector. A subsequent rotation by - cZ brings the 1 axis parallel to the meridonal plane of the scattered light beam. The above rotations are mathematically equivalent to pre- and post-multiplying the scattering matrix by rotation matrices. Various papers outline this operation ,3,2’hence only the final result will be given here. If the rotation operator is denoted by L(a),for a rotation of U, then the phase matrix is given by (7) In the above equation, note that p is the cosine of 13and that (Tcan be found by using spherical trigonometry. The relevant formulae for (T,0 and L(a) are listed in appendix A. So, if C, represents cos 2a,, C, equals cos 2a,, S, equals sin 20, and S, equals sin 2a,, c, PA(@)
Pl,W w2,hw#4
GP2,W
C2P,,WG
GP2lW
S2P,,@>G
-
-
0
SIP2,W
~2~2dW,
-
C~P~@)S
-
+ C2~2,vvs
-
f32p,,WS
+ C2~2,WG
s232,WG
W2,W
=
0
Wd@)
GM@)
s24,(@) s2,w
’ I (8)
For Rayleigh scattering, the d2, term is zero, which results in the decoupling of the circular polarized component. I4 The phase matrix is the vector analog of the phase function used in the scalar radiative transfer theory. 3. MULTIPLE
SCATTERING
OF POLARIZED
LIGHT
In this section, the transport equation will be described along with a solution in terms of the Successive Order of Scattering Series (SOSS). This series will later be evaluated by the Monte Carlo
Multiple scattering of a polarized laser beam
-Rayleigh
_______.~&43~ . . .-. . . x=0.813 10’ : ;
m=l.l97+Oi
..X=1.5?2 ------x=2.379 ---x=3.264
,
0.5.
175
m=l.l97+Oi.
-Royleigh _**.430 -------x&813
.. ..__
X=1.592
------x=2.379 ----x=3.264
‘*r-----_____-__~4).3,3 rd.592 ------1~2.379 -x=3.264
8
-I m=l.l97+Oi
-Rayleigh -------xdl.430
10
i 6 K
:
4
:
i
P
CL
Fig. 2. Elements of the Mie phase matrix for m = 1.197 + Oi.
technique. The problem considered in this paper is to describe the multiple scattering of a polarized narrow beam of light (such as a laser) incident on a plane-parallel scattering medium. Consider a polarized incident light beam with a spatial distribution given by
I,“C(X,Y,W+#J) = LWPc_YP(P - lP(dJ).
(9)
The term I, indicates the Stokes vector of the incident light beam, and the terms involving the delta functions indicate that the light beam is concentrated at the origin of the coordinate system and is incident normal to &he XY plane. The above stated problem is illustrated in Fig. 3. Mathematically, the backward Monte Carlo technique is equivalent to Monte Carlo integration of the SOSS.‘7~26 In general, the Stokes vector I, at any location and in any direction of propagation in a scattering medium can be represented by an SOSS relationship like
I, = f1. = I, + I , a JQSRT 58 ?-~B
(10)
A. Ambirajan and D. C. Look
176
where the subscript t denotes the total (inclusive of all orders of scatter) Stokes vector at some location which is the sum of the diffuse and collimated components. The diffuse component of the light field is denoted without the subscript “t” and is given by
(11) Further, I, is the Stokes vector for photons which have been scattered n times and is represented by a set of multidimensional integrals. The zeroth-order Stokes vector, I,,, refers to the contribution of the incident polarized searchlight beam to the total Stokes vector and is given by’* = WX)&.Y)&~ - 1)@4) e-’ .
I&,YJ;&)
(12)
The Monte Carlo method is used to evaluate the integrals associated with each order of scatter. Thus the purpose of this section is to present these integrals in a general form. The derivation of these integrals follows a similar procedure used for the unpolarized searchlight beam problem.‘8 In the subsequent sections, a 4 x 4 matrix called the “Stokes matrix” (M) will be evaluated. The Stokes matrix is distinct from the scattering matrix, Z, and the phase matrix, P. The latter two matrices refer to the response of a single scatterer to an incident light beam whereas the Stokes matrix refers to the response of a multiple scattering medium to an incident light beam. The scattering matrix (Z), phase matrix (P) and the Stokes matrix (M) all have the same properties as a Mueller matrix.2’ The scattered beam is thus described by a Stokes matrix defined by the equation: I = MI,.
(13)
Once M is found, Eq. (13) can be used to obtain the scattered Stokes vectors, I. The advantage of this method is that once M is found, one can operate on M with I, (the incident polarization state) to find the resultant Stokes vector for any incident state of polarization. Further, similar to the definitions for the various orders of Stokes vectors (I”), various orders of Stokes matrices (M,) can be defined as M,=iM,=M,+M 0
I,
w
Wl) wi9
S(Y) I
cI L
Scattering layer Fig. 3. Geometry under consideration.
(14)
Multiple scattering of a polarized laser beam
177
Fig. 4. Coordinate system.
Note that,
as in Eq. (13)
I, = ;, I,, = i M,I, = MJ, , n=ll
( 15a)
,I=0
I = f I, = i MJ, = MI,, “=I PI= 1
(15b)
1, = MA, .
(15c)
and
In the subsequent
subsections,
formulae
will be presented
for M, M,,, M,, Mz and M,,, respectively.
3.1. Optical coordinate system In this paper, all coordinates will be specified in optical units. An optical coordinate is defined as the product of the extinction coefficient, be, and the geometric coordinate. The optical coordinate in the X direction is defined to be x, (i.e., x = &I’where fie is constant). The same convention will apply for the Y and Z directions. The rectangular optical coordinate system chosen for this problem is shown in Fig. 4, which depicts a beam of light, originating at the point (x,,y,,z,,), passing through point (xJ,z). The azimuthal angle of the direction of propagation of the beam of light is 4, the polar angle is 19(p is the cosine of this polar angle). The azimuthal and polar angles are specified with respect to a local coordinate system (X/J’,; ‘) with its origin at (XJJ). A useful relationship will be obtained later by noting now that x0 = x + (zO - z) tan (3 cos C#I,
(164
y. = y + (z,, - z) tan 0 sin C$
(16b)
178
A. Ambirajan and D. C. Look
and tan
e = J1 - j.?/p
.
(16~)
Equations such as these will be used to specify the straight lines along which beams of light travel. In general, the radiation field is a function of both the direction of propagation (P,c$) and the location coordinates (x,y,z). Thus the Stokes vector (I) and the Stokes matrix (M) are usually represented as I(x,y,z;p,t$) and M(x,y,z;&), respectively. 3.2. Vector radiative transport equation Because the intensity at any location can be broken into the sum of a diffuse component and a collimated component,‘8 the intensities without subscripts represent diffuse components only. The matrix form of the vector radiative transport equation follows as Eq. (17). Note that the order of the source terms is indicated.
P
$f +M(x,y,z;wb)
=
z P(P,A~,W(X)~(Y)e-’ (17)
The relevant boundary conditions are for ,U> 0 WV,O;P,+)
= 0
(18a)
M(s,y,&,#~)
= 0.
(18b)
and for p < 0
The bold faced 0 of Eqs 18a and 18b is the zero or null matrix and is given by 0
o=
0
0
0000. 0 0 0 0 0 i 0
0
0 0
I
(19)
Thus, as seen in Eq. (14) the total Stokes matrix at any location is the sum of the diffuse and collimated components, CoIlmated
DlfTUX
where the first term on the right-hand side represents the attenuated but unscattered Stokes matrix (M,). The second term on the right-hand side represents the diffuse part of the light field (M), thus
M = f M,,(x,y,wc,~) II =
(21)
I
In this paper, Eq. (17) will be solved for the diffuse Stokes matrix (M). Note that K is the 4 x 4 identity matrix [i.e., K = diag(l,l, l,l)]. 3.3. The first order of scatter The derivation of the expression for the first order of scatter is similar to that for the scalar searchlight problem,” and will be omitted here for brevity. Thus the first-order diffuse Stokes matrix (M,) is given by
Multiple scattering of a polarized laser beam
[ jf$
x exp -
179
x2 + y2 6[+M - arctan(y/x)]
yI < p < y2 ,
(22)
d---l
where YI=
L-z
-
(234
x2 + y* + (L - z)*
and
y2=J-&q.
Wb)
3.4. The second order of scatter
As with the preceding section, Ref. 18 presents details of the derivation. The formula of the second order of scatter for the diffuse Stokes matrix (M,), for p > 0, is
M&,)~,r;/4)
= ( E >‘“-;_>“’
[[arcsin
- arcsm(y,)]
where
z, = &
ln[a(e”~)-
l)+
11.
(25)
Similarly, for the p < 0 case, the double scattered diffuse Stokes matrix is Ze-: -~lL~‘Lm””
l[arcsin(y?)
- arcsin(
where
(27) Note that, in the above equations, j_4= arcsin{ /I[arcsin(y,) - arcsin( x, = x + (z, - 2)
Jl-P2 P
+ arcsin(
cos 4
9
,
(284 (28b)
180
A. Ambirajan and D. C. Look
Yl = y +
YI=
(2, - z)
Jl -CL2sin 4
,
P
L - z,
-
(284
x: + y: + (L - zJ2 ’
Pe)
and 4, = arctanQ,/x,)
.
(280
3.5. Higher orders of scatter
As mentioned previously, the equations for M, (n > 2) can be derived in a similar manner to that done in an earlier paper” and only the resultant form will be presented here. Thus, for p > 0 M,(x,y,z;p&)
=
~(1 - e-““)
x M,-I(xn-I,~n-l,~n-l + w-~~L~)
dLI @-, da,-, . (294
Similarly, for the case of p < 0, M,(x,Y,z;~,~)
=
041 - e-(‘-L”r)
x Mn-,(xn-,a-,,$,-,
+w,-,AJdL,
dyn-, da,-,
W’b)
where PL,-I =2/L, 4”-, X,-I =x+*,_,
Wd
- 1,
= 27v-, JI-Tci P
W’b)
9
(3Oc)
cosl#iJ
and (304
For p < 0, the expression for
rc/,,_ , is $.
, = p ln[a, _ ,( 1 - e -‘i,‘) + e - Z,l’l]
(30e)
and for ~1< 0, the expression for tin _, is $, _ , = p ln[g _ ,( 1 _ e - (; - L)!tl)+ e - CT - L)rLl]
(300
3.6. Arbitrary incident spatial distribution In many practical situations, the incident light beam may not be concentrated at one spatial location as presented in Eq. (9), instead the incident beam may be spread out spatially.‘7~28This situation can be handled in a manner analogous to that used in Ref. 18. Thus all orders of the diffuse Stokes vector field for such a situation (spread out though collimated incident light beam) can be computed by using the theory developed for the searchlight problem. So, the situation previously discussed for the case of a normally incident delta function loading [see Eq. (9)] can be generalized to that of a normally incident beam of arbitrary incident spatial distribution of intensity and polarization. This can be achieved by using the superposition principle.29,30 The incident beam is assumed to have the form 1inc(x~Y~z;P~4)
=
lof(x#P(x)G(Y)G(P
-
1)6(4)
7
(31)
Multiple
scattering
of a polarized
181
laser beam
wheref(x,y) is the incident spatial distribution of energy. Let all of the total (inclusive of all orders) diffuse Stokes vector field due to this spatial distribution be represented by I,(.u,y,z;~,~).The corresponding Stokes matrix for this case is M,(x,y,z;p,t$).If the searchlight loading is visualized as an impulse load, then the total (inclusive of all orders) diffuse Stokes matrix for this case, M(.ug,z;p,c#~),can be considered as the impulse response. Hence, by superposition, the multiply-scattered Stokes matrix for an incident distribution given by Eq. (31) is obtained as f(s’,f)M(x
- X’J - y’,q&)
dx’ dy’ .
(32)
The above equation is analagous to conoo1z+zg2ythe impulse response with the incident load distribution. 4.
MONTE
CARLO
ESTIMATOR
Before outlining the Monte Carlo estimator characteristics of the problem being presented, a brief outline of the basis behind the Monte Cario estimator used in this paper will be given and then the various estimators will be introduced. The goal of this section is to present estimators for the first-, second- and nth-order Stokes matrices, namely M,, M, and M,, respectively. In this section, the symbol (I) N implies the mean value of Z after N simulations. The symbol E[M] is the estimated value of the integral involved in the matrix M. Consider the integral
s
’wfx)g(x) dx ,
Q=
(33)
0
where w(x) is the weighting factor such that w(x) > 0
(34)
and
s
w(x) dx = 1 .
(35)
0
Thus w(x) has the characteristics of a probability density function” while g(x) is an arbitrary function. If x, is a random variable from a distribution w(x,), then an estimator of Q is given by
edPI = &xl)
(36)
and the mean value of the estimator after N events is
(Q>N = (l/W~eJQl. ,=I
Alternatively,
(37)
it is possible to sample xi from a uniform distribution, and define the estimator as
edQ1= dO+fx,) .
(38)
The latter estimator [i.e., Eq. (38)] will, in general, lead to a larger variance in the solution when compared with use of the former estimator [i.e., Eq. (36)]. Thus the selection of Eq. (36) to represent the estimator will result in a smaller value and is therefore termed a variance reduction technique.” Consider the cumulative distribution of w(x)
u = F(x)
=
s
wb)dy
0
Olxll
(39)
182
A. Ambirajan and D. C. Look
and the inverse transformation x=F-‘(U).
(40)
Thus Eq. (33) can be rewritten, by changing variables, as Q= [g[F-‘(U)]dU.
(41)
If U is sampled from a uniform distribution, then this is equivalent to sampling x from the distribution w(x).” Thus the estimator has the form
4Ql = gF’- ‘(UJI9
(42)
where U, is the Chsample from a uniform distribution and is equivalent to the estimater given by Eq. (36). This is the basis for O’Brien’s26 assertion that most variance reduction techniques are essentially equivalent to integral transformations. This is also the reason for the rather lengthy transformations carried out in the derivation of the integral forms for M,, M, and M, in Ref. 18. Now, from Eq. (6) and the fact that pI, > 0, it can be seen that the element p,, of the scattering matrix [Eq. (5)] has the characteristics of a probability density function. It would thus seem reasonable to sample the direction of a scattered photon from this “distribution”. Remember that p,, is the scalar Mie phase function.23 In the stochastic simulation of scalar (unpolarized) photon trajectories, the direction of the photon after each scattering event is obtained by sampling the scalar phase function32 (p,,). The same procedure will be used in the present paper to simulate the trajectory of a polarized photon. In essence, the trajectory of the polarized photon will be identical to that of a scalar photon. However, the weights associated with the trajectory will be different from the scalar case.” Another algorithm for the simulation of polarized photon trajectories, proposed recently, is the method of symmetrized trajectories. 33This method apparently produces reduced variances for problems involving spheroidal shell atmospheres. It is not used here. For the second order of scatter, the estimators are obtained based on Eqns (24) and (26). Thus the estimator for the double scattered Stokes matrix when ,u > 0 is e,[M,(x,y,z;p&)]
= (z
)I ‘~1~~“’
[arcsin(y2)
- arcsin(
where tl and p are sampled from uniform distributions. Note that x,, y,, 4,, y, and yz are obtained from Eqs. (28). Similarly, for ~1< 0 the following estimator is obtained 2e-:_e
1 ILi’L-““’ [arcsin
- arcsin(
exp[ - js
JG] . (44)
Since the direction of the sampled photon is obtained from p,,, Eq. 29a indicates that for p > 0
where Ii/b_, is obtained by selecting a random number from a uniform distribution and x;,_, and JJ;_, are derived from I&-, by means of Eqs. (30). pb _, and #, _, are obtained by sampling the phase function distribution pII by using an algorithm suggested by Barkstrom.-‘? The process of dividing the estimator in Eq. (45) by p,, is a consequence of the fact that the sampling is obtained
183
Multiple scattering of a polarized laser beam
from this (p,,) distribution. estimator for p < 0 is
This process is often called Importance
Sampling.3’ Similarly, the
(46) In Eqns (45) and (46) the angle 0 is given by cos 0 = /yin_, - (1 -/1)“2(1 -p,-,)‘~cos(~
- f#&_,).
(47)
Thus the estimator for higher order scatters is derived recursively from lower order scatters. 5. RESULTS In this section, the theory outlined above will be used to analyze the multiple scattering of a narrow, circularly polarized light beam by a plane-parallel medium composed of spherical dielectric scatterers. The backscattered component of the radiation will be studied. The spherical scatterers were assumed to have a relative refractive index of m = 1.197 + O.OOOi,which corresponds to latex particles suspended in water. A variety of situations were studied as shown in Table 1. For each size parameter, the backscattered Stokes vector was obtained for four optical thicknesses (,C) of the slab. The various elements of the Mie phase matrix are shown plotted in Fig. 2 for the size parameters of interest in this study. Further, all results presented in this paper are for an albedo (w) of one. Results for smaller albedos can easily be obtained since the intensity at each order of scatter is stored.” It should be noted that the numerical technique used computes the entire backscatter reflection matrix (M) of the medium. As pointed out earlier, M operates on the incident Stokes vector (IJ as shown in Eq. (13). If the incident state of polarization is circular, the diffuse radiation field is symetric about the laser beam.“4 In this study, the analysis will focus on the case of a circularly polarized searchlight beam incident on the medium. The Stokes vector of the incident light beam is
(48)
where 1, is the intensity of the incident beam. The aim of this study is to examine the physical nature of the scattering of polarized light and to see how multiple scattering affects the degree of polarization. Knowledge of this fact would aid in determining the validity of using scalar theories in analyzing the multiple scattering of polarized light. The quantities of interest are the elements of the Stokes vector, the degree of linear polarization (P,), and the degree of circular polarization (PC). Another quantity of interest is the mean number of scattering events (I) undergone by a photon. Since Z,(x,y,z;p,~$) is the intensity (note that this represents the first element of the Stokes vector) of light which has undergone n scattering events, Table I. Range of values of the size parameter (x) and optical thickness (L). The number of orders of scatter are the numbers in the body of the table. L Y
0.430 I .592 2.379 3.264
0.1
0.5
I.0
5.0
25 25 25 25
50 50 50 50
50 50 50 50
100 100 100 100
A. Ambirajan
184
and D. C. Look
and X:= ,Z,(x,y,z;&) is the intensity of light which has been scattered at least once (i.e., n 2 l), the probability that a photon gets scattered n times (where n 2 l), pnr is
Thus the mean number of scattering events [I(x,y,z;p,+)] undergone by a photon at a particular location of interest is the estimated value of n, and is evaluated with the equation
In the following sections, the results for multiple scattering of a circularly polarized laser beam will be presented. Comparisons with the scalar theory will also be presented. A FORTRAN code was written to implement the above theory and was executed on a HP-9000 work-station. The random numbers were generated by a multiplicative congruential generator with
(b) x=0.430 -+-~~=I.592 -x=2.379 -x=3.264
IO”
(Osqll (0.241~) (0.36Opj (0.494fi)
-
GO.430 ~I.592 x=2.379 ~~3.264
IO”
IO”
Optical radius, r
(0.065p (0.241~ (0.36Op (0.494p
lo”
IO”
Optical radius, r
(4
(cl ---cGo.430 --a-x=1.592 -x=2.379 ~~3.264
(O.Ot+) (0.241p) (0.3flOp) (0.494@
-t --o-c ---A---
Optical radius, r Fig. 5. Normalized
intensities
for a number
GO.430 x=1.592 x=2.379 ~~3.264
(0.065# (0.241~) (0.36Op) (0.494v)
Optical radius, r of size parameters
and optical
thicknesses.
18.5
Multiple scattering of a polarized laser beam
shuffling, and a multiplier of 16 807.35 For each case considered in this paper, 1 400 000 photons were traced. The computer code used in this study was validated in two ways. The first method of validation was to assume a phase matrix such that only the p,, term was non-zero. Thus the intensity predicted by the code should yield results which are identical to the scalar case for the Z term. Further, the Q, U, and V terms should all be zero. This indeed was the case. Secondly, to check the validity of the polarized part of the code, the data were compared with the polarized one-dimensional code of Evans? This was done by convolving [see Eq. (32)] the polarized searchlight solution with a uniform distribution [f(x,y) = l] over the incident surface for the phase matrices of interest in this study. The results from the polarized searchlight code compared well with the one-dimensional data. 5.1. ZntensitJ
The radial (see Fig. 3) variation of the backscattered polarized intensity appears to have the same form as that observed for the scalar problem. 27.28 This can be seen in Fig. 5, where the normalized intensity increases and then decreases with an increase in the optical radius (r). Some interesting features can be noted. For increasing particle size parameters, the backscattered intensity decreases. This can be attributed to the fact that larger particles have a larger asymmetry factorz3 which implies that a greater portion of the energy of the incident beam is scattered in the forward direction. For larger optical depths [see Fig. 5(d)], the backscattered intensity increases because a larger proportion of the energy of the incident beam is backscattered.
. A v x
40 -
5
IA.5
x=0.430(0.065&L) x=1.592(0.241~) ~2.379 (0.36Op) ~3.264 (0.494~)
CJe1.0
I”=( 1.0,O.l)
z ee F “, ‘:
30.4
1 201, 1
=
A
*
A
A:.
n
: v
I.. IO-
,
‘X ov IO“
%
.
I
,xX y
Y”
x
;
-
v )i
+ 7,
v IO”
IO”
ld
Optical radius, r 50
. A V q
40 :
E %
30
.
..,,....,
. . . . . ...
. . . . .. IA.0
x=Oo.430 (0.065p) x=1.592(0.241~) x=2.379(0.36Q1) x=3.264(0.49+)
.,. 6Fl.O
I,=( l,O,O,I1
.
1
3
F 7 20‘:I A;:
;::
A ’
t’
b t
=
IO-
, w IO”
“XI .w . . - v lo”
I x ; .r...J IO’
7
‘I
_,, ,,
t
,:.
loo
Optical radius. r
Fig. 6. Percentage difference between the vector and scalar radiative transfer solutions for the diffuse light field.
186
A. Ambirajan and D. C. Look
-
-x=0).430 -x=1.592 -x=2.379 -x=3.264
x=0).430(0.065 x=1.592(0.241 x=2.379(0. x=3.264(0.
(0.065&l) (0.2418) (0.3m) (0.494N)
Optical radius, r
r t
0.8
0.6
L=l.a=l.O
1
Io=ll,o,o,l)
j
-c
-c x=0.430(0.065~) --a-.-x=1.592 (0.241~) --c x=2.379 (0.36@) -----b x=3.264(0.499)
-
x=0.430(o.Ofjsp) x=1.592(0.241~) x=2.379(0.36Op) x=3.264(0.494~)
d
0.0: lo”
IO”
IO”
Optical radius, r
IO0
lo”
IO”
I@’
IO”
Optical radius, r
Fig. 7. Degree of linear polarization for a number of size parameters and optical thicknesses.
In Fig. 6 the vector and scalar” estimates of the intensity for two optical thicknesses (L) are compared. Clearly, for smaller size parameters, the difference between the vector and scalar theories is significant whereas it is smaller for larger size parameters. Also, it appears that for very large optical radii the difference between vector and scalar solutions decreases. However, this does not necessarily imply that the diffuse light field at large optical radii or small size parameters is unpolarized. Quite the contrary is observed in the following presentation. Thus even though the vector and scalar solutions approach each other as far as the intensity is concerned, the diffuse light field may still be polarized. In other words, the Q, CJ, and V terms may be significant with respect to the intensity I. This indeed was the case in our numerical simulations. 5.2. Degree of linear polarization In Fig. 7 the degree of linear polarization for the above cases is outlined. The behavior here is very interesting. The degree of linear polarization at small optical radii is the same irrespective of optical depth. This is because at small optical depths, the multiple scattering is primarily second-order (see Fig. 8). When the optical radius is small (i.e., x’ + JJ* is small), the quantities y, and y2 in Eqs 23a and 23b tend to one and zero, respectively, when z = 0 (upper surface) and x2 + y2+0 (small optical radii). Thus the double scattered Stokes matrix for backscatter [Eq. (26)]
187
Multiple scattering of a polarized laser beam
becomes essentially independent of optical depth. However, the variation with respect to particle size is harder to explain, particularly at small radii. One clue to this variation is the sharp variations in the functions comprising the phase matrix for Mie scatterers as shown in Fig. 2. This type of oscilation with respect to particle size for side scatter in a dilute suspension was first observed experimentally by Look and Chen.j4 For an optical depth of 0.1, the degree of linear polarization at large optical radii exhibits a behavior that is monotonic with respect to particle size. The degree of linear polarization decreases with an increase in particle size. The case of x = 3.264, however, shows an anomalous behavior in the sense that the data are highly oscillatory at even small optical radii. This is because the intensity is very low in this region for x = 3.264. Thus when the degree of linear polarization is computed by dividing by the intensity, the errors in this quantity lead to unpredictable results. Similar trends are observed for larger optical thicknesses. However, the degree of linear polarization decreases for a given size parameter with an increase in the optical thickness. The interesting feature in Fig. 7 is due to the fact that, even at reasonably large optical radii, the light field exhibits a high degree of linear polarization. In Fig. 8, the mean number of scattering events (I) for a medium with an optical thickness of 5.0 at an optical radius of 1.0 is greater than 5. Thus even for quite a large number of scattering events (see Fig. 8) the light field retains its polarized character (Fig. 7).
25.0
I
,,
1
111..(
‘....‘,
,....,.,
,?-
M.1, 20.0 -
__ .-...----------
-
uF1.0
.
lo=1 1,O,O.l I
. ----...-. ---
GO.430 (0.065p) x=l.592(0.241~) x=2.379 (0.3fX&, x=3.264 (0.494P)
15.0 -
x=0 430 (O.Oql) n=I.J92(0.241p) x=2.379 (0.360)~) x=3.264 (0.494~)
15.0 -
L
L 10.0 -
. 10.0 -
5.0 -
0.0 lo”
, “...,J
, ‘...*d 10-l
I@*
0.0’
“Lei,J’ Id
“11,,1’
lo”
1V’
Optical radius, r
‘ll(l.f lo”
‘.J
‘l(((ii’ IO”
Optical radius, r
25.0
. . . . ..I
.,.,..,
v...,,
r
Ls5.0, co=l.O _ I,=(l.O.O,l) 20.0
15.0
0.0’ lo”
20.0 GO.430 (0.065~) ._------. x=1.592 (0.241~) -----.x=2.379 (0.3cXy~) .. x=3.264 (0.494~)
’
,....A lo”
.
..‘,,,I
......--. p
.
x=0.430 (0.065p) x=1.592 (0.241~) x=2.379 (0.3w) x=3.264 (0.4%)
IS.0 -
I
lo“
Optical radius, r
,..1,d loo
.
,I
I 10“
. . .,..I
. . . ,,,, I
IO“
loo
Optical radius, r
Fig. 8. Mean number of scattering events (r) for a number of size parameters and optical thicknesses.
188
A. Ambirajan
and D. C. Look
x=0.430(0. -x=1.592 (0. --c ~22.379(0. --b ~3.264 (0.
-
0.5
-----o-
-
x=0.430(0.065p) x=1.592(0.241~) x=2.379(0.3m) ~3.264 (0.494~)
.I..-n
lo”
lo”
IO’
lo” Optical radius, r
lo“
l@’
IO”
Optical radius, r
“OV] xd.430 -x=1.592 ~2.379 ~3.264
lo”
-t -
(O.O+&) (0.241~) (0.3q) (0.494~)
100
10”
Optical radius, I Fig. 9. Degree
of circular
polarization
x=0.430(O.OLqL) x=1.592(0.241~) x=2.379(0.36Op) x=3.264(0.4%)
lo”
10-l Optical radius,
for a number
of size parameters
IO”
I
and optical
thicknesses.
5.3. Degree of circular polarization
In Fig. 9, the degree of circular polarization is shown. In this case PCis identical for all optical depths at small optical radii. The anomalous behavior for a size parameter of s = 3.264 can be explained in a similar manner to the previous paragraph. At large optical radii, the trend is an increase in the degree of circular polarization for larger size parameters. This can be explained as follows: for small size parameters, the dz, term in the scattering matrix [see Eq. (5)] is negligible as shown in Fig. 2. This implies that the circular component of the vector radiative transfer equation is essentially decoupled.14 Thus the circular component (V) of the vector radiative transfer equation behaves as a scalar intensity with an effective albedo of less than one. It thus approaches zero much faster than the actual intensity (I). However, with an increase in the size parameter, the quantity d2, becomes more significant. The physical effect of d?, is the coupling of the circular and linear components of the light field, as pointed out by Zege and Chaikovskya.‘5 Thus, for a larger size parameter, an increase in the conversion of linear light to circular light can be expected. Another interesting feature to be observed is that the sign of the degree of circular polarization goes negative as the optical radius is increased. In other words, the direction of rotation of the electric vector changes at larger optical radii. However, the degree of circular polarization does not seem to decrease appreciably at large optical radii for increased optical depths of the medium. Again, the interesting physical effect is that the light field is polarized at large optical radii.
189
Multiple scattering of a polarized laser beam
5.4. Analysis for various orders of scatter In the above analysis, it is seen that the light field retains its polarized character at quite significant optical radii from the incident light beam. This is a very interesting phenomenon because one would expect that, at large optical radii from the beam, the increased number of scattering events would increasingly depolarize the light field. This does not seem to be the case. To understand this phenomenon further, the Stokes vector will be analyzed with respect to the order of scatter. Plots of some polarization parameters for two optical depths of the medium are shown in Figs 10 and 1I. The data in these figures are plotted against the order of scatter (N). The quantity I, [see Fig. IO(a) and Fig. 1l(a)] is the backscattered intensity for photons that have been scattered N times and the quantity ,/m [see Fig. IO(b) and Fig. 1l(b)] is the backscattered linearly polarized light that has been scattered N times. It can be observed that IN and dm decrease after N = 8 exponentially with respect to the order of scatter. A similar trend was observed for scalar multiple scattering in one-dimensional problems” and for scalar multiple scattering by a searchlight beam.‘” The linearly polarized component (JmV) and the on the same graph in Fig. 10(c) and Fig. 1l(c). The ratio of these two quantities (
0
’
“I’,
”
10
20
s ‘a
30
Order of scatter, N
1.01 L 0.9 :
I ,
Cd) I
I ,
k0.S ~1.0 x=2.379(0.36%) r&O01
0.8
loa
0.6
. .
lo”
. .
0.5 :
7 lo* r
b
. v
. ’
lo“
. .
101
1, w (Q;+U5’” l
.
0
’ 5 Order
.
8 ’ 10
of scatter, N
.
0.3 :
. .:
loa IV0
:
0.4 : . .
0.2 : . . v . . . . ’
.
_) . 8 . 3
.
15
0
. l m. ‘L”“““‘._ s
10
15
Order of scatter, N
Fig. IO. Various polarization parameters plotted as a function of the order of scatter for an optical depth of L = 0.5.
A. Ambirajan and D. C. Look
190
. . .,.
.
I,..
1
.I....,.
.
.
W .
L.&o ~1.0 x=2.379(0.360~)
Order of scatter, N
Order of scatter, N
(cl
1 ” ” I.. ” IA.0 o=l .O x=2.379(0.360~) r=O.Ool
100 . IO”
1.0.
0.9
7
(4 1
,
I
7
1
L=S.O~1.0 x=2.379(0.360~) PO.001
0.8 :
. 10-Z
0.7 y .
lo-’
..
I
. 10” lo” r lo1 r
. . .
.
. .
’ .
v.7.
8 ‘, . (Q;+U$‘”
0.6 ; 0.5 :
. . .
. S . . 1
lo” r
Order of scatter, N
Order of scatter, N
Fig. 11. Various polarization parameters plotted as a function of the order of scatter for an optical depth of L = 5.0.
I,,,) is P,N [see Fig. 10(d) and Fig. 1l(d)], which is the degree of linear polarization of a photon that has been scattered N times. For all the cases considered here, the incident light beam is circularly polarized and the intensities for higher orders converge for all optical radii. This was observed in an earlier paper by the present authors.‘* The same is observed here for the linearly polarized component. The curious feature illustrated in Fig. IO(d) and Fig. 1l(d) is that the degree of linear polarization initially decreases with an increase in the order of scatter and then actually increases after approximately the fifth order for both the optical thicknesses considered. This was observed for the other size parameters also. This would seem to go against intuition, which suggests that a higher number of scattering events should lead to a decrease in the polarization of the light field. It should be noted here that since most of the energy is concentrated at lower orders for smaller optical radii, the net (inclusive of all orders) light field will still be largely determined by the lower orders. The idea behind these plots, however, is to illustrate the fact that higher orders are polarized and would thus influence the degree of polarization of the diffuse light field. For a given order of scatter, the degree of linear polarization decreases with an increase in optical thickness. Another noteworthy point is that, with an increase in the optical radius, the higher orders of scatter gain
I91
Multiple scattering of a polarized laser beam
increasing significance in specifying the diffuse light field as shown in Fig. 10(a) and Fig. 1l(a). The increase in the polarization of higher order photons would tend to cause the the net light field to stay polarized to significant optical radii. The large scatter in the higher order terms is indicative of the decreased accuracy in the evaluation of these terms. 6. CONCLUSIONS
In this paper a backward Monte Carlo estimator was developed for the multiple scattering of a narrow, polarized light beam. It was used to analyze the multiple scattering of a circularly polarized light beam. Results from this study were compared with data in which polarization effects were ignored (scalar theory). Certain other polarization parameters were studied. Data were compared for a variety of size parameters (using Mie theory), optical thicknesses, and optical radii. It should be emphasized that the numerical data presented in this study are for the backscattered component of the diffuse light field. Significant differences were observed between the intensities computed with the vector radiative transfer theory and the scalar theory. Differences between the two theories as great as 25% were observed for the intensity and were greater for smaller optical thicknesses and smaller size parameters. The diffuse light field was polarized even at significant optical radii from the incident light source. The degree of polarization (both linear and circular) was strongly dependent on the size parameter. At smaller optical radii, where double scattering was predominant, the various polarization parameters were essentially independent of the optical thickness of the medium and depended only on the size parameters of the scatterers. The influence of the order of scatter on the various polarization parameters was presented. The intensity and other components of the Stokes vector were observed to decrease exponentially with respect to the order of scatter for the optical thicknesses considered here after N = 8. Further, it was observed that light that has been scattered more than five times is still polarized. Acknow[edgemPnrs-The authors of this paper wish to thank Dr F. K. Evans of the University making his one-dimensional nit.colarado.edu/ - evans).
polarized
radiative
transfer
code
freely
available
on
REFERENCES I. Kienle, A., Lilge, L., Patterson, M. S., Hibst, R., Steiner, R. and Wilson,
the
of Colorado for kindly world wide web (http:‘/
B. C., Appl. Opt., 1996, 35,
2304. 2. Wang, L. and Jacques, S. L., J. Opt. Sot. Am. A, 1993, 10, 1746. 3. Groenhuis, R. A. J., Ferwada, H. A. and Ten Bosch, J. J., Appl. Opt., 1983, 22, 2456. 4. Wilson, B. C. and Adams, G., Med. Phys., 1983, 10, 824.
5. Bissonnette, L. R., Bruscaglioni, P., Ismaelli, A., Zaccanti, G., Cohen, A., Benyahu, Y., Kleiman, M., Egert, S., Flesia, C., Schwendimann, P., Starkov, A. V., Noormohammadian, M., Oppel, I_. G., Winkler. D. M., Zege, E. P., Katsev, 1. L. and Polonsky, I. N., Appl. Phys., 1995, B60, 355. 6. Yamazaki, H., Tojo, K. and Miyanami, K., Powder Technol., 1992, 70, 93. 7. Carswell, A. I. and Pal, S. R., Appl. Opt., 1985, 24, 3464. 8. Look, D. C. Jr. Opt. Eng., 1989, 28, 160. 9. Look, D. C. Jr. and Chen, Y. R., J. Thermophys. Heat Transfer, 1993, 7, 631. 10. Look, D. C. Jr. and Chen, Y. R., Appl. Opt., 1995, 34, 144. II. Dogariu, M. and Asakura, T., J. Opt. (Paris), 1993, 24, 271. 12. Dogariu, M. and Asakura, T., Waves in Random Media, 1994, 4, 429. 13. Mueller, D. W. Jr. and Crosbie, A. L., submitted to J. &ant. Spectrosc. Radiat. Tram/k. 14. Mueller, D. W. Jr. and Crosbie, A. L., submitted to J. Quant. Spectrosc. Radiat. Transfer. 15. Zege, E. P. and Chaikovskaya, L. I., J. Quant. Spectrosc. Radiat. Tramfer, 1996, 55, 19. 16. Collins, D. G., Bllttner, W. G., Wells, M. B. and Horak, H. G., Appl. Opt., 1972, 11, 2684. 17. Marchuk, G. I., Mikhailov, G. A., Nazaraliev, M. A., Darbinjan, R. A., Kargin, B. A. and Elepov, B. S., The Monte Carlo Methods in Atmospheric Optics. Springer-Verlag, Berlin, 1980. 18. Ambirajan, A. and Look, D. C. Jr. J. Quant. Spectrosc. Radiat. Transfer, 1996, 56, 317. 19. Ishimaru, A., Part. Part. Syst. Charact., 1994, 11, 183. 20. Ishimaru, A., Wave Propagation and Scattering in Random Media, Vol. I. Academic Press, New York, 1978. 21. Hovenier, J. W. and van der Mee, C. V. M., Astron. Astrophys., 1983, 128, I. 22. Voss, K. J. and Fry, E. S., Appl. Opt., 1984, 23, 4427. 23. van de Hulst, H. C., Light Scattering by Small Particles. Dover, New York, 1981. 24. Dave, J. V., Appl. Opt., 1970, 9, 1888.
192
A. Ambirajan and D. C. Look
Wiscombe, W. J., Appl. Opt., 1980, 19, 1505. O’Brien, D. M., J. Quant. Spectrosc. Radiat. Transfer; 1992, 48, 41. Crosbie, A. L. and Dougherty, R. L., J. Quant. Spectrosc. Radiat. Transfer, 1982, 27, 149. Crosbie, A. L. and Dougherty, R. L., J. Q uant. Spectrosc. Radiat. Transfer, 1982, 28, 233. Hecht, E., Optics, 2nd edn. Addison-Wesley Publishing Company, New York, 1993. Spanier, J. and Gelbard, E. M., Monte Carlo Principles and Neutron Transport Problems. Addison-Wesley Publishing Co., Reading MA, 1969. 31. Kales, M. H. and Whitlock, P. A., Monte Carlo Methods, Vol. I: Basics. John Wiley & Sons, New York, 1986. 32. Barkstrom, B. R., J. Quant. Spectrosc. Radiat. Transfer, 1992, 53, 23. 33. Voshchinnikov, N. V. and Karjukin, V. V., Astron. Astrophys., 1994, 288, 1. 34. Look, D. C. Jr. & Chen, Y. R., paper presented at 6th AIAA/ASME Joint Thermophysics and Heat 25. 26. 27. 28. 29. 30.
Transfer Conference, Colarado Springs, Colarado,
1994. AIAA Paper No. 94-2093.
Users Manual, Version 1.1, Vol. 3. Houston, January 1989, Ch. 8. 35. IMSL MATH/LIBRARY, 36. Evans, F. K. and Stephens, G. L., J. Quant. Spectrosc. Radiat. Transfer, 1991, 46, 413. 37. van de Hulst, H. C., Multiple Light Scattering. Academic Press, New York, 1980.
APPENDIX
A
In this appendix formulae for the rotation angles used in the radiative transfer of polarized light will be given. Note that 0 is the included angle between the incident and scattered light beams and that cos 0, and cos 19~are represented as p, and p>, respectively. The formula for cos 0 can be obtained by using trigonometry and is given by cos 0 = cos 8, cos 92 - (1 - co? &)‘:‘(I - co? e,)’ 2cos(fjz- 4,)
(AlI
By using spherical trigonometry, the following expressions can be obtained cos u, =
cos
e2 - cos 8, cos 0 sin 8, sin 0 ’
cos u2= cos 8, -
cos
e2 cos 0
sin e2 sin 0
sin u, =
sin e2 sin(@, - &) sin 0
sin (r2=
sin 8, sin(f$, - &) sin 0
’
64’4 643)
Thus we can obtain cos 2a,, cos 20, sin 2a, and sin 2~~ by using the standard trigonometric identities cos 20 = 2 cos’ 0 - I
W)
sin 2a = 2 sin D cos u
647)
and Further, the rotation matrix L(u) is given by I L(o) =
; 0
0
0
0
sin 2u 0 cos 2u _ sin 2u cos 20 0 0 0 1I
W9