Recursive T matrix algorithm for resonant multiple scattering: Applications to localized plasmon excitations
aa r X i v : . [ phy s i c s . op ti c s ] M a y Recursive T matrix algorithm for resonant multiple scattering :Applications to localized plasmon excitations Brian Stout and J.C. Auger and Alexis Devilez Case 161, Institut Fresnel, Facult´e des Sciences et Techniques de St. J´erˆome13397 Marseille cedex 20, France Center for Laser Diagnostics, Dept. of Applied Physics andYale University, New Haven, CT 06520 USA ∗ (Dated: July 2008) A matrix balanced version of the Recursive Centered T Matrix Al-gorithm (RCTMA) applicable to systems possessing resonant inter-particle couplings is presented. Possible domains of application in-clude systems containing interacting localized plasmon resonances,surface resonances, and photonic jet phenomena. This method is ofparticular interest when considering modifications to complex sys-tems. The numerical accuracy of this technique is demonstratedin a study of particles with strongly interacting localized plasmonresonances.c (cid:13)
OCIS codes:
Keywords: Multiple scattering, T matrices, morphology dependent resonances,plasmon resonances
1. Introduction
It has been well established that certain kinds of recursive T matrix algorithms (knownas RCTMA) are numerically stable and can be used to solve the Foldy-Lax multiple-scattering equations for particles exhibiting “modest” inter-particle couplings. By “modestcouplings”, we refer to situations in which the order of orbital number of the Vector Spherical ∗ Electronic address: [email protected], [email protected], [email protected]
Wave Functions (VSWFs) necessary to describe the field scattered by each particle in anaggregate of particles are not too much larger than that necessary for describing isolatedparticles. The “modest coupling” criteria apply to a host of multiple scattering situations,including systems of dielectric particles comparable in size to the wavelength and for mostpacking fractions including dense packing. The modest coupling criteria can also apply tometallic particles under certain conditions.Like any multiple scattering technique not employing matrix balancing, the RCTMA canencounter numerical difficulties in certain extreme situations of strongly coupled resonantphenomenon. In this work, we present a matrix balanced form of the Recursive Centered T Matrix Algorithm (or RCTMA) that can readily be employed even in the presence of strong( i.e. resonant) inter-particle couplings. The rather extreme situation of “strong couplings”studied here will generally require carefully micro-scaled engineered systems where high Q -factor resonances can occur for particles illuminated in isolation, and in which the parti-cles are sufficiently closely spaced that neighboring particles modify the resonance responseproperties. Examples of strong inter-particle couplings can be found in particles exhibitingplasmon resonances, surface resonances, or even photonic jet phenomenon.In section 2, the notation is introduced in a brief review of the relevant multiple-scatteringtheory. Section 3 describes an analytic matrix balancing procedure used to ‘well-condition’the multiple scattering system of equations. A matrix balanced RCTMA is derived in section4. Essential formulas for applications are summarized in section 5. Their applications arethen demonstrated by applying matrix balanced RCTMA calculations to study systems ofinteracting localized plasmon excitations. Some known and novel aspects of interactinglocalized plasmon excitations are presented herein.
2. Multiple-scattering theory - VSWF approach
Let us consider an arbitrary incident electromagnetic field incident on a collection of three-dimensional particles (as shown in fig.1). The particles are considered as ‘individual’ scat-terers if they can be placed in a circumscribing sphere lying entirely within the homogeneousmedium (actually this constraint can frequently be relaxed, cf. ).The electromagnetic field incident on an N -particle system, E i , is developed in termsof the transverse regular VSWFs developed about some point O arbitrarily chosen as the x x x N O x yz x R R R N R E inc Fig. 1.
Schematic of an field incident on a collection of scatterers centered on x , x , ..., x N . The radii ofthe respective circumscribing spheres are denoted R , R , ..., R N . system ‘origin’: E i ( r ) = E ∞ X n =1 n X m = − n { Rg [ M nm ( k r )] a ,n,m + Rg [ N nm ( k r )] a ,n,m } = E X q =1 ∞ X p =1 Rg [ Ψ q,p ( k r )] a q,p ≡ E Rg (cid:2) Ψ t ( k r ) (cid:3) a (1)where E is a real parameter determining the incident field amplitude. Eq.(1), introduces acondensed notation for the VSWFs, M nm and N nm : Ψ ,p ( k r ) ≡ M n,m ( k r ) and Ψ ,p ( k r ) ≡ N n,m ( k r ). The notation Rg [ ] stands for “the regular part of” and distinguishes theseregular VSWFs from the “irregular” scattered VSWFs (cf. appendix A). In the second lineof eq.(1), the two subscripts ( n, m ) are replaced by a single subscript p defined such that p ( n, m ) ≡ n ( n + 1) − m and has the inverse relations : n ( p ) = Int √ p m ( p ) = − p + n ( n + 1) . (2)The last line of eq.(1), adopts the compact matrix notation allowing the suppression of thesummation symbols. The superscripted ( t ) stands for the transpose of a column ‘matrix’of composed of VSWFs into a row ‘matrix’ of these functions.For points external to all individual circumscribing spheres, the total field, E t ( r ) canbe written as the sum of the incident field, and a set of ‘individual’ scattered fields, E ( j )s ,centered respectively on each of the particle centers: E t ( r ) = E i ( r ) + N X j =1 E ( j )s ( r j )= E Rg (cid:2) Ψ t ( k r ) (cid:3) a + E N X j =1 Ψ t ( k r j ) f ( j ) N (3)where each scattered field, E ( j )s , is developed, with coefficients f ( j ) N , on a basis of outgoingVSWFs defined with respect to the associated particle center, denoted x j . The sphericalcoordinates relative to each scatterer are denoted r j ≡ r − x j .The crucial idea of Foldy-Lax multiple-scattering theory is that there exists an e xcitationfield, E ( j )exc ( r ), associated with each particle which is the superposition of the incident fieldand the field scattered by all the other particles in the system (excluding the field scatteredby the particle itself). From this definition, the excitation field of the j th particle can bewritten E ( j )exc ( r j ) ≡ E Rg (cid:2) Ψ t ( k r j ) (cid:3) e ( j ) N ≡ E i ( r ) + N X l =1 ,l = j E ( l )s ( r l )= E Rg (cid:2) Ψ t ( k r j ) (cid:3) " J ( j, a + N X l =1 ,l = j H ( j,l ) f ( l ) N (4)where e ( j ) N are the coefficients of the excitation field in a regular VSWF basis centered on the j th particle. In the second line of eq.(4), we have used the translation-addition theorem )and introduced the notation where J ( j, ≡ J ( k x j ) is a regular translation matrix and H ( j,l ) ≡ H [ k ( x j − x l )] is an irregular translation matrix. Analytical expressions for thematrix elements of J ( k x j ) and H ( k x j ) are given in refs. .The other key idea of multiple scattering theory is that the field scattered by the object, f ( j ) N , is obtained from the excitation field e ( j ) N via the 1-body T matrix, T ( j )1 , derived whenone considers the particle to be immersed in an infinite homogeneous medium. This relationis then expressed as f ( j ) N = T ( j )1 e ( j ) N (5)[The index 1 on the T ( j )1 indicates that this T matrix concerns an isolated particle, henceforthreferred to as a ‘1-body’ T matrix.] Employing eq.(5) in eq.(4), one obtains a Foldy-Lax setof equations for the excitation field coefficients : e ( j ) N = J ( j, a + N X l =1 ,l = j H ( j,l ) T ( l )1 e ( l ) N j = 1 , ..., N (6)For numerical applications where one is obliged to solve the equations on a truncatedVSWF basis, it is advantageous to work with a set of formally equivalent equations involvingthe scattering coefficients f ( j ) N . This set of equations is derived by multiplying each of eqs.(6)from the left by T ( j )1 and using eq.(5) to obtain f ( j ) N = T ( j )1 J ( j, a + T ( j )1 N X l =1 ,l = j H ( j,l ) f ( l ) N j = 1 , ..., N (7)In the RCTMA, one calculates the centered multiple scattering transition matrices, T ( j,k ) N ,which directly yield the scattered field coefficients in terms of the field incident on the systemthrough the expression f ( j ) N = N X k =1 T ( j,k ) N a ( k ) with a ( k ) ≡ J ( k, a (8)In this equation, we have introduced the column matrix a ( j ) which contains the coefficientsof the field incident on the entire system developed on a VSWF basis centered on the j th particle.
3. Basis set truncation and matrix balancing
Although the multiple scattering formulas of the previous section are expressed as matrixequations on VSWF basis sets of infinite dimension, the finite size of the scatterers natu-rally restricts the dimension of the dominate VSWF contributions. In order to discuss thisphenomenon analytically, we consider the case of spherical scatterers. For non-sphericalscatterers, the matrix balancing procedure described below should be applied to the circum-scribing spheres of the particles.The Mie solution for a sphere of radius R j immersed in a homogeneous host medium, canbe cast in the form of a 1-body T matrix that is diagonal in a VSWF basis centered on theparticle: h T ( j )1 i q,p ; q ′ ,p ′ = δ q,q ′ δ p,p ′ T ( j, n ( p ) , q ) (9)where the T ( j )1 ( n ( p ) , q ) correspond to the Mie coefficients and depend on q and n (cf. eq.(1)).With the objective of matrix balancing, it is helpful to express the Mie coefficients of thescatterers in terms of the Ricatti Bessel and Hankel functions, respectively ψ n ( z ) ≡ zj n ( z )and ξ n ( z ) ≡ zh n ( z ), and their logarithmic derivativesΦ n ( z ) ≡ ψ ′ n ( z ) ψ n ( z ) Ψ n ( z ) ≡ ξ ′ n ( z ) ξ n ( z ) (10)The T matrix elements of eq.(9) for a sphere of dielectric contrast ρ j ≡ k j /k can then becast in the convenient form : T ( j, n,
1) = ψ n ( kR j ) ξ n ( kR j ) µ j µ Φ n ( kR j ) − ρ j Φ n ( ρ j kR j ) ρ j Φ n ( ρ j kR j ) − µ j µ Ψ n ( kR j ) ≡ ψ n ( kR j ) ξ n ( kR j ) T ( j, n, T ( j, n,
2) = ψ n ( kR j ) ξ n ( kR j ) µ j µ Φ n (cid:0) ρ j kR j (cid:1) − ρ j Φ n ( kR j ) ρ j Ψ n ( kR j ) − µ j µ Φ n ( ρ j kR j ) ≡ ψ n ( kR j ) ξ n ( kR j ) T ( j, n,
2) (11)where k is the wavenumber in the external medium. The n ormalized T matrix coefficients, T ( j, n, q ), of eq.(MieT) contain a rich resonant structure. The ratios ψ n ( kR j ) /ξ n ( kR j )on the other hand have an exponentially decreasing behavior for large, n ≫ kR j as isdemonstrated in fig.2 for kR = 10. One can remark from figure 2 that | ψ n ( kR ) /ξ n ( kR ) | become quite small beyond n max = kR + 3 and its value at n = 14 is ∼ − . Althoughthese factors permit an appropriately truncated VSWF basis set to contain essentially allthe physical information necessary for accurate calculations, they also tend to produce ill-conditioned linear systems when one is obliged to enlarge the VSWF space far beyond ≈ kR + 3 in order to account for strong coupling phenomenon.A solution to the above problem is to ‘balance’ the matrix manipulations in section 4below by defining “normalized” scattering and incident coefficients: h f ( j ) i q,p ≡ ξ n ( p ) ( kR j ) (cid:2) f ( j ) (cid:3) q,p (cid:2) a ( j ) (cid:3) q,p ≡ ψ n ( p ) ( kR j ) (cid:2) a ( j ) (cid:3) q,p (12)For notational purposes, it is convenient to define diagonal matrices (cid:2) ψ ( j ) (cid:3) and (cid:2) ξ ( j ) (cid:3) withRicatti-Bessel functions along their diagonals, namely (cid:2) ψ ( j ) (cid:3) q ′ ,p ′ ; q,p ≡ δ q,q ′ δ p,p ′ ψ n ( p ) ( kR j ) and (cid:2) ξ ( j ) (cid:3) q ′ ,p ′ ; q,p ≡ δ q,q ′ δ p,p ′ ξ n ( p ) ( kR j ). This notation allows normalized or ‘balanced’ versions ofthe one-body and many-body T -matrices to be defined respectively as T ( j )1 ≡ (cid:2) ξ ( j ) (cid:3) T ( j )1 (cid:2) ψ ( j ) (cid:3) − and T ( j,k ) N ≡ (cid:2) ξ ( j ) (cid:3) T ( j,k ) N (cid:2) ψ ( k ) (cid:3) − (13) n Log ( ψ n ( ) / ξ n ( )) Fig. 2.
Plot of the spherical Bessel to Hankel function ratio, ψ n ( kR ) /ξ n ( kR ) occurring in the Miecoefficients when kR = 10. In terms of these normalized quantities, eq.(8) then reads f ( j ) N = N X k =1 T ( j,k ) N a ( k ) j = 1 , ..., N (14)In the next section, these ‘normalized’ T ( j )1 and T ( j,k ) N are used to derive a matrix balancedversion of the recursive T matrix algorithm.
4. Derivation of a matrix balanced recursive algorithm
In this section, we derive a matrix balanced version of the Recursive Centered T MatrixAlgorithm (RCTMA) using purely algebraic manipulations. The recursive algorithm can beinvoked once we have a solution for the T ( j,k ) N − matrices of a N ≥ T (1 , ≡ T (1)1 .One then considers an arbitrarily positioned particle being added to the system. Theexcitation field on a particle N added to the system can be expressed as the superpositionof three fields. The first contribution is simply the field incident on the system, the secondcontribution results from the scattering of the incident field by the N − N , and finally the third contribution comes from field scattered by theparticle N onto the N − N th particle as an excitationfield. Invoking the translation-addition theorem and eq.(8), these three contributions canbe expressed in matrix form as e ( N ) N = a ( N ) + N − X j,k =1 H ( N,j ) T ( j,k ) N − a ( k ) + N − X j,k =1 H ( N,j ) T ( j,k ) N − H ( k,N ) f ( N ) N (15)Defining now the normalized irregular translation matrices and excitation coefficientsrespectively as H ( j,k ) ≡ (cid:2) ψ ( j ) (cid:3) H ( j,k ) (cid:2) ξ ( k ) (cid:3) − and e ( j ) N ≡ (cid:2) ψ ( j ) (cid:3) e ( j ) N (16)the normalized form of eq.(15) can be written e ( N ) N = a ( N ) + N − X j,k =1 H ( N,j ) T ( j,k ) N − a ( k ) + N − X j,k =1 H ( N,j ) T ( j,k ) N − H ( k,N ) f ( N ) N (17)where we also used the definitions in eqs.(12) and (13).Recalling that the excitation field is linked to the scattered field by the 1-body T -matrixvia eq.(5), and invoking the definitions of eq.(12) and (16) we can write e ( j ) N = h T ( j )1 i − f ( j ) N (18)Employing this relation for particle N on the LHS of eq.(17) and rearranging we obtain (h T ( N )1 i − − N − X j,k =1 H ( N,j ) T ( j,k ) N − H ( k,N ) ) f ( N ) N = a ( N ) + N − X j,k =1 H ( N,j ) T ( j,k ) N − a ( k ) . (19)We now take the normalized T ( N,N ) N matrix to be expressed as T ( N,N ) N = (h T ( N )1 i − − N − X j,k =1 H ( N,j ) T ( j,k ) N − H ( k,N ) ) − (20)With this assignment, we multiply both sides of eq.(19) by T ( N,N ) N and obtain an expressionconsistent with equation (14): f ( N ) N = T ( N,N ) N a ( N ) + T ( N,N ) N N − X j,l =1 H ( N,j ) T ( j,k ) N − a ( k ) = T ( N,N ) N a ( N ) + N − X l =1 T ( N,k ) N a ( k ) = N X l =1 T ( N,k ) N a ( k ) (21)where we have assigned the matrix T ( N,k ) N , k = N as T ( N,k ) N = T ( N,N ) N N − X j =1 H ( N,j ) T ( j,k ) N − (22)One completes the description of the scattering by the system by remarking that the fieldscattered by the other particles in the system are the superposition of the field that wouldbe scattered by the N − N th particle, plus the field scatteredfrom the N − N th particle.Using again the translation-addition theorem, the field coefficients of f ( j ) N can in turn beexpressed in a form consistent with equation (14) as f ( j ) N = N − X l = k T ( j,k ) N − a ( k ) + N − X k =1 T ( j,k ) N − H ( k,N ) f ( N ) N = N − X l =1 T ( j,k ) N − a ( k ) + N − X k =1 T ( j,k ) N − H ( k,N ) T ( N,N ) N a ( N ) + N − X l =1 N − X k =1 T ( j,l ) N − H ( l,N ) T ( N,k ) N a ( k ) = T ( j,N ) N a ( N ) + N − X k =1 T ( j,k ) N a ( k ) = N X k =1 T ( j,k ) N a ( k ) (23)where we invoked eq.(21). In the second and third lines we have defined the T ( j,N ) N and T ( j,l ) N matrices such that T ( j,N ) N = N − X k =1 T ( j,k ) N − H ( k,N ) T ( N,N ) N (24a) T ( j,k ) N = T ( j,k ) N − + N − X l =1 T ( j,l ) N − H ( l,N ) T ( N,k ) N (24b)At this point, all the T ( j,k ) N matrices have been obtained and the matrix manipulationsin eqs.(20), (22) and (24) can then be repeated to add as many particles to the system asdesired. A. Relationship with system matrix inversions
Although the recursive algorithm is quite efficient for systems with relatively small numbersof particles, for systems with many particles, one may prefer to try and solve an entire N -particle system directly. A balanced linear system for the entire system corresponding to0our recursive algorithm can be obtained by applying the relation of eq.(5) to the left handside of eq.(7), then multiplying both sides of the resulting equations by the (cid:2) ψ ( j ) (cid:3) matrixand finally rearranging to obtain a system of balanced linear equations for the unknownscattering coefficients: h T ( j )1 i − f ( j ) N − N X k =1 ,k = j H ( j,k ) f ( k ) N = a ( j ) j = 1 , ..., N (25)where we used eqs.(12) and (13). The system of linear equations in eq.(25) can in principlebe directly solved by inverting the balanced system matrix: f (1) N f (2) N ... f ( N ) N = h T (1)1 i − − H (1 , · · · − H (1 ,N ) − H (2 , h T (2)1 i − · · · − H (2 ,N ) ... ... . . . ... − H ( N, − H ( N, · · · h T ( N )1 i − − a (1) a (2) ... a ( N ) (26)Once we have inverted this system, one can associate each block with a corresponding T ( j,k ) N matrix as shown below as f (1) N f (2) N ... f ( N ) N = T (1 , N T (1 , N · · · T (1 ,N ) N T (2 , N T (2 , N · · · T (2 ,N ) N ... ... . . . ... T ( N, N T ( N, N · · · T ( N,N ) N a (1) a (2) ... a ( N ) (27)which is the same form as the desired solutions given in eq.(14).
5. Summary and applications to localized plasmon excitations
In this section, we will apply the RCTMA to solve systems exhibiting strong interactionsbetween localized plasma resonances. We begin this section by summarizing the balancedrecursive algorithm. We then recall some useful formulas for extracting physical quanti-ties from the T matrix. Finally, we carry out some illustrative calculations for stronglyinteracting systems. A. Summary of the balanced RCTMA algorithm
In order to implement the RCTMA, one must first solve the 1-body T -matrices, T (1)1 , T (2)1 , ... , T ( N tot )1 , for all the particles in the system. Normalized versions of the 1-body T matrices1and the irregular translation matrices , H ( j,k ) , are then calculated via T ( j )1 ≡ (cid:2) ξ ( j ) (cid:3) T ( j )1 (cid:2) ψ ( j ) (cid:3) − H ( j,k ) ≡ (cid:2) ψ ( j ) (cid:3) H ( j,k ) (cid:2) ξ ( k ) (cid:3) − (28)where the diagonal matrices (cid:2) ψ ( j ) (cid:3) q,q ′ ,p,p ′ = δ q,q ′ δ p,p ′ ψ n ( p ) ( kR j ) and (cid:2) ξ ( j ) (cid:3) q,q ′ ,p,p ′ = δ q,q ′ δ p,p ′ ξ n ( p ) ( kR j ) respectively have Ricatti-Bessel and Ricatti-Hankel functions on the di-agonal. ( R j being the radius of the circumscribing sphere of the jth scatterer).The balanced recursive algorithm is that the solution for the T ( N,N ) N matrix is obtainedfrom the T matrices of the N − T ( j,k ) N − , via the matrix inversion in eq.(20). Allthe other matrices T ( j,k ) N with j = N or k = N are then obtained via matrix multiplicationsand additions via equations (22) and (24). This process is then repeated as many times asdesired. B. Physical quantities
When the incident field is a plane wave, it is convenient to express physical quantities interms of cross sections. Appealing to the far-field approximation of the field, the extinctionand scattering cross sections of clusters of N objects can be respectively expressed σ ext = − k Re " N X k =1 a ( j ) , † f ( j ) N and σ scat = 1 k N X j,k =1 f ( j ) , † N J ( j,k ) f ( k ) N (29)It is also possible to produce analytical expression for local field quantities like indi-vidual absorption cross sections. For lossy scatterers in a lossless host medium, one canobtain individual particle absorption cross sections by integrating the Poynting vector on acircumscribing sphere surrounding the particle to obtain the formula as σ ( j )a = − k Re n f ( j ) , † N e ( j ) N o − k (cid:12)(cid:12)(cid:12) f ( j ) N (cid:12)(cid:12)(cid:12) (30)In an analogous fashion, optical forces on the particles can be calculated by integratingthe Maxwell tensor on a circumscribing sphere surrounding the particle . It is frequentlyconvenient to characterize the optical force by vector ‘cross sections’, −→ σ opt , defined suchthat the time averaged optical force on particles immersed in a liquid dielectric of refractionindex n med can be expressed as F opt = k S inc k n med c −→ σ opt (31)2where k S inc k = (cid:13)(cid:13) Re { E ∗ inc × H inc } (cid:13)(cid:13) is the incident irradiance. The binding force and itsassociated cross section, σ b , between two particles separated by a relative position vector r pos ≡ r − r , can be defined as F b ≡
12 ( F − F ) · b r pos ≡ k S inc k n med c σ b (32) C. Interacting localized plasmon excitations
For those conductors, such as the noble metals, that support surface plasmon resonances,one can usually observe localized plasmon resonances in sufficiently small particles. Theseresonances are typically dominated by absorption if the particles are sufficiently small withrespect to the incident wavelength and by scattering for larger particles. We chose to studysilver spheres 50 nm in diameter immersed in air (for which both scattering and absorptionare non-negligible).The study is carried out for wavelengths ranging from the near ultra-violet through thevisible (300 to 850 nm). We ignore the relatively modest finite size corrections to damping and simply adopt the bulk dielectric constant of silver from ref. and extrapolate betweenthe experimental values. The extinction, scattering and absorption cross sections for theseparticles are readily obtained from Mie theory and are displayed in figure 3 as a function offrequency. These spheres are quite small with respect to visible wavelengths, (size parametersin the 300 ↔
800 nm wavelength range go through kR = 0 . ↔ .
20) and the isolated particlecross sections are obtained to high precision with n max = 4. One can also see from figure 3that the strength of the plasmon resonance for these particles is about half due to absorptionand about half due to scattering.One of the principal sources of interest of the localized plasmon resonances is their ca-pacity to produce large field enhancements in regions much smaller than the incident fieldwavelength. This property is demonstrated in fig.4a) with a 2D and 1D plot of the electricfield intensity in and near an isolated 50 nm diameter silver sphere illuminated near itsresonance peak ( λ = 365 nm with N Ag = 0 .
077 + 1 . i ). The plots in Fig.4 are performedin a plane containing the center of the sphere and perpendicular to k inc (the polarizationlies along the horizontal axis). The dimensionless extinction and scattering ‘efficiencies’, Q = σ/ ( πR ), at this frequency are are respectively Q ext = 14 .
48 and Q scat = 6 .
200 300 400 500 600051015 λ (nm) σ / ( π R ) Q ext Q abs Q scat Fig. 3.
Total cross section ‘efficiencies’, Q ≡ σ/ ( πR ) for an isolated 50 nm diameter sphere. the T matrix calculated by RCTMA contains information for arbitrary incident fields, westudy the physically interesting case of a plane wave perpendicular to the axis separatingthe particles. As is widely known, the response then depends strongly on the polarizationof the incident light. In figures 5a) and 5c), the extinction and scattering cross sectionsper particle are presented when the polarization is respectively perpendicular and parallelto the symmetry axis. From figure 5c), one sees that the cross section for the parallel toaxis polarization presents a two sphere coupled resonance that is strongly red-shifted withrespect to the isolated particle resonance. The polarization perpendicular to this axis onthe other hand presents only slight modifications with respect to an isolated sphere.The optical binding force cross sections for these same polarizations are respectively plot-ted in figures 5b) and 5d). While the binding force for the polarization perpendicular to theparticle axis (cf.5c)) is slightly repulsive, the force for polarization parallel to the resonancecan be highly attractive with the dimensionless | Q b | attaining amplitudes of three ordersof magnitude. There has already been experimental and theoretical evidence supportingthe existence of optical force couplings in particles with plasmon excitations although suchhigh precision calculations at such small separations seems not to have been presented beforenow.This dimer system dramatically illustrates the ‘strong’ coupling category since correctcalculations require that the VSWF space be enlarged far beyond the predominantly dipolarresponse characterizing the particles in isolation. The normalized cross sections per particleare given in the table 1 for different values of the VSWF space truncation. Although it was4Fig. 4. Electric field intensity || E t || / || E inc || in an isolated 50 nm diameter sphere ( λ = 365 nm, N Ag = 0 .
077 + 1 . i ). In a) is presented a 2D (hot) plot of the electric field intensity in a plane perpendicularto the wavevector and containing the origin of the sphere (the horizontal axis lies along the polarizationdirection). Fig b) is a 1 D plot of the field intensity along the line in this plane containing the direction ofelectric field polarization. necessary to go to ∼ n max = 30 to achieve 4 digit precision in all the cross sections, the tableindicates that results were already quite good at n max = 20.A base 10 logarithmic intensity field map for a two silver sphere dimer illuminated withlight polarized along the symmetry axis (frequency near the coupled sphere resonance maxi-mum ( λ = 467 nm and N Ag = 0 .
048 + 2 . i ) is presented in figs.6a) and figs.6b) which arerespectively a 2D plot (in the same plane as figure 4) and a 1D plot along the symmetry axis.The size parameter of the individual spheres is kR = 0 .
34 and the isolated cross sections atthis frequency are Q ext = 0 .
136 and Q scat = 0 .
300 400 5000246810 a) ⊥ polarization σ / ( π R ) Q ext Q scat
300 400 50005 b) λ (nm) σ b / ( π R ) Q b
400 500 600 700051015 || polarization c) Q ext Q scat
400 500 600 700−6000−30000 d) λ (nm) Q b Fig. 5.
Dimensionless cross section ‘efficiencies’ per particle, Q = σ/ (2 πR ) and binding force ‘efficiencies’for a dimer of 50 nm diameter spheres (1 nm separation). In a) and b) the polarization is perpendicular tothe symmetry axis and in c) and d) it is parallel to the symmetry axis. n max Q ext / Q scat / Q b -417 -3639 -5530 -5918 -6000 -6015 -6018 -6018 Table 1.
Dimensionless cross section ‘efficiencies’ per particle in function of the VSWF truncation, n max . Q ext = σ ext / (2 πR ), Q scat = σ scat / (2 πR ) and Q b = σ b / ( πR ). The system is a dimer composed of D = 50 nm diameter silver spheres (1 nm separation) at ( λ = 467 nm and N Ag = 0 .
048 + 2 . i ) An important word of caution should be made at this point. Although 1 nm separationmay appear to ‘nearly’ touching, the coupled resonance is in fact quite sensitive to exactseparation details when resonant particles are so closely separated. For example, at a sep-aration distance of 0 . λ ≃
516 nm as compared with λ ≃
467 nm for a 1 nm separation, and the multipoleorder has to be pushed to n max ≃
50 to achieve four digit accuracy in the cross sections.Nanometer scale separations are not necessarily theoretical idealizations however as recent6Fig. 6.
Logarithmic scale plots of the field intensity for a two sphere dimer with ( λ = 467 nm and incidentlight polarized along the sphere axis ( N Ag = 0 .
048 + 2 . i ). In a) is a 2D plot in the plane containing thecenters of the spheres and the polarization vector while b) is a 1D Logarithmic plot along the symmetryaxis of the spheres. Figures. c) and d) are the same as a) and b) respectively but for a 5 sphere chain ofspheres at its resonance maximum ( λ = 561 nm and N Ag = 0 . . i ). (cf. figure 7). . Nevertheless, in applications likeDNA separators, one may well have to consider the strong optical forces between theseparticles on account of the exceptionally strong attractive optical forces efficiencies of theseresonances. For instance, the binding force efficiency at 0 . mathrmnm separation was calculated at Q b = − λ ≃ . Q b = − λ ≃
467 nm) . The question of perfect spheres exactly in contacthowever seems untenable from an experimental standpoint and quite difficult from a the-oretical standpoint on account of the singular behavior of the contact point. Theoretical‘separations’ of 1 ˚ A for instance require multipole truncations of the order of n max & Q b , which is the binding optical force between outermost spheres andtheir nearest neighbor and Q b , which is the binding force between central sphere and each ofits nearest neighbors. It is interesting to remark that addition of other spheres in the chaindramatically lessens the strong binding force interactions between spheres even though thefields between the spheres (cf. fig.5) can still be almost as high as the dimer case.We remark that the interactions have continued to red-shift and widen the coupled“chain” resonance. This chain resonance peaks at ≈
561 nm and N Ag ≃ . . i ( Q ext / .
416 and Q scat / . n max = 20 sincecalculations at n max = 24 produced negligible differences on this scale.Despite the dominance of scattering, a considerable amount of absorption is still presentin the 5 sphere chain. Furthermore, from the field maps in figures 6c) and 6d), one can seethat the central sphere has the highest field internal field intensities, and one consequentlyexpects increased absorption in the central sphere. This supposition can readily be confirmedquantitatively by using eq.(30) to calculate the absorption in each individual sphere. The8
300 400 50002468 a) ⊥ polarization σ / ( π R ) Q ext Q scat
300 400 500−10−505 b) λ (nm) σ b / ( π R ) Q b,1 Q b,2
400 500 600 700051015 || polarization c) Q ext Q scat
400 500 600 700−800−4000400 d) λ (nm) Q b,1 Q b,2 Fig. 7.
Total cross section and binding ‘efficiencies’ for a line of 5 ‘touching’ silver spheres 50 nm in diameter(1 nm separation). In a) and b) the polarization is perpendicular to the symmetry axis and in c) and d) itis parallel to the symmetry axis. results are given in table 2. Q a , Q a , Q a , Q a , Q a , . .
333 3 .
030 2 .
333 0 . Table 2.
Individual absorption efficiencies Q a ,j ≡ σ a ,j / ( πR ) in a five sphere chain at λ = 561 nm and N Ag = 0 . . i We conclude this section with some calculations systems concerning larger chains ofparticles. One can remark that the chain coupled resonance continued to red-shift andwiden when passing from the dimer to the five particle chain. Results for the extinction andscattering cross sections chains of 10 and 20 sphere chains are presented in figure 8 for thesame polarizations and incident directions as considered previously.One readily sees that ultraviolet and perpendicular responses per particle seem to havestabilized for large chains. The collective chain response on the other hand continued to9broaden and slightly red shift as one passed from 10 to 20 sphere chains and it is an interestingpoint for future studies to examine the evolution of this phenomenon for even longer chainsand to study the impact of defaults in the chains.
300 400 5000246 a) ⊥ polarization Q Q Q Q
400 600 800024681012 || polarization b) λ (nm) Q Q Q Q Fig. 8.
Total extinction and scattering cross section efficiencies per particle in chains of 10 and 20 particles.In a) the polarization is perpendicular to the symmetry axis and in b) it is parallel to the symmetry axis.
6. Conclusions
The balanced recursive algorithm can give useful and highly accurate information in systemswith large numbers of strongly interacting resonances. This has been demonstrated hereinfor the case of localized plasmon resonances and the studies presented here suggest thatchains of closely spaced localized plasmons can have potentially interesting applicationswith respect to frequency shifting and broadening. Although not demonstrated here, thistechnique also proves useful for treating closely spaced systems possessing surface resonancesof ‘whispering gallery’ type.It is worth remarking that matrix balancing seems to be a useful method to employ inalmost any Foldy-Lax equation solution scheme, be that for direct system matrix inversion,iterative techniques or linear system solutions. In fact, some modern matrix inversion pro-grams actually integrate numerical matrix balancing into their algorithms. Nevertheless,since the matrix balancing in Foldy-Lax equations can be obtained analytically at relatively0low computational cost, it seems beneficial to carry out this balancing explicitly rather thanrelying on purely numerical treatments.The matrix balanced RCTMA has potentially interesting applications for other kinds ofresonance phenomenon, notably whispering Gallery modes. Such studies are currently un-derway. Furthermore, the ability of the matrix balanced RCTMA to study defaults and smallmodifications in large complicated systems is particularly promising and will be employedin subsequent studies.Brian Stout and Alexis Devilez would like to thank Ross McPhedran, Evgeny Popovand Nicolas Bonod for helpful discussions. This work was funded in part by the grantANR-07-PNANO-006-03 ANTARES of the French National Research Agency.
Appendix A: Vector spherical wave functions
The vector spherical wave functions can be readily written in terms of the Vector sphericalharmonics (VSHs) and outgoing spherical
Hankel functions : Ψ ,p ( k r ) ≡ M nm ( k r ) ≡ h + n ( kr ) X nm ( θ, φ ) Ψ ,p ( k r ) ≡ N nm ( k r ) ≡ kr hp n ( n + 1) h + n ( kr ) Y nm ( θ, φ ) + (cid:2) krh + n ( kr ) (cid:3) ′ Z nm ( θ, φ ) i (A1)In the same manner, the regular VSWFs are obtained by replacing the spherical Hankelfunctions in eq.(A1) by spherical Bessel functions. Our adopted definition of the VSHs is Y nm ( θ, φ ) ≡ b r Y nm ( θ, φ ) Z nm ( θ, φ ) ≡ r ∇ Y nm ( θ, φ ) p n ( n + 1) X nm ( θ, φ ) ≡ Z nm ( θ, φ ) ∧ b r (A2)where the Y nm ( θ, φ ) are the scalar spherical harmonics. References
1. B. Stout, J.-C. Auger, and J. Lafait, “A transfer matrix approach to local field calcu-lations in multiple scattering problems,” J. Mod. Opt. , 2129–2152 (2002).2. J.-C. Auger and B. Stout, “A recursive centered T-Matrix algorithm to solve the multi-ple scattering equation : numerical validation,” J. Quant. Spect. & Rad. Trans. ,533–547 (2003).3. A. Doicu and T. Wriedt, Light Scattering by Systems of Particles (Springer, 2006).4. L. Tsang, J. A. Kong, and R. T. Shin,
Theory of Microwave Remote Sensing (JohnWiley & Sons, 1985).15. W. C. Chew,
Waves and Fields in Inhomogeneous Media (IEEE Press, New York,1994).6. M. Lax, “Multiple Scattering of Waves,” Rev. Mod. Phys. , 287–310 (1951).7. B. Stout, C. Andraud, S. Stout, and J. Lafait, “Absorption in multiple scatteringsystems of coated spheres,” J. Opt. Soc. Am. A , 1050–1059 (2003).8. B. Stout, J.-C. Auger, and J. Lafait, “Individual and aggregate scattering matrices andcross sections : conservation laws and reciprocity,” J. Mod. Opt. , 2105–2128 (2001).9. D. W. Mackowski, “Calculation of total cross sections in multiple-sphere clusters,” J.Opt. Soc. Am. A , 2851–2861 (1994).10. O. Moine and B. Stout, “Optical force calculations in arbitrary beams by use of thevector addition theorem,” J. Opt. Soc. Am. B , 1620–1631 (2005).11. M. I. Mishchenko, L. D. Travis, and A. Lacis, Scattering, Absorption and Emission ofLight by Small Particles (Cambridge University Press, 2002).12. C. F. Bohren and D. R. Huffman,
Absorption and Scattering of Light by Small Particles (Wiley-Interscience, New York, 1983).13. E. . D. E. Gray,
American Institute of Physics Handbook 3rd edition (Mcgraw-Hill Tx,1972).14. Z. Li, M. K all, and H. Xu, “Optical forces on interacting plasmonic nanoparticles in afocused Gaussian beam,” Phys. Rev. B , 085,412 (2008).15. S. Bidault, F. J. G. Abajo, and A. Polman, “Plasmon-Based Nanolenses Assembled ona Well-Defined DNA Template,” J. AM. Chem. S.130