# Diffraction of Light by a Two-Dimensional Lattice of Spheres

Authors: Bernard de Dormale; Vo-Van Truong
Journal: International Journal of Optics
 (2012)

---

## Abstract

Two-dimensional arrays of particles are of great interest because of their very characteristic optical properties and numerous potential applications. Although a variety of theoretical approaches are available for the description of their properties, methods that are accurate and convenient for computational procedures are always sought. In this work, a new technique to study the diffraction of a monochromatic electromagnetic field by a two-dimensional lattice of spheres is presented. The method, based on Fourier series, can take into account an arbitrary number of terms in the multipole expansion of the field scattered by each sphere. This method has the advantage of leading to simple formulas that can be readily programmed and used as a powerful tool for nanostructure characterization.


---

## Body

## 1. Introduction

The selective absorption in the visible spectral region observed in granular noble-metal films or nanoparticles has been the subject of numerous investigations in the past [1]. This phenomenon is attributed to plasma resonance in the particles forming the granular films, distinguishing the particle optical properties from the bulk optical properties of the metal itself. Early applications of these granular films included solar absorbers [2] and selective optical filters [3]. With recent advances in particle synthesis and the capacity to form periodic nanoparticle arrays by various methods, the promise of application of these nanoparticles has increased manyfold, covering areas of sensing, optical, electrooptical, as well as solar cells and medical diagnostics [4].To characterize the two-dimensional distribution of nanoparticles on a substrate, it was frequent to use an effective-medium theory that approximated the granular films to an effective film having an effective optical thickness and effective optical properties [5–7]. The particles were considered very small compared to the light wavelength so that the dipole approximation could be used. This approach has served well the characterization of films in many studies, giving a good correlation between theory and experiments [7, 8]. When particles get slightly larger, the dipole approximation may no longer hold and retardation effects have to be taken into account, especially when nanoparticle-substrate interaction is also considered [9]. A few attempts at a comprehensive treatment are faced with conceptual difficulties and prohibitive computing times [10]. Other methods for modeling the optical properties of nanoparticle arrays include the discrete dipole approximation (DDA) [11], the T-matrix method [12], and the finite-difference time-domain (FDTD) method [13]. Each method has its own advantages and shortcomings but generally they required elaborate numerical computational procedures. Among other things, our method, based on Fourier series, can take into account an arbitrary number of terms in the multipole expansion of the field scattered by each spherical nanoparticle. It has the advantage of leading to simple formulas that can be readily programmed and used as a powerful tool for nanostructure characterization.In a previous work, we introduced a new way to solve the problem of particle arrays consistently and efficiently, when only the electric dipole radiation of the particles was considered [14, 15]. The theoretical treatment of the scattered field was based on a variant of the electromagnetic Green function method for a two-dimensional lattice of spheres. We will now show how this approach can be extended to take into account an arbitrary number of multipole terms. The present paper will deal with the principles of the method and general formulas, whenever possible, will be developed. In a separate recent publication, examples of numerical calculations have been given, showing interesting features both for dielectric nanoparticle arrays as well as metallic nanoparticle arrays, together with comparison with some experimental data [16].In the following, we will make use of the basis of vector harmonics [17, 18]:(1)Xlm-1(θ,ϕ)=l2l+1nYlm(θ,ϕ)-il+12l+1n×Xlm(θ,ϕ),Xlm0(θ,ϕ)=Xlm(θ,ϕ),Xlm1(θ,ϕ)=-l+12l+1nYlm(θ,ϕ)-il2l+1n×Xlm(θ,ϕ),
where n=x/r, Ylm(θ,ϕ) are the usual spherical harmonics, and
(2)Xlm(θ,ϕ)=1l(l+1)LYlm(θ,ϕ),L=1ix×∇
are the commonly known vector harmonics [19]. The more important properties of the Xlms(θ,ϕ)'s are listed in the appendix of [18]. We will also use the related vector-valued functions [18]:
(3)Hlms(x)=rl+sXlms(θ,ϕ),
whose components are either zero or homogeneous polynomials of degree l+s in x, y and z, which means that Hlms(λx)=λl+sHlms(x). Hlm-1(x)andHlm0(x) will appear in many formulas; the most convenient way we know to compute them is to first use [18] the following:
(4)Hlm-1(x)=1l(2l+1)∇(rlYlm(θ,ϕ)),
then
(5)Hlm0(x)=-i2l+1l+1x×Hlm-1(x).
Calculations can be made in Cartesian coordinates, since the quantities rlYlm(θ,ϕ) appearing in (4) are polynomials in x, y, and z.


## 2. Outline of the Method

The problem at hand is the diffraction of a monochromatic electromagnetic field of angular frequencyω and wave vector k by a structure made of identical, nonoverlapping spheres of radius ρ, whose centers xi=(xi,yi,0) form a hexagonal lattice in the xy-plane. Assuming for now that no substrate is present, the traditional approach to this subject can be summarized in the following four conceptual steps [20] (we will use Gaussian units throughout this paper, and a harmonic time dependence through a factor  exp(-iωt) is understood for all time dependent quantities).Step 1.
We will callEi,Bi the field radiated by the sphere centered at xi, and by Ei′,Bi′ the total electromagnetic field incident on that sphere, which is the superposition of the incoming field Ein,Bin and of the fields scattered by all the other spheres:
(6)Ei′=Ein+∑l≠iEl,Bi′=Bin+∑l≠iBl.
Because of the symmetry of the problem, all the Ei's and Bi's are merely copies of each other, differing only by a phase factor, and thus, assuming that x0 is the origin:
(7)Ei(x)=exp⁡(ik⋅xi)E0(x-xi),Bi(x)=exp⁡(ik⋅xi)B0(x-xi).
Therefore,
(8)E0′(x)=Ein(x)+∑i≠0exp⁡(ik⋅xi)E0(x-xi),B0′(x)=Bin(x)+∑i≠0exp⁡(ik⋅xi)B0(x-xi).Step 2.
E
0
′
,
B
0
′ and E0,B0appearing in (8), have the respective multipole expansions:
(9)E0′(x)=∑l,m[ikaE′(l,m)∇×(jl(kr)Xlm(θ,ϕ))E′0+aM′(l,m)jl(kr)Xlm(θ,ϕ)ik],B0′(x)=∑l,m[ikaE′(l,m)jl(kr)Xlm(θ,ϕ)E′0-ikaM′(l,m)∇×(jl(kr)Xlm(θ,ϕ))],
where jl(r) are the spherical Bessel functions of the first kind, and
(10)E0(x)=∑l,m[ikaE(l,m)∇×(hl(1)(kr)Xlm(θ,ϕ))E′0+aM(l,m)hl(1)(kr)Xlm(θ,ϕ)ik],B0(x)=∑l,m[aE(l,m)hl(1)(kr)Xlm(θ,ϕ)E′0-ikaM(l,m)∇×(hl(1)(kr)Xlm(θ,ϕ))],
where hl(1)(r) are the first Hankel functions [19]. The multipole coefficients aT′(l,m), aT(l,m) (with T=E for electric or M for magnetic) in (9) and (10) must be related by Mie's formulas:
(11)aE(l,m)=αE(l)aE′(l,m),aM(l,m)=αM(l)aM′(l,m),
where for a given value of k=|k|, the polarisabilities αT(l) are known functions of the common size and dielectric constant of the spheres [21].Step 3.
In view of (8), it should be possible to use techniques like addition theorems to express the aT′(l,m)'s in terms of the aT(l,m)'s [17]. Therefore, if we call aTin(l,m) the multipole coefficients of the incoming field, (11) could be put into a form the following:
(12)aE(l,m)=αE(l)(aEin(l,m)+∑l′m′A(l,m;l′,m′)aE(l′,m′)+∑l′m′B(l,m;l′,m′)aM(l′,m′)),aM(l,m)=αM(l)(aMin(l,m)+∑l′m′C(l,m;l′,m′)aE(l′,m′)+∑l′m′D(l,m;l′,m′)aM(l′,m′)),
from which, theoretically, the aT(l,m)'s can be computed, hence E0andB0. In practice, however, one has to reduce (12) to an equation involving only a finite number of unknowns, which is equivalent to approximating E0,B0 by a field with only a finite number of nonvanishing multipole coefficients.Step 4.
Finally, we have for the scattered field the following:(13)Esc(x)=∑iEi(x)=∑iexp⁡(ik⋅xi)E0(x-xi),Bsc(x)=∑iBi(x)=∑iexp⁡(ik⋅xi)B0(x-xi).
Again, addition theorems must be used to express the fields in (13) in terms of the aT(l,m)'s.There is a fundamental difficulty involved in the approach just described: although the series in (13) or (8) may very well converge, this seems to be no longer true when the approximation at the end of the third step is made [20]. Even though this problem can be dealt with by some averaging process [20], the resulting calculations are not only exceedingly complex, they are also far too inefficient to be of practical use. There is thus a real need for improvements, both in terms of consistency and efficiency, and we propose here a new calculation technique addressing both concerns.Our approach distinguishes itself from the traditional scheme in two fashions: first, we will use Fourier series, instead of (13), to express EscandBsc. Then, since we can write
(14)∑i≠0Ei=∑iEi-E0,∑i≠0Bi=∑iBi-B0,
we will use, instead of (8), the equivalent “subtractive” formula:
(15)E0′=Ein+Esc-E0,B0′=Bin+Bsc-B0,
where, again, EscandBsc are expressed in terms of Fourier series. To avoid dealing with singularities in the operation, the current densities of the multipoles, at this stage, are replaced by narrow-shaped Gaussian distributions. For the dipole approximation [14, 15], this method not only proved to be fully consistent, it also led to accurate results hundreds of times faster than with the usual procedure.The calculation of the Fourier series ofEscandBsc in terms of the multipole coefficients aT(l,m) of the field E0,B0 is presented in Section 3, while in Section 4, we show how to retrieve the multipole coefficients aT′(l,m) of E0′,B0′ from the expressions for Esc,Bsc,  and E0,B0. Some auxiliary results are relegated to the Appendix.


## 3. Calculation of the Scattered Field

Let us callJi the current density induced in the ith sphere (i.e., the source of Ei,Bi): we will compute the  scattered field as the field radiated by the current density:
(16)J=∑iJi.
Just as in (7), we must have the following:
(17)Ji(x)=exp⁡(ik⋅xi)J0(x-xi),
where k=(kx,ky,kz) is the incoming wave vector. Thus, (16) can be rewritten as
(18)J(x)=∑iexp⁡(ik⋅xi)J0(x-xi),
which has the following form: exp(ik·x) times a periodic function. This means that J can be developed in Fourier series [14]:
(19)J(x)=∑gJg(z)exp⁡(i(kx+gx)x+i(ky+gy)y),
where the summation is extended to all vectors g=(gx,gy) of the two-dimensional reciprocal lattice. The Fourier coefficients are given by
(20)Jg(z)=1S∬exp⁡(-i(kx+gx)x-i(ky+gy)y)J(x)dxdy,
where the integral extends to the “primitive cell” of area S in the xy-plane (note that in (20) we could have used J0 instead of J, because  J(x)=J0(x) in the primitive cell (to avoid unnecessarily technical arguments, we will assume that spheres do not extend beyond their respective cell boundaries).Our task is to compute the scattered field in terms of the multipole moments ofJ0. To achieve that, we first solve Maxwell's equations:
(21)∇×Bsc(x)+ikEsc(x)=4πcJ(x),∇×Esc(x)-ikBsc(x)=0,
using a vector potential Asc, solution of
(22)(Δ+k2)Asc(x)=-4πcJ(x),
in terms of which the field is given by [14]
(23)Esc(x)=ik(∇(∇⋅Asc(x))+k2Asc(x)),Bsc(x)=∇×Asc(x).Next, we develop the vector potential and the scattered field in Fourier series, with respective Fourier coefficientsAgsc,Egsc,andBgsc, using formulas similar to (19) and (20). Substitution of the series into (22) and (23) gives the following:
(24)(d2dz2-(kx+gx)2-(ky+gy)2+k2)Agsc(z)=-4πcJg(z),(25)Egsc(z)=ik(Dg(Dg⋅Agsc(z))+k2Agsc(z)),Bgsc(z)=Dg×Agsc(z),
where the vector operator Dg is defined by the following:
(26)Dgf(z)=(i(kx+gx)f(z),i(ky+gy)f(z),df(z)dz).
If we also define the following:
(27)νg=k2-(kx+gx)2-(ky+gy)2(0≤arg(νg)<π),kg±=(kx+gx,ky+gy,±νg),
we can write the retarded solution of (24) as
(28)Agsc(z)=2πicνg∫-∞+∞exp⁡(iνg|z-z′|)Jg(z′)dz′,
which gives the scattered field:
(29)Egsc(z)=-2πωνg{∫-∞zexp⁡(iνg(z-z′))-2πωνg×[k2Jg(z′)-(kg+⋅Jg(z′))kg+]dz′-2πωνg+∫z+∞exp⁡(iνg(z′-z))-2πωνg×[k2Jg(z′)-(kg-⋅Jg(z′))kg-]dz′∫}-4πiω(z^⋅Jg(z))z^,Bgsc(z)=-2πcνg{2πcνg∫-∞zexp⁡(iνg(z-z′))(kg+×Jg(z′))dz′-2πcνg=+∫z+∞exp⁡(iνg(z′-z))(kg-×Jg(z′))dz′}.The last step consists in expressing (29) in terms of the multipole coefficients aT(l,m) of E0,B0. Any current density with the aT(l,m)'s as multipole moments will, outside the source, radiate the same field as E0,B0. Thus, if we are only interested in the value of Esc,Bsc outside the layer of spheres, we can use a point source with the same multipole moments [18]:
(30)J0(x)=∑l=1∞c(-1)l+1kl+12l+1l+1×∑m=-ll(aE(l,m)Hlm-1(∇)-ikaM(l,m)∇×Hlm-1(∇))δ(x),
the differential operators Hlm-1(∇) in (30) are obtained by replacing x,y,andz in Hlm-1(x) by the corresponding partial derivative. This leads to the Fourier coefficients:
(31)Jg(z)=∑l=1∞c(-1)l+1Skl+12l+1l+1×∑m=-ll(aE(l,m)Hlm-1(Dg)-ikaM(l,m)Dg×Hlm-1(Dg))×δ(z),
and, after some manipulations, to the following:
(32)Egsc(z)=-exp⁡(iνg|z|)2πSνg∑l=1∞(-i)l-1kl+12l+1l+1∑m=-llkg±×(aM(l,m)Hlm-1(kg±)-1kaE(l,m)kg±×Hlm-1(kg±)),Bgsc(z)=-exp⁡(iνg|z|)2πSνg∑l=1∞(-i)l-1kl+12l+1l+1∑m=-llkg±×(aE(l,m)Hlm-1(kg±)+1kaM(l,m)kg±×Hlm-1(kg±)),
where kg+ should be used for z>0 and kg- for z<0. Using (5) and (32) can be more compactly rewritten as
(33)Egsc(z)=exp⁡(iνg|z|)2πSνg∑l=1∞(-i)lkl+1×∑m=-ll(1kaM(l,m)Hlm0(kg±)-1kaE(l,m)kg±×Hlm0(kg±)1k),Bgsc(z)=exp⁡(iνg|z|)2πSνg∑l=1∞(-i)lkl+1×∑m=-ll(1kaE(l,m)Hlm0(kg±)+1kaM(l,m)kg±×Hlm0(kg±)).


## 4. Subtraction Process

Using a point source, like we just did, is not suitable to computeE0′,B0′ through formula (15): this gives formulas for Esc,Bsc and E0,B0 that are both singular at the origin. To avoid such singularities, we will use a variant of Ewald's method [22], and replace the δ-function in (30) by a suitably normalized Gaussian distribution:
(34)Φ(x)=1(2πσ2)3/2exp⁡(k2σ22-r22σ2),
where σ is assumed to be small, compared with the distance between neighboring dipoles, so that we can presume Ji(0)=0fori≠0. The current density (30) and its Fourier coefficients (31) will thus be replaced by the following:
(35)J0(x)=∑l=1∞c(-1)l+1kl+12l+1l+1×∑m=-ll(ikaE(l,m)Hlm-1(∇)-ikaM(l,m)∇×Hlm-1(∇))Φ(x),(36)Jg(z)=∑l=1∞c(-1)l+1kl+12l+1l+1×∑m=-ll(ikaE(l,m)Hlm-1(Dg)-ikaM(l,m)Dg×Hlm-1(Dg))Φg(z),
with
(37)Φg(z)=1S(2πσ2)1/2exp⁡(νg2σ22-z22σ2).
When these coefficients are substituted in (29), the corresponding formulas for Egsc,Bgsc are more complicated than (33), because it is now necessary to use integrations by parts to compute the scattered field and it is difficult to express the results in a general manner.We will now show how to compute the multipole coefficientsaT′(l,m) of  E0′,B0′ using formula (15): the technique is based on results of [18], which show that if a vector field assumes the following form:
(38)F(x)=∑l,m,sflms(r)Xlms(θ,ϕ),
one has
(39)[(Hlms(∇))*⋅F(x)]x=0=1(2π)3/2lim⁡η→01η2(l+s)+3×∫0∞exp⁡(-r22η2)flms(r)rl+s+2dr.
The change of variable r=η2u in formula (39) allows to rewrite it as
(40)[(Hlms(∇))*⋅F(x)]x=0=2l+s-1π3/2∫0∞(lim⁡η→0flms(η2u)(η2u)l+s)exp⁡(-u)ul+s+1/2du,
that is,
(41)[(Hlms(∇))*⋅F(x)]x=0=2l+s-1Γ(l+s+3/2)π3/2×lim⁡r→0flms(r)rl+s,
where Γ(x) is the gamma function.Assuming thatσ is small enough, E0′,B0′ will assume the general form for a free electromagnetic field [18]:
(42)E0′(x)=∑l,m(-aE′(l,m)l+12l+1jl-1(kr)Xlm-1(θ,ϕ)+aM′(l,m)jl(kr)Xlm0(θ,ϕ)+aE′(l,m)l2l+1jl+1(kr)Xlm1(θ,ϕ)),B0′(x)=∑l,m(aM′(l,m)l+12l+1jl-1(kr)Xlm-1(θ,ϕ)+aE′(l,m)jl(kr)Xlm0(θ,ϕ)-aM′(l,m)l2l+1jl+1(kr)Xlm1(θ,ϕ)).
Since
(43)lim⁡r→0jl+s(kr)rl+s=π1/22l+s+1Γ(l+s+3/2)kl+s,
this leads to the following formulas [18]:
(44)aE′(l,m)=-4πkl-12l+1l+1[(Hlm-1(∇))*⋅E0′(x)]x=0,aM′(l,m)=4πkl-12l+1l+1[(Hlm-1(∇))*⋅B0′(x)]x=0.
For Ein,Bin, we will just have the following:
(45)[(Hlm-1(∇))*⋅Ein(x)]x=0=(Hlm-1(ik))*⋅Ein(0),[(Hlm-1(∇))*⋅Bin(x)]x=0=(Hlm-1(ik))*⋅Bin(0)
and, quite similarly, for Esc,Bsc,
(46)[(Hlm-1(∇))*⋅Esc(x)]x=0=∑g[(Hlm-1(Dg))*⋅Egsc(z)]z=0,[(Hlm-1(∇))*⋅Bsc(x)]x=0=∑g[(Hlm-1(Dg))*⋅Bgsc(z)]z=0.The task is a bit more complicated forE0,B0. First, we will use (5) to rewrite (35) as
(47)J0(x)=∑l=1∞c(-1)l+1kl+1×∑m=-ll(aE(l,m)2l+1l+1Hlm-1(∇)+1kaM(l,m)Hlm0(∇)2l+1l+1)Φ(x).
Using then the fact that [18]
(48)Hlms(∇)Φ(x)=(-1)l+sσ2(l+s)Φ(x)Hlms(x)
leads to
(49)J0(x)=∑l=1∞ckl+1σ2l-2×∑m=-ll(aE(l,m)2l+1l+1Hlm-1(x)-aM(l,m)kσ2Hlm0(x)2l+1l+1)Φ(x).
If we now develop E0,B0 and J0 in vector harmonics:
(50)E0(x)=∑l,m,sElms(r)Xlms(θ,ϕ),B0(x)=∑l,m,sBlms(x)Xlms(θ,ϕ),J0(x)=∑l,m,sJlms(r)Xlms(θ,φ),
we will have
(51)Jlm-1(r)=aE0(l,m)c(2π)3/2σ2l+1kl+12l+1l+1rl-1exp⁡(k2σ22-r22σ2),Jlm0(r)=-aM0(l,m)c(2π)3/2σ2l+3kl+2rlexp⁡(k2σ22-r22σ2),Jlm1(r)=0.
According to formulas (A.3) and (A.6) of the Appendix,
(52)Elm-1(r)=-4πk2cl+12l+1×(hl-1(kr)∫0rjl-1(kr′)Jlm-1(r′)r′2dr′-jl-1(kr)∫0rhl-1(kr′)Jlm-1(r′)r′2dr′+jl-1(kr)∫0∞hl-1(kr′)Jlm-1(r′)r′2dr′)-4πikcl2l+1Jlm-1(r),Blm-1(r)=-4πk2cl+12l+1×(hl-1(kr)∫0rjl(kr′)Jlm0(r′)r′2dr′-jl-1(kr)∫0rhl(kr′)Jlm0(r′)r′2dr′+jl-1(kr)∫0∞hl(kr′)Jlm0(r′)r′2dr′).
Since Jlm-1(r)=O(rl-1) and Jlm0(r)=O(rl)(we say that f(x)=O(xn) if the limit limx→0(f(x)/xn)  exists), (52) takes the following form:
(53)Elm-1(r)=-4πk2cl+12l+1jl-1(kr)×∫0∞hl-1(kr′)Jlm-1(r′)r′2dr′-4πikcl2l+1Jlm-1(r)+O(rl+1),Blm-1(r)=-4πk2cl+12l+1jl-1(kr)×∫0∞hl(kr′)Jlm0(r′)r′2dr′+O(rl+1),
so that, according to (41),
(54)[(Hlm-1(∇))*⋅E0(x)]x=0=-aE(l,m)(2π)3/2σ2l+1l+12l+1exp⁡(k2σ22)×[∫0∞hl-1(kr)exp⁡(-r22σ2)rl+1dr+ll+12lΓ(l+1/2)π1/2ikl+2],[(Hlm-1(∇))*⋅B0(x)]x=0=aM(l,m)(2π)3/2σ2l+3kl+12l+1exp⁡(k2σ22)×∫0∞hl(kr)exp⁡(-r22σ2)rl+2dr.
Combining (44)–(46) and (54) gives the final result:
(55)aE′(l,m)=-4πkl-12l+1l+1×{(Hlm-1(ik))*⋅Ein(0)+∑g[(Hlm-1(Dg))*⋅Egsc(z)]z=0}-2πaE(l,m)(kσ)2l+1exp⁡(k2σ22)(Il-1(kσ)+il(l+1)(2l-1)!!),aM′(l,m)=4πkl-12l+1l+1×{(Hlm-1(ik))*⋅Bin(0)+∑g[(Hlm-1(Dg))*⋅Bgsc(z)]z=0}-2πaM(l,m)(kσ)2l+3exp⁡(k2σ22)Il(kσ),
where
(56)(2l-1)!!=1⋅3⋅5⋅7⋯(2l-1),Il(α)=∫0∞hl(z)exp⁡(-z22α2)zl+2dzIl(α) can easily be computed using the recursion formula:
(57)Il(α)=α2(lim⁡z→0zl+1hl(z)+Il-1)=α2(Il-1-i(2l-1)!!),
which can be deduced from an integration by parts, and using the fact that, for l=0,
(58)I0(α)=∫0∞exp⁡(iz)izexp⁡(-z22α2)z2dz=α3exp⁡(-α22)[π2+i∫0αexp⁡(t22)dt]-iα2.


## 5. Numerical Applications

In the following, we will show a few examples of calculations for illustration purpose. In particular, we would like to compare calculations based uniquely on the electric dipole terms and those obtained by including multipolar terms for better accuracy. The spherical particles are distributed in a two-dimensional hexagonal array as described previously and the metal composing the particles is gold. As particle deposits are usually done on glass, the optical refractive index supposed in the wavelength range considered (from 400 nm to 800 nm) is 1.5. The optical constants for gold have been determined by ellipsometric studies using thin films of this metal (we acknowledge with thanks Professor Georges Bader of Université de Moncton for kindly providing us with the gold optical constants).From the calculations performed with small particles, such as those with a diameter smaller than 50 nm, it is observed that the dipolar approximation with retardation effects is quite similar to more accurate calculations using higher multipolar terms. Figure1 shows the reflectance for a two-dimensional array of 50 nm diameter spherical particles on glass with an interparticle distance of 60 nm. The dipolar approximation and the multipolar calculations both have a reflectance peaking at 570 nm with only a slight variation in amplitude in the wavelength range considered.  Variations in the interparticle distance will not alter this conclusion.Figure 1
Reflectance of a two-dimensional array of 50 nm gold nanoparticles deposited on glass separated by an interparticle distance of 60 nm. Solid line: dipolar approximation, dotted line: multipolar calculations.For larger particles, the difference between the dipolar and multipolar calculations becomes considerable. Figure3 shows the case of 200 nm particles separated by a distance of 400 nm. The dipolar calculations gave a reflectance that peaks near 540 nm, whereas the reflectance given by the multipolar calculations exhibits a maximum near 570 nm. A broadening of the dipolar reflectance spectra is equally observed when multipolar terms are included. Putting the particles at a closer distance, that is, 300 nm instead of 400 nm, will accentuate the interactions between the particles and the difference between the dipolar and multipolar calculations for the reflectance (Figure 3). Compared to results in Figure 2, the curves in this latter figure show a further broadening accompanied by a red shift in the peak wavelength.Figure 2
Reflectance of a two-dimensional array of 200 nm gold nanoparticles deposited on glass separated by an inter-particle distance of 400 nm. Solid line: dipolar approximation, dotted line: multipolar calculations.Figure 3
Reflectance of a two-dimensional array of 200 nm gold nanoparticles deposited on glass separated by an inter-particle distance of 300 nm. Solid line: dipolar approximation, dotted line: multipolar calculations.It should be noted that calculations done for other metals such as silver and copper will give similar qualitative features.


## 6. Conclusion

The formulas we have obtained are relatively simple, considering the complexity of the problem encountered. When the spheres are suspended in vacuo, all that remains to do, to completely solve the problem, is to calculate explicitly the Fourier coefficients of the scattered field. This task might be a little tedious, but the results will involve nothing more elaborate than the complementary error function. From a numerical perspective, this means that in terms of special functions, the only procedures needed are those involving the computation of the erfc(x) function and the Dawson's integral which appears in (58).As for the case where a substrate has to be taken into account, a technique analogous to the one introduced in [15] can be used.



---

*Source: 901970-2012-08-14.xml*