Comparison of numerical methods in near-field computation for metallic nanoparticles

Four widely used electromagnetic field solvers are applied to the problem of scattering by a spherical or spheroidal silver nanoparticle in glass. The solvers are tested in a frequency range where the imaginary part of the scatterer refractive index is relatively large. The scattering efficiencies and near-field results obtained by the different methods are compared to each other, as well as to recent experiments on laser-induced shape transformation of silver nanoparticles in glass. © 2011 Optical Society of America OCIS codes: (050.1755) Computational electromagnetic methods; (290.0290) Scattering. References and links 1. G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys.-Leipzig 330, 377–445 (1908). 2. A. Taflove and S. C. Hagness, Computational Electrodynamics, Third Edition. (Artech House 2005). 3. J. Niegemann, M. Konig, K. Stannigel, and K. Busch, “Higher order time-domain methods for the analysis of nano-photonic systems,” Photon. Nanostructures 7(1), 2–11 (2009). 4. P. Monk, Finite Element Method for Maxwell’s Equations, (Oxford, 2006). 5. A. C. Cangellaris and D. B. Wright, “Analysis of the numerical error caused by stair-stepped approximation of a conduction boundary in FDTD simulations of electromagnetic phenomena,” IEEE Trans. Antenn. Propag. 39(10), 1518–1525 (1991). 6. D. W. Lynch and W. R. Hunter, “Comments on the Optical Constants of Metals and an Introduction to the Data for Several Metals,” in Handbook of Optical Constants of Solids, vol. 1, E. D. Palik, ed (Academic, San Diego, 1985). 7. A. Vial and T. Laroche, “Description of dispersion properties of metals by means of the critical points model and application to the study of resonant structures using the FDTD method,” J. Phys. D Appl. Phys. 40(22), 7152– 7158 (2007). 8. A. Doicu, T. Wriedt and Yu. Eremin, Light Scattering by Systems of Particles, Null-Field Method with Discrete Sources: Theory and Programs (Springer 2006). 9. A. Doicu and T. Wriedt, “Near-field computation using the null-field method,” J. Quant. Spectrosc. Radiat. Transf. 111(3), 466–473 (2010). 10. B. T. Draine and P. J. Flatau, “Discrete dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). 11. B. T. Draine and P. J. Flatau, “User Guide to the Discrete Dipole Approximation Code DDSCAT 7.1”, http://arXiv.org/abs/1002.1505v1 (2010). 12. COMSOL Multiphysics demonstration CD-ROM can be requested at http://www.comsol.com 13. C. Hafner and L. Bomholt, The 3D electromagnetic wave simulator (Wiley 1993). 14. A. Stalmashonak, A. Podlipensky, G. Seifert, and H. Graener, “Intensity-driven, laser induced transformation of Ag nanospheres to anisotropic shapes,” Appl. Phys. B 94(3), 459–465 (2009). 15. A. Stalmashonak, G. Seifert, and H. Graener, “Spectral range extension of laser-induced dichroism in composite glass with silver nanoparticles,” J. Opt. A, Pure Appl. Opt. 11(6), 065001 (2009). 16. A. Stalmashonak, C. Matyssek, O. Kiriyenko, W. Hergert, H. Graener, and G. Seifert, “Preparing large-aspectratio prolate metal nanoparticles in glass by simultaneous femtosecond multicolor irradiation,” Opt. Lett. 35(10), 1671–1673 (2010). #140555 $15.00 USD Received 4 Jan 2011; revised 11 Mar 2011; accepted 16 Mar 2011; published 22 Apr 2011 (C) 2011 OSA 25 April 2011 / Vol. 19, No. 9 / OPTICS EXPRESS 8939 17. W. J. Wiscombe, “Improved Mie scattering algorithms,” Appl. Opt. 19(9), 1505–1509 (1980). 18. D. Gutkowicz-Krusin and B. T. Draine, “Propagation of electromagnetic waves on a rectangular lattice of polarizable points”, http://xxx.arxiv.org/abs/astro-ph/0403082 (2004). 19. M. A. Yurkin, “Discrete dipole simulations of light scattering by blood cells”, Dissertation (2007), ISBN 90– 5776–169–6


Introduction
Design and study of modern metallic nanostructures often require accurate near-field computation, a task that can be tackled by a number of computational methods.In this paper, we compare different methods regarding their accuracy and computational time and also use them for the interpretation of experimental data.Considering general nanostructure shapes quickly rules out purely analytic methods such as Mie theory [1], which is restricted to simple shapes like spheres.Numerical time domain methods as Finite-Difference Time Domain (FDTD) [2], the Discontinuous Galerkin Time Domain (DGTD) [3] and the Finite Element Method (FEM) [4] allow better representation of arbitrarily shaped structures.The latter two methods, in contrast to FDTD, rely on an unstructured mesh, which allows a more accurate representation of curved surfaces.This not only increases efficiency, but also avoids unphysical behavior of the fields due to the staircase approximation in FDTD [5].The representation of material properties, especially of a frequency-dependent permittivity as in the case of metals, also depends on the methods.Comparing numerical results with experiments suggests the use of measured values for the permittivity [6].Because time domain methods usually require an approximation of the permittivity with analytic models like the Drude-Lorentz (see e.g [2].) or Critical Points model [7], this rules out time domain and restricts our choice to frequency domain methods.Among these are semi-analytical methods, which express the electric field in terms of basic functions.T-Matrix type methods usually rely on a single vector spherical harmonics expansion, Multiple Multipole (MMP) methods make use of several expansions with freely chosen origins, and the Discrete Dipole Approximation (DDA) employs a large number of dipoles to represent the scatterer.The more general FEM in frequency domain approximates the electric field in terms of low order polynomials on small volumes.We compare the following implementations: the Null Field Method with Discrete Sources (NFM-DS) of T-Matrix type [8,9], DDSCAT 7.0 [10,11], which is a DDA implementation, the FEM provided by COMSOL Multiphysics® 3.5 [12] and MMP 3D [13] as the representative of MMP methods.This comparison is done not only to validate the commonly used numerical methods, but also to informally indicate how the proper choice of the method can drastically reduce the required computation time.The chosen test scattering problem involves a silver spherical or spheroidal nanoparticle (NP) embedded in glass and illuminated in the wavelength range 250nm-900nm.This scattering problem is challenging because of the expected field enhancement near the poles of the NP, and because the imaginary part of the refractive index of silver is relatively large in the above wavelength range, see Fig. 1a.Apart from the purely computational interest, the test scattering problem is highly relevant in the interpretation of recent experimental findings on laser-induced shaping of silver NPs, as explained in Section 3. The test scattering problem is described in detail in Section 2. The experimental setup for the shaping of silver NPs is described in Section 3. Section 4 contains a very brief description of the four considered electromagnetic field solvers.Numerical results are presented and discussed in Section 5, and conclusions are stated in Section 6.

The test scattering problem
Figure 1b shows an instance of the considered scattering problem, which in fact consists of two special cases.In both cases, the scatterer is of silver and has the volume of a ball of radius 15nm.In the first case, the boundary of the scatterer is spherical, and in the second case it is prolate spheroidal with the ratio of the longer semi-diameter to the shorter given by 2.4: 1.In the coordinate system of Fig. 1b, the boundary of the scatterer is described by the implicit equation (x/a) 2 + (y/a) 2 + (z/b) 2 = 1, with a = b = 15nm in the first case and a 11.2nm, b 26.9nm in the second case.The scatterer is immersed in glass with refractive index 1.5, and illuminated by a time-harmonic, z-polarized, uniform plane wave of amplitude 1 V/m and propagating in the positive x-direction.The free-space wavelength λ 0 of the incident field is allowed in the range 250nm-900nm.Figure 1a shows the real and the imaginary part of the refractive index n Ag of silver for these wavelengths, obtained by spline interpolation of the data from Lynch and Hunter [6].

Laser-induced shaping of silver nanoparticles
Two of the authors (G.S. and A.S.) have recently extensively studied the interaction of linearly polarized ultrashort femtosecond laser pulses with silver NPs embedded in glass.The typical samples, corresponding largely to the case studied here numerically, consist of sodalime glass (n = 1.5) and contain spherical silver NPs of approximately 15nm radius in a layer near the sample surface.The average filling factor is f = 0.001, so that interaction between individual NPs are negligible.It was discovered that these NPs experience a persisting transformation into spheroidal shape upon irradiation with several hundred intense femtosecond laser pulses at a wavelength near to the surface plasmon (SP) resonance at 400nm wavelength [14,15].Macroscopically, the shape change into uniformly oriented spheroids is observed as optical dichroism.Depending on the applied intensity and number of pulses, spheroids with maximal aspect ratio of 1.8 could be obtained by near-resonant excitation (λ 400nm).If under these conditions the intensity or number of pulses applied to one position is getting too high, one observes partial destruction (being reflected in a bleaching of the SP bands) instead of further elongation of the NPs [14].
To achieve larger aspect ratios, different irradiation wavelengths have to be used.In particular, even rather strongly red-shifted, "off-resonant" excitation (with respect to the initial SP band) can induce a very effective shape transformation of the NPs, in spite of the relatively low remaining SP absorption [15].Moving towards longer wavelengths, a maximum aspect ratio of 2.4 is being achieved around λ = 600nm; then, increasing the irradiation wavelength further, the efficiency of the nanoparticle shape transformation decreases more and more, and finally (for λ > 670nm) the laser pulses do not evoke any measurable extinction changes any longer.Even the very strong laser fundamental at 800nm does not cause dichroism when irradiating a sample with original, spherical silver NPs.
We found, however, a very instructive technique to obtain even larger aspect ratios by using simultaneous irradiation with two different wavelengths.For example, using a strong pulse at 532nm and a much weaker one at 800nm wavelength, aspect ratios of up to 4 were obtained [16].This rather unexpected result appears to be a crucial point to test and refine our model for the physical processes behind the step-by-step shape transformation of Ag NPs under irradiation with several hundred intense fs laser pulses.The qualitative picture, being suitable to explain all experimental findings has been proposed in [14] and is based on the electric field enhancement at the particle-glass interface.Key feature are two emission processes: (i) anisotropic emission of electrons due to the electric field of the laser pulse, strengthened by the field enhancement at the particle surface, and (ii) isotropic emission of electrons and (later) silver ions due to the heating of the particle.In the further course of events, emitted electrons and silver ions recombine and diffuse back towards the remaining particle, leading to material deposition favorably at the poles.
Therefore, anisotropic emission of electrons appears to be the driving force for anisotropic shape transformation, and the electric field distribution around the nanoparticle seems to be the decisive point governing the maximum anisotropy achievable by the transformation process.Hence, we will study the field enhancement for the two shapes outlined in section 2, i.e. for the cases of a sphere and a spheroid with an aspect ratio of 2.4.The latter refers to the maximum aspect ratio, which can be achieved by fs laser-induced shape transformation using a single wavelength irradiation.

Numerical methods
Four field solvers are here applied to the test scattering problem: the Null Field Method with Discrete Sources (NFM-DS) [8,9], DDSCAT 7.0 [10,11], COMSOL Multiphysics Suite 3.5 [4,12], and MMP 3D [13].All the programs except COMSOL are available free of charge.NFM-DS is an implementation and extension of the well-known T-Matrix Method [8], and it is applicable to dielectric isotropic, anisotropic, chiral and perfectly electrically conducting particles of arbitrary shape.It is also applicable to layered/inhomogeneous particles and particles on surfaces, as well as to aggregates of particles.The T-Matrix Method is based on expressing the incident field E i , the field in the scatterer interior E int and the scattered field E s in terms of truncated series of spherical vector wave functions (SVWF), Here 1 M  , 1 N  are regular SVWF, and The transition matrices are computed for the considered particle (scatterer) using null-field equations [8].The so-called T-Matrix int T s int  T i maps the coefficients of the series for the incident field to the coefficients of the series for the corresponding scattered field.According to the Wiscombe criterion [17] a spherical scatterer of radius 15nm and refractive index 1.5 requires the parameter N of 3 or 4. For non-spherical particles of the same volume, or when the electromagnetic field needs to be evaluated close to the scatterer surface, the number of multipoles should be increased for accuracy.We use the maximal order N = 8 of SVWF in the field expressions.For the near-field data of Section 5, the radiating integrals are computed by discretizing the generatrix of the particle surface by 25 up to 400 points.DDSCAT 7.0 is an electromagnetic field solver based on the Discrete Dipole Approximation (DDA).It is applicable to scattering and absorption by isotropic and anisotropic dielectric single or multiple particles of arbitrary shape.A particle is modeled by an array of N polarizable point dipoles.The dipoles are placed at points r μ in a cubic lattice inside the particle volume and with lattice spacing d.For isotropic materials, the polarization vector  P of the dipole at lattice point r μ is assumed equal to E exc   , where the complex number   is the polarizability of the dipole, and where is the sum be written as the system of linear equations 1 , 1, , in the polarization vectors P  .After the system is inverted, the near-field and far-field quantities of the given scattering problem can be computed using the above expression for  A and summing up the field contributions from the individual dipoles.The expression for   depends on the chosen polarizability model for the point dipoles.Several models are commonly used for DDA; DDSCAT employs the Lattice Dispersion Relation as treated in [18].Chapter 3 of [19] contains an overview and a detailed discussion of the DDA formulation of scattering problems, including polarizability models.Two validity criteria for DDSCAT are stated in [10,11]: that the interdipole separation d be small compared to all structural lengths of the particle (so that the dipole array models the shape of the particle with sufficient accuracy), and that d be sufficiently small compared to the plane-wave wavelength in the particle material.The second criterion may be satisfied [10] if dn , where n is the complex refractive index of the particle material and 0  is the free-space wavelength of the incident field.For noble-metal particles with relatively large imaginary part of the refractive index, the interdipole separation must be further reduced.For the near-field and farfield data of Section 5, arrays of at least 20 × 20 × 20 dipoles and at least 20 × 20 × 48 dipoles are used to model the spherical and the prolate spheroidal particle, respectively, giving the largest used lattice spacing of d = 1.5nm.The structural lengths of the considered particles are thus at least 20 times larger than the used interdipole separation, and, since COMSOL Multiphysics is a commercially available integrated environment for modeling, solution and postprocessing of a wide variety of boundary problems.It employs the Finite Element Method (FEM) [4], where the system under consideration is represented by a mesh consisting of tetrahedral elements.This mesh allows accurate and efficient geometric representation of the scatterer.On each tetrahedral element, the solution, in our case the components of the electric field, are approximated by an expansion in terms of basis functions, usually polynomials up to third or fourth degree.The elements are coupled to their neighbors by the boundary conditions, which leads to a system of linear equations.The solution of this system of equations yields the expansion coefficients and thus the electric field is obtained.As the details of this method are very involved and depend on the specific implementation, we refer the reader to a classic textbook [4] or, e.g. the COMSOL manual.The results shown in Section 5 are computed using the iterative GMRES solver.Due to symmetry, it suffices to discretize only ¼ of the NP and use a mesh with 43244 tetrahedral elements.This results in a minimum element quality of 0.1387 and an element volume ratio of 3.09•10 5 .Finally, MMP 3D is a semi-analytic code based on the Multiple Multipole (MMP) expansion.This solver is applicable to scattering problems involving piecewise linear, homogeneous and isotropic domains.In MMP 3D, the fields are approximated by linear combinations of elementary solutions of the Maxwell system, such as plane waves, waveguide modes and 2D and 3D multipoles.The availability of different types of elementary solutions allows efficient formulation of practical problems, tailored to the given geometry.In a typical scattering formulation, the scattered field and the field in the interior of the scatterer are expressed in terms of sums of multipoles with different origins of expansion, ( ) (The superscript (j) signifies the j"th origin of expansion.)The coefficients in the sums are found by generalized point matching at the scatterer boundary.The number of testing points where the boundary condition is enforced is usually much larger than the number of used multipoles.For the numerical results presented below, the boundary condition is tested at 255 points on the particle surface, and a total of n = 6 multipole solutions are used in the representation of E s and int E .The origins of expansion of the multipoles are distributed along the z-axis in the interior of the particle.The symmetry of the scattering problem with respect to the xy-plane is exploited to halve the actual number of unknown coefficients.The radiating and regular multipoles have coinciding origins of expansion, and they employ SVWFs of maximal order 20, 11 and 6, and maximal degree 7, 11 and 6, respectively.The maximal order of the SVWFs could be reduced without a significant increase in the solution error, but, as seen in Section 5, the near-field execution time for MMP 3D is relatively short even for the above choice of the number of SVWFs used.

Numerical results
The comparison of the numerical results in this section is done for the solver parameters fixed in Section 4. The execution times stated below are included for illustration only, since the solvers were run on different computers.The NFM-DS was run on a single core of a 2.33GHz Intel® Xeon® processor with 16Gb RAM, DDSCAT and MMP 3D used a single core of a 2.66GHz Intel® Core2 Duo CPU with 4Gb RAM, and COMSOL was run on a 2.4GHz Intel® Core2 Quad processor with 8Gb RAM.For the near-field comparisons, the free-space incident wavelengths λ 0 = 410nm and λ 0 = 688nm are used, since they correspond to off-resonance and on-resonance values, respectively.The computed total electric near fields are sampled on subsets of a fixed rectangular grid G in the xz-plane, i.e., in the plane spanned by the direction of propagation of the incident field and by the incident polarization.The 121 × 86 grid consists of points (x j , z j ) spaced equidistantly between 60nm and 60nm along the x-axis and between 0nm and 85nm along the z-axis (the spacing is 1nm).Since, as it turns out, the quality of the computed near field can deteriorate significantly as the observation point approaches the NP surface, it is convenient to have a measure of discrepancy between the near fields that disregards the field values in a well-defined, variable neighborhood of the NP.Thus, for the spherical scatterer the 2  -norm of a field E is in the following defined by  around the NP in the xz-plane.In both cases, the closer the parameter ε is to 1 the more the field in the immediate vicinity of the scatterer is taken into account when calculating the field norms and distances.Figure 2 shows the relative 2  near-field error  For the shown values of ε, the relative near-field error is less than 10% for DDSCAT and less than 1.8% for NFM-DS and COMSOL.The quality of the near field decreases significantly for all three methods as observation points closer to the scatterer boundary are included (that is, as ε approaches 1).In particular, while NFM-DS gives by far the overall best near field, the decrease in quality as ε  1 is so large for this method that the relative nearfield error for NFM-DS is close to that produced by COMSOL if all sampling points in the grid G outside the scatterer are taken into account.Next, Fig. 3 shows the scattering efficiency of the prolate spheroidal scatterer as function of the incident free-space wavelength (from 250nm to 900nm), as calculated by all four solvers.NFM-DS and COMSOL agree well with the reference MMP 3D solution over almost the entire shown wavelength spectrum.There is very good agreement between the DDSCAT and MMP results in the range 250nm-400nm.However, DDSCAT with a 20 × 20 × 48 array shows a relatively large error in the range 550nm-650nm, and it overestimates the scattering efficiency for longer wavelengths.It turns out that modeling the scatterer with an array of 40 × 40 × 96 points shows only a relatively small improvement in the computed scattering efficiency, and that the number of polarisable points needs to be much larger (e.g., 80 × 80 × 192, corresponding to the lattice spacing d = 0.28nm and


) to obtain satisfactory far-field results over the whole range of wavelengths.For each wavelength, the computation of the scattering efficiency takes on average about 0.2 sec for NFM-DS, 18 sec for COMSOL Multiphysics, 20 sec for DDSCAT 20 × 20 × 48, and 1750 sec for DDSCAT 80 × 80 × 192.At 318nm and 688nm, the two wavelengths considered in detail in the following, MMP 3D takes about 44 sec to calculate the scattering efficiency.The DDSCAT far-field computation takes longer for wavelengths where the imaginary part of the refractive index of silver is relatively large. near-field error for NFM-DS, DDSCAT and COMSOL as function of ε for the prolate spheroidal scatterer at λ 0 = 318nm and λ 0 = 688nm, respectively.Again, the near-field error is relative to the MMP 3D result, and it is plotted on a logarithmic scale.For NFM-DS, results with 50 and 400 generatrix discretization points are shown, and for DDSCAT results with arrays of 20 × 20 × 48 and 40 × 40 × 96 point dipoles are given.For the shown values of ε, and for methods other than NFM-DS 50, the relative near-field error is less than 10% at 318nm and less than 18% at 688nm.For large ε, at λ 0 = 318nm, the relative error seems to settle at 10 6 for NFM-DS 50 and NFM-DS 400, at 10 5 for DDSCAT 40 × 40 × 96 and at 10 3 for COMSOL.At λ 0 = 688nm, the relative error seems to settle at 10 -4.5 0.003% for NFM-DS 50 and NFM-DS 400, at 10 2 for DDSCAT 40 × 40 × 96 and at 10 -1.65 2.2% for COMSOL.The difference between the error levels in Figs. 4 and 5 is likely due to the fact that n Ag = 0.93 + i0.51 (with time dependence it e   ) at 318nm and n Ag = 0.14 + i4.44 at 688nm, that is, Im n Ag 0.55 Re n Ag at 318nm while Im n Ag 32 Re n Ag at 688nm.The quality of the near field decreases significantly for NFM-DS, DDSCAT and COMSOL as ε approaches 1.While NFM-DS produces by far the smallest error overall, COMSOL appears to give the best-quality near field closest to the scatterer boundary.For ε close to 1, the relative near-field error produced by NFM-DS 400 is approximately 3% (worse than DDSCAT) at 318nm and approximately 10% (better than DDSCAT) at 688nm.In contrast, the corresponding COMSOL error is approximately 0.2% at 318nm and approximately 8% at 688nm.Note that, for both incident wavelengths, the COMSOL relative error varies the least with respect to the parameter ε.This, as well as the relatively good quality of the COMSOL near field close to the particle, could be due to the fact that the COMSOL mesh density increases as the distance to the particle decreases, while the other methods rely on radiating sources at the particle surface or inside the particle volume and do not compensate for the proximity of the observation point to the particle.Figures 6 and 7 show the relative near-field error for NFM-DS with 25 up to 400 generatrix discretization points, for λ 0 = 318nm and λ 0 = 688nm, respectively.As expected, in both cases the near-field quality increases with better discretization.The calculation of the near field over the grid G took (for both incident wavelengths) approx.16     Figure 8 shows the logarithm of the total near-field amplitude in the xz-plane (see the coordinate system shown in Fig. 1b) obtained at λ 0 = 318nm by the NFM-DS, DDSCAT and MMP 3D, respectively.Figure 9 shows these data in the yz-plane.The near field is evaluated on an 0.5nm rectangular mesh with points spaced equidistantly along the x and z axes.At 318nm incident wavelength, n Ag = 0.93 + i0.51 and the imaginary part of n Ag attains a global minimum in the  chosen wavelength range (see Fig. 1a).As attested by the above analysis in terms of 2  distances, the near fields exterior to the scatterer are in good agreement; in particular, they show local minima near the poles of the prolate spheroid.The interior fields of NFM-DS, COMSOL and MMP agree, while the interior field produced by DDSCAT has a significantly  lower amplitude and appears more homogeneous in comparison.Figures 10 and 11 show the near-field results in the xz-plane and the yz-plane, respectively, obtained at λ 0 = 688nm incident wavelength.Again, the log 10 -values of the total near-field amplitudes are shown.For the incident wavelength of 688nm, n Ag = 0.14 + i4.44 and the imaginary part of n Ag is substantialabout 32 times larger than the real part.The exterior fields are in good agreement, with significant field enhancement near the poles of the scatterer and local minima in its equatorial plane.NFM-DS, COMSOL and MMP 3D agree in the interior of the scatterer, while the DDSCAT field shows unphysical features there.In the remainder of this section, we will be concerned with calculations related to the experiments outlined in Section 3. Our choice falls to COMSOL because it has shown good results in comparison to the other methods plus it comes with several postprocessing abilities that will be useful in the following.Our aim is to study the enhancement of the normal component of the electric field averaged over the particle surface

  nE n F = dA
This quantity is a measure for the strength of the anisotropic emission of electrons from the particle and therefore indicates the "efficiency" of the transformation process.In Fig. 12 we examine the enhancement F n as a function of the aspect ratio.The calculations were done for the two wavelengths used simultaneously in the experiment [16], 532nm and 800nm, and for wavelengths 10nm below and above these values.All curves are normalized to the enhancement value of the sphere at the respective wavelength.Anticipating that the experiment causes a step-by-step shape transformation towards larger aspect ratio with growing number of applied laser pulses, we can read the x axis of Fig. 12 also as being roughly proportional to processing time.Following first the wavelengths around 532nm, we notice that the field enhancement of the sphere is sufficient to start the process.Then, with growing pulse number a peak enhancement is reached around an aspect ratio of 1.75; further on, the field enhancement decreases and for aspect ratios of 2.4 -2.5 it drops below the initial value.This is in good agreement with single wavelength irradiation experiments, where maximum aspect ratios at typically these values were obtained [14].It is obvious that continued irradiation would only further decrease the anisotropic electron emission process, so that isotropic thermal emission finally prevails, preventing any further NP elongation.In the case of simultaneous two-color irradiation, however, the 800nm component apparently starts to take control over the anisotropic emission from this point on: while for wavelengths around 800nm there is no significant enhancement for aspect ratios below 1.75 (Fig. 12), it increases for larger values and already around 2.3 exceeds the value for wavelengths around 532nm.We therefore conclude that the field enhancement is indeed the essential parameter providing sufficient anisotropic electron emission, and as such can be considered the driving force of the transformation process.Using two wavelength irradiation it is possible to obtain larger aspect ratios, because the long wavelength pulse yields a good field enhancement, whereas the short wavelength pulse fails to do so for larger aspect ratios.So the calculations presented here are a strong confirmation of the previously developed model for the NP shape transformation.

Conclusion
We compared the far-field and near-field performance of four electromagnetic field solvers for silver nanoparticles in glass.
In the far field, the NFM-DS and COMSOL show very good correspondence with the reference MMP solution over a wide range of wavelengths and scatterer refractive indices.DDSCAT shows relatively large far-field error in the incident wavelength range 550nm -650nm, but is otherwise in relatively good agreement with the MMP data.There is, in general, good correspondence in the external near field obtained by the four numerical methods.In particular, the behavior of the field close to the poles of the spheroidal scatterer is computed/predicted consistently.However, it is also shown that the quality of the computed near field decreases near the scatterer boundary for NFM-DS, COMSOL and DDSCAT.Thus, caution should be exercised when computing near fields of silver nanoparticles very near the scatterer surface.Comparison of the interior near fields computed by DDSCAT 7.0 at wavelengths 318nm and 688nm indicates that the method is sensitive to large values of the imaginary part of scatterer refractive index.The results of the calculations were also successfully tested against recent experimental results of creating large aspect ratio silver nanoparticles in glass by simultaneous two-color fs laser irradiation.The field enhancement obtained by COMSOL as a function of aspect ratio and irradiation wavelength is in good agreement with the idea that the permanent shape transformation of the NPs from spheres to spheroids is controlled by the amount of anisotropic, field-driven electron emission.Insofar experiments and calculations give a self-consistent picture, which nicely explains many previously unresolved aspects of laser-induced shape transformation of initially spherical metallic Ag nanoparticles.

3  M , 3 N
are radiating SVWF.The coefficients in the above series are linked by the so-called transition matrices int 

2  2 
-norm disregards the near field in the closed annulus around the spherical NP in the xz-plane.Similarly, for the prolate spheroidal scatterer the 2 -norm of a field E is defined for 1 < ε < -norm disregards the near field in the closed set DS, DDSCAT and COMSOL as function of ε for the spherical scatterer at λ 0 = 410nm, where we have chosen the semi-analytic MMP as reference solution.(MMP 3D is chosen as reference because of the absence of discretized surface integrals in that method.Analytic expressions for scattering by spheroidal particles exist, but, in contrast to the spherical case, the evaluation of the near field is not easily available.)The near-field error is plotted on a logarithmic scale.In NFM-DS, 400 discretization points are used along the generatrix of the particle surface, and the DDSCAT results are shown for arrays of 20 × 20 × 20 and 40 × 40 × 40 point dipoles modeling the particle.

Figures 4
Figures 4 and 5 show the relative 2  near-field error for NFM-DS, DDSCAT and COMSOL as function of ε for the prolate spheroidal scatterer at λ 0 = 318nm and λ 0 = 688nm, respectively.Again, the near-field error is relative to the MMP 3D result, and it is plotted on a logarithmic scale.For NFM-DS, results with 50 and 400 generatrix discretization points are shown, and for DDSCAT results with arrays of 20 × 20 × 48 and 40 × 40 × 96 point dipoles sec for NFM-DS 25, 31 sec for NFM-DS 50, 60 sec for NFM-DS 100, 115 sec for NFM-DS 200 and 238 sec for NFM-DS 400.For DDSCAT 40 × 40 × 96, this calculation took approx.297 sec at 318nm and 619 sec at 688nm.COMSOL took 126 sec at 318 nm and 156 sec at 688nm, and MMP 3D took 6 sec.The COMSOL times include mesh generation and assembly.

Fig. 12 .
Fig. 12. Surface average of the normal field component enhancement as a function of the aspect ratio for wavelengths around those used in the experiment.Enhancements are normalized to the values for a sphere at the respective wavelengths.
at the lattice point r μ ) of the incident electric field and the properly retarded electric fields radiated by all the other dipoles in the array.It holds with the time dependence