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

2.7 Composite Particles

139

If A = A−1, and Alp, l, p = 1, 2, . . . , N(k) + 1, are the block-matrix components of A, the recurrence relation for T -matrix calculation read as

T (k) =

N(k)+1,N(k)+1T (k−1) +

N(k)

N(k)+1,lT

11 k r(k)

A

pT

A

− s

0l

l=1

+

N(k)

11 k

N(k)

lpT

k

r(k)

l,N(k)+1T (k−1) +

pT

11

r(k)

.

T

s

0l A

A

− s

0p

l=1

p=1

The procedure is repeated until all N particles are exhausted. At each iteration step, only a [N(k) + 1]-scatterer problem needs to be solved, where N(k) usually is much smaller than N. Thus, we can keep the dimension of the problem manageable even when N is very large. The geometric constraint of the T -matrix method requires that the N(k) particles completely reside inside the spherical shell between Rcs(k − 1) and Rcs(k). However, for a large number of particles and small size parameters, numerical simulations certify the accuracy of the recursive scheme even if this geometric constraint is violated.

As mentioned before, the T 11(−ksr0l) matrices are used to translate the incident field from the global coordinate system to the lth particle coordinate system, while the T 11(ksr0l) matrices are used to translate the scattered fields from the basis of the lth particle to the global coordinate system. A recursive T -matrix algorithm using phase shift terms instead of translation matrices has been proposed by Auger and Stout [4].

2.7 Composite Particles

A composite particle consists of several nonenclosing parts, each characterized by arbitrary but constant values of electric permittivity and magnetic permeability. The null-field analysis of composite particles using the addition theorem for regular and radiating vector spherical wave functions is equivalent to the multiple scattering formalism [189]. The translation addition theorem for radiating vector spherical wave functions introduces geometric constraints which are not fulfilled for composite particles. Several alternative expressions for the transition matrix have been derived by Str¨om and Zheng [219] using the Q matrices for open surfaces (the interfaces between the di erent parts of a composite particle). In the present analysis we avoid the use of any local origin translation and consider a formalism based on closed-surface Q matrices.

2.7.1 General Formulation

To simplify our presentation we first consider a composite particle with two homogeneous parts as shown in Fig. 2.8. We assume that the surfaces S1c = S1 S12 and S2c = S2 S12 are star-shaped with respect to O1 and O2,

140 2 Null-Field Method

P

n1

M1 r1 r19

	r		O1		r29	S1	O		S1
	r2			r 9				1	n21
		O					O
	Di,1				S12
Ds							S12
	Di,2			r99	r199 S2
		O					O		S2
			2					2
				r299	M2				n12
									n2

Fig. 2.8. Geometry of a composite particle

respectively. The main reason for introducing the star-shapedness restriction is that we want to keep the individual constituents reasonably simple. A novel feature is the appearance of edges, which are allowed by the basic regularity assumptions of the T -matrix formalism [256]. A composite particle can be treated by the multiple scattering formalism under appropriate geometrical conditions. The boundary-value problem for the composite particle depicted in Fig. 2.8 has the following formulation.

Given the external excitation Ee, He as an entire solution to the Maxwell equations, find the scattered field Es, Hs and the internal fields Ei,1, Hi,1 and Ei,2, Hi,2 satisfying the Maxwell equations

× Es = jk0µsHs ,	× Hs = −jk0εsEs in	Ds ,		(2.151)
× Ei,1 = jk0µi,1Hi,1 ,	× Hi,1 = −jk0εi,1Ei,1		in Di,1 , (2.152)
and
× Ei,2 = jk0µi,2Hi,2 , × Hi,2 = −jk0εi,2Ei,2		in	Di,2 ,	(2.153)
the boundary conditions
n1 × Ei,1 − n1 × Es = n1 × Ee ,
n1 × Hi,1 − n1 × Hs = n1 × He				(2.154)
on S1,
n2 × Ei,2 − n2 × Es = n2 × Ee ,
n2 × Hi,2 − n2 × Hs = n2 × He				(2.155)

2.7 Composite Particles

141

on S2, and

n12 × Ei,1 + n21 × Ei,2 = 0 ,
n12 × Hi,1 + n21 × Hi,2 = 0	(2.156)

on S12, and the Silver–M¨uller radiation condition for the scattered field (2.3).

The Stratton–Chu representation theorem for the incident and scattered fields in Di, where Di = Di,1 Di,2 S12 = R3−Ds, together with the boundary conditions (2.154) and (2.155), lead to the general null-field equation

× S1 [ei,1 (r ) − ee (r )] g (ks, r, r ) dS (r )
+	j
		× × S1 [hi,1 (r ) − he (r )] g (ks, r, r ) dS (r )
	k0εs
+ × S2			[ei,2 (r ) − ee (r )] g (ks, r, r ) dS (r )	(2.157)
+	j
		× × S2 [hi,2 (r ) − he (r )] g (ks, r, r ) dS (r ) = 0 , r Di .
	k0εs
Passing from the origin O to the origin O1, using the identities

g (ks, r, r ) = g (ks, r1, r1) , g (ks, r, r ) = g (ks, r1, r1 ) ,

restricting r1 to lie on a sphere enclosed in Di,1 and taking into account the integral form of the boundary conditions (2.156)

$ N

3

%

$ N

3

%

ν

ν

ei,1 ·

M

3

dS +

ei,2 ·

M

3

dS = 0 ,

S12

ν

S12

ν

$ M

3

%

$ M

3

%

ν

dS +

ν

dS = 0 ,

hi,1 ·

N

3

S12

hi,2 ·

N

3

S12

ν

ν

yields

jks2

#

$ N

3

(ksr1) %

ν

[ei,1 (r1) − ee (r1)] ·

3

(ksr )

π S

M

1c

ν

1

$%+

M

3

(ksr1)

µs

ν

dS (r

+ j

[h

(r

)

−

h

(r

)]

·

)

εs

i,1

1

e

1

N

3

(ksr )

1

ν

1

2

#

$

3

(ksr1 )

%

+

jks

[e

(r )

e (r )]

N

(2.158)

−

·

ν

π

S

i,2

2

e

1

M

3

(ksr )

2c

ν

1

$%+

M

3

(ksr1 )

µs

ν

+ j

[h

(r )

h (r )]

dS (r ) = 0 , ν = 1, 2, . . .

εs

i,2

2

− e

1

·

N

3

(ksr )

2

ν

1

142 2 Null-Field Method

Sr	P1
	r1

O1

S2c

Di,2	r199
	O2

M2

Fig. 2.9. Illustration of the auxiliary sphere in the exterior of S2c

The physical fact that the tangential fields are continuous across the intersecting surface is the key to the structure of the null-field equations, and it is not used explicitly latter on. To simplify the null-field equations (2.158) we use the Stratton–Chu representation theorem for the incident field in the interior and exterior of the closed surface S2c

$ %

−Ee(r1) 0

=

×

e

(r ) g (k

, r

, r ) dS(r )

S2c

e

1

s

1

1

1

+

j

h

(r ) g (k , r

, r ) dS(r ) ,

$

k0εs × ×

S2c

e

1

s 1

1

1

r1 Di,

2

%

.

r1 R3 − Di,2

Choosing an auxiliary sphere in the exterior of S2c as in Fig. 2.9, and restricting r1 to lie on this sphere, we derive

jks2					$ N	3			(ksr1 ) %
				ee (r1 ) ·		ν
π		S			M		3		(ksr )
			2c				ν
								1

$%

							M		3	(ksr1 )
		µs
									ν
+ j				h	(r )	·					dS (r ) = 0 , ν = 1, 2, ... ,
		εs		e	1		N	3		(ksr )	1
								ν
										1

whence, using the identities ee(r ) = ee(r1) and he(r ) = he(r1), we obtain

jks2		#
						[ei,1 (r1) − ee
π		S1c
				µs
+ j						[hi,1 (r1) − he
				εs

$%

N 3 (ksr1)

(r )] · ν

M 3ν (ksr1)

$%+

M 3 (ksr1)

(r )] · ν dS (r1)

N 3ν (ksr1)

+ jks2				e	(r )		$ N	3		(ksr1 ) %
						·		ν
				i,2	2				3		(ksr )
	π S						M
			2c						ν
											1

$%

M

3

(ksr1 )

µs

ν

dS (r ) = 0 ,

+ j

h

i,2

(r )

·

εs

2

N

3

(ksr )

2

ν

1

(2.159)

ν = 1, 2, . . .

2.7 Composite Particles

143

Similarly, passing from origin O to origin O2 and restricting r2 to lie on a sphere enclosed in Di,2, gives

jks2					$ N	3		(ksr2) %
				ei,1 (r1) ·		ν
π		S			M		3		(ksr	)
			1c				ν
								2

$%

M

3

(ksr2)

µs

ν

+ j

h

i,1

(r

)

·

dS (r )

εs

1

N

3

(ksr

)

1

ν

2

2

#

$

3

(ksr2 )

%

+

jks

[e

(r )

e (r )]

N

−

·

ν

2

S

i,2

2

e

M

3

(ksr )

2c

ν

2

$ M

3

(ksr2 ) %+

µs

ν

+ j

[h

(r )

−

h (r )]

·

dS (r ) = 0 ,

εs

i,2

2

e

N

3

(ksr )

2

ν

2

(2.160)

ν = 1, 2, . . .

The surface fields ei,1, hi,1 and ei,2, hi,2 are approximated by finite expansions in terms of regular vector spherical wave functions as in (2.132) and (2.133) respectively. Inserting these expansions into the null-field equations (2.159) and (2.160), yields the system of matrix equations

Q31

(k

, 1, k

, 1)i

1

+ Q31

(k

, 1, k

, 2)i

2

= Q31

(k

, 1, k , 0)e ,

1

s

i,1

2

s

i,2

1

s

s

Q31

(k

, 2, k

, 1)i

1

+ Q31

(k

, 2, k

, 2)i

2

= Q31

(k

, 2, k , 0)e , (2.161)

1

s

i,1

2

s

i,2

2

s

s

where the significance of the vectors i1, i2 and e is as in Sect. 2.6. The Q matrices involving vector spherical wave functions with arguments referring to di erent origins are given by (2.84)–(2.87).

To compute the scattered-field expansion centered at the origin O we use the Stratton–Chu representation theorem for the scattered field Es in Ds. The expressions of the scattered field coe cients are given by (2.136) with S1c and S2c in place of S1 and S2, respectively, and in matrix form, we have

s = Q11	(k , 0, k	i,1	, 1)i	1	+ Q11	(k , 0, k	i,2	, 2)i	2	.	(2.162)
1	s				2	s

The null-field equations (2.161) and the scattered field equation (2.162) completely describe the T -matrix calculation and correspond to equations (2.134) and (2.137) of Sect. 2.6.

2.7.2 Formulation for a Particle with N Constituents

For a composite particle with N constituents, the system of matrix equations consists in the null-field equations in the interior of Slc

N

Q31p (ks, l, ki,p, p)ip = Q31l (ks, l, ks, 0)e , for l = 1, 2, . . . , N (2.163)

p=1

and the matrix equation corresponding to the scattered field representation

N
s = Ql11(ks, 0, ki,l, l)il .	(2.164)

l=1

144 2 Null-Field Method

Considering the global matrix A with block-matrix components

	Alp = Qp31(ks, l, ki,p, p) ,	l, p = 1, 2, . . . , N ,	(2.165)
we express the solution to the system of matrix equations (2.163) as
	$ N	%
il =	AlpQp31(ks, p, ks, 0) e , l = 1, 2, . . . , N ,
	p=1

where A stay for A−1, and Alp, l, p = 1, 2, . . . , N, are the block-matrix components of A. In view of (2.164), the transition matrix of the composite particle becomes

NN

T = Q11l (ks, 0, ki,l, l)AlpQ31p (ks, p, ks, 0) .

l=1 p=1

For an axisymmetric composite particle, the transition matrix is blockdiagonal with each block matrix corresponding to a di erent azimuthal mode m. Since there is no coupling between the di erent m indices, each block matrix can be computed separately. To derive the dimension of the transition matrix, we consider an axisymmetric particle and set

Nmax(l) =	! Nrank(l) ,	m = 0	,
	Nrank(l) − |m| + 1 , m = 0

where Nrank(l) is the maximum expansion order of the region l and m is the azimuthal mode. Then, we have

dim Q31p = 2Nmax(p) × 2Nmax ,

dim Alp = 2Nmax(l) × 2Nmax(p) ,

dim Q11l = 2Nmax × 2Nmax(l),

and therefore

dim (T ) = 2Nmax × 2Nmax .

Thus, Nmax gives the dimension of the transition matrix and

	Nmax =	! Nrank ,	m = 0,
		Nrank − |m| + 1 , m = 0,

where Nrank is determined by the size parameter of the composite particle. It should be mentioned that in our computer code, Nrank and Nrank(l), l = 1, 2, . . . , N, are independent parameters controlling the T -matrix calculations.

The contributions of Zheng [276, 277], Zheng and Str¨om [279], and Zheng and Shao [278] demonstrated the applicability of the null-field method to particles with complex geometries and optical properties. The above formalism can be extended to nonaxisymmetric particles in which case, the domain of analysis can be discretized into many volume elements with di erent shapes and optical constants.

2.7 Composite Particles

145

2.7.3 Formulation with Discrete Sources

For a composite particle with two constituents, the null-field equations formulated in terms of distributed vector spherical wave functions

jks2 #		$	N	3		(ksr1) %
				ν
		[ei,1 (r1) − ee (r )] ·			3		(k r	)
π
	S1c		Mν				s 1

$%+

3

(k

r )

+ [e

(r

)

−

e

(r )]

·

Nν

s

1

dS (r

)

i,1

1

e

M

3

(ksr1)

1

ν

+

jks2

e

(r )

$ N

3

(ksr1 ) %

(2.166)

i,2

·

ν

π

S2c

2

3

(k r )

Mν

s

1

$%

3

(ksr1 )

µs

ν

+ j

h

i,2

(r )

·

M

dS (r ) = 0 , ν = 1, 2, . . .

εs

2

N

3

(ksr1 )

2

ν

and

jks2

$

3

(ksr2) %

ν

ei,1 (r1) ·

N

3

(k

r )

π

S1c

Mν

s

2

$%

3

(ksr2

)

µs

ν

+ j

h

i,1

(r

)

·

M

dS (r )

εs

1

N

3

(ksr2)

1

ν

2

#

$

3

(ksr2 )

%

+

jks

[e

(r )

e (r )]

N

(2.167)

−

·

ν

2

i,2

2

e

3

(k r )

S2c

Mν

s 2

$%+

3

(ksr2 )

µs

ν

+ j

[h

(r )

−

h

(r )]

·

M

dS (r ) = 0 , ν = 1, 2, . . .

εs

i,2

2

e

N

3

(ksr2 )

2

ν

are equivalent to the general null-field equation (2.157). The distributed vector spherical wave functions in (2.166)–(2.167) are defined as

1,3	1,3	[k(r1	− z1,nez )] ,
Mmn	(kr1) = M m,|m|+l
1,3	1,3	[k(r1	− z1,nez )] ,
Nmn	(kr1) = N m,|m|+l
and
1,3	1,3	[k(r2	− z2,nez )] ,
Mmn	(kr2) = M m,|m|+l
1,3	1,3	[k(r2	− z2,nez )] ,
Nmn	(kr2) = N m,|m|+l

where {z1,n}∞n=1 is a dense set of points situated on the z-axis and in the interior of S1c, while {z2,n}∞n=1 is a dense set of points situated on the z- axis and in the interior of S2c. In the null-field method with discrete sources,

146 2 Null-Field Method

the approximate surface fields eNi,1, hNi,1 and eNi,2, hNi,2 are expressed as linear combinations of distributed vector spherical wave functions as in (2.132) and

(2.133) but with M1µ(ki,1r1), Nµ1(ki,1r1) in place of M 1µ(ki,1r1), N 1µ(ki,1r1) and M1µ(ki,2r2 ), Nµ1(ki,2r2 ) in place of M 1µ(ki,2r2 ), N 1µ(ki,2r2 ), respectively.

For a composite particle with N constituents, the block-matrix elements

31

Alp are given by (2.165) with Qp in place of Q31p , and the expressions of

31

the elements of the Qp matrix are given by (2.84)–(2.87) with the localized vector spherical wave functions replaced by the distributed vector spherical wave functions. The transition matrix is

N N

11

31

lp

T =

Ql (ks, 0, ki,l, l)A

Qp (ks, p, ks, 0) ,

l=1 p=1

11

where the Ql

matrices contains as rows the vectors M

1

, N

1

and as columns

ν

ν

the vectors

µ1 ,

µ1

31

, while the Qp matrices contains as rows and columns

M3

N3

1

1

the vectors M

,

N

and M

µ, N

µ, respectively.

ν

ν

For axisymmetric particles, the z-axis of the particle coordinate system is the axis of rotation and the use of distributed vector spherical wave functions improves the numerical stability of the T -matrix calculations, especially for highly elongated and flattened composite particles [51].

2.8 Complex Particles

The transition matrix for complex particles consisting of a number of arbitrary constituents can be derived by combining the results established in the previous sections. As an example, we consider the particle geometry depicted in Fig. 2.10. The scatterer is a composite particle with two inhomogeneous constituents and each inhomogeneous part is a layered particle. We denote by i1 and i1 the vectors containing the expansion coe cients of the internal fields in Di,1, by i2 and i2 the vectors containing the expansion coe cients of the internal fields in Di,2 and by i3 and i4 the vectors containing the expansion coe cients of the internal fields in Di,3 and Di,4, respectively. i1 andi2 refer to the expansion coe cients corresponding to radiating vector spherical wave functions, while i1, i2, i3 and i4 refer to the expansion coe cients corresponding to regular vector spherical wave functions.

Taking into account the results established for composite particles, we see that the null-field equations in the interior of S1c = S1 S12 and S2c = S2 S12 are

Q331 (ks, 1, ki,1, 1)i1 + Q311 (ks, 1, ki,1, 1)i1

+Q332 (ks, 1, ki,2, 2)i2 + Q312 (ks, 1, ki,2, 2)i2 = Q311 (ks, 1, ks, 0)e ,

and

			2.8 Complex Particles		147
			Ds
	O3	i3	S1
	O1
	Di,3	S	3
	O
Di,1 S12		i1, i1
Di,2	O2		i2, i2

O4	i4	S2
Di,4
S4

Fig. 2.10. Geometry of a complex particle

Q331 (ks, 2, ki,1, 1)i1 + Q311 (ks, 2, ki,1, 1)i1

+Q332 (ks, 2, ki,2, 2)i2 + Q312 (ks, 2, ki,2, 2)i2 = Q312 (ks, 2, ks, 0)e ,

respectively. Using the results established for inhomogeneous particles, we deduce that the null-field equations in the exterior of S1c and the interior of S3 are

−i1 + Q113 (ki,1, 1, ki,3, 3)i3 = 0 ,

−Q311 (ki,1, 3, ki,1, 1)i1 + Q313 (ki,1, 3, ki,3, 3)i3 = 0 ,

while the null-field equations in the exterior of S2c and the interior of S4 take the forms

−i2 + Q114 (ki,2, 2, ki,4, 4)i4 = 0 ,

−Q312 (ki,2, 4, ki,2, 2)i2 + Q314 (ki,2, 4, ki,4, 4)i4 = 0 .

To derive the coe cients of the scattered-field expansion centered at O we proceed as in the case of composite particles and obtain

s = Q131 (ks, 0, ki,1, 1)i1 + Q111 (ks, 0, ki,1, 1)i1

+Q132 (ks, 0, ki,2, 2)i2 + Q112 (ks, 0, ki,2, 2)i2 .

The above equations form a system of seven equations with seven unknowns. Eliminating the vectors corresponding to the internal fields it is possible to obtain a relation between the scattered and incident field coe cients and thus, to derive the expression of the transition matrix.

The T -matrix description is also suitable for scattering configurations with several layered particles, Fig. 2.11, or inhomogeneous particles with several enclosures, Fig. 2.12. In the first case, we compute the transition matrices

148 2 Null-Field Method

S12	Ds		S22
O1		O2
Di,12		Di,22
S11	O	Di,21	S21
Di,11

Fig. 2.11. Geometry of two layered particles

Ds

S1

S2 S3

O2 O3

Di,2 O Di,3

Di,1

Fig. 2.12. Geometry of an inhomogeneous particle with two enclosures

of each layered particle with the formalism presented in Sect. 2.5, and then use the multiple scattering formalism to derive the transition matrix of the ensemble. In the second case, we compute the transition matrix of the twoparticle system, and then calculate the transition matrix of the inhomogeneous particle using the formalism presented in Sect. 2.4. In fact, we can consider any combination of separate and consecutively enclosing surfaces and for each case immediately write down the expression of the transition matrix according to the prescriptions given in the earlier sections.

2.9 E ective Medium Model

The optical properties of heterogeneous materials are of considerable interest in atmospheric science, astronomy and optical particle sizing. Random dispersions of inclusions in a homogeneous host medium can be equivalently considered as e ectively homogeneous after using various homogenization formalisms [36,210]. The e ective permittivity of a heterogeneous material relates