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

2.6 Multiple Particles

129

where i1 = [cN1,µ, dN1,µ]T, i2 = [cN2,µ, dN2,µ]T, and as usually, e = [aν , bν ]T is the vector containing the expansion coe cients of the incident field in the global

coordinate system. Further, defining the scattered field coe cients

s1 = Q111 (ks, ki,1)i1 ,

s2 = Q112 (ks, ki,2)i2 ,

and introducing the individual transition matrices,

T1 = −Q111 (ks, ki,1) Q311 (ks, ki,1) −1 ,

T2 = −Q112 (ks, ki,2) Q312 (ks, ki,2) −1 ,

we rewrite the matrix system (2.134) as

s1

−

S

rtr

= T 1

S

rt

T

1

10e ,

12 s2

s2

−

S

rtr

= T 2

S

rt

T

2

20e ,

21 s1

and find the solutions

rtr

T 2

rtr

−1

rt

+

rtr

T 2

rt

s1 = T 1 I

12

21

T 1

S

10

12

S

20

,

− S

S

S

rtr

rtr

T 2

−1

rt

+

rtr

T 1

rt

(2.135)

s2 = T 2 I

21

T 1 12

S

20

21

10

.

− S

S

S

S

To compute the T matrix of the two-particles system and to derive a scatteredfield expansion centered at the origin O of the global coordinate system we use the Stratton–Chu representation theorem for the scattered field Es in Ds. In the exterior of a sphere enclosing the particles, the expansion of the approximate scattered field ENs in terms of radiating vector spherical wave functions reads as

N

ENs (r) = fνN M 3ν (ksr) + gνN N 3ν (ksr) ,

ν=1

where the expansion coe cients are given by

$ fνN %

jks2

N

$ N

1

(ksr ) %

ν

=

e

(r )

gνN

π S1

i,1

1

M

1

(ksr )

ν

$ M

1

(ksr ) %

µs

N

ν

+ j

hi,1

(r1)

1

dS (r1)

εs

N

(ksr )

ν

+

jks2

eN

(r ) $ N

1

(ksr ) %

ν

π

S2

i,2

2

M

1

(ksr )

ν

$ M

1

(ksr ) %

µs

N

ν

+ j

hi,2

(r2 )

1

dS (r2 ) .

(2.136)

εs

N

(ksr )

ν

130 2 Null-Field Method

Finally, using the addition theorem for the regular vector spherical wave functions

M

1

(k r

)

M

1

(k

r

)

ν

tr

s

µ

s

1

1

(ksr )

=

S01

N

1

(ksr

)

N

ν

µ

ν

µ

1

with

S01tr = T 11 (ksr01) R (α1, β1, γ1) ,

and

M

1

(k

r

)

M

1

(k r )

ν

s

tr

µ

s 2

1

(ksr )

=

S02

N

1

(ksr )

N

ν

µ

ν

µ

2

with

S02tr = T 11 (ksr02) R (α2, β2, γ2) ,

we obtain

s = S01tr Q111(ks, ki,1)i1 + S02tr Q211(ks, ki,2)i2
= S01tr s1 + S02tr s2 ,	(2.137)

where s = [fνN , gνN ]T is the vector containing the expansion coe cients of the scattered field in the global coordinate system. Combining (2.135) and (2.137), and using the identities S20rt (S10rt )−1 = S21rtr and S10rt (S20rt )−1 = S12rtr, yields [187]

T =

S

tr

− S

rtr

S

rtr

−1

S

rtr

S

rtr

S

rt

01T 1

I

T

2

T

1

2

21

10

12

21

I + 12 T

+

S

tr

I

− S

rtr

S

rtr

−1

I +

S

rtr

S

rtr

S

rt

02T 2

T 2

T 1

12

20

,

21

T 1 12

21

where the explicit expressions of the transformation matrices S21rtr are given by

S12rtr = R (−γ1, −β1, −α1) T 11 (ksr12) R (α2, β2, γ2) ,

(2.138)

and S12rtr

and

S21rtr = R (−γ2, −β2, −α2) T 11 (−ksr12) R (α1, β1, γ1) ,

respectively. Equation (2.138) gives the system transition matrix T in terms of the individual transition matrices T 1 and T 2, and the transformation matrices S and S. S and S involve translations of the regular and radiating vector spherical wave functions, respectively, and geometric constraints are introduced by the S matrices. Obviously, the geometric restrictions r12 > r2

and r12

2.6 Multiple Particles

131

> r1 introduced by the S matrices are fulfilled if the smallest circumscribing spheres of the particles do not overlap. The following feature of

equation (2.138) is apparent: if T 2 = 0, then T = S01tr T 1S10rt , and if T 1 = 0, then T = S02tr T 2S20rt , as it should. In general, T is a sum of two terms, each

of which is a modification of these limiting values. The various terms in a formal expansion of the inverses occurring in (2.138) can be considered as multiple-scattering contributions. The terms involving only T 1 and T 2 repre-

sent reflections at S1	and S2, respectively, the terms involving only T 1		rtr	T 2
		12
rtr		S

and T 2S21 T 1 represent consecutive reflections at S2 and S1, and S1 and S2, respectively, etc.

2.6.2 Formulation for a System with N Particles

The generalization of the T -matrix relation to a system with more than two constituents is straightforward. The system of matrix equations consists in the null-field equations in the interior of Sl

			N
	−			rtr		S	rt	N
sl		T l	S		sp = T l		l0 e for l = 1, 2, . . . ,		(2.139)
			lp

p=l

and the matrix equation corresponding to the scattered field representation

N
s = S0trl sl .	(2.140)

l=1

In practical computer calculations it is convenient to consider the global matrix A with block-matrix components

All = I , l = 1, 2, . . . , N ,

A	lp	−	S	rtr	,		l, p = 1, 2, . . . ,	N	,
		= T l lp				l = p ,

and to express the solution to the system of matrix equations (2.139) as

	$ N	%
sl =	AlpT pSprt0	e , l = 1, 2, . . . , N ,	(2.141)
	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.140), the system T -matrix becomes

NN

T = S0trl AlpT pSprt0

l=1 p=1

and this transition matrix can be used to compute the scattering characteristics for fixed or random orientations of the ensemble [153, 165].

132 2 Null-Field Method

In order to estimate the computer memory requirement it is important to know the dimensions of the matrices involved in the calculation. If Nrank(l) and Mrank(l) are the maximum expansion order and the number of azimuthal modes for the lth particle, then the dimension of the transition matrix T l is dim(T l) = 2Nmax(l) × 2Nmax(l), where

Nmax(l) = Nrank (l) + Mrank (l) [2Nrank(l) − Mrank(l) + 1] .

The dimension of the global matrix A is

N	N
dim (A) = 2	Nmax(l) × 2 Nmax(l) ,
l=1	l=1
since

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

If Nmax gives the dimension of the system T -matrix dim (T ) = 2Nmax × 2Nmax ,

we have

dim Sprt0 = 2Nmax(p) × 2Nmax , dim S0trl = 2Nmax × 2Nmax(l) .

In our computer code, the parameters Nmax and Nmax(l), l = 1, 2, . . . , N, are independent. Nmax(l) is given by the size parameter of the lth particle, while Nmax is given by the size parameter of a sphere centered at O and enclosing the particles.

2.6.3 Superposition T -matrix Method

For a system of N particles with αl = βl = γl = 0, l = 1, 2, . . . , N, the transformations of the vector spherical vector wave functions involve only the addition theorem under coordinate translations, i.e.,

S0trl = T 11 (ksr0l) ,

Sprt0 = T 11 (−ksr0p) ,

and

Srtr = T 31 (ksrlp)

lp

for l, p = 1, 2, . . . , N. Consequently, the matrix equation (2.139) takes the form

2.6 Multiple Particles

133

	N
sl = T l T 11 (−ksr0l) e + T 31 (ksrlp) sp ,			(2.142)
	p=l
while (2.140) becomes	N
	s = T 11 (ksr0l) sl .		(2.143)

l=1

The precedent equations can be written in explicit form by indicating the vector and matrix indices

N

11

31

[(sl)ν ] = [(Tl)νν ]

T

(−ksr0l) [eµ]

T

(ksrlp) (sp)µ

,

ν

ν

µ

µ

p=l

(2.144)

N

[sν ] = T

11

(ksr0l) (sl)µ

,

(2.145)

νµ

l=1

and can also be derived by using the so-called superposition T -matrix method [169]. The superposition T -matrix method reproduces the two cooperative effects characterizing aggregate scattering: interaction between particles and far-field interference [273]. The interaction e ect take into account that each particle is excited by the initial incident field and the fields scattered by all other particles, while the far-field interference is a result of the incident and scattered phase di erences from di erent particles. The superposition T -matrix method involves the following general steps:

1.The expansions of the incident and scattered fields in each particle coordinate system.

2.The representation of the field exciting a particle by a single-field expansion (which includes the incident field and the scattered fields from all other particles).

3.The solution of the transmission boundary-value problem for each particle.

4.The derivation of the scattered-field expansion in the global coordinate system.

The interaction e ect is taken into account in steps 2 and 3, where the scattered fields from other particles are transformed and included in the incident field on each particle, while the interference e ect is taken into account in steps 1 and 4 that address the incident and scattered path di erences, respectively. For a system of N spheres, the individual component T matrices are diagonal with the standard Lorenz–Mie coe cients along their main diagonal, and in this case, the superposition T -matrix method is also known as the multisphere separation of variables technique or the multisphere superposition method [21, 24, 29, 72, 150]. Solution of (2.144) have been obtained using direct matrix inversion, method of successive orders of scattering, conjugate

134 2 Null-Field Method

gradient methods and iterative approaches [74, 95, 198, 244, 271]. Originally all methods were implemented for spherical particles, but Xu [274] recently extended his computer code for axisymmetric particles.

In the following analysis, we derive (2.144) and (2.145) by using the superposition T -matrix method. The field exciting the lth particle can be expressed as

N

Eexc,l (rl) = Ee (rl) + Es,p (rl) ,

p=l

where Ee is the incident field and Es,p is the field scattered by the pth particle. In the null-field method, transformation rules between the expansion coe cients of the incident and scattered fields in di erent coordinate systems have been derived by using the integral representations for the expansion co- e cients. In the superposition T -matrix method, these transformations are obtained by using the series representations for the electromagnetic fields in di erent coordinate systems. For the external excitation, we consider the vector spherical wave expansion

Ee(r) = aµM 1µ (ksr) + bµN 1µ (ksr)

µ

and use the addition theorem

M µ1 (ksr) =		11 (ksr0l) M ν1 (ksrl) ,
N µ1 (ksr)	Tµν		N ν1 (ksrl)
to obtain
Ee (rl) = al,ν M ν1 (ksrl) + bl,ν N ν1 (ksrl)
ν
with
al,ν =		11	(ksr0l) aµ .
bl,ν		Tµν	bµ

Similarly, for the field scattered by the pth particle, we consider the series representation

Es,p (rp) = fp,µM 3µ (ksrp) + gp,µN 3µ (ksrp)

µ

and use the addition theorem

M 3µ (ksrp) N 3µ (ksrp)

=	31	(ksrpl) M ν1 (ksrl)	for rl < rpl ,
	Tµν	N ν1 (ksrl)

2.6 Multiple Particles

135

to derive

Es,p (rl) = flp,ν M 1ν (ksrl) + glp,ν N 1ν (ksrl)

ν

with

flp,ν =	31	(ksrpl) fp,µ .
glp,ν	Tµν	gp,µ

Thus, the field exciting the lth particle can be expressed in terms of regular vector spherical wave functions centered at the origin Ol

Eexc,l (rl) = al,ν M ν1		(ksrl) + bl,ν N ν1 (ksrl) ,
ν

where the expansion coe cients are given by

						N
[(el)		] =	11	(ksr0l)	[eµ] +		31		(sp)
	ν		Tµν				Tµν	(ksrpl)		µ
						p=l

and as usually, el = [al,ν , bl,ν ]T. Further, using the T -matrix equation sl = T lel, we obtain

N

[(sl)

] = [(Tl)

]

11

31

(ksrpl)

, (2.146)

ν

νν

(ksr0l) [eµ] +

Tµν

(sp)

µ

Tµν

p=l

and since (cf. (B.74) and (B.75))

Tν31µ (ksrlp) = Tµν31 (ksrpl)

and

Tν11µ (−ksr0l) = Tν11µ (ksr

we see that (2.144) and (2.146) coincide. The scattered field is a superposition

individual particles,

l0) = Tµν11 (ksr0l) ,

of fields that are scattered from the

N

Es(r) = Es,l (rl)

l=1

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

whence, using the addition theorem

136 2 Null-Field Method

M µ3 (ksrl)	= 33	(		ksr0l) M ν3 (ksr)					for	r > r0l ,
N µ3 (ksrl)	Tµν		−			N ν3 (ksr)
we derive the scattered-field expansion centered at O
Es(r) = fν M ν3 (ksr) + gν N ν3 (ksr)
	ν
with		N
					33	( ksr0l)				(2.148)
	[sν ] =				Tµν		(sl)	µ	.
						−

l=1

Since,

Tνµ11 (ksr0l) = Tµν11 (−ksr0l) = Tµν33 (−ksr0l) ,

we see that (2.145) and (2.148) are identical.

2.6.4 Formulation with Phase Shift Terms

For an ensemble of N particles with αl = βl = γl = 0, l = 1, 2, . . . , N, the system T -matrix is given by

N N
T = T 11 (ksr0l) AlpT pT 11 (−ksr0p) .	(2.149)

l=1 p=1

In (2.149), the translation matrices T 11(−ksr0p) and T 11(ksr0l) give the relations between the expansion coe cients of the incident and scattered fields in di erent coordinate systems. In the present analysis, these relations are derived by using the direct phase di erences between the electromagnetic fields in di erent coordinate systems.

For a vector plane wave of unit amplitude and wave vector ke, ke = ksek , we have

Ee (rl) = epolejksek ·r0l ejksek ·(r−r0l )

= ejksek ·r0l aν M 1ν (ksrl) + bν N 1ν (ksrl)

ν

=al,ν M 1ν (ksrl) + bl,ν N 1ν (ksrl)

ν

and therefore

al,ν		aν
bl,ν	= ejksek ·r0l	bν	.

Further, using the far-field representation of the field scattered by the lth particle,

2.6 Multiple Particles 137

E

(r

) =

ejksrl

!E

(e ) + O

1

"

rl

s∞,l

s,l

l

r

rl

and the approximation

ejksrl

=

ejksr e−jkser ·r0l

1 + O

1

,

rl

r

r

we see that the angular-dependent vector of scattering coe cients

N

s (er ) = e−jkser ·r0l sl

l=1

approximates the s vector in the far-field region, and that the angulardependent transition matrix

N N
T (er ) = ejks(ek ·r0p −er ·r0l )AlpT p	(2.150)

l=1 p=1

approximates the T matrix in the far-field region. For spherical particles, this method is known as the generalized multiparticle Mie-solution [272, 273].

Equation (2.150) shows that we have to impose

dim (T ) = dim AlpT p = dim (T lp) ,

where T lp = AlpT p. Assuming dim(T ) = 2Nmax × 2Nmax, we can set Nmax(p) = Nmax for all p = 1, 2, . . . , N, and in this case, dim(T ) = dim(T p).

Alternatively, we can set Nmax				= maxp Nmax(p)		and use the convention
(T	)	0 whenever ν > 2N			(l) and µ > 2Nmax(p). This formulation
	lp νµ ≡		max			11
avoids the computation of the translation matrices T11 and involves only the
inversion of the global matrix A. The dimensions of T						depend on translation

distances and are proportional to the overall dimension of the cluster, while the dimension of A is determined by the size parameters of the individual particles. Therefore, the generalized multiparticle Mie-solution is limited by the largest possible individual particle size in a cluster and not by the overall cluster size. However, the angular-dependent transition matrix can not be used in the analytical orientation-averaging procedure described in Sect. 1.5 and in view of our computer implementation, this method is not e ective for T -matrix calculations.

2.6.5 Recursive Aggregate T -matrix Algorithm

If the number of particles increases, the dimension of the global matrix A becomes excessively large. Wang and Chew [250, 251] proposed a recursive T - matrix algorithm, which computes the T matrix of a system of n components by using the transition matrices of the newly added q components and the

,Ncs p=1

138 2 Null-Field Method

T matrix of the previous system of n − q components. In this section we use the recursive T -matrix algorithm to analyze electromagnetic scattering by a system of identical particles randomly distributed inside an “imaginary” spherical surface.

Let us consider Ncs circumscribing spheres with radii Rcs(k), k = 1, 2, . . . , Ncs, in increasing order. Inside the sphere of radius Rcs(1) there are N(1) particles having their centers located at r(1)0l , l = 1, 2, . . . , N(1), and in each spherical shell bounded by the radii Rcs(k − 1) and Rcs(k) there are N(k) particles having their centers located at r(0kl ), l = 1, 2, . . . , N(k). Obviously, N = N(p) is the total number of particles and Rcs(Ncs) is the radius of the circumscribing sphere containing the particles. For simplicity, we assume that the particles are identical and have the same spatial orientation. The scattering geometry for the problem under examination is depicted in Fig. 2.7.

At the first iteration step, we compute the system T -matrix T (1) of all particles situated inside the sphere of radius Rcs(1). At the iteration step k, we compute the system T -matrix T (k) of all particles situated inside the sphere of radius Rcs(k) by considering a [N(k) + 1]-scatterer problem. In fact, the transition matrix T (k) is computed by using the system T -matrix T (k−1) of the previous ,kp=1−1 N(p) particles situated inside the sphere of radius Rcs(k − 1) and the individual transition matrix T p of all N(k) particles situated in the spherical shell between Rcs(k − 1) and Rcs(k). In this specific case, the blockmatrix components of the global matrix A are given by

All

= I ,

l = 1, 2, . . . , N(k) + 1 ,

Alp = −T pT 31

l = p ,

l, p = 1, 2, . . . , N(k) ,

ksrlp(k)

,

Al,N(k)+1

=

−

T

pT

31

k

r(k)

=

−

T

pT

31

k

r(k)

, l = 1, 2, . . . ,

N

(k) ,

s

l0

− s

0l

AN(k)+1,l = −T (k−1)T 31 ksr(0kl ) , l = 1, 2, . . . , N(k) .

Rcs(1)	Rcs(k-1)
O	rol	(k)
		O
	Rcs(Ncs)	R	(k)
			cs

Fig. 2.7. Illustration of the recursive T -matrix algorithm for many spheres