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.
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]: where , are the usual spherical harmonics, and are the commonly known vector harmonics [19]. The more important properties of the 's are listed in the appendix of [18]. We will also use the related vector-valued functions [18]: whose components are either zero or homogeneous polynomials of degree in x, y and z, which means that . will appear in many formulas; the most convenient way we know to compute them is to first use [18] the following: then Calculations can be made in Cartesian coordinates, since the quantities 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 by a structure made of identical, nonoverlapping spheres of radius , whose centers form a hexagonal lattice in the -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 is understood for all time dependent quantities).
Step 1. We will call the field radiated by the sphere centered at , and by the total electromagnetic field incident on that sphere, which is the superposition of the incoming field and of the fields scattered by all the other spheres: Because of the symmetry of the problem, all the 's and 's are merely copies of each other, differing only by a phase factor, and thus, assuming that is the origin: Therefore,
Step 2. and appearing in (8), have the respective multipole expansions: where are the spherical Bessel functions of the first kind, and where are the first Hankel functions [19]. The multipole coefficients , (with for electric or M for magnetic) in (9) and (10) must be related by Mie's formulas: where for a given value of , the polarisabilities 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 's in terms of the 's [17]. Therefore, if we call the multipole coefficients of the incoming field, (11) could be put into a form the following: from which, theoretically, the 's can be computed, hence . In practice, however, one has to reduce (12) to an equation involving only a finite number of unknowns, which is equivalent to approximating by a field with only a finite number of nonvanishing multipole coefficients.
Step 4. Finally, we have for the scattered field the following: Again, addition theorems must be used to express the fields in (13) in terms of the '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 . Then, since we can write we will use, instead of (8), the equivalent “subtractive” formula: where, again, 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 of in terms of the multipole coefficients of the field is presented in Section 3, while in Section 4, we show how to retrieve the multipole coefficients of from the expressions for , and . Some auxiliary results are relegated to the Appendix.
3. Calculation of the Scattered Field
Let us call the current density induced in the th sphere (i.e., the source of ): we will compute the scattered field as the field radiated by the current density: Just as in (7), we must have the following: where is the incoming wave vector. Thus, (16) can be rewritten as which has the following form: times a periodic function. This means that can be developed in Fourier series [14]: where the summation is extended to all vectors of the two-dimensional reciprocal lattice. The Fourier coefficients are given by where the integral extends to the “primitive cell” of area in the -plane (note that in (20) we could have used instead of , because 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 of . To achieve that, we first solve Maxwell's equations: using a vector potential , solution of in terms of which the field is given by [14]
Next, we develop the vector potential and the scattered field in Fourier series, with respective Fourier coefficients , using formulas similar to (19) and (20). Substitution of the series into (22) and (23) gives the following: where the vector operator is defined by the following: If we also define the following: we can write the retarded solution of (24) as which gives the scattered field:
The last step consists in expressing (29) in terms of the multipole coefficients of . Any current density with the 's as multipole moments will, outside the source, radiate the same field as . Thus, if we are only interested in the value of outside the layer of spheres, we can use a point source with the same multipole moments [18]: the differential operators in (30) are obtained by replacing in by the corresponding partial derivative. This leads to the Fourier coefficients: and, after some manipulations, to the following: where should be used for and for . Using (5) and (32) can be more compactly rewritten as
4. Subtraction Process
Using a point source, like we just did, is not suitable to compute through formula (15): this gives formulas for and 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: where is assumed to be small, compared with the distance between neighboring dipoles, so that we can presume . The current density (30) and its Fourier coefficients (31) will thus be replaced by the following: with When these coefficients are substituted in (29), the corresponding formulas for 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 coefficients of using formula (15): the technique is based on results of [18], which show that if a vector field assumes the following form: one has The change of variable in formula (39) allows to rewrite it as that is, where is the gamma function.
Assuming that is small enough, will assume the general form for a free electromagnetic field [18]: Since this leads to the following formulas [18]: For , we will just have the following: and, quite similarly, for ,
The task is a bit more complicated for . First, we will use (5) to rewrite (35) as Using then the fact that [18] leads to If we now develop and in vector harmonics: we will have According to formulas (A.3) and (A.6) of the Appendix, Since and (we say that if the limit exists), (52) takes the following form: so that, according to (41), Combining (44)–(46) and (54) gives the final result: where can easily be computed using the recursion formula: which can be deduced from an integration by parts, and using the fact that, for ,
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. Figure 1 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.
For larger particles, the difference between the dipolar and multipolar calculations becomes considerable. Figure 3 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.
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() 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.
Appendix
A. Multipole Expansion of Retarded Fields inside the Source
When solving Maxwell's equations, the electrical field , magnetic induction , current density , and vector potential can all be expended using the vector harmonics : Substituting (A.1) into (22) gives first, for the retarded field the following: and when (A.2) is itself substituted into (25), the final result is
Acknowledgments
Financial support from Faculté des Études Supérieures at Université de Moncton and the Natural Science and Engineering Research Council of Canada (NSERC) is gratefully acknowledged.