MULTIPLE SCATTERING OF RADIATION IN CLUSTERS OF DIELECTRICS

A fast, accurate, and general technique for solving Maxwell{close_quote}s equations in the presence of a finite cluster of arbitrarily disposed dielectric objects is presented. The electromagnetic field is first decomposed into multipoles with respect to centers close to each of the objects of the cluster and multiple scattering is carried out until convergence is achieved. Radiation scattering cross sections are obtained using this method for clusters formed by homogeneous spheres and coated spheres made of different materials (Al, Si, and SiO{sub 2}), including magnetic ones. Near- and far-field distributions are offered as well. {copyright} {ital 1999} {ital The American Physical Society}


I. INTRODUCTION
The electromagnetic response of small structures to external sources constitutes a field of intense research due to the promising technological applications of photonic devices. 1 These are in essence materials tailored on the micrometer scale to exhibit unusual properties of transmission 2 and reflection 3 of radiation that can be exploited in the design of filters and mirrors for light. Special emphasis has been placed in the study of periodic structures in connection to photonic bands, which are formed in a similar way as electronic bands in solids. [4][5][6][7][8][9][10][11]1,12,13 The photonic properties of these materials have been tested in macroscopic structures using radiation in the GHz regime, where absolute band gaps have already been obtained. 9,10 On the theory side, different methods have been developed to solve Maxwell's equations in the presence of periodic arrays, including extensions of those used in low-energy electron diffraction, 4,5,13 plane-wave 7 and Bloch wave 8 expansions, and the transfermatrix approach. 11,14 Actual photonic microstructures have been recently employed to confine light in a finite region of space, 15,16 resulting in well-defined narrow modes that could be eventually used in laser design. This opens the field of photonic chemistry, where arbitrarily distributed micrometer elements are combined to confine, scatter, or emit light. Therefore, theoretical methods suited to solve the electromagnetic problem near clusters of dielectric objects are needed to explore the possibilities offered by these structures as they become increasingly complex. The present work constitutes an attempt to advance in this direction.
The solution of Maxwell's equations finds application in a number of spectroscopy techniques, including spontaneous emission, 17 light emission from scanning tunneling microscopes, 18 scanning near-field optical microscopy ͑SNOM͒, 19,20 and electron-energy-loss spectroscopy ͑EELS͒. 21,14 For instance, SNOM permits one to obtain spatial resolution on a nanometer scale by using an external probe to bring the light of a laser of much larger wavelength to the area of the specimen under examination, 19 posing the problem of understanding the measured far-field induced by local interaction of laser light with structures whose size is much smaller than the wavelength. This technique combines spatial localization and radiation of energy appropriate to sample valence-band features. 20 The simulation of EELS relies on the availability of methods to calculate the field induced by an external electron interacting with complex nanostructures. Beyond some nonrelativistic [22][23][24] and relativistic [25][26][27] analytical solutions for simple geometries, numerical methods have been successfully used in this respect. 14 In particular, the relativistic boundary element method 28 consists of representing the interfaces of a given heterostructure by means of charges and currents that are solved self-consistently in the presence of an external field. 29 However, the magnitude of the numerical problem becomes unaffordable for targets composed of a large number of elements. In an attempt to circumvent this issue, granular structures have been studied using effectivemedium theories, taking advantage of averaging over random distributions of objects. [30][31][32][33][34] The formalism introduced in this paper allows one to reduce the solution of Maxwell's equations in the presence of a cluster of arbitrarily distributed dielectric objects to the knowledge of the individual scattering properties of each of the objects. For that purpose, the electromagnetic field is first decomposed into multipoles around each constituent of the cluster, and multiple elastic scattering of the multipole expansions ͑MESME͒ is carried out until convergence is achieved. The time required to numerically solve this problem is proportional to the square of the number of objects in the cluster. This permits one to compute radiation scattering cross sections for a cluster formed by a large number of objects within any desired degree of accuracy. 35 The present formalism shares many features with other theories employed to study electron diffraction in solids. [36][37][38][39][40][41][42] This is a consequence of the fact that electrons moving in the interstitial region between the atoms of a solid, where the potential is nearly flat, are governed by the same Helmholtz wave equation that rules the propagation of free photons, the only difference lying in the matching conditions satisfied by either the electron wave functions or the electromagnetic fields. 1 The theory of MESME is presented in Sec. II, where a general, computationally efficient technique is derived that allows us to solve Maxwell's equations in clusters of arbitrarily disposed dielectrics. Section III describes further computational details. The application to the simulation of radiation scattering is given in Sec. IV. Some numerical examples

II. MULTIPLE ELASTIC SCATTERING APPROACH TO THE ELECTROMAGNETIC PROBLEM
The electromagnetic field induced by interaction of an external source with a cluster of dielectric objects will be solved in terms of MESME around the objects of the cluster. Some of the underlying ideas involved in this formalism have been borrowed from cluster models for the simulation of electron diffraction in solids, [39][40][41][42] and adapted to deal with the electromagnetic field rather than the electron wave function, exploiting the well-known analogy between electrons in solids and light in nanostructures. 1,12 A similar extension from electrons to photons was carried out by Ohtaka and co-workers 4,5,13 allowing methods developed for the simulation of low-energy electron diffraction to be employed in the study of photonic band structures for periodic arrays of spheres, and, more recently, for two-dimensional lattices of cylinders. 43,44 Instead, finite clusters of arbitrarily distributed objects are considered here, and the photonic properties of a variety of systems are investigated. A full curved-wave ͑mul-tipole͒ treatment of the electromagnetic field is used, in combination with improved iterative multiple-scattering techniques. The electron wave functions employed in the electron-diffraction analogy are replaced here by scalar functions that represent the electromagnetic field. 45 The analysis carried out in this section combines the following elements. 35 ͑a͒ The external electromagnetic field ͑e.g., the field set up by incoming radiation or by an external fast electron͒ is decomposed into multipoles around each object of the cluster. This decomposition is reviewed in Sec. II A using a notation appropriate for its application to MESME.
͑b͒ The scattering on each individual dielectric object is represented by complex scattering matrices, which permit one to obtain the contribution of each object to the scattered field, as discussed in Sec. II B.
͑c͒ The field resulting from scattering on a given object needs to be propagated to other objects of the cluster where further scattering can take place. This involves nontrivial translations of multipoles, which will be examined in Sec. II C.
͑d͒ Finally, the solution of Maxwell's equations is expressed in terms of self-consistently-calculated scattered fields. This leads to the MESME secular equation presented in Sec. II D.
Therefore, this formalism permits one to express the solution of Maxwell's equations in the presence of a cluster of arbitrarily disposed dielectric objects in terms of the individual scattering properties of the constituents of the cluster. Spherical objects, for which the scattering properties are collected in analytical expressions given in Sec. II B, will be considered here for simplicity. The details of the formalism are given next and a more schematic picture is offered in Fig.  1 and its caption.

A. Electromagnetic multipoles and scalar functions
Let us start by expressing the electromagnetic field in terms of multipoles with respect to a given origin r ␣ within a homogeneous region of space free of charges and currents. The electric field E is by necessity transversal in that region, and thus, it can be expressed in frequency space as 45,27 where kϭ/c, L ␣ ϭϪi(rϪr ␣ )ϫ" is the orbital angular momentum operator relative to the position r ␣ , and ␣ M and FIG. 1. Diagrammatic representation of the elements involved in the solution of Maxwell's equations in the presence of a cluster of dielectric objects using multiple elastic scattering of multipole expansions ͓MESME; see Eq. ͑18͒, reproduced in this figure as well͔. The electromagnetic field is expressed in terms of scalar functions ␣ , made up of multipoles relative to dielectric objects ␣, ␤, etc. The external field acting on object ␣ ͑upper right corner of this figure͒, ␣ ext , is a superposition of spherical plane waves with no net energy flux around the object ͓double-arrow line; see Eq. ͑6͔͒. Its scattering at ␣, represented by the scattering matrix t ␣ , gives rise to a superposition of scattered outgoing waves ͓t ␣ ␣ ext ; see Eqs. ͑7͒ and ͑8͔͒ which adds to the induced part of the total field produced by self-consistent scattering on object ␣, ␣ ind . The remaining contribution to ␣ ind is provided by the scattering of the self-consistent induced field coming from each other object ␤ ␣. The latter is calculated in various steps ͑starting from the upper left corner͒: first, the system is rotated by means of a rotation matrix R ␣␤ such that the bond vector r ␣ Ϫr ␤ is made to point along the positive z axis; the resulting rotated outgoing waves centered at r ␤ (R ␣␤ ␤ ind ) are then expressed in terms of spherical plane waves centered at r ␣ by linearly translating spherical harmonics (G ␣␤ z is the operator for translations along the z axis, as explained in Sec. II C͒; an additional translation operator T ␣␤ z is necessary to compensate for the lack of invariance of multipoles under translation of the origin; the result is then rotated back to the original position and scattered at object ␣. Finally, summation over ␤ ␣ yields Eq. ͑18͒. ␣ E are the so-called magnetic and electric scalar functions, respectively. If the region under consideration is filled with medium j and described by its frequency-dependent local dielectric function ⑀ j and magnetic permeability j , one finds, upon insertion of Eq. ͑1͒ into Maxwell's equations, that the scalar functions satisfy the wave equation where k j ϭkͱ⑀ j j ͑the square root is chosen here to have a non-negative imaginary part͒. Besides, the magnetic field is found to be and the scalar functions can be determined from E using the identities 45,27 Notice that the longitudinal modes are explicitly left out of this formalism, 45,27 preventing from possible numerical problems. Equation ͑2͒ implies that the multipole expansion of the electromagnetic field in the homogeneous region under consideration can be constructed as a sum of free spherical waves. The absence of sources in that region indicates that these waves cannot lead to a net energy flux through any closed surface contained within it, so that they can be termed spherical plane waves by analogy to conventional plane waves. In particular, the scalar functions that describe the external field can be expanded in terms of spherical harmonics Y L as where Lϭ(l,m), j L (u)ϭi l j l (͉u͉)Y L (û ) represents one of the noted spherical plane waves, j l is a spherical Bessel function, and groups both magnetic and electric components.

B. Single scattering on a dielectric object
Now let us focus on a dielectric object located near r ␣ and surrounded by a homogeneous medium jϭ0 of dielectric function ⑀ 0 and magnetic permeability 0 , so that, for a given frequency component , the momentum of the electromagnetic field is k 0 ϭkͱ⑀ 0 0 ͑if jϭ0 refers to the vacuum, one has ⑀ 0 ϭ 0 ϭ1 and k 0 ϭkϭ/c).
The total electric field is the superposition of the external field and the field induced by scattering on the object under consideration, that is, EϭE ext ϩE ind . In addition, the electromagnetic field in the homogeneous medium jϭ0 is a combination of outgoing and incoming spherical waves, represented by spherical Hankel functions h l (ϩ) (k 0 r) and h l (Ϫ) (k 0 r), respectively. 46 In particular, E ind finds its sources in the charges and currents induced by the external field in the dielectric object, and therefore, it has to be a combination of only outgoing waves. In other words, the corresponding scalar functions can be written where h L (ϩ) (u)ϭi l h l (ϩ) (͉u͉)Y L (û ). The superscript ss refers to the fact that the field induced in the presence of just one object can be considered to be the result of single scattering by comparison to the case of a cluster of several objects, where multiple scattering becomes relevant, as shown latter in this work. Equation ͑7͒ is valid for r outside a sphere centered at r ␣ and fully containing the dielectric object ͑i.e., containing the sources of the induced field͒.
Within the linear-response approximation, the components of the scattered field have to be proportional to those of the external field. Therefore, one can write, in terms of the coefficients of expressions ͑6͒ and ͑7͒, where t ␣,LL Ј is the so-called scattering matrix. This is schematically shown in the upper right corner of Fig. 1, where the external field acting on ␣, represented by a double-arrow line to emphasize the fact that it is made up of spherical plane waves j L , gives rise to outgoing scattered waves h L (ϩ) that generate ␣ ss , represented by outgoing arrows. The elements of the scattering matrix can be determined for each L by solving Maxwell's equations in the presence of the dielectric object with the asymptotic condition implicitly defining t ␣,L Ј L , and the requirement that be finite everywhere. If r ␣ lies within a homogeneous region of space j, possibly inside the dielectric object, the finiteness of means that only spherical plane waves j L Љ ͓k j (rϪr ␣ )͔ contribute at that point, since both outgoing and incoming spherical waves diverge at their origin.
Although the formalism presented in this work can be applied to arbitrarily shaped objects, whose scattering matrices are generally dense and can be obtained numerically, 29 we will only consider for simplicity spherically symmetric objects, for which the matching conditions satisfied by the fields ͑i.e., the continuity of the normal displacement, the tangential electric field, the normal magnetic induction, and the tangential magnetic field͒ reduce, after using Eqs. ͑1͒ and ͑3͒, to the continuity of ␣ M , ⑀ ␣ E , (1/)(1ϩr‫ץ/ץ‬r) ␣ M , and (1ϩr‫ץ/ץ‬r) ␣ E , where local response of the materials involved in the system is assumed and ⑀ and are the frequency-dependent response functions that take values ⑀ j and j inside medium j. Thus magnetic and electric scalar 3 ϭk 2 a, and 0 ϭka. This expression converges smoothly to the homogeneous sphere limit when both b˜0 or aϭb. It also reproduces the polarizability coefficients in the limit c˜ϱ. 48 Radiation scattering cross sections are presented for isolated spheres in Figs. 2 and 3, and they are discussed in Sec. V.

C. Translation of electromagnetic multipoles
When more that one dielectric object is considered, the results of scattering from each object ͑e.g., ␤) need to be propagated to each of the rest of the objects of the cluster ͑e.g., ␣), where further scattering events can take place. This involves both ͑i͒ translation of the spherical harmonics appearing in the multipole expansions discussed above, and ͑ii͒ translation of the origin of multipoles from r ␤ to r ␣ along the bond vector d ␣␤ ϭr ␣ Ϫr ␤ .
͑i͒ Translation of spherical harmonics. The spherical harmonics and spherical Hankel functions that made up the multipole expansion of the results of scattering at r ␤ are given in coordinates relative to r ␤ . They can be expressed with respect to a new origin r ␣ by using the formula for translation of spherical harmonics, [36][37][38] FIG. 2. Total radiation-scattering cross section for an isolated sphere of SiO 2 as a function of incoming photon energy for various values of the sphere radius a, as indicated by means of labels. The cross section, calculated by using Eq. ͑30͒, has been normalized to the projected area a 2 .
FIG. 3. Dependence of the elastic cross section on scattering angle for isolated SiO 2 spheres of radius ͑a͒ aϭ20 nm and ͑b͒ aϭ160 nm, as given by the integral of Eq. ͑29͒ over azimuthal directions. The differential cross section has been multiplied by sin and it is given per degree and normalized to the projected area a 2 . The contour curves limiting white areas correspond to values of 0.001, 0.011, and 0.06 in the left part of figure ͑a͒, the right part of figure ͑a͒, and figure ͑b͒, respectively. The distance between consecutive contour curves corresponds to fixed values of 0.0002, 0.001, and 0.005, respectively.
is the Green function of Eq. ͑2͒, written in the basis set of spherical harmonics attached to r ␣ and r ␤ , and is a Gaunt integral. When d ␣␤ is directed along the positive z axis ͑i.e., the quantization direction͒, the Green function reduces to A recurrence relation has been reported that permits us to evaluate G ␣␤,LL Ј z efficiently. 41 Notice that G ␣␤ propagates magnetic and electric components separately.
͑ii͒ Translation of the origin of multipoles. Unfortunately, the scalar functions are not invariant under translations of the origin of coordinates. This comes from a lack of invariance of the angular momentum operator involved in Eq. ͑1͒, which transforms as when changing the origin from r ␤ to r ␣ . Therefore, the electric field induced by scattering at r ␤ can be expressed as in terms of operators relative to r ␣ . Following the discussion of Sec. II A, E ␤ ind can be written where ␣␤ ind accounts for the contribution to the induced field coming from scattering at ␤ and expressed in terms of spherical plane waves around ␣. Now magnetic and electric components of ␣␤ ind can be obtained from Eqs. ͑4͒ and ͑5͒ after substituting E by E ␤ ind as given by Eq. ͑12͒, leading to the following rule of transformation of the scalar functions under translation of the origin along the vector d ␣␤ : 49 where the substitution ٌ 2˜Ϫ k 0 2 has been performed. Now a linear translation operator T can be defined to represent this rule. T acts on the multipole components of the above scalar functions and takes a particularly simple form when d ␣␤ is directed along the positive z axis. More explicitly, as shown in the Appendix. Notice that ␤ ind is made up of spherical outgoing waves (h L (ϩ) ) centered at r ␤ ͑this is represented by outgoing arrows centered around ␤ in the upper left corner of Fig. 1͒. The translation of spherical harmonics discussed in point ͑i͒ above allows one to express ␤ ind in terms of spherical plane waves ( j L ) centered at a new origin r ␣ ͑see the diagram on the left part of Fig. 1͒. On the other hand, the translation of the origin of multiples discussed in point ͑ii͒ operates on ␤ ind ͓already translated as explained in point ͑i͔͒ to produce ␣␤ ind , which is also made up of spherical plane waves centered at r ␣ .
In brief, the combination of these two operators allows us to obtain the coefficients of the multipole expansion of ␤ ind around ␣ as T ␣␤ G ␣␤ ␤ ind , where ␤ ind represents the vector formed by the coefficients ␤,L ind of the multiple expansion of ␤ ind around ␤.

D. Multiple scattering
Now let us focus on a cluster of dielectric objects labeled by coordinate vectors r ␣ . Our aim is to obtain the scalar functions appropriate for the electromagnetic field that satisfies Maxwell's equations in the presence of the cluster, and to relate these functions to the scattering properties of each individual component of the cluster. The induced part of the self-consistent scattered field will be expressed as the sum of contributions coming from the different cluster components ␣, that is, Since ␣ ind finds its origin in the charges and currents induced in object ␣, it has to be a combination of spherical outgoing waves centered around r ␣ : is valid for r outside a sphere centered at r ␣ and fully containing object ␣.
Within the single-scattering ͑ss͒ approach, ␣ ind is given by the scattered field ␣ ss discussed in Sec. II B. This results from single scattering of the external field as expressed in Eq. ͑8͒, conveniently written using matrix notation as where the vector ␣ ss(ext) has components ␣,L ss(ext) and the matrix t ␣ has components t ␣,LL Ј .
The singly scattered field coming from a certain object ␤ can in turn suffer scattering on every other object of the cluster ␣ ␤. That is, ␣ ind receives contributions coming from previous scattering on every other object ␤ ␣. This leads to a self-consistent relation for the induced field that can be written as where the first term on the right-hand side is the result of direct single scattering of the external field at ␣; the second term describes both the propagation of the self-consistently scattered field originating in every other object ␤ ␣ from ␤ to ␣ and the subsequent scattering of this field at ␣; ϭ1 is introduced for convenience; and the operator H ␣␤ accounts for the propagation just noted. H ␣␤ can be conveniently constructed in four steps as follows ͑see the left part of Fig. 1 for a schematic representation of this procedure͒. 35 ͑i͒ Following previous work of electron diffraction in solids, 36,39 the bond vector d ␣␤ is rotated onto the z axis by using a rotation matrix R ␣␤ , 46 which acts on the spherical harmonics of the multipole expansion of ␤ ind . The Euler angles corresponding to this rotation can be chosen (0,, Ϫ) if (,) are the polar angles of d ␣␤ .
͑ii͒ The resulting rotated scalar functions are then propagated a distance d ␣␤ along the positive direction of the z axis. This is accomplished by multiplying by the Green function G ␣␤ z given by Eq. ͑11͒.
͑iii͒ The lack of invariance of multipoles under translations of the origin needs to be overcome by multiplying by the linear translation operator T ␣␤ z described in Sec. III C ͓Eq. ͑15͔͒.
͑iv͒ Finally, the z axis has to be rotated back onto the d ␣␤ direction using R ␣␤ Ϫ1 , and one has The analogy with electron diffraction in solids 38,39 is complete, except for the lack of invariance just pointed out. The rotation matrices can be in turn decomposed into azimuthal and polar rotations as 46 and this helps to soften the computational demand in the evaluation of Eq. ͑18͒, as shown in Sec. III.

III. COMPUTATIONAL PROCEDURE
An efficient scheme is presented in this section that allows us to solve Eq. ͑18͒ within affordable limits in both computation time and storage demand. The operators involved in Eqs. ͑18͒-͑20͒ are approximated by finite matrices of dimension ͓(l max ϩ1) 2 ͔ 2 , where l max is the maximum orbital angular momentum number under consideration ͑this is for each component, electric and magnetic; note that Green functions and rotation matrices act independently on each of these components͒. In the calculations that follow, convergence has been achieved for l max р12 in most cases.
The direct inversion of Eq. ͑18͒ is computationally prohibitive for large clusters. A procedure for solving it that mimics the multiple-scattering expansion consists in starting with ␣ ind,1 ϭ ␣ ss as a guess for ␣ ind , and then calculating the result of scattering up to an order nϾ2 by using the iterative relation This is nothing but the Taylor expansion of ␣ 's in powers of t ␣ 's. This procedure has been found computationally convenient in many cases where ␣ ind,n˜ ␣ ind as n˜ϱ, particularly in electron diffraction in solids, 39,42 but it may lead to divergences when any of the eigenvalues of tH has a modulus larger than 1.
Here Eq. ͑18͒ has been solved by using the recursion method, 50,51 where plays the same role as the energy in former electronic band-structure calculations. Although one is only interested in the value ϭ1, the recursion method results advantageous because it is fully convergent even in situations where Eq. ͑21͒ leads to divergences. 38 This is the case in many of the examples offered below, where strong scatterers are placed relatively close to each other.
The computational costs of both the iteration procedure and the recursion method is basically coming from the multiplications H ␣␤ ␤ ind . These two methods require the same number of those multiplications per iteration, namely, N(N Ϫ1), where N is the number of objects in the cluster. The factorization of H ␣␤ given in Eqs. ͑19͒ and ͑20͒ has the virtue of reducing both ͑a͒ the storage capacity required to evaluate those multiplications and ͑b͒ the computational effort.
͑a͒ A significant reduction in memory requirements can be accomplished if the coefficients of each polar rotation R mm Ј (l) (0,,0) ͓see Eq. ͑20͔͒, each azimuthal rotation (Ϫ1) mЈ e imЈ , each translation of spherical harmonics G ␣␤ z , and each translation of the origin of multipoles T ␣␤ z are computed and stored only the first time that they are encountered along the entire calculation. In addition, for any arbitrary cluster, the points r ␣ can be chosen in such a way that they form a highly symmetrical mesh, where many bond distances and angles are repeated ͑this might require that the points r ␣ do not lie necessarily in the geometrical centers of the objects that they label͒. If this is the case, the total number of different bond distances and bond polar angles is considerably smaller than the number of bond vectors d ␣␤ . To illustrate this, let us put forward the example of a simple-cubiclattice cube of side p in units of the lattice constant; this cluster contains p 3 nodes and (2pϪ1) 3 Ϫ1 different bond vectors, a number that has to be compared with at most 3p 2 bond distances, since the square of the distance between any pair of nodes has to equal an integral number, and the distance between opposite corners is ͱ3p ͑a better estimate for this case results in Ϸ1.8p 2 different bond distances for large p values͒.
͑b͒ For a given maximum value of the orbital angularmomentum number l max , the dimension of each vector ␣ is (l max ϩ1) 2 ͑for each component, electric and magnetic͒, so that every matrix-vector product H ␣␤ ␤ involves 4(l max ϩ1) 4 complex multiplications. However, all of the matrices that appear on the right-hand side of Eq. ͑19͒ are sparse, as one can see from Eqs. ͑11͒, ͑15͒, and ͑20͒. A detailed inspection leads to the conclusion that only Ϸ(20/3)(l max ϩ1) 3 complex multiplications are needed to evaluate the product H ␣␤ ␤ when H ␣␤ is decomposed as shown in Eq. ͑19͒. This is a factor of Ϸl max /2 smaller than the direct matrix-vector product.
Further reduction in computational and storage demand can be achieved if symmetry relations for the Green functions and the rotation matrices 46 are used ͑e.g., G ␣␤,lm,l Ј m z ϭG ␣␤,l Ј Ϫm,lϪm z ). Notice that multiplications by scattering matrices t ␣ do not affect significantly the total computational costs for relatively large clusters ͑e.g., above ten objects͒, since they take place just outside the summation over cluster objects ␤ in Eq. ͑18͒, so that only N of those multiplications are needed per interation.
A fully automated implementation of these ideas has been performed, resulting in a new code ͑MESME͒ that provides the solution of Eq. ͑18͒ for clusters of arbitrarily distributed constituents by investing a computation time where A is a constant (Aϳ10 Ϫ4 s on a Pentium at 333 MHz͒.

IV. APPLICATION TO THE SCATTERING OF RADIATION
The formalism presented in Sec. II will be applied here to the study of scattering of radiation in a cluster of arbitrarily distributed dielectric objects. Section IV A will be devoted to derive analytical expressions for the multipole expansion coefficients of an incoming plane wave. The far field produced by the charges and currents induced in the cluster will be discussed in Sec. IV B.

A. Multipoles of the incoming radiation
When the cluster is illuminated by a plane wave, the external electric field can be written 52 where ⑀ ជ is the ͑complex͒ polarization vector, which is assumed to be normalized as ͉⑀ ជ ͉ϭ1, and K i is the momentum of the incoming light. E ext must satisfy the wave equation in the surrounding medium jϭ0, and therefore, ͉K i ͉ϭk 0 ϭkͱ⑀ 0 0 , where kϭ/c. 52 The multipole coefficients of the scalar functions introduced in Sec. II A can be obtained by expanding Eq. ͑22͒ in spherical plane waves and inserting the resulting expression into Eqs. ͑4͒ and ͑5͒. 45 One finds where ⍀ i refers to the polar angles of K i ,

B. Induced electromagnetic field in the far-field limit and scattering cross section
The induced electric field can be obtained by inserting into Eq. ͑1͒ the self-consistently calculated induced part of the scalar functions, made up of spherical outgoing waves h L (ϩ) ͓k 0 (rϪr ␣ )͔ in the interstitial medium jϭ0 just outside the cluster objects. At very large distances from the cluster, the spherical outgoing waves behave like where K f ϭk 0 r. Furthermore, " can be substituted by iK f in Eq. ͑1͒ ͑this represents the leading term in the r˜ϱ limit͒. Using Eqs. ͑1͒, ͑3͒, and ͑17͒, the induced electric and magnetic far fields are found to be and is the scattering amplitude and ⍀ denotes the polar angles of r.
The radiated energy associated to the induced far field can be derived from the energy flux across an arbitrarily large sphere of radius r centered at the cluster, and this flux can be in turn calculated from the integral of the normal Poynting vector as 52 where the integral over the time has been included. Expressing the fields in terms of their frequency components, one finds can be interpreted as the probability of radiating a photon of energy per unit of energy range and unit of solid angle around the direction ⍀. Inserting Eqs. ͑25͒ and ͑26͒ into Eq. ͑28͒, and noticing that f-rϭ0, 53 one finds represents the differential cross section for elastically scattered photons, which can be computed from the coefficients of the induced scalar functions using Eqs. ͑24͒ and ͑27͒. The total interaction cross section can be divided as tot ϭ el ϩ inel , where inel accounts for inelastic losses, connected to absorption of the incoming radiation by objects of the cluster. The total cross section can be expressed in terms of the imaginary part of the scattering factor f along the forward direction by means of the optical theorem for radiation, which reads 52

V. RESULTS AND DISCUSSION
Clusters formed by spheres made of Al, SiO 2 , and Si are considered next. Their scattering matrices have been calculated using Eq. ͑9͒, where the frequency-dependent dielectric functions of these materials are employed. In particular, the response of Al has been approximated by a Drude dielectric function with bulk plasma energy p ϭ15 eV and damping ϭ1.06 eV. The dielectric functions of Si and SiO 2 have been taken from optical data. 54 Figure 2 shows the total cross section of isolated SiO 2 spheres, as calculated from Eq. ͑30͒. Different sphere radius a have been considered and the results have been normalized to the projected area a 2 . Substantial variations in the cross section can be observed in the region below 10 eV, where the dielectric function has a smooth structure and is basically real. As the sphere size increases, the low-energy peak in the cross section shifts toward lower energies. This is connected to retardation effects in the electromagnetic signal when the sphere radius is comparable in size to the wavelength of the radiation ͑see the upper scale in Fig. 2͒. In a simplified picture where one considers a frequency-independent dielectric constant ͑actually, this is nearly the case in SiO 2 in the energy range 3-8 eV, where Re͕⑀͖Ϸ2.1Ϫ2.9ӷIm͕⑀͖), one FIG. 4. Total radiation-scattering cross section for two touching SiO 2 spheres of radius aϭ160 nm as a function of incoming photon energy for different orientations of the cluster ͑see the inset͒. The incoming radiation is right circularly polarized ͑RCP͒. The cross section has been normalized to the projected area of two spheres 2a 2 . The result for an isolated sphere ͑multiplied by a factor of 2͒ is also shown for comparison. The wavelength of the radiation is given on the upper scale, normalized to a.
can argue that the frequency of the eigenmodes will scale with a in such a way that a/c is a dimensionless constant number, and hence, the larger the radius a, the smaller the frequency . For energies above 12 eV the normalized cross section depends very weakly on a.
Angular distributions of elastically scattered photons are represented in Fig. 3 for two different radii of SiO 2 spheres, as calculated from Eq. ͑29͒. The angular distribution of scattered radiation has a marked anisotropic character. Notice in particular that forward scattering appears to be dominant over a wide energy range for the larger sphere under consideration (aϭ160 nm͒.
Actually, this anisotropy in scattering is translated into an orientational dependence of the total cross section in the case of the two-sphere SiO 2 cluster of Fig. 4, where right circularly polarized radiation 52 ͑RCP͒ has been considered. The result derived from isolated spheres ͑thick solid curve, multiplied by a factor of 2͒ is very close to the one obtained when the dimer is oriented perpendicular to the direction of the incoming radiation. If the dimer is aligned with the di-rection of incidence of the radiation, forward scattering on the lower sphere focuses scattered radiation on the second sphere, and this contributes to enhance the role of multiple scattering in this case ͑thin solid curve͒. 55 It should be pointed out that the single-scattering approach results in a total cross section that equals the sum of cross sections of the components of the cluster ͓see Eqs. ͑27͒ and ͑30͔͒.
In Fig. 5, the near field has been represented for the two extreme orientations of the SiO 2 dimer considered above, and also for the isolated sphere. The magnitude represented in the figure is the square of the induced electric field, which is discontinuous on the sphere surfaces, represented by thick circumferences. Notice that the variations of the field in each sphere relative to the case of an isolated sphere are very small. Nevertheless, the maximum of the induced field, located on the side of the sphere opposite to the incoming radiation, enhances the interaction between the two spheres of the dimer when this is oriented along the direction of incidence of the radiation. Elastic and total cross sections, given by Eqs. ͑29͒ and ͑30͒, respectively, have been compared in Fig. 6͑a͒ for clusters of Al spheres of radius aϭ19 nm. Elastic-scattering cross sections have been obtained by integrating the differential cross section over angles of scattering. The elastic cross section lies always below the total one, and the difference between the two of them accounts for absorption of photons by the spheres. The different features of the energy dependence of the cross section for the isolated sphere are shifted and split when one considers clusters of spheres. A significant orientational dependence can be observed in the case of the dimer, where low-energy modes of different energy are excited depending on the orientation of the cluster relative to the incoming radiation. In particular, peak A stays at nearly the same energy for clusters of three and four spheres, indicating that this feature is connected to spheresphere interaction in dimers oriented perpendicularly with respect to the direction of the incoming beam. RCP radiation is considered in Fig. 6͑a͒, but a strong dependence on the polarization is observed for oriented clusters, like the dimer of Fig. 6͑b͒. Only the electric-field component of the external radiation parallel to the dimer axis contributes to excite the low-energy peak at 5.2 eV. Notice the triple-crossing point E, which occurs due to both the symmetry of the cluster, implying that RCP ϭ LCP ͑i.e., right and left circularly polarized radiation lead to the same cross section͒, and the general identity RCP ϩ LCP ϭ LPʈx ϩ pЌx , where the notation of Fig. 6͑b͒ has been adopted.
The scattering on Al spheres has a strong electric character for the sphere radius under consideration. Actually, if one sets t l M ϭ0 in the calculations presented in Fig. 6, the results are nearly indistinguishable on the scale of the figure. This was expected from previous results on the magnitude of the magnetic components of the scattering matrix for Drude spheres. 27 However, for the SiO 2 spheres of Figs. 2 and 4 the magnetic components play a significant role even for the smallest spheres (aϭ20 nm͒, so that t l M cannot be dismissed in that case.
The scattering-angle distributions of some of the elastic cross sections discussed in Fig. 6͑a͒ are analyzed in Fig. 7. The dimer oriented along the direction of incidence of the radiation ͓Fig. 7͑b͔͒ has a focusing effect on the scattered radiation, enhancing forward scattering as compared to the case of an isolated sphere 55 ͓see Fig. 7͑a͔͒. The two orientations of the dimer contemplated in Fig. 7 showed certain similarities in the integrated cross section for Ͼ7 eV ͓see Fig. 6͑a͔͒, but their angular distributions reveal notorious differences. Lower-energy peaks are observed when the dimer is oriented perpendicular to the radiation direction ͓Fig. 7͑c͔͒. Notice the symmetry between forward and backward scattering in the low-energy region in this case. Instead, the high-energy region has a dominant forward-scattering character. When a third sphere is added to form a triangle oriented normal to the incoming radiation ͓Fig. 7͑d͔͒, there are not substantial changes as compared to the dimer case. However, a fourth sphere on top of the triangle ͓see Fig.  7͑e͔͒ and a fifth underneath ͓Fig. 7͑f͔͒ produce an effect of focusing of the scattered radiation along the forward direction.
Even more dramatic differences between different clusters are observed in the doubly-differential angular distributions FIG. 6. ͑a͒ Total radiation scattering cross section ͑solid curves͒ and elastic scattering cross section ͑broken curves͒ for clusters of 1-5 Al spheres of radius aϭ19 nm. The separation between sphere surfaces is 2 nm. The cross section has been normalized to the number of spheres N in each case and the projected area of one sphere a 2 . The radiation is RCP, and is moving upwards, as shown in the upper inset. The sphere clusters are disposed as illustrated in the rest of the insets. Consecutive curves have been shifted a value of 3 upwards to improve readability. The wavelength of the radiation is given on the upper scale, normalized to a. ͑b͒ Dependence of the total scattering cross section on the polarization of the incoming radiation for the horizontal two-sphere cluster considered in ͑a͒. The cases of linear polarization ͑LP͒ parallel and perpendicular to the bond vector have been considered, as well as circular polarization.
FIG. 7. Dependence of the elastic cross section on scattering angle for clusters of 1-5 Al spheres of radius aϭ19 nm, disposed as shown in the insets. The separation between sphere surfaces is 2 nm. The incoming radiation is RCP, and it is incident as shown by white arrows. The differential cross section has been integrated over azimuthal directions and multiplied by sin and it is given per degree and normalized to both the number of spheres in each case and the projected area of one sphere a 2 . Contour curves limiting white areas correspond to values of 0.035 in ͑a͒ and ͑b͒, 0.022 in ͑c͒-͑e͒, and 0.026 in ͑f͒. Consecutive contour curves differ by a value of 0.005 in ͑a͒ and ͑b͒, and 0.002 in ͑c͒-͑f͒.
represented in Fig. 8, where radiation of energy corresponding to features A and B of Fig. 6͑a͒ has been contemplated. The figure illustrates how different polarizations of the incoming radiation lead also to very different diffraction patterns. The symmetry of the cluster is reflected in that of the far field, as one can see by comparing the two-sphere case with the rest of the clusters. However, a rotation from the directions of symmetry of the cluster is observed is the case of circularly polarized radiation. The intensity of scattered photons is maximum along the forward direction (ϭ0) in all cases.
Near-field distributions ͑i.e., the square of the induced electric field͒ have been represented in Fig. 9 for isolated Al spheres and clusters formed by two and three Al spheres. The parameters of the spheres are the same as in Fig. 6. The energy of the radiation corresponds to the features marked in Fig. 6͑a͒. Circularly polarized radiation has been considered. Notice that the two modes illustrated in the figure for the isolated sphere show very different distributions of induced electric field, which can be ascribed to dipole and quadrupole modes in Figs. 9͑a͒ and 9͑b͒, respectively. For the dimer, a large enhancement of the induced field can be observed in the region between the two spheres in Figs. 9͑c͒ and 9͑d͒, indicating a strong sphere-sphere interaction. The lack of specular symmetry in the plane normal to the radiation in Fig. 9͑d͒ can be attributed to the polarization of the radiation. For the same energy but using radiation aligned with the dimer, the 5.2-eV mode does not show up in the cross section, as can be seen in Fig. 6͑a͒, and this is translated into the very weak induced field shown in Fig. 9͑e͒. Finally, the case of three spheres ͓Fig. 9͑f͔͒ shows a strong localization of the FIG. 8. Distribution of the elastic cross section ͓Eq. ͑29͔͒ over directions of scattering for various clusters of 2-4 Al spheres of radius aϭ19 nm whose surfaces are separated by 2 nm. Each row represents results obtained from the same cluster ͑see insets͒, with the radiation coming normal to the plane of representation from underneath ͑i.e., ʈz). Various polarizations of the incoming light have been considered ͑the polarization employed in each column is indicated in the upper insets͒. The energy of the radiation corresponds to that of points A and B in Fig. 6͑a͒ ͑i.e., A ϭ5.2 eV and B ϭ8.2 eV͒. ͑a͒-͑p͒ represent forward-scattering results ͑i.e., angles of scattering ϭ0Ϫ90°, with ϭ0 in the center of each figure͒, whereas ͑q͒-͑t͒ correspond to backward scattering ͑i.e., angles of scattering ϭ90°-180°, with ϭ180°in the center of each figure͒. The doubly differential cross section is given per stereoradian, and normalized to both the number of spheres in each case and the projected area of one sphere a 2 . Contour curves limiting white areas correspond to values of 0.3 in ͑a͒ and ͑d͒, 0.6 in ͑b͒, 0.15 in ͑c͒, 0.4 in ͑e͒-͑p͒, and 0.09 in ͑q͒-͑t͒. Consecutive contour curves differ by a value of 0.05 in ͑a͒, ͑d͒, and ͑e͒-͑p͒, 0.1 in ͑b͒, 0.025 in ͑c͒, and 0.01 in ͑q͒-͑t͒.
FIG. 9. Contour maps of the square of the induced-electric-field strength in the vicinity of Al spheres of radius aϭ19 nm illuminated by RCP radiation of energy as indicated by labels ͓see corresponding points in Fig. 6͑a͔͒. The cases of an isolated sphere ͓͑a͒ and ͑b͔͒, a two-sphere cluster with different orientations relative to the incoming radiation ͓͑c͒-͑e͔͒, and a three-sphere cluster ͑f͒ have been considered. The planes of representation are schematically shown in the insets relative to both the incoming beam direction and the position of the spheres. The sphere surfaces are represented by thick circumferences on the figures. The separation between the surfaces of the spheres is 2 nm in ͑c͒-͑f͒. Contour curves limiting white areas correspond to a value of 5, 5.5, 8, 8, 3.5, and 9 a.u. in ͑a͒-͑f͒, respectively. The difference in value between consecutive contour curves is 0.5 a.u. in all cases. The strength of the external field has been normalized to 1 a.u. induced field in the interstitial region between the spheres, as well as a certain asymmetry also connected to the polarization of the radiation.
The near field induced by linearly polarized light incident along different angles with respect to the normal to a triangle formed by three Al spheres of radius aϭ19 nm, like those considered in Fig. 9͑f͒, is represented in Fig. 10. When the polarization vector is contained in the plane determined by the normal to the triangle and the direction of incidence of the radiation ͓s polarization; Figs. 10͑a͒-10͑c͔͒, the near field is very sensitive to the angle of incidence. For normal incidence, the polarization direction has a large projection on two of the bond vectors of the cluster, producing an intense induced electric field in those regions. For in-plane incidence, the polarization vector is normal to all bond vectors, and, as a result, the induced field is very weak. This is in qualitative agreement with the conclusions extracted from Figs. 6͑b͒ and 9͑c͒-9͑e͒, where it was shown that only the electric field component parallel to the bond vector contributes to excite a mode at the energy A ϭ5.2 eV under consideration. For p polarization ͓Figs. 10͑d͒-10͑f͔͒ the near field is rather insensitive to the incidence angle, and this can be attributed to the fact that the projections of the polarization vector on the bond vectors of the cluster are independent of the angle on incidence in this case. Also, the wavelength of the radiation is over 17 times larger than the sphere radius ͑see upper scales in Fig. 6͒, and, hence, phase differences in the interaction of external radiation with the spheres for various incidence angles are small.
The evolution of the induced electric field on its way out from the target is represented in Fig. 11. The near field is shown on different planes parallel to the same three-Alsphere cluster considered above ͓Figs. 11͑a͒-11͑f͔͒, and also on concentric hemispheres centered at the cluster ͓Figs. 11͑g͒-11͑j͔͒. The electric field evolves smoothly, forming different patterns that cannot be directly identified with the far field up to relatively large distances from the cluster. For hemispheres at distances of 150 and 200 nm, like those contemplated in Figs. 11͑i͒ and 11͑j͒, the induced field starts resembling the far field represented in Fig. 8͑h͒ for the same geometry. Upon inspection of Eq. ͑17͒, one concludes that the far-field approximation is valid at distances R from the cluster much larger than the cluster size and such that k 0 R Ӷͱl max . In the present case, this condition becomes R Ӷ125 nm.
The maxima of the induced electric field are located just outside the spheres in the case of Al ͑see Figs. 9-11͒, whereas for SiO 2 the induced field is larger inside the spheres ͑see Fig. 5͒. This is related to the fact that the wave equation ͑2͒, which rules the propagation of free scalar functions in Al ͑SiO 2 ), is the same as the Schrödinger equation for free electrons moving in an repulsive ͑attractive͒ potential, since the real part of the dielectric constant is negative ͑positive͒ for the radiation energies under consideration.
Single-scattering ͑ss͒ results have been compared with full multiple elastic scattering of multipole expansions ͑MESME͒ in Fig. 12 for four-Al-sphere clusters of different sizes. Multiple scattering effects are shown to be important, and they lead to dramatic modifications in the cross section. When the dimensions of the cluster are increased by scaling its geometrical parameters, a general shift of the maxima in the cross sections toward lower energies is observed as a result of retardation effects similar to those discussed above for SiO 2 spheres.
A larger target has been considered in Fig. 13, where the effect of multiple scattering of external radiation has been studied for a cluster of 60 spheres with the same structure as carbons in a C 60 molecule. For SiO 2 spheres ͓Fig. 13͑a͔͒, the total scattering cross section ͑solid curves͒, as calculated from Eq. ͑30͒, follows quite closely the elastic-scattering cross section ͑dashed curves͒ in the region below Ϸ8.5 eV, where SiO 2 is transparent and absorption is negligible. For Al spheres ͓Fig. 13͑b͔͒, both cross sections are relatively featureless. The cross section of the Al cluster lies below the result of single scattering ͑dotted curves, representing 60 times the total cross section of an isolated sphere͒ in most of the energy range under consideration. In silica, the prominent feature obtained for single scattering at around 9 eV is converted into a dip when multiple scattering is switched on.
The case of Si spheres coated with SiO 2 has been illustrated in Fig. 14. Isolated spheres of radius aϭ20 nm have been considered in Fig. 14͑a͒ for different radii of the inner Si core, b. A prominent feature dominates the spectrum for the case of pure Si ͑i.e., bϭa). The exciton of SiO 2 shows up at around 10.5 eV when the sphere is coated, and it remains as a dominant feature for the case of pure SiO 2 . For three touching spheres, multiple scattering is observed to play a significant role when the coating layer is thin. Notice that a coating layer of just 2 nm of SiO 2 is able to produce sizable variations in the cross section. A shift of the modes FIG. 11. Contour maps of the square of the induced-electric-field strength in the vicinity of a cluster of three Al spheres of radius aϭ19 nm illuminated by RCP radiation. The energy of the radiation is A ϭ5.2 eV. The surfaces of representation are schematically shown in the insets relative to both the incoming beam direction and the position of the spheres. Different planes perpendicular to the direction of the radiation have been considered in ͑a͒-͑f͒. The distance from the planes to the centers of the spheres is 5 ͓ϩ5͔ 30 nm in ͑a͒-͑f͒, respectively. The near field has been also represented for points lying on concentric hemispheres whose distance from the cluster center is 50 ͓ϩ50͔ 200 nm in ͑g͒-͑j͒, respectively. In this case, the representation is based upon polar angles with respect to the forward direction determined by the center of the hemispheres, in a similar way as in Fig. 8͑h͒. The sphere surfaces are represented by thick circumferences in ͑a͒-͑c͒. The separation between the surfaces of the spheres is 2 nm. Contour curves limiting white areas correspond to a value of 4 a.u. in ͑a͒-͑d͒, 2.2 in ͑e͒, and 1.2 in ͑f͒. The difference in value between consecutive contour curves is 0.5 a.u. in ͑a͒-͑d͒, and 0.2 in ͑e͒ and ͑f͒.
The maximum values in ͑g͒-͑j͒ are 2.5, 0.36, 0.24, and 0.19 a.u., respectively. The strength of the external field has been normalized to 1 a.u. toward lower frequencies is observed here again when multiple scattering is switched on.
Metallic spheres of frequency-independent negative dielectric function are contemplated in Fig. 15. For nonmagnetic isolated spheres ͓Fig. 15͑a͔͒, a general trend toward lower energies is observed in the features of the cross-section spectrum when the magnitude of the dielectric constant increases. Also, the cross section increases in magnitude in the same direction, and this is consistent with the fact that negative dielectric functions act like repulsive potentials for light, so that, the larger the magnitude of ⑀, the stronger the interaction. For a cluster of three nonmagnetic spheres ͓Fig. 15͑b͔͒, multiple-scattering effects are important and they are enhanced when the magnitude of the dielectric function increases ͑notice in particular the dramatic change that occurs when ⑀ goes from Ϫ20 to Ϫ40͒.
Finally, the lower part of Fig. 15 is devoted to studying magnetic spheres as a function of the magnetic permeability for ⑀ϭϪ10. The attenuation of the electric field inside the spheres increases with following the exponential law exp(ϪrͱϪ⑀), so that and ͉⑀͉ play a similar role in this respect. Actually, the modes of both isolated spheres ͓Fig. 15͑c͔͒ and three-sphere clusters ͓Fig. 15͑d͔͒ are also shifted toward lower energies when the permeability increases.

VI. CONCLUDING REMARKS
In summary, a general, computationally efficient technique has been presented for solving Maxwell's equations in a cluster of arbitrarily disposed dielectric objects. This method permits one to express the solution of Maxwell's equations in terms of the scattering properties of the constituents of the cluster. The computation time scales with the square of the number of elements of the cluster.
The generality of this method relies on the fact that there is no restriction on the shape or internal structure of the constituents of the clusters under consideration, apart from the condition imposed by the translation formula of spherical harmonics ͓see Eq. ͑10͔͒, which requires that the objects of the cluster can be embedded inside nonoverlapping spheres. Work to overcome this restriction is in progress.
The present formalism has been applied to the study of scattering of radiation and different examples of simulation for clusters formed by spheres have been offered, including clusters of magnetic spheres and coated spheres. The effect of multiple scattering of radiation is dramatic in some of the clusters under consideration. In general, metallic clusters are able to produce larger multiple scattering effects than clusters formed by insulators.
Near-field distributions have been presented for various clusters. Notice that there is a connection between the elec-FIG. 12. Total scattering cross section for RCP radiation incident on tetrahedral clusters of Al spheres of various sizes. Single scattering results ͑ss; dashed curves͒ and full multiple-scattering results ͑MESME; solid curves͒ are shown for clusters of spheres of radius aϭ0.475d, whose centers are separated a distance d ͑see labels͒. The radiation is coming along a direction perpendicular to one of the faces of the tetrahedron, as shown in the inset. The cross section has been normalized to both the number of spheres and the projected area of one sphere a 2 . Consecutive curves have been shifted a value of 5 upwards to improve readability.
FIG. 13. ͑a͒ Radiation scattering cross section for a cluster formed by 60 SiO 2 spheres with the same structure as carbons in a C 60 molecule. The radius of the spheres is 49.5 nm. The nearestneighbor bond distance is 100 nm. The elastic cross section ͓dashed curve; integral of Eq. ͑29͔͒ is compared with the total cross section ͓solid curve; Eq. ͑30͔͒, and 60 times the total cross section of an isolated sphere ͑dotted curve͒. ͑b͒ The same as ͑a͒ for Al spheres.
tric field induced at a position near a cluster illuminated by an external plane wave incident along a certain direction, and the radiation emitted along that direction due to a localized source at the position under consideration. 13 This reciprocity theorem permits us to extract conclusions on the information contained in SNOM images by analyzing the near field induced in the cluster by external radiation. In this sense, the method presented here can be applied to the analysis of SNOM in complex geometries.
The present theory can deal with different external sources. For instance, if the external field is provided by a fast electron, the induced field acting back on the electron can be used to simulate electron energy losses in scanning transmission electron microscopes.
Finally, an interesting possibility is offered by the extension of quasicrystals 56,57 to their photonic counterpart, where large clusters of dielectrics would be distributed to form a nonperiodic structure able to reflect light forming fivefold symmetry patterns. Well-defined electronic bands are observed in quasicrystals. 56 Similarly, photonic band gaps have been predicted in two-dimensional photonic quasicrystals, 58 and one would expect similar results in three-dimensional photonic structures, where cluster models like the one developed here are ideal to study local properties.

APPENDIX: TRANSLATION OF THE ORIGIN OF ELECTROMAGNETIC MULTIPOLES
This appendix is devoted to deriving Eq. ͑15͒, which represents the explicit form of the transformation rule of multipole coefficients of scalar functions under translations of the origin of multipoles from r ␤ to r ␣ along the vector d ␣␤ ϭr ␣ Ϫr ␤ , when the latter is directed along the positive z FIG. 14. Total scattering cross section for RCP radiation incident on both ͑a͒ an isolated Si sphere coated with SiO 2 and ͑b͒ a cluster of three touching spheres. The sphere radius is aϭ20 nm in all cases and the radius of the Si core b is changed as shown by labels. In ͑b͒, the radiation is coming normal to the plane defined by the sphere centers. The cross section has been normalized to the projected area of one and three spheres in ͑a͒ and ͑b͒, respectively. Consecutive curves have been shifted 0.2 upwards to improve readability.
FIG. 15. Total scattering cross section for RCP radiation incident on both an isolated sphere ͓͑a͒ and ͑c͔͒ and a cluster of three identical spheres along the direction normal to the plane of the sphere centers ͓͑b͒ and ͑d͔͒, as a function of frequency . Nonmagnetic spheres (ϭ1) have been considered in ͑a͒ and ͑b͒ with different values of the dielectric constant (⑀ϭϪ4, Ϫ5, Ϫ6, Ϫ8, Ϫ10, Ϫ15, Ϫ20, and Ϫ40). Results for magnetic spheres are shown in ͑c͒ and ͑d͒ for ⑀ϭϪ10 in all cases and different values of the magnetic permeability (ϭ1, 2, 4, 6, 8, 10, 12, 14, and 16). The results have been normalized to the projected area of either one sphere or three spheres, in each case. The separation between the surfaces of the spheres is a/10 in the three-sphere cluster, where a is the sphere radius. The wavelength of the radiation ϭ2c/ is given on the upper scale, normalized to a. Consecutive curves have been shifted 0.2 upwards to improve readability in ͑c͒ and ͑d͒.
axis. This can be accomplished by projecting Eqs. ͑13͒ and ͑14͒ onto multipole components, noticing that each multipole under consideration is actually a spherical plane wave j L ͓see comments immediately after Eq. ͑6͒ and those at the end of Sec. II C͔. In particular, the ␦ function in Eq. ͑15͒ comes from the first term on the right-hand side of Eqs. ͑13͒ and ͑14͒. Moreover, the last term in these equations can be projected using the identity where u ␣ ϭk 0 (rϪr ␣ ), Lϭ(l,m), and the substitution (1/L ␣ 2 ) j L˜j L /l(lϩ1) has been performed. This gives rise to crossed electric-magnetic terms in Eq. ͑15͒.
The remaining second term requires a more careful analysis. One can proceed by considering the operator Using the identity ٌ 2 j L ϭϪk 0 2 j L ͓see Eq. ͑2͒ and notice that the medium filling the interstitial region between objects of the cluster corresponds to jϭ0͔, Eq. ͑A2͒ leads to ẑ •͑L ␣ ϫ" ͒ j L ͑ u ␣ ͒ϭϪi lϩ1 k 0 ͭ u ␣ ͓ j l ͑ u ␣ ͒ϩ j l Љ͑u ␣ ͔͒Y L ͑ u ␣ ͒ ϩ ͫ j l Ј͑u ␣ ͒Ϫ j l ͑ u ␣ ͒ u ␣ ͬ ϫ͑1Ϫ 2 ‫ץ͒‬ Y L ͑ u ␣ ͒ ͮ , where ϭ(zϪz ␣ )/͉rϪr ␣ ͉. Now, expressing Y L and (1 Ϫ 2 ‫ץ)‬ Y L in terms of Y lϩ1,m and Y lϪ1,m ͑see for instance the appendices of Messiah 46 ͒ and using the recurrence relations of the spherical Bessel functions j l , one finds where D L is defined by Eq. ͑16͒. Finally, using Eqs. ͑A1͒ and ͑A3͒ in the multipole expansion of Eqs. ͑13͒ and ͑14͒, and rearranging the summation index l, one obtains Eq. ͑15͒.