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

2.5 Layered Particles

119

we express i1 and i1 as

i1 = −A11e , i1 = −A21e ,

and use (2.100) to obtain

T = − Q131 (ks, 1, ki,1, 1) A11 + Q111 (ks, 1, ki,1, 1) A21 .

For axisymmetric layers and axial positions of the origins Ol (along the z-axis of rotation), the scattering problem decouples over the azimuthal modes and the transition matrix can be computed separately for each m. Specifically, for each layer l, we compute the Ql matrices and assemble these matrices into the global matrix A. The matrix A is inverted, and the blocks 11 and 21 of the inverse matrix are used for T -matrix calculation. Because A is a sparse matrix, appropriate LU–factorization routines (for sparse systems of equations) can be employed.

An important feature of this solution method is that the expansion orders of the surface field approximations can be di erent. To derive the dimension of the global matrix A, we consider an axisymmetric particle. If Nrank(l) is the maximum expansion order of the layer l and, for a given azimuthal mode m, 2Nmax(l) × 2Nmax(l) is the dimension of the corresponding Q matrices, where

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

then, the dimension of the global matrix A is given by

dim (A) = 2Nmax × 2Nmax ,

with

N−1

Nmax = Nmax (N) + 2 Nmax(l) .

l=1

The dimension and occupation of the matrix A is shown in Table 2.2 for three layers.

Since

dim A11 = dim A21 = dim Q131 = dim Q111 = 2Nmax(1) × 2Nmax(1) , it follows that

dim (T ) = 2Nmax(1) × 2Nmax(1) .

Thus, the dimension of the transition matrix is given by the maximum expansion order corresponding to the first layer, while the maximum expansion

120 2 Null-Field Method

Table 2.2. Occupation of the global matrix

	2Nmax(1)	2Nmax(1)	2Nmax(2)	2Nmax(2)	2Nmax(3)
2Nmax(1)	Q133	Q131	0	0	0
2Nmax(1)	−I(−Q113)	0	Q213	Q211	0
2Nmax(2)	0	−Q131	Q233	Q231	0
2Nmax(2)	0	0	−I(−Q213)	0	Q311
2Nmax(3)	0	0	0	−Q231	Q331

For distributed sources, the identity matrix I is replaced by the Q13 matrix

orders corresponding to the subsequent layers are in descending order. In contrast to these prescriptions, the solution method using a recurrence relation for T -matrix calculation requires all matrices to be of the same order, i.e.,

	dim (T ) = dim T l,l+1,...,		=		1	dim (	Q	l) = 2Nmax(1)	×	2Nmax(1) .
				2
		N

This requirement implies that the same number of basis functions must be used to approximate the surface fields on each layer. For concentrically layered spheres, this requirement is not problematic because the basis functions are orthogonal on spherical surfaces. For nonspherical layered particles, we approximate the surface fields by a complete system of vector functions and it is natural to use fewer basis functions for smaller layer surfaces. However, for convergence tests it is simpler to consider a single truncation index [181, 248].

2.5.3 Formulation with Discrete Sources

For a two-layered particle as shown in Fig. 2.2, the null-field equations formulated in terms of distributed vector spherical wave functions (compare to (2.70), (2.71) and (2.72))

jks2

#

$

3

(ksr1) %

ν

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

N

3

(k r )

π

S1

Mν

s 1

$

3

(ksr1) %+

µs

ν

+ j

[h

(r

)

h

(r

)]

M

dS (r

) = 0 ,

εs

i,1

1

−

e

1

·

N

3

(ksr1)

1

ν

jki2,1

$

1

(ki,1r1) %

ν

−

ei,1 (r1) ·

N

1

(k

r

)

π

S1

Mν

i,1

1

+ j

$ M

1

(ki,1r1) % dS (r

µi,1

(r

)

)

h

i,1

·

ν

εi,1

1

N

1

(ki,1r1)

1

ν

(2.113)

ν = 1, 2, . . . ,

2.5 Layered Particles

121

jki2,1

$

1

(ki,1r1 ) %

ν

+

ei,2 (r2 ) ·

N 1

(k

r )

(2.114)

π

S2

M

i,1

1

ν

$

1

(ki,1r1 ) %

µi,1

ν

+ j

h

i,2

(r )

·

M

dS (r ) = 0 ,

ν = 1, 2, . . . ,

εi,1

2

N

1

(ki,1r1 )

2

ν

and

jki2,1

$

3

(ki,1r2) %

ν

−

ei,1 (r1) ·

N

3

(k

r

)

π

S1

Mν

i,1

2

$%

+ j

M

3

(ki,1r2)

µi,1

(r

)

dS (r

)

h

i,1

·

ν

εi,1

1

N

3

(ki,1r2)

1

ν

jki2,1

$

3

(ki,1r2 ) %

ν

+

ei,2 (r2 ) ·

N 3 (k

r )

(2.115)

π

S2

M

ν

i,1

2

+ j

$ M

3

(ki,1r2 ) % dS (r ) = 0 , ν = 1, 2, . . .

µi,1

(r )

h

i,2

·

ν

εi,1

2

N

3

(ki,1r2 )

2

ν

are equivalent to the general null-field equations (2.68) and (2.69). The distributed vector spherical wave functions in (2.113)–(2.115) 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 S1, while {z2,n}∞n=1 is a dense set of points situated on the z-axis and in the interior of S2. Due to their completeness property, the distributed vector spherical wave functions can be used to approximate the surface fields as in

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

For a multilayered particle, it is apparent that the solution methods with distributed sources use essentially the same matrix equations as the solution methods with localized sources. The matrices A1, Al,l−1 and AN are given

pq pq

by (2.109), (2.110) and (2.112), respectively, with Ql in place of Ql , while the matrix Al−1,l is

l,l+1,...,N

122 2 Null-Field Method

		13
Al−1,l =		−Ql−1(ki,l−1, l − 1, ki,l−1, l − 1)	0		.
			31

0−Ql−1(ki,l−1, l, ki,l−1, l − 1)

pq

The expressions of the elements of the Ql 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

13		11	11	21
T = Q1 (ks, 1, ki,1, 1) A			+ Q1 (ks, 1, ki,1, 1) A
31
×Q1 (ks, 1, ks, 1) ,
11	13

where the Q1 (ks, 1, ki,1, 1) and Q1 (ks, 1, ki,1, 1) matrices contain as rows

the

vectors M

1

(ksr1),

N

1

(ksr1)

and as

columns the vectors

Mµ1

(ki,1r1),

ν

ν

(k

r

) and

3 (k

r

),

3(k

r

), respectively, while the Q

31

, 1)

1

i,1

i,1

i,1

(k , 1, k

Nµ

1

Mµ

1

Nµ

1

s

s

matrix contains as rows and columns the vectors M3ν (ksr1), Nν3(ksr1) and M 1µ(ksr1), N 1µ(ksr1), respectively.

The use of distributed vector spherical wave functions improves the numerical stability of the null-field method for highly elongated and flattened layered particles. Although the above formalism is valid for nonaxisymmetric particles, the method is most e ective for axisymmetric particles, in which case the z-axis of the particle coordinate system is the axis of rotation. Applications of the null-field method with distributed sources to axisymmetric layered spheroids with large aspect ratios have been given by Doicu and Wriedt [50].

2.5.4 Concentrically Layered Spheres

For a concentrically layered sphere, the precedent relations simplify to those obtained in the framework of the Lorenz–Mie theory. In this specific case, we use the recurrence relation (2.107) with Ql given by (2.106). All matrices are diagonal and denoting by (Tl,l+1,...,N )1n and (Tl,l+1,...,N )2n the elements of the matrix T , we rewrite the recurrence relation as

(T		l,l+1,...,N					)1
							n
	=		−		jn (xl) [mr,lxljn (mr,lxl)] rn − jn (mr,lxl) [xljn (xl)] pn					,
					(1)			(1)
					hn	(xl) [mr,lxljn (mr,lxl)] rn − jn (mr,lxl) xlhn			(xl) pn	(2.116)
(T	l,l+1,...,N					)2
						n
= −					jn (xl) [mr,lxljn (mr,lxl)] qn − mr2,ljn (mr,lxl) [xljn (xl)] sn						,
				hn(1) (xl) [mr,lxljn (mr,lxl)] qn − mr2,ljn (mr,lxl) xlhn(1) (xl) sn

(2.117)

2.5 Layered Particles

123

where

			= 1 + (T		)	1	hn(1) (mr,lxl)
	p	n		l+1,l+2,...,N				,
						n	jn (mr,lxl)
	s		= 1 + (T		)	2	hn(1) (mr,lxl)
		n		l+1,l+2,...,N				,
						n	jn (mr,lxl)
	r		= 1 + (T		)1		mr,lxlhn(1) (mr,lxl)		,
		n		l+1,l+2,...,N
							[mr,lxljn (mr,lxl)]
						n
	q		= 1 + (T		)2		mr,lxlhn(1) (mr,lxl)		,
		n		l+1,l+2,...,N
							[mr,lxljn (mr,lxl)]
						n

xl = ki,lRl = k0mlRl is the size parameter of the layer l and

ml

mr,l = ml−1

(2.118)

(2.119)

(2.120)

(2.121)

is the relative refractive index of the layer l with respect to the layer l − 1. To obtain a stable scheme for computing the T matrix, we express the above recurrence relation in terms of the logarithmic derivatives An and Bn:

n

mr,lAn (mr,lxl) rn +

pn jn (xl) − jn−1 (xl) pn

1

xl

(Tl,l+1,...,N )n

= −

,

n

(1)

(1)

mr,lAn (mr,lxl) rn +

pn hn

(xl) − hn−1

(xl) pn

xl

(2.122)

An (mr,l xl )

n

qn +

sn jn (xl) − jn−1 (xl) sn

2

mr,l

xl

(Tl,l+1,...,N )n

= −

, (2.123)

An (mr,l xl )

n

(1)

(1)

qn +

sn

hn

(xl) − hn−1 (xl) sn

mr,l

xl

where pn and sn are given by (2.118) and (2.119), respectively, and rn and qn are now given by

1 Bn (mr,lxl) h(1)n (mr,lxl)

rn = 1 + (Tl+1,l+2,...,N )n An (mr,lxl) jn (mr,lxl)

and

2 Bn (mr,lxl) h(1)n (mr,lxl)

qn = 1 + (Tl+1,l+2,...,N )n An (mr,lxl) jn (mr,lxl) ,

respectively.

The computation of Lorenz–Mie coe cients for concentrically layered spheres has been considered by Kerker [115], Toon and Ackerman [222], and Fuller [76], while recursive algorithms for multilayered spheres have been developed by Bhandari [13] and Mackowski et al. [154].

124 2 Null-Field Method

2.6 Multiple Particles

In this section we extend the null-field method to the case of an arbitrary number of particles by using the translation properties of the vector spherical wave functions. Taking into account the geometric restriction that the particles do not overlap in space, we derive the expression of the transition matrix for restricted values of translations. Our treatment closely follows the original derivation given by Peterson and Str¨om [187, 188].

2.6.1 General Formulation

For reasons of clarity of the presentation we first consider the generic case of two homogeneous particles immersed in a homogeneous medium with a relative permittivity εs and a relative permeability µs. The scattering geometry is depicted in Fig. 2.5. The surfaces S1 and S2 are defined with respect to the particle coordinate systems O1x1y1z1 and O2x2y2z2, respectively, while the coordinate system of the ensemble or the global coordinate system is denoted by Oxyz. The coordinate system O1x1y1z1 is obtained by translating the coordinate system Oxyz through r01 and by rotating the translated coordinate system through the Euler angles α1, β1 and γ1. Similarly, the coordinate system O2x2y2z2 is obtained by translating the coordinate system Oxyz through r02 and by rotating the translated coordinate system through the Euler angles α2, β2 and γ2. The main assumption of our analysis is that the smallest circumscribing spheres of the particles centered at O1 and O2, respectively, do not overlap. The boundary-value problem for the two scattering particles depicted in Fig. 2.5 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

			P
		r1		r2
				r
n1	M1			M2 n2
			r29	r199
	r19	r9		r99 r299
				r12

S1	O1		r01
		D
			i,1

Ds

r02	Di,2	O2	S2
O

Fig. 2.5. Geometry of two scattering particles

2.6 Multiple Particles

125

and Ei,2, Hi,2 satisfying the Maxwell equations
× Es = jk0µsHs ,	× Hs = −jk0εsEs in	Ds ,		(2.124)
× Ei,1 = jk0µi,1Hi,1 ,	× Hi,1 = −jk0εi,1Ei,1	in	Di,1 ,	(2.125)
and
× Ei,2 = jk0µi,2Hi,2 , quad × Hi,2 = −jk0εi,2Ei,2		in	Di,2 ,	(2.126)
the boundary conditions
n1 × Ei,1 − n1 × Es = n1 × Ee ,
n1 × Hi,1 − n1 × Hs = n1 × He				(2.127)
on S1 and
n2 × Ei,2 − n2 × Es = n2 × Ee ,
n2 × Hi,2 − n2 × Hs = n2 × He				(2.128)

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

The Stratton–Chu representation theorem for the scattered field Es in Di,1 and Di,2 together with the boundary conditions (2.127) and (2.128) yield the general null-field equation

Ee(r) + × S1 ei,1 (r ) g (ks, r, r ) dS (r )
+	j
		× × S1 hi,1 (r ) g (ks, r, r ) dS (r )
	k0εs
+ × S2 ei,2 (r ) g (ks, r, r ) dS (r )
+	j
		× × S2 hi,2 (r ) g (ks, r, r ) dS (r ) = 0 , r Di,1 Di,2 .
	k0εs

Before we derive the null-field equations, we seek to find a relation between the expansion coe cients of the incident field in the global coordinate system

Oxyz,

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

ν

and the expansion coe cients of the incident field in the particle coordinate system O1x1y1z1,

Ee (r1) = a1,ν M 1ν (ksr1) + b1,ν N 1ν (ksr1) .

ν

126 2 Null-Field Method

M

r1		r
		S
O1	r01	Ds
		O

Fig. 2.6. Auxiliary surface S

For this purpose we choose a su ciently large auxiliary surface S enclosing O and O1 (Fig. 2.6) and in each coordinate system we use the Stratton–Chu representation theorem for the incident field in the interior of S. We obtain

$ %

aν bν

and

$ %

a1,ν b1,ν

jks2

$ N

3

(ksr ) %

ν

= −

ee (r )

3

(ksr )

π

S

M

$

ν

%

µs

M

3

(ksr )

ν

dS (r ) ,

+ j

he (r )

3

(ksr )

εs

N

ν

=

jks2

e

(r

) $ N

3

(ksr1) %

−

ν

π

S

e

1

M

3

(ksr

)

ν

1

$%

					M		3		(ksr1)
	µs
							ν
+ j			he (r1)			3				dS (r1) ,
		εs			N				(ksr1	)
						ν

respectively. Using the addition theorem for radiating vector spherical wave functions

M

3

(k r

1

)

M

3

(k

r )

ν

s

rt

µ

s

3

(ksr )

=

S10

N

3

(ksr )

,

N

ν

µ

ν

1

µ

where

S10rt = R (−γ1, −β1, −α1) T 33

(−ksr01) ,

for

r > r01

and taking into account that S10rt is a block-symmetric matrix, yields

a1,ν

= S10rt

aµ

(2.129)

b1,ν

bµ

.

ν

µ

2.6 Multiple Particles

127

The condition r > r01 can always be satisfied in practice by an appropriate choice of the auxiliary surface S, whence, using the identity T 33(−ksr01) = T 11(−ksr01), we see that

S10rt = R (−γ1, −β1, −α1) T 11 (−ksr01) .

We proceed now to derive the set of null-field equations. Passing from the origin O to the origin O1, using the relations

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

and restricting r1 to lie on a sphere enclosed in Di,1, gives

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

$%

M

3

(ksr1)

µs

ν

dS (r

+ j

h

(r

)

·

)

εs

i,1

1

N

3

(ksr )

1

ν

1

+

jks2

e

(r )

$ N

3

(ksr1 ) %

(2.130)

i,2

·

ν

π

S

2

2

M

3

(ksr )

ν

1

$ M

3

(ksr1 ) %

$ a1,ν

%

µs

ν

+ j

h

(r )

·

dS (r ) =

−

, ν = 1, 2, ... ,

εs

i,2

2

N

3

(ksr )

2

b1,ν

ν

1

where the identities ei,2(r2 ) = ei,2(r1 ) and hi,2(r2 ) = hi,2(r1 ) have been used. For the general null-field equation in Di,2 we proceed analogously but restrict r2 to lie on a sphere enclosed in Di,2. We obtain

jks2

$ N

3

(ksr2) %

ν

ei,1 (r1) ·

3

(ksr )

π S

1

M

ν

2

$ M

3

(ksr2) %

µs

ν

+ j

hi,1

(r1) ·

3

(ksr )

dS

εs

N

ν

2

+

jks2

e

(r )

$ N

3

(ksr2 ) %

i,2

·

ν

π

S

2

2

M

3

(ksr )

ν

2

$ M

3

(ksr2 ) %

µs

ν

+ j

hi,2

(r2 ) ·

3

(ksr )

dS

εs

N

ν

2

(r1)

(2.131)

$%

(r ) =	−	a2,ν	, ν = 1, 2, ... ,
2		b2,ν

where, as before, we have taken into account that ei,1(r1) = ei,1(r2) and

hi,1(r1) = hi,1(r2).

µi,1

128 2 Null-Field Method

The surface fields ei,1, hi,1 and ei,2, hi,2 are the tangential components of the electric and magnetic fields in the domains Di,1 and Di,2, respectively, and the surface fields approximations can be expressed as linear combinations of regular vector spherical wave functions,

$

eNi,1(r1) hNi,1(r1)

% N
= c1N,µ	j
µ=1	−

+dN

1,µ −j

n1(r1) × M 1µ(ki,1r1)

εi,1 n1(r1) × N 1µ(ki,1r1)

n1(r1) × N 1µ(ki,1r1)

				(r	)		M 1 (k		r	) (2.132)
	εi,1 n					×
	µi,1	1		1			µ	i,1	1

and

$

eNi,2(r2 ) hNi,2(r2 )

% N
N
= µ=1 c2,µ −j
+d2N,µ
		j
	−

n2(r2 ) × M µ1 (ki,2r2 )

εi,2

n

(r )

×

N

1 (k

i,2

r )

µi,2

2

2

µ

2

n2(r2 ) × N µ1 (ki,2r2 )

(r )

M 1

(k

r ) . (2.133)

εi,2 n

×

µi,2

2

2

µ

i,2

2

Inserting (2.132) and (2.133) into (2.130) and (2.131), using the addition theorem for vector spherical wave functions

M

3

(ksr1 )

rtr

M

1

(k r

2

)

ν

µ

s

N

3

(ksr )

=

12

N

1

(ksr )

ν

µ

ν

1

µ

2

with

rtr =

(

γ

,

β ,

α )

31 (k

s

r

12

)

(α

, β , γ ) ,

for

r

< r

,

12

R −

1

−

1

−

1

T

R

2

2 2

2

12

S

and

M

3

M

1

(k

r

)

ν

(ksr2)

rtr

µ

s

1

=

,

3

21

1

N

(ksr

)

ν

µ

N

(ksr

)

ν

2

µ

1

with

S

rtr =

(

γ ,

−

β ,

−

α )

T

31 (

−

k r

)

R

(α , β , γ )

for

r

< r

,

R −

2

2

2

s 12

1 1 1

1

12

21

and taking into account the transformation rule for the incident field coe - cients (2.129), yields the system of matrix equations

	31		S	rtr	11	(ks, ki,2)i2 =	−S	rt
					Q2			10e ,
Q1 (ks, ki,1)i1 + 12
S	rtr	11			31		−S	rt
		Q1	(ks, ki,1)i1 + Q2			(ks, ki,2)i2 =		20e ,	(2.134)
21