Tryfon Antonakakis

Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy

John Desanto

Gérard Granet

Editor Evgeny Popov

Boris Gralak

Sébastien Guenneau

Daniel Maystre

André Nicolet

Brian Stout
email: brian.stout@fresnel.fr

Fréderic Zolla

Benjamin Vial

GratinGs: theory and numeric applications

Introduction to particulate gratings

Lattice sums of spherical harmonic functions are well suited for modeling gratings composed of periodic arrays of identical discrete particles, henceforth referred to as particulate gratings (cf. fig. (6.1)). By discrete particles, we mean that the particles have a physical boundary such that there exists a region between the individual particles that is governed by the host material's constitutive relations. This feature makes particulate gratings somewhat different from most of the other diffraction grating problems studied in this book which are usually characterized by a substrate and a superstrate with distinct constitutive parameters. The techniques of this chapter can be extended to include the effects of a nearby planar interface, [START_REF] Stefanou | Scattering of light from a two-dimensional array of spherical particles on a substrate[END_REF][START_REF] Stefanou | Scattering of electromagnetic waves by a disordered twodimensional array of spheres[END_REF][START_REF] Stefanou | Heterostructures of photonic crystals: frequency bands and transmission coefficients[END_REF] but such considerations complicate the problem somewhat and this chapter therefore concentrates on substrate-free particulate gratings. Theoretical analysis of the particulate grating problem can draw on both single-particle scattering theory and techniques originally developed for solid state physics. The solid state analogy is clear from the similarity of this problem to the scattering of waves by crystal lattices, particularly in the "muffin tin" approximation [START_REF] Slater | Wave functions in a periodic potential[END_REF]. Summations of the spherical harmonic fields scattered by the (infinite) number of particles in the lattice involve semi-convergent series and will generally go under the name of "lattice sums". By lattice sum, we mean sums of the form ∑ Λ Φ r r r j where Λ refers to the ensemble of points, r r r j , in a periodic lattice, and Φ is a given function.

Lattice sums have applications in many fields and their study dates back to the 19th century treating conditionally convergent sums of solutions to the Laplace equations (most notably in the Madelund constant of ionic crystals). Nevertheless, they were not always recognized as a specific branch of study, and their derivations tended to be scattered throughout the literature. This situation is changing however with the appearance of extensive reviews in recent years [START_REF] Linton | Lattice sums for the helmholtz equation[END_REF][START_REF] Linton | One-and two-dimensional lattice sums for the threedimensional helmholtz equation[END_REF][START_REF] Moroz | Exponentially convergent lattice sums[END_REF][START_REF] Moroz | Quasi-periodic green's functions of the helmholtz and laplace equations[END_REF]. Furthermore, another monograph, dedicated entirely to lattice sums, is appearing at the same time as this one. [START_REF] Borwein | Lattice Sums : Then and Now[END_REF] As developed in detail in the aforementioned monograph, the study of the (scale invariant) Laplace equation lattice sums have generated a number of important analytic results. The grating problem on the other hand involves propagating waves and consequently requires lattice sums of Helmholtz-type solutions. Although there are fewer fully analytic results for the (scale dependent) Helmholtz lattice sums than for the Laplace equation case, analytic manipulations remain essential for regularizing and accelerating a numerical analysis.

In solid-state physics, Helmholtz (i.e. Schrödinger) equation lattice sums are a key aspect of the Korringa-Kohn-Rostoker (KKR) methods for band-structure calculations in crystals. [START_REF] Papanikolaou | Conceptual improvements of the kkr method[END_REF][START_REF] Korringa | On the calculation of the energy of a bloch wave in a metal[END_REF][START_REF] Kohn | Solution of the schrödinger equation in periodic lattices with an application to metallic lithium[END_REF] In KKR theory, lattice sums intervene in the calculation of the "structure constants" of the lattice Green function and their regularization generally goes under the name of Ewald sum techniques. The Ewald sum method is quite intricate, but its basic principle can be viewed as separating a semi-convergent sum into slowly and rapidly convergent parts and then to transform the slowly convergent part into reciprocal space via the Poisson sum formula where it becomes a rapidly convergent series.

Although Ewald sum methods are proven to be quite efficient for most of the problems encountered in solid state physics, there utility has been repeatedly criticized for gratingtype applications (requiring numerically unwieldy evaluations of incomplete Gamma functions with negative real arguments [START_REF] Winitzki | Computing the incomplete gamma function to arbitrary precision[END_REF], poor numerical properties for high multipole orders or large wavenumber, k ). A number of authors have consequently looked for alternative lattice sum techniques since the pioneering work of Kerker over 30 years ago. In this chapter, we simply discuss and compare some of our preferred methods in the appendices. Our emphasis will instead be placed on painting a complete gratings-picture analysis capable of describing both near and far-field phenomenon in particulate gratings.

The matrix elements of the Ω propagation matrix introduced in section 6.3.3 correspond to the "structure constants" of a KKR theory. More precisely, due to the differences between the Schrödinger equation and Maxwell's propagation equation, the Ω matrix elements can be written as a superposition of the KKR structure constants. In both KKR and particulate grating theory, one desires to calculate the lattice Green function. A fundamental method choice in this chapter is to use the language of T-matrices. Notably, we will see that the quasi-periodic Green function can be expressed as a lattice sum of multiple-scattering T-matrices. The multiplescattering T-matrices themselves are calculated in terms of the single-particle T-matrices, and the Ω matrix.

The T-matrix manipulations are carried out on a basis set of solutions to the Helmholtz equation, which we generally refer to as partial waves, (PWs), also commonly referred to as spherical wave functions (SWFs). This T-matrix approach is also generally adopted in the KKR calculations [START_REF] Papanikolaou | Conceptual improvements of the kkr method[END_REF], but in the light scattering community, the terminology "T-matrix method" is often considered to be synonymous with extended boundary condition technique (also called Null-field methods), but the T-matrix is a general theoretical construct that relates the field incident on a particle to the field scattered by the particle. As such, it can be seen as providing a complete solution to the single-particle scattering problem. In practice, the T-matrix can be generated by a wide variety of techniques including DDA, method of moments, and fictitious source techniques.

The T-matrix of an individual particle depends on the shape and the constitutive parameters of the particles, both of which can be quite arbitrary as long as the particle response is linear (including anisotropic constitutive parameters, magnetic permeability contrast etc.). However, since the T-matrix can be viewed as being the complete solution of a 1-body problem, its determination can be viewed as being separate from the grating problem. In this chapter, we simplify the T-matrix part of the problem by considering only isotropic spherical scatterers. The T-matrix of such scatterers is diagonal in the partial wave basis with its elements being determined analytically from Mie theory. We insist however, that for particulate gratings composed of more exotic scatterers, it generally suffices to insert the appropriate T-matrix to obtain the response of lattices composed of such scatterers. We refer the interested to reader to reviews of the T-matrix methods. [START_REF] Mischenko | T-matrix theory of electromagnetic scattering by particles and its applications: a comprehensive reference database[END_REF][START_REF] Mishchenko | T-matrix computations of light scattering by nonspherical particles: A review[END_REF] The methods developed in this chapter can be adapted to the study gratings composed of periodic infinite cylinders. However, there are fundamental differences in the mathematics, since this problem is usually addressed by solving 2-dimensional Helmholtz equations. We therefore neglect this problem in order to concentrate on the fully 3-dimensional problem of particulate gratings like those of figure 6.1.

The first five sections constitute the heart of this chapter since they describe the general mathematical analysis of gratings using spherical harmonic lattice sums. The last four sections treat numerical methods for calculating lattice sums and special functions. Support material, corrections, and erratum will be made available and updated during the year 2013 at www.fresnel.fr/perso/stout/index.htm.

Waves and partial waves

A fundamental aspect of the particulate gratings is that they can be viewed as a multiple-scatting phenomenon with light propagating through the host medium between individual scatterings events. The wave equation for light in this homogeneous isotropic medium is:

∇ ∇ ∇ × ∇ ∇ ∇ × E E E + k 2 E E E = 0 0 0 , (6.1)
where k = √ ε r µ r ω c is the wavenumber of the host medium. Solutions of eq.( 6.1) satisfy both the vector Helmholtz equation,

∆E E E + k 2 E E E = 0 0 0 , (6.2)
and the additional constraint that the longitudinal field components are null:

∇ ∇ ∇ • E E E = 0 . (6.3)
A basis set for solutions to the vector Helmholtz equation of eq.( 6.2) can be readily constructed starting from the scalar Helmholtz equation:

∆ϕ + k 2 ϕ = 0 . (6.4)
As well established in textbooks, eq.( 6.4) can be solved by separation of variables in spherical coordinates with 'regular' solutions taking the form of spherical harmonics, Y n,m multiplied by spherical Bessel functions, j n (kr), that are regular for all values of r. There are also 'irregular' solutions spherical Neumann functions, y n (kr), that have essential singularities as kr → 0. Details concerning the properties and calculation of the Y n,m (θ , φ ) are given in section 6.8.1.

The spherical coordinate solutions to eq.( 6.4) will henceforth be referred to as scalar partial waves and will be defined as:

J n,m (kr r r) ≡ j n (kr)Y n,m ( r r r) , and Y n,m (kr r r) ≡ y n (kr)Y n,m ( r r r) .

(6.5)

The regular partial waves, J n,m , can serve as a basis set for any source free incident field solution to eq.(6.4). Outgoing partial waves solutions of the Helmholtz equation, denoted H n,m , will be of primary interest in grating theory since they will be used to describe fields scattered by the grating. They are defined as a superposition of the regular and irregular partial waves:

H n,m (kr r r) ≡ h n (kr)Y n,m ( r r r) = J n,m (kr r r) + iY n,m (kr r r) . (6.6)
Incident field solutions to the vector Helmholtz equation of eq.( 6.2) can be expressed as scalar partial waves associated with unit vectors along each axis i.e.:

E E E inc (r r r) = x x x ∑ n,m α (x) n,m J n,m (kr r r) + y y y ∑ n ,m α (y) n ,m J n ,m (kr r r) + z z z ∑ n ,m α (z) n ,m J n ,m (kr r r) . (6.7)
The field scattered field scattered by a particle in the context of the vector Helmholtz equation can likewise be developed in terms of the outgoing spherical waves:

E E E scat (r r r) = x x x ∑ n,m β (x) n,m H n,m (kr r r) + y y y ∑ n ,m β (y) n ,m H n ,m (kr r r) + z z z ∑ n ,m β (z)
n ,m H n ,m (kr r r) , (6.8) provided that the coordinate system origin is chosen to lie inside the particle and that the field description is applied only to regions lying outside the particle. The field in eq.(6.8) represents the field scattered by a single scatterer, so the grating problem in terms of partial waves must sum over the field scattered by all the particles in the lattice. Finding efficient ways for calculating the lattice sum will therefore figure prominently in the subsequent sections of this chapter.

Before studying T-matrices in the next section, we first address an important technical issue. The field expansions in eq.(6.7) and eq.(6.8) have both transverse and longitudinal components and therefore are not generally solutions of the light propagation problem of eq.(6.1). The transverse wave condition of eq.( 6.3) can be satisfied by requiring that the Cartesian field coefficients, α (x,y,z) n,m , (and respectively β (x,y,z) n,m ) satisfy certain relations amongst themselves. The important point is to remark that the constraint conditions, although somewhat complex in spherical coordinates, only affect the partial wave coefficients, and not the partial waves themselves. Consequently, in the rest of this chapter, we can generally restrain our attention to lattice sums of scalar partial wave sums even though the end goal is describe electromagnetic field scattering.

Expressing the transverse vector partial waves in terms of the Cartesian axis partial waves of eq.(6.7) or eq.(6.8) is a relatively complex affair involving angular momentum coupling formalism, coordinate transformations, and recurrence relations. [START_REF] Newton | Scattering Theory of Waves and Particles[END_REF] Consequently, it is more common to invoke one of the various methods that have been devised over the years to directly generate transverse partial waves: Debye potentials, Hertz potentials, the Boulenkamp-Casimir approach [START_REF] Jackson | Classical Electrodynamics : Third Edition[END_REF], pilot vector techniques [START_REF] Cho | Waves and Fields in Inhomogeneous Media[END_REF], etc. Whatever one's "preferred" technique and notation, the two types of transverse partial waves, Ψ Ψ Ψ J ,q,p (often denoted M M M J ,p and N N N J ,p in the literature), can be expressed:

Ψ Ψ Ψ J ,1,p (kr r r) ≡ j n (kr) X X X n,m ( r r r) Ψ Ψ Ψ J ,2,p (kr r r) ≡ 1 kr n (n + 1) j n (kr)Y Y Y n,m ( r r r) + [kr j n (kr)] Z Z Z n,m ( r r r) , (6.9)
where X X X n,m , Y Y Y n,m , and Z Z Z n,m are the vector spherical harmonics (VSHs), (described in section eq.( 6.9.3)). The first subscript, J , on Ψ Ψ Ψ J ,q,p serves to indicate that the radial dependence is governed by spherical Bessel functions. A value of q = 1 in the second subscript indicates a "transverse electric" (TE) type wave (i.e. possessing no radial electric field component), while q = 2 indicates transverse magnetic (TM) waves. In order to minimize the number of subscripts, we adopt the common procedure that the third subscript p of Ψ Ψ Ψ J ,q,p , replaces the two multipole subscripts n and m by defining its value such that [START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF]:

p (n, m) ≡ n(n + 1) -m . (6.10)
With the notation of eq.( 6.9), one can express any incident field satisfying equation (6.1) in terms of the transverse vector partial waves:

E E E inc (r r r) = E ∑ q=1,2 ∞ ∑ p=1
Ψ Ψ Ψ J ,q,p (kr r r) a q,p , (

where a q,p are (dimensionless) field coefficients, and E is a constant with the dimension of electric field and which can be used to adjust the field strength. With this notation, the field scattered by a particle whose circumscribing sphere is centered at a position x x x j can be written:

E E E scat r r r j = E ∑ q=1,2 ∞ ∑ p=1 Ψ Ψ Ψ H ,q,p kr r r j f ( j) q,p , (6.12)
where r r r j ≡ r r rx x x j , and f

q,p are the scattering coefficients of the particle j. The index, H . on the Ψ Ψ Ψ H ,q,p indicates that the radial dependence should be governed by spherical Hankel functions, h n (kr), rather than the spherical Bessel functions, j n (kr), found in the Ψ Ψ Ψ J ,q,p functions of eq.(6.9).

T-matrix theory

Green functions and T-matrices

The fundamental object that one would like to calculate in a multiple-scattering system (like a particulate grating) is the system Green's function. However, the dyadic Green's function for a homogeneous medium has a strongly singular behavior and needs to be defined in the context of distributions. [START_REF] Cho | Waves and Fields in Inhomogeneous Media[END_REF] The T-matrix formalism allows us to largely circumvent this singular behavior, and also to work directly in terms of fields which is often more manageable than the relatively intricate dyadic Green's function formalism. For instance, the operator form of the Green function of a single object in a homogeneous medium can be written: G G G = g g g + g g gt t tg g g , (

where t t t is the isolated particle (or 1-body) T-matrix operator, and g g g is the Green function operator of the homogeneous medium (sometimes called a propagator). In this formalism, the singular behavior is relegated to the propagator, g g g, leaving the (non-singular) scattering response due to the object being described by t t t. Furthermore, when considering excitations outside the scatterer, the homogeneous Green function, g g g to the right of the t t t operator acting on the sources generates the incident field, while the g g g to the left of it generates the scattered field. [START_REF] Sheng | Introduction to Wave Scattering, Localization and Mesoscopic Phenomena[END_REF] In the partial wave basis, t t t then truly takes the form of a matrix, henceforth denoted, t, that relates the incident field coefficients to the scattered field coefficients:

f = ta , (6.14)
where a and f are column matrices composed respectively of the incident field and scattered field coefficients (cf. eqs. (6.11 and (6.12)). [START_REF] Stout | Recursive t-matrix algorithm for resonant multiple scattering: Applications to localized plasmon excitations[END_REF] The 1-body T-matrix, t, in this expression is now truly a matrix relating field coefficients of partial wave field decompositions. This T-matrix formalism can be extended to include systems containing N particles. The system Green function can be written,

G G G = g g g + g g g N ∑ j=1 T T T ( j) g g g , (6.15)
where the multiple-scattering (or N-body) T-matrix operators, T T T ( j) , are associated with each particle and which incorporates all the multiple-scattering effects due to the presence of the N -1 other particles in the system. Passing once again to a partial wave field description, the multiple-scattering T-matrix, T ( j) , generates the field scattered by each particle in terms of the field incident on the system:

f ( j) = T ( j) a ( j) , (6.16)
where a ( j) indicates the incident field developed on a coordinate system centered on the j th particle. All multiple-scattering phenomenon and some rather subtle technical difficulties have all been incorporated into the definition of T ( j) , but nowadays they can be calculated rather readily for systems with a finite number of particles starting from the 1-body T-matrices, t ( j) , of the individual particles. [START_REF] Stout | Recursive t-matrix algorithm for resonant multiple scattering: Applications to localized plasmon excitations[END_REF] The number of particles in a grating problem is infinite (from the ideal mathematical standpoint), which given their physical content would render exact calculations of T ( j) impossible. Nevertheless, the fact that the system is identical when viewed from any lattice site allows the T ( j) matrices to be calculated as a lattice sum as we shall see in section 6.3.3. We first rapidly review below our notation and terminology for lattices.

Direct and reciprocal lattices

A lattice, Λ, of dimension d Λ , is invariant under a coordinate system translation along any vector, r r r j j j , that can be expressed r r r j j j =

d Λ ∑ i=1 j i a a a i , (6.17)
where a a a i are the primitive lattice vectors, and j j j ≡ ( j 1 , ..., j d Λ ) is a shorthand notation for a set of d Λ relative integers, j i ∈ Z. In order to diminish the number of subscripts, we will sometimes employ an alternative notation for the primitive lattice vectors: a a a ≡ ≡ ≡ a a a 1 , b b b ≡ ≡ ≡ a a a 2 , and c c c ≡ ≡ ≡ a a a 3 . It also proves convenient to define the x and y axis of the system so that the primitive lattice vectors can be expressed: When d Λ = 3, the r r r j j j ensemble defines a crystalline type lattice, henceforth denoted (L), as frequently encountered in photonic crystals and "meta-materials". A two-dimensional grating, or mono-layer lattice (ML), like that of figure 6.1, occurs when the system invariance only occurs for 2D displacements of r r r j j j = j a a a a + j b b b b. Finally, linear chains (C) are only invariant with respect to translations of r r r j = ja a a.

The reciprocal lattice, Λ * , is defined in terms of lattice 'wave-vectors', p p p g g g , defined in terms of the primitive reciprocal lattice vectors, a a a i :

p p p g g g = 2π d Λ ∑ i=1 g i a a a i , (6.18)
where g i ∈ Z. The primitive reciprocal vectors, a a a j , are defined such that their scalar products with respect to a a a j satisfy: a a a i • a a a j = δ i j i, j = 1, ..., d Λ . (

From eqs.(6.17) (6.18) and (6.19), one readily finds that the reciprocal lattice vectors, p p p g g g , of eq.(6.18), have the property r r r j j j • p p p g g g = 2πN , (

where N is some integer (which results in exp ir r r j j j • p p p g g g = 1 for all r r r j j j and p p p g g g ).

The unit cell "volume", A , appears repeatedly in theories of particulate lattices. For lattice dimensions of d Λ = 1, 2, and 3, the corresponding A 1,2,3 is given by:

   A 1 = |a a a| d Λ = 1 A 2 = |a a a × × × b b b| d Λ = 2 A 3 = |(a a a × × × b b b) • c c c| d Λ = 3 , (6.21)
with dimensions of "length" for d Λ = 1, "area" for d Λ = 2 and "volume" for d Λ = 3. The corresponding "volume" of the reciprocal space lattice sites are given by A -1 .

Grating T-matrices

Each site of a lattice is identical to all the others so that the multiple scattering T-matrices of eq.(6.16) are all equal, i.e. T ( j) = T . The scattering coefficients f ( j) are then given by:

f ( j) = Ta ( j) . (6.22)
The trouble with this equation is that the coefficients f ( j) and a ( j) are expressed on localized partial wave basis, but the long range nature of scattered fields would require the T -matrices to act on very high multipole orders in order to account for these long-rang interactions. Since manipulating high multipole orders is numerically inefficient, one considerably simplifies this problem by only calculating the multiple-scattering T -matrices for incident fields satisfying a quasi-periodic condition. Quasi-periodicity can be viewed as requiring the partialwave decomposition of the incident field on each lattice site, r r r j j j , to satisfy, a ( j) = e iβ β β •r r r j j j a , (

where 'a' corresponds to the incident field coefficients at the origin, and β β β , the 'on-shell' quasiperiodicity vector. The term 'on-shell' indicates that the quasi-periodic vector satisfy |β β β | = k, and although this considerably restricts the type of incident fields (and modes) that one can treat, it does correspond to many situations of physical interest including both homogeneous and inhomogeneous plane-waves incident on diffraction gratings.

The quasi-periodic condition can be viewed as a partial Fourier transform description in that the overall field behavior is of an oscillatory nature, while the quasi-periodic T-matrix, T β β β , describes local-field perturbations due to the presence of the particles. Consequently, one expects the T β β β matrices to be well approximated on a truncated (i.e. finite) partial-wave basis (similar to the behavior of the isolated particle T-matrices [START_REF] Stout | Recursive t-matrix algorithm for resonant multiple scattering: Applications to localized plasmon excitations[END_REF]). The quasi-periodic condition allows the Foldy-Lax equations for the multiple-scattering T-matrices to take the form:

T β β β = t + tΩ β β β T β β β , (6.24)
where the Ω β β β matrix designates a quasi-periodic lattice sum of the irregular translation-addition matrices:

Ω β β β (k, a a a i≤d ) = ∑
r r r j j j ∈Λ r r r j j j =0 0 0 e iβ β β •r r r j j j H kr r r j j j . (

The analytical properties of the irregular translation-addition matrix, H (x x x), are described in section 6.8.3 where one also gives expressions for its matrix elements. The exclusion of the 'origin' lattice site, r r r j j j = 0 0 0, from the sum in Ω β β β has a physical significance in that it accounts for propagation of the light scattered by all the other particles in the lattice onto the particle at the origin (the light 'scattered' by the particle onto itself has already been included in the individual T-matrix, t). One finds in section 6.8.3 that each matrix element of the translation-addition matrix, H kr r r j j j , can be written: H kr r r j j j p,q;p ,q = ∑ l,m C l,m p, q; p , q h l kr j j j Y l,m r r r j j j , (

where the sum over the multipole indices, (l, m) is finite. Expressions for the C l,m (p, q; p , q ) coefficients [START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF][START_REF] Cho | Waves and Fields in Inhomogeneous Media[END_REF][START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF] are given in the section 6.8. Inserting eq.( 6.26) into eq.( 6.25) and rearranging the summations, we find

Ω β β β p,q;p ,q = ∑ l,m
C l,m p, q; p , q ∑ r r r j j j ∈Λ r r r j j j =0 0 0 e iβ β β •r r r j j j h l kr j j j Y l,m r r r j j j

≡ ∑ l,m
C l,m p, q; p , q S l,m (k, β β β ) . (6.27)

where we have defined S l,m (k, β β β ) as a Hankel function lattice sum such that:

S n,m (k, β β β ) ≡ S H n,m (k, β β β ) ≡ ∑
r r r j j j ∈Λ r r r j j j =0 0 0 e iβ β β •r r r j j j H n,m kr r r j j j . (6.28)

We recall that H n,m (x x x) was defined in eq.( 6.6) as a partial wave of the spherical Hankel function type.

It will occasionally prove useful to calculate the analogous lattice sums over the partial waves of the Bessel or Neumann types, denoted respectively, S J n,m (k, β β β ) and S Y n,m (k, β β β ). Since we will principally be concerned with the partial wave lattice sums of the Hankel function type, S n,m without a superscript will always indicate a lattice sum of the Hankel function type. We also remark that the exclusion of the origin position from the lattice sum is important from a mathematical standpoint since the Hankel functions have an essential singularity at their origin.

The solution to eq.( 6.24) for the multiple-scattering T-matrix is readily formulated in terms of matrix inversion:

T β β β = t -1 -Ω β β β -1 . (6.29)
Once the T β β β matrix is known, the scattering field coefficients for any particle, j, in the lattice is the same as the coefficients at the origin but multiplied by a e iβ β β •r r r j j j phase factor. In the matrix notation, this is simply expressed:

f ( j) β β β = e iβ β β •r r r j j j f β β β = e iβ β β •r r r j j j T β β β a , (6.30)
where a is the column matrix composed of the incident field coefficients developed around the origin.

Far-fields

The field 'scattered' by the grating (i.e. 'transmitted' and 'reflected' diffraction orders) can be determined by inserting eq.(6.30) into eq.(6.12) wherein the scattered field takes the form of a lattice sum of the transverse-outgoing-vector partial waves, Ψ Ψ Ψ H ,q,p , described in eq.(6.9) of section 6.2:

E E E s,Λ (r r r) = E ∑ r r r j j j ∈Λ e iβ β β •r r r j j j Ψ Ψ Ψ H kr r r j j j f β β β ≡ E ∑ q=1,2 ∞ ∑ p=1 ∑ r r r j j j ∈Λ
Ψ Ψ Ψ H ,q,p kr r r j j j e iβ β β •r r r j j j f β β β q,p . (6.31)

We will see in eq.( 6.62) of section 6.4.4 that for each multipole order, p = 1, ...∞, and transverse wave type q = 1, 2; the term in brackets can be re-expressed as an infinite sum of plane waves. Only a finite subset of these waves are of the propagative type however (the rest are all of an evanescent nature). Consequently, the multipole summation of eq.( 6.31) allows one to calculate the efficiency of each reflected or transmitted order in the far field.

Near-fields

Another quantity of physical interest is that of near fields in a particulate grating (non-linear effects, SERS, etc.). The plane wave expansion discussed above for far fields could be invoked in principal, but for near fields one must also calculate the (infinite) evanescent orders that one could neglect in the far field. The convergence of the plane wave expansion will generally be poor near the grating, which renders this approach unattractive.

As long as the incident field is quasi-periodic with respect to the grating, one needs only to determine the near fields in a single Brillouin zone around a given lattice site (the site at the origin being the most practical). In this case, it seems clear that the localized multipolar field developments are well adapted to the development of the local field in the Brillouin zone. In multiple scattering theory, the f β β β coefficients give the field scattered by the particle at the origin, while the excitation field corresponds to the field at that was 'incident' on this particle, i.e. the superposition of the field incident on the grating and the field scattered by all the other particles in the system. This excitation field can be developed on the regular partial waves and its coefficients, e β β β , related to the scattering coefficients via the 1-body T-matrix via the relation:

e β β β = t -1 f β β β . (6.32)
The total field in the Brillouin zone is simply a superposition of the scattered and excitation field:

E E E (B.z.) t (r r r) = E Ψ Ψ Ψ H (kr r r) f β β β + Ψ Ψ Ψ J (kr r r)t -1 f β β β ≡ E ∑ q=1,2 ∞ ∑ p=1
Ψ Ψ Ψ H ,q,p (kr r r) f q,p + Ψ Ψ Ψ J ,q,p (kr r r) e q,p . (6.33)

Propagating modes

The response of a grating or chain of particles will often be dominated by the excitation of nearly guided modes in the structure (commonly referred to as leaky modes). Unlike, lossless 3D lattices which have true propagating modes described by real values of β β β , the leaky modes in a grating or chain of particles are so-named because they generally radiate energy to the far field even in the absence of intrinsic losses. Consequently, leaky modes can be described by a complex valued β β β -vector (or complex frequency). The determination of a leaky mode thus involves searching for a complex pole in the determinant of the multiple-scattering T-matrix, T β β β . Since matrix inversions are numerically expensive, one may prefer to look for zero eigenvectors, ν α , of the matrix t

-1 -Ω β β β α , i.e. t -1 -Ω β β β α ν α = 0 . (6.34)
However, the search for zero eigenvalues can limit the implantation of complex analysis methods that have proven useful in determining the position of poles in the complex plane. The Floquet mode associated with the eigenvector, ν α , can be constructed from eq.( 6.34) coupled with the plane-wave development the terms in eigenvector, ν α :

E E E (F.m.) α (r r r) = E ∑ q=1,2 ∞ ∑ p=1 ∑ r r r j j j ∈Λ
Ψ Ψ Ψ H ,q,p kr r r j j j e iβ β β α •r r r j j j [ν α ] q,p .

(6.35)

Before finishing this section, it should be pointed out that matrix inversion solutions to the multiple scattering problem (like that given in eq.(6.29)) were disregarded for a long time in favor of iterative solutions to the T-matrix or underlying linear system of equations. The reason for this is that the matrix t -1 -Ω β β β α is generally ill-conditioned. This difficulty can be generally overcome by analytical matrix balancing as described in the next section 6.3.4.

Matrix balancing

Although not necessary from a formal standpoint, analytical matrix balancing improves the conditioning of the matrices occurring in multiple-scattering calculations for both matrix inversion and eigenvalue resolution. [START_REF] Stout | Recursive t-matrix algorithm for resonant multiple scattering: Applications to localized plasmon excitations[END_REF] Analytical matrix balancing can be achieved by multiplying a matrix from both the right and left by diagonal matrices, [ξ ] and [ψ] -1 , whose matrix elements are given by:

[ψ] q,q ,p,p = δ q,q δ p,p ψ n (kR) ,

[ξ ] q,q ,p,p = δ q,q δ p,p ξ n (kR) , (6.36)

where ψ n (kR) and ξ n (kR) are respectively the regular and irregular spherical Ricatti-Bessel functions (cf. 6.168) and R the radius of the minimal circumscribing sphere surrounding the scatterers.

Matrix balancing can be readily formulated by defining normalized incident and scattering coefficients, a and f respectively such that:

a ≡ [ψ] a , f β β β ≡ [ξ ] f β β β . (6.37)
The associated normalized or 'balanced' matrices are defined [START_REF] Stout | Recursive t-matrix algorithm for resonant multiple scattering: Applications to localized plasmon excitations[END_REF]:

t ≡ [ξ ]t [ψ] -1 , T β β β ≡ [ξ ] T β β β [ψ] -1 , Ω β β β ≡ [ψ] Ω β β β [ξ ] -1 . (6.38)
The above definitions were chosen such that eqs.(6.22) and (6.24) respectively take the form: The normalized T-matrix, T , is then obtained via the generally well-conditioned matrix inversion:

f ( j) ≡ T β β β a ( j) , (6.39)
T β β β = t -1 -Ω β β β -1 . (6.41)
Since we generally want the non-normalized T-matrix for applications, we reconstruct, T β β β via a final multiplication by our diagonal matrices:

T β β β = [ξ ] -1 T β β β [ψ] . (6.42)

Mathematical relations for lattice sums

This section is dedicated to reviewing the mathematical relations that allow one to treat lattice sums for lattices of dimensions d Λ = 1, 2, 3 (i.e. particulate chains, gratings, and crystals). They will notably allow us to evaluate the lattice sum in eq.(6.27) which is used to calculate far-field response from gratings. These relations were derived (and often rederived) in many places, and we refer the reader to refs. [START_REF] Linton | Lattice sums for the helmholtz equation[END_REF][START_REF] Linton | One-and two-dimensional lattice sums for the threedimensional helmholtz equation[END_REF][START_REF] Moroz | Exponentially convergent lattice sums[END_REF][START_REF] Moroz | Quasi-periodic green's functions of the helmholtz and laplace equations[END_REF][START_REF] Enoch | Sums of spherical waves for lattices, layers and lines[END_REF] for additional details and perspectives.

The most difficult mathematical problem to address will be the evaluation of Hankel function lattice sum, S H n,m that was introduced in eq.(6.28) for the calculation of the Ω β β β matrix of eq. (6.27).

S H n,m ≡ ∑ r r r j ∈Λ r r r j =0 0 0 e iβ β β •r r r j j j H n,m kr r r j j j (6.43)

When not otherwise specified, partial waves lattice sums will always be assumed to be of the Hankel functions type. The underlying reason for this appears in the translation-addition where Hankel functions allow one to re-express waves scattered by a given lattice site as waves incident on a different lattice site. We will occasionally consider Bessel and Neumann types of scalar partial waves:

S Y n,m ≡ ∑ r r r j ∈Λ r r r j =0 0 0
e iβ β β •r r r j j j Y n,m kr r r j j j S J n,m ≡ ∑ r r r j ∈Λ r r r j =0 0 0 e iβ β β •r r r j j j J n,m kr r r j j j , (6.44)

The interest of these sums in part is due to the fact that

S H n,m = S J n,m + iS Y n,m but also because S J n,m
can potentially prove useful in certain applications. We will see that the relations developed in this chapter permit the S J n,m sum to be evaluated in closed form, but unfortunately the Neumann partial wave sum, S Y n,m , appears to be as difficult to evaluate as the Hankel partial wave sum, and a closed form expression does not appear to be possible.

Lattice reduction

Although Ewald sums are a time honored technique in solid state physics, a considerable amount of effort has recently been devoted to what has come to be called Lattice reduction techniques. The basic idea turns around the fact that lattice sums tend to be more practical for d Λ = 1 and d Λ = 3 than for the grating dimension of d Λ = 2.

First one chooses a coordinate system such that a preferred axis (like the z axis) will along a given lattice vector. For example a a a = (0, 0, a), and b b b = (0, b 2 , b 1 ). With this basis the 2D Mono-Layer lattice sum, S ML n,m can then be written expressed as a superposition of a Chain sum, S C n,m , containing the origin (z-axis), and a superposition of all the chain sums 'above' the central chain (z>0), denoted S ML+ n,m or 'below' the central chain, S ML- n,m (z<0). The central chain can be readily be evaluated using one of the techniques described in this chapter, while the integral expression for Hankel functions described in section 6.4.4 allows one to derive efficient expressions for S ML+ n,m and S ML- n,m . Lattice reduction can also be applied in the reverse direction, with one expressing the 3D crystalline lattice sum, S L n,m , as the superposition of a monolayer sum in the z = 0 plane, with sums of all monolayers with z > 0, S L+ n,m and all monolayers with z < 0, S L- n,m . There exists efficient techniques for calculating the crystalline lattice sum, S L n,m while one can determine efficient expressions for S L,+ n,m and S L,- n,m using again the integral expressions of section 6.4.4. In this reverse lattice sum method, the monolayer lattice sum is expressed:

S ML n,m = S L n,m -S ML,+ n,m -S L,- n,m . (6.45)
One should that the choice of orientation of the coordinate axis will not be the same in general for different lattice sum techniques, but these differences can be be compensated for by using the rotation matrices of section 6.8.4

Lattice reduction is based on the idea that it can prove numerically efficient to carry out lattice sums for a lattice dimensions other than that desired. To construct a 3D periodic media, we rotate the 2D lattice of the preceding section back to an orientation in the xOy plane with a a a = (a, 0, 0), and b b b = (b 1 , b 2 , 0), and now c c c = (0, c 2 , c 3 ). The quasi-periodic vector is given by β β β = (β 1 , β 2 , β 3 ). The 3D lattice sum can then be written:

S L n,m = S ML n,m + S L+ n,m + S L- n,m . (6.46)
The S L+ n,m denotes all the z > 0 planes, while S L- n,m sums all the z < 0 planes. The lattice reduction technique breaks the sum down into elements which tend to have a decreasing difficulties for divergence. An interest of the lattice reduction technique is that it can be adapted to partial lattices. For instance, large but finite chains, a finite number of infinite chains or finally a finite number of infinite planes.

Plane wave expansion

The expansion of a plane wave in terms of partial waves allows one to transform between partial wave and Fourier transforms. It reads:

e ik k k•r r r = 4π ∞ ∑ ν=0 µ=ν ∑ µ=-ν i ν j ν (kr)Y * ν,µ k k k Y ν,µ ( r r r) = ∞ ∑ ν=0 µ=ν ∑ µ=-ν p ν,µ Ψ ν,µ (r r r) , (6.47)
where Ψ ν,µ (r r r) are the scalar partial wave functions discussed in section 6.2, and p ν,µ the coefficients in the development of a scalar plane wave on a partial wave basis i.e. :

Ψ ν,µ (r r r) ≡ j ν (kr)Y ν,µ ( r r r) , p n,m = 4πi n Y * n,m k k k . (6.48)
One can produce an integral expression of j n (kr)Y n,m ( r r r) by multiplying both sides of eq.( 6.47) by Y n,m k k k and integrating over all directions of k k k.

dΩ k k k e ik k k•r r r Y n,m k k k = 4π dΩ k k k ∞ ∑ ν=0 µ=ν ∑ µ=-ν i ν j ν (kr)Y ν,µ ( r r r)Y * ν,µ k k k Y n,m k k k = 4πi n j n (kr)Y n,m ( r r r) . (6.49)
We have thus found that regular partial waves are an angular Fourier transform of the spherical harmonics:

Ψ J ,n,m ≡ j n (kr)Y n,m ( r r r) = 1 4πi n dΩ k k k e ik k k•r r r Y n,m k k k . (6.50)
Likewise, the transverse regular partial waves, Ψ Ψ Ψ J can be expressed as an angular Fourier transform of the vector spherical harmonics:

Ψ Ψ Ψ J ,q=1,n,m (kr r r) = i -n 4π dΩ k k k e ik k k• • •r r r X X X n,m ( k k k) Ψ Ψ Ψ J ,q=2,n,m (kr r r) = i 1-n 4π dΩ k k k e ik k k• • •r r r Z Z Z n,m (k k k) . ( 6
.51)

Poisson summation formula

The Poisson summation formula is a crucial mathematical tool for evaluating lattice sums. It allows one to pass from a sum over the real lattice vectors to a sum over the reciprocal lattice vectors. Formally, it can be written:

∞ ∑ r r r j j j ∈Λ e ik k k•r r r j j j = (2π) d Λ A d Λ ∑ p p p g g g ∈Λ * δ k k k -p p p g g g , (6.52)
where A d Λ is the "volume" of the reciprocal lattice cell. Since long and short range interactions can both be strong for lattice problems the Poisson summation formula often does not directly accelerate the lattice sum, but it nevertheless proves invaluable for a number of useful that we will derive in the rest of this chapter.

For the 1-D sum in eq.( 6.90), this can be written:

∞ ∑ j=-∞ e i(k k k+β β β )•( z z z ja) = 2π a ∞ ∑ g=-∞ δ k z + β z - 2π a g . (6.53)
We then write this relation in a dimensionless form:

∞ ∑ j=-∞ e i(k k k+β β β )•( z z z ja) = 2π kai n ∞ ∑ g=-∞ δ k z k + β z k -g 2π ka .
(6.54)

Integral expressions for outgoing partial waves

The Weyl identity expressed the Hankel function of order 0 as an integral of plane waves:

h 0 (kr) = 1 2πk ∞ -∞ dk x dk y exp (±ik k k • • • r r r) k z = 1 2πk ∞ -∞ dk x dk y exp [±i (k x x + ik y y + ik z z)] k z z ≷ 0 , (6.55)
where the plus sign is taken for z > 0 and the minus sign is used when z < 0. The k z component is fixed by the constraint that

k 2 x + k 2 y + k 2 z = k 2 , namely k z = k 2 -k 2 x -k 2 y .
It is interesting to remark that the spherical Bessel function is a superposition of plane waves that are constrained to satisfy k k k = k. Since the reciprocal space integration in eq.( 6.55) is carried out in the xOy plane, it is convenient to define a specific symbol for the wavevector in the xOy plane, K K K = k x x x x + k y y y y, and the full wavevector is then, k k k = K K K + k z z z z. It is also convenient to define dimensionless or normalized quantities:

K ≡ |K K K| /k γ z ≡ k z /k = 1 -K 2 . (6.56)
If we take the position vector r r r in eq.( 6.55) to lie along the z axis, r r r = = = r z z z, then we can integrate over the azimuthal angle to obtain a single integral expression for Hankel functions that can be used in lattice sums:

h 0 (kr) = ∞ 0 dKK exp [±iγ z kr] γ z z ≷ 0 . (6.57)
Wittmann pointed out that the above Weyl identity of eq.( 6.55) can be generalized to all partial waves of the Hankel function type [START_REF] Wittmann | Spherical wave operators and the translation formulas[END_REF] :

Ψ H ,n,m ≡ h n (kr)Y n,m ( r r r) = i -n 2π ∞ -∞ dk x dk y k x + ik y m P m n (γ z ) exp (±i (k x x + k y y + k z z)) γ z z ≷ 0 . (6.58)
If we take r r r again to lie along the z axis, we find an integral expression for spherical Hankel functions:

h n (kr) = i -n ∞ 0 dKK P n (γ z ) exp [iγ z kr] γ z . (6.59)
The integral relation of eq.( 6.58) can also be extended to the outgoing vector partial waves:

Ψ Ψ Ψ H ,q=1,n,m (kr r r) = i -n 2π ∞ -∞ dk x dk y exp(±i(k x x+k y y+k z z)) γ z X X X n,m ( k k k) Ψ Ψ Ψ H ,q=2,n,m (kr r r) = i 1-n 2π ∞ -∞ dk x dk y exp(±i(k x x+k y y+k z z)) γ z Z Z Z n,m (k k k) z ≷ 0 . (6.60)
The Poisson sum rule allows one to express quasi-periodic 2D lattice sum in terms of 2D reciprocal lattice vectors. For the scalar partial waves, one has: ∑ r r r j j j ∈Λ exp iβ β β • • • r r r j j j Ψ H ,n,m kr j j j = ∑ p p p g g g ∈Λ * 2πi -n kk + g g g,z A 2 Y n,m k k k ± g g g exp ik k k ± g g g • r r r z ≷ 0 , (6.61) while for the vector partial waves, ∑ r r r j j j ∈Λ exp iβ β β • • • r r r j j j Ψ Ψ Ψ H ,1,n,m (kr r r j j j ) = ∑ p p p g g g ∈Λ *

2πi -n kk + g g g,z A 2 X X X n,m k k k ± g g g exp ik k k ± g g g • r r r ∑ r r r j j j ∈Λ exp iβ β β • • • r r r j j j Ψ Ψ Ψ H ,2,n,m (kr r r j j j ) = ∑ p p p g g g ∈Λ * 2πi 1-n kk + g g g,z A 2 Z Z Z n,m k k k ± g g g exp ik k k ± g g g • r r r z ≷ 0 . (6.62)
In the partial wave lattice sums of eqs.(6.61) and (6.62), the wavevector k k k ± g g g is given by:

k k k ± g g g ≡ β β β + p p p g g g ± z z z k 2 -β β β + p p p g g g 2 , (6.63)
and k + g g g z is its z component:

k + g g g,z ≡ k k k ± g g g • z z z = k 2 -β β β + p p p g g g 2 .
(6.64)

One remarks that for a real Bloch vector β β β , the wavevector k k k ± g g g is real i.e. propagative in nature only for those lattice vectors for which k > β β β + p p p g g g .

(6.65)

Partial wave rotation

Let us consider once again a row matrix, Ψ t J ,Y ,H , composed of one of the types of partial waves in eq.(6.5). The fact that the choice of the orientation of the coordinate system is arbitrary imposes transformation relations amongst the partial waves. Let us consider a position M given by the vector r r r in our chosen coordinate system. We next consider another coordinate system with the same origin, but rotated by the 3 Euler angles, α, β , and γ in which the same point M is now designated by a vector r r r (n.b. |r r r | = |r r r| = r). The linear relationship between the row matrix in these 2 coordinate systems is then:

Ψ t (r r r) = Ψ t r r r D (α, β , γ) .
(6.66)

If the rotated coordinate systems is taken such that r r r lies along the z axis in the rotated coordinate n this relation takes the form:

Ψ t (r r r) = Ψ t (r z z z) D (φ , θ , 0) . (6.67)
In component form this reads for Hankel function sums:

h n (kr)Y n,m ( r r r) = h n (kr)Y n0 (0, 0) D n,0;n,m (φ r r r , θ r r r , 0) = 2n + 1 4π h n (kr) D n,0;n,m (φ r r r , θ r r r , 0) , (6.68)
where we used eq.( 6.75) for an expression of the Y n0 (0, 0).

Numerical Examples

This section will be expanded considerably in the second edition.

Far and near field response from gratings

As discussed in section 6.3, once the lattice sums have been determined for all the Ω β β β matrix elements, and the lattice T-matrix of eq.(6.29) obtained, one has ready access to both the far and near field response of the system. However, the quasi-periodic lattice for all the multipole orders must be calculated anew whenever one looks for response to a different quasi-periodicity vector, β β β or wavenumber k (i.e. frequency).

Modes for particulate chains

There is considerable interest in calculating and characterizing the 'propagating' modes of periodic chains, gratings, and finite stacks of particulate gratings. A major difficulty is that these systems are open, so propagating 'modes' are necessarily lossy which entails that they don't exist for real values of frequency and quasi-periodicity vector β β β . This situation was in fact already present in the T-matrices of the individual particles which have do have poles at any real frequencies, but one can find poles at complex frequencies. For planar surfaces, Greffet has recently argued [START_REF] Archambault | Surface plasmon fourier optics[END_REF] that the Leaky-modes can be described by letting either frequency or wavevector be described by a complex number. This idea has recently been employed by several authors for calculating modes in infinite particulate chains where the component of β β β along the chain axis is allowed to be a complex number. [START_REF] Rolly | Dispersion relations in metal nanoparticle chains: necessity of the multipole approach[END_REF][START_REF] Campione | Complex bound and leaky modes in chains of plasmonic nanospheres[END_REF][START_REF] Conforti | Dispersive properties of linear chains of lossy metal nanoparticles[END_REF][START_REF] Koenderink | Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains[END_REF][START_REF] Koenderink | Plasmon nanoparticle array waveguides for single photon and single plasmon sources[END_REF] There are still some technical issues under discussion considering modes in gratings due to ambiguities in the best manner to define a complex propagation vector in the grating plane. [START_REF] Fructos | Characterization of complex plasmonic modes in two-dimensional periodic arrays of metal nanospheres[END_REF]. This issue will be treated in more depth in the second edition.

Typically, one has looked for propagating modes in particulate arrays of sub-wavelength particles metallic particles. There is however an increased in interest in high index dielectrics. Due to the complexity of the full multipole approach, most works searches for modes in the complex plane have adopted what amounts to be a electric dipole approximation to eq.(6.29). [START_REF] Rolly | Dispersion relations in metal nanoparticle chains: necessity of the multipole approach[END_REF] We have recently argued that electric dipole is insufficient in the presence of strong interactions that are provoked by resonances. [START_REF] Rolly | Dispersion relations in metal nanoparticle chains: necessity of the multipole approach[END_REF] These results and conclusions are reviewed here.

We adopt the same parameters for a plasmonic chain as Conforti and Guasoni. [START_REF] Conforti | Dispersive properties of linear chains of lossy metal nanoparticles[END_REF] Namely, we consider an infinite chain of identical 50nm diameter silver particles separated by d = 75nm (center-to-center). The system is immersed in a non-magnetic medium with relative permittivity ε = 2.25 (n = 1.5).

The figures are plotted with normalized frequencies and wave-vectors:

ω ≡ ωd 2πc = d λ v β ≡ β d 2π (6.69)
where λ v is the vacuum wavelength. The light line for these parameters is given by ω = β n med . The dispersion relations of the principal propagating modes calculated in the electric dipole approximation are plotted in figure 6.2 (dashed curves). They are then compared with fully converged n max = 10 calculations of these dispersion relations (solid line) in this same figure by solving eq.(6.34). The imaginary part of the dispersion relations for dipolar and converged multipole calculations are given in figure 6.2a). and 6.2b) merit some commentary. It is immediately clear that the dipole approximation provides a moderately accurate prediction of dispersion relations only over a narrow range of frequencies for which the imaginary part of the propagating wavevector is rather small, and the real part is near the light line. One should also recall that symmetry dictates that if a given value of β corresponds to mode at a given frequency, then by symmetry, -β is also a solution to these equations. For the sake of clarity, these symmetric modes are not presented in these figures.

Like Conforti and Guasoni [START_REF] Conforti | Dispersive properties of linear chains of lossy metal nanoparticles[END_REF], we find a transverse mode, labeled "T2" whose imaginary part of β is opposite in sign with the real part of β . It may prove physically relevant to think of this T2 mode as a backscattering mode, or to interpret this in terms of negative effective index. It is interesting to remark that the T2 mode tends toward the edge of the Brillouin zone at low frequencies. It has recently been argued that these Our dipole approximation predictions for the longitudinal mode are quite similar to that of ref. [START_REF] Conforti | Dispersive properties of linear chains of lossy metal nanoparticles[END_REF] wherein the dipole prediction is that the mode "folds back" before reaching the edge of the Brillouin zone. The full multipole calculations on the other hand predict that the longitudinal mode goes to the edge of the Brillouin zone, and that the "fold back" only occurs after it has gone "beyond" the edge of the Brillouin zone. In our calculations, the "T1" mode is quite close to the light line, and henceforth rather poorly confined by the plasmon chain so its importance in applications seems limited. In our calculations, the dipole approximation for the "T1" mode is quite similar to the multipole solution except that we only found that the full multipole solution predicted both extremities of the T1 mode to lie on the light line. Due to the system design (sub-wavelength resonance particles) one expected to find significant guiding of modes appear only in the frequency domains where the scattering cross section of the individual particles is non-negligible. To illustrate this point, we plot the extinction and scattering cross section for an individual particle in the chain in figure 6.3. We remark in particular that near individual particle resonance maxima, all the guided modes of figure 6.2 lie near the light line, and as seen in figure 6.2b) it is here also that their imaginary parts are smallest. Furthermore, with the exception of the 'backscattering' mode T2, all guided modes apparently cease to exist when one moves sufficiently far away from the scattering resonance frequency.

The reader has probably remarked some strange behavior of the modes in the electric dipole approximation at high frequencies. For instance, at around ω = 0.225 a "kink" appears in the longitudinal mode, and a spurious T2 solution emerges from the light line. We carried out mode calculations with various multipole cutoffs and found that such kinks and spurious solutions were relatively commonplace (at high or low frequencies) when low numbers of multipoles are used in the simulations and such behavior disappears when higher multipole orders are used. It is also worth remarking that for high order simulations, the Re[β ] of the modes terminate at either the light line, or the edge of the Brillouin zone, but modes can terminate at undiscriminating positions in β space when calculations are carried out at low order. The mode diagrams of figures 6.2a) and 6.2b) were somewhat unconventional since they did not display symmetric, -β , modes, and allowed the dispersion relation of the longitudinal mode to move outside the Brillouin zone. A more conventional representation of the dispersion relations is given in figure 6.4 which includes the symmetric modes, but only displays modes when Re[β ] has positive values lying within Brillouin zone (here we display only the converged multipole calculations). Transverse modes with Im[β ] > 0 modes are given by solid blue(black) lines, while transverse modes with Im[β ] < 0 are dashed dashed blue(black) lines. Longitudinal modes with Im[β ] > 0 are given by solid cyan(gray) lines while longitudinal modes with negative Im[β ] are in dashed cyan(gray). It is interesting to note that the longitudinal modes extend to the positive edge of the Brillouin zone and that the "fold back" only occurs when Im[β ] of the longitudinal mode is negative. One can also remark that transverse modes, T2, with both positive and negative Im[β ] exist above the light line, but that their imaginary parts are quite large. Longitudinal modes above the light line also exist at frequencies below the particle resonance maximum, but these modes remain quite close to the light line.

Chain sums

Hankel function chain sums

A periodic chain of wave scattering is defined by a lattice vector a a a, such that there is an elementary 'scatterer' at all positions r r r j i.e.: r r r j ≡ ja a a j ∈ Z j = -∞, ... -2, -1, 0, 1, 2, ..., ∞ . (6.70)

A chain sum for a quasi-periodicity vector β β β is defined:

S C n,m (k, a, β ; a a a) ≡ S C,H n,m (k, a, β ; a a a) ≡ ∑ j =0 j∈Z
H n,m ( ja a a) e i jaβ β β • a a a , (6.71)

One can remark that the chain sum, S C n,m , depends on the amplitude of the lattice vector a = |a a a|, and its direction, and the scalar product between a a a and another vector β β β which we will call the 'incident' or 'quasi-periodicity' vector. The chain sum in fact only depends on the scalar product between β β β and the periodicity vector:

β ≡ β β β • a a a .
(6.72)

One remarks that the direction of a a a depends on the orientation of the coordinate system. We can take advantage of this fact to define the z axis as the direction of a a a such that: r r r j = ja z z z , (6.73) but one must keep in mind that the expression for S C n,m is reference frame dependent. In this coordinate system, the chain sum, S C n,m (k, a, β ; z z z), takes the form :

S C n,m (k, a, β ; z z z) ≡ ∑ j =0 j∈Z H n,m ( ja z z z) e iβ β β • z z z ja = ∑ j =0 j∈Z h n (k | j| a)Y n,m j | j| z z z e iβ a j = δ m,0 λ n,0 ∞ ∑ j=1
h n ( jka) e i jβ a + (-1) n e -i jβ a , (

where we used the fact that only the m = 0 scalar spherical harmonics are non null at θ = 0, π:

Y n,m (0, 0) = δ m,0 λ n,0 = 2n + 1 4π Y n,0 (π, 0) = δ m,0 (-) n Y n,m (0, 0) . (6.75)
The analytical expressions for the first few Hankel function are:

h 0 (x) = - i x e ix h 1 (x) = e ix -1 x - i x 2 h 2 (x) = e ix i x - 3 x 2 - 3i x 3 , (6.76)
which are readily obtained from the general analytic expression for Hankel functions of arbitrary order:

h n (x) = (-i) n+1 n ∑ s=0 i s 2 s s! (n + s)! (n -s)!
e ix x s+1 . (6.77)

Integral technique for Hankel lattice sums

Generalizing the Weyl integral to produce an integral expression for spherical Hankel functions gives us the integral:

h n (kr) = 1 i n ∞ 0 dK K P n (k z /k) exp (iγ z kr) γ z , (6.78)
where we recall that

k z = √ k 2 -K 2 and K = k 2 x -k 2
y is the wavevector component in the xOy plane. In the last line we defined the dimensionless quantities, K ≡ K/k and γ z = 1 -K 2 so that this integral could be calculated using dimensionless variables. Now if we actually try to evaluate this integral, we will have numerical problems when we go past the point where k z = 0. Since the singularity coming from the k z denominator lies just above the real axis, we can analytically continue the integration into the fourth quadrant of the complex plane. Any angle will do as long as the resulting line integral is sufficiently far from the positive real axis or the negative imaginary axis. We will generally take an angle of 45 • as a reasonable compromise. We will find that the integrand will decrease exponentially for large |K| in the complex plane so that we don't have much problem with the integral extending to infinity.

Thanks to the integral expression for Hankel functions of eq.( 6.78), we are now ready to treat an infinite chain sum for a chain oriented along the z axis:

S C n,m (k, a, β ; z z z) = ∑ j∈Z * exp (iβ a j) h n (k | j| a)Y n,m j | j| z z z = ∞ ∑ j=1 exp (iβ a j) h n (ka j)Y n,m (0, 0) + ∞ ∑ j=1 exp (-iβ a j) h n (ka j)Y n,m (π, 0) = δ m,0 2n + 1 4π ∞ ∑ j=1 exp (iβ a j) h n (ka j) + (-) n ∞ ∑ j=1 exp (-iβ a j) h n (ka j) (6.79)
where we used:

Y n,m (0, 0) = δ m,0 2n + 1 4π P n (1) , Y n,m (π, 0) = δ m,0 2n + 1 4π P n (-1) , (6.80)
and P n (1) = 1 , P n (-1) = (-1) n . (6.81)

Using the integral relation of eq.( 6.78), we have:

S C n,m (k, a, β ; z z z) = δ m,0 2n + 1 4π 1 i n k ∞ 0 dK K P n (k z /k) k z × ∞ ∑ j=1 exp [i (k z + β ) jd] + (-) n ∞ ∑ j=1 exp [i (k z -β ) ja] . (6.82)
We have finally an integral expression for the chain sums:

S C n,m (k, a, β ; z z z) = ∑ j∈Z * exp (iβ jd) h n kr j Y n,m r r r j = δ m,0 2n + 1 4π 1 i n ∞ 0 dK K P n (γ z ) γ z ×   1 exp -i k z + β ka -1 + (-) n exp -i k z -β ka -1   . ( 6
.83)

Polylog approach to Hankel chain sums

Inspection of eqs.(6.79) and (6.77) shows that all terms in the chain sum can be expressed in terms of polylogarithm functions which are defined by [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]:

Li n (z) = ∞ ∑ j=1 z j j n . (6.84)
The chain sum expressed in terms of polylogarithms is then:

S C n,m (k, a, β ; z z z) = δ m,0 2n + 1 4π n ∑ s=0 (-i) n+1 i s 2 s s! (n + s)! (n -s)! × (Li s+1 exp [i (k + β ) a] + (-) n Li s+1 exp [i (k -β ) a]) (ka) s+1 .
(6.85)

This was the chain sum for outgoing Hankel functions, but we will also sometimes be interested in incoming Hankel functions, or Bessel functions of the fourth kind. These are expressed:

h (4) n (x) ≡ h - n (x) = j n (x) -iy n (x) , (6.86)
and their chain sums are:

S C n,m (k, a, β ; z z z) = δ m,0 2n + 1 4π n ∑ s=0 (-i) n+1 i s 2 s s! (n + s)! (n -s)! × (Li s+1 exp [-i (k -β ) a] + (-) n Li s+1 exp [-i (k + β ) a]) (ka) s+1 .
(6.87)

Bessel function chain sums

Although fully analytic expressions for Hankel function chain and lattice sums do not seem to exist currently, the Bessel functions lattice and chain sums do have analytic expressions. These Bessel function sums are useful in their own right for certain applications:

S C,J n (k, a, β ; z z z) = ∞ ∑ j=-∞, j =0
Y n,m r r r j j n ( jka) e i jβ a

= ∞ ∑ j=-∞ Y n,m r r r j j n ( jka) e i jβ a - ∞ ∑ j=-∞ Y 0,0 ( r r r 0 ) j 0 ( jka) = ∞ ∑ j=-∞ Y n,m r r r j j n ( jka) e i jβ a - 1 √ 4π δ n,0 . (6.88)
Using the integral expression for Y n,m r r r j j n ( jka) as an integral over the directions a wavenumber k k k as derived in eq.( 6.50) allows us to write:

j n kR j Y n,m r r r j e iβ β β •( z z z ja) = 1 4πi n Y n,m k k k e iβ β β • z z z ja e ik k k• z z z ja dΩ k k k . (6.89)
The lattice sum of the Bessel type then can be written:

S C,J n (k, a, β ; z z z) = 1 4πi n Y n,m k k k ∞ ∑ j=-∞ e i(k k k+β β β )•( z z z ja) dΩ k k k - 1 √ 4π δ n,0 . (6.90)
At this point, one invokes the Poisson summation formula which can be written formally as:

∞ ∑ r r r j j j ∈Λ e ik k k•r r r j = (2π) d Λ A ∑ p p p g ∈Λ * δ k k k -p p p g , (6.91)
where A is the "volume" of the reciprocal cell. For the 1-D sum in eq.( 6.90), this can be written:

∞ ∑ j=-∞ e i(k k k+β β β )•( z z z ja) = 2π a ∞ ∑ g=-∞ δ k z + β z - 2π a g . (6.92)
We then write this relation in a dimensionless form:

∞ ∑ j=-∞ e i(k k k+β β β )•( z z z ja) = 2π kai n ∞ ∑ g=-∞ δ k z k + β z k -g 2π ka . (6.93)
Putting this relation into the k k k integral of eq.( 6.90), we then obtain a finite sum expression for S C,J n :

S C,J n (k, a, β ; z z z) = - 1 √ 4π δ n,0 + πi n ka g max ∑ g=g min Y n0 cos β z,q , (6.94)
where since -1 < k z /k < 1we only sum over those values of g for which

-1 < ℜ β z a + 2πg ka < 1 . (6.95)
The values g min and g max are:

g min = - β z a -ka 2π + 1 g max = - β z a + ka 2π . (6.96)
The angle cos β z,g in eq.( 6.94) is given by:

cos β z,g ≡ β z a + 2πg ka = ℜ [β z /k] + iℑ [β z /k] + g 2π ka , (6.97)
where we used the parity relation:

Y n,m (-r r r) = (-1) n Y n,m (-r r r) . (6.98)
We recall that the sin and cosines for a complex angle, θ k = θ + θ , are given by: cos θ k = e iθ e -θ + e -iθ e θ 2 = cos θ cosh θi sin θ sinh θ , (6.99) and sin θ k = e iθ e -θe -iθ e θ 2i = sin θ cosh θ + i cos θ sinh θ .

(6.100)

Chain sum rotation

The chain sums expressions given in eqs.(6.83), (6.87), and (6.94) all took advantage of the facilities presented by orienting the chain of particles along the z axis. When performing lattice reduction techniques, it is necessary to have chain sums in other orientations. The chain sum is obtained by applying:

D (n) 0,m (θ , φ ) = 4π 2n + 1 Y n,m (θ , φ ) . (6.102)
Thus is trivial to write chain sums of any type (J ,H ,Y ) in an arbitrary orientation, r r r, in terms of the chain sum in the direction z z z by the simple relation:

S C n,m (β ; r r r) = S C n (β ; z z z) 4π 2n + 1
Y n,m (θ r r r , φ r r r ) . (6.103)

An orientation along the x axis is for example:

S C n,m (β ; x x x) = S C n (β ; z z z) 4π 2n + 1 Y n,m π 2 , 0 , (6.104)
is useful when carrying out a monolayer sum in the next section. An orientation along the y axis is for example:

S C n,m (β ; y y y) = S C n (β ; z z z) 4π 2n + 1 Y n,m π 2 , π 2

Integral technique

For the integral technique, it is useful to adopt a coordinate system where the a a a lattice vector lies along the y axis. In this coordinate system, the basis vectors are a a a = (0, a, 0 In this case, fixing j b = 0 and summing over j a corresponds to a chain sum along the y axis, and j b ≶ 0, corresponds to a term in the z ≶ 0 half plane respectively. The Mono-Layer (ML) lattice sum, S ML n,m , can be written:

S ML n,m = S C n,m (β 1 ; y y y) + S ML+ n,m + S ML- n,m , (6.107)
where S C n,m (β 1 ; y y y) is the chain sum along the y axis, and S ML n,m is the sum of all the sites with z = 0 The lattice sum for all j b > 0 (z > 0) sites can be expressed as an integral:

S ML+ n,m = - (-) m i n ka ∞ ∑ g=-∞ ∞ 0 dk x γ z P (m) n (γ) k g -ik x m + k g + ik x m × 1 1 -exp -i kb 2 β 2 -k g + kb 1 β 1 + γ , (6.108)
while the lattice sum for all j b < 0, can be expressed:

S ML- n,m = - 1 (-i) n ka ∞ ∑ g=-∞ ∞ 0 dk x γ z P (m) n (γ) k g -ik x m + k g + ik x m × 1 1 -exp -i kb 2 -β 2 + k g + kb 1 -β 1 + γ . (6.109)
In both of these expressions, we defined k g as the reciprocal lattice vector along the y axis:

k g ≡ β 2 + 2πg/a = kk g , (6.110)
and normalized the components of the quasi-periodic vector β as:

β 1 ≡ β 1 /k β 2 ≡ β 2 /k , (6.111)
and γ is reciprocal lattice vector along the z axis:

γ ≡ 1 -k 2 g -k 2 x . (6.112)
At the end of this calculation, one should keep in mind that lattice sum was carried out in a system where the lattice was in the yOz plane. One can obtain the expression for the lattice sum in the xOy plane by rotating the lattice sums by 90 • around the y axis in a clockwise manner, and then 90 • around the new y axis, and finally 90 • around the new z axis.

Modified Bessel function sums

Although the integral technique is rather efficient, one may prefer to obtain do a little more analytic work and obtain the lattice sum in a manner which takes the from of a lattice sum in a 2D host space. As for the integral it proves convenient to place the lattice in the yOz plane, but this time, one places the a a a lattice vector along the z axis, so that the lattice vectors can be written: The lattice sites in this system can be expressed: The reciprocal lattice is given by:

K K K g g g=(g a ,g b ) = g a a a a + g b b b b , (6.114)
where the reciprocal lattice vectors are:

a a a = 1 ab 2 (0, -b 1 , b 2 ) b b b = 0, 1 b 2 , 0 . (6.115)
This time, lattice reduction is performed by carrying out the lattice sum on the z axis of the working coordinate system which is to say that one sets j b = 0, and sums over all j a . The lattice sum is then achieved by: S ML n,m = S C n,m (β 1 ; z z z) + S ML,+ n,m + S ML,- n,m , (

where S ML,± n,m is the sum of all sites except those along the z axis:

S ML,± n,m = ∑ j b ∈Z * e i j b (β 1 b 1 +β 2 b 2 ) ∞ ∑ j a =-∞
e i j a β 1 a H n,m kr r r j j j . (6.117)

Since j b = 0 in this sum and all lattice sites are in the yOz plane, the azimuthal angle of r r r j j j is either π/2 or -π/2 for j b > 0 or j b < 0 respectively. This allows us to conclude in this plane, H n,m doesn't depend on the sign of m:

H n,-m kr r r j j j = H n,m kr r r j j j . (6.118)

We then appeal to an integral representation:

h n (kr) P m n (cos θ ) = (-i) n+1 π ∞ -∞ e ikzt K m (-ikργ (t)) P m n (t) dt (6.119) = (-i) n-m π ∞ -∞ e ikzt H m (kργ (t)) P m n (t) dt , (6.120)
where z = r cos θ , kρ = (kr) 2 -(kz) 2 > 0, and K m (z) is a modified Bessel function defined by: K m (z) ≡ i m+1 H m (iz) , (6.121) and finally γ (t) is defined such that:

γ (t) = i √ t 2 -1 |t| ≥ 1 √ 1 -t 2 t < 1 . (6.122)
The modified Bessel functions, K m (z). The j a sum can then be written:

∞ ∑ j a =-∞ e i j a β 1 a H n,m kr r r j j j = (-i) n π (-) m [sgn ( j b )] m ∞ ∑ j a =-∞ e i j a β 1 a ∞ -∞
e ik( j a a+ j b b 1 )t H m (kργ (t)) P m n (t) dt (6.123)

= (-i) n π (-) m [sgn ( j b )] m ∞ -∞ ∞ ∑ j a =-∞ e i j a (β 1 +kt)a e ik j b b 1 t H m (kργ (t)) P m n (t) dt . (6.124)
Using the 1D Poisson sum formula, we have:

∞ ∑ j a =-∞ e i j a (kt+β 1 )a = 2π a ∞ ∑ g=-∞ δ kt + β 1 + g 2π a = 2π ka ∞ ∑ g=-∞ δ t + β 1 k + g 2π ka . (6.125) ∞ ∑ j a =-∞ e i j a β 1 a H n,m kr r r j j j = (-i) n π (-) m [sgn ( j b )] m ∞ -∞ e ik j b b 1 t ∞ ∑ j a =-∞ e i j a (β 1 +kt)a H m (kργ z (t)) P m n (t) dt = 2 (-i) n ka (-) m [sgn ( j b )] m ∞ ∑ g=-∞ e -iβ 1p j b b 1 H m (kb 2 | j b | γ g ) P m n -β 1,g = 2i n ka [sgn ( j b )] m ∞ ∑ g=-∞ e -iβ 1,g j b b 1 H m (kb 2 | j b | γ g ) P m n β 1,g , (6.126)
where β 1,g and β 1,g are defined:

β 1,g ≡ β 1 + g 2π a β 1,g ≡ β 1 k + g 2π ka (6.127) and γ g ≡ γ β 1,g . (6.128)
We have therefore the monolayer sum:

S ML,± n,m = 2i n ka ∞ ∑ g=-∞ P m n β 1,g ∑ j b ∈Z * e i j b( β 1 -β 1,g) b 1 e i j b β 2 b 2 [sgn ( j b )] m H m (kb 2 | j b | γ g ) = 2i n ka ∞ ∑ g=-∞ P m n β 1,g ∞ ∑ j b =1 e i j b w g + (-) m e -i j b w g H m (kb 2 j b γ g ) = 2i n ka ∞ ∑ g=-∞ P m n β 1,g S m (w g , kb 2 γ g ) , (6.129)
where we defined:

w g ≡ β 2 b 2 - 2pπb 1 a . ( 6

.130)

One can rotate the lattice sum to put it in the desired coordinate system in the xOy plane. This can be obtained simply by a rotation of 90 • about the y axis. The expression of eq.( 6.129) is still in the form of the 2 double infinite sums like our initial expression. However, only a finite number of g values will correspond to propagating modes, i.e. β 1,g < 1, and the other values for g correspond to evanescent modes and are exponentially convergent. therefore infinite series sums. The j b sum in eq.( 6.129) is known as a Schlömilch series, and it can be expressed as a finite sum of Bernoulli polynomials.

Schlömilch series

The Schlömilch series can be expressed

S m (λ , µ) ≡ ∞ ∑ j=1 e iλ j + (-) m e -iλ j H m (µ j) . ( 6

.131)

The zero order sum is

S 0 (λ , µ) = -1 - 2i π C + log µ 4π + 2 Θ 0 + ∑ g∈Z * 2 Θ g + i π |g| , (6.132)
where C 0.5772 is Euler's constant and

Θ g = µ 2 -λ 2 g 1/2 λ g = λ + 2gπ . ( 6
Y n,m (θ , φ ) = 2n + 1 4π (n -m)! (n + m)! 1 2 P m n (cos θ ) exp(imφ ) ≡ P m n (cos θ ) exp(imφ ) , (6.134)
where in the second line we have introduced the normalized associated Legendre functions, P m n (cos θ k ) ≡ λ n,m P m n (cos θ k ), where the λ n,m normalization factor is defined:

λ n,m ≡ 2n + 1 4π (n -m)! (n + m)! 1 2
. (6.135)

These scalar spherical harmonics are normalized with respect to an integration over the solid angles :

4π 0 dΩY * ν,µ (θ , φ )Y n,m (θ , φ ) ≡ (-1) µ π 0 sin θ dθ 2π 0 dφ Y ν,-µ (θ , φ )Y n,m (θ , φ ) = δ n,ν δ m,µ . (6.136)
In principal, the Legendre polynomials, P n (x) = P 0 n (x), can be obtained from Rodrigues' formula:

P n (x) = 1 2 n n! d n dx n x 2 -1 n , (6.137)
but in practice we will calculate them with recurrence relations. Likewise, the associated Legendre functions could be obtained for from the expression:

P m n (x) = (-1) m 1 -x 2 m/2 d m dx m P n (x) . (6.138)
Their calculation is simplified by noting that the normalized associated Legendre functions have the convenient parity property that:

P -m n (x) = (-1) m P m n (x) . (6.139)
There are alternative ways of calculating the scalar spherical harmonics that are better for formulating lattice sums and reflections from a physical interface. In lattice sums and reflections from surfaces, the spherical harmonics will be evaluated in terms of the direction of the incident or reflected wavevectors, k k k:

Y n,m (θ k , φ k ) = Y n,m k k k = Y n,m (k k k/ / /k) = Y n,m (k x /k, k y /k, k z /k) , (6.140)
where we recall that:

k z k = cos θ k k x /k = sin θ k cos φ k k y /k = sin θ k sin φ k , (6.141)
and we keep in mind that P m n are functions of cos θ k = k z /k. Since x = cos θ , and the P n (x) are polynomials in x, the d m dx m P n (x) are functions of cos θ . The factor 1x 2 m/2 corresponds to sin m θ k with no ambiguity in sign since Re {θ k } is ∈ (0, π). One should remark that the 1x 2 m/2 is non-polynomial so that is why one refers to them as associated Legendre functions. For applications involving reciprocal space and/or integrations in the complex plane it proves useful to explicitly extract this factor, and define associated Legendre polynomials, which we shall denote, P m n (not to be confused with the normalized associated Legendre functions).

For positive m we have then:

Y n,m (k x /k, k y /k, k z /k) = P m n (cos θ k ) exp (imφ k ) = λ n,m (-1) m sin m θ k (cos φ k + i sin φ k ) m d m dx m P n ( k z k ) = (sin θ k cos φ k + i sin θ k sin φ k ) m (-1) m λ n,m d m dx m P n ( k z k ) = (k x /k + ik y /k) m (-1) m λ n,m d m dx m P n ( k z k ) = (k x /k + ik y /k) m P m n ( k z k ) = K k |m| exp (imφ k ) P m n ( k z k ) . (6.142)
where we have defined the normalized associated Legendre polynomials, P m n , such that:

P m n ( k z k ) ≡ (-1) m λ n,m d m dx m P n ( k z k ) . (6.143)
The parameter, K K K ≡ k x x x x + k y y y y , (6.144) corresponds to the momentum space vector in the x-y plane.

The wave vector components k x , k y , and k z are related to the possibly complex angles, θ k and φ k , via the relations:

k x = k sin θ k cos φ k k y = k sin θ k sin φ k k 2 z = k 2 -K 2 = k 2 cos 2 θ k , (6.145)
This relation of eq.( 6.142) for Y n,m can be extended to negative m by writing :

Y n,m (k x /k, k y /k, k z /k) = K k |m| exp (imφ k ) P m n ( k z k ) m ≷ 0 , (6.146)
as long as we define P -m n such that :

P -m n ≡ (-1) m P m n . (6.147)
The objective of the above procedure was to define P m n (x) that are always polynomials of x for both positive and negative m. This is in contrast to the associated Legendre functions P m n (x) which are not polynomials in terms of x.

Translation-addition theorem for scalar partial waves

Let us consider a point M in a system using spherical coordinates. We consider a second system of spherical coordinates centered on the position r r r 0 . The position of M in this second system centered on r r r 0 is: r r r = r r r --r r r 0 . (6.148)

We take the usual convention of outgoing scalar partial waves as products of spherical Hankel function and scalar spherical harmonics:

Ψ H ,n,m (kr r r) ≡ Ψ (3) n,m (kr r r) ≡ h n (kr)Y n,m (θ , φ ) , (6.149)
while the regular scalar partial waves replace the spherical Hankel functions with spherical Bessel functions: Ψ J ,n,m (kr r r) ≡ Ψ

n,m (kr r r) ≡ j n (kr)Y n,m (θ , φ ) . (6.150)

One can construct a row 'matrices' composed of the Ψ J ,n,m (kr r r) or Ψ H ,n,m (kr r r) functions respectively then the translation-addition theorem for scalar partial waves can be compactly expressed in matrix form:

Ψ t H (kr r r) = Ψ t H kr r r . α(kr r r 0 ) r > r 0 Ψ t H (kr r r) = Ψ t J kr r r . β (kr r r 0 ) r < r 0 Ψ t J (r r r) = Ψ t J r r r . β (kr r r 0 ) ∀ r r r j , (6.151)

where the elements of the irregular translation-addition matrix, α have extremely simple expressions in terms of the 3Y coefficients : stands for the Wigner 3J coefficients. It is worth remarking that the '3Y' coefficients of eq.(6.154) are closely related to the Gaunt coefficients developed in quantum mechanics for treating the helium atom (mostly differing on account of different normalization conditions).

Vector translation-addition theorem

The vector translation-addition theorem are the vector analogue of the scalar translation theorem discussed above in section6.8.2. This would be an almost trivial extension of the scalar addition theorem if we were working with solutions of the vector Helmholtz equation in a Cartesian basis like those of eqs.(6.7) and (6.8). The additional complication is due to the fact we want to the vector translation-addition theorem to act on the purely transverse waves like those of eq.(6.9). Defining column matrices, Ψ Ψ Ψ , composed of the transverse vector partial waves, the vector translation-addition theorem is written:

Ψ Ψ Ψ t H (kr r r) = Ψ Ψ Ψ t H kr r r J(kr r r 0 ) r > r 0 Ψ Ψ Ψ t H (kr r r) = Ψ Ψ Ψ t J kr r r H(kr r r 0 ) r < r 0 Ψ Ψ Ψ t J (r r r) = Ψ Ψ Ψ t J r r r J(kr r r 0 ) ∀ |r 0 | , (6.156)

where the matrix J(kr r r 0 ) matrix can be expressed in terms of spherical scalar β (kr r r 0 ) matrices of eq.(6.153) (expressed in terms of spherical Hankel functions) while the H(kr r r 0 ) matrices can be expressed in terms of the α(kr r r 0 ) matrices of eq.(6.152).

Explicitly, the H(kr r r 0 ) matrix can be expressed:

H(kr r r 0 ) = A ν,µ;n,m (kr 0 , θ 0 , φ 0 ) B ν,µ;n,m (kr 0 , θ 0 , φ 0 ) B ν,µ;n,m (kr 0 , θ 0 , φ 0 ) A ν,µ;n,m (kr 0 , θ 0 , φ 0 ) . (6.157)

The vector coefficients A ν,µ;n,m are then calculated using : (cos β ) . (6.161) 6.9 Recurrence relations for special functions Partial wave descriptions are composed of products of spherical harmonic and spherical Bessel types special functions. For a numerical analysis, it is important to calculate these functions rapidly and accurately. Recurrence relations prove to be a good manner to achieve this goal. Multipole expansions must be truncated to a given order n max , which determines the the strength of spatial field variations. Inspection of the translation-addition theorem formulas show us that we will need to evaluate spherical harmonic and spherical Hankel functions up to order 2n max .

starting with the initial value χ 0 (z) =cos z.

Having calculated ψ n (z) via eq.( 6.178) and χ n (z) via eq.( 6.186), one can finally construct h n (z) by h n (z) ≡ j n (z) + iy n (z) , (6.187)

with the real and imaginary parts of ξ n (z) now both being calculated up to machine accuracy.

The ratio of spherical Bessel functions to spherical Hankel functions also occurs frequently in Mie theory, and I found it convenient and more accurate to calculate these ratios directly using upward recurrence relations, notably:

j n (z) h n (z) = j n-1 (z) h n-1 (z) n -ϕ n-1 (z) n -ϕ (3)
n-1 (z) For coated spheres, it is can also useful to use the analogous recurrence relation:

, ( 6
j n (z) y n (z) = ψ n-1 (z) y n-1 (z) n -ϕ n-1 (z) n -ϕ (2)
n-1 (z)

, ( 6

Figure 6 . 1 :

 61 Figure 6.1: Particulate grating with lattice vectors a a a and b b b.

  a a a = (a, 0, 0), b b b = (b x , b y , 0), and c c c = (c x , c y , c z ).

  andT β β β a = ta + tΩ β β β T β β β a . (6.40)

Figure 6 . 2 :

 62 Figure 6.2: Real and imaginary parts of the dispersion relations in the dipole approximation (dashed curves), and fully converged multipole calculations with n max = 10 (full lines). The Longitudinal mode with positive imaginary part is in cyan(gray) the "T1" mode with positive imaginary part is in blue(black line). The

  Figures 6.2a) and 6.2b) merit some commentary. It is immediately clear that the dipole approximation provides a moderately accurate prediction of dispersion relations only over a narrow range of frequencies for which the imaginary part of the propagating wavevector is

Figure 6 . 3 :

 63 Figure 6.3: Normalized extinction is a solid blue (line) and scattering cross section given by a dashed green line of a silver monomer in terms of frequency (a = 25nm). reproduced with permission : http://dx.doi.org/10.1364/JOSAB.29.001012

Figure 6 . 4 :

 64 Figure 6.4: Positive and imaginary parts of the dispersion relations in the Re[β ] > 0 part of the Brillouin zone. Transverse modes with: Im[β ] > 0 modes are solid blue(black) lines, while that with Im[β ] < 0 is given by a dashed blue(black) line. Longitudinal modes with Im[β ] > 0 are solid cyan(gray) lines, while those with Im[β ] < 0 is a dashed cyan(gray) line. reproduced with permission : http://dx.doi.org/10.1364/JOSAB.29.001012

  Y n,m ( r r r) = Y n,0 ( z z z) D (n) 0,m (θ r r r , φ r r r ) = (θ r r r , φ r r r ) , (6.101) which just translates the relation derived in Edmonds (eq.(4.1.25) page 59) that:

. ( 6

 6 .105)6.7 2D Grating lattice sumsA 2D periodic media is characterized by two basic lattice vectors a a a, b b b. Although, we want to describe a system with lattice vectors a a a = (a, 0, 0), and b b b = (b 1 , b 2 , 0), we are going to work in a rotated coordinate systems in which the lattice will be placed in the xOy plane.

  ), b b b = (0, b 2 , b 1 ), then lattice sites are given by: r r r j j j=( j a , j b ) = j a a a a + j b b b b = (0, j a a + j b b 2 , j b b 1 ) .(6.106)

  a a a = (0, 0, a) , and b b b = (0, b 2 , b 1 ) .

r r r j

  j j=( j a , j b ) = j a a a a + j b b b b = (0, j b b 2 , j a a + j b b 1 ) = j b b 2 y y y + ( j a a + j b b 1 ) z z z . (6.113)

  α ν,µ,nm (kr r r 0 ) = 4πi ν-n n+ν ∑ p=|n-ν| i p 3Y (n, m; ν, µ; p) h p (kr 0 )Y p,m-µ (θ 0 , φ 0 ) ,(6.152) while the the elements of the regular translation-addition matrix, β ν mu,;n,m , are the same as the α ν,µ;n,m coefficients but with the j p (kr 0 ) replacing the h p (kr 0 ) function i.e.:β ν,µ;n,m (kr r r 0 ) ≡ 4πi ν-n n+ν ∑ p=|n-ν| i p 3Y (n, m; ν, µ; p) j p (kr 0 )Y p,m-µ (θ 0 , φ 0 ) ,(6.153)where 3Y (n, m; ν, µ; p) are the 3Y coefficients defined by the angular integral of three scalar spherical harmonics:3Y (n, m; ν, µ; p) ≡ m (θ , φ )Y * ν,µ (θ , φ )Y * p,m-µ (θ , φ ) sin θ dθ dφ = (-) µ (-) m-µ π 0 2π 0 Y n,m (θ , φ )Y ν-µ (θ,φ )Y p,µ-m (θ , φ ) sin θ dθ dφ = (-) m (2n + 1)(2ν + 1)(2p + 1

A ν,µ;n,m = 1 2 1 ν 2 2ν + 1 2ν - 1 1ν

 121 (ν + 1) n (n + 1) 2µmα ν,µ;n,m + + (nm) (n + m + 1) (νµ) (ν + µ + 1)α ν,µ+1;n,m+1+ (n + m) (nm + 1) (ν + µ) (νµ + 1)α ν,µ-1;n,m-1 . (6.158)When filling up a matrix, with the A ν µ,nm coefficient, we should fill them up column by column (calculate all the ν, µ) elements for a fixed n, m. Then, for each n, m, the B ν µ,nm coefficients can be calculated from the previous (i.e. ν -1) scalar coefficients :B ν,µ,n,m = -i 1 (ν + 1) n (n + 1) 2m (νµ) (ν + µ)α ν-1,µ;n,m + (nm) (n + m + 1) (νµ) (νµ -1)α ν-1,µ+1;n,m+1 -(n + m) (nm + 1) (ν + µ) (ν + µ -1)α ν-1,µ-1;n,m-1 .(6.159)6.8.4 Rotation matricesUnder rotation, each of the four blocks of a vector translation-addition matrix transform following the rotation matrix, D (α, β , γ), which is expressed in terms of the 3 Euler angles, α, β , and γ. The D (α, β , γ) matrix elements are described in detail in ref.[START_REF] Edmonds | Angular Momentum in Quantum Mechanics[END_REF], and are block diagonal in the orbital (multipole) 'quantum' number, n :[D (α, β , γ)] ν,µ,n,m = δ n,ν exp (iµα) d (n) µ m (β ) exp (imγ) .(6.160)The elements d(n)µm are standard,[START_REF] Edmonds | Angular Momentum in Quantum Mechanics[END_REF] and the d(n)µ,m term in the rotation matrices can be expressed in terms of the Jacobi polynomials[START_REF] Edmonds | Angular Momentum in Quantum Mechanics[END_REF] :d (n) µm (β ) = (n + µ)! (nµ)! (n + m)! (nm)

1 -

 1 ize -2izize -2iz i + z

  zsin z z sin z + cos z . (6.194) Other useful relations are obtained from the classic spherical Bessel function recurrence relations:f n (z) = z f n-1 (z) + z f n+1 (z) 2n + 1 f n (z) = n f n-1 (z) -(n + 1) f n+1 (z) 2n + 1 ,(6.195)

Recurrence relations for associated Legendre polynomials

The recurrence relations for the P m n polynomials are are initialized with:

We can then calculate all the P n n up to 2n max with the recurrence relation:

.163)

The P m n with m = n -1 are calculated via:

All the remaining P m n with m = 1, ..., n -2 can be calculated for each n = 3, ..., 2n max using the relation :

Although we obtain the P 0 n in the above scheme, it can sometimes prove useful to obtain the normalized Legendre Polynomials through the recurrence relation :

2n + 1 2n -3 P 0 n-2 (u) . (6.166)

The P m n with negative values of m are calculated using :

.167)

Logarithmic Bessel functions

The spherical Bessel, Neumann and Hankel functions of the Ricatti form, are simply these functions multiplied by their argument. The advantage of this form is that they have better limit properties for small arguments. Their definitions are respectively:

Multiplying the logarithmic derivatives of these functions by their argument defines the functions, ϕ (1) , ϕ (2) , ϕ (3) :

The ϕ

n can also be generated by particularly efficient numerical recursion relations. They also provide particularly symmetric expression for the Mie coefficients of spherical scatterers and formulations of matrix balancing. The Wronskian relation for Ricatti-Bessel functions:

takes a simple form in terms of the ϕ

(1,2,3) n functions as:

. (6.171)

The ϕ (3) can be reliably calculated numerically from the upward recurrence relation:

starting with an initialization of

The ξ n (z) functions in most situations can then be readily calculated numerically from the recurrence relation:

n-1 (z) , (6.174) starting from the initial value ξ 0 (z) = -ie iz . One should note that analytical expressions exist for the spherical Ricatti-Hankel functions, ξ n (z), and these can be useful at low multipole order:

When one deals with high multipole orders however, it usually is more practical to exploit the recurrence relations of eq.( 6.172) and (6.174).

The regular Ricatti Bessel functions, ϕ

n (z), obey the same recurrence relations as the ϕ

(3) n (z) functions. If one calculates them by via upward recurrence, things may work fine for the first few recurrence calculations, but at some point, the recurrence relations can go completely off course. The solution to this problem has been known for quite some time is that the ϕ n (z) functions should be calculated starting from high values of n in a reverse recurrence relation. I start usually with n equal to at least n max + 20 where n max is the largest value that I want to use in calculations, and start simply with ϕ n max +20 (z) = 0. The ϕ n (z) functions so obtained have always been the correct ones up to machine precision. The reverse recurrence relation is:

.176)

One can check calculations by verifying that the ϕ 0 (z) obtained by backward recurrence is equal to the analytical result:

Once the ϕ n (z) functions have been calculated, one can readily generate the ψ n (z) functions with the upward recurrence relation which is the direct analogue of eq.( 6.174)

starting with the initial value ψ 0 (z) = sin z. Analytical expressions for the lowest ψ n (z) are:

It is perhaps worth remarking that, there is one potential numerical problem with using eq.( 6.174) to calculate spherical Hankel functions. For some values of z the real and imaginary parts of the spherical Hankel functions can have extremely different absolute values. For concreteness, let us assume that |Im (h n (z))| |Re (h n (z))|, then using eq.( 6.174), the calculated value of Im (h n (z)) will usually be quite inaccurate if its absolute value is less than last significant figure in the calculation of Re (h n (z)). This problem can be circumvented (for real values of z at least) by calculating the spherical Neumann functions (denoted here by y n (z) but some authors denote it n n (z)). We recall that the Neumann functions are real-valued provided that z is real valued.

The Ricatti Neumann functions are defined :

The first few Ricatti Neumann functions, χ n (z), are

We define a ϕ (2) 'logarithmic derivative' Neumann function as :

.183)

We can calculate the ϕ

n from the upward recurrence relation:

with an initialization of ϕ

(2)

Once the ϕ

n functions have been calculated, one can readily generate the χ n (z) functions with the upward recurrence relation which is the direct analogue of eq.( 6.174)

n-1 (z) , (

with f n (z) = j n (z) , h n (z) to obtain :

with ψ n (z) ≡ z j n (z) + j n (z) (6.197) we obtain a convenient expression for the derivative of Ricatti-Bessel functions: .198) or the expression:

(6.199)

Vector Spherical Harmonics

There is no universally accepted notation for the Vector Spherical Harmonics (VSHs). Our notation for their normalized forms is X X X n,m , Y Y Y n,m , and Z Z Z n,m where they are respectively defined by

The scalar spherical harmonics, Y n,m (θ , φ ), do have a nearly universal convention for their definitions [START_REF] Jackson | Classical Electrodynamics : Third Edition[END_REF] which we recalled in eq.(6.134).

For numerical calculations of the VSHs, it is convenient to introduce the normalized functions u m n and s m n defined as:

where the P m n are the Legendre functions defined in eqs.(6.137), (6.138), and (6.139), and γ n,m a normalization factor given by

.203)

The transverse VSHs have compact expressions in terms of u m n and s m n : The normalized u m n functions can be readily calculated with recurrence relations: