Expansion of Arbitrary Electromagnetic Fields in Terms of Vector Spherical Wave Functions
W. L. Moreira, A. A. R. Neves, M. K. Garbos, T. G. Euser, P. St. J. Russell, C. L. Cesar
aa r X i v : . [ phy s i c s . op ti c s ] M a y Expansion of Arbitrary Electromagnetic Fields in Terms of Vector Spherical Wave Functions ∗ W. L. Moreira , A. A. R. Neves , M. K. Garbos , T. G. Euser , P. St. J. Russell , † and C. L. Cesar ‡ Instituto de F´ısica Gleb Wataghin, State University of Campinas, 13.083-970 Campinas, S˜ao Paulo, Brazil National Nanotechnology Laboratory, Istituto di Nanoscienze of CNR,Universit`a del Salento, via Arnesano, 73100, Lecce, Italy and Max Planck Institute for the Science of Light, 91058 Erlangen, Germany (Dated: May 10, 2012)Since 1908, when Mie reported analytical expressions for the fields scattered by a spherical particleupon incidence of an electromagnetic plane-wave, generalizing his analysis to the case of an arbitraryincident wave has proved elusive. This is due to the presence of certain radially-dependent terms in theequation for the beam-shape coefficients of the expansion of the electromagnetic fields in terms of vectorspherical wave functions. Here we show for the first time how these terms can be canceled out, allowinganalytical expressions for the beam shape coefficients to be found for a completely arbitrary incident field.We give several examples of how this new method, which is well suited to numerical calculation, can beused. Analytical expressions are found for Bessel beams and the modes of rectangular and cylindricalmetallic waveguides. The results are highly relevant for speeding up calculation of the radiation forcesacting on spherical particles placed in an arbitrary electromagnetic field, such as in optical tweezers.
PACS numbers: 03.50.De,41.20.-q,42.25.-p,42.25.Bs
Gustav Mie, in his celebrated 1908 paper [1], used thevector spherical wave function (VSWF), or partial waveexpansion (PWE), of a linear polarized plane-wave to gen-eralize scattering theories to spherical particles of any size,from geometrical optics to the Rayleigh regime, and thuswas able to clarify many phenomena, for example in atmo-spheric physics. He obtained analytical expressions for theexpansion coefficients based on special mathematical iden-tities related to a plane-wave. This beam expansion wasnecessary for applying boundary conditions at a sphericalinterface. Since then, with the arrival of lasers and opticalwaveguides, the diversity and complexity of possible inci-dent fields has become enormous so that the restriction toan incident plane-wave has become unrealistic.Different experiments, ranging from particle levitationand trapping [2, 3], to the ultrahigh-Q microcavities usedin cavity QED experiments [4, 5], use different beams. Forexample, very high numerical aperture beams are used inoptical tweezers and confocal microscopy [6–8], evanes-cent fields in near-field microscopy [9, 10], and the waveg-uide modes of a fiber taper are employed to couple lightto the whispering gallery modes of spherical microcavi-ties [11]. Optical forces, absorption, Raman scattering andfluorescence can be greatly enhanced inside spherical mi-crocavities at Mie resonances [12–15]. Laguerre-Gaussian,Hermite-Gaussian and Bessel beams [16, 17], and the in-ternal electromagnetic field of hollow core photonic crys-tal fibers [18, 19], are used to trap and transport particles.The understanding of all these phenomena requires a pre-cise knowledge of the VSWF coefficients of the incidentbeams. A generalized Lorenz-Mie theory was developed tohandle the many variants of beams beyond classical plane-waves, and the expansion coefficients in these cases areknown as beam shape coefficients (BSC) [20, 21]. More- over, because the VSWFs form an orthogonal complete ba-sis, they can be used to study scattering and forces [22]on non-spherical particles, and are the starting point of thepowerful T-matrix methods [23].The calculation of BSCs for an arbitrary beam has al-ways been a complicated task, requiring significant effort.Furthermore, there is a fundamental problem with thesecalculations: an expansion of any function in some basisis complete only when the expansion coefficients can bewritten in terms of scalar products, or integrals, with de-fined numerical values. This task has actually never beenaccomplished for the VSWFs of an arbitrary beam becausethe integral over the solid angle does not explicitly elimi-nate the radial dependence, at least up until now. So far aswe are aware, the current literature lacks any mathematicalproof that this radial function, which appears after integra-tion over all solid angles for any type of beam that satisfiesMaxwell’s equations, can exactly cancel out the sphericalBessel function that appears on the other side of the BSCequation. If this is not true, then the BSC could not be aconstant independent of the radial coordinate – as requiredfor a successful expansion.This non-radial dependence of the BSC has been provenonly for the case of plane-waves and for a high numeri-cal aperture focused Gaussian beam [24]. Working with anelectromagnetic mode inside a hollow cylindrical waveg-uide, we also have been able to obtain analytical expres-sions for constant BSCs that depend only on the positionof the reference frame. This raises the fundamental ques-tion, whether it would be possible to prove that the spheri-cal Bessel function would naturally emerge from the solidangle integral for any type of electromagnetic field. Thepurpose of this letter is to show, we believe for the firsttime, that this is indeed possible. The implications of thisresult for computational light scattering is very noteworthy.We show how the new method can be used to calculate theBSCs for plane-waves, cylindrical and rectangular waveg-uide modes and Bessel beams.The dimensionless BSCs G T E/T Mlm for an incident field E = E ( r ) , H = H ( r ) are defined in the equations [25] (cid:20) E Z H (cid:21) = E X p,q (cid:20) G T Elm G T Mlm (cid:21) M lm + (cid:20) G T Mlm − G T Elm (cid:21) N lm , (1) E is an electric field dimension constant, k N lm = i ∇ × M lm , M lm = j l ( kr ) X lm ( ˆr ) , j l ( kr ) are spherical Besselfunctions, X ( ˆr ) = L Y lm ( ˆr ) / p l ( l + 1) are the sphericalharmonics, Z = p µ/ε , k = ω √ µε and L = − i r × d/d r , d/d r = ˆx ∂ x + ˆy ∂ y + ˆz ∂ z is the gradient operator in thecoordinates (direct) space. Throughout this paper we usethe convention that the terms inside [] are parts of separateequations, i.e., (1) contains two equations, the first relating E to G T Elm and G T Mlm and the second, Z H to G T Mlm and − G T Elm .The usual procedure for obtaining the BSC’s involvesmultiply both sides of (1) by X ∗ l ′ m ′ , take scalar productswith the fields and integrate over the solid angle Ω . Dueto the orthogonality properties of the vector spherical har-monics [25, 26], one can easily show that E j l ( kr ) (cid:20) G T Elm G T Mlm (cid:21) = Z d Ω( ˆr ) X ∗ lm ( ˆr ) · (cid:20) E ( r ) Z H ( r ) (cid:21) , (2) Ω( ˆr ) is the solid angle with respect to an arbitray originnot related to any particular point of the incident beam, ex-plicit said to be in the direct space. Equation (2) does notyield explicit expressions for the BSCs because the LHSstill contains the radially-dependent spherical Bessel func-tion. Our goal is to extract this function from the RHS andcancel it out with the one on the LHS, for any general inci-dent electromagnetic field. To accomplish this we use theFourier transform F of the fields (cid:20) E ( r ) H ( r ) (cid:21) = 1(2 π ) / Z d k ′ (cid:20) E ( k ′ ) H ( k ′ ) (cid:21) e i k ′ · r . (3)By this definition, one can show that F { L ψ ( r ) } = L Ψ( k ) and that the angular momentum operator in re-ciprocal k-space L has the same form as in real r-space (is Hermitian) and is given by L = − i k × d/d k , i.e. d/d k = ˆx ∂ k x + ˆy ∂ k y + ˆz ∂ k z is the gradient operator in the Fourier (reciprocal)space. Using this property and the Rayleigh expan-sion e i k ′ · r = 4 π P ∞ l =0 i l j l ( k ′ r ) P lm = − l Y lm ( ˆr ) Y ∗ lm ( ˆk ′ ) and making M lm ( k ) = j l ( kr ) X lm ( ˆk ) , X lm ( ˆk ) = L Y lm ( ˆk ) / p l ( l + 1) , we obtain j l ( kr ) (cid:20) G T Elm G T Mlm (cid:21) = i l E r π Z d k ′ M ∗ lm · (cid:20) E Z H (cid:21) . (4)Now, only Fourier transforms of the form k F ( k ′ ) = δ ( k ′ − k ) F k ( ˆk ′ ) will represent a field F ( r ) that satisfies the wave equation ∇ F + k F = 0 in three dimensions.So we define the ˆk -only dependent fields (cid:20) E ( k ′ ) H ( k ′ ) (cid:21) = δ ( k ′ − k ) k ′ " E k ( ˆk ′ ) H k ( ˆk ′ ) . (5)Imposing this restriction one finally obtains that (cid:20) G T Elm G T Mlm (cid:21) = i l E r π Z d Ω k ′ X ∗ lm ( ˆk ′ ) · " E k ( ˆk ′ ) Z H k ( ˆk ′ ) . (6)To calculate the coefficients placed at arbitrary position r independent of the coordinate system of the fields onecan use the translation property of Fourier transform [26].Since this form is free of any radially-dependent function,our goal has been achieved.To calculate these Fourier transforms numerically itcould be convenient to use Laplace series [26], alsoknown as spherical harmonic transforms [27]. Us-ing the Laplace series expansion [ E ( ˆk ′ ) , Z H ( ˆk ′ )] = P p ′ ,q ′ [ e p ′ ,q ′ , h p ′ ,q ′ ] Y ∗ p ′ ,q ′ ( ˆk ′ ) one obtains (cid:20) G T Elm G T Mlm (cid:21) = i l E r π X p ′ ,q ′ (cid:20) e p ′ ,q ′ h p ′ ,q ′ (cid:21) · h p ′ , q ′ | L | p, q i p l ( l + 1) . (7)The e l ′ ,m ′ / h l ′ ,m ′ coefficients can be calculated by sev-eral algorithms freely available in internet and the matrix h p ′ , q ′ | L | p, q i is shown in most quantum mechanics books.From now on we shall use the parameter p = ± (1) and will write the components of the fiels in the complexcircular basis ( ˆe − , ˆz , ˆe + ) , so that the components are inthe form a = [ a − , a z , a + ] in which a ± = ˆe ± · a and a z = ˆz · a . We also have that ˆe ± = ( ˆx ± i ˆy ) / √ . Plane wave:
The fields of a arbitrarily polarized planewave are given by [ E ( r ) , Z H ( r )] = E [ ˆ ǫ , ˆk × ˆ ǫ ] e i k · r , ˆ ǫ is the unit polarization vector, therefore the angular Fouriertransform is easily calculated, one obtains from equation(6) that (cid:20) G T Elm G T Mlm (cid:21) = i l π X ∗ lm ( ˆk ) · (cid:20) ˆ ǫ ˆk × ˆ ǫ (cid:21) . (8)For the special case ˆk = ˆz and for a circularly polarizedwave with ˆ ǫ = ˆe p we have (cid:20) G T Elm G T Mlm (cid:21) = i l q π (2 l + 1) δ lp (cid:20) ∓ i (cid:21) , (9) δ n,m is the Kronecker delta function. This is the resultshown in Jackson’s book [25] at less a √ factor due toour choice of complex basis. General hollow waveguide mode:
The TM and TEmodes of a hollow, cylindrical waveguide of arbitrarycross-sectional shape are given as function of g ( r ) = g ρ ( ρ ) e ik z z , g ρ ( ρ ) is the scalar solution of the transversewave equation ( ∂ x + ∂ y + γ ) g ρ ( ρ ) = 0 satisfying theboundary conditions at the waveguide surfaces [25], (cid:20) E T M Z H T E (cid:21) = k z k ˆz × (cid:20) − Z H T M E T E (cid:21) (10) = E (cid:20) ˆz + i k z γ dd ρ (cid:21) g ( r ) , (11) k z is the wavevector in the z direction, k = k z + γ , γ isthe transverse wavevector and d/d ρ = ˆx ∂ x + ˆy ∂ y is thetransverse gradient operator.It should be useful to introduce the cylindrical coordinatesystem in the r-space and k-space. In the r-space we have ρ = x ˆx + y ˆy , ρ = ρ · ρ and ˆ ρ = ρ /ρ . We have also φ = − y ˆx + x ˆy and ˆ φ = φ /ρ . In the same way, in k-space we change ρ → γ and φ → ζ , and for sphericalcoordinates we also change θ → ξ obtaining an equivalentsystem of coordinates. So, the Fourier transforms of fields(10) and (11) are given by (cid:20) E T M Z H T E (cid:21) = E G ( k ′ ) (cid:20) ˆz − k z γ γ ′ (cid:21) (12) (cid:20) E T E − Z H T M (cid:21) = E kγ ′ γ G ( k ′ ) ˆ ζ ′ . (13)Now, similar to the previous case, only G ρ ( k ′ ) = G γ ( ζ ′ ) δ ( γ ′ − γ ) /γ can represent in the Fourier domaina scalar function g ρ ( ρ ) that satisfies the transverse waveequation ( ∂ x + ∂ y + γ ) g ρ ( ρ ) = 0 . So, we have G ( k ′ ) = G γ ( ζ ′ ) δ ( γ ′ − γ ) δ ( k ′ z − k z ) /γ ′ = G γ ( ζ ′ ) δ ( k ′ − k ) δ ( ζ ′ − ζ ) / (sin ξ ′ k ′ ) . Applying these results to (12) and (13) andusing (6) we obtain (cid:20) G T Elm [ T M ] G T Mlm [ T E ] (cid:21) = 2 i l m p l ( l + 1) k γ Q ml (cid:18) k z k (cid:19) G mγ (14) (cid:20) G T Mlm [ T M ] − G T Elm [ T E ] (cid:21) = 2 i l − p l ( l + 1) Q ml ′ (cid:18) k z k (cid:19) G mγ , (15) Q ml ( x ) is the normalized Associated Legendre Polynomialin the way that Y ml ( θ, φ ) = Q ml (cos θ ) e imφ and cos ξ = k z /k , Q ml ′ ( x ) is the derivative of Q ml ( x ) , G mγ = Z dζ ′ G γ ( ζ ′ ) e − imζ ′ . (16) Rectangular hollow metallic waveguide mode:
Inthis case there are two scalar functions g ρ , given by g T Eρ ( ρ ) = cos( k ma x ′ ) cos( k nb y ′ ) and g T Mρ ( ρ ) =sin( k ma x ′ ) sin( k nb y ′ ) , k ma = mπ/a and k nb = nπ/b , ζ mn = arctan( k nb /k ma ) = arctan( na/mb ) and γ mn = π p m /a + n /b , m , n integers and the origin at thelower left corner of the waveguide [25]. If the origin isplaced at ( x , y , then x ′ = x + x and y ′ = y + y . TheFourier transforms of these fields can be easily calculated,and the result is a sum of four Dirac delta functions a thepoints ( k ma , k nb ) , ( − k ma , k nb ) , ( k ma , − k nb ) , ( − k ma , − k nb ) . These delta functions can be written in cylindrical coor-dinates to give us the function G γ ( ζ ′ ) . The integral over ζ ′ given by G qγ in (16) – the plus sign is for the TE-waveguidemode and the minus sign is for the TM-waveguide mode –is G mγ = α (cid:0) e imζ Re { i m e iφ − } ± e − imζ Re { i m e iφ + } (cid:1) (17)in which α = πi − m e ik z z and φ ± = k x x ± k y y . TheBSCs can now be found substituting (17) in (14) and (15). Cylindrical hollow metallic waveguide mode:
Thescalar solution for the electromagnetic fields in termsof cylindrical coordinates with the origin on axisfor the metallic waveguide are given by g ρ ( ρ ) = J m ( γ mn ρ ) e ± imφ , J m ( x ) are order m Bessel functions, J ′ m ( x ) is the derivative of J m ( x ) , γ m,n = χ m,n /R or γ m,n = χ ′ m,n /R , with R being the radius of the cylinder, χ m,n the n -th root of J m ( x ) for the TM mode and χ ′ m,n being the n -th root of J ′ m ( x ) for the TE mode [25]. Fromnow on we define the functions ψ m ( k ; r ) = J m ( γρ ) e ismφ e ik z z (18) Ψ m ( k ; ˆk ′ ) = √ π ( − i ) m e ismζ ′ δ ( ξ ′ − ξ )sin ξ ′ (19)where again s = ± (1) and remember the addition theo-rem [28] to express the scalar function in terms of a newcoordinate system r = r ′ + r , written as a convolution ψ m ( k ; r ′ + r ) = ∞ X j = −∞ ψ m − j ( k ; r ) ψ j ( k ; r ′ ) . (20)The Fourier transform of ψ m ( k ; r ) can be written as Ψ m ( k ; ˆk ′ ) δ ( k ′ − k ) /k ′ . So, we can say that for a cylin-drical waveguide g ρ ( ρ ) = ψ M ( k ; r ) , where M denotesthe propagating mode and this way G γ ( ζ ′ ) = √ π ∞ X j = −∞ ψ M − j ( k ; r )( − i ) j e isjζ ′ (21) G mγ = 2 π ( − i ) sm ψ M − sm ( k ; r ) . (22)The BSC’s can now be found substituting (22) in (14) and(15). The on-axis case can be obtained by setting r = 0 ,which implies that J M − sm ( γ mn ρ ) = δ m,sM , and there-fore G mγ = 2 π ( − i ) M δ m,sM . (23)As we can see the sum over m in (1) will disappear. Bessel Beams:
We have calculated BSCs for two kindBessel beams electromagnetic fields that obey Maxwellequations [29]. They are derived using vector potential andLorentz gauge [30] and are given as function of ψ alreadydefined in (18) and (20) and p = ± (1) . (cid:20) E z Z H z (cid:21) = E (cid:20) ψ M ˆz + ∇∇ · [ ψ M ˆz ] /k − i ∇ × [ ψ M ˆz ] /k (cid:21) (24) (cid:20) E p Z H p (cid:21) = E (cid:20) ψ m − ˆe p + ∇∇ · [ ψ M − ˆe p ] /k − i ∇ × [ ψ M − ˆe p ] /k (cid:21) . (25)The Fourier transform of these fields are given by (cid:20) E z Z H z (cid:21) = E Ψ M " ˆz − ˆk ′ ( ˆk ′ · ˆz ) ˆk ′ × ˆz (26) (cid:20) E p Z H p (cid:21) = E Ψ M − " ˆ e p − ˆk ′ ( ˆk ′ · ˆ e p ) ˆk ′ × ˆ e p . (27)The k ′ component is obviously null by orthogonality withangular momentum. Using again the addition theorem (20)and making c lm ± = p l ( l + 1) − q ( q ± one can showthat the coefficients are given by (28) and (29).In conclusion, we have shown that radially-independentamplitudes (the BSCs) of a complete set of vector spheri- cal wave functions can be calculated explicitly for an arbi-trary electromagnetic field. We have shown how this resultcan be used to determine the BSCs for several beam-typescommonly employed in photonics, although of course themethod is not restricted to applications within the field ofoptics. This new-found ability to evaluate the BSCs ofthe VSWFs analytically makes it much easier to explorerapidly the influence of experimental parameters in prac-tical field scattering and optical tweezer systems. Withthis analytical breakthrough, the long-standing problem ofevaluating the BSCs for an arbitrary field has been solvedand the non-radial dependence of the BSCs proven, allow-ing one to avoid unnecessary approximations in the numer-ical evaluation of these quantities. (cid:20) G T Elm G T Mlm (cid:21) z = 4 πi l − sm p l ( l + 1) ψ M − sm ( k ; r ) (cid:20) mQ ml ( k z /k ) i ( γ/k ) Q ml ′ ( k z /k ) (cid:21) (28) (cid:20) G T Elm G T Mlm + ip ( k z /k ) G T Elm (cid:21) p = 4 πi l − s ( m − p ) p l ( l + 1) ψ M − − s ( m − p ) ( k ; r ) (cid:20) c lm − p Q m − pl ipm ( γ/k ) Q ml (cid:21) (29) ∗ W. L. Moreira and C. L. Cesar belong to CEPOF-FAPESPand INFABIC-CNPq and are grateful to FAPESP, CNPq,CAPES and PETROBRAS for financial assistance. We alsothank Rodrigo G. Pereira at Kavli Institute for TheoreticalPhysics (UCSB) for theoretical discussions about the angu-lar momentum operator in k-space. † [email protected] ‡ lenz@ifi.unicamp.br[1] G. Mie, Ann. Phys. , 377 (1908)[2] A. Ashkin, Phys. Rev. Lett. , 156 (1970)[3] A. Ashkin, J. M. Dziedzic, and T. Yamane, Nature , 769(1987)[4] L. Collot, V. Lefevreseguin, M. Brune, J. M. Raimond, andS. Haroche, Europhys. Lett. , 327 (1993)[5] D. W. Vernooy, A. Furusawa, N. P. Georgiades, V. S.Ilchenko, and H. J. Kimble, Phys. Rev. A , R2293 (1998)[6] K. Svoboda and S. M. Block, Ann. Rev. Biophys. Biomol.Struct. , 247 (1994)[7] S. Hell, G. Reiner, C. Cremer, and E. H. K. Stelzer, J. Mi-crosc. , 391 (1993)[8] G. J. Brakenhoff, P. Blom, and P. Barends, J. Microsc. ,219 (1979)[9] E. Betzig and J. K. Trautman, Science , 189 (1992)[10] E. J. Sanchez, L. Novotny, and X. S. Xie, Phys. Rev. Lett. , 4014 (1999)[11] M. Cai, O. Painter, and K. J. Vahala, Phys. Rev. Lett. , 74(2000)[12] A. A. R. Neves, A. Fontes, L. Y. Pozzo, A. A. de Thomaz,E. Chillce, E. Rodriguez, L. C. Barbosa, and C. L. Cesar,Opt. Express , 13101 (2006) [13] F. Vollmer, D. Braun, A. Libchaber, M. Khoshsima,I. Teraoka, and S. Arnold, Appl. Phys. Lett. , 4057 (2002)[14] S. M. Spillane, T. J. Kippenberg, and K. J. Vahala, Nature , 621 (2002)[15] J. Ng, C. Chan, P. Sheng, and Z. Lin, Opt. Lett. , 1956(2005)[16] J. Arlt and K. Dholakia, Opt. Comm. , 297 (2000)[17] L. Novotny, E. J. Sanchez, and X. S. Xie, Ultramicroscopy , 21 (1998)[18] P. S. J. Russell, Science , 358 (2003)[19] T. G. Euser, M. K. Garbos, J. S. Y. Chen, and P. S. J. Russell,Opt. Lett. , 3674 (2009)[20] G. Gouesbet, B. Maheu, and G. Grehan, JOSA-A , 1427(1988)[21] G. Gouesbet, J. Quant. Spectrosc. Radiat. Transfer ,1223 (2009)[22] T. A. Nieminen, H. Rubinsztein-Dunlop, and N. R. Hecken-berg, J. Quant. Spectrosc. Rad. Transfer , 627 (2001)[23] M. I. Mishchenko, L. D. Travis, and D. W. Mackowski, J.Quant. Spectrosc. Rad. Transfer , 535 (1996)[24] A. A. R. Neves, A. Fontes, L. A. Padilha, E. Rodriguez,C. H. B. Cruz, L. C. Barbosa, and C. L. Cesar, Opt. Lett. ,2477 (2006)[25] J. D. Jackson, Classical Electrodynamics , 3rd ed. (Wiley,1999)[26] G. B. Arfken and H. J. Weber,
Mathematical Methods forPhysicists , international ed. (Elsevier, 2005)[27] R. Suda and M. Takami, Math. Comp. , 703 (2002)[28] M. Abramowitz and I. A. Stegun, Handbook of mathemati-cal functions , 9th ed. (Dover, 1970)[29] K. T. McDonald(1988), arXiv:physics/0006046v1[30] L. W. Davis, Phys. Rev. A19