1. Introduction

Gustav Mie, in his celebrated 1908 paper [1], used vector spherical wave functions (VSWF), or partial wave expansion (PWE), of a linear polarized plane-wave to generalize scattering theories applied to spherical particles of any size, ranging from geometrical optics to the Rayleigh regime. Mie’s contribution to light scattering had broad applications, clarified many phenomena, such as, for example, the composition of clouds in Venus [2]. To different experiments, ranging from particle levitation and trapping [3, 4], to QED ultrahigh-Q microcavity experiments [5, 6], different beams are used. For example, very high numerical aperture beams are used in optical tweezers and confocal microscopy [7, 8, 9], evanescent fields in near-field microscopy [10, 11], and the waveguide modes of a fiber taper are employed to couple light to the whispering gallery modes of spherical microcavities [12, 13]. Optical forces, absorption, Raman scattering, fluorescence and nonlinear optics can be greatly enhanced inside spherical microcavities at Mie resonances [14, 15, 16, 17, 18]. Laguerre-Gaussian, Hermite-Gaussian and Bessel beams [19, 20], and the internal electromagnetic field of hollow-core photonic crystal fibers [21, 22] can be used to trap and transport particles. To appreciate all these phenomena, a precise knowledge of the VSWF coefficients of the incident beams is required.

The obtained analytical expressions for the expansion coefficients were based on special mathematical identities related to a plane wave. This beam expansion was necessary to apply the boundary conditions to a spherical interface. With the arrival of lasers and optical waveguides, the complexity and diversity of possible incident fields has become enormous so that the restriction to a description based on an incident plane-wave approximation became unrealistic. A generalized Lorenz-Mie theory was developed to handle the many variants of beams beyond the approximation of plane-waves, where the expansion coefficients in these cases are known as beam shape coefficients (BSC) [23, 24, 25, 26]. Literature shows different definitions of beam shape coefficients, therefore we are using the uppercase letter “G” for the BSC to avoid confusion with lowercase letter “g” first used by Gouesbet’s group, that extracted the plane wave prefactor. Moreover, because the VSWFs form an orthogonal complete basis, they can be used to study scattering and optical forces [27] on non-spherical particles, and are the starting point of the powerful T-matrix method [28]. The calculation of BSCs for an arbitrary beam, however, has always been a complicated task, requiring significant theoretical and computational effort. Furthermore, there is a fundamental problem with these calculations: an expansion of any function in some basis is complete only when the expansion coefficients can be written in terms of scalar products, or integrals, with defined numerical values. This task has been previously accomplished only for highly focused beams used for optical trapping, [29, 30], but not for an arbitrary Maxwellian beam, because the integral over the solid angle did not explicitly eliminate the spherical Bessel function radial dependence of RHS. The successful expansion for the VSWFs requires that the BSC is a constant independent of the radial coordinate.

This non-radial dependence of the BSC is trivial for plane-waves, and has been proven for the case of high numerical aperture focused Gaussian beam [29]. Working with an electromagnetic mode inside a hollow cylindrical waveguide, we have also been able to obtain analytical expressions for constant BSCs that depend only on the position of the reference frame. This raised the fundamental question if it would be possible to prove that the spherical Bessel function would naturally emerge from the solid angle integral (see Eq. (2)) for any type of electromagnetic field? The purpose of this work is to first demonstrate that this is indeed possible. The implication of this result for computational light scattering is very noteworthy. Secondly, we apply this new method to calculate the BSCs for plane- waves, cylindrical and rectangular waveguide modes and Bessel beams.

2. Beam shape coefficients determination

The dimensionless BSCs $G_{lm}^{TE/TM}$ for incident harmonic electromagnetic fields E = E(r), H = H(r) are defined in the equations of ref. [31]

The usual procedure to obtain the BSC’s involves taking the scalar product of ${\mathbf{X}}_{l^{\prime }{m}^{\prime }}^{\ast }(\widehat{\mathbf{r}})$, with the fields in Eq. (1) and integrating over the solid angle $\mathrm{\Omega }(\widehat{\mathbf{r}})$. Due to the orthogonality properties of the vector spherical harmonics [31, 32], one can easily show that

Our goal is to extract this function from the RHS and cancel it out with the one on the LHS, for any general incident electromagnetic field that satisfies Maxwell’s equations. The most elegant way to do this is to take the Fourier transform of Eq. (1), including the functions M_lm(r) and N_lm(r), which will result in vector spherical harmonics in the Fourier domain. We then use the orthogonality properties of vector spherical harmonics again, as done in the direct space, to obtain an expression for BSC’s which is independent of the spherical Bessel function, for any kind of field that satisfies the Maxwell equations and the wave equation. The two key steps are:

To accomplish this, we use the Fourier transform ℱ of the fields, defined as,

One can show that $ℱ\{\mathbf{L}\psi (\mathbf{r})\}=ℱ\left\{\mathbf{L}\right\}ℱ\{\psi (\mathbf{r})\}=ℒ\mathrm{\Psi }(\mathbf{k})$ and that the angular momentum operator in reciprocal k-space $ℒ=ℱ\left\{\mathbf{L}\right\}$ has the same form as in real r-space (is Hermitian) and is given by $ℒ=-i\mathbf{k}\times d/d\mathbf{k}$, i.e. $d/d\mathbf{k}=\widehat{\mathbf{x}}{\partial }_{k_x}+\widehat{\mathbf{y}}{\partial }_{k_y}+\widehat{\mathbf{z}}{\partial }_{k_z}$ is the gradient operator in the Fourier (reciprocal) space. To do this we apply the angular momentum operator over a generic function ψ(r) written by means of the Fourier spectrum Ψ(k), as in Eq. (3). So the operator will act over the kernel exp{ik·r}, and obtain

The next step is to write down the vector product by its components and integrate by parts,

As an example, we can take l = x, m = y e n = z to obtain the x component of the integral. As we know that ψ has a Fourier transform, it implies that it vanishes on the extremes of integration. For the x component we have

Finally, we obtain

Using this property, the Rayleigh expansion $e^{i\mathbf{k}\cdot \mathbf{r}}=4\pi {\displaystyle {\sum }_{l=0}^{\infty }{i}^l{j}_l(k^{\prime }r){\displaystyle {\sum }_{m=-l}^l{Y}_{lm}(\widehat{\mathbf{r}})Y_{lm}^{\ast }}({\widehat{\mathbf{k}}}^{\prime })}$ and making ${ℳ}_{lm}(\mathbf{k})=j_l(kr){\mathcal{X}}_{lm}(\widehat{\mathbf{k}})$, in which ${\mathcal{X}}_{lm}(\widehat{\mathbf{k}})=ℒY_{lm}(\widehat{\mathbf{k}})/\sqrt{l(l+1)}$, we obtain

Noting that only Fourier transforms of the form $k^2ℱ({\mathbf{k}}^{\prime })=\delta (k^{\prime }-k){ℱ}_{\mathrm{k}}({\widehat{\mathbf{k}}}^{\prime })$ will represent a field F(r) that satisfies the wave equation ∇²F + k²F = 0 in three dimensions, we define the $\widehat{\mathbf{k}}$-only dependent fields

The Dirac Delta function now is used to extract the spherical Bessel function from the integral in k′, leaving only a solid angle integral, and Eq. (8), becomes

Note that in Eq. (10), with these BSCs, our goal has been achieved since they are now free from any radially-dependent function. To calculate the coefficients placed at arbitrary position r₀ independent of the coordinate system of the fields one can use the translation property of Fourier transform [32]. The translation properties of the BSC’s in the conventional way are quite complicated [33], but it is just a phase term in the Fourier domain. In the examples we used, the translation was done in every case without any complication.

It may be useful to mention that to calculate these Fourier transforms numerically it could be convenient to use Laplace series [32, 34], also known as spherical harmonic transforms [35]. Using the Laplace series expansion $[{\epsilon }_{\mathrm{k}}({\widehat{\mathbf{k}}}^{\prime }),Z{\mathscr{H}}_{\mathrm{k}}({\widehat{\mathbf{k}}}^{\prime })]={\displaystyle {\sum }_{p^{\prime },q^{\prime }}[{\mathbf{E}}_{p^{\prime },q^{\prime }},{\mathbf{H}}_{p^{\prime },q^{\prime }}]Y_{p^{\prime },q^{\prime }}}({\widehat{\mathbf{k}}}^{\prime })$ one obtains

The E_l′,_m′ /H_l′,_m′ coefficients can be calculated by several algorithms freely available in internet and the matrix $〈l,m\left|ℒ\right|p^{\prime },q^{\prime }〉$ is shown in most quantum mechanics books.

3. Gauge invariance and determination of BSCs from the vector potential

We have shown the calculation of BSC’s using the vector electromagnetic fields. But the fields can be written in terms of the scalar Φ and vector A potentials by the equations:

Here, for convenience, we adjusted the dimensions of A₀ to E₀, in a way that A becomes dimensionless, and ZH₀ = E₀ = kcA₀. In the Fourier domain the gradient −∇Φ becomes $-i\mathbf{k}\tilde{\mathrm{\Phi }}$ which tells us that Φ is not taken into account on the calculations of BSCs. This is a very interesting result since the angular momentum does not have a component in the radial direction $\widehat{\mathbf{k}}$. Therefore, there is no dependence on any gauge transformation and only the vector potential is necessary to obtain the BSCs.

The scalar and vector potential also obeys wave equation, so in Fourier space we have $\mathcal{A}={\mathcal{A}}_{\mathrm{k}}\delta (k^{\prime }-k)/k^2$. We can immediately write the electromagnetic fields as

Using these relations in Eq. (10) we obtain

In this way, we can obtain the BSCs directly, without having to explicitly calculate the electromagnetic beams. These results also tell us that one can easily extend the idea to calculate the BSCs using the polarization potentials (Hertz vectors), or scalar potentials like Whittaker, Debye and specially Bromwich, which were used in the original formulation of the GLMT [36, 37, 38]. Another interesting approach is the method used by Davis [39]. Davis looked for solutions of electromagnetic wave equation whose vector potential A has only one component, in some system of coordinates (taking care about the components of the vector laplacian [40]). An example of the use of Davis method is shown in the calculation of the components of the Bessel Beams [41, 42]. Both solutions shown in this work were obtained using the vector potential. It is also interesting to note that the explicit expressions for the Bessel Beams depend on the scalar potential in the Lorenz gauge. Nevertheless, the beam shape coefficients do not care about scalar potential or any gauge (i.e. gauge invariance). This result also opens up the possibility to generate new kinds of beams by directly proposing any new valid vector potential.

4. Derivation of typical BSCs

In this section we shall make use of our main result to obtain BSC for known wavefields. The first test is to obtain the known results for a plane wave. Afterwards, we calculate BSCs for unknown results and test them numerically. We have an expression for generic waveguides, and use them for two specific cases, the rectangular and the cylindrical waveguides.

From now on we shall use the parameter p = ±(1) in order to avoid misunderstandings with other ±’s of different meanings. We write the components of the fields in the complex circular basis $(\widehat{\mathbf{e}}-,\widehat{\mathbf{z}},{\widehat{\mathbf{e}}}_{+})$, in which ${\widehat{\mathbf{e}}}_p=(\widehat{\mathbf{x}}+ip\widehat{\mathbf{y}})/\sqrt{2}$.

4.1. Plane wave

The fields of an arbitrarily polarized plane wave are given by

For the special case $\widehat{\mathbf{k}}=\widehat{\mathbf{\epsilon }}$ and for a circularly polarized wave with $\widehat{\mathbf{\epsilon }}={\widehat{\mathbf{e}}}_p$ we have, using the properties of the scalar spherical harmonics

4.2. General hollow waveguide mode

The TM and TE modes of a hollow waveguide of arbitrary cross-sectional shape are given as function of $g(\mathbf{r})=g_{\gamma }(\mathit{\rho })e^{i{k}_{z}z}$, g_γ(ρ) is the scalar solution of the transverse wave equation

Similar to the previous case, only G_γ(γ′) = G_γ(ζ′)δ(γ′ − γ)/γ can represent in the Fourier domain a scalar function g_γ(ρ) that satisfies the transverse wave equation $({\partial }_x^2+{\partial }_y^2+{\gamma }^2)g_{\gamma }(\mathit{\rho })=0$. Therefore, we have $G({\mathbf{k}}^{\prime })=G_{\gamma }({\zeta }^{\prime })\delta ({\gamma }^{\prime }-\gamma )\delta ({k^{\prime }}_2-k_z)/{\gamma }^{\prime }=G_{\gamma }({\zeta }^{\prime })\delta (k^{\prime }-k)({\xi }^{\prime }-\xi )/(\sin {\xi }^{\prime }{k^{\prime }}^2)$. Applying these results to Eq. (22) and Eq. (23) and using Eq. (10) we obtain

$Q_l^m(x)$ is the normalized Associated Legendre Function in the way that $Y_l^m(\theta ,\varphi )=Q_l^m(\cos \theta )e^{im\varphi }$ and cosξ = kz/k, $Q_l^m{}^{\prime }(x)$ is the derivative of $Q_l^m(x)$. We made this choice instead of the classical $P_l^m$ given in [32] because the $P_l^m$ make the equations unnecessarily big and confusing. Moreover, the Associated Legendre Functions are numerically determined by recurrence, making $Q_l^m$ more stable than $P_l^m$. The Fourier coefficient of the azimuthal component is given by

4.2.1. Rectangular hollow metallic waveguide mode

For this case, the rectangular waveguide (RWG), there are two scalar functions g_γ, given by $g_{\gamma }^{TE}(\mathit{\rho })=\cos (k_a^m{x}^{\prime })\cos (k_b^n{y}^{\prime })$ and $g_{\gamma }^{TM}(\mathit{\rho })=\sin (k_a^m{x}^{\prime })\sin (k_b^n{y}^{\prime })$, $k_a^m=m\pi /a$ and $k_b^n=n\pi /b$, ${\zeta }^{mn}=\arctan (k_b^n/k_a^m)=\arctan (na/mb)$ and ${\gamma }_{mn}=\pi \sqrt{m^2/a^2+n^2/b^2}$, for m and n integers, and the origin is at the lower left corner of the waveguide [31]. For an origin placed at (x₀, y₀, 0) we find x′ = x + x₀ and y′ = y + y₀. The Fourier transforms of these fields can be easily calculated, and the result is a sum of four Dirac delta functions at points $(k_a^m,k_b^n)$, $(-k_a^m,k_b^n)$, $(k_a^m,-k_b^n)$, $(-k_a^m,-k_b^n)$. These delta functions can be written in cylindrical coordinates to give us the function G_γ(ζ′). The integral over ζ′ given by $G_{\gamma }^q$ in Eq. (26) (the plus sign is for the TE-waveguide mode and the minus sign is for the TM-waveguide mode) is

4.2.2. Cylindrical hollow metallic waveguide mode

The scalar solution for the electromagnetic fields in terms of cylindrical coordinates with the origin on axis for the metallic waveguide (shortly CWG) are given by g_γ(ρ) = J_m (γ_mnρ)e^±imϕ, where J_m(x) are Bessel functions of order m, ${J^{\prime }}_m(x)$ is the derivative of J_m(x), γ_m,n = χ_m,n/R or ${\gamma }_{m,n}={{\chi }^{\prime }}_{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 ${{\chi }^{\prime }}_{m,n}$ being the n-th root of ${J^{\prime }}_m(x)$ for the TE mode [31]. From now on we define the function

This function is the solution of the scalar wave equation in cylindrical coordinates and will be largely used from now on. The Fourier transform of ψ_m(k; r) can be written as ${\mathrm{\Psi }}_m(\mathbf{k};{\widehat{\mathbf{k}}}^{\prime })\delta (k^{\prime }-k)/{k^{\prime }}^2$, where

We use s = ±(1) to specify the polarity of the cylindrical beam and distinguish from the sign p of the complex circular base.

To account for a translation in the reference system we express the scalar function in terms of a new coordinate system r = r′ + r₀, remembering the addition theorem [43] written as a convolution

For a cylindrical waveguide we have g_γ(ρ) = ψ_M(k; ρ), where M denotes the propagating mode. The azimuthal Fourier component becomes

The BSC’s can now be found substituting Eq. (32) in Eq. (24) and Eq. (25). The on-axis case can be obtained by setting ρ₀ = 0, which implies that J_M−sm(γ_mnρ₀) = δ_m,sM, and therefore

As we can see the sum over m in Eq. (1), will contribute only for terms m = sM.

4.3. Bessel beams

We have calculated BSCs for two kind of Bessel Beams (BB) electromagnetic fields that obey Maxwell equations [41]. They are derived using vector potential and Lorenz gauge [39] and are given as function of ψ_M already defined in Eq. (28) and Eq. (30) and p = ±(1). The two cases are distinguished by the subscript in the square brackets, where z says that the adimensional vector potential A has only the z component $(\mathbf{A}={\psi }_M\widehat{\mathbf{z}})$, shortly referred as BBZ, and the subscript p says that the vector potential A is circularly polarized $(\mathbf{A}={\psi }_M{\widehat{\mathbf{e}}}_p)$, shortly referred as BBP. The complete electromagnetic fields are written as

Using ψ_M = ψ_M(k;r′+r₀), and the the addition theorem Eq. (30) and making $c_p^{lm}=\sqrt{l\left(l+1\right)-m\left(m+p\right)}$ and W_l,s,m = 4πi^l−sml^−1/2(l + 1)^−1/2 one can show that the coefficients that describe the fields in Eq. (34) and Eq. (35) are given by

It’s easy to verify that if we make $\widehat{\mathbf{k}}=\widehat{\mathbf{z}}$ and M = 0 on the definitions of ψ_M on the Eq. (35) we shall get a circularly polarized plane wave, given by doing $\widehat{\epsilon }={\widehat{\mathbf{e}}}_p$ in the Eq. (17), and the BSC’s in Eq. (37) will become the same given by Eq. (18) when we do $\widehat{\epsilon }={\widehat{\mathbf{e}}}_p$.

5. Validation

The partial wave expansion in Eq. (1) is performed using the vector spherical wave functions M_lm and N_lm, which depends on spherical Bessel functions j_l(kr), that become null for kr ≪ l and show a fast oscillation decay for kr ≫ l. These spherical Bessel functions give to the expansion an ”areal” convergence that follows the rule of thumb l_max ≈ kr_max. Suppose we want to evaluate the expansion in a position r away from the origin, we need to sum from 1 to l^*, in the way that the functions j_l>l^* will give a small contribution to the point. In our examples we expand over a maximum radius smaller than 10 wavelengths, which implies that kr ≈ 60. We used the value 64 for all cases, which allow us to conclude that the truncation with l_MAX = 64 can handle the cases shown with very high precision. These results can be seen in Fig. 1 where we show the increase of the “areal” convergence of the expansion with l_MAX for two distinct directions. The black curve is the plot of the analytical expression. In Fig. 2 we have a reference geometry of the beams calculated, the directions, the profiles of z-component of the Poynting vectors and vectors representing the other components. In Fig. 1 we show the z component of the Poynting vector along the 45 degrees direction (x = y), passing both by origin of the fields on its own system of coordinates and by the origin of the expansion. The vertical cyan lines in the Fig. 1(b) are the limits of the cylindrical waveguide. On the upper axis of each figure we have the value of kr, and we can see that the truncation on l_MAX gives an good approximation for all the points that kr < l_MAX.

 

Fig. 1 z-Component of the Poynting vectors with respect to truncation in l_MAX along the diagonal (y = x).

Download Full Size | PDF

 

Fig. 2 Plots of the components of the Poynting vector of the rectangular and cylindrical TM waveguides and the two cases Bessel Beams: in density is the real part of the z component (normalized); the vector field in blue is the imaginary part of transverse component and in magenta its the real part. If we treat the TE case, the vectors represent the imaginary part (magenta) will invert its directions. All the other values remain unchanged, the same happens if we change p→−p. The green lines indicate the origin of the expansion.

Download Full Size | PDF

We truncate the infinite summation in l_MAX in Eq. (1), which gives 2l +1 terms, remembering that for every l we have a summation in −l ≤ m ≤ l. Therefore, the total number of terms in the summation is l_MAX (l_MAX + 2). The total number of points is N = n_xn_yn_zl_MAX(l_MAX + 2). In order to evaluate the best l_MAX truncation value, and the calculation speed, we wrote a code in C language using parallel computation by OpenMP, and ported this code into R environment [44, 45], which allowed us to run the calculations in a few minutes. The main parameters are n_x = n_y = 200, n_z = 1 and l_MAX = 2,4,8,16,32,64. As basic unit we take a wavelength of λ = 500 nm, which gives a wavenumber of k = 2π/λ. We use k_x = k_y = mπ/a, with m = 4, $\gamma =\sqrt{2}{k}_x$. We used this value of γ for all the cases. For the cylindrical waveguides we also used M = N = 2, and the radius for the TM case obtained by R = χ_MN/γ. We also used p = s = +1, and performed an expansion over a six wavelengths sized square area (a = b = 6λ). The origin of the expansions is r₀ = (1,1,0)a/4 for the rectangular waveguide whose own origin is on bottom right corner and r₀ = −(1,1,0)a/4 for the cylindrical cases, whose own origin is on the center, as shown Fig 2. It is convenient to use Poynting vectors because the electromagnetic TE/TM fields can be obtained from each other by the rule ZH^TE → E^TM and E^TE → −ZH^TM, which implies that S^TE → S^TM*. So, the electromagnetic flux, given by the real part of the z component of S, do not change by the interchanging between TE and TM.

Although BSC’s do not depend on the scalar potential, the gradient term in Eq. (34) and Eq. (35) plays an important role in the components of the fields. We just have to think about the z-component of the Poynting vector: $2S_z=E_x{H}_y^{*}-E_y{H}_x^{8*}$. The transversal components of the electric field are given by the term corresponding to the scalar potential ∇∇·A, that, if dropped out, the z component of the Poynting vector vanishes. Nonetheless this term is not used for the calculation of the BSC’s. To conclude, we cite again J. D. Jackson[46]: “whatever propagation or non propagation characteristics are exhibited by the potentials in a particular gauge, the electric and magnetic fields are always the same and display the experimentally verified properties of causality and propagation at the speed of light”. This enforces the robustness of the vector spherical wave function expansion and our method to calculate the beam shape coefficients.

6. Conclusion

We developed a novel approach, applicable to any vectorial electromagnetic field that satisfies the vector wave equation and Maxwells equations, to obtain radially-independent amplitudes (the BSCs) of a complete set of vector spherical wave functions, without any approximation. We have applied this result to determine the BSCs for several beam-types commonly employed in photonics. Note that the presented approach is not restricted to fields with radial or azimuthal polarization, but that it can quite easily handle the common linear polarized fields not examined here. The new result in this paper is the definition of the BSC, through the use of Eq. (10) or Eq. (11). We also show the gauge invariance of our procedure and call attention that our result could be used to engineer new kinds of beams with some desirable BSCs or to propose a beam by providing its potential vector. We emphasize that there are no approximations in the calculation of the BSCs obtained.

This analytical result has been suggested as a path to solve some open problems in generalized Lorenz–Mie theories and related topics [47]. We note that this method is not restricted to applications within the field of optics. This new-found ability to analytically evaluate the BSCs of the VSWFs makes it much more effortless to rapidly explore the influence of applications. With this contribution, the problem of evaluating the BSCs for an arbitrary field has been solved and the non-radial dependence of the BSCs proven, allowing one to avoid unnecessary approximations in the numerical evaluation of these quantities that results in some residual dependence on the radial coordinate.