Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.astro.spbu.ru/DOP/8-GLIB/op3/JQSRT_110_1356.pdf Äàòà èçìåíåíèÿ: Fri Nov 19 12:11:23 2010 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 01:08:10 2012 Êîäèðîâêà: ISO8859-5

ARTICLE IN PRESS
Journal of Quantitative Spectroscopy & Radiative Transfer 110 (20 09) 1356 -1368

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer
jo urnal h ome pag e: w ww.else vier.c om/locate/jqsrt

Separation of variables method for multilayered nonspherical particles
Alexander Vinokurov a,Ó, Victor Farafonov a, Vladimir Il'in
a b c

b,c

St. Petersburg University of Aerospace Instrumentation, Bol. Morskaya 67, St. Petersburg 190000, Russia Astronomical Institute, St. Petersburg State University, Universitetskij pr. 28, St. Petersburg 198504, Russia Pulkovo Astronomical Observatory, Pulkovskoe shosse 65/1, St. Petersburg 196140, Russia

article info
Article history: Received 29 December 20 08 Received in revised form 27 February 20 09 Accepted 27 February 20 09 PACS: 42.25Fx Keywords: Light scattering Separation of variables' method Nonspherical inhomogeneous scatterers

abstract
A separation of variables method based on expansions of the electromagnetic fields in terms of spherical wave functions is expanded at nonspherical (axisymmetric) particles with a rather large number of layers. Commonly used alternative approaches to systems of linear algebraic equations relative to unknown field expansion coefficients for layered particles are considered in some detail. The SVM code developed is compared with the EBCM, GMT and DDA codes designed for multilayered scatterers and some numerical results obtained for nonspherical scatterers with up to 10 0 layers are presented as illustrations. & 20 09 Elsevier Ltd. All rights reserved.

1. Introduction Light scattering by layered particles is of growing interest for modelling optical properties of disperse media studied in astrophysics, physics of atmosphere and oceans, ecology, biophysics as well as in various applications. Such particles are applied to represent natural scatterers of both quasilayered and more complex structures. A large number of methods have been used to treat layered particles, but this field has been never reviewed in sufficient detail. Therefore, a state of the art review is required. The separation of variables method (SVM) was most often applied to treat light scattering by layered scatterers (see [1,2]). There are numerous works devoted to layered spheres as the SVM solution (theory) of Mie [3] is easily extended to such spheres (see for details [4]). After the first paper of Aden and Kerker [5] on core-mantle spheres, there were dozens papers where this theory was refined and applied to different cases (see reviews in [6,7]). New works on the subject appear till now (see, e.g., [8-11]). The SVM using a cylindrical basis was utilized to obtain solutions for layered infinite cylinders which find some useful applications. Pioneer works were made on core-mantle circular cylinders by Kerker and Matijevich [12] and Shah [13]. Recent studies are presented in [14-16]. The SVM with a spheroidal basis was applied first to core-mantle spheroids [17] and later to multilayered ones (see a review [18] and recent papers [19-23]). Note that in all the works the particles had spheroidal layers with confocal

Ó Corresponding author. Tel.: +7 812 7832236; fax: +7 812 4287129.

E-mail address: alexander.a.vinokurov@gmail.com (A. Vinokurov). 0 022-4073/$ - see front matter & 20 09 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.20 09.02.031


ARTICLE IN PRESS
A. Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368 1357

boundaries. A breakthrough was made by Han et al. [24] who obtained an exact solution for spheroids with nonconfocal layers. Layered scatterers of arbitrary (but usually axisymmetric) shape were treated by the extended boundary condition method (EBCM) using a spherical basis and the T -matrix method (TMM). The classic works are a theoretical paper of ? Peterson and Strom [25] and a work of Wang and Barber [26] with numerical results. Many papers that apply this approach to layered scatterers are cited in comprehensive reviews of Mishchenko et al. [27-29]. We note two interesting modifications of the approach. Kyurkchan et al. [31] have extended their modification of the EBCM called pattern equation method to coated scatterers. The approach deals with scattering patterns representing the fields in the far field zone and is easily generalized on multilayered objects. Petrov et al. [30] have calculated the T -matrix for a scatterer with layered structure from a Sh-matrix that depends only on the shape of the layers (no dependence on their size parameters and refractive indices). The authors expand their approach at continuously varying layer parameters. Two methods mainly distinguished by use of different basis--the (generalized) point-matching method (PMM) and the generalized multipole technique (GMT) are close to the EBCM/TMM approach. Al-Rizzo and Tranquilla [32] have demonstrated that the generalized PMM can be efficiently applied to core-mantle axisymmetric scatterers. Doicu and Wriedt [33,34] have combined their null-field method with discrete sources (actually the GMT) with the TMM to treat layered three-dimensional scatterers of rather large eccentricity. Other methods--the coupled dipole method (CDM), the method of moments (MoM), the finite difference time domain method (FDTD), etc. [1,2] are more universal and hence much more computational time consuming than the SVM, EBCM, PMM and GMT when applied to simple shape scatterers. Therefore, the CDM, MoM, FDTD, etc. are mainly used for complicated objects. The discrete dipole approximation (DDA) is a typical representative of the CDM methods based on the volume integral formulation of the scattering problem. Core-mantle ellipsoids are involved in the standard DDA code of Draine and Flatau [35]. An example of application of the DDA to complex core-mantle particles can be found in [36]. The FDTD method was applied to quite different layered objects: antennas [37], biological tissues [39], human head [38], etc. This method is often used also to consider optical properties of layered media (see, e.g., [41]). Some other methods recently modified to treat layered particles are presented in [42-44] (see, e.g., [45] for references to earlier works). Besides the methods mentioned, some approximations can be also useful in study of layered scatterers. Multilayered ellipsoids in the Rayleigh (RA) and quasistatic (QSA) approximations were considered, for instance, in [46]. Layered spheres have been recently studied in the Rayleigh-Gans approximation in [47]. The anomalous diffraction (van de Hulst's) and geometrical optics approximations are applied to multilayered particles in [48]. Another way to approximate the optical properties of layered scatterers is the effective medium theory (EMT--see for more details [49]). A special rule of this theory for layered particles was described in [50]. The EMT is very efficient, but may be inaccurate in particular for layered scatterers (see discussion in [51]). Naturally, any of the approximations has an essentially limited accuracy level and applicability range. Three notes should be added. First, particles with an (large) inclusion can be considered as core-mantle ones. However, we skipped references to most papers on such particles as the main attention is paid here to multilayered scatterers. Second, it should be mentioned that there are many not works on scattering of acoustic waves by layered obstacles (e.g., [52,53]). Third, because of application demands essential efforts were directed at solution of the inverse problem for layered scatterers and some interesting results have been obtained (see, e.g., [54,55]). Numerous works recently done on wave scattering by layered particles reflects the fact that they provide a very useful model for various inhomogeneous natural scatterers. Though many methods were applied or specially developed to treat nonspherical particles with layers, it is still hardly possible to derive sufficiently accurate results when the number of (nonspheroidal) layers exceeds 5-10. In this paper, the SVM based on expansions of the fields in terms of spherical wave functions is applied to really multilayered nonspherical (axisymmetric) particles. In Section 2 main features of our approach (separation of the fields, special scalar potentials used, etc.) are described. Two alternatives used to solve systems of linear algebraic equations relative to the scattered field expansion coefficients are discussed in some detail. In Section 3 we consider convergence, accuracy and speed of our SVM approach. It is compared with the EBCM, GMT, DDA methods applied to multilayered scatterers. Some numerical results for multilayered scatterers are also presented. Main conclusions are drawn in Section 4.

2. Basic equations 2.1. Formulation of the problem We consider an L-layered scatterer embedded in a homogeneous medium. All layer boundaries Si ?i Ì 1; 2; ... ; Lî are assumed to be axisymmetric and to have the same symmetry axis. The boundaries divide the particle and the medium into domains Di ?i Ì 1; 2; ... ; L ? 1î. They are characterized by dielectric permittivity i and magnetic permeability mi pffiffiffiffiffiffiffiffi (see Fig. 1), the corresponding wave number is ki Ì i mi k0 , where k0 is the wave number in vacuum.


ARTICLE IN PRESS
1358 A . Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368

Fig. 1. An axisymmetric layered scatterer (z-axis coincides with the symmetry one) and notation used.

The electromagnetic fields in the domain Di are denoted by E?iî ; H ?iî. The fields satisfy the Maxwell equations 8 > E?iî ?r î Ì Ð 1 = Ò H ?iî ?r î; r 2 D ; > > i < ii k0 1 > ?iî > H ?r î Ì = Ò E?iî ?rî; r 2 Di ; > : imi k0 where i Ì 1; 2; ... ; L ? 1. The boundary conditions at all the layer surfaces are 8 < E?iî ?r îÒ ni ?r î Ì E?i?1î ?r îÒ ni ?r îjr2S ; i : H ?iî ?r îÒ ni ?r î Ì H
?i?1î

(1)

?r îÒ ni ?r îj

r 2S

i

;

(2)

where ni ?r î is the outer normal to the surface Si and i Ì 1; 2; .. . ; L. The fields E?iî ; H ?iî in the domains Di with ipL should be represented by sums of incoming and outgoing fields E
?iî

ÌE

?iî 1

? ? E2iî ;

H

?iî

ÌH

?iî 1

? ? H 2iî . ?L?1î ?L?1î 1 ?L?1î ?L?1î 1

(3) ÌE ,H ÌH . As incident plane wave can be ?1î ?1î 1 ; H 1 when considering the with the boundary conditions

In the innermost layer (a particle core) like inside a homogeneous particle one has just E the scattered field must satisfy the Sommerfeld radiation condition at infinity and an represented by incoming waves, we can denote the former by E?1î ; H ?1î and the latter by E 2 2 fields outside the particle. ?1î ?1î Thus, the problem is to find the unknown scattered field (E2 ; H 2 ) as a solution to Eqs. (1) (2) for the given incident plane wave ?E?1î ; H ?1î î. 1 1 2.2. Features of the approach used

To solve the problem we apply the SVM with a spherical basis. The following generalization of this method on nonspherical scatterers is used. All the fields are expanded in terms of spherical wave functions, the expansions are substituted in boundary conditions, and then integration of these conditions over the boundaries gives a system of linear equations relative to unknown scattered field expansion coefficients (see for more details, e.g., [56,57]). We extend this approach to layered scatterers and modify it by extracting an axisymmetric part of the fields and utilizing special scalar potentials for different parts of the fields. ? ? For each domain Di , the fields Esiî ; H siî (s Ì 1; 2; i Ì 1; 2; .. . ; L ? 1) are divided in two parts: E
?iî s

ÌE

?iî s;A

?E

?iî s;N

;

H

?iî s

ÌH

?iî s;A

?H

?iî s;N

,

(4)

on the azimuthal angle j and are hereafter called axisymmetric. Averaging of the where the fields angle should give zero. It is simple to prove that the light scattering problem can nonaxisymmetric be solved independently for axisymmetric ?EA Ì E1;A ? E2;A î and nonaxysimmetric ?EN Ì E1;N ? E2;N î parts of the fields [58]. Such representation of the fields has advantages and disadvantages discussed in [59] and can be easily skipped. Special scalar potentials are introduced for each of the field parts. For the axisymmetric parts, we employ p
?iî s ? Ì Esi;î ;j cos j; A ?iî s;A;j ? qsiî Ì H ?iî s;A;j

E ; H ?i;î do not depend sA ? ? parts Esi;î ; H si;î over this N N

?iî s;A

cos j,
?iî s;A ?iî s;A

(5)

,H are j-components of the fields E and H . Other components of the fields are derived from the where Maxwell equations. The potentials p?iî ; q?iî resemble the Abraham ones, but in contrast to them satisfy the corresponding Helmholtz equations [60]. It is possible to show that the boundary conditions for p?iî and q?iî can be separated, with p being related to the TE mode task and q to the TM mode one (see for more details [59]). A definition of the TE and TM modes is

? Esi;î ;j A


ARTICLE IN PRESS
A. Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368 1359

given, e.g., in [49]. Spherical coordinates ?r ; y; jî related with the axisymmetric particle are introduced in such a way that all layer surface equations take the form r Ì r i ?yî. In the case of the nonaxisymmetric parts, for the TE mode we use E
?iî s;N ? ? Ì = Ò?U siî iz ? V siî r î;

H

?iî s;N

Ì

1 ? = Ò = Ò?U ?iî iz ? V siî rî, s imi k0

(6)

and for the TM mode E
?iî s;N

Ì

1 ? ? = Ò = Ò?U siî iz ? V siî rî; ii k0

H

?iî s;N

? ? Ì = Ò?U siî iz ? V siî r î,

(7)

where iz is the unit vector along the particle symmetry axis. As the scalar potentials introduced satisfy Helmholtz equations they can be expanded in terms of spherical functions as follows: p
?iî 1 î 1 Xa l Ì1 ?iî 1;l ?iî l

? q1i

Ì

b1;

jl ?ki r îP1 ?cos yî cos j; l

U V

?iî 1 ?iî 1

Ì
?iî 2 ?iî 2

11 XX a mÌ1 lÌm

?iî 1;lm

b1

jl ?iî ;lm

?ki r îP m ?cos yî cos mj, l

(8)

p

?iî 2 î

? q2i

Ì

1 Xa

?1î h ?iî l lÌ1 b2;l

?iî 2;l

?ki r îP ?cos yî cos j;

1 l

U V

Ì

?iî 11 X X a2;lm mÌ1 lÌm

b2

?iî ;lm

hl ?ki r îP m ?cos yî cos mj, l

?1î

(9)

where hl ?ki r î and jl ?ki r î are the first kind Hankel and Bessel functions, P m ?cos yî the associated Legendre functions. So, to l ?1î ?1î î find the scattered field at any point r we need to determine the expansion coefficients a?1lî , b2;l and a?1lm , b2;lm . Naturally, the 2; 2; optical properties of a scatterer such as cross-sections, scattering matrix, etc. are expressed through these coefficients (see, e.g., [61]). One should realize that our selection of the scalar potentials U ; V is equivalent to use of the corresponding vector wave functions in the field expansions. For instance, for the TE mode and the electric field, instead of Eqs. (6) and (8) we could write E?iî ?r î Ì E 1
?iî 1; A

?1î

?r î?

11 XX ?a mÌ1 lÌm

?iî 1;lm

M z ?r î? b1 lm

?iî ;lm

M r ?r îî; lm

r 2 Di ,

(10)

where the vector wave functions are M z ?r î Ì = Ò?iz clm ?r îî; lm M r ?r î Ì = Ò?r clm ?r îî. lm (11)

Here clm ?r î is a solution to the scalar Helmholtz equation Dc ? ki c Ì 0. ?iî ? It is easy to prove that the coefficients asi;î ; bs;lm in Eqs. (8)-(10) are the same. Note that instead of our potentials U ; V lm one can use two usual Debye potentials V e;m and correspondingly expand the fields in the functions M r and lm N r Ì = Ò Mr . lm lm ? ? It should be noted that the axisymmetric task is actually scalar--one scalar potential (psiî or qsiî ) defines the fields ? Esi;î ; H ?i;î for each i; s. Hence solution of this task is simple and very fast. Our large experience of dealing with the A sA axisymmetric and nonaxisymmetric tasks indicates that their properties (convergence behaviour, accuracy, etc.) are very close. It makes the axisymmetric task very useful in extensive studies of applicability ranges of different approaches. Solution of this task also allows one quickly to derive proper values of such technical parameters as the number of knots in quadrature formula and of terms kept in the expansions necessary to reach given accuracy of results. A back side of our extracting the axisymmetric part of the fields is that generally we need to keep 1-2 additional terms in the field (potential) expansions to have as high same accuracy as that obtained without extracting. If necessary, one can skip this extraction simply by taking m Ì 0 in Eq. (10). 2.3. Linear systems to be solved The systems of linear algebraic equations relative to the scattered field expansion coefficients are derived as follows. We write the boundary conditions (2) separately for the axisymmetric (EA ; H A ) and nonaxisymmetric (EN ; H N ) parts of the fields and replace these field parts with the corresponding scalar potentials, using Eqs. (5)-(7). The potential expansions (8)-(9) are substituted in the reformulated boundary conditions. The conditions for each boundary are then multiplied by functions Pm ?cos yî sin my and integrated over the corresponding surface Si . Completeness of the spherical functions allows n one to get systems of linear algebraic equations relative to the expansion coefficients of the scattered field and all internal ones. In the axisymmetric task we consider the first N terms in the potential expansions and multiply the conditions written for p; q by N different spherical harmonics with m Ì 1. In the nonaxisymmetric task N 2 =2 terms are kept in the U ; V expansions and multiplication is made for N Ð m ? 1 spherical harmonics with m Ì 1; 2; ... ; N.


ARTICLE IN PRESS
1360 A . Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368

For an axisymmetric scatterer, one can solve the systems arising for different m independently, and hereafter we skip the index m. So, for each m one gets finite systems which look in the matrix form as follows: ! ! ! ! x?iî x?i?1î E?i?1î F ?i?1î A?iî B?iî Ì ; i Ì 1; 2; ... ; L. (12) y?iî y?i?1î G?i?1î H?i?1î C ?iî D?iî For each m the vectors x?iî respectively. The elements of and the associated Legendre Zp ?i jl ?ki r i ?yîîP m ? Amî;nl Ì l
0 ? ? and y?iî contain the expansion coefficients ?fa1i;îlm gN m ; fb1;lm gN m î and ?fa2i;îlm gN m ; fb2;lm gN m î, lÌ lÌ lÌ lÌ matrices A?iî ; B?iî ; C ?iî ; D?iî are integrals of products of the Bessel or first kind Hankel functions functions (and their first derivatives) over the surface Si , for instance ?iî ?iî

cos yîP m ?cos yî sin y dy, n

(13)

where r Ì r i ?yî is the innermost layer EL?1 2 ! ?Lî ?Lî A B C ?Lî D?Lî

surface equation of Si . Expressions for other matrix elements are presented in [62]. Note that in the Ì 0, i.e., y?L?1î Ì 0 and one has ! ! x?Lî E?L?1î Ì x?L?1î . (14) y?Lî G?L?1î

It should be pointed out that there are two alternative approaches to Eqs. (12)--one either composes one large system (see, e.g., [21,31,34] and most works on core-mantle particles) or avoids this by using recursive [25,26,33] or iterative [20,22,30,50] relations. 2.3.1. Single system approach It is easy to bring the known vector x?1î into the right-hand side of system (12) for i Ì 1 and to join partly related other systems for i Ì 2; 3; ... ; L. This gives the following system relative to the unknown coefficients of expansions of the scattered and internal fields: 10 ?1î 1 0 ?1î 1 0 y A ÐB?1î E?2î F ?2î 0 ÑÑ Ñ 0 CB ?2î C B C ?1î C B CB C H?2î 0 ÑÑ Ñ 0 CB x B ÐD?1î G?2î C CB ?2î C B B C B0C B B0 CB C 0 CB y ÐA?2î ÐB?2î E?3î ÑÑ Ñ CB B CB C C B ?3î C Ì B 0 Cx?1î . (15) B0 0 CB x ÐC ?2î ÐD?2î G?3î ÑÑ Ñ CB . C B . C B B . C B.C B Ñ ÑÑ ÑÑÑ ÑÑ Ñ ÑÑÑ Ñ ÑÑ ÑÑ Ñ CB . C B . C CB B CB C C B @0 0 0 0 ÑÑ Ñ E?L?1î AB y?Lî C B 0 C A@ A @ 0 0 0 0 ÑÑ Ñ G?L?1î x?L?1î 0 The dimension of this system is 2N L in the axisymmetric task and 4?N ? 1 Ð mîL for m Ì 1; 2; ... ; M in the nonaxisymmetric one. Usually, the number of azimuthal terms to be considered ?M î can be essentially smaller than the number of radial ones ?Nî. An advantage of this approach is flexibility in choosing N as its different values can be selected in layers. This may be useful for scatterers with the layer surfaces of different types as convergence of methods mainly depends on surfaces' shape. However, the approach meets an obvious problem for particles with a large number of layers as the system size and computational time quickly grow with an increase of L. This problem is solved by use of sparse matrix inversion algorithms [34]. 2.3.2. Iteration/recursion approach This approach to system (12) can be presented in different ways. For example, we introduce matrices P 1 ; P 2 so that ! ! ! P1 x?1î A?1î B?1î Ì x?L?1î . (16) P2 y?1î C ?1î D?1î The matrices P 1 ; i Ì L; L Ð 1; ... ! P1 Ì P2 P2 are easily derived by considering Eq. (12) with the inverted left-hand side matrix consequently for ; 2 and in its normal form for i Ì 1. This gives !2 !Ð1 !3 !Ð1 ! Y E?L?1î E?2î F ?2î 4LÐ1 A?iî B?iî E?i?1î F ?i?1î 5 A?Lî B?Lî . (17) G?2î H?2î G?L?1î G?i?1î H?i?1î C ?Lî D?Lî C ?iî D?iî iÌ2

We solve the system that is obtained from Eq. (16) ! ! ! ÐB?1î P 1 y?1î A?1î Ì x?1î . ÐD?1î P 2 x?L?1î C ?1î

(18)

Thus, here we make L Ð 1 inversions of 2N Ò 2N matrices in the axisymmetric part and 4?N ? 1 Ð mîÒ 4?N ? 1 Ð mî matrices for each m to be considered in the nonaxisymmetric one. So, the algorithm is more robust to an increase of the number of layers than the single matrix scheme presented above.


ARTICLE IN PRESS
A. Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368 1361

Note that one can rewrite Eq. (17) as follows: !2 ! !Ð1 3 ! L ?L?1î Y E?iî F ?iî P1 A?iî B?iî 4 5E Ì . ?iî P2 G?L?1î H?iî C ?iî D?iî i Ì2 G This rather iterative relation allows one easily to formulate a recursive relation for the matrices P 1 ; P ! ! ! !Ð1 P 1 ?Lî P1 ?L Ð 1î E?2î F ?2î A?2î B?2î Ì , P 2 ?Lî P2 ?L Ð 1î G?2î H?2î C ?2î D?2î
2

(19)

(20)

where Ps ?Lî means the matrix P s for a particle with L layers (s Ì 1; 2) and the innermost layer of an L-layered particle is characterized by the matrices A?2î ; B?2î ; ... ; H?2î . Instead of solution of a linear system one often calculates a T -matrix that relates the unknown and known field expansion coefficients in the form y?1î Ì T x?1î . In our case the T -matrix is equal to T Ì ?A
?1î

Ð P1 P

Ð1 2

C

?1î Ð1

î

?B

?1î

Ð P1 P

Ð1 2

D?1î î.

(21)

Note that in an earlier paper [62] we used other auxiliary matrices Q 1 ; Q 2 defined by relation ! ! Q1 x?1î Ì x?L?1î . Q2 y?1î This gave a more simple expression for the T -matrix T Ì Q 2Q
Ð1 1

(22)

.

(23)

However, our use of Q 1 ; Q 2 and Eq. (22) instead of P1 ; P2 and Eq. (18) surprisingly led to several order lower accuracy that could be reached. ? It should be noted that Peterson and Strom [25] suggested a simple recursive relation for T-matrices of layered scatterers. It may be written as T ?Lî Ì ?Q
11

ÐQ

13

T ?L Ð 1îî?Q

31

ÐQ

33

T ?L Ð 1îîÐ1 ,

(24)

where T ?Lî is the T -matrix of an L-layered particle and Q ij are some matrices described, e.g., in [26]. Obviously, Eq. (24) is applicable to the T -matrices arisen in our axisymmetric and nonaxisymmetric tasks with T ?1î Ì G
?L?1î

?E

?L?1î Ð1

î

.

(25)

However, our paper [50] demonstrates that an iterative scheme for the T -matrices being generally equivalent to the recursive relation (24) should be more computationally efficient 0 1 ! ! ?iî ?iî L Y Q 31 ÐQ 33 Q 1 ?1î Q 1 ?Lî @ A Ì . (26) Q 2 ?Lî Q 2 ?1î Q ?iî ÐQ ?iî
iÌ2 11 13 ?i When the EBCM method is applied, elements of the matrices Q jkî and Q j ?1î are some surface integrals and use of the ~ iterative scheme (26) allows one to avoid L inversions of the matrix Q 1 Ì ?Q 31 Ð Q 33 T î made in the recursion relation (24). ~ In the SVM method an extra inversion of Q 1 to get T -matrix at each step in the recursion scheme is less essential as a matrix including A; B; C ; D is anyway inverted for each i to get the matrix with Q ?iî . jk So, we see that the iterative and recursive schemes are tightly related and present an alternative to the single matrix scheme described above. Use of the matrices P 1 ; P 2 or Q 1 ; Q 2 looks to be preferable to that of the T -matrices and then there is practically no difference between iterative and recursive schemes.

3. Numerical results and discussion We have implemented the suggested SVM approach to layered nonspherical scatterers as a Fortran 77 code. Both alternative schemes of solution of the equation system were included. All computations were performed with an Intel 1.8 Hz processor. 3.1. Convergence, accuracy and computational time These three related aspects of our code in the case of layered spheroids are reflected in Fig. 2, where we consider how a measure of relative accuracy of results--relative difference of the extinction and scattering cross-sections for nonabsorbing particles

dÌ

jC C

ext ext

ÐC ?C

sca sca

j

(27)


ARTICLE IN PRESS
1362 A . Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368

Fig. 2. Dependence of an error measure d (left panels) and computational time t (right panels) on the number of terms N kept in the field expansions for 2 and 6-layered prolate spheroids. The aspect ratio of the layer surfaces is the same ?ai =bi Ì 1:5î, the size parameters xV;1 Ì 3, xV;i Ì xV;iÐ1 Ð 0:25, refractive index outside the particle m1 Ì 1 and in even and odd layers mi Ì 1:33 and 1.7, respectively, the incident wave propagation angle a Ì 45 . Results were obtained with our SVM code for the single matrix (SVMsm) and iterative (SVMit) schemes and with the EBCM code from [50].

and computational time t depend on the number of terms kept in the field (potential) expansions N for prolate spheroids with cyclically repeating layers of several materials. The behaviour of accuracy reflected in the dependence d?N î is typical of the SVM, EBCM and PMM methods--when N grows, accuracy first increases but after some value of N rapidly decreases [59]. Similar figures were obtained for another accuracy measure--relative difference in a cross-section obtained when N and N Ð 1 terms are considered

dN Ì

jC N Ð C N C N Ð1

Ð1

j

.

(28)

The dependence of computational time on N may be approximated as t $Na with a Ì 2:5 Ð 3, which is also typical of the methods. The data given for the EBCM are discussed in the next subsection. Comparing two alternative schemes, one can conclude from Fig. 2 that use of a big system allows one to reach better accuracy than application of an iterative/recursive relation. The convergence speed (here the slope of d?N î) is nearly the same, but the former approach is 2-10 0 times slower than the latter one and this difference grows with N and in particular with the number of layers L. Note that standard Gauss-Jordan matrix elimination subroutines used to crash when the number of equations exceeds about 500-10 0 0 (as seen in Fig. 2) and further one should apply sparse matrix inversion procedures. The upper left panel also demonstrates results from the paper [62] where we used Eqs. (22) instead of Eqs. (18). As a result, maximum accuracy reached was lower by a factor of 10 0, though with an increase of L this difference decreased. One can also see that the convergence speed does not depend on the number of layers L, while maximum accuracy reached does. It is still very high (about 10Ð10 ) for a spheroid with six layers. So, even for a scatterer with 10 0 layers one should get several correct digits in results. Required computational time nearly linearly grows with the number of layers when the iteration scheme is applied. 3.2. Comparison with other methods We have compared our SVM code for layered scatterers with available codes based on other methods, namely the EBCM, GMT and DDA. Fig. 3 shows the phase function computed for 2 to 10 0-layered spheroids by our code, a GMT code from [34] and (for Lp10) the standard DDA code [35]. The GMT method described in [34] is often called the null field method with discrete sources (and sometimes a TMM). The book [34] well explains the method where the fields are expanded in terms of spherical functions at different points


ARTICLE IN PRESS
A. Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368 1363

Fig. 3. Phase function computed by the SVM, GMT and DDA methods for a multilayered prolate spheroid. Equivolume layers consist of vacuum ?m Ì 1î and ice ?m Ì 1:33î with their volume fractions being equal to 0.33 and 0.67, respectively. The aspect ratio of the layer surfaces is the same ?ai =bi Ì 1:5î, the particle size parameter xV Ì 3, the axial incidence ?a Ì 0 î.

Table 1 Normalized cross-sections obtained with the SVM and GMT codes.

L

SVM N Q
sca

GMT Q
ext

N 10 8 10 8 8 8

Q

sca

Q

ext

2 6 8 10 50 10 0

38 34 30 30 20 20

0.832175034371 0.86357148727 0.87266563639 0.87911779316 0.9065868 0.911002

0.832175034371 0.86357148726 0.87266563640 0.87911779317 0.9065863 0.9110 04

0.8322 0.8636 0.8722 0.878 0.88 0.9

0.8323 0.8630 0.8723 0.876 0.89 0.8

(sources) fr n gNÌ1 , which can be presented in our notation like the following (cf. Eq. (10)): n E?iî ?r î Ì 1
Mrank Nrank XX mÌ0 nÌ1 ?iî ?anm M r m;m?l ? ?r Ð r niî î? bnm N ?iî r m;m?l ? ?r Ð r niî îî,

(29)

where l Ì 1 for m Ì 0 and l Ì 0 otherwise. The GMT code was configured to choose sources automatically, the number of integration points was N int Ì 1000, the number of expansions Nrank was varying from 4 to 40 but kept the same for all layers, the number of azimuthal expansions M rank was determined by the code automatically. The figure demonstrates that the results obtained with the GMT and our SVM codes match very well for particles with up to 10 layers. Essential difference at large scattering angles can be observed only for a very large number of layers. A deeper insight in accuracy of the results is provided by Table 1 where we give normalized extinction and scattering cross-sections for the spheroids


ARTICLE IN PRESS
1364 A . Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368

considered in Fig. 3. Computational time required by the GMT and our SVM codes to obtain the data in the table was comparable. This may be related only with a large number of integration points used to be taken in the GMT code as generally this method should be faster than the SVM. Table 1 indicates a good match of the SVM results with the law of conservation of energy (Q ext is equal to Q sca for dielectric particles). This fact shows that the proposed SVM method is very robust to an increase of the layer number and should provide reliable results even for particles with more than a hundred of layers. The DDA code applied is a very popular CDM tool employed to solve the light scattering problem for nonspherical inhomogeneous particles. We added a simple subroutine modelling dipole locations for a multilayered particle (target) and used the basic code described in [35]. The total number of dipoles was N dip Ì 1:5 Ò 106 , but similar results were obtained also for Ndip Ì 105. There is a problem of the DDA method when applied to scatterers with a large (over about 10) number of layers--it becomes hardly possible to keep the correct volume fractions of materials involved as the layers become too narrow for appropriate distribution of dipoles even when N dip exceeds millions. Hence in Fig. 3 we present only data obtained for 4 and 10-layered particles. One can well see how small difference between the DDA and other methods for L Ì 4 becomes essential for L Ì 10. Nevertheless, the DDA being also a rather slow method is still very useful in treatment of complex scatterers with a few layers. The EBCM code used was described in [50]. It includes an iterative scheme based on Eq. (22) that could be improved as discussed above. Further improvement of performance of the code can be expected after use of the approach utilized in the well-known EBCM/TMM code of Mishchenko [1]. However, Fig. 2 has clearly demonstrated that the EBCM method is not quite suitable for treatment of layered scatterers. The reason may be in the fact observed for homogeneous particles for which the applicability range of the EBCM (in contrast to the PMM and SVM) is strictly limited in space of particle shape parameters (see discussion in [59]). When more shapes (of layers) are involved, this range quickly decreases. An illustration is Fig. 4 for core-mantle spheroids. It is well seen from the figure that the EBCM provides acceptable accuracy only for the pffiffiffi aspect ratio a1 =b1 o 2 (while the SVM works well beyond this limit). Note that for homogeneous spheroids the EBCM gives accurate cross-sections for a=b44210. When more layers are considered, convergence of the EBCM becomes worse (see Fig. 2) and one can hardly obtain reliable results for the number of layers larger than 5. 3.3. Illustrative examples To illustrate capacity of the method and code suggested, we present some more numerical results obtained for porous particles. Figs. 5 and 6 show normalized scattering cross-section and phase function, respectively, for spheroids with

Fig. 4. Accuracy of cross-sections of core-mantle spheroids achieved by the SVM and EBCM methods in dependence on the particle aspect ratio a1 =b1 and the size parameter xV;1. Other parameters are as in Fig. 2.


ARTICLE IN PRESS
A. Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368 1365

Fig. 5. Normalized scattering cross-section Q sca in dependence on the size parameter xV for 2, 6, 10, 20, and 100-layered prolate spheroids. Materials of equivolume layers are vacuum ?m Ì 1î and dirty ice ?m Ì 1:33 ? i0:01î. The volume fraction of vacuum (porosity) is 0.33 (left panel) and 0.9 (right panel). The aspect ratio of all layers is the same ai =bi Ì 1:5, the axial incidence ?a Ì 0 î. For homogeneous spheroids, refractive index was obtained using the Maxwell-Garnett rule.

Fig. 6. Phase function for 2, 6, 10, 20 and 10 0-layered spheroids. The outermost layer size parameter xV;1 Ì 3, other parameters are as in Fig. 5.

cyclically repeating layers of ice and vacuum. The number of layers L changes from 2 to 10 0 for two values of particle porosity (the volume fraction of vacuum) P Ì 0:33 and 0.9. Note that the (integral) cross-section for P Ì 0:9 and the phase function (differential cross-section) for P Ì 0:33 rather weakly depend on L when L44. In contrast, this cross-section for P Ì 0:33 and the phase function for P Ì 0:9 essentially change with an increase of L even when L$100. We believe that there should be some limit values of the cross-section and the phase function as the number of layers tends to infinity. It took place for layered spheres considered in [63], where its author needed, however, about several thousands of layers to approach this limit. The behaviour of Figs. 5-6 confirms this conclusion for spheroids. The results presented in this section make more understandable the following Table 2 where we very roughly estimate properties of the methods mentioned in Section 1 when they are applied to layered scatterers. The table shows the type of scatterers to which the given method is applicable, a maximum reachable size parameter xV Ì 2pr V =l, where rV is the radius of a sphere which volume is equal to that of a nonspherical particle, a maximum aspect ratio a=b when the method is


ARTICLE IN PRESS
1366 A . Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368

Table 2 Methods applied to layered nonspherical scatterers. Methoda EMT RA GMT SVMe EBCM SVMsphh CDM FDTDk
a b c d e f g h i j k

Particle type Any Confocal ellipsoidsc Any Axisymmetric Axisymmetricf Confocal spheroids Any Any
b

xV Medium to large Small Medium Medium Mediumg Large Medium Medium(?)
b

a=b Medium to large Medium to larged Medium to large Medium Mediumg Large Large Large(?)
b

L
Any Any Large Large Small Large Mediumj Medium(?)

Accuracy Low Medium to low Medium to high High High High Medium Medium(?)

Speed High High Medium Medium Medium Medium Low Low

i

See the names of the methods in Section 1. Depends on the method used for homogeneous particles with an effective refractive index. Properties of an extension to nonconfocal ellipsoids in [46] are not clear. Large for so called quasistatic approximation [46]. The method developed in this work. For nonaxisymmetric scatterers the speed becomes rather low. Medium to large when extended precision is used. SVM with a spheroidal basis. Properties of an extension to nonconfocal spheroids in [24] are not clear. When dipoles form a cubic lattice. Properties of the FDTD method when applied to layered scatterers are not quite clear.

applied to spheroidal scatterers, a possible number of layers L, accuracy of results and speed of the method relative to others. The table is not a verdict for the methods, but just a background to better realize the place of the suggested approach (labelled by SVM in the table) among others. We see that our method can be useful when one needs rather efficiently to treat simple shape nonspherical scatterers with a (very) large number of layers if these particles are more complex than spheroids with the confocal layer boundaries. Note that such spheroids form a very special set of particles as any layer shape is strictly determined by its volume and the particle aspect ratio. Our method can also be applied when one requires (very) high accuracy of results for layered scatterers of a simple shape. It may be in particular useful as a testing tool for other codes. 4. Conclusions We have developed a special version of the SVM method for axisymmetric multilayered scatterers. Numerical tests show that our approach provides high accuracy results for simple shape particles with up to 10 0 layers. A consideration of two general computational schemes used for layered scatterers has been performed. We find that the single system scheme allows one usually to reach more accurate results than the iterative or recursive scheme. However, the former requires an order of magnitude longer computational time than the latter and is hardly applicable to scatterers with 10 and more layers without use of sparse matrix inversion algorithms. Comparison of the developed SVM code with available codes based on other (EBCM, GMT, DDA) methods applicable to layered particles allows one to see the field of applicability of our code. It should be most useful for treatment of layered particles being more complex than spheroids with strictly confocal layer boundaries and having a rather large number of layers or any number of layers if high accuracy of results is required. Note that solution of these tasks is hardly possible by other light scattering methods.

Acknowledgements The authors are thankful to B. Draine and Th. Wriedt for providing the DDA and GMT codes, respectively, and M. Prokopjeva for performing calculations with the DDA code. The work was partly supported by the Grants RFFI 07-02-00831, NTP 2.1.1/665 and NSh 1318.20 08.2. References
[1] [2] [3] [4] Mishchenko MI, Hovenier JW, Travis LD. Light scattering by nonspherical particles. San Diego: Academic Press; 20 00. Kahnert FM. Numerical methods in electromagnetic scattering theory. JQSRT 20 03;79-80:775-824. ? ? ? Mie G. Beitrage zur Optik Truber Medien, speziell kolloidaler Metallosungen. Ann Phys 1908;25:377-445. Fuller KA, Mackowski DW. Electromagnetic scattering by compounded spherical particles. In: Mishchenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press; 20 0 0. p. 225-72. [5] Aden AL, Kerker M. Scattering of electromagnetic waves from two concentric spheres. J Appl Phys 1951;22:1242-6.


ARTICLE IN PRESS
A. Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368 1367

[6] Gonzales F, Moreno F. Introduction to light scattering from microstructures. In: Moreno F, Gonzales F, editors. Light scattering from microstructures. Lectures notes in physics, vol. 534. Berlin: Springer; 20 0 0. p. 1-19. [7] Babenko VA, Astafyeva LG, Kuzmin VN. Electromagnetic scattering by disperse media. London: Springer-Praxis; 2003. [8] Gurwich I, Shiloah N, Kleiman M. Calculations of the Mie scattering coefficients for multilayered particles with large size parameters. JQSRT 20 01;70:433-40. [9] Yang W. Improved recursive algorithm for light scattering by a multilayered sphere. Appl Opt 20 03;42:1710-20. [10] Li R, Han X, Jiang H, Ren KF. Debye series for light scattering by a multilayered sphere. Appl Opt 20 06;45:1260-70. [11] Li R, Han X, Shi L, Ren KF, Jiang H. Debye series for Gaussian beam scattering by a multilayered sphere. Appl Opt 20 07;46:4804-12. [12] Kerker M, Matijevich E. Scattering of electromagnetic waves from concentric infinite cylinders. J Opt Soc Am 1961;51:506-8. [13] Shah GA. Scattering of plane electromagnetic waves by infinite concentric circular cylinders at oblique incidence. Mon Not R Astron Soc 1970;148:93-102. [14] Yousif HA, Elsherbeni AZ. Electromagnetic scattering from an eccentric multilayered cylinder at oblique incidence. J Electrmagn Wave Appl 1999;13:325-36. [15] Gurwich I, Shiloah N, Kleiman M. The recursive algorithm for electromagnetic scattering by titled infinite circular multi-layered cylinder. JQSRT 1999;63:217-29. [16] Li R, Han X, Jiang H, Ren KF. Debye series of normally incident plane-wave scattering by an infinite multilayered cylinder. Appl Opt 2006;45:6255-62. [17] Onaka T. Light scattering by spheroidal grains. Ann Tokyo Astron Obs 1980;18:1-54. [18] Ciric IR, Cooray FR. Separation of variables for electromagnetic scattering by spheroidal particles. In: Mishchenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press; 20 00. p. 89-130. [19] Gurwich I, Kleiman M, Shiloah N, Cohen A. Scattering of electromagnetic radiation by multilayered spheroidal particles: recursive procedure. Appl Opt 20 00;39:470-7. [20] Gurwich I, Kleiman M, Shiloah N, Oaknin D. Scattering by an arbitrary multi-layered spheroid: theory and numerical results. JQSRT 20 03;79-80:649-60. [21] Barton JP. Internal, near-surface and scattered electromagnetic fields for a layered spheroid with arbitrary illumination. Appl Opt 20 01;40:3598-607. [22] Farafonov VG. New recursive solution of the problem of scattering of electromagnetic radiation by multilayered spheroidal particles. Opt Spectra 20 01;90:743-52. [23] Li LW, Kang XK, Leong MS. Spheroidal wave functions in electromagnetic theory. New York: Wiley; 20 02. [24] Han Y, Zhang H, Sun X. Scattering of shaped beam by an arbitrarily oriented spheroid having layers with non-confocal boundaries. Appl Phys B 20 06;84:485-92. ? [25] Peterson B, Strom S. T -matrix formulation of electromagnetic scattering multilayered scatterers. Phys Rev D 1974;10:2670-84. [26] Wang DS, Barber PW. Scattering by inhomogeneous nonspherical objects. Appl Opt 1979;18:1190-7. [27] Mishchenko MI, Videen G, Babenko VA, Khlebtsov NG, Wriedt T. T -matrix theory of electromagnetic scattering by particles and its applications: a comprehensive reference database. JQSRT 20 04;88:357-406. [28] Mishchenko MI, Videen G, Babenko VA, Khlebtsov NG, Wriedt T. Comprehensive T -matrix reference database: a 20 04-20 06 update. JQSRT 20 07;106:304-24. [29] Mishchenko MI, Videen G, Babenko VA, Khlebtsov NG, Wriedt T, Zakharova NT. Comprehensive T -matrix reference database: a 20 06-20 07 update. JQSRT 20 08;109:1447-60. [30] Petrov D, Shkuratov Y, Zubko E, Videen G. Sh-matrices method as applied to scattering by particles with layered structure. JQSRT 20 07;106:437-54. [31] Kyurkchan AG, Demin DB, Orlova NI. Solution based on the pattern equation method for scattering of electromagnetic waves by objects coated with dielectric materials. JQSRT 20 07;106:192-202. [32] Al-Rizzo HM, Tranquilla JM. Electromagnetic scattering from dielectrically coated axisymmetric objects using the generalized point-matching technique II. Numerical results and comparison. J Comput Phys 1995;119:356-73. [33] Doicu A, Wriedt T. Null-field method with discrete sources to electromagnetic scattering from layered scatterers. Comput Phys Commun 20 01;138:136-42. [34] Doicu A, Wriedt T, Eremin Y. Light scattering by systems of particles. Berlin: Springer; 20 06. [35] Draine BT, Flatau PJ. User guide for the discrete dipole approximation code DDSCAT 6.1. arXiv preprint: astro-ph/0409262; 20 04. [36] Gordon HR, Du T. Light scattering by nonspherical particles: application to coccoliths detached from Emiliania huxleyi. Limnol Oceanogr 20 01;46:1438-54. [37] Somasiri NP, Chen X, Robertson ID, Rezazadeh AA. Comparison of numerical techniques for full wave analysis of multi-layered microstrip antenna. Int J Infrared Millim Waves 20 02;23:1777-85. [38] Kim J, Rahmat-Samii Y. Simulations, design and characterization of implanted antennas inside a human body. IEEE Trans Microwave Theory Tech 20 04;52:1934-43. [39] Karlowski A. Improving accuracy of FDTD simulations in layered biological tissues. IEEE Microwave Wireless Components Lett 20 04;14:151-2. [41] Winton SC, Kosmas P, Rappaport CM. FDTD simulations of TE and TM plane waves at nonzero incidence in arbitrary layered media. IEEE Trans Antennas Propag 20 06;53:1721-8. [42] Choi MK. Numerical calculation of light scattering from a layered sphere by boundary-element method. J Opt Soc Am 20 01;18:577-83. [43] Rysakov W, Ston M. The light scattering by core-mantle sphere. JQSRT 20 01;69:121-9. ? [44] Seydou F, Duraiswami R, Gumerov NA, Seppanen T. TM electromagnetic scattering from multilayered dielectric bodies--numerical solution. ACES J 2004;19:10 0-7. [45] Medgyesi-Mitschang LN, Eftimiu C. Scattering from axisymmetric obstacles embedded in axisymmetric dielectrics: the method of moments. Appl Phys 1979;19:275-85. [46] Posselt B, Farafonov VG, Il'in VB, Prokopjeva MS. Light scattering by multilayered ellipsoidal particles in the quasistatic approximation. Meas Sci Technol 20 02;13:256-62. [47] Sloot PMA, Fidgor CG. Elastic light scattering from nucleated blood cells: rapid numerical analysis. Appl Opt 1986;25:3559-65. [48] Kokhanovsky AA. Optics of light scattering media: problems and solutions. Berlin: Wiley-Praxis; 1999. [49] Bohren CF, Huffman DR. Absorption and scattering of light by small particles. New York: Wiley; 1983. [50] Farafonov VG, Il'in VB, Prokopjeva MS. Light scattering by multilayered nonspherical particles: a set of methods. JQSRT 20 03;79-80:599-626. [51] Voshchinnikov VB, Il'in VB, Henning T. Modelling the optical properties of composite and porous interstellar grains. Astron Astrophys 20 05;429:371-81. [52] Athanasiadis C. On the acoustic scattering amplitude for a multi-layered scatterer. J Aust Math Soc Ser B 1998;39:431-48. [53] Charalambopoulos A, Dassios G, Fotiadis DI, Massalas CV. Scattering of a point generated field by a multilayered spheroid. Acta Mech 20 01;150:4107-19. [54] Gutman S. Identification of multilayered particles from scattering data by clustering method. J Comput Phys 20 00;163:529-46. [55] Yan G, Zhao H. The far field operator for a multilayered scatterer. Comput Math Appl 2002;43:631-9. [56] Barton JP. Electromagnetic field calculations for an irregularly shaped near-spheroidal particle with arbitrary illumination. J Opt Soc Am A 20 02;40:3598-607. [57] Farafonov VG, Vinokurov AA, Il'in VB. Comparison of the light-scattering methods using the spherical basis. Opt Spektrosk 20 07;102:1006-16. [58] Farafonov VG. Diffraction of a plane electromagnetic wave by a dielectric spheroid. Diff Equat (Sov) 1983;19:1765-77. [59] Farafonov VG, Il'in VB. In: Kokhanovsky AA, editor. Light scattering reviews. Berlin: Springer-Praxis; 20 06. p. 125-77.


ARTICLE IN PRESS
1368 A . Vinokurov et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1356 -1368

[60] Farafonov VG, Il'in VB. Single light scattering by simple model nonspherical heterogeneous particles. Berlin: Springer-Praxis; 20 09, in preparation. [61] Voshchinnikov NV, Farafonov VG. Optical properties of spheroidal particles. Astrophys Space Sci 1993;204:19-86. [62] Farafonov VG, Vinokurov AA. Light scattering by multilayered axisymmetric particles: solution of the problem by the separation of variables method. Opt Spektrosk 20 08;105:318-31. [63] Johnson BR. Exact theory of electromagnetic scattering by a heterogeneous multilayered sphere in the infinite-layer limit: effective-media approach. J Opt Soc Am 1999;16:845-52.