Evaluation of E. M. Fields and Energy Transport in Metallic Nanoparticles with Near Field Excitation

We compare two ways of calculating the optical response of metallic nanoparticles illuminated by near field dipole sources. We develop tests to determine the accuracy of the calculations of internal and scattered fields of metallic nanoparticles at the boundary of the particles and in the far field. We verify the correct transport of energy by checking that the evaluation of the energy flux agrees at the surface of the particles and in the far field. A new test is introduced to check that the surface fields fulfill Maxwell's equations allowing evaluation of the validity of the internal field. Calculations of the scattering cross section show a faster rate of convergence for the principal mode theory. We show that for metallic particles the internal field is the most significant source of error.


INTRODUCTION
Many of the recent developments in Nanophotonics imaging and sensingare based on the interaction of metallic particles with sources ofradiation located at sub wavelength distances from the particles [1].The necessity to understand howto optimize experimental set-ups and to extract the opticalpropertiesof the nanoparticles from the experimental results has been a strongdriver of the demand for accurate modeling of the interaction betweenincoming light and nanoparticles [2].InNear Field Optical Microscopy in illumination mode one is interestedin calculating the light originating from near field interactionsafter it passes into the far field region, where the detector isplaced [3].For other forms of microscopy and for surface enhancedspectroscopy and sensing, one needs to find the energy flux near thesurface of the nanoparticles.Because fully analytical calculationsare possible only for the few shapes for which the Maxwell's equationsadmit separation of variables, it is important to develop tests forassessing the ability of different methods to calculate quantities ofinterest such as cross sections, field intensities and energy fluxesthat have different convergence rates with respect to computationalparameters.Several efficient techniqueshave beendeveloped to study scattering in non-spherical particles [4][5][6][7][8].In this paper we compare two implementations of therecently developedtheory of the principal modes (TPM) for internal andscattering fields [9][10][11][12], with the Discrete Sources Method(DSM) [13], which is very fast and able to calculate the fields at any point in space.All these algorithms are able to treat non-sphericalparticles, and are based on the decomposition of internal andscattered fields into sums of fields produced by electric and magnetic multipoles distributed inside and/or outside the particles [13][14][15].The principal mode theory issemi-analytical and based on the decomposition of internal andexternal fields into orthogonal modes which are the generalization ofMie's solutions [16] tonon-spherical particles and whose amplitudes are found by projectingthe incident fields on the modes themselves.The explicitdeterminations of internal and scattering modes are used to findresonances [9] and develop control methods [12].On the contrary, theDSM determines the amplitude of internal and scattered multipoles bysolving an overdetermined system of equations.For these methods wecompare the error in the boundary conditions at the surface and wefind how well these methods satisfy exact relations for the scatteredenergy flux and for internal and scattered fields.

METHODOLOGY
We compare the performances of TPM and DSM using the same number and distribution of electric and magnetic multipoles inside and outside the particles and the same set of points on the surface of the particles.We have used a grid of 8,000 points on the generatrix of the surface [10,17] for all the results shown in this paper and we have checked the numerical convergence, see thediscussion in the following section.
Theinternal and scattered fields excited by a given incident field are determined by satisfying the boundary conditions, where p is a point on the surface and n is the unit vector normal to the surface, are sets of  electric and magnetic multipolessummed to represent the internal and scattered fields, indexed i and s respectively,at that point.Sampling the multipoles on the generatrix of the surface leads to the   m n  matrix L , where m is the number of sampling points and n the number of multipoles, which is the same starting point for all the methods used here.Gauss-Legendre quadrature points are used both for the sampling of the fields and also for surface integration along the generatrix line [10].
The DSMdirectly solves for the expansion coefficients of the multipoles, / i s c  , in a least squares sense typically by using Gaussian elimination.Alternatively the over-determined set of linear equations Lx f  can also be solved by using the decomposition, L QR  , where Q is a square matrix whose columns are orthogonal and R is an upper triangular matrix [18].The number of columns of Q is the rank of R , i.e. the dimension of the largest invertible minor of R .Theoretically the multipoles used are linearly independent, so Q should have n columns if using exact numerical precision.In practice, some of the multipoles give rise to columns of L that appear linearly dependent when usingfinite numerical precision; QR algorithms where the number of columns of Q are determined by a user defined upper bound on the ratio between the largest and the smallest eigenvalues (the estimated condition number) of R , are available [18].The eigenvalues of R that would give rise to a poorer condition number are removed, and the corresponding columns of Q are eliminated.This procedure effectively reduces the number of the functions used to span the solution by eliminating the functions that are most effected by numerical noise.
The TPM method instead constructs n pairs of internal and scattered modes that are orthogonal on the surface, each consisting of a linear combination of themultipoles.This is achieved by consideringsubmatrices   n n n   ) formed by sampling the internal and scattered modes, respectively.We then find internal and scattered orthonormal modes using either the decompositions U U to find the principal modes, i.e. two sets of internal and scattered modes that are correlated pairwise on the surface of the particle.For each incident field, the internal and scattered fields are found by projecting the incident fields on to the principal modes [9].The amplitude of then th internal principal modeis obtained by using the expression, Similarly the amplitudes of the modes in the scattered space can be obtained by changing the sign of Eqn.(2) and exchanging the sets of principal internal and scattered modes, n i and n s respectively, where the projection of the incident field is given by .Note that unlike the DSM, in the principal mode theory one can control separately the numerical solutions for the subspaces of the internal and scattering multipoles.

RESULTS AND DISCUSSION
We investigate the validity of the numerical solutions to the scattering problem calculated via three different methods; DSM using QR decomposition hereafter referred to as QR, TPM using solely Singular Value Decomposition (SVD) and also a TPM combination of both algorithms (QR+SVD).To provide a fair comparison between algorithms we limit the rank of the output for each method via regularization to be the same for all methods and study the effect of incrementing that limit.Simulations were run for two distinct particle types, a nanodisc of radius 400nm and depth 35nm and a nanorod of length 400nm and diameter 35nm, with rounded edges.Other than geometry, the two particles differ in the type of sources used to represent the fields.For the rod multipole sources are distributed along the symmetry axis in the real space whereas for the disc the sources are located in the complex space effectively making them ring sources distributed concentrically along the particle radius [14].The particles were illuminated by a near field source of wavelength 720nm comprised of a combination of electric and magnetic point dipoles located 50nm from the particle surface.The approximate locations of the near field source, moved to obtain average values for some tests by using different locations and polarisations, are highlighted in Fig. 1.Firstly, we compare the convergence of the solutions by plotting the differential scattering cross sections (DSCS), the angular variation of the electric field intensity in the far field [13], of each of the three methods by increasing the rank from an effective minimum.These results were obtained by calculating the light scattered by the excited particles into the far field along the generatrix line, 0   , sampling θ at equal intervals between the poles of the particle's symmetry axis at 0 and π, shown in Fig. 2. We observe that for minimal rank there is an obvious advantage to the TPM methods, which while not fully converged show the main features of the spectrum at the correct angles.The QR solution however, for both the rod and disc particles, fails to even approximately produce these features of the solution when the rank is minimal.As the rank is increased both TPM methods converge more rapidly than the pure QR solution which requires the maximum rank achievable with the TPM methods to show full convergence, for the disc, and approximate convergence for the rod.Note that with these particular source configurations the upper bound on the rank obtainable for SVD and QR+SVD when no limit is imposed is almost half that observed for the QR algorithm.
As we are solving the scattering problem by using a surface method we test the numerical validity of the surface fields primarily through the fractional L 2 surface error, where the norm of the surface field residual is calculated in terms of each input field [10], Where the tangential components of the incident, internal and scattered fields projected onto the particle surface are represented by 0 f , i f , and s f respectively.Due to the cylindrical symmetry of the particles we are able to separate components of the fields according to their angular phase dependence

 
exp im , where m is the component of the optical angular momentumalong the symmetry axis.Convergence of the fractional L 2 error test was checked with varying sized grids, ranging from 6000 to 12000 points, and a small oscillation of the computed value was observed but with a maximum deviation of ca.6% of the results shown in Fig. 3.As we would expect from the results of the DSCS, SVD and QR+SVD perform much better with minimal rank and produce an acceptably small (less than 4%) L 2 error on the disc for 2 m  , as shown in Fig. 3.As we increase the rank, we find that QR catches up with the TPM methods and that we have an error of less than 10% of the incident field for 6 m  .Increasing the rank further for the pure QR case does produce an even lower L 2 error and it does begin to outperform the other methods at high m .The rod particle is much easier to integrate and we observe a very low residual up to 7 m  however due to the limited radius of the particle only the fields for 2 m  are non-negligible.For this type of particle, QR must retain a much higher rank of the composite matrices to perform as well as the TPM methods and so with limited rank the L 2 error fluctuates strongly as the near field source is scanned across the particle.comprised of a combination of electric and magnetic point dipoles located 50nm from the particle surface.The approximate locations of the near field source, moved to obtain average values for some tests by using different locations and polarisations, are highlighted in Fig. 1.Firstly, we compare the convergence of the solutions by plotting the differential scattering cross sections (DSCS), the angular variation of the electric field intensity in the far field [13], of each of the three methods by increasing the rank from an effective minimum.These results were obtained by calculating the light scattered by the excited particles into the far field along the generatrix line, 0   , sampling θ at equal intervals between the poles of the particle's symmetry axis at 0 and π, shown in Fig. 2. We observe that for minimal rank there is an obvious advantage to the TPM methods, which while not fully converged show the main features of the spectrum at the correct angles.The QR solution however, for both the rod and disc particles, fails to even approximately produce these features of the solution when the rank is minimal.As the rank is increased both TPM methods converge more rapidly than the pure QR solution which requires the maximum rank achievable with the TPM methods to show full convergence, for the disc, and approximate convergence for the rod.Note that with these particular source configurations the upper bound on the rank obtainable for SVD and QR+SVD when no limit is imposed is almost half that observed for the QR algorithm.
As we are solving the scattering problem by using a surface method we test the numerical validity of the surface fields primarily through the fractional L 2 surface error, where the norm of the surface field residual is calculated in terms of each input field [10], Where the tangential components of the incident, internal and scattered fields projected onto the particle surface are represented by 0 f , i f , and s f respectively.Due to the cylindrical symmetry of the particles we are able to separate components of the fields according to their angular phase dependence

 
exp im , where m is the component of the optical angular momentumalong the symmetry axis.Convergence of the fractional L 2 error test was checked with varying sized grids, ranging from 6000 to 12000 points, and a small oscillation of the computed value was observed but with a maximum deviation of ca.6% of the results shown in Fig. 3.As we would expect from the results of the DSCS, SVD and QR+SVD perform much better with minimal rank and produce an acceptably small (less than 4%) L 2 error on the disc for 2 m  , as shown in Fig. 3.As we increase the rank, we find that QR catches up with the TPM methods and that we have an error of less than 10% of the incident field for 6 m  .Increasing the rank further for the pure QR case does produce an even lower L 2 error and it does begin to outperform the other methods at high m .The rod particle is much easier to integrate and we observe a very low residual up to 7 m  however due to the limited radius of the particle only the fields for 2 m  are non-negligible.For this type of particle, QR must retain a much higher rank of the composite matrices to perform as well as the TPM methods and so with limited rank the L 2 error fluctuates strongly as the near field source is scanned across the particle.comprised of a combination of electric and magnetic point dipoles located 50nm from the particle surface.The approximate locations of the near field source, moved to obtain average values for some tests by using different locations and polarisations, are highlighted in Fig. 1.Firstly, we compare the convergence of the solutions by plotting the differential scattering cross sections (DSCS), the angular variation of the electric field intensity in the far field [13], of each of the three methods by increasing the rank from an effective minimum.These results were obtained by calculating the light scattered by the excited particles into the far field along the generatrix line, 0   , sampling θ at equal intervals between the poles of the particle's symmetry axis at 0 and π, shown in Fig. 2. We observe that for minimal rank there is an obvious advantage to the TPM methods, which while not fully converged show the main features of the spectrum at the correct angles.The QR solution however, for both the rod and disc particles, fails to even approximately produce these features of the solution when the rank is minimal.As the rank is increased both TPM methods converge more rapidly than the pure QR solution which requires the maximum rank achievable with the TPM methods to show full convergence, for the disc, and approximate convergence for the rod.Note that with these particular source configurations the upper bound on the rank obtainable for SVD and QR+SVD when no limit is imposed is almost half that observed for the QR algorithm.
As we are solving the scattering problem by using a surface method we test the numerical validity of the surface fields primarily through the fractional L 2 surface error, where the norm of the surface field residual is calculated in terms of each input field [10], Where the tangential components of the incident, internal and scattered fields projected onto the particle surface are represented by 0 f , i f , and s f respectively.Due to the cylindrical symmetry of the particles we are able to separate components of the fields according to their angular phase dependence

 
exp im , where m is the component of the optical angular momentumalong the symmetry axis.Convergence of the fractional L 2 error test was checked with varying sized grids, ranging from 6000 to 12000 points, and a small oscillation of the computed value was observed but with a maximum deviation of ca.6% of the results shown in Fig. 3.As we would expect from the results of the DSCS, SVD and QR+SVD perform much better with minimal rank and produce an acceptably small (less than 4%) L 2 error on the disc for 2 m  , as shown in Fig. 3.As we increase the rank, we find that QR catches up with the TPM methods and that we have an error of less than 10% of the incident field for 6 m  .Increasing the rank further for the pure QR case does produce an even lower L 2 error and it does begin to outperform the other methods at high m .The rod particle is much easier to integrate and we observe a very low residual up to 7 m  however due to the limited radius of the particle only the fields for 2 m  are non-negligible.For this type of particle, QR must retain a much higher rank of the composite matrices to perform as well as the TPM methods and so with limited rank the L 2 error fluctuates strongly as the near field source is scanned across the particle.The error in the propagation of the scattered fields can be determined by comparing the integral of the Poynting vectors on the surface and also at infinity [13].As with the L 2 error this test was checked for convergence by varying the number of grid points and the maximum deviation from the results reported in Fig. 3 was ca.0.05%.This flux ratio gives an indication as to the quality of the scattered field produced by evaluating the error in the propagation of the special functions and for the disc particle all methods perform similarly despite the large difference in DSCS results and L 2 error, particularly with minimal rank.For the rod particle however, we only expect valid results from the flux ratio where the scattered field is non-negligible.Again, the QR method cannot compete with the TPM methods particularly with minimal rank however it does eventually perform as well when the rank is steadily increased beyond what is shown in Fig 3.
A further test of the validity of the calculated fields is the introduction of the Stratton-Chu test at infinity, where the scattered field can be compared with an exact solution of Maxwell's equations, and the Internal field should be exactly zero [13]; e e e r e dS r C  We use the asymptotic form of the multipole sources [13], A E , to evaluate the field calculated using from the scattered field as   . The convergence of the Stratton-Chu test for the scattered field was again tested by varyingthe grid size and showed a maximum discrepancy of the order 1E-6 when compared with the asymptotic values.A much larger fluctuation in results at wide angles was observed for the test of the internal field, which uses a different kernel to the scattered.In Fig. 4 we again see evidence that the disc particle is particularly difficult to integrate when compared to the rod which gives excellent results over all the m channels for each method.For the disc only up to 4 m  at low rank and 6 m  at higher rank give a value close to the value of the asymptotic sources for the scattered field, a result similar to that observed for the flux ratio.In fact the Stratton-Chu test proves to be a more stringent test of the scattered field than the flux ratio.The Stratton-Chu test on the internal field of the disc appears to indicate that there is a problem with the field.The field is evaluated along a line in the far field from zero to  , and is non-zero for wide angles around / 2  .This is due to the fact that the grid along the curved edge of the particle is particularly difficult to integrate and has not converged for this number of grid points.To indicate this more clearly, in Fig. 4c we plot also the Stratton-Chu test for the internal field, calculated using QR+SVD on the surface of the disc, expanded out onto a sphere with radius equal to that of the disc, where it convincingly passes this test for all m channels at all angles.The QR algorithm appears to perform better for the internal field than the other methods but this is due to the fact that it assigns the sources significantly smaller amplitudes when solving for the fields at the particle boundary, as such the values calculated at infinity also appear smaller.To the best of our knowledge the Stratton-Chu test is the first procedure developed to evaluate the quality of the internal field.
We have observed that for low rank that there is a clear advantage to using a method which splits the space into two subspaces not only for the extra information about the system which is obtained but also for the accuracy in the calculations performed.There is also another advantage to using the TPM methods, due to the sequential way in which the surface fields are calculated using SVD they can be written out to be used again for a different excitation of the same particle.While, for the initial calculation QR proves to be slightly quicker, as shown in Table 1., for multiple calculations SVD and QR+SVD need only calculate the fields once and the subsequent calculations are significantly faster, by a factor of ~5 for QR+SVD and ~7 for SVD.Table 1.Total computational time for a full solution of the scattering problem for the disc particle using an AMD Opteron Processor 6344 2.6 GHz system averaged over 5 runs.For QR+SVD and pure SVD we highlight the time taken for the initial calculation and also subsequent calculations for the same particle where the fields are read back in.

CONCLUSION
In conclusion we have shown that the TPM -by separately considering the internal and scattered subspaces of the electromagnetic fields -hasa faster convergence with the number offunctions used to expand these solutions than the Discrete Source Method implemented through QR decomposition.In addition to evaluating thefractional error in the calculated solutions at the particle surface, we have also quantitatively tested the quality of the numerical evaluation of the transport of energy away from the particle.However this does not reveal errors in the field inside of the particle, hencewe havedemonstrated that the Stratton-Chu relations offer an excellent metric for validating the reliability of both the internal and scattered fields, providing the first general test for the internal field.For metallic particles these tests reveal that the internal field is the most significant source of numerical error in these calculations.We have also shown that for particles with large aspect ratio, such as those considered here, the accuracy of the surface quadrature is extremely important, this suggests that integration via adaptive grids may be beneficial in improving the accuracy of the calculations.

Fig. 1 .
Fig. 1.Sampling points of near field excitation for a rounded gold nanodisc and nanorod.The red points indicate the approximate location of the near field source as it was scanned over the gold nanoparticles at a height of 50nm, for a rod wi th dimensions (l=400nm, d=35nm) and a disc (d=800nm, z=35nm).The blue circles indicate the location of the near field source for the differential scattering cross section and Stratton-Chu measurements.There are 15 sampling points for each particle (the centre of the disc was sampled with 3 different polarisations.)All of the following simulations were performed using these particles.

Fig. 1 .
Fig. 1.Sampling points of near field excitation for a rounded gold nanodisc and nanorod.The red points indicate the approximate location of the near field source as it was scanned over the gold nanoparticles at a height of 50nm, for a rod wi th dimensions (l=400nm, d=35nm) and a disc (d=800nm, z=35nm).The blue circles indicate the location of the near field source for the differential scattering cross section and Stratton-Chu measurements.There are 15 sampling points for each particle (the centre of the disc was sampled with 3 different polarisations.)All of the following simulations were performed using these particles.

Fig. 1 .
Fig. 1.Sampling points of near field excitation for a rounded gold nanodisc and nanorod.The red points indicate the approximate location of the near field source as it was scanned over the gold nanoparticles at a height of 50nm, for a rod wi th dimensions (l=400nm, d=35nm) and a disc (d=800nm, z=35nm).The blue circles indicate the location of the near field source for the differential scattering cross section and Stratton-Chu measurements.There are 15 sampling points for each particle (the centre of the disc was sampled with 3 different polarisations.)All of the following simulations were performed using these particles.

Fig. 2 .
Fig. 2.Convergence of differential scattering cross sections (DSCS) along the generatrix line with increasing rank.The DSCS, in arbitrary units, for the three different algorithms plotted against the far field angleθ, varied incrementally between 0 and π between the poles of the particle's symmetry axis showing convergence with increasing rank of the solution matrices for a (a) disc and (b) rod.

Fig. 3 ..
Fig. 3. Fractional L2 surface error and flux ratio between the particle surface and the far field for a gold (a)-(d) nanodisc and (e)-(h) nanorod for the three solution methods plotted against the absolute value of the index of optical angular momentum, m, of the field.All points shown are the averages  of the 15 sampling points indicated in Fig. 1 and the error bars show the standard deviation

Fig. 4 .Fig. 4 .Fig. 4 .
Fig. 4. Evaluation of the internal and scattered fields using the Stratton-Chu test at infinity for the disc (a)-(d) and rod (e)-(h) particles.Plotted, for each of the three algorithms, are the average fields, ES-C, calculated at 15 different points along θ from 0 toπ,φ=0 for different scattering channels m.The scattered fields are compared with exact solutions of the Maxwell equations EA.Fig 4c) shows an additional plot, QR+SVD*, where the internal field calculated on the surface of the disc was expanded onto the surface of a sphere of equal radius are the tangential components of the surface electric and magnetic fields, r e is the unit vector towards the evaluation point in the far field, r  is the point on the surface S and s D is the scattering domain.