Methods for streamlining expensive function calls

The focus of this thesis has been developing and validating new methods of streamlining expensive function calls. These new methods of streamlining are applied to the problem of modeling electrodynamic response in in nite two-dimensional nanoparticle (NP) arrays. Using these methods, an exact analytical solution for electromagnetic interaction with two-dimensional spherical periodic in nite arrays has been developed. The new solution o ers fundamental insight into the singularities responsible for dramatic sensitivity in- creases seen in sensory applications. The new analytical solution also provides increased computational time with 60x improvement demonstrated in this work. In addition, a new electroless deposition method followed by thermal annealing was developed by Ahn et al. for spherical Au array fabrication. To complement the model and fabrication technique, additional improvements in NP array sample analysis have been adapted and developed. Algorithms for NP particle density, size, and shape analysis are described using Atomic Force Microscope (AFM) and Scanning Electron Microscope (SEM) with a comparison between the two microscopy methods. Using fabrication and analysis, the methods of streamlining expensive function calls can be tested and improved.

METHODS FOR STREAMLINING EXPENSIVE FUNCTION CALLS by Benjamin Jeppson Taylor A thesis submitted to the faculty of The University of Utah in partial ful llment of the requirements for the degree of Master of Science Department of Chemical Engineering The University of Utah August 2011 Copyright c Benjamin Jeppson Taylor 2011 All Rights Reserved Th e Un i v e r s i t y o f Ut ah Gr a d u a t e S c h o o l STATEMENT OF THESIS APPROVAL The thesis of Benjamin Jeppson Taylor has been approved by the following supervisory committee members: Donald Keith Roper , Chair 12/23/10 Date Approved James C. Sutherland , Member Date Approved Edward M. Trujillo , Member 01/12/11 Date Approved and by JoAnn S. Lighty , Chair of the Department of Chemical Engineering and by Charles A. Wight, Dean of The Graduate School. ABSTRACT The focus of this thesis has been developing and validating new methods of streamlining expensive function calls. These new methods of streamlining are applied to the problem of modeling electrodynamic response in in nite two-dimensional nanoparticle (NP) arrays. Using these methods, an exact analytical solution for electromagnetic interaction with two-dimensional spherical periodic in nite arrays has been developed. The new solution o ers fundamental insight into the singularities responsible for dramatic sensitivity in- creases seen in sensory applications. The new analytical solution also provides increased computational time with 60x improvement demonstrated in this work. In addition, a new electroless deposition method followed by thermal annealing was developed by Ahn et al. for spherical Au array fabrication. To complement the model and fabrication technique, additional improvements in NP array sample analysis have been adapted and developed. Algorithms for NP particle density, size, and shape analysis are described using Atomic Force Microscope (AFM) and Scanning Electron Microscope (SEM) with a comparison between the two microscopy methods. Using fabrication and analysis, the methods of streamlining expensive function calls can be tested and improved. CONTENTS ABSTRACT : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : iii LIST OF FIGURES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : vi ACKNOWLEDGMENTS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : x CHAPTERS 1. INTRODUCTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1 1.1 Nanoscience Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2. THEORY AND MODELING : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 13 2.1 PNAS Submission [20] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Supporting PNAS Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Design Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3. FABRICATION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 33 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Heating Au Island Films Forms Nanoparticles . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Multiple Experiments: Drop Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.4 Experimental Results Compared to Theoretical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4. DATA ACQUISITION AND ANALYSIS : : : : : : : : : : : : : : : : : : : : : : : : : : : 40 4.1 Data Aquistion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 Watershed-Transform and Shape Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.3 Applying Watershed to SEM Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.4 Particle Diameters, SEM vs AFM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.5 Particle Densities, SEM vs AFM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5. CONCLUSION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 46 5.1 Concept to Reality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 REFERENCES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 49 LIST OF FIGURES 1.1 Coulomb's law states that the force experienced between two particles q1 and q2 is inversely proportional to the distance between them, r2, and directly proportional to the magnitude of the charges. If charges have equal signs, the force will be a repelling force, and if charges are opposite, the forces will be attractive. Coulomb's law was derived from experimental data compiled by Coulomb. Figure copyright is owned by author . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Gauss's law states that the electric ux through any closed surface is propor- tional to the enclosed electric charge. If no charges are contained within the closed surface, no electric ux can move through the surface. If no electric ux moves through the surface, the internal electric eld is zero everywhere. Figure copyright is owned by author . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Faraday's law of induction states that the electromagnetic force (EMF) around a closed loop is directly proportional to the rate of change of the magnetic ux. The gure demonstrates a simple example where a bar magnet moving through a closed loop induces a current and causes a local magnetic eld on the wire loop. This principle is responsible for electricity generated from generators which are capable of harvesting mechanical motion to produce electricity. Figure copyright is owned by author. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 SEM image of a randomly oriented Au spherical nano particle array. The array is periodic, and therefore continues randomly across the entire chip. The chemical methods used to produce this sample are outlined in detail in Chapter 3. Reprinted with permission from Langmuir [23] . . . . . . . . . . . . . . . 8 1.5 Scanning Electron Microscope (SEM) diagram showing the basic elements of a SEM. The electron gun bombards the sample with an electron beam and is focused and moved with magnetic lenses and elds. The sample is typically placed within a high vacuum environment in order to limit the additional noise that an atmosphere induces. The primary electrons from the beam are scattered according to the topography and material of the sample. An electron collector collects the scattered electrons and ampli es the signal with a photomultiplier tube. The resulting signal is converted into a raster image as the electron beam scans the sample side-to-side. [23] Figure copyright is owned by author . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.6 Atomic Force Microscope (AFM) diagram showing basic elements of an AFM. A cantilever tip moves along the top of the sample in order to detect sample topography. As the cantilever moves up and down, a laser re ects o the cantilever arm to detect movement via a photo diode. A raster image is made from the electrical signal detected by the photo diode as the arm is scanned side-to-side. Reprinted with permission from Wikipedia [3] . . . . . . . . . . . . . . . . 11 2.1 Submitted to Proceedings of the National Academy of Sciences of the United States of America Grid layout for a 5x5 square grid where rij represents the vector between the center point, i, and the point, j. Dotted and solid circles are shown to represent the nested grid concept used to obtain the analytical solution for the retarded dipole sum, S, for a square grid array. . . . . . . . . . . . . 22 2.2 Submitted to Proceedings of the National Academy of Sciences of the United States of America The retarded dipole sum, S , may also be visualized as the three contributing equations (2.30,2.31,2.32) . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3 Submitted to Proceedings of the National Academy of Sciences of the United States of America The retarded dipole sum, S, contribution can be split into static (r1), intermediate (r2), and far- eld (r3) contributions for greater understanding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Submitted to Proceedings of the National Academy of Sciences of the United States of America The singularities in the retarded dipole sum, S, can be seen at the expected locations where A = p2 and A = 2 . An in nite number of singularities exist for an in nite array but due to the D3 factor in (2.29) singularities resulting in larger A are less signi cant . . . . . . . . . . . . . . . . . . . . . 24 2.5 Submitted to Proceedings of the National Academy of Sciences of the United States of America Computed values of the retarded dipole sum in (27) for nite arrays of size L converge to exact solution for retarded dipole sum in (2.29) for in nite array. As L increases, real and imaginary parts of (2.27) red and blue shift, respectively. For values of L > 20, both real and imaginary parts exhibit narrowed lineshapes. The exact solution is computed using a smaller grid size (400) than in Figure 4 (1000), which increases the oscillation in the computed values. Inset: Computation time for nite arrays increases exponentially with L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Submitted to Proceedings of the National Academy of Sciences of the United States of America The behavior of the SIII double sum is shown with the interpolation points. This simulation used 10002 NPs which produces oscillation 29 2.7 (Table and caption reproduced with permission from IEEE) Near- eld plasmon features (<400nm) with far- eld extinction peaks (>450 nm) for 50 nm radius Ag NP array embedded in air. The di erent colored lines represent di erent interparticle spacings and the black line is spectra for an in nite spacing. In nite spacing can also be thought of as a single particle spectra. . . . . . . . . . 31 vi 2.8 (Table and caption reproduced with permission from IEEE) Sensitivity of AU NP extinction e ciency (EE) with respect to change in local refractive index, o, where sensitivity is inversly proportional to the peak full width half max (FWHM) and de ned as the maximum extinction e ciency delta. A near eld subset from the inset plot shows e ects of changing particle radius, r, and interparticle spacing d. The colorscale shows sensitivity values up to 188. Dashed red lines show where relative separation distance, d=2r, is constant as r and d are varied. The insert plot shows extraordinary far- eld e ects predicted by the new algorithm. Dashed lines in the insert plot show dis- continuities predicted using the new algorithm with adjacent high sensitivity. This plot allows the search space to be viewed for the optimization problem, and demonstrates the predictability of the singularities. When optimizing this problem in practice for higher dimensions the results are not visualized, but optimized using a metaheuristic global search algorithm. . . . . . . . . . . . . . . . . . 32 3.1 (Figure and caption reproduced with permission from Langmuir [23]) SEM images of fuzzed quartz slides after immersion ( rst row) in Au plating solution for (A) 0.33 min, (B) 1 min, and (C) 4 min; after respective immersed slides are annealed (second row, images D, E, and F) for 3 h at 250 C; and after respective immersed, annealed slides are heated (third row, images G, H, and I) for 20 min at 800 C. Inset are transmission UV-vis spectra corresponding to each image. Percent transmission is shown in images A-C, E, and F. Extinction in arbitrary units (au) is shown in images D, G-I. All images are the same magni cation; see scalebar on the left. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 (Figure and caption reproduced with permission from Langmuir [23]) SEM images of fused quartz slides annealed (250 C; 3h) and heated (800 C; 20 min) after immersion in Au plating solution for (a) 8 min, (b) 30 min, and (C) 60 min. Weblike strands are NaCl-lie crystal structure removed after rinsing with DDD-H2O which were identi ed using X-Ray Di raction (XRD). . . . . . . 35 3.3 (Figure and caption reproduced with permission from Langmuir [23]) SEM images after annealing (250 C, 3 h) fused quartz slides exposed to microliter droplets of EL plating solution exposed for (A) 0.33 min, (B) 1 min, and (C) 2 min. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Low vacuum SEM image showing increased Au deposits at the drop boundary as the drop evaporated from N2 convection during the deposition time . . . . . . 36 3.5 Relative pixel intensity from SEM image of droplet edge shows increased Au deposits near the boundary and the change in Au coverage uniformity moving towards the center of the particle. The vertical dashed line represents the best approximation of the outside droplet boundary. Relative slopes at the beginning and end of the plot should not be taken literally since charging may be e ecting perimeter intensity in the Field Of View (FOV). Figure copyright is owned by author . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 vii 3.6 Theoretical results deviate from the experimental results for several reasons. The experimental results are created with randomly oriented arrays with a particle size distribution. This violates the equal polarization assumption used in the theoretical derivation (2.15). In order to perform a more rigorous appraisal of the theoretical results, a lithography method of fabrication would be needed to ensure equal spacing and particle size. Figure copyright is owned by author . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1 (Figure and caption reproduced with permission from Langmuir [23]) The watershed method distinguishes particle characteristics more realistically than simple threshold methods when particles are grouped in clusters. (A) Original SEM image of Au-coating on quartz slide. (B) Threshold method of segmen- tation results in particle "fusion" with size distribution skewed by particles at the border. (C) Watershed method with border removal quantitatively distinguishes particle size, shape, and number density. . . . . . . . . . . . . . . . . . . . 41 4.2 Each particle was centered by nding the individual centers of mass and aligned by using the golden search method to maximize the moment about a particular axis. The resulting 3D histogram allows for mean shape determi- nation and shape con dence plots analysis. A shape con dence plot would be an outline shape that contained a quantile of the particle population as shown by the two black contour traces. Figure copyright is owned by author . . . . . . . 42 4.3 (Figure and caption reproduced with permission from Langmuir [23]) Mean diameter and density of particles obtained after annealing at 250 C (3 h) and heating at 800 C (20 min) shift to larger scales as E.L. deposition time increases from (A) 0.33 to 1 to 4 min. (B) Mean diameter and particle density for annealed and heated lms deposited at 8, 30, and 60 min. . . . . . . . . . . . . . 43 4.4 (Figure and caption reproduced with permission from Langmuir [23]) Particle densities estimated by watershed analysis of AFM ( lled circles) and SEM (open circles) images for immersion times at 0.33, 1, 4, 8, 30, and 60 min. The 0.33 min sample was damaged prior to AFM imaging. . . . . . . . . . . . . . . . . . . . 45 viii ACKNOWLEDGMENTS The author would like to acknowledge Dr. Randy Polson at the Dixon Laser Institute of the University of Utah, Salt Lake City, Utah, USA, for technical assistance with AFM, and also Analia G. Dall' Asen and Wonmi Ahn for their close collaboration in the lab and data analysis, and Dr. Keith Roper for his support and guidance. CHAPTER 1 INTRODUCTION 1.1 Nanoscience Introduction As scienti c research continuous to nd bottlenecks in computing power, care must be taken to ensure maximum e ciency within the algorithms themselves. This work demonstrates a novel method of streamlining computationally expensive function calls to improve simulation results. This concept is demonstrated with an algorithm used for sim- ulating electrodynamic response of two-dimensional nanoparticle (NP) arrays. Accurately simulating the electric eld associated with these dipole array sensors allows for sensor design optimization and increased sensitivity. By increasing the sensor sensitivity, antigen detection on the order of a single molecule may be possible. In order to validate the new results, NP array fabrication as well as metrology techniques were used. For metrology, this work covered UV-VIS spectroscopy, Scanning Electron Microscopy (SEM), and Atomic Force Microscopy (AFM) as well as new image processing algorithms that were needed to generate summary statistics. The new analytical solution that was derived using this redundancy elimination technique o ers more fundamental insight. Now the exact loca- tions of solution singularities can be identi ed without numerical techniques or expensive simulations as was previously done. In addition to more insight, speed improvements on the order of 60x are also realized which allow more particles, and therefore more accuracy per simulation. 1.2 Modeling 1.2.1 History of Electromagnetism Mankind has always been aware of electromagnetism with the dramatic e ects of light- ning and only during the past millennium, have we nally begun to gain a complete understanding of the physics behind this fascinating phenomenon. An uncon rmed Greek legend dates the discovery of magnetism back to 900 B.C., when a Greek shepherd named Magnus experienced a force pulling on the iron nails in his sandals while crossing a bed 2 of black rock [18]. The region was named Magnesia and the black rock became known as magnetite. Electrostatic potential e ects were rst documented by a Greek philosopher by the name of Thales who noticed that amber could pick up a feather after being rubbed by cat fur. The rst widespread use of electromagnetism came around 1100 AD with the development of the compass which relied on using a magnetized needle and the earth's magnetic eld for direction. 1.2.2 Coulomb's Law During the past millennium, hundreds of scientists have contributed to the develop- ment of electromagnetic understanding with a plethora of experimentation. One such experimenter, Girolamo Cardano, published De Substilitate (1550) [18] which distinguished for the rst time the di erences between magnetic and electric forces. One of the rst quantitative expressions was later developed by a French physicist, Charles Augustin de Coulomb. Coulomb developed an expression from experimental data that allowed him to predict the force, F, between two charged electrostatic objects [18]. F = ke q1q2 r2 Coulomb's law (1.1) The magnitude of the force, F, was directly proportional to the magnitude of a charge, q1, in the presence of a second charge, q2, as seen in Fig 1.1. As the two charges q1 and q2 were brought closer together, Coulomb noticed that the force was inversely proportional to the distance. The proportionality constant ke also known as Coulomb's constant was derived from the experimental data but was later shown to have more signi cance than a simple constant. 1.2.3 Gauss's Law Despite the usefulness of Coulomb's law, it was not a general expression yet, but valid only for stationary point charges. It was not until 1835 that Carl Friedrich Gauss derived a gerneral expression called Gauss's law [18]. The new expression relied on the concept of an electric eld which was rst introduced by Michael Faraday. The electric eld is de ned as the space surrounding an electric charge or moving magnetic eld. An electric eld represents the force that a hypothetical point charge would experience in a given location. With the new concept of an electric eld, Gauss's law stated that the electric ux through any closed surface was proportional to the amount of electric charge enclosed by that surface. 3 +q2 +q1 +q2  -q1      r Figure 1.1. Coulomb's law states that the force experienced between two particles q1 and q2 is inversely proportional to the distance between them, r2, and directly proportional to the magnitude of the charges. If charges have equal signs, the force will be a repelling force, and if charges are opposite, the forces will be attractive. Coulomb's law was derived from experimental data compiled by Coulomb. Figure copyright is owned by author The electric ux can be de ned as the number of electric eld lines moving through an arbitrary surface. Two stronger charges will experience a higher electric ux between the charges since there will be a higher number of equipotential lines. This law can be explained by the expression Integral form: I S E @A = Q 0 (1.2) Di erential form: ~r E = 0 (1.3) where 0 is the material speci c permittivity, Q is the total charge, and is the charge density. The di erential element, @A, is a vector with the magnitude of an in nitesimal surface area element with direction perpendicular to the surface, as seen in Fig 1.2. The divergence1, ~r , is used to represent the same expression in di erential form allowing the bene ts of vector calculus to be used. 1.2.4 Gauss's Law of Magnetism Similar to Gauss's law, an analogous expression was developed for the magnetic eld which stated that the total magnetic ux through a closed three-dimensional surface was zero 1The divergence is a vector calculus operator which determines the magnitude of a vector eld's source or sink. For example, the divergence of the magnetic eld, ~B , can be de ned by: ~r ~B = ~X @Bx @x +~Y @By @y +~Z @Bz @z for Cartesian coordinates 4          E E E ∫sĖdA = εQ0 dAθ Figure 1.2. Gauss's law states that the electric ux through any closed surface is proportional to the enclosed electric charge. If no charges are contained within the closed surface, no electric ux can move through the surface. If no electric ux moves through the surface, the internal electric eld is zero everywhere. Figure copyright is owned by author Integral form: I S B @A = 0 (1.4) Di erential form: ~r B = 0, (1.5) where B is the magnetic eld and S represents any closed surface enclosing the magnetic eld. Unlike Gauss's law for the electric eld, magnetic charges are assumed to come in pairs, also known as dipoles. The theoretical existence of a single magnetic point charge, referred to as a magnetic monopole, has been debated with no proof on either side of the argument until recently. New research published in Nature September 2009 on spin ice reports quasiparticles resembling the behavior of magnetic monopoles [4]. Gauss's law in the standard set of Maxwell equations (1.8-1.11) excludes the possibility of such a particle existing. Despite this recent discovery, the assumption of no magnetic monopoles has held 5 S N motion wire loop EMF B B B B Figure 1.3. Faraday's law of induction states that the electromagnetic force (EMF) around a closed loop is directly proportional to the rate of change of the magnetic ux. The gure demonstrates a simple example where a bar magnet moving through a closed loop induces a current and causes a local magnetic eld on the wire loop. This principle is responsible for electricity generated from generators which are capable of harvesting mechanical motion to produce electricity. Figure copyright is owned by author. up experimentally for nearly a century and is considered a valid assumption for the vast majority of cases and therefore is excluded in this work. 1.2.5 Faraday's Law An additional breakthrough for electrodynamic understanding came from Faraday's law of induction which was also developed from experimental data. Faraday's law of induction states [17] that the induced electromotive force or EMF in any closed circuit is equal to the time rate of change of the magnetic ux ( @B @t ) through the circuit, as shown in Fig 1.3. Faraday's law: ~r E = 􀀀 @B @t (1.6) A similar but opposite law was developed called Ampere's law which related a magnetic eld to its current source. In Faraday's law, a moving magnetic ux created a current, 6 where in Ampere's law, a moving current created a magnetic eld. 1.2.6 Ampere's Law Ampere's law relates the magnetic eld outside a closed loop to the current passing through it. The law was derived in 1826 by Andr-Marie Ampre [18]; however, the orig- inal formulation did not account for time varying charges. Maxwell added the missing displacement current correction to Ampere's law in 1861 to get Ampere's law: ~r B = J + @E @t (1.7) where and are material-speci c permittivity and permeability, respectively, and ~r is the curl2. In equation (1.8-1.11) J, and t represent total current density, charge density, and time, respectively. 1.2.7 Maxwell's Equation In 1864, the greatest single step in the world's knowledge of electricity was made by James Clerk Maxwell [12]. Maxwell, who had compiled the current knowledge of electromagnetic phenomena into twenty di erential equations, hypothesized that electric elds traveling at the speed of light was not coincidence. Maxwell's hypothesis that light and electromagnetic elds were related led to his mathematical derivation that light was an electromagnetic wave with oscillating magnetic and electric elds. Maxwell was also able to calculate the speed of light using constants from electrostatic experiments to be 3:107x108m=s (actual speed is: 2:997x108). Despite the majority of all current electromag- netic understanding stemming from Maxwell's equations, the general equations that are used today were simpli ed by Oliver Heaviside [18] to four coupled di erential equations: 2The curl is a vector calculus operator which de nes the rotation of a vector eld. For example, the curl of the electric eld, ~E, in Cartesian coordinates can be de ned by: ~r E = ~X @Ez @y 􀀀 @Ey @z +~Y 􀀀 @Ex @z 􀀀 @Ez @x + ~Z @Ey @x 􀀀 @Ex @y 7 Gauss's law: ~r E = (1.8) Faraday's law: ~r E = 􀀀 @B @t (1.9) Ampere's law: ~r B = J + @E @t (1.10) Gauss's law of magnetism: ~r B = 0 (1.11) Even though Maxwell's equations can be solved directly for simple situations, it is often convenient to introduce the idea of potentials in order to obtain a smaller number of second- order equations. From Maxwell's equations, the electromagnetic expression for light can be derived[7] to be: ~r 2 􀀀 1 c2 @2 @t2 ~E = ~r 2 􀀀 1 c2 @2 @t2 ~B = 0 (1.12) where c is the speed of light. This expression for light allows us to understand and predict the dramatic e ects of light interacting on the nano-scale. 1.2.8 Nanoparticle Arrays Recent attention in nanoscience has focused on the potential uses found in periodic arrays of metal nanoparticles. Silver (Ag) and Gold (Au) arrays have been shown to have dramatic sensitivities to environmental refractive index changes making them ideal candidates for biosensors. This work focuses speci cally on the experimental and theoretical results of spherical Au arrays, as shown in Fig 1.4. Using similar nanoparticle arrays, single cell analyte detection is possible due to rapid changes in the elctromagnetic behavior of such arrays. 1.2.9 Nanoparticle Array Modeling In x2.1-2.4, theoretical methods are discussed for predicting the behavior of light with the nanoparticle (NP) arrays. Numerical methods are available for solving the coupled partial di erential methods which make up Maxwell's equation. However, high computational costs usually limit a direct approach due to the number of particles and the standard numerical methods available which require spacial discretization. Other numerical methods which are faster and more commonly used to model the behavior of light interacting with NP arrays include Finite Di erence Time Domain (FDTD) [19], T-matrix [14], and 8 Figure 1.4. SEM image of a randomly oriented Au spherical nano particle array. The array is periodic, and therefore continues randomly across the entire chip. The chemical methods used to produce this sample are outlined in detail in Chapter 3. Reprinted with permission from Langmuir [23] Coupled Dipole Approximation (CDA) [21]. FDTD and T-matrix are considered exact multipole methods because they account for all signi cant poles. The CDA, which is only a dipole approximation, can be valid for smaller NP sizes 100 nm. Despite the relative computational speed improvement using the CDA approach outlined by Schatz [22], optimization of design parameters is still limited. A new analytical solution derived in this work can now provide precise locations for all singularities for a square grid array for particle sizes 100 nm for nonmagnetic particles. Singularity design criteria are of particular interest because they are usually related with enhanced sensitivity characteristics. The new solution, derived in Section x2.2, provides more fundamental insight and the computational e ciency necessary to allow true multivariable optimization of all design parameters: NP spacing, NP size, NP material, medium material, and wavelength used. 9 1.3 Fabrication 1.3.1 Nanoparticle Array Fabrication Chapter 3 introduces fabrication of Au NP arrays on Silicon (Si) substrate. Various methods of Au NP array fabrication have been developed including ion implantation, surface etching [16], laser ablation [9], inductively coupled thermal plasmas [11], thermodynamic self organization [8], and thermal evaporation of thin Au lm followed by annealing [6]. In addition to NP array fabrication, Au thin lms and island lms (thin lms with voids) have been created using electrodeposition, electroless deposition (EL) [15], sputtering [10], and vacuum evaporation [5]. Current methods for fabrication tend to be limited by geometry, material, and cost. In collaboration with Ahn et al. [2] a new EL method has been developed, as outlined in Chapter 3, which allows for tunable NP size distributions as well as the ability to coat internal 3D surfaces at ambient conditions. This new EL method is also preferred since it does not require expensive equipment and has the ability for rapid coating and fabrication making it a potential economical candidate for large scale NP array fabrication. 1.4 Analysis 1.4.1 Scanning Electron Microscope A Back Scattered Electron (BSE) SEM is popular for nano-scale analysis. It works by using a high-powered electron beam that emits electrons down onto the sample. The sample is typically placed in a high-vacuum environment in order to limit additional noise that an atmosphere would induce. The electron beam is moved with angstrom (1x10􀀀10 meters) precision using magnetic elds to focus and scan the beam in a raster (side-to-side) pattern. As shown in Fig 1.5, while the beam is scanning, electrons from the beam are scattered and collected by a BSE detector. Within the detector, the electrons enter a photomultiplier tube to amplify the signal which is then converted into an image. The SEM is susceptible to image noise and therefore, additional steps are needed which include averaging multiple scans and applying digital lters [3]. BSE SEM o ers a powerful tool capable of obtaining topographical sample dimensions, material composition, and electrical conductivity. 1.4.2 Atomic Force Microscopy In addition, the AFM [1] has been developed to get topography information on samples using a mechanical cantilever arm that is moved across the surface in a raster pattern; see 10 Electron Gun Electron Beam Computer Generated Image Anode Magnetic Lens Scanning Coils Stage Specimen BSE Detector Photo-Multiplier Tube Figure 1.5. Scanning Electron Microscope (SEM) diagram showing the basic elements of a SEM. The electron gun bombards the sample with an electron beam and is focused and moved with magnetic lenses and elds. The sample is typically placed within a high vacuum environment in order to limit the additional noise that an atmosphere induces. The primary electrons from the beam are scattered according to the topography and material of the sample. An electron collector collects the scattered electrons and ampli es the signal with a photomultiplier tube. The resulting signal is converted into a raster image as the electron beam scans the sample side-to-side. [23] Figure copyright is owned by author Fig 1.6. As the microscopic lever moves up and down over the particles or changes in the sample topography, an electrical signal is detected via a laser and converted into an image. The AFM is considerably slower than the SEM due to the mechanical movement of the arm. The AFM also lacks resolution relative to the SEM. However, the AFM does o er topographic data that can complement the SEM analysis. A detailed comparison between the data sets o ered with SEM and AFM is presented later in this work in Chapter 4. 11 Figure 1.6. Atomic Force Microscope (AFM) diagram showing basic elements of an AFM. A cantilever tip moves along the top of the sample in order to detect sample topography. As the cantilever moves up and down, a laser re ects o the cantilever arm to detect movement via a photo diode. A raster image is made from the electrical signal detected by the photo diode as the arm is scanned side-to-side. Reprinted with permission from Wikipedia [3] 1.4.3 Nanoparticle Sample Analysis Chapter 4 reviews common methods for NP data acquisition and related image pro- cessing algorithms used for analysis. Popular tools for NP data acquisition include AFM and SEM images which have been previously mentioned. However, AFM tip convolution produced by the AFM probe tip can cause exaggerated particle diameters if improper algorithms are used. A detailed comparison between the data produced by AFM and SEM is given in Chapter 4 to provide a better understanding between these two methods. Once the data have been collected, it can be analyzed using various image processing algorithms. Commonly used algorithms include threshold particle segmentation, which uses a particular gray scale shade to designate particles from silicon. The threshold method is simple but can also create exaggerated particle sizes due to clumping multiple particles into one. Superior methods include using the watershed transform [13] which is capable of segmenting adjacent Sample Surface , , 12 particles. In addition to particle segmentation, an additional algorithm has been developed which allows particle shape analysis (see 4.1.1). 1.5 Conclusion Making use of thousands of years of electromagnetic understanding leading to Maxwell's equation, all macroscopic electrodynamic phenomena can be explored. Individual equations and laws previously derived have been scrutinized and re-derived to produce a sound set of equations that have held up to experimental data. Exciting enhancements at the nano scale can be modeled using computational approximation and in some cases analytical solutions derived from Maxwell's equations for a greater understanding. In addition to modeling, Au NP arrays can be fabricated and analyzed to provide greater understanding, theoretical validation, and feedback on model assumptions. The Greek pre x nano is becoming more well known as researchers and scientists break new ground in the nanoscience area using powerful tools of function call reduction discussed in this work. CHAPTER 2 THEORY AND MODELING 2.1 PNAS Submission This chapter was submitted to the Proceedings of the National Academy of Science. 14 Exact solution to the coupled dipole approximation of Maxwell's equations for electromagnetic interactions in periodic nanostructures Ben J. Taylor¤, Braden M. Harbiny, and D. Keith Roper¤ y ¤University of Utah Department of Chemical Engineering, and yUniversity of Arkansas Ralph E. Martin Department of Chemical Engineering Submitted to Proceedings of the National Academy of Sciences of the United States of America An exact solution to the coupled dipole approximation of Maxwell's equations for electromagnetic interactions in periodic subwavelength nanostructures has been derived. Contributions to the retarded dipole sum from near-¯eld, intermediate, and far-¯eld radiative interactions are distinguished. The solution precisely identi¯es the source, location, and magnitude of periodic divergences in real and imaginary components of the retarded dipole sum which arise due to the wavenumber and grid size, independent of particle size and dielectric constants of particle and environment. Far-¯eld radiative contributions are signi¯cant, outweighing near-¯eld and intermediate e®ects. Real and imaginary parts computed for retarded dipole sum in ¯nite arrays converge toward the exact solution as array size increases. The exact solution permits validation of approximations to retarded dipole sum using ¯nite di®erence time domain and discrete dipole approaches. The new solution allows rapid optimization and design of novel optoelectronic, spectroscopic, biomedical and sensing systems which use nanoparticle arrays as well as integrated circuit waveguides, optical ¯lters, photovoltaics and antennae. Section 1: Introduction Periodic subwavelength nanostructures such as one- and two-dimensional arrays of holes or particles exhibit extraordinary spectral resonances in theory [1, 2, 3, 4, 5, 6, 7, 8] and experiment [9, 10, 11, 12, 13, 14, 15]. Periodic nanostructures allow superluminal propagation [16] and could improve surface enhanced spectroscopies [17, 18, 19], sensors [20, 21], photodynamic therapies [22], subdi®ractive mapping of biomolecules [23], integrated circuit waveguides [24], optical ¯lters, photovoltaics and antennae. Lattice dimension [7] and material size, shape, and dielectric functions can be selected to phase-match incident and di®ractive coupled ¯elds and induce large polarizabilities that enhanced local surface ¯elds and produce narrow, high-amplitude resonances at tunable frequencies. Control of strong local EM ¯elds has recently permitted nanoparticle (NP) tuning of magnetic resonance frequency and lineshape of split ring resonator metamaterials [25] and modulation of polarization of light scattered from NP [26]. Rapid development of spectroscopically-active periodic subwavelength nanostructures has focused recent attention on solving Maxwell's equations for electromagnetic (EM) interactions in nanoscale lattices. The complete problem has been analyzed computationally using discrete dipole approximation (DDA), ¯nite di®erence time domain (FDTD), and T-matrix. Some limiting cases of the problem have been solved analytically to obtain exact expressions. However, existing exact expressions for two- dimensional (2D) arrays neglect contributions to local electromagnetic ¯eld from far-¯eld radiation. Solutions to classical EM theory are necessary, to describe optical properties of multi-NP systems such as surface enhanced Raman spectroscopy (SERS) [27] especially in combination with quantum mechanical methods for electric ¯elds very close to NP surfaces [28]. Exact solutions improve matching of atomistic and continuum-scale methods to allow multiscale characterization of such phenomena. The DDA [29, 30] allows scattering and absorption of an incident EM ¯eld to be computed for dielectric structures of arbitrary geometry and composition from induced dipole moments that arise from interactions in a set of point dipoles representing a scatterer. DDA results for nanoscale lattices have been validated using exact solutions to Maxwell's equations for two contiguous spheres [31]. DDA results were comparable with T-matrix (i.e. coupled multipole) simulation [24] of scattering from arrays based on Mie theory using a truncated multipolar expansion around each particle and the translation addition theorem for vector harmonics [3]. DDA accurately predicted spectral wavelength shifts of Ag cylinders [11]. Compared with FDTD, DDA computations are generally applicable, easier to use, faster, more accurate and more sensitive to local dielectric constant. But computational cost becomes prohibitive when solving the number of equations required for large arrays, and when varying multiple parameters (e.g. particle size, spacing, and dielectric constant) [31]. DDA computation using the package DDSCAT [17, 18, 19] inaccurately predicted values of polarization in large particle arrays for interparticle spacing equal to the wavelength of light ¸ [8]. While some numerical methods could take advantage of symmetry in the scatterer, symmetry has not been used to simplify DDA computation. The coupled dipole approximation (CDA) is formally similar to DDA, except that each scatterer is represented by a single dipole whose polarizability is equal to that of the particle. A lattice dispersion formula may be used to make the dipole array response identical to that of the bulk continuum. T-matrix provides exact solutions to the CDA that allow attribution and quanti¯cation of resonant-energy transmission and absorption features. These features are determined by material dielectric properties, nanostructure dimensions and periodicity on the basis of limiting condition sets. Computational results that used Reserved for Publication Footnotes www.pnas.org/cgi/doi/10.1073/pnas.0709640104 PNAS Issue Date Volume Issue Number 1{8 15 CDA for in¯nite arrays of NP were shown to produce spectral features consistent with those from T-matrix for interparticle spacings, D, greater than 2R in NP of radius, R, equal to 30nm and 5 nm, respectively [3].The adequacy of CDA may decrease in describing NP > 100 nm in which multipolar dynamic depolarization and radiative damping decrease scattering. Exact expressions for induced dipole ¯elds in periodic arrays of in¯nite dimensions have been derived for one-dimensional (1D) dipole chains [2] and for 1D chains of cylinders [8] in the long-wavelength limit (rij ¿ ¸) and for 2D dipoles neglecting far-¯eld radiation, where rij is the spacing between between two particles i and j. Merchiers [32] derived expressions for far- ¯eld cross sections of 2D arrays as functions of induced electric and magnetic dipole moments. Comparing analytic dipole and computational multipole approaches shows multipole e®ects produce large deviations in spectral peakshifts only for jrij j=R < 2:5 [33]. This validates development of exact expressions that neglect multipole e®ects in this range. Experimental studies of periodic subwavelength nanostructures has typically been limited to disk shaped NP arrays made by electron-beam lithography [34]. This is due to di±culties in fabricating stable arrays of spheres. These studies show increas- ing interparticle spacing (i.e. grating constant, jrij j) increases transverse-mode resonance peak and scattering intensity while decreasing Lorentzian peakwidth. This increases Q-factor due to ohmic loss from electron-phonon scattering that occurs instead of radiative loss. The purpose of this work has been to derive an exact solution to the CDA of Maxwell's equations for electromagnetic inter- actions in periodic dipole arrays that include contributions due to far-¯eld radiation. This exact solution provides a benchmark for resource-intensive computational analysis. It is applicable to analyze superluminal propagation [2]. An exact solution makes it possible to extend results recently computed from Maxwell's equations via generalized Mie theory (GMT) [26, 27], DDA [28], and FDTD [22] in order to describe modulation of polarization, enhanced EM ¯elds (SERS), and photodynamic therapy, respectively, and concepts such as NP sensing networks [20, 21], light-induced self assembly [35], and subdi®ractive mapping of biomolecules [23] examined in small, random NP clusters to regular, two-dimensional NP systems. It allows interpretation of phenomena observed in physical subwavelength nanostructures of increasing complexity. It can be used to guide parameter selection in design and manufacture of integrated devices. Examples summarized in this and previous paragraphs represent some of many possible speci¯c applications for the exact solution. Section 2: Model The phenomena of classical electricity and magnetism interacting with physical objects can be completely described on a macro- scopic scale using Maxwell's equations and the Lorentz Force law, F = q(E + v £ B), where F is the force experienced by a particle, q is the charge on the particle, E is the electric ¯eld, v is the particle velocity, and B is the magnetic ¯eld. In a homogeneous, isotropic, linear medium in which magnetic monopoles are absent the Maxwell equations become ~r¢E = ½ ² [1] ~r£E = ¡ @B @t [2] ~r£B = ¹J + ¹² @E @t [3] ~r¢B = 0 [4] where ² and ¹ are material-speci¯c permittivity and permeability, respectively, and ~r£ is the curl. In (1-4) J, ½ and t represent total current density, charge density, and time respectively. Even though Maxwell's equations can be solved directly for simple situations it is often convenient to introduce the idea of potentials in order to obtain a smaller number of second-order equations. For the case of nonmagnetic material we can neglect1 the magnetic ¯eld Maxwell equations. The electric ¯eld may be de¯ned in terms of a vector potential, A, when there are no sources (½ = 0 and J = 0) using E = ic k ~r £ (~r £ A) [5] where k = !p¹²=c = 2¼=¸ is the wavenumber, ! represents angular frequency, ¸ is wavelength, c is the speed of light, and i = p¡1. In Jackson [37] the sinusoidal time varying retarded vector potential is de¯ned as A(x) = ¹ 4¼ eikr r 1X m=0 (¡ik)m m! Z J(x0)(n ¢ x0)md3x0 [6] where x is the point where the ¯eld is measured, x0 corresponds to the source of the E ¯eld, and n is the unit vector between the points x and x0. The vector potential for a dipole can be derived by keeping only the 0th term (m=0) in (6) to produce A(r) = ¡ ¹ 4¼ eikr r i!p [7] 1The e®ective ¹ may deviate from unity [36] when the second Mie coe±cient, b1, is no longer negligible for nonmagnetic materials (e.g. Au, Ag, ...). 2 www.pnas.org/cgi/doi/10.1073/pnas.0709640104 Footline Author 16 where r = jx¡x0j is the displacement vector and p is the electric dipole moment de¯ned by p = R x0½(x0)d3x0, and r = jrj. The dipole vector potential can be substituted into the electric ¯eld de¯ned in terms of a vector potential to produce [see Supporting Information: Appendix A] E(r; p) = eikr 4¼² ½µ 1 r3 ¡ ik r2 ¶ [3(^r ¢ p)^r ¡ p] + k2 r [(^r £ p) £^r] ¾ [8] which is the general vector form of the electric ¯eld of a retarded oscillating dipole where ^r is the unit vector in the r direction, de¯ned by ^r = r=jrj. We now proceed to use the de¯nition of the dipole electric ¯eld (8) to model the behavior of light interacting with a nanoparticle (NP) array. The local electric ¯eld, Eloc, experienced by an individual NP upon which light is incident is made up of a contribution due to the oscillating electric ¯eld, E, from neighboring NPs given by (8) and a contribution due to incident light, Einc. The electric ¯eld due to the incident light can be derived from Maxwell's equations when there are no free sources and ² and ¹ are assumed constant to produce a wave equation, ~r 2E ¡ ¹² @2E @(ct)2 = 0 [9] whose solution is Einc = Eoeik¢r0¡i!t [10] where r0 is a position in the wave¯eld and k is the wavevector. Using the principle of superposition the local electric ¯eld experienced by any individual NP can be expressed as Eloc = Einc + X j6=i E(rij ; pj) [11] where the subscripts i and j denote a vector between the origin, i, and the neighboring particle, j. The expression for the local electric ¯eld can be used to relate polarization, pi, to the local ¯eld using a scalar polarizability, ®i, viz. pi = ®iEloc. [12] For the case when the incident light is perpendicular to the array in the ^x-^y plane and polarized in the ^x direction and the array is symmetric about both ^x and ^y (8) can be used to produce [see Supporting Information: Appendix A] X j6=i E(rij ; pj) = X j6=i eikrij 4¼² ½µ 1 r3 ij ¡ ik r2 ij ¶ (3 cos2(µij) ¡ 1) + k2 sin2(µij) rij ¾ pj [13] where µij is the angle between the position vector, ^rij, and the polarization, pj, and rij = jrijj. Now after combining (13), (11) into (12) we obtain pi = ®iEinc + ®i 4¼² X j6=i eikrij 8>>>>< >>>>: static (SZ) z }| { 3 cos2(µij) ¡ 1 r3 ij ¡ intermediate (IZ) z }| { ik(3 cos2(µij) ¡ 1) r2 ij + radiation (RZ) z }| { k2 sin2(µij) rij 9>>>>= >>>>; pj. [14] In (14) contributions to the oscillating dipole electric ¯eld arise from near-¯eld static zone (SZ;1=r3), induction zone (IZ;1=r2), and far-¯eld radiation zone (RZ;1=r). In the case of an in¯nite array, pi = pj = p and (14) can be solved for p yielding p = ®Eloc = Einc 1 ® ¡ S(k) = Eoeik¢r¡i!t 1 ® ¡ S(k) [15] where S(k) is the retarded dipole sum given by: S = 1 4¼² 1X j=1 j6=i fi;j(k) = 1 4¼² 1X j=1 j6=i 8>>>< >>>: SZ z }| { 3 cos2(µij) ¡ 1 r3 ij ¡ IZ z }| { ik(3 cos2(µij) ¡ 1) r2 ij + RZ z }| { k2 sin2(µij) rij 9>>>= >>>; eikrij . [16] The retarded dipole sum, S(k), corresponds to the sum of the electric dipole ¯elds at any array node divided by the polarization, p. The local electric ¯eld, Eloc, in the in¯nite array is equivalent at every array node. Its magnitude and direction varies periodically from node to node throughout the in¯nite array. The function, fi;j(k), in (16) can be transformed to a function fi;j(k;D) in rectilinear coordinates using rij = D p x2 + y2 and angle µij = tan¡1(y=x). In this transformation, variables x and y are integer multiples of intraparticle spacing, D, for a square array along the corresponding coordinate axes ^x and ^y (See Fig 1). Footline Author PNAS Issue Date Volume Issue Number 3 17 The transformation for µij and rij re°ects a case when the polarization is parallel to the rows of particles in the ^x direction and k is perpendicular to the array. After substitution into (16) we obtain fx;y(k;D) = 8>>>< >>>: SZ z }| { (2x2 ¡ y2) D3 p x2 + y2 5 + IZ z }| { ik(y2 ¡ 2x2) D2(x2 + y2)2 + RZ z }| { k2y2 D p x2 + y2 3 9>>>= >>>; eikDpx2+y2 . [17] For selected values of variable k and parameter D, the function fx;y exhibits symmetry around both ^x and ^y axes. Figure 1 shows that by choosing to center the reference frame within the in¯nite lattice we may take advantage of mirror symmetry about the ^x and ^y axes. We select this coordinate system to reduce the number of function calls to compute S. The computation of fx;y for any dipole in quadrant II, III, and IV thereby reduces to computation of the corresponding symmetric point of fx;y in quadrant I. We later introduce an algorithm to explicitly identify these symmetric points. To ensure a single center point exists we increment the grid size, L, by odd numbers. The retarded dipole sum, S, may be approximated using a ¯nite grid size L = 2N + 1, N = 1; 2; 3 : : :1, to obtain SL£L. To identify symmetric fx;y points which are redundant computations in the truncated SLxL summation we visualize individual summation terms for a given k and D as matrix elements. For N = 1 and N = 2, respectively, we obtain S3£3 = 1 4¼² X 2 4 f1;1 f0;1 f1;1 f1;0 0 f1;0 f1;1 f0;1 f1;1 3 5 and [18] S5£5 = 1 4¼² X 2 66664 f2;2 f1;2 f0;2 f1;2 f2;2 f2;1 f1;1 f0;1 f1;1 f2;1 f2;0 f1;0 0 f1;0 f2;0 f2;1 f1;1 f0;1 f1;1 f2;1 f2;2 f1;2 f0;2 f1;2 f2;2 3 77775 [19] where each element in the respective L £ L matrix of (18) and (19) are summed to obtain the truncated estimator, SL£L, of the dipole sum. Computation of SL£L is simpli¯ed by considering the symmetry, which reduces the number of function calls, (N + 1)2 ¡ 1, relative to the number of elements in the matrix, 2N + 1, by a fraction of (N + 1)2 ¡ 1 (2N + 1)2 [20] to provide a ¯rst estimate to the improved speed of the new algorithm. For large grid sizes, i.e. N ! 1, (20) approaches 1=4. Once redundant matrix elements have been identi¯ed they can be factored out and grouped. This can be done by writing the subscript variables x and y of fx;y into two M-dimensional vectors X and Y, where M = 8 PN i=1 i. This process may begin at the top center position of the outer shell in an L £ L matrix and move clockwise around the perimeter, repeating the process for each consecutive nested shell in the matrix to ¯ll the vectors. For example, in (18) the subscript variables of the top center function call f0;1 are transformed into initial elements of the respective vectors, X1 = 0, and Y1 = 1. Repeating this process for the entire matrix in (18) yields · X Y ¸ = · 0 1 1 1 0 1 1 1 1 1 0 1 1 1 0 1 ¸ [21] With the input variables of the 3£3 array expressed in vector form (21) it is clear that only three function calls are required to produce results for all eight elements in the array. For larger non-square or multidimensional arrays we introduce an algorithmic approach to group the redundant function calls. Such an algorithm would be helpful to eliminate redundant function calls for 3D clusters of NPs. An algebraic grouping algorithm is helpful since visual grouping in more complex situations becomes tedious. The truncated estimator, SL£L, can thus be represented in a linear algebraic format after grouping identical function calls in (21), viz., S3£3 = 1 4¼² (4f1;1 + 2(f1;0 + f0;1)) . [22] The resulting expression, (22), factors out all apparent redundant function calls of (18) thereby speeding computation. For larger arrays the redundant function calls can likewise be factored out to yield for N=2 and N=3 S5£5 = S3£3 + 1 4¼² (4(f2;1 + f2;2 + f1;2) + 2(f2;0 + f0;2)) [23] S7£7 = S5£5 + 1 4¼² (4(f3;1 + f3;2 + f3;3 + f1;3 + f2;3) + 2(f3;0 + f0;3)) [24] respectively. We group like terms in (22), (23), and (24) to reveal a trend which allows the function calls to be represented generally for any value of L: 4 www.pnas.org/cgi/doi/10.1073/pnas.0709640104 Footline Author 18 SL£L = S(L¡2)£(L¡2) + 1 4¼² ( 2(fN;0 + f0;N ) + 4 Ã XN j=1 fN;j + NX¡1 j=1 fj;N !) . [25] The new expression, (25), may be simpli¯ed further by combining the summation terms into a single summation. This is done by adding an additional term, fN;N , to the overall expression in (25) resulting in: SL£L = interior z }| { S(L¡2)£(L¡2) + 1 4¼² 8>>>>< >>>>: perimeter z }| { 2(fN;0 + f0;N ) + 4fN;N + 4 NX¡1 j=1 [fN;j + fj;N ] 9>>>>= >>>>; . [26] The leftmost term on the right hand side of (26) represents the interior points on any L £ L grid. The bracket terms in (26) represent the perimeter points that remain to be calculated. Now (26) can be rewritten as a double summation to account for the nesting of consecutive shells in the S(L¡2)£(L¡2) term that make up the interior points. SL£L = 1 4¼² XN m=1 2 64 axial z }| { 2(fm;0 + f0;m)+ diagonal z }| { 4fm;m 3 75 + off¡axial=diagonal(OAD) z }| { 1 ¼² XN m=2 mX¡1 n=1 [fm;n + fn;m] [27] where m and n are integers used to increment inner and outer summations. The double summation term in (27) includes all of the o®-axial/diagonal (OAD) elements. In (27) the axial elements are represented by fm;0 and f0;m and the diagonal elements are represented by fm;m. The remaining terms fm;n and fn;m represent the remaining elements that are o®-axial and o®-diagonal. With the new algorithm, an expression similar to (20) can be derived to estimate the overall improvement in computing speed. The number of function calls, ((N + 1)(N + 2) ¡ 1)=2, relative to the number of elements in the matrix is represented by the fraction (N + 1)(N + 2) ¡ 1 2(2N + 1)2 [28] where for large grid sizes (N ! 1), (28) approaches 1/8 for the new algorithm. The symbolic functions in (17) are now substituted into the new algorithm in (27), which consists of nested double-sum function calls, to allow algebraic and trigonometric simpli¯cation of function terms, using a reciprocal scaled wavelength de¯ned as A = kD = 2¼D ¸ . Its use ultimately permits single variable interpolation. This reduces (16) when N = 1, to the following solution [see Supporting Information: Appendix B]: S(A;D) = 1 4¼²D3 (2SI (A) + SII (A) + 4SIII (A)), [29] SI (A) = axial SZ z }| { P3(eiA)¡ axial IZ z }| { iAP2(eiA)¡ axial RZ z }| { A2 ln(1 ¡ cos(A) ¡ i sin(A)) [30] SII (A) = diag SZ z }| { P3(eip2A) p2 ¡ diag IZ z }| { iAP2(eip2A)¡ diag RZ z }| { p2A2 ln(1 ¡ cos(p2A) ¡ i sin(p2A)) [31] SIII (A) = OAD SZ z }| { Q3(eiA)¡ OAD IZ z }| { iAQ2(eiA)+ OAD RZ z }| { A2Q1(eiA) [32] where Pm is a polylogarithm de¯ned by Pm(Z) = 1X n=1 Zn nm = Z Z 0 Pm¡1(t) t dt [33] and Qz(x) is a double sum de¯ned as Qz(x) = 1X m=2 mX¡1 n=1 xpm2+n2 (m2 + n2)z=2 [34] and can be used to understand and simplify (32). The polylogarithms are special functions that can quickly be calculated within machine precision using the integral de¯nition, and for which analytical solutions exist for various cases. The complete solution, S(A;D), includes near- and far-¯eld contributions due to interactions between axial (SI (A)), diagonal (SII (A)), and o®- axial/o®-diagonal (SIII (A)) elements. The latter semi-analytical component requires further numerical reduction [see Supporting Information: Appendix C]. Footline Author PNAS Issue Date Volume Issue Number 5 19 Section 3: Discussion Divergences: existence and location. The new solution in (29) - (32) provides fundamental insight into the extinction spectra and increases computational speed and accuracy. The product S(A;D)D3 is a function of only the reciprocal scaled wavelength, A = kD. Singularities and divergences in exist when A = kD = 2¼n q x2f + y2 f [35] where the subscript f represents the ¯rst element along axial and diagonal rays of particles emanating from i at a constant angle, µif , and n is an integer multiple. For example, in the case of (30) when xf = 0 and yf = 1, (35) produces A = 2¼n which when substituted into (30) produces divergences when the natural log goes to zero for all n. Likewise, in the case of (31) when xf = 1 and yf = 1, (35) produces A = 2¼n=p2 which when substituted into (31) produces divergences for all n. The expression in (35) is equivalent to saying divergences exist when the wavelength is equal to the distance between the center particle i and the ¯rst element f, ¸ = rif n . [36] Real and imaginary solution components: axial, diagonal, and o®-axial/o®-diagonoal elements. The exact solution to the CDA for Maxwell's equations in (29) - (32) allows individual contributions from each set of contributing dipoles to be distinguished. Figure 2 shows real and imaginary components of contributions from axial (SI ,(30)), diagonal (SII ,(31)), and o®-axial o®-diagonal (SIII ,(32)) terms to the retarded dipole sum. The y-axis scale from -100 to +100 and the x-axis scale from 0 to 4.5 illustrate the behavior of the terms over an extended range of A=¼ to show the periodic behavior of the divergences. The SII and SIII components introduce a divergence at A = p2=¼. When added together, as in (29), the sum of these contributions produces a singularity in S at A = p2=¼ (as in Figure 5). Additional divergences observed to occur in SII and SIII components at integer multiples of A = p2=¼ do not produce such a singularity. Divergences occur in the SI term at integer multiples of A = 2¼ as indicated by (35). Divergences in the SIII term also appear at multiples of A = 2¼, although the source of these is not transparent in (32). Adding contributions from all points in the 2D array smooths the periodic behavior of the complete retarded dipole sum. Real and imaginary solution components: near-, induction- and far-¯eld contributions. Figure 3 shows real and imaginary components of contributions from static (r¡1, SZ), induction (r¡2,IZ), and far-¯eld (r¡3,RZ) radiative contributions to the retarded dipole sum that arise in SI , SII , and SIII . The magnitude of the near-¯eld contribution in the range shown is less than 5. In contrast to assumptions made by previous exact solutions which neglected e®ects of far-¯eld radiation, the magnitude of radiative contribution is clearly non-negligible: it dwarfs the intermediate and near-¯eld contributions. Except for the singularity at A = p2¼, the same divergences occur in Figure 3 as were seen in Figure 2. Grouping the contributions into static, intermediate, and far-¯eld terms simpli¯es illustration of the retarded dipole sum. Figure 4 shows real and imaginary components of the complete retarded dipole sum as de¯ned in (29). The large divergences that occur periodically as shown in (35) and (36) arise primarily due to the contributions from far-¯eld radiation. As grid size becomes very large the magnitude of these divergences becomes in¯nite. The behavior of the retarded dipole sum is a strong function of grid size near the divergence points. Larger grid sizes approximate behavior near these points more accurately. Convergence of the exact solution. To con¯rm the exact solution, Figure 5 compares real and imaginary parts of S for an in¯nite periodic lattice, given by (29), with values of S computed for a center particle in a ¯nite array, in the array size range L = [5; 40], by substituting (17) into (27). As L increases, peak values of the computed real part red shift and narrow until at L > 20 (400 nanoparticles), the ¯nite S begins to exhibit a narrow lineshape similar to the exact solution for an in¯nite array. For L = 100 the peakwidth of the ¯nite S is essentially equal to the in¯nite solution; only its amplitude remains smaller. The vertical displacement in Figure 5 between the ¯nite S sum and the in¯nite solution in (29) arises from constituent sinusoidal functions of the in¯nite solution, that are not present in (27) after performing the rectilinear transformation on the function call fi;j(k) in (16). The exact solution is computed using a smaller grid size (400) than in Figure 4 (1000), which increases the oscillation in the computed values. Inset in Figure 6 are the times (in seconds) required to compute values of the ¯nite retarded dipole sum at integer values of ¸ as L increases. Computation was performed using an AMD Turion dual-core CPU and computation time was recorded by onboard time analysis routines within the MATLAB TM software package. Computation time on the order of 60 hours is required to compute SLxL for an 80 x 80 array of nanoparticles for 600 distinct wavelengths. This compares with ca. 60 seconds required to compute the in¯nite solution. Convergence of computed ¯nite (i.e L < 1) retarded dipole sum values SLxL toward the in¯nite solution S(A;D) obtained by allowing N ! 1 in (27) is a semi-inductive method for validating the accuracy of the development in (49)-(68), which provides su±cient accuracy at reasonable computational expense. Computational methods such as FDTD and DDA which approximate (29) for one-dimensional arrays similary show that SL£L approaches the 1-D analog of (29) as N ! 1 [4]. However, the computational complexity of these methods is commonly greater than that of calculating SL£L using (27). To illustrate, computation times on the order of 10 hours per wavelength were reported to estimate the retarded dipole sum for an array of 1000 nanoparticles [3]. Validation of approximate solutions. The exact solution permits validating results from approximations to the retarded dipole sum obtained by ¯nite di®erence time domain (FDTD) or discrete dipole approaches. As examples, prior computation [38] suggested that divergences disappear as R < 30nm and a new divergence appears around A = ¼. This is clari¯ed by the exact solution, which shows divergences persist (albeit at small magnitudes) at all values of R, and the apparent divergence at A = p2¼ ' 4:44 disappears into a singularity as the grid size increases. Computations by DDA, T-matrix, FDTD require search algorithms to identify locations of the divergences, while the exact solution provides the precise location of divergences in the retarded dipole sum in (35). Beyond results from FDTD and DDA, the exact solution indicates there are an in¯nite number of 6 www.pnas.org/cgi/doi/10.1073/pnas.0709640104 Footline Author 20 divergences for an in¯nite array. Validation of (29) using DDA is precluded by inaccuracies reported to arise when interparticle spacing D is equal to incident electromagnetic wavelength ¸ [8]. Application of the exact solution. The exact solution in (29) - (32) shows the value of ¸ at which divergences occur is independent of other system parameters including size and dielectric constant of the dipolar inclusion, and dielectric constant of the surrounding medium. This speci¯city, and rapid calculation, make application of (29) preferable for multivariate char- acterization or optimization, relative to either evaluating large arrays like (27) using the ¯nite retarded dipole sum, or using approximations like FDTD. For example, S(A;D) may be employed to determine cross sections for absorption, scattering, and extinction due to dipole-dipole interactions of inclusions in a homogenous host, so-called extended Maxwell-Garnett composites in which inclusions are not small compared to wavelength described by e®ective-medium theory [39]. The exact solution to the CDA is being compared with results from T-matrix, CDA computations, and experiments to assess validity of CDA for describing arrays of NP with R > 50nm in which multipolar dynamic depolarization and radiative damping may decrease scattering. The exact solution enables rapid, precise estimations of the e®ect of changes in ¸ and D on optical cross sections, across large ranges in values of these parameters. With the exact solution in hand, values of k and D may be selected to produce divergences at a desired location (e.g. a semiconductor bandgap). Section 4: Conclusion An exact solution to the coupled dipole approximation of Maxwell's equations for electromagnetic interactions in periodic dipole arrays that include contributions due to far-¯eld radiation has been derived. This solution precisely identi¯es the source, location, and magnitude of periodic divergences in real and imaginary components of the retarded dipole sum which arise solely due to the wavenumber and grid size independent of particle size and dielectric constants of particle and environment. Contributions to the retarded dipole sum from near-¯eld, intermediate, and far-¯eld radiative interactions are distinguished. Far-¯eld radiative contributions are signi¯cant, outweighing near-¯eld and intermediate e®ects. Computed values of real and imaginary components of the retarded dipole sum in ¯nite arrays converge toward the exact solution as array size increases. Deriving the new solution produces a new method to optimize design of periodic subwavelength nanostructures for a wide range of applications. The improved computational speed of the exact solution allows for multivariate optimization of NP arrays. A single variable interpolation technique has been introduced which allows optimizing design of a sensor using a single retarded dipole array of su±cient grid size. The new solution allows straightforward implementation and design of NP arrays for a wide range of applications. It is a useful benchmark for analyzing performance of computational analysis. Phenomena observed in physical subwavelength nanostructures of increasing complexity can be interpreted with the solution. Divergences identi¯ed by the solution correspond to extraordinary spectral resonances in arrays of holes or particles. Applications to improve surface enhanced spectroscopies, sensors, integrated circuit waveguides, optical ¯lters, photovoltaics and antennae are envisioned. ACKNOWLEDGEMENTS. This work was supported in part by NSF CMMI-0909749. The authors would like to acknowledge valuable conversations with Phillip Blake and Wonmi Ahn. AUTHOR CONTRIBUTIONS. D.K.R. organized and directed the work and the preparation of the manuscript. B.J.T. performed computational derivation, analysis, and prepared drafts of the model section and ¯rst three paragraphs of the discussion section. B.M.H. performed computations and analysis to assist in revising the manuscript and prepared drafts of the remaining paragraphs of the discussion section.. Footline Author PNAS Issue Date Volume Issue Number 7 21 1. Meier M, Liao PF, Wokaun A (1985) Enhanced Fields on rough surfaces: dipolar interactions among particles of sizes exceeding the Rayleigh limit. J Opt Soc Am B 2:931-949. 2. Markel VA (1993) Coupled-dipole approach to scattering of light from a one-dimensional periodic dipole structure. J Mod Opt 40:2281-2291. 3. Zhao LL, Schatz GC, Kelly L (2003) The extinction spectra of silver nanoparticle arrays: inuence of array structure on plasmon resonance wavelength and widths. J Phys Chem B 107:7343-7350. 4. Zou S, Schatz GC, Janel N (2004) Silver nanoparticle array structures that produce remarkably narrow plasmon lineshapes. J Chem Phys 120:10871-10875. 5. Erin M et al. (2004) Investigating narrow plasmons in nanoparticle arrays fabricated using electron beam lithography. SPIE Proceedings 5513:22-29. 6. Zou S, Schatz GC (2004) Narrow plasmonic/photonic extinction and scattering line shapes for one and two dimensional silver nanoparticle arrays.J Phys Chem 121:12606-12612. 7. Fung KH, Chan CT (2007) Plasmonic modes in periodic metal nanoparticle chains: a direct dynamic eigenmode analysis. Opt. Lett. 32:973-975 8. Yin S, Luo X, Du C, Zhang Y, Deng Q (2008) The coupled electric ¯eld e®ects on localized surface plasmon resonance in nanoparticle arrays. J Appl Phys 104:024308. 9. Lamprecht B, Lechner RT, Ditlbacher H, Krenn JR, Leitner A, Aussenegg FR, Schider G (2000) Metal Nanoparticle Gratings: Inuence of Dipolar Particle Interaction on the Plasmon Resonance. Phys Rev Lett 84:4721-4725. 10. Haynes CL, Zhao LL, Van Duyne RP, Schatz GC, Gunnarsson L, Prikulis J, Kasemo B Kall M, McFarland Ad (2003) Nanoparticle Optics: The Importance of Radiative Dipole Coupling in Two-Dimensional Nanoparticle Arrays. J Phys Chem B 107:7337-7342. 11. Chamberlin DR, Sultana KA, Chow EK, Sigalas MM, Liu M, Grot AC, Fan S, Wang Z (2005) SERS and plasmon resonance of engineered nanoparticle arrays. Preceedings of SPIE 592708:1-8. 12. Felidj N, Aubard J, Levi G, Hohenau A, Krenn JR, Aussenegg FR, Laurent G (2005) Grating-induced plasmon mode in gold nanoparticle arrays. J Chem Phys 123:221103. 13. Hicks EM, Gunnarsson L, Zou S, Kasemo B, Kall M, Spears KG, Schatz GC, Van Duyne RP, Rinpnas.bstdzevicius T (2005) Controlling plasmon line shapes through diractive coupling in linear arrays of cylindrical nanoparticles fabricated by electron beam lithography. Nano Letters 5:1065-1070. 14. Chumanov G, Malynych SZ (2006) Nanoparticle Assemblies and Superstructures In Nanoparticle Assemblies, Superstructures pp. CRC press Boca Raton, FL. 15. McLellan E, Rindzevicius T, Kall M, Zou S, Spears K, Schatz G, Van Duyne R, Gunnarsson L (2005) Plasmonic and Di®ractive Coupling in 2D Arrays of Nanoparticles produced by Electron Beam Lithography. Mater Res Soc Symp Proc 951:E09-20. 16. Govyadinov AA, Markel VA (2008) From slow to superluminal propagation: Dispersive properties of surface plasmon polaritons in linear chains of metallic nanospheroids. Phys Rev B 78:035403. 17. Carron KT, Meier M, Wokaun A, Fluhr W (1986) Resonances of two-dimensional particle gratings in surface-enhanced Raman scattering. J Opt Soc Am B 3, 430-440. 18. Felidj N, Levi G, Krenn JR, Salerno M, Schider G, Lamprecht B, Leitner A, Aussenegg FR, Aubard J (2002) Controlling the optical response of regular arrays of gold particles for surface-enhanced Raman scattering. Phys Rev B 65:075419. 19. Genov DA, Shalaev VM, Wei A, Sarychev AK (2004) Resonant Field Enhancements from Metal Nanoparticle Arrays. Nano Letters 4, 153-158. 20. Li H, Rothberg L (2004) Colorimetric detection of DNA sequences based on electrostatic interactions with unmodied gold nanoparticles. PNAS 101:14036-14039. 21. Souza GR, Staquicini FI, Ozawa MG, Snyder EY, Sidman RL, Miller JH, Arap W, Pasqualini R, Christianson DR (2005) Networks of gold nanoparticles and bacteriophage as biological sensors and cell-targeting agents. PNAS 103:1215-1220. 22. Zhang Y, Previte MJ, Geddes CD, Aslan K (2008) Plasmonic engineering of singlet oxygen generation. PNAS 105:1798-1802. 23. Agrawal A, Wang GD, Wang MD, Nie S, Deo R (2008) Nanometer-scale mapping and single-molecule detection with color-coded nanoparticle probes. PNAS 105, 3298-3303. 24. Alu A, Engheta N (2006) Theory of linear chains of metamaterial/plasmonic particles as subdiraction optical nanotransmission lines. Phys Rev B 74:205436. 25. Driscoll T, Bason DN, Palit S, Cho SY, Jokerst NM, Smith DR (2007) Tuned permeability in terahertz split-ring resonators for devices and sensors. Applied Physics Letters 91:062511. 26. Shegai T, Dadosh T, Zhang Z, Xu H, Haran G, Li Z (2008) Managing light polarization via plasmon-molecule interactions within an asymmetric metal nanoparticle trimer. PNAS 105, 16448-16453. 27. Jackson JB, Halas NJ (2004) Surface-enhanced Raman scattering on tunable plasmonic nanoparticle substrates. PNAS 101:17930-17935. 28. GC, Schatz (2007) Using theory and computation to model nanoscale properties. PNAS 104:6885-6892. 29. Purcell EM, Pennpacker CR (1973) Scattering and absorption of light by nonspherical dielectric grains. Astrohphys J 186:705-714. 30. Yurkin MA, Hoekstra AG (2007) The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength. JQSRT 106:546-557. 31. Yurkin MA, Hoekstra AG (2007) The discrete dipole approximation: an overview and recent developments. JQSRT 106:558-559. 32. Merchiers O, Gonzalez F, Saiz JM, Moreno F (2007) Light scattering by an ensemble of interacting dipolar particles with both electric and magnetic polarizabilities. Phys Rev A 76, 043834. 33. Khlebstov B, Zharov B, Khlebstov N, Melnikov A (2006) Absorption and scattering of light by a dimer of metal nanospheres: comparison of dipole and multipole approaches. Nanotech 17, 1437-1445. 34. Biring S, Wang J, Wang Y, Wang H (2008) Light scattering from 2D arrays of monodispersed Ag-nanoparticles separated by tunable nano-gaps: spectral evolution and analytical analysis of plasmonic coupling. Optics Express 16:15312. 35. Klajn R, Grzybowski BA, Bishop KJ (2007) Light-controlled self-assembly of reversible and irreversible nanoparticle suprastructures. PNAS 104:10305-10309. 36. Ruppin R (2000) Evaluation of extended Maxwell-Garnett theories. Optics Communication 182:273-279. 37. Jackson J (1999) Classical Electrodynamics. (John Wiley Sons, INC.), Third edition, pp. 243-245. 38. Zou S, Schatz GC (2006) Theoretical studies of plasmon resonances in one-dimensional nanoparticle chains: narrow lineshapes with tunable widths. Nanotechnology 17:2813-2820. 39. Roper D.K, Ahn W, Taylor B, Y. D'Asen Enhanced spectral sensing by electromagnetic coupling with localized surface plasmons on subwavelength structures IEEE Sensors Accepted Feb. 17, 2009 8 www.pnas.org/cgi/doi/10.1073/pnas.0709640104 Footline Author 22 Figure 2.1. Submitted to Proceedings of the National Academy of Sciences of the United States of America Grid layout for a 5x5 square grid where rij represents the vector between the center point, i, and the point, j. Dotted and solid circles are shown to represent the nested grid concept used to obtain the analytical solution for the retarded dipole sum, S, for a square grid array. 2.2 Supporting PNAS Figures ,, ~ ~ J : e, ', , .,0 , '0 I 0 : .... ~ ... ~ 'I I I I I e e : ~ _____ .J o ..... r ·· . lJ :... ...} ····0 @ .. \ x~ lJ ~ ~ ._._.() .. _._IE) - - - ~ o o 0 000 Uniqlle calls: e el1ter point: Redllndant calls: '4._ 23 Figure 2.2. Submitted to Proceedings of the National Academy of Sciences of the United States of America The retarded dipole sum, S , may also be visualized as the three contributing equations (2.30,2.31,2.32) Figure 2.3. Submitted to Proceedings of the National Academy of Sciences of the United States of America The retarded dipole sum, S, contribution can be split into static (r1), intermediate (r2), and far- eld (r3) contributions for greater understanding 100 0 +-100 0 -100 0 100 o +-100 o -100 o - - 0.5 0.5 - --./ ,/ \ / \ / / - ./ \ ~- === '..:::: "- " I " 1 1.5 / / - - - - .L _ - -- - "- "- "' "- "- 1 1.5 1 \ / \ / \ / V / \ --- 2 2.5 Afrr I I I I / I I i - -- _ I ~ I I I / I I / I / ..- r \ I I I \ I \ 1 ..- -l "- -I' \ I \ \ / "- \ ~ 2 2.5 Afrr 1/ I ( / / // / / 1/ __ Re(SI) - - Re(SII) ___ Re(SIII) I __ lm(S,) I __ __ _ Im(SII) __ _ lm(SIII) ~ I 3 3.5 4 I I I / 1 I 1 I 1 I / I / / / f ) - I ' - --Re(r-1) "- I -2 I - - Re(r ) I -3 I - - - Re(r ) I '\ --lm(r-1) I I I \ -2 \ - - - -Im(r ) I ..- l - -3 - -v I - - - Im(r ) I 1 I \ "I I 'I I \ I I I I I I I 3 3.5 4 24 Figure 2.4. Submitted to Proceedings of the National Academy of Sciences of the United States of America The singularities in the retarded dipole sum, S, can be seen at the expected locations where A = p2 and A = 2 . An in nite number of singularities exist for an in nite array but due to the D3 factor in (2.29) singularities resulting in larger A are less signi cant Figure 2.5. Submitted to Proceedings of the National Academy of Sciences of the United States of America Computed values of the retarded dipole sum in (27) for nite arrays of size L converge to exact solution for retarded dipole sum in (2.29) for in nite array. As L increases, real and imaginary parts of (2.27) red and blue shift, respectively. For values of L > 20, both real and imaginary parts exhibit narrowed lineshapes. The exact solution is computed using a smaller grid size (400) than in Figure 4 (1000), which increases the oscillation in the computed values. Inset: Computation time for nite arrays increases exponentially with L 25 Supporting Information : Exact solution to the coupled dipole approximation of Maxwell's equations for electromagnetic interactions in periodic nanostructures Ben J. Taylor¤, Braden M. Harbiny, and D. Keith Roper¤ y ¤University of Utah Department of Chemical Engineering, and yUniversity of Arkansas Ralph E. Martin Department of Chemical Engineering Submitted to Proceedings of the National Academy of Sciences of the United States of America Appendix Section A: Electric field derivation The electric dipole field can be calculated by keeping only the 0th term (m = 0) in (6) which gives A(x) = ¹ 4¼ exp(ikr) r Z J(x0)d3x0. [37] Integration by parts and combining the relationship, i!½ = ~r ¢ J, from the continuity equation for charge, reduces (37) to A(x) = ¡ ¹ 4¼ exp(ikr) r i! Z x0½(x0)d3x0. [38] From the definition of the electric dipole moment, p = R x0½(x0)d3x0, we can further simplify (38) to get A(r; p) = ¡ ¹ 4¼ exp(ikr) r i!p. [39] The electric field can be found by substituing (39) into (5): E(r; p) = ic k ~r £ (~r £ A) [40] = c k !¹ 4¼ ½ ~r µ ikeikr r ¡ eikr r2 ¶ £ (^r £ p) + µ ikeikr r ¡ eikr r2 ¶ ~r £ (^r £ p) ¾ [41] = 1 4¼² ½µ 2 r3 ¡ 2ik r2 ¡ k2 r ¶ ^r £ (^r £ p) + µ ik r ¡ 1 r2 ¶ ~r £ (^r £ p) ¾ eikr [42] where ^r = r=r is a unit vector. The vector identity ~r £ (^r £ p) = ¡p=r ¡ (p ¢ ^r)^r=r can be used to further simplify (42) to obtain [1] E(r; p) = 1 4¼² ½µ 2 r3 ¡ 2ik r2 ¡ k2 r ¶ ^r £ (^r £ p) ¡ µ ik r2 ¡ 1 r3 ¶ (p + (p ¢ ^r)^r) ¾ eikr [43] which can be reduced by substituting in Lagrange's formula ^r £ (^r £ p) = (^r ¢ p)^r ¡ (^r ¢ ^r)p and factoring like terms to produce E(r; p) = 1 4¼² ½µ 1 r3 ¡ ik r2 ¶ [3(^r ¢ p)^r ¡ p] ¡ k2 r [(^r ¢ p)^r ¡ p] ¾ eikr [44] = 1 4¼² ½µ 1 r3 ¡ ik r2 ¶ [3(^r ¢ p)^r ¡ p] + k2 r [(^r £ p) £^r] ¾ eikr [45] which is the general form of the electric field for an oscillating dipole which corresponds to (6). Reserved for Publication Footnotes www.pnas.org/cgi/doi/10.1073/pnas.0709640104 PNAS Issue Date Volume Issue Number 1-5 26 In the ^x-^y plane the unit vector has cartesian coordinates ^r = cos(µ)^x + sin(µ)^y. If the polarization vector is parallel with ^x then µ is the angle between the two vectors ^r and p and the polarization can be defined as p = Px^x. Substituting the cartesian expression for the polarization and position vector we obtain a two component vector field defined by Ex = eikr 4¼² ½µ 1 r3 ¡ ik r2 ¶ (3 cos2(µ) ¡ 1) + k2 sin2(µ) r ¾ Px^x [46] Ey = eikr 4¼² ½µ 1 r3 ¡ ik r2 ¶ 3 cos(µ) sin(µ) + k2 cos(µ) sin(µ) r ¾ Px^y. [47] In the case of an array of particles symmetric about both ^x and ^y, the Ey, contributions sum to zero and the contribution from neighboring dipole fields can be written as X j6=i E(rij ; pj) = X j6=i eikrij 4¼² ½µ 1 r3 ij ¡ ik r2 ij ¶ (3 cos2(µij) ¡ 1) + k2 sin2(µij) rij ¾ pj [48] as seen in (13). Section B: Analytical Solution Derivation Beginning with the final result from the derivation for SL£L (27) SL£L = 1 4¼² XN m=1 2 64 axial z }| { 2(fm;0 + f0;m)+ diagonal z }| { 4fm;m 3 75 + off¡axial=diagonal(OAD) z }| { 1 ¼² XN m=2 mX¡1 n=1 [fm;n + fn;m] [49] we take the axial term when N = 1 Saxial = 1 4¼² 1X m=1 axial z }| { 2(fm;0 + f0;m) [50] and substitute in the expression for fx;y(k;D) in (50) and simplify to get Saxial(k;D) = 1 4¼² 1X m=1 eikDm 8>>>< >>>: Static (SZ) z }| { 2 µ (2) + (¡1) D3m3 ¶ + Intermediate (IZ) z }| { 2 µ (¡i2k) + (ik) D2m2 ¶ + radiation (RZ) z }| { 2 µ (0) + (k2) Dm ¶ 9>>>= >>>; [51] = 1 2¼² 1X m=1 eikDm ½ 1 D3m3 ¡ ik D2m2 + k2 Dm ¾ . [52] We introduce A = kD which yields Saxial(A;D) = 1 2¼²D3 ( 1X m=1 ¡ eiA¢m m3 ¡ iA 1X m=1 ¡ eiA¢m m2 + A2 1X m=1 ¡ eiA¢m m ) [53] = 1 2¼²D3 8>>< >>: SZ z }| { P3 ³ eiA ´ ¡ IZ z }| { iAP2 ³ eiA ´ ¡ RZ z }| { A2 ln ³ 1 ¡ eiA ´ 9>>= >>; [54] = 1 2¼²D3 fSI (A)g . [55] where P3 and P2 are polylogarithms defined in (33) and the RZ term contains the anlytical solution to P1, P1 = ¡ln(1 ¡ x). Likewise, the diagonal term is derived beginning with Sdiagonal = 1 4¼² 1X m=1 diagonal z }| { 4fm;m [56] we substitute (17) and simplify to get 2 www.pnas.org/cgi/doi/10.1073/pnas.0709640104 Footline Author 27 Sdiagonal(k;D) = 1 ¼² 1X m=1 eikDp2m2 8>>>< >>>: SZ z }| { (2m2 ¡ m2) D3p2m2 5 ¡ IZ z }| { ik(m2 ¡ 2m2) D2p2m2 4 + RZ z }| { k2m2 Dp2m2 3 9>>>= >>>; [57] = 1 ¼² 1X m=1 eikDp2m ½ p2 8D3m3 ¡ ik 4D2m2 + p2k2 4Dm ¾ . [58] Substituting in A = kD provides Sdiagonal(A;D) = 1 ¼²D3 1X m=1 8< : p2 8 ³ eiAp2 ´m m3 ¡ iA 4 ³ eiAp2 ´m m2 + p2A2 4 ³ eiAp2 ´m m 9= ; [59] = 1 4¼²D3 8>>>< >>>: SZ z }| { p2 2 P3 ³ eiAp2 ´ ¡ IZ z }| { iAP2 ³ eiAp2 ´ ¡ RZ z }| { p2A2 ln ³ 1 ¡ eiAp2 ´ 9>>>= >>>; [60] = 1 4¼²D3 fSII (A)g . [61] The remaining off-axial and off-diagonal terms can be derived beginning with SOAD = off¡axial=diagonal(OAD) z }| { 1 ¼² XN m=2 mX¡1 n=1 [fm;n + fn;m] [62] we substitute (17) to produce SOAD(k;D) = 1 ¼² 1X m=2 mX¡1 n=1 8>>>>< >>>>: SZ z }| { eikDpm2+n2 D3pm2 + n2 3 ¡ IZ z }| { ikeikDpm2+n2 D2pm2 + n2 2 + RZ z }| { k2eikDpm2+n2 Dpm2 + n2 9>>>>= >>>>; . [63] Substituting in A = kD provides SOAD(k;D) = 1 ¼²D3 1X m=2 mX¡1 n=1 8>>>>< >>>>: SZ z }| { (eiA)pm2+n2 pm2 + n2 3 ¡ IZ z }| { iA (eiA)pm2+n2 m2 + n2 + RZ z }| { A2 (eiA)pm2+n2 pm2 + n2 9>>>>= >>>>; [64] = 1 ¼²D3 8>< >: SZ z }| { Q3(eiA)¡ IZ z }| { iAQ2(eiA)+ RZ z }| { A2Q1(eiA) 9>= >; [65] = 1 ¼²D3 fSIII (A)g [66] where Qz(x) is a double sum defined as Qz(x) = 1X m=2 mX¡1 n=1 xpm2+n2 (m2 + n2)z=2 . [67] We combine results from (55, 61, 66) to obtain the final expression for S S(A;D) = 1 4¼²D3 (2SI (A) + SII (A) + 4SIII (A)) [68] which is the final form of the new analytical solution beginning in (29). Footline Author PNAS Issue Date Volume Issue Number 3 28 Section C: Evaluation of off­axial and off­diagonal terms in S In every continuous span of (32)1, i.e. A=[0,p2¼],[p2¼,2¼],. . . , a finite number of interpolation points may be used to fit (32) using piecewise cubic interpolation. Factoring out 1=D3 in (29) has allowed (32) to be represented as a function of A alone, which allows one dimensional interpolation as opposed to bicubic interpolation in variables k and D. To optimize the spacing between the finite number of points, a Monte Carlo simulation was used. The Monte Carlo simulation was chosen due to its simplicity to code with respect to other optimization routines, such as the genetic algorithm or game theory. The Monte Carlo simulation was able to reproduce identical results when sufficiently large simulations were run. Table 1 identifies interpolation points obtained across the interval A = [0; 2¼]. Over this interval, a value of R2¸ 0.99 was obtained using piecewise cubic interpolation. The number and/or spacing of interpolation points in the Monte Carlo simulation could be increased to obtain higher accuracy. New interpolation points are required to extend computation of the exact solution to values A=¼ > 4:2 (e.g. into the infrared region). Figure 1 shows the interpolation points together with the behavior of real and imaginary components of (32) - off-axis and off-diagonal contributions to the retarded dipole sum - plotted versus A=¼ up to a value just larger than 2. The simulation shown used 10002 NPs. Oscillations observed in Figure 2 dampen out with larger grid sizes. The divergence at A=¼ = p2 becomes a singularity when contributions from (30) and (31) are included to obtain the complete retarded dipole sum (shown in Figure 5), whereas the divergence at A=¼ = 2 remains a divergence. The y-axis scale from -200 to +200 and the x-axis scale from 0 to 2.1 were selected to show the concentration of the interpolation points near the divergences. It is useful to note that the retarded dipole sum need not be recalculated when another variable besides k and D is adusted during multivari-able optimization on any mathematical derivative of the retarded dipole, like polarization, optical cross section, or electric field in (15). As an example, when examining effects of varying particle diameter or dielectric constant on optical cross section, only a single set of calculated or interpolated values of the retarded dipole sum for a given matrix of k and D is required. 1 Equations (1) ­( 36) refer to equations within the main document. 4 www.pnas.org/cgi/doi/10.1073/pnas.0709640104 Footline Author 29 Figure 2.6. Submitted to Proceedings of the National Academy of Sciences of the United States of America The behavior of the SIII double sum is shown with the interpolation points. This simulation used 10002 NPs which produces oscillation 200 -U)- Qc:::) 0 -200 -U)- 0 I 200 I 0 0.5 1 A/n (kDln) 1.5 2 30 2.3 Design Optimization To illustrate the usefulness of the new analytical solution, an Ag NP array sensor application will be optimized for sensitivity. The NP array sensor is most sensitive when sharp spectral peaks exist where slight shifts in local refractive index cause dramatic changes in the peak wavelength, as seen in Figure 2.7. With the new model, the existence of extraordinary peaks can be predicted with the knowledge of the location of singularities prior to computation. Knowing the precise location of the singularities limits the time needed to locate them as was previously done using brute force search methods [22]. The sensitivity, which is the maximum change in extinction e ciency, can be seen in Figure 2.7 for various interparticle spacings. For an in nite spacing, or single particle spectra, the sensitivity is the lowest because a change in refractive index will have a minimal change in extinction e ciency. However, for the blue curve which is an interparticle spacing of 480nm the peak is very sharp and therefore a small change in refractive index will lead to a large change in extinction e ciency. So a spacing of 480nm would be preferred for Table 1: Points of interpolation for computing real and imaginary components of (32) using one dimensional piecewise cubic interpolation in A. Interpolation points were determined by using a Monte Carlo simulation to determine optimal spacing. The values reported are the result of a 100,0002 NP simulation. 1 A=kD Re(SIII) Im(SIII) 0.00 0.3150 0 0.50 -0.1267 0.3467 1.40 -0.4922 -0.4190 2.10 0.2066 -0.8707 2.75 1.4565 -0.5822 3.55 2.0348 2.7089 4.00 -0.5870 5.9506 4.25 -6.0960 8.6120 4.35 -11.3767 9.8267 4.40 -16.9787 10.6193 J27r -oo( Grid --+ (0) 10.6193 R"2 0.9981 0.9997 J27r -00 0 4.463 -23.09 -10.67 5.423 -13.85 -10.95 5.583 -4.34 -10.86 5.963 27.45 -2.93 31 Figure 2.7. (Table and caption reproduced with permission from IEEE) Near- eld plasmon features (<400nm) with far- eld extinction peaks (>450 nm) for 50 nm radius Ag NP array embedded in air. The di erent colored lines represent di erent interparticle spacings and the black line is spectra for an in nite spacing. In nite spacing can also be thought of as a single particle spectra. sensitivity for silver NP radius of 50nm in a medium with a local refractive index of 1. Figure 2.8 shows a continuation of an Ag sensor design optimization similar to Figure 2.7 but with higher step-resolution combined with the optimization of particle size. The colormap is linked to sensitivity which is proportional to the maximum slope of the spectral response and the desired optimization tness parameter. With the new algorithm, Figure 2.8 was created in less than a day and quickly identi es ideal spacing and particle size for the desired target wavelength. These results are signi cant because they allow for rapid NP sensor design optimization and therefore maximum sensitivity. In addition to optimizing the NP size and interparticle spacing, the new model can optimize refractive index and even NP material allowing new, unstudied materials to be noticed as potential replacements for Ag and Au NP arrays. The optimization results from the model can be used to match the manufacturing process and minimize the number of experiments needed to achieve the desired output with the aid of the advanced analysis techniques discussed in x4.3-4.5. 32 Figure 2.8. (Table and caption reproduced with permission from IEEE) Sensitivity of AU NP extinction e ciency (EE) with respect to change in local refractive index, o, where sensitivity is inversly proportional to the peak full width half max (FWHM) and de ned as the maximum extinction e ciency delta. A near eld subset from the inset plot shows e ects of changing particle radius, r, and interparticle spacing d. The colorscale shows sensitivity values up to 188. Dashed red lines show where relative separation distance, d=2r, is constant as r and d are varied. The insert plot shows extraordinary far- eld e ects predicted by the new algorithm. Dashed lines in the insert plot show discontinuities predicted using the new algorithm with adjacent high sensitivity. This plot allows the search space to be viewed for the optimization problem, and demonstrates the predictability of the singularities. When optimizing this problem in practice for higher dimensions the results are not visualized, but optimized using a metaheuristic global search algorithm. , 6S E 2B S 0 12 , 2 0.1 CHAPTER 3 FABRICATION 3.1 Introduction Nanoparticle (NP) fabrication is the nal step of producing sensors. After electromag- netic theory has been used to optimize the theoretical NP size, spacing, and material, a method must be developed to create such a device. Using advanced electroless deposition (EL), Ahn et al. [2] have developed a low cost technique capable of producing random arrays of Au NPs with varying size and spacings. 3.2 Heating Au Island Films Forms Nanoparticles Nanoparticles (NP) have been fabricated on Silicon (Si) substrate using various methods including ion implantation and surface etching [16], laser ablation [9], inductively coupled thermal plasmas [11], thermodynamic self organisation [8], and thermal evaporation of thin Au lm followed by annealing [6]. Island lms and thin lms have been formed using electrodeposition, electroless de- position (EL) [15], sputtering [10], and vacuum evaporation [5]. A new NP fabrication technique has been developed which uses EL to create island lms and then thermally annealed the island lms to create NPs of controllable size distributions. The new EL process has an advantage over other methods for being able to rapidly coat fragile, 3D, or internal surfaces at ambient conditions without requiring expensive specialized equipment or conductive surfaces in the case of electrodepostion. The EL process is also relatively inexpensive for mass fabrication and scales well making it a promising economical solution for NP array applications in the market. Ahn et al. outlined the EL plating method [2] which involved three steps. To summarize: rst, sensitization of quartz slides (GE 124 fused quartz, Chemglass Inc., Vineland, NJ) by immersion in SnCl2 and tri uoroacetic acid followed by rinsing with distilled deonized, degassed (DDD) DDD-H2O and drying with N2 gas; second, silver (Ag) deposition was accomplished using AgNO3 followed by rinsing with DDD-H2O and drying with N2 gas; 34 nally, immersion of Ag-derivatized quartz slides into the gold (Au) solution Na3Au(SO3)2 for the desired deposition time and rinsing and drying to complete. For a complete recipe detailing concentrations, times, and additional solutions and steps needed for the new EL method along with more experimental results refer to the full paper [2]. Once the Au island lms have been created, the transformation from island lms to NPs can be seen in Fig 3.1. The rst row shows EL island lms for (A) 0.33 min (B) 1 min, and (C) 4 min Na3Au(SO3)2 deposition times. The second row are the same samples after annealing for three hours at 250 C. The nal step involves heating the annealed samples to 800 C for 20 min to create spherical NPs, as shown in the bottom row in Fig 3.1. The transmission and extinction spectra were collected by Ahn [2] who assisted in the collection of SEM images as well in Fig 3.1. Notice as the samples become more spherical Figure 3.1. (Figure and caption reproduced with permission from Langmuir [23]) SEM images of fuzzed quartz slides after immersion ( rst row) in Au plating solution for (A) 0.33 min, (B) 1 min, and (C) 4 min; after respective immersed slides are annealed (second row, images D, E, and F) for 3 h at 250 C; and after respective immersed, annealed slides are heated (third row, images G, H, and I) for 20 min at 800 C. Inset are transmission UV-vis spectra corresponding to each image. Percent transmission is shown in images A-C, E, and F. Extinction in arbitrary units (au) is shown in images D, G-I. All images are the same magni cation; see scalebar on the left. 35 the inset UV-vis spectra show a sharpening resonance peak that is anticipated with NP array theoretical modeling. Longer deposition times were also studied as well with Fig 3.2 showing (A) 8 min, (B) 30 min, and (C) 60 min Na3Au(SO3)2 deposition times. The sharper contrast between Fig 3.1 and Fig 3.2 is a result of the decreased magni cation and also a higher quality low vacuum SEM used at Brigham Young University. The weblike structures seen in Fig 3.2 are NaCl-like crystal structures determined by X-Ray Di raction (XRD) that were later removed after prolonged rinsing with DDD-H2O. Manual masking had to be performed to allow data analysis due to the image processing interference created by the crystal structures. The data and statistics from the NP samples in Fig 3.2 and Fig 3.1 are discussed in the following data analysis chapter. 3.3 Multiple Experiments: Drop Method To reduce the amount of Na3Au(SO3)2 needed to plate a thin lm and to increase the number of deposition time experiments, an EL drop method was developed with the help of Ahn. After the proposed procedures of Ag and Au deposition on quartz slides, a drop of 1-5 L Na3Au(SO3)2 solution was placed on a slide for various amounts of time. Images (see Fig 3.3) and spectra resulting from drop volumes in this range were similar. To maintain the necessary N2 to prevent Na3Au(SO3)2 from oxidizing, a polypropylene container was used during the deposition period to enclose the sample. Afterward, the Au drop-plated slide was removed and washed with DDD-H2O and dried with N2 to create the samples. The NP diameters size distributions were similar throughout the drop; however, exces- sive Au deposition was noticed around the droplet perimeter. The excessive Au deposits can Figure 3.2. (Figure and caption reproduced with permission from Langmuir [23]) SEM images of fused quartz slides annealed (250 C; 3h) and heated (800 C; 20 min) after immersion in Au plating solution for (a) 8 min, (b) 30 min, and (C) 60 min. Weblike strands are NaCl-lie crystal structure removed after rinsing with DDD-H2O which were identi ed using X-Ray Di raction (XRD). 36 Figure 3.3. (Figure and caption reproduced with permission from Langmuir [23]) SEM images after annealing (250 C, 3 h) fused quartz slides exposed to microliter droplets of EL plating solution exposed for (A) 0.33 min, (B) 1 min, and (C) 2 min. be explained by the increased local concentration of Na3Au(SO3)2 as the solvent evaporated at the droplet boundary during the deposition time. A low vacuum SEM image can be seen in Fig 3.4 along with an average intensity plot created from Fig 3.4 in Fig 3.5 showing the relative changes in SEM intensity as a function of droplet radius. The charging e ects observed near the edge caused varying intensity in the eld of view (FOV); therefore, the trending intensity seen towards the edges of the image should not be interpreted as changes in Au concentration. Figure 3.4. Low vacuum SEM image showing increased Au deposits at the drop boundary as the drop evaporated from N2 convection during the deposition time 37 5 10 15 20 25 0 20 40 60 80 100 μm Relative Intensity Figure 3.5. Relative pixel intensity from SEM image of droplet edge shows increased Au deposits near the boundary and the change in Au coverage uniformity moving towards the center of the particle. The vertical dashed line represents the best approximation of the outside droplet boundary. Relative slopes at the beginning and end of the plot should not be taken literally since charging may be e ecting perimeter intensity in the Field Of View (FOV). Figure copyright is owned by author 3.4 Experimental Results Compared to Theoretical As shown in Fig 3.6, the theoretical results produce a higher extinction e ciency due to the ordered grid. A randomly oriented array of NPs cancels out far- eld interactions and therefore would dampen the extinction. The spread between the experimental results is most likely caused by a di erence in particle size distributions between the samples which was an artifact from varying spacing. Experimental results re ect relative shifts but not overall t. To more accurately test the validity of the theoretical solution, an ordered array with xed NP size would be required. This could be achieved using lithography, and is the next logical step in this research because theoretically, ordered arrays produce orders of magnitude sensitivity improvements. 3.5 Conclusion Despite the novelty of the new EL method and its low cost appeal, dramatic sensitivities can only be realized with an evenly spaced array and equal particle sizes. If spacing and particle size are not equal, then the equal polarization assumption (2.15) used in the derivation is not valid. The theoretical work in Section 2.2 shows that far- eld e ects are additive and therefore exact spacing is preferred to achieve maximum sensitivity. 38 Figure 3.6. Theoretical results deviate from the experimental results for several reasons. The experimental results are created with randomly oriented arrays with a particle size distribution. This violates the equal polarization assumption used in the theoretical derivation (2.15). In order to perform a more rigorous appraisal of the theoretical results, a lithography method of fabrication would be needed to ensure equal spacing and particle size. Figure copyright is owned by author 39 Future work will most likely involve some form of lithography similar to what is used for semiconductor manufacturing. E-beam lithography would be preferred for initial research but mass production of sensors would require a more economical solution such as pho- tolithography. Photolithography methods would require an expensive reticle mask, but would save money long-term. CHAPTER 4 DATA ACQUISITION AND ANALYSIS 4.1 Data Aquistion The NP arrays were analyzed using Atomic Force Microscopy (AFM) and a Scanning Electron Microscope (SEM). The SEM images were obtained using a Philips XL30 ESEM FEG (FEI, Hillsboro, Oregon) using a backscattered electron detector mode operating with a 15 kV electron accelerating voltage. The low vacuum setting was used and the pressure was varied independently for each image based on charging and desired resolution over a low vacuum range of 0.4-1.2 Torr. Acquisition preprocessing software automatically stretches to accommodate low- and high-intensity in an image making intensity comparisons between images susceptible to error. The low vacuum setting helped reduce the charging e ects that were observed on high vacuum SEM images taken previously which improved the image contrast and accuracy in particle size. 4.2 Watershed-Transform and Shape Analysis SEM particle analysis is typically performed using threshold algorithms for isolating particles. Fig 4.1 shows that this method can produce inaccurate particle size and shape measurements when adjacent particles are close and when poor image resolution and quality result in clumping. A watershed transformation [13] was used to segment the clumps and provide more accurate size and shape characterization. Brie y, the watershed transform segments the image by relying on the image topology. The original image is transformed into the complement, making local maxima into local minima and vice versa. The watershed transform then begins to ll the now local minima volume, hence the name. While lling, if the local minima volume contacts an adjacent local minima volume, a division is made and the particle becomes segmented. Two algorithms have been developed to complement the adapted watershed transfor- mation. The rst algorithm strips out particles that are touching the border of an image to exclude incomplete particles from the measurements of particle size, shape, and number 41 Figure 4.1. (Figure and caption reproduced with permission from Langmuir [23]) The wa- tershed method distinguishes particle characteristics more realistically than simple threshold methods when particles are grouped in clusters. (A) Original SEM image of Au-coating on quartz slide. (B) Threshold method of segmentation results in particle "fusion" with size distribution skewed by particles at the border. (C) Watershed method with border removal quantitatively distinguishes particle size, shape, and number density. density, as seen in Fig 4.1:C. The sodium crystal-like structures in Fig 3.2 that disappeared when samples were rinsed copiously with DDD-H2O after EL were manually corrected after analysis by the watershed method to di erentiate these structures from Au NPs. Crystals selected as particles were manually deselected and missed particles were manually reselected using the open source computer program GIMP (v. 2.2.11; www.gimp.org). The image processing was accomplished using image processing tools (Image Processing Toolbox 6.0) in MATLAB (v. 7.4; MathWorks). The second algorithm has been developed to analyze particle shape as shown in Fig 4.2. The algorithm selects individual particles and then optimizes the theta, , of rotation to minimize the di erence in particle shape between the particles. The resulting two- dimensional shapes are then added on top of each other to produce a 3D histogram (Fig 4.3). The resulting 3D histogram can then be used to estimate the probability of nding a particle with a particular shape and easily determine the average sphericity of the Au NPs and other shape-related metrics. 4.3 Applying Watershed to SEM Images Fig 4.3 shows particle size distributions determined by watershed analysis from SEM images in Fig 3.1 (row 3) and Fig 3.2 for heated Au island lms immersed at di erent times. Fig 4.3:A and Fig 4.3:B show that as immersion time increases from 8 to 30 to 60 min, average particle diameters produced by successive annealing and heating changed from 169.0 2.5 to 266 22 to 97.0 9.5 nm, particle density changed from 7.2 to 4.3 to 12 counts per cm2 x 10􀀀8, and Full Width Half Max (FWHM) changed from 174 to 60.7 to 113 nm. 42 Figure 4.2. Each particle was centered by nding the individual centers of mass and aligned by using the golden search method to maximize the moment about a particular axis. The resulting 3D histogram allows for mean shape determination and shape con dence plots analysis. A shape con dence plot would be an outline shape that contained a quantile of the particle population as shown by the two black contour traces. Figure copyright is owned by author 43 Figure 4.3. (Figure and caption reproduced with permission from Langmuir [23]) Mean diameter and density of particles obtained after annealing at 250 C (3 h) and heating at 800 C (20 min) shift to larger scales as E.L. deposition time increases from (A) 0.33 to 1 to 4 min. (B) Mean diameter and particle density for annealed and heated lms deposited at 8, 30, and 60 min. ~ '"E 8 -:-i U1 ~ 6 C ::l 0 ~u 4 >- :~ U1 c 2 (]) ~ c 0 10' i 1.5 -:-i El c 1 :::J o U . : .. ··0.33 min 1 min - 4min (Al Circle equivale nt diameter (nm) : " ...:.'.:, . , . : " 102 -- Bmin -- 30 min ·· .. ·· 60 min (8) Circle equivalent diameter (nm) 44 Increasing deposition time enlarges Au nanoparticles and broadens the size distribution. EL samples are skewed to smaller particle values by a large number of ne particles identi ed by watershed analysis in the 30 and 60 min samples. 4.4 Particle Diameters, SEM vs AFM Table 4.1 shows particle diameters obtained by watershed analysis of SEM images and AFM. The nonlinear behavior of deposition time versus particle size may have been related to the surface area over volume fraction. With a constant deposition ux relative to the surface area, the behavior of deposition would be inversely proportional to the diameter. Therefore, for smaller diameters a faster particle growth rate is expected. AFM diameters exceed SEM diameters in each case except for 4-min deposition, due in part to convolution between probe tip and particle shape. The tip convolution is due to the angle of the tip that prevents the tip of the probe from immediately contacting the substrate once it leaves the particle, therefore leading to an exaggerated particle size. A second computer algorithm has been developed, arbitrarily called "Peak Height Method (PHM)," to determine particle heights from the local minima of particle gradient obtained from inverted AFM images. Table 4.1 shows the mean value of AFM particle heights obtained by this method compared to SEM particle diameters obtained by watershed analysis. As expected, the PHM more closely matches the SEM values than using watershed on AFM data due to the elimination of AFM tip convolution. Table 4.1. (Table and caption reproduced with permission from Langmuir [23] ) The AFM (PHM) ts the SEM (WM) data better than the AFM (WM), which overpredicts the particle size, due to tip convolution. time (min) particle diameter (nm) SEM(WM) AFM(WM) AFM(PHM) 1 67.7 1:1 111.9 3:9 42.3 1:1 4 119.6 4:2 114.1 4:4 49.98 0:82 8 169.4 2:5 241.3 4:6 116.7 2:6 30 266 22 764 80 183 15 60 97 13 580 42 136.9 8:7 45 Figure 4.4. (Figure and caption reproduced with permission from Langmuir [23]) Particle densities estimated by watershed analysis of AFM ( lled circles) and SEM (open circles) images for immersion times at 0.33, 1, 4, 8, 30, and 60 min. The 0.33 min sample was damaged prior to AFM imaging. 4.5 Particle Densities, SEM vs AFM Fig 4.4 compares particle densities obtained by watershed analysis of AFM and SEM images. To nd the particle densities on the Si surface, software counted the number of particles and divided by the total analyzed surface area. Compared with watershed analysis of SEM, AFM images returned smaller particle counts per square micrometer by a factor of 5 at 1 min immersion and by a factor of 14 at 60 min immersion as shown in Fig 4.4. Watershed analysis of SEM images appears to identify small ( 5 nm) particles more often than AFM images. Measured particle densities are most similar at 8 min, where SEM and AFM images each show monomodal distributions of 169 nm particles. 4.6 Conclusion For proper model validation and troubleshooting, accurate image processing methods must be used to analyze SEM and AFM images. The watershed method and PHM provided a more accurate analysis of our NP samples. In addition to developing the algorithms, the comparison between the two most common acquisition methods, SEM and AFM, showed shortcomings of AFM. Result con dence was increased using automated image processing algorithms for SEM which allowed for large sample statistical analysis (over 1000 particles). CHAPTER 5 CONCLUSION 5.1 Concept to Reality This work has covered a wide variety of topics including NP array modeling, fabrication, and model veri cation through NP analysis. Some of the more noteworthy points of this work as well as novel discoveries have been outlined in bulleted form below. A new analytical solution to Maxwell's equation for in nite two-dimensional arrays for retarded dipole interactions between spherical nanoparticles was derived. This new solution is valid for nonmagnetic spheres that are small enough (< 100 nm) that multipole interactions can be ignored and large enough that quantum e ects can also be ignored (2.29). The new Monte Carlo interpolation method allowed for a decreased number of function evaluations. Using this new method to replace a xed step size in independent variables allows for larger array simulations. Using the Monte Carlo interpolation technique that was developed, 100; 000; 000; 000 NP 2D array simulations were demon- strated in this work and multivariable design and optimization was improved by increasing the algorithm e ciency (2.C). This new variable step size method is only accurate when the exact knowledge of the singularities are known in order to determine the proper interpolation ranges, i.e. A = [0p2 ]; [p2 2 ], etc. An expression (2.35) capable of de ning all singularities within the new analytical solution (2.29) was also obtained. The knowledge of these singularity locations and behavior o ers new fundamental insight into NP array sensitivity enhancement. In addition, singularities which were previously searched for numerically can be found analytically without using costly direct search methods. Two-dimensional NP array interparticle spacing, NP diameter (valid for NP diameters < 100nm but su ciently large to ignore quantum e ects), surrounding medium, and NP material can be tuned using the speed of the new algorithm to optimize sensor design for target wavelengths (2.8). 47 In coordination with Ahn et al. [2], a new fabrication technique was developed which involves electroless deposition (EL) followed by thermal annealing to produce random gold nanoparticle arrays (3.1). A comparison between AFM and SEM postprocessing particle size and spacial particle density analysis was conducted. A solution was also developed to counter AFM tip- convolution using the Peak Height algorithm (Table 4.1) to produce more closely matched datasets between AFM and SEM tools. A particle shape analysis algorithm was invented allowing randomly oriented nanopar- ticles to be analyzed for shape. Particles were rotated and aligned to produce a 3D histogram of stacked particle shapes. From the 3D histogram, particle shapes could be produced that represented any desired percentile of the population allowing particle shape con dence plots. The watershed transform was adapted to show the particle selection bene ts for touching particles compared to simpler threshold alternatives which tend to exaggerate particle size due to clumping (4.2). Combining the new analytical solution Monte Carlo technique, EL annealing method, and analysis methods provides a total concept to reality procedure for producing optimized NP designs. NP composition, spacing, size, wavelength, and refractive index can be opti- mized for every scenario within the model validity range. Once the parameters have been determined for an Au array, EL annealing experimentation (Fig 4.1) can be performed. The watershed transform can then be used on SEM images to acquire particle size distributions, count densities, and shape characteristics. Automated SEMs and scatterometers used in the semiconductor industry can be used likewise to achieve real-time process control and array analysis in a manufacturing environment. 5.2 Future Work The analytical solution is valid for a nite range of spherical NP diameters for 2D arrays due to multipole and quantum e ects. Future work to improve the current model validity range by eliminating assumptions or maximize its use have been listed below. A general expression valid for spherical NP diameter >100 nm could be derived by replacing (2.13) with the higher multipole expression for a retarded electric eld and redoing the remainder of the derivation. This new expression would then be valid for all nonmagnetic spherical 2D arrays where NP diameter size was large enough to 48 ignore quantum e ects. Au and Ag are common NP array materials for sensor applications. Using the speed improvements of the new analytical solution, NP material could be optimized to search for a new theoretical material with superior optical properties. Once an ideal material was identi ed, the NP 2D array fabrication technique would then have to be developed as well. The redundancy elimination procedure for the 2D in nite NP array could be expanded to include that of a 3D in nite NP array. The derivation would be modi ed where the 3x3 matrix in (2.18) would be replaced by a 3x3x3 matrix and (2.19) by 5x5x5 before continuing on with the derivation. The dilogarithm, P2(Z) = P 1n=1 Zn n2 , o ers a signi cant roadblock to the continued simpli cation of the analytical solution. Progress towards obtaining a dilogarithm analytical solution or symbolic regression approximation could o er additional insight and improvement. Adding the substrate electric eld contributions to (2.13) for a more accurate solution using the principle of superposition. An other alternative is using a surrounding medium with an identical refractive index as the substrate to eliminate substrate electromagnetic interference. Despite the bene ts of the current model, future work is necessary to expand the model validity range and usefulness. With the existing model, direct search methods could be used to inversely solve the necessary material optical properties for maximum sensitivity. The desired optical properties could be compared to a database of known material, or using physical chemistry nonexistent compounds/alloys could be suggested. REFERENCES [1] Wikipedia, the free encyclopedia, Nov 2009. [2] Wonmi Ahn, Benjamin Taylor, Analia G. Dall' Asen, and D. Keith Roper. Electroless gold island thin lms: Photoluminescence and thermal transformation to nanoparticle ensembles. J Opt Soc Am B, 2:931{949, 2008. [3] Hitachi vendor Andy. Hitachi vendor training. In IMFT Hitachi vendor training, July 2009. [4] S. T. Bramwell, S. R. Giblin, S. Calder, R. Aldus, D. Prabhakaran, and T. Fennell. Measurement of the charge and current of magnetic monopoles in spin ice. Nature, 461:7266, 2009. [5] Sockalingum G D, Beljebbar A, Morjani H, Angiboust J F, and Manfait M. Title not added. Biospectrosc., 4:S71{S78, 1998. [6] Spadavecchia J, Prete P, Lovergine N, Tapfer L, and Rella R. Title not added. Phys. Chem. B, 109:17347{17349, 2005. [7] J Jackson. Classical Electrodynamics, pages 243{245. John Wiley Sons, INC., third edition, 1999. [8] Beomseok K, Tripp S, and Wei A. Self-organization of large gold nanoparticle arrays. J. Am. Chem. Soc., 123:7955{7956, 2001. [9] Hata K, Yoshida S, Fujita M, Yasuda S, Makimura T, Murakami K, and Shigekawa H. Self-assembled monolayer as a template to deposit silicon nanoparticles fabricated by laser ablation. J. Phys. Chem. B, 105:10842{10846, 2001. [10] Khriachtchev L, Heikkila L, and Kuusela T. Title not added. Appl. Phys. Lett., 78:1994{1996, 2001. [11] Shigeta M and Watanabe T. Growth mechanism of silicon-based functional nanopar- ticles fabricated by inductively coupled thermal plasmas. Appl. Phys., 40:2407{2419, 2007. [12] James Clerk Maxwell and Joseph John Thomson. A Treatise on Electricity and Magnetism, volume 2. Clarendon press, 3 edition, 1904. [13] F Meyer. Topographic distance and watershed lines. Sig. Proc., 38:113{125, 1994. [14] Michael I. Mishchenkoa, Gorden Videenb, Victor A. Babenkoc, Nikolai G. Khlebtsovd, and Thomas Wriedte. T-matrix theory of electromagnetic scattering by particles and its applications: a comprehensive reference database. Journal of Quantitative 50 Spectroscopy Radiative Transfer, 88:357{406, 2004. [15] Mallory G O and Hajdu J B. Electroless Plating: Fundamentals and Applications, chapter 1. American Electroplaters and Surface Finishers Society: Orlando, FL, 1990. [16] Nakano S, Ishizuka M, Wang Q, and Ogiso H. Nanoparticle fabrication for micro total analysis system by using ion implantation and surface etching. Surface and Coatings Technology, 187:167{171, 2004. [17] M N O Sadiku. Elements of Electromagnetics, page 386. Oxford University Press, fourth edition, 2007. [18] Tapan K. Sarkar, Robert Mailloux, Arthur A. Oliner, Magdalena Salazar-Palma, and Dipak L. Sengupta. History of Wireless, page 22. Wiley-IEEE Press, rst edition, 2006. [19] A. Ta ove, K. R. Umashankar, B. Beker, F. A. Harfoush, and K. S. Yee. Detailed fdtd analysis of electromagnetic elds penetrating narrow slots and lapped joints in thick conducting screens. Antennas and Propagation, IEEE Transactions, 32:247{257, 1988. [20] Ben Taylor, Braden Harbin, and Keith Roper. Exact solution to the coupled dipole approximation of maxwell's equation for electromagnetic interactions in periodic nanos- tructures. PNAS, Submitted 2011. [21] Maxim A. Yurkin and Alfons G. Hoekstra. The discrete dipole approximation: an overview and recent developments. J.Quant.Spectrosc.Radiat.Transf., 106:558{589, 2007. [22] Schatz GC Zou S. Theoretical studies of plasmon resonances in one-dimensional nanoparticle chains: narrow lineshapes with tunable widths. Nanotechnology, 17:2813{ 2820, 2006.