We extend a simple dipole approximation model to predict nonlinear scattering from small particles. This numerical method is known as Discrete Dipole Approximation (DDA) and has been extensively used to model linear scattering by small particles of various shapes and sizes. We show here that DDA can be used to efficiently model second harmonic scattering by small particles. Our results are compared with experimental data and other computational methods.
©2010 Optical Society of America
The scattering of light by small particles has been an active area of research for over a century. Scattering by a particle is a function of various factors such as the particle’s shape and size, the distribution (intensity, direction and polarization) of the incident light field, and the orientation of the particle relative to the excitation field. Developing a computational model that can readily account for these factors has become all the more important today because of the wide variety of nanoparticles synthesized and used . Nonlinear optical properties, specifically second harmonic generation (SHG) from small particles are highly dependent on the particle geometry. It is well known that under the dipolar approximation, SHG is forbidden from a centrosymmetric medium. At an interface where this symmetry is broken, SHG can be observed. However, in the case of spherical particles which are small relative to the excitation wavelength, SHG generated from one part of the particle surface is cancelled out by SHG generated from the other parts of the particle surface. An expression for second order polarization in a centrosymmetric medium, proposed by Adler , formed the basis of most bulk models of SHG from small spheres. Agarwal and Jha  proposed the first model for SHG from small metal spheres, predicting a dipolar bulk response and a quadrupolar surface response from the particles. Others have taken into account the dynamics of electrons within a small sphere [4, 5] to predict its second harmonic scattering properties. SHG from an array of quantum dots was calculated by modeling each quantum dot as a particle in a box which is being perturbed by excitation light . All the above works consider the bulk of the particle as a source of SHG. There is SHG even from the surface of a spherical particle if the particle is not small as compared to the excitation wavelength or if the excitation field is not uniform across the particle. The theory of SHG from small spheres under inhomogeneous illumination has been reported [7, 8]. For spherical particles whose size is comparable to the excitation wavelength, the charge induced on the surface by incident light can give rise to higher harmonics . Polar molecules immobilized on the surface of small polystyrene beads when excited, scatter second harmonic light which is polarized and has a specific angular distribution [10, 11]. Second harmonic response from small spheres with arbitrary surface response was calculated by Dadap et. al [12, 13]. All these analytical models give us an insight into the interaction of light with small particles. While these models are important, they are restricted to spherical particles in relatively homogeneous environment. With the increasing importance of non-spherical metal nanoparticles being recognized [14–16], it is important to develop a computational frame work that can efficiently predict the SHG properties of these nonspherical particles. In addition to particle geometry, the illumination conditions may also be far from simple with significant field gradient or complex polarization distribution. To add to this complexity, the local environment of the nanoparticles may also be heterogeneous, especially for biomedical applications. Numerical methods are better equipped to solve these problems.
Here we use a numerical method called the discrete dipole approximation (DDA)  to calculate second harmonic scattering properties of small particles. This method was formulated to simulate linear scattering and absorption by interstellar dust particles [17, 18]. In DDA, a scatterer is assumed to be made up of small polarizable dipoles (Fig. 1 ) which interact among themselves and with the external field. The optical properties of the scatterer are given by the summation optical properties of all the constituent dipoles. Draine and associates developed an efficient Fourier transforms based algorithm for DDA to facilitate fast computation of scattered fields by particles of various geometries [19, 20]. The algorithm is implemented in FORTRAN language and it is available as a free program - DDSCAT . Hoekstra et. al. developed a parallel computing version of DDA which made it possible to simulate light scattering by large particles [22, 23]. DDA has been used in different areas of research but with different names, like coupled dipole method (CDM)  and polarizable dipole model (PDM) . One of the main advantages of DDA is its ability to handle arbitrary shapes of scatterers. The method is easy to implement and it is computationally undemanding. Although DDA was developed for materials with moderate magnitudes of refractive index, it can be used to calculate scattering properties of metal nanoparticles which have high magnitudes of refractive index . The results predicted in the far field and the highly sensitive near field agree well with results from other numerical methods. As such, it has been extensively used to calculate linear optical properties of noble metal nanoparticles, especially for non-spherical shapes, clusters of nanoparticles and nanoparticles in heterogeneous local environments [27–30]. DDA has also been adapted to calculate scattering by a wide variety of objects like periodic scatterers , photonic crystals , magnetic nanoparticles , metamaterials  and blood cells . Apart from scattering, this method has also been used to calculate optical forces on nanoparticles . All these applications of DDA consider linear interaction of light with matter. We have shown here that by taking into account nonlinear interaction between light and matter, DDA can be extended to predict nonlinear optical scattering from small particles. To demonstrate our method we have calculated the spatial distribution of scattered second harmonic light from different types of particles.
When excited by light, the strong near field around nanoparticles suggests enhanced nonlinear scattering from the surface of these particles but calculation of the scattered nonlinear light by a nanoparticle using DDA has not been reported yet. We achieve this by extending the DDA to calculate induced nonlinear dipoles within small particles when excited by light. We briefly describe the linear model and then extend it to the nonlinear regime. To begin with, a scatterer can be assumed to be made up of N small sub-volumes which are arranged on a cubic lattice. If E inc,i is the incident field and P i (1) is the first order polarization induced at the center of the ith sub-volume, then these quantities are related by
The superscript in braces refers to the order of the process. α i,ω is polarizability of the ith dipole at angular frequency ω. A ij is the term defining interaction between ith and jth dipoles. If r i is the position vector of the ith dipole, thenEq. (1) in a compact form
Finding exact solution to Eq. (3) is a challenging task because of the large size of matrix A. But approximate solutions to this equation can be calculated by various forms of conjugate gradient methods which are fast in finding solutions to reasonable accuracy [37,38]. Once the induced dipole moments are known, the scattered field is given by superposition of the individual dipoles fields. One can also calculate extinction cross-section (Cext) and absorption cross-section (Cabs) of the scatterer .
All the above equations describe linear phenomena and we have reproduced them here for the purpose of clarity. We extend this dipole approximation from Eq. (1) to calculate nonlinear scattering by the scatterer. If we carefully look at Eq. (1), we realize that the linear local field is given by
This local field is the actual driving force for the various orders of polarization. Typically, we neglect the higher orders of local field relative to the first order local field because they are comparatively much smaller in magnitude. Since we are interested in second harmonic scattering, we focus on the second order polarization, P (2). If β i is the first hyperpolarizability of the ith sub-volume, then
Equation (8) differs from Eq. (3) mainly in the driving force for the induced polarizations, which appears on the right hand side of these equations. For the second harmonic polarization, the driving force is proportional to the square of the linear local field weighted by first hyperpolarizability tensor. It should also be noted that the interaction term A ij in Eqs. (7) & (8) will be a function of the second harmonic wave number and not the excitation wave number. Incidentally, similar models for SHG from surfaces were reported by Wijers et. al.  and Poliakov et. al. . Wijers et. al. assumed a silicon wafer surface to be made up of 20-80 stacked layers where each layer acts as a dipole. In Wijers’ model, the geometry is simple – a linear arrangement of dipoles. On the other hand, Poliakov et. al.  used their model to calculate various types of induced nonlinear polarization in molecules adsorbed on rough metallic surfaces. The present work, however, differs from these two earlier reports because here we focus on far field second harmonic scattering properties of nanoparticles of different types and geometries. Shape and size of a nanoparticle play important roles in its second order optical properties and we demonstrate this effect by means of simple examples. We also show here that our method can be used for nanoparticle composites. Another dipole model [41, 42] worth mentioning here considers the dipoles to be independent of each other and driven only by the incident field. This model can be termed as uncoupled dipole model and it is an approximation to our nonlinear coupled dipole model. The coupled dipole model is a two step improvement of the uncoupled dipole model. The coupled dipole model first takes into account the corrected linear local field rather than the incident field itself. Then it calculates the corrected second order polarization rather than obtaining it directly by squaring the linear local field.
Since our model is an extension of DDA, it can be used for scatterers of random shapes. The criteria for choosing the number of dipoles in a scatterer or the size of the sub-volumes are the same as described by Draine et. al.  except that one has to consider the smallest wavelength involved in the entire process. In our case, the second harmonic wavelength is the smallest wavelength. The refractive index of gold nanoparticles was obtained from Blanchard et. al. . The exact value of hyperpolarizability (β) of gold is not known. However we can say that the surface of the particle is asymmetric and hence hyperpolarizability for dipoles on the surface is non-zero. We have considered only one component (β surf ,⊥ ⊥ ⊥) of the hyperpolarizability tensor which acts on the electric field normal to the surface (E loc,surf ,⊥) and yields a second order polarization in the same direction .
Results & discussion
Based on the algorithm proposed by Draine and associates , we developed a MATLAB code to implement the nonlinear DDA model. Firstly, we calculated distribution of scattered second harmonic light from a 14 nm gold nanosphere. The excitation source is a plane polarized light of wavelength 800 nm, polarized along x-direction and propagating along positive z-direction. The nanosphere is located at the origin of the coordinate axes and it is assumed to be made up of dipoles which are 1 nm apart on a cubic lattice. We calculated the induced second order dipole moments (P i (2)’s) and the scattered second harmonic field was calculated from these dipoles. The three polarization components (x, y and z) and the total intensity of the scattered second harmonic field are represented on a polar plot (Fig. 2 ). The distance from the origin represents the magnitude in space. The magnitudes have been normalized with the maximum intensity value to give an estimate of relative strengths of field components. These results are in good agreement with the theoretical analysis of second harmonic scattering from small spherical particles reported earlier .
Extending our work to a less trivial case, we apply this DDA formalism to calculate the SHG from non-spherical scatterers. Although DDA models a scatterer as a collection of dipoles, the overall scattering depends on the geometry of the scatterer among other factors. The geometry of the scatterer influences the interaction between the dipoles and therefore the overall scattering profile is a superposition of various multipoles such as dipole, quadrupole and so on. This is clearly evident in the hyper Rayleigh scattering (HRS) experiments on gold and silver nanospheres [45, 46]. In this experimental setup incident light propagates along positive the z-axis. The polarization of incident light is rotated in the xy plane and the polarization angle is measured with respect to the positive x-axis. At each incident polarization angle, the intensity of x-polarized scattered second harmonic light from a suspension of nanospheres was measured along y-axis. In reality, small metal nanospheres are not perfectly spherical. Due to this asymmetry in shape, the scattered second harmonic light is predominantly dipolar is nature as long as the nanospheres are small. With increase in size of a nanosphere, retardation effects make the quadrupolar scattering stronger than the dipolar scattering. A finite element analysis of this HRS experiment showed that even small metal nanoparticles which are perfectly spherical in shape show quadrupolar scattering but a slight deviation from spherical shape in these small metal nanoparticles gives rise to a dipolar scattering . The perfect spherical shape leads to a cancellation of dipolar scattering. Therefore, quadrupolar scattering, though weak in small nanoparticles, shows up. In metal nanoparticles which are not perfectly spherical, the dipolar behavior dominates as long as the nanoparticle is small. With increasing size, the quadrupolar behavior takes over the dipolar behavior. We simulated these HRS experiments using our nonlinear DDA model. In Fig. 3 we show the effect of shape of a nanoparticle on the scattered second harmonic light. The samples were gold nanosphere, 16 nm in diameter and two types of nanorods, each 12 nm in diameter but with lengths of 14 nm and 16 nm. To account for the free rotation of the particle in a suspension, we averaged the intensity of scattered second harmonic light over six different orientations of the particle for each incident polarization angle. These six directions correspond to axes joining the six diametrically opposite pairs of vertices in an icosahedron. Since an icosahedron is one of the five Platonic solids, its vertices sample the entire solid angle of 4π equally. Figure 3 shows the strength of x-polarized second harmonic light, measured along positive y-axis, as we change the polarization angle of incident light in the xy-plane. The second harmonic scattering from nanospheres is predominantly quadrupolar in nature. On the other hand, second harmonic scattering from the gold nanorods tends to be dipolar with increasing aspect ratio of the rods.
To compare our results with experimental results, we simulated HRS experiments on silver nanospheres . The excitation wavelength was set at 780 nm. Scattering samples were silver nanospheres of diameters 40 nm, 60 nm and 80 nm. To match the experimental conditions closely our simulated nanospheres should have all forms of surface and bulk asymmetries which occur in silver nanospheres used for experiments. This, however, is not possible because of the large number of asymmetries. Therefore we used just one asymmetry , elongation in one direction, for all three sizes of nanospheres. With increase in size of the spheres, the degree of asymmetry decreases and the quadrupolar behavior dominates as can be seen in second harmonic scattering patterns (Fig. 4 ). All these results are in good agreement with the experimental results  and numerical simulations .
So far we compared the results from our model with theoretical and experimental results for metal nanoparticles. But our model can also predict nonlinear scattering from dielectric particles. Yang et. al.  reported HRS from a water suspension of polystyrene beads which were coated with malachite green molecules. They used polystyrene beads of 510 nm, 700 nm and 980 nm in diameter and measured position dependent p-polarized second harmonic scattering for s- and p-polarized incident light. In Yang’s experiments, the detector was rotated about the sample from positive z-direction (+) to negative z-direction (-) to measure the spatial distribution of HRS.
To test our model, we simulated these experiments using polystyrene beads of diameters 510 nm, 680 nm and 986 nm. The incident wavelength was set at 840 nm to match the experimental conditions. Owing to the low refractive index of polystyrene, we set the distance between dipoles as 17nm. The excitation light can be either s- or p-polarized, and the second harmonic scattering is a combination of s- and p-polarizations. The experimental results  show that for p-polarized excitation, the p-polarized HRS becomes more forward directed with increase in size of the beads. It was also reported that the p-polarized HRS becomes weaker and spatially less localized when the polarization of excitation light is changed from p to s. Both these results are well captured by our model (Fig. 5 ). The angular localization of the main lobes and the first side lobes, and relative strengths of the HRS intensities match well with the experimental results . For example, the magnitudes of p-polarized HRS from 686nm beads when excited with s-polarized light and that from 510 nm beads when excited with p-polarized light are comparable. One aspect of our results which does not match with the Yang’s results is size of side lobes. This mismatch is partly due to the fact that the polarizability of malachite green molecules is not known quantitatively. We observed that as we increased the polarizability of the surface bound molecules, the side lobes decreased in size. But higher refractive index implies higher computational costs, especially when we have to find the right polarizability by iteration. It is possible to match the experimental results by iteratively choosing the right polarizability but due to the marginal improvement in results, we avoided this extra computational cost.
In summary, we have demonstrated that DDA is a very simple model to predict second harmonic scattering from nanoparticles of various kinds. We have shown how shape and size of nanoparticles influence their nonlinear optical properties. With increasing applications of nanoparticles in optical experiments, this method should be helpful in providing good prediction of nonlinear scattering properties. The model does not restrict itself to specific shapes of scatterers or illumination conditions. We can model a single particle as well as a collection of particles by this method. The way we presented the model here, one can also use it for other nonlinear scattering processes such as third harmonic, CARS, etc. This is just a preliminary report exploring DDA for examples of nonlinear scattering applications.
The authors would like to thank Prof. Jerome Mertz for his valuable comments and Singapore MIT Alliance for funding this work.
References and Links
1. C. L. Nehl and J. H. Hafner, “Shape-dependent plasmon resonances of gold nanoparticles,” J. Mater. Chem. 18(21), 2415–2419 (2008). [CrossRef]
2. E. Adler, “Nonlinear Optical Frequency Polarization in a Dielectric,” Phys. Rev. 134(3A), A728–A733 (1964). [CrossRef]
3. G. S. Agarwal and S. S. Jha, “Theory of second harmonic generation at a metal surface with surface plasmon excitation,” Solid State Commun. 41(6), 499–501 (1982). [CrossRef]
4. X. M. Hua and J. I. Gersten, “Theory of second-harmonic generation by small metal spheres,” Phys. Rev. B 33(6), 3756–3764 (1986). [CrossRef]
5. O. A. Aktsipetrov, P. V. Elyutin, A. A. Nikulin, and E. A. Ostrovskaya, “Size effects in optical second-harmonic generation by metallic nanocrystals and semiconductor quantum dots: The role of quantum chaotic dynamics,” Phys. Rev. B 51(24), 17591–17599 (1995). [CrossRef]
6. A. Guerrero and B. S. Mendoza, “Model for great enhancement of second-harmonic generation in quantum dots,” J. Opt. Soc. Am. B 12(4), 559–569 (1995). [CrossRef]
7. V. L. Brudny, B. S. Mendoza, and W. Luis Mochán, “Second-harmonic generation from spherical particles,” Phys. Rev. B 62(16), 11152–11162 (2000). [CrossRef]
8. W. L. Mochán, J. A. Maytorena, B. S. Mendoza, and V. L. Brudny, “Second-harmonic generation in arrays of spherical particles,” Phys. Rev. B 68(8), 085318 (2003). [CrossRef]
9. D. Östling, P. Stampfli, and K. H. Bennemann, “Theory of nonlinear optical properties of small metallic spheres,” Z. Phys. D At. Mol. Clust. 28, 169–175 (1993). [CrossRef]
10. J. Martorell, R. Vilaseca, and R. Corbalán, “Scattering of second-harmonic light from small spherical particles ordered in a crystalline lattice,” Phys. Rev. A 55(6), 4520–4525 (1997). [CrossRef]
12. J. I. Dadap, J. Shan, K. B. Eisenthal, and T. F. Heinz, “Second-Harmonic Rayleigh Scattering from a Sphere of Centrosymmetric Material,” Phys. Rev. Lett. 83(20), 4045–4048 (1999). [CrossRef]
13. J. I. Dadap, J. Shan, and T. F. Heinz, “Theory of optical second-harmonic generation from a sphere of centrosymmetric material: small-particle limit,” J. Opt. Soc. Am. B 21(7), 1328–1347 (2004). [CrossRef]
14. N. J. Durr, T. Larson, D. K. Smith, B. A. Korgel, K. Sokolov, and A. Ben-Yakar, “Two-photon luminescence imaging of cancer cells using molecularly targeted gold nanorods,” Nano Lett. 7(4), 941–945 (2007). [CrossRef] [PubMed]
15. Y. Jung, H. Chen, L. Tong, and J.-X. Cheng, “Imaging Gold Nanorods by Plasmon-Resonance-Enhanced Four Wave Mixing,” J. Phys. Chem. C 113(7), 2657–2663 (2009). [CrossRef]
17. E. M. Purcell and C. R. Pennypacker, “Scattering and Absorption of Light by Nonspherical Dielectric Grains,” Astrophys. J. 186, 705–714 (1973). [CrossRef]
18. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]
20. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]
21. B. T. Draine and P. J. Flatau, “User Guide for the Discrete Dipole Approximation Code DDSCAT 7.0,” (2008).
22. A. G. Hoekstra and P. M. A. Sloot, “Coupled dipole simulations of elastic light scattering on parallel systems,” Int. J. Mod. Phys. C 6(5), 663–679 (1995). [CrossRef]
23. A. G. Hoekstra, M. D. Grimminck, and P. M. A. Sloot, “Large Scale Simulations of Elastic Light Scattering by a Fast Discrete Dipole Approximation,” Int. J. Mod. Phys. C 9(1), 87–102 (1998). [CrossRef]
24. A. Lakhtakia, “General theory of the Purcell-Pennypacker scattering approach and its extension to bianisotropic scatterers,” Astrophys. J. 394, 494–499 (1992). [CrossRef]
25. N. Arzate, B. S. Mendoza, and R. A. Vázquez-Nava, “Polarizable dipole models for reflectance anisotropy spectroscopy: a review,” J. Phys. Condens. Matter 16(39), S4259–S4278 (2004). [CrossRef]
26. W.-H. Yang, G. C. Schatz, and R. P. Van Duyne, “Discrete dipole approximation for calculating extinction and Raman intensities for small particles with arbitrary shapes,” J. Chem. Phys. 103(3), 869–875 (1995). [CrossRef]
27. N. Félidj, J. Aubard, and G. Levi, “Discrete dipole approximation for ultraviolet–visible extinction spectra simulation of silver and gold colloids,” J. Chem. Phys. 111(3), 1195–1208 (1999). [CrossRef]
28. K. L. Kelly, E. Coronado, L. L. Zhao, and G. C. Schatz, “The Optical Properties of Metal Nanoparticles: The Influence of Size, Shape, and Dielectric Environment,” J. Phys. Chem. B 107(3), 668–677 (2003). [CrossRef]
29. P. K. Jain, K. S. Lee, I. H. El-Sayed, and M. A. El-Sayed, “Calculated absorption and scattering properties of gold nanoparticles of different size, shape, and composition: applications in biological imaging and biomedicine,” J. Phys. Chem. B 110(14), 7238–7248 (2006). [CrossRef] [PubMed]
30. S. W. Prescott and P. Mulvaney, “Gold nanorod extinction spectra,” J. Appl. Phys. 99(12), 123504–123507 (2006). [CrossRef]
31. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for periodic targets: theory and tests,” J. Opt. Soc. Am. A 25(11), 2693–2703 (2008). [CrossRef]
32. F. Bordas, N. Louvion, S. Callard, P. C. Chaumet, and A. Rahmani, “Coupled dipole method for radiation dynamics in finite photonic crystal structures,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 73(5), 056601 (2006). [CrossRef] [PubMed]
35. M. A. Yurkin, K. A. Semyanov, P. A. Tarasov, A. V. Chernyshev, A. G. Hoekstra, and V. P. Maltsev, “Experimental and theoretical study of light scattering by individual mature red blood cells by use of scanning flow cytometry and a discrete dipole approximation,” Appl. Opt. 44(25), 5249–5256 (2005). [CrossRef] [PubMed]
36. P. C. Chaumet, A. Rahmani, A. Sentenac, and G. W. Bryant, “Efficient computation of optical forces with the coupled dipole method,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 72(4), 046708 (2005). [CrossRef] [PubMed]
39. C. M. J. Wijers, T. Rasing, and R. W. J. Hollering, “Second harmonic generation from thin slabs in the discrete dipole approach,” Solid State Commun. 85(3), 233–237 (1993). [CrossRef]
40. E. Y. Poliakov, V. A. Markel, V. M. Shalaev, and R. Botet, “Nonlinear optical phenomena on rough surfaces of metal thin films,” Phys. Rev. B 57(23), 14901–14913 (1998). [CrossRef]
41. L. Moreaux, O. Sandre, and J. Mertz, “Membrane imaging by second-harmonic generation microscopy,” J. Opt. Soc. Am. B 17(10), 1685–1694 (2000). [CrossRef]
42. J.-X. Cheng and X. S. Xie, “Green's function formulation for third-harmonic generation microscopy,” J. Opt. Soc. Am. B 19(7), 1604–1610 (2002). [CrossRef]
43. N. P. Blanchard, C. Smith, D. S. Martin, D. J. Hayton, T. E. Jenkins, and P. Weightman, “High-resolution measurements of the bulk dielectric constants of single crystal gold with application to reflection anisotropy spectroscopy,” Phys. Status Solidi 0(8c), 2931–2937 (2003). [CrossRef]
44. G. Bachelier, I. Russier-Antoine, E. Benichou, C. Jonin, and P.-F. Brevet, “Multipolar second-harmonic generation in noble metal nanoparticles,” J. Opt. Soc. Am. B 25(6), 955–960 (2008). [CrossRef]
45. J. Nappa, G. Revillod, I. Russier-Antoine, E. Benichou, C. Jonin, and P. F. Brevet, “Electric dipole origin of the second harmonic generation of small metallic particles,” Phys. Rev. B 71(16), 165407 (2005). [CrossRef]
46. I. Russier-Antoine, E. Benichou, G. Bachelier, C. Jonin, and P. F. Brevet, “Multipolar Contributions of the Second Harmonic Generation from Silver and Gold Nanoparticles,” J. Phys. Chem. C 111(26), 9044–9048 (2007). [CrossRef]