Doicu A., Wriedt T., Eremin Y.A. Light scattering by systems of particles (OS 124, Springer, 2006

2.9 E ective Medium Model

149

the average electric displacement to the average electric field [137],

D = εe E ,

where D = D DdV , E = D E dV and D is the domain occupied by the heterogeneous material. The e ective-medium theories operate with various assumptions which determine their domain of applicability. The usual derivations of the e ective-medium approximations are based on electrostatic considerations and assume that the inclusions are much smaller than the wavelength and that the interactions between the inclusions can be neglected. The average electric displacement is expressed in terms of the average electric and polarization fields

D = εs E + P ,

while the average polarization field P is related to the dipole moment of a single inclusion p by the relation

P = n0p ,

where εs is the relative permittivity of the host medium and n0 is the number of particles per unit volume. Defining the exciting field Eexc as the local field which exists within a fictitious cavity that has the shape of an inclusion

1

Eexc = E + 3εs P

and assuming that the dipole moment can be expressed in terms of the exciting field as

p = αEexc

with α being the polarizability scalar, yields

εe = εs	3εs + 2n0α	.	(2.168)
	3εs − n0α

The main issue in the e ective-medium approximations is to relate the polarizability α to the relative permittivity εr, where εr = εi/εs, and εi is the relative permittivity of the inclusion. For spherical particles and at frequencies for which the inclusions can be considered very small, the Maxwell–Garnett formula uses the relation

α =	4	πR3εs			3 (εr − 1)			,	(2.169)
	3
						εr + 2
which yields			1 + 2c εr−1
εe = εs						εr+2	,		(2.170)
			1	−		c εr−1
						εr+2

150 2 Null-Field Method

where R is the radius of the spherical particles and c is the fractional volume (or the fractional concentration) of the particles, c = (4/3)πR3n0. Note that size-dependent Maxwell–Garnett formulas incorporating finite-size effects have been derived by using the Lorenz–Mie theory [52, 57], the volume integral equation method [129] and variational approaches [211]. However, in microscopic view, waves propagating in the composite do undergo multiple correlated scattering between the particles. When the concentration of the inclusions is high and the size of the inclusion becomes comparable to the wavelength, more complete models which takes into account the inclusion size, statistical correlations, and multiple-scattering e ects, are needed.

In this section, we consider the scattering of a vector plane wave incident on a half-space with randomly distributed particles. Our treatment follows the procedure described by Varadan et al. [236], Varadan and Varadan [235], and Bringi et al. [28] and concerns with the analysis of the coherent scattered field and in particular, of the e ective propagation constants in the case of normal incidence. This approach is based on the T -matrix method and on spherical statistics, even though the particles may be nonspherical. Lax’s quasi-crystalline approximation [138] is employed to truncate the hierarchy of equations relating the di erent orders of correlations between the particles, and the hole correction and the Percus–Yevick [186] approximation are used to model the pair distribution function. The integral equation of the quasicrystalline approximation gives rise to the generalized Lorentz–Lorenz law and the generalized Ewald–Oseen extinction theorem. The generalized Lorentz– Lorenz law consists of a homogeneous system of equations and the resulting dispersion relation gives the e ective propagation constant for the coherent wave. The generalized Ewald–Oseen extinction theorem is an inhomogeneous equation that relates the transmitted coherent field to the amplitude, polarization, and direction of propagation of the incident field.

Other contributions which are related to the present analysis includes the work of Foldy [69], Waterman and Truell [258], Fikioris and Waterman [65], Twerski [230, 231], Tsang and Kong [223, 224] and Tsang et al. [227]. The case of a plane electromagnetic wave obliquely incident on a half-space with densely distributed particles and the computation of the incoherent field with the distorted Born approximation have been considered Tsang and Kong [225, 226].

2.9.1 T -matrix Formulation

The scattering geometry is shown in Fig. 2.13. A linearly polarized plane wave is incident on a half-space with N identical particles centered at r0l, l = 1, 2, . . . , N, and each particle is assumed to have an imaginary circumscribing spherical shell of radius R. The particles have the same orientation and we choose αl = βl = γl = 0 for l = 1, 2, . . . , N. The relative permittivity and relative permeability of the homogeneous particles are εi and µi, and εs and µs are the material constants of the background medium (which also occupies

			2.9 E ective Medium Model 151
		P
	z
	rp		rl
Op	r
	rlp
			Ol
O1	r0p		R	ON
		r0l
ε i, µi			εs, µs
				y
x	O

ke eβ eα

Fig. 2.13. Geometry of identical, homogeneous particles distributed in the halfspace z ≥ 0

the domain z < 0). The vector plane wave is of unit amplitude:

Ee(r) = epolejke·r ,

and propagates along the z-axis of the global coordinate system, ke = ksez . Our analysis is entirely based on the superposition T -matrix method which

expresses the scattered field Es in terms of the fields Es,l scattered from each individual particle l (cf. (2.147)),

N

Es(r) = Es,l (rl)

l=1

N
= fl,µM µ3 (ksrl) + gl,µN µ3 (ksrl) .	(2.171)
l=1 µ

Taking into account that the particle coordinate systems have the same spatial orientation (no rotations are involved), we rewrite the expression of the scattered field coe cients sl, sl = [fl,µ, gl,µ]T, given by (2.142) as

	N
sl = T ejksz0l e +	T 31 (ksrlp) sp ,		(2.172)

p=l

where e is the vector of the incident field coe cients in the global coordinate system and T is the transition matrix of a particle.

152 2 Null-Field Method

Before going any further we make some remarks on statistically averaged wave fields. We define the probability density function of finding the first particle at r01, the second particle at r02, and so forth by p(r01, r02, . . . , r0N ). The probability density function can be written as [236]

p (r01, r02, . . . , r0N ) = p (r0l) p (r01, r02, . . . , , . . . , r0N | r0l)

= p (r0l, r0p) p (r01, r02, . . . , , . . . , , . . . , r0N | r0l, r0p) ,

where p(r0l) is the probability density of finding the particle l at r0l, p(r0l, r0p) is the joint density function of the particles l and p and the superscript indicates that the term is absent. From Bayes’ rule we have

p (r0l, r0p) = p (r0l) p (r0p | r0l) ,

where p(r0p | r0l) is the conditional probability of finding the particle p at r0p if the particle l is known to be at r0l. The conditional densities with one and two particles held fixed are related by the relation

p (r01, r02, . . . , , . . . , r0N | r0l) = p (r0p | r0l)

×p (r01, r02, . . . , , . . . , , . . . , r0N | r0l, r0p) .

The configurational average of a statistical quantity f is defined as

E{f } =	· · · f p (r01, r02, . . . , r0N )
D	D

×dV (r01) dV (r02) , . . . , dV (r0N ) ,

while the conditional averages with one and two particles held fixed are

El{f } =	· · · f p (r01, r02, . . . , , . . . , r0N | r0l)
D	D

×dV (r01) dV (r02) , . . . , , . . . , dV (r0N )

and

Elp{f } =	· · · f p (r01, r02, . . . , , . . . , , . . . , r0N | r0l, r0p)
D	D

×dV (r01) dV (r02) , . . . , , . . . , , . . . , dV (r0N ) ,

respectively.

Averaging (2.172) over the positions of all particles excepting the lth, we obtain the conditional configurational average of the scattered field coe cients

jk z	N	1			31
El {sl} = T e s	0l e +			Dp T		(ksrlp) Elp {sp} g (rlp) dV (r0p) ,
		V
	p=l

(2.173)

2.9 E ective Medium Model

153

where g is the pair distribution function of a statistically homogeneous medium,

p (r0p | r0l) =	1		g (r0p − r0l) =	1	g (rlp)
		V		V

and V is the volume accessible to the particles. In our analysis we assume that both N and V are very large, but the number of particles per unit volume n0,

n0 = N/V,

is finite. The integration domain Dp is the half-space z0p ≥ 0, excluding a spherical volume of radius 2R centered at r0l, i.e., |r0p −r0l| ≥ 2R. The above equation indicates that the conditional average with one particle fixed El{sl} is given in terms of the conditional average with two particles fixed Elp{sp}. In order to close the system we use the quasi-crystalline approximation

Elp {sp} = Ep {sp} for l = p ,

(2.174)

which implies that there is no correlation between the lth and pth particles other than there should be no interpenetration of any two particles.

To solve (2.173) we seek the plane wave solution

El {sl} = sejKsz0l ,

(2.175)

where s is the unknown vector of the scattered field coe cients and the effective wave vector Ke is assumed to be parallel to the wave vector of the incident field, i.e., Ke = Ksez . Substituting (2.174) and (2.175) into (2.173),

,N

and taking into account that for identical particles the sum p=l can be replaced by N − 1 and for a su ciently large N, n0 ≈ (N − 1)/V , we obtain

sejKsz0l = T ejksz0l e + n0 T 31 (ksrlp) g (rlp) sejKsz0p dV (r0p) .

Dp

(2.176) The integral in (2.176) can be written as (excepting the unknown vector s)

			Jl = J1,l + J2,l
where
	Dp T	31		jKsz0p
J1,l =			(ksrlp) e		dV (r0p) ,

J2,l = T 31 (ksrlp) [g (rlp) − 1] ejKsz0p dV (r0p) ,

Dp

and the rest of our analysis concerns with the computation of the matrix terms

J1,l and J2,l.

154 2 Null-Field Method

In view of (2.144), J1,l can be expressed as

A3 mn,

−

m n

(ksrlp) B3 mn,

−

m n

(ksrlp)

jK z

0p

J1,l =

−

−

e

s

dV (r0p) ,

B3

(k r

lp

) A3

(k r

lp

)

Dp

−mn,−m n

s

−mn,−m n

s

whence, taking into account the series representations for the translation co- e cients (cf. (B.72) and (B.73))

A−3 mn,−m n (ksrlp) =

2jn −n

nn (n + 1) (n + 1)

jn a (

−

m,

−

m

|

n , n, n ) u3

(k r

lp

) ,

×

n

1

m −mn

s

B−3 mn,−m n (ksrlp) =

2jn −n

nn (n + 1) (n + 1)

×

jn b

1

(

−

m,

−

m

|

n , n, n ) u3

(k r

lp

)

n

m −mn

s

with (cf. (B.68) and (B.69))

a1

(

−

m,

−

m

|

n , n, n ) = π mm π|m| (β) π|m | (β) + τ |m| (β) τ

|m |

(β)

0

n

n

n

n

×Pn|m −m| (cos β) sin βdβ ,

b1 (

−

m, m

|

n , n, n ) =

π mπ|m| (β) τ |m |

(β) + m τ |m| (β) π|m | (β)

−

−

0

n

n

n

n

×Pn|m −m| (cos β) sin βdβ ,

we see that the calculation of J1,l requires the computation of the integral

Imm1 n = ejKsz0p u3m −mn (ksrlp) dV (r0p) .

Dp

The procedure to calculate the integral Imm1 n is described in Appendix C and we have (cf. (C.2)),

16πR3

2n + 1

Imm1 n =

jn

ejKsz0l δmm Fn (Ks, ks, R)

2

2

2

(ksR)

− (KsR)

2π

2n + 1

+

jn +1

ejksz0l δmm

2

ks2 (Ks − ks)

		2.9 E ective Medium Model 155
with
	(1)
Fn (Ks, ks, R) = (ksR) hn		(2ksR) jn (2KsR)

− (KsR) h(1)n (2ksR) (jn (2KsR)) .

Relying on this result we find that the matrix J1,l can be decomposed into the following manner:

J1,l = J1,RejKsz0l + J1,z ejksz0l .

(2.177)

The elements of the matrix J1,R,

(J1,R)11−mn,−m n (J1,R)12−mn,−m n J1,R = (J1,R)21−mn,−m n (J1,R)22−mn,−m n ,

are given by

(J1,R)−11mn,−m n = (J1,R)−22mn,−m n

=

32πR3

jn −n

Rmnn11 δmm ,

(ksR)2 − (KsR)2

nn (n + 1) (n + 1)

(J1,R)−12mn,−m n = (J1,R)−21mn,−m n

= −

32πR3

jn −n

Rmnn12 δmm ,

(ksR)2 − (KsR)2

nn (n + 1) (n + 1)

where Rmnn11 and Rmnn12 are defined as

∞

11,12

2n

+ 1

1,2

n

Rmnn =

(−1)

Fn (Ks, ks, R) Imnn n ,

2

n =0

and

I1

= π m2π|m| (β) π|m| (β) + τ |m| (β) τ |m| (β) P

n

(cos β) sin βdβ ,

mnn n

0

n

n

n

n

I2	= m	π π|m|	(β) τ		|m|	(β) + τ |m| (β) π|m|		(β) P	n
mnn n		n			n	n	n
		0
The matrix J1,z is of the form:
		(			)11		( )12
		J1,z =		J1,z −mn,−m n			J1,z −mn,−m n
			(J1,z )−21mn,−m n				(J1,z )−22mn,−m n

(cos β) sin βdβ .

(2.178)

,

156 2 Null-Field Method

and the matrix elements are given by

(J1,z )−11mn,−m n = (J1,z )−22mn,−m n

=

4π

jn −n+1

1

δmm ,

mnn

ks2 (Ks

−

ks) nn (n + 1) (n + 1) S

(J1,z )−12mn,−m n = (J1,z )−21mn,−m n

=

4π

jn −n+1

2

δmm

,

mnn

−ks2 (Ks

−

ks) nn (n + 1) (n + 1) S

where the expressions of Smnn1 and Smnn2 are similar to those of R11mnn and R12mnn excepting the factor Fn ,

∞

1,2

2n

+ 1

1,2

n

Smnn =

(−1)

2

Imnn n .

n =0

For m = 1 we have

S11nn

= S−1 1nn = snn ,

(2.179)

S12nn

= −S−2 1nn = −snn ,

(2.180)

where the coe cient snn is given by (cf. (C.6))

snn =

(−1)n+n

.

(2.181)

n(n + 1) (2n + 1)

n (n

+ 1) (2n + 1)

4

To derive the expression of J2,l, we need to compute the integral

Imm2 n = Dp [g (rlp) − 1] ejKsz0p um3 −mn (ksrlp) dV (r0p)

and the result is (cf. (C.3))

2n + 1

Imm2 n = 32πR3jn

ejKsz0l δmm Gn (Ks, ks, R) ,

2

where

∞

(1)

2

Gn (Ks, ks, R) =

1

[g (2Rx) −

1] hn

(2ksRx) jn (2KsRx) x

dx .

The matrix term J2,l then becomes

J2,l = J2,RejKsz0l

(2.182)

2.9 E ective Medium Model

157

and the elements of the matrix J2,R,

(J2,R)11−mn,−m n (J2,R)12−mn,−m n J2,R = (J2,R)21−mn,−m n (J2,R)22−mn,−m n

are given by

(J2,R)11−mn,−m n = (J2,R)22−mn,−m n

= 64πR3

jn −n

Rmnn21

δmm ,

nn (n + 1) (n + 1)

(J2,R)−12mn,−m n = (J2,R)−21mn,−m n

= −64πR3

jn −n

Rmnn22 δmm ,

nn (n + 1) (n + 1)

where Rmnn21 and Rmnn22 are defined as

∞

2n + 1

1,2

21,22

n

Rmnn =

(−1)

Gn (Ks, ks, R) Imnn n .

2

n =0

Taking into account the expressions of the integral terms I

1,2

(cf. (2.178))

it is apparent that

mnn n

Rmnn11

= R−11mnn

,

Rmnn12 = −R−12mnn

(2.183)

and

Rmnn22 = −R−22mnn .

Rmnn21 = R−21mnn

,

(2.184)

Substituting (2.177) and (2.182) into (2.176) gives two types of terms. One type of terms has a exp(jksz0l) dependence and corresponds to waves traveling with the wave number of the incident wave, while the other type of terms has a exp(jKsz0l) dependence and corresponds to waves traveling with the wave number of the e ective medium. The terms with wave number ks should balance each other giving the generalized Ewald–Oseen extinction theorem,

n0J1,z s + e = 0 ,

(2.185)

while balancing the terms with wave number Ks gives the generalized Lorentz– Lorenz law

[I − n0T (J1,R + J2,R)] s = 0 .

(2.186)

For an ensemble of particles with di erent orientations we proceed analogously. If the configurational distribution and the orientation distribution are assumed to be statistically independent, the total probability density

158 2 Null-Field Method

function can be expressed as a product of two functions, one of which p(r01, r02, . . . , r0N ) describes the configurational distribution and the other of which p(Ω1, Ω2, . . . , ΩN ) describes the distribution of particles orientations, i.e.,

p(r01, r02, . . . , r0N ; Ω1, Ω2, . . . , ΩN )

=p (r01, r02, . . . , r0N ) p (Ω1, Ω2, . . . , ΩN ) ,

where Ωl is the set of Euler angles specifying the orientation of the lth particle and Ωl = (αl, βl, γl). Each of these probability density functions is normalized to unity and the problem of computing the configurational and orientation averages are separated. Assuming that the particle orientations are statistically independent

p (Ω1, Ω2, . . . , ΩN ) = p (Ω1) p (Ω2) ...p (ΩN ) ,

and that the individual orientation distribution is uniform

1

p (Ωl) = 8π2 , l = 1, 2, . . . , N ,

we average (2.172) over all configurations for which the lth particle is held fixed and over all orientations. We obtain

	jk z		0l ' (		N		1			31	(
		s						'
El {sl} = e			T	e +			V	Dp	T T (ksrlp) Ep (sp) g (rlp) dV (r0p) ,
					p=l
where . denotes the orientation average
			1
El {sl} =					Ω · · ·			D · · ·	D slp (r01, r02, . . . , , . . . , r0N | r0l)
		(8π2)N				Ω

×dV (r01) dV (r02) , . . . , , . . . , dV (r0N ) dΩ1dΩ2, . . . , dΩN ,

sl = sl(Ω1, Ω2, . . . , ΩN ), Ω is the unit sphere, dΩl = sin βldβldαldγl and

' (		1
T	=		Ω T (αl, βl, γl) dΩl
		8π2

is the orientation-averaged transition matrix. Making the assumption

'	31	( ' (	31	(ksrlp) Ep (sp) ,
T T		(ksrlp) Ep (sp) = T T

which implies that Ep{sp} depends only on Ωp [227], we obtain the generalized Lorentz–Lorenz law in terms of the orientation-averaged transition matrixT instead of the transition matrix T of a particle in a fixed orientation. We note that the submatrices of the orientation-averaged transition matrix are diagonal and their elements do not depend on the azimuthal indices (see Sect. 1.5).