Simple algorithm for partial wave expansion of plasmonic and evanescent fields

Xinning Yu; Qian Ye; Huajin Chen; Shiyang Liu; Zhifang Lin

doi:10.1364/OE.25.004201

1. Introduction

The T-matrix approach [1–3] rates certainly among the most popular and physically transparent methods in the field of light scattering. Besides the evaluation of the T-matrix itself, another important ingredient in this approach lies in the derivation or (numerical) calculation of the partial wave expansion coefficients (PWECs) that depicts the incident field in terms of the regular vector spherical wave functions (VSWFs) [2–8]. The calculation of the PWECs of a variety of optical beams actually forms the base of the generalized Lorenz-Mie theory [9, 12, 13], where the PWECs are known as beam shape coefficients.

The PWECs for general incident field are generally computed by taking the inner product of the incident field with the VSWFs [9, 12–16], termed projection method, which requires numerically computing integrals of various oscillating functions. As one goes to higher orders of VWSFs for large particle, the highly oscillatory integrands inevitably results in slow convergence, imposing a limitation on efficiency and particle size in solving the light scattering problem. Explicit expressions for the PWECs, together with a stable numerical algorithm, are therefore highly desirable. Among a few cases where the PWECs have been worked out analytically are a homogeneous plane wave impinging in arbitrary direction [2,6–8], Bessel beams of zero [17] and arbitrary order [18], and the evanescent waves [10–13]. The difference between a homogeneous plane wave and an evanescent wave lies actually in the fact that the former has three real Cartesian components of the wave vector, while the latter possesses one complex plus two real Cartesian components of the wave vector. It is this difference that defers the derivation of explicit PWECs by decades, meanwhile an alternative algorithm, the complex-angle approach [19, 20], was proposed, where the conventional Mie theory for homogeneous plane wave can be adopted to tackle the scattering of evanescent wave through a rotation by a complex angle. The explicit PWECs are, however, not yet available for plasmonics field [14, 21], which, differing from the evanescent wave, is characterized by a wave vector with more than one complex Cartesian components.

In this paper, we first work out the explicit PWECs for a complex wave vector (complex-k) field with three Cartesian components k_x, k_y, k_z being arbitrarily constant complex, except for limiting to the case of non-absorptive medium so that the sum of the squares of its three Cartesian components, $k_{x}^{2} + k_{y}^{2} + k_{z}^{2} = ω^{2} ε μ$ , is positive. Here ε and μ are permittivity and permeability of medium where the optical field exists. The complex-k field contains the evanescent waves and simple plasmonic fields as its special cases. It actually serves as the basis functions to model any electromagnetic fields. Next, we present a recurrence algorithm to evaluate the PWECs numerically up to arbitrary order of expansion in VSWFs, thus offering an accurate and powerful tool for the scattering of generic electromagnetic fields that can be modelled by a superposition of the complex-k fields such as the evanescent and plasmonic waves. The approach is especially useful for the large particle where numerical integration in the conventional projection method [14–16] becomes more and more numerically challenging and time consuming.

As applications of our approach, we calculate the optical force and torque acting on a spherical particle immersed in the plasmonic field, up to particle radius R as large as over 30 times of the incident wavelength λ and with the cut-off order n_c of VWSF expansion close to 500. Our approach remains stable while the conventional projection method is challenged due to the numerical integration of highly oscillating integrands. It is found that, as the particle radius R increases, the optical torque on an absorptive sphere keeps parallel to the orientation of spin of optical field for some field polarizations, and it changes its direction for some other polarizations.

2. Explicit partial wave expansion

In the source-less non-absorptive medium, the incident electric field of a time harmonic electromagnetic wave with complex-k is given by

E_{inc} = E e^{i k \cdot r} = E_{0} (E_{x} e_{x} + E_{y} e_{y} + E_{z} e_{z}) e^{i (k_{x} r_{x} + k_{y} r_{y} + k_{z} r_{z})},

where e_x, e_y and e_z are the unit vectors in the Cartesian coordinates system, E₀ > 0 denotes the amplitude, and three complex Cartesian components, E_x, E_y, and E_z, characterize the polarization of the complex-k field. The position vector r and the complex wave vector k read

r = r_{x} e_{x} + r_{y} e_{y} + r_{z} e_{z}, k = k_{x} e_{x} + k_{y} e_{y} + k_{z} e_{z},

where the three complex Cartesian components of k can take arbitrary complex numbers, except for being subject to the dispersion relation in non-absorptive medium

k \cdot k = k_{x}^{2} + k_{y}^{2} + k_{z}^{2} = k^{2} = ω^{2} ε μ > 0,

where k = 2πn_d/λ, and λ denoting the wavelength in vacuum with wave number k₀ = 2π/λ. In our paper the surrounding medium is nonmagnetic, so we denote the medium’s relative permittivity ε_d = ε/ε₀ and relative permeability μ_d = μ/μ₀ = 1 with the permittivity ε₀ and permeability μ₀ in vacuum. Then

n_{d} = \sqrt{ε_{d}}

is the relative index of refraction of the medium. It follows from ∇ · E = 0 that

k_{x} E_{x} + k_{y} E_{y} + k_{z} E_{z} = 0 .

In the T-matrix approach, any monochromatic incident field is expanded in terms of VSWFs. For the complex-k field given by Eq. (1), the expansion reads

E_{inc} = - i E_{0} \sum_{n = 1}^{\infty} \sum_{m = - n}^{n} E_{m n} [p_{m n} N_{m n} (k, r) + q_{m n} M_{m n} (k, r)],

where p_mn and q_mn are the PWECs to be derived, and

E_{m n} = (2 n + 1) i^{n} \frac{(n - m)!}{(n + m)!} .

The regular VSWFs [2,7,8] are given by

\begin{array}{l} L_{m n} (k, r) = [τ_{m n} (\cos θ) e_{θ} + i τ_{m n} (\cos θ) e_{ϕ}] \frac{j_{n} (k r)}{k r} e^{i m ϕ} + e_{r} P_{n}^{m} (\cos θ) \frac{1}{k} \frac{d [j_{n} (k r)]}{d r} e^{i m ϕ} \\ M_{m n} (k, r) = [i τ_{m n} (\cos θ) e_{θ} - τ_{m n} (\cos θ) e_{ϕ}] j_{n} (k r) e^{i m ϕ}, \\ N_{m n} (k, r) = [τ_{m n} (\cos θ) e_{θ} + i τ_{m n} (\cos θ) e_{ϕ}] \frac{1}{k r} \frac{d [r j_{n} (k r)]}{d r} e^{i m ϕ} + e_{r} n (n + 1) P_{n}^{m} (\cos θ) \frac{j_{n} (k r)}{k r} e^{i m ϕ}, \end{array}

where e_r, e_θ and e_ϕ are the unit vectors in the spherical coordinates system

r = \sqrt{r_{x}^{2} + r_{y}^{2} + r_{z}^{2}}

, and j_n(x) is the spherical Bessel functions. Two auxiliary functions are

π_{m n} (z) = \frac{m}{\sqrt{1 - z^{2}}} P_{n}^{m} (z), τ_{m n} (z) = - \sqrt{1 - z^{2}} \frac{d P_{n}^{m} (z)}{d z}

where

P_{n}^{m} (z)

is the associated Legendre function of the first kind of type 2 [22], defined by

P_{n}^{m} (z) = \frac{{(- 1)}^{m}}{Γ (1 - m)} \frac{{(1 + z)}^{m / 2}}{{(1 - z)}^{m / 2}}_{2} F_{1} (- n, n + 1; 1 - m; \frac{1 - z}{2}),

with the integers n and m satisfying n ≥ 1 and −n ≤ m ≤ n. Here ₂F₁(a, b; c; z) denotes the hypergeometric functions [22,23]. It is remarked that the functions defined in Eq. (9) represents an analytic continuation of the conventional associated Legendre functions within the unit circle |z| < 1 in the cut complex plane [22]. They have branch cuts extending from −∞ to −1 then from +∞ to +1 in the complex z plane and the arguments of (1 ± z) falling inside (−π, π], namely, −π < arg (1 ± z) ≤ π [22], so they are single-valued functions. When z is a real number in the range −1 ≤ Re(z) = z ≤ 1, with arg(1 ± z) = 0, Eq. (9) reduces identically to the conventional associated Legendre functions of the first kind

P_{n}^{m} (x)

[4, 5] that follows the Rodrigues formula

P_{n}^{m} (x) = \frac{{(1 - x^{2})}^{m / 2}}{2^{n} n!} \frac{d^{n + m}}{d x^{n + m}} {(x^{2} - 1)}^{n} .

To derive the PWECs p_mn and q_mn, we first rewrite the electric field given by Eq. (1) into a form showing explicit the polarization

E = E_{0} (E_{x} e_{x} + E_{y} e_{y} + E_{z} e_{z}) = E_{0} (p_{1} {\hat{θ}}_{k} + p_{2} {\hat{ϕ}}_{k}),

where we have defined the complex unit vectors to describe polar “direction”

{\hat{θ}}_{k}

, azimuthal “direction”

{\hat{ϕ}}_{k}

, as well as the radial direction

\hat{k}

\begin{array}{l} \hat{k} = \sin α \cos β e_{x} + \sin α \sin β e_{y} + \cos α e_{z}, \\ {\hat{θ}}_{k} = \cos α \cos β e_{x} + \cos α \sin β e_{y} - \sin α e_{z}, \\ {\hat{ϕ}}_{k} = - \sin β e_{x} + \cos β e_{y}, \end{array}

which form an orthogonal set satisfying

\hat{k} \cdot {\hat{θ}}_{k} = {\hat{θ}}_{k} \cdot {\hat{ϕ}}_{k} = \hat{k} \cdot {\hat{ϕ}}_{k} = 0, \hat{k} \cdot \hat{k} = {\hat{θ}}_{k} \cdot {\hat{θ}}_{k} = {\hat{ϕ}}_{k} \cdot {\hat{ϕ}}_{k} = 1 .

The complex direction cosines of $\hat{k}$ , viz, (sin α cos β, sin α sin β, cos α), are governed by

\cos α = \frac{k_{z}}{k}, \sin α = \frac{k_{ρ}}{k}, \cos β = \frac{k_{x}}{k_{ρ}}, \sin β = \frac{k_{y}}{k_{ρ}},

here

k_{ρ} = k {(1 - k_{z}^{2} / k^{2})}^{1 / 2}

. It is remarked that k_x, k_y and k_z take complex values subject to Eq. (3), and the argument of the complex k_ρ lies between

- \frac{π}{2}

and

\frac{π}{2}

, since the square root function here has the branch cuts from −∞ to 1 and from +1 to +∞ in the complex plane. A comparison between Eqs. (1) and (10) yields p₁ and p₂ that depict the complex amplitudes of the p and s polarizations,

p_{1} = \frac{k E_{z}}{k_{ρ}}, p_{2} = \frac{k_{x} E_{y} - k_{y} E_{x}}{k_{ρ}} .

Next we take advantage of the following expansion of a unit dyadic in terms of VSWFs [2],

\hat{I} e^{i k \cdot r} = \sum_{n = 0}^{+ \infty} \sum_{m = - n}^{+ n} [A_{m n}^{(0)} L_{m n} (k, r) + B_{m n}^{(0)} M_{m n} (k, r) + C_{m n}^{(0)} N_{m n} (k, r)],

where Î is the unit dyad, and the vector coefficients read

A_{m n}^{(0)} = i E_{m n} P_{n}^{m} (\cos α) \hat{k} e^{- i m β},

B_{m n}^{(0)} = \frac{- i E_{m n}}{n (n + 1)} [π_{m n} (\cos α) {\hat{θ}}_{k} - i τ_{m n} (\cos α) {\hat{ϕ}}_{k}] e^{- i m β},

C_{m n}^{(0)} = \frac{- i E_{m n}}{n (n + 1)} [τ_{m n} (\cos α) {\hat{θ}}_{k} - i π_{m n} (\cos α) {\hat{ϕ}}_{k}] e^{- i m β} .

The proof of Eq. (15) is given in Appendix A.

Then it follows from Eqs. (1) and (10) that

\begin{array}{l} E_{inc} = E e^{i k \cdot r} = E \cdot \hat{I} e^{i k \cdot r} \\ = E_{0} (p_{1} {\hat{θ}}_{k} + p_{2} {\hat{ϕ}}_{k}) \cdot \sum_{n = 0}^{+ \infty} \sum_{m = - n}^{+ n} [A_{m n}^{(0)} L_{m n} (k, r) + B_{m n}^{(0)} M_{m n} (k, r) + C_{m n}^{(0)} N_{m n} (k, r)], \end{array}

which, together with Eqs. (5), (12), and (16), leads to the explicit expressions for the PWECs p_mn and q_mn,

\begin{array}{l} p_{m n} = \frac{1}{n (n + 1)} [p_{1} τ_{m n} (\cos α) - i p_{2} π_{m n} (\cos α)] {(\cos β - i \sin β)}^{m}, \\ q_{m n} = \frac{1}{n (n + 1)} [p_{1} τ_{m n} (\cos α) - i p_{2} π_{m n} (\cos α)] {(\cos β - i \sin β)}^{m}, \end{array}

As a special case, for the wave propagating in the xz plane (k_y = 0, sin β = 0 and cos β = 1), the p polarization of incident field corresponds to p₁ = 1 and p₂ = 0. In Appendix B, we also present an alternative proof of p polarization case of Eq. (17) by working out analytically the integration in the conventional projection method.

3. Recurrence algorithm for partial wave expansion coefficients

Although the explicit PWECs are given by Eq. (17), their numerical evaluation is nontrivial since the arguments of two auxiliary functions π_mn and τ_mn are complex, see, Eq. (13). In our calculation, we adopt the recurrences due to Xu [7] for computing π_mn and τ_mn with real cos α, see also appendix A of Ref. [11]. It is found that the recurrences keep stable and give correct values for π_mn and τ_mn even when cos α takes complex value z, provided that appropriate branch cuts are drawn.

The recurrence algorithm starts by introducing two renormalized auxiliary functions, and ${\tilde{τ}}_{m, n}$ , defined by

{\tilde{π}}_{m, n} = c_{m n} π_{m n}, {\tilde{τ}}_{m, n} = c_{m n} τ_{m n},

with

c_{m n} = \sqrt{\frac{(2 n + 1) (n - m)!}{n (n + 1) (n + m)!}} .

The renormalized auxiliary functions satisfy [7]

{\tilde{π}}_{- m, n} = - {(- 1)}^{m} {\tilde{π}}_{m, n}, {\tilde{τ}}_{- m, n} = {(- 1)}^{m} {\tilde{τ}}_{m, n},

so only m ≥ 0 cases need to be evaluated. The recurrence relations for their evaluation are

{\tilde{π}}_{n, n} = {[\frac{n (2 n + 1) (1 - z^{2})}{2 (n + 1) (n - 1)}]}^{\frac{1}{2}} {\tilde{π}}_{n - 1, n - 1}, {\tilde{π}}_{n - 1, n} = {[\frac{(n - 1) (2 n + 1)}{(n + 1)}]}^{\frac{1}{2}} z {\tilde{π}}_{n - 1, n - 1},

\begin{array}{l} {\tilde{π}}_{m, n} = - {[\frac{(2 n + 1) (n - 1) (n + m - 1) (n - m - 1) (n - 2)}{(n + 1) (n - m) (n + m) n (2 n - 3)}]}^{\frac{1}{2}} {\tilde{π}}_{m, n - 2} \\ + {[\frac{(2 n + 1) (n - 1) (2 n - 1)}{(n + 1) (n - m) (n + m)}]}^{\frac{1}{2}} z {\tilde{π}}_{m, n - 1}, 0 < m \leq n - 2, \end{array}

{\tilde{τ}}_{0, n} = {[\frac{(2 n + 1) (2 n - 1)}{(n - 1) (n + 1)}]}^{\frac{1}{2}} z {\tilde{τ}}_{0, n - 1} - {[\frac{(2 n + 1) n (n - 2)}{(n - 1) (n + 1) (2 n - 3)}]}^{\frac{1}{2}} {\tilde{τ}}_{0, n - 2}, n > 2,

{\tilde{τ}}_{m, n} = \frac{n z}{m} {\tilde{π}}_{m, n} - \frac{1}{m} {[\frac{(n + m) (n - m) (2 n + 1) (n - 1)}{(2 n - 1) (n + 1)}]}^{\frac{1}{2}} {\tilde{π}}_{m, n - 1}, m > 0 .

The recurrence algorithm proceeds as follows. Suppose one has got the values of ${\tilde{π}}_{m, n}$ and ${\tilde{τ}}_{m, n}$ for n ≤ l and 0 ≤ m ≤ n for each n. Then one can evaluate ${\tilde{π}}_{l + 1, l + 1}$ and ${\tilde{π}}_{l, l + 1}$ based on Eqs. (20a). With ${\tilde{π}}_{m, l + 1}$ for 1 ≤ m ≤ l − 1 obtained by Eq. (20b) and ${\tilde{π}}_{0, n} = 0$ by Eq. (8), one has all ${\tilde{π}}_{m, l + 1}$ for 0 ≤ m ≤ l + 1. Next, the values of ${\tilde{τ}}_{m, l + 1}$ for m = 0 and m > 0 are calculated based on Eqs. (20c) and (20d), respectively. The algorithm starts with the initial values

\begin{array}{l} {\tilde{π}}_{0, n} (z) = 0, {\tilde{π}}_{1, 1} (z) = \frac{\sqrt{3}}{2}, & {\tilde{π}}_{1, 2} (z) = \frac{\sqrt{5}}{2} z, & {\tilde{π}}_{2, 2} (z) = \frac{\sqrt{5}}{2} \sqrt{1 - z^{2}}, \\ {\tilde{τ}}_{0, 1} (z) = - \frac{\sqrt{6}}{2} \sqrt{1 - z^{2}}, & {\tilde{τ}}_{1, 1} (z) = \frac{\sqrt{3}}{2} z, & {\tilde{τ}}_{0, 2} (z) = - \frac{\sqrt{30}}{2} z \sqrt{1 - z^{2}}, \\ {\tilde{τ}}_{1, 2} (z) = \frac{\sqrt{5}}{2} (2 z^{2} - 1), & {\tilde{τ}}_{2, 2} = \frac{\sqrt{5}}{2} z \sqrt{1 - z^{2}} . \end{array}

In all calculation, the branch cuts for the square root function are drawn along the real axis from −∞ to −1 and from +1 to +∞, and $\sqrt{1 - z^{2}} \equiv \sqrt{1 - z} \sqrt{1 + z}$ with −π < arg(1 ± z) ≤ π, so all numerical calculation can be done in any computer language directly.

4. Numerical simulation

To corroborate our approach, we proceed to calculate optical force and optical torque acting on a metallic spherical particle immersed in a plasmonic field [14]. The force and torque are calculated based on Maxwell stress tensor method, with explicit expressions in terms of the PWECs of incident waves and the Mie coefficient of a spherical particle given in Refs. [14–16, 24–29]. The plasmonic field is launched at an interface between a metal plate of permittivity ε_m and a non-dissipative fluid of permittivity ε_d > 0, through total internal reflection using a prism in a Kretschmann configuration [14, 21]. The metallic spherical particle resides in the fluid. The multiple scattering occurring between the particle and the interface is omitted. All parameters are the same as in Ref. [14] for the purpose of comparison, including x axis taken as the propagation direction and z axis perpendicular to the metal plate. The plasmonic field has p polarization characterized by E_x = k_z/k, E_y = 0, and E_z = −k_x/k, with three Cartesian components of the wave vector given by [14]

\frac{k_{x}}{k} = \sqrt{\frac{ε_{m}}{ε_{d} + ε_{m}}}, \frac{k_{y}}{k} = 0, \frac{k_{z}}{k} = - \sqrt{\frac{ε_{d}}{ε_{d} + ε_{m}}},

which satisfy Eq. (3), and p₁ = 1, p₂ = 0 are determined by Eq. (14). So the PWECs given by Eq. (17) are readily evaluated with Eq. (13). We adopt n_d = 1.33 for water as in [14]. The operating wavelength is fixed at λ = 594 nm so ε_m = −8.767+1.535i [30] for gold particle and plate. Our results for non-vanishing components of optical force and torque, normalized by the intensity

I_{0} = ε_{0} E_{0}^{2} n_{d} c / 2

, as a function of particle size are shown in Fig. 1, together with results from Ref. [14]. Perfect agreement is reached, validating our approach. In Fig. 1 we also display optical force and torque on a larger sphere where the direct numerical integration for calculating PWECs in conventional projection method is difficult to converge, manifesting the robustness of our approach for optical field that can be modelled by a superposition of complex-k fields.

Fig. 1 Normalized optical force components f_x, f_z and torque $T_{y}$ acting on a gold sphere in a plasmonic field versus sphere radius R. (a) The force components f_x (left axis, black line) and f_z (right axis, blue line) as a function of R. The red circles and squares denote f_x and f_z from Ref. [14], while the continuous curves represent results based on our approach; (b) The same as panel (a) except for optical torque $T_{y}$ . The open circles and continuous curve exhibit, respectively, results from Ref. [14] and the calculation based on our approach; (c) The force components f_x (left axis, black line) and f_z (right axis, blue line) evaluated based on our approach for larger range of R including the radius regime (R > 2μm) where numerical integration for PWECs is difficult to converge; (d) The same as panel (c) except for optical torque $T_{y}$ .

Download Full Size | PDF

It is remarked that for the expansion of any incident field in terms of VWSFs, see Eq. (5), the summation over index n is practically cut off at n_c, with an empirical formula for n_c given by

n_{c} = x + 4.05 x^{1 / 3} + 2,

for the plane wave [31]. Here x = 2πn_d R/λ is the size parameter [5, 31]. It is noted that n_c proposed in Eq. (23) is actually not large enough for the scattering of arbitrary complex-k fields. In our simulation, we set the cut-off index equal to 1.5 times of n_c given in Eq. (23), to ensure an accurate description of complex-k field and convergence in simulation results. It is this large n_c for large radius that challenges the convergence of numerical integration and results in a difficulty in evaluating PWECs by conventional projection method.

It is observed that the force and torque display oscillations every δR ≈ λ/10, as shown in Figs. 1(a) and 1(b). Such oscillations originate from multipolar resonance effects [14]. It seems that the oscillation behavior disappears for R > 6 μm, as displayed in Figs. 1(c) and 1(d). The reason is that as the radius increases, the oscillation peaks become narrow for non-dissipative particle and characterized by sharp morphology-dependent resonances for a large sphere [32, 33]. The dissipation of gold ruins these peaks, leaving us with a smooth monotonic curve.

In typical (electric) plasmonic field that has p polarization, where the electric (magnetic) field is parallel (perpendicular) to the plane of incidence in Kretschmann configuration (namely, with E_y = 0), the optical torque exerting on a spherical particle remains parallel to the direction of spin [14, 20], irrespective of the particle radius. The direction of spin of optical field is determined by [34]

s = κ \times η / | κ \times η |,

where k = κ+iη is the complex wave vector, with κ and η denoting the real and imaginary parts of k. In the case corresponding to Fig. 1, the spin direction is determined by the wave vector given in Eq. (22), which means s = −e_y. In Figs. 1(b) and 1(d), the optical torque is seen to point towards −e_y, as expected.

As examples of the applications of our approach, we calculate the optical torque acting on a gold sphere residing in some complex-k fields. Typical results are shown in Fig. 2, where the operating wavelength λ = 594 nm, the index of fluid medium n_d = 1.33 and the permittivity of sphere ε_m = −8.767 + 1.535i remain unchanged as in Fig. 1. The dimensionless complex-k is taken to be k_x/k₀ = −0.40 + 0.30i, k_y/k₀ = 0.90 − 0.20i and k_z/k₀ determined by Eq. (3), so that the spin direction s = 0.71e_x + 0.64e_y − 0.29e_z is based on Eq. (24). The only difference between Figs. 2(a) and 2(b) is the polarizations of electric fields: E_x = −0.68 − 0.10i, E_y = 0.45 + 0.22i for Fig. 2(a) and E_x = 0.87 − 0.36i, E_y = 0.19 − 0.53i for Fig. 2(b), the corresponding E_z is calculated by Eq. (4) separately. The chosen values for the incident fields actually represent some plasmonic fields observed in a rotated coordinate system. In Fig. 2(a), the direction of optical torque acting on the gold sphere is seen to keeps parallel to the spin direction of the complex-k field, with the ratio of three Cartesian components to total torque $T_{t} = \sqrt{T_{x}^{2} + T_{y}^{2} + T_{z}^{2}}$ , i.e. $T_{x} / T_{t} = 0.71$ , $T_{y} / T_{t} = 0.64$ and $T_{z} / T_{t} = - 0.29$ , regardless of the particle size. On the other hand, in Fig. 2(b) the direction of torque is subject to an obvious change as the particle size varies. Our preliminary results suggest that for any complex wave vector k, we can cast k_y = 0 by a rotation of coordinate system. If, in this new coordinate system, we have either E_y = 0 or E_x = E_z = 0, namely, the optical field exhibits either p or s polarization, as understood based on Kretschmann configuration, then the optical torque will keep pointing towards the direction of spin determined by Eq. (24). On the other hand, if in the new coordinate system, the field is a superposition of p and s polarization, then the optical torque will change its direction as the particle radius varies. Maybe the reason is the unbalanced light pressure forces acting on different parts of the particle [35, 36]. A more in-depth analysis that reveals the underlying physical origin will be given in a separate paper.

Fig. 2 The direction cosines of optical torque ( $T_{x} / T_{t}$ , $T_{y} / T_{t}$ , $T_{z} / T_{t}$ ), as well as the magnitude of torque, as a function of particle radius R when the a dissipative gold particle is immersed in a complex-k field characterized by the same wave vector components k_x/k₀ = −0.40 + 0.30i, k_y/k₀ = 0.90 − 0.20i, and k_z/k₀ determined by Eq. (3) with (a) E_x = −0.68−0.10i, E_y = 0.45+0.22i, with E_z calculated by Eq. (4); (b) E_x = 0.87 − 0.36i, E_y = 0.19 − 0.53i, with E_z calculated by Eq. (4) as well. The direction of torque may keep fixed to the spin direction Eq. (24) of the optical field, or change as the R increases, dependent on the polarization of the complex-k field.

Download Full Size | PDF

5. Conclusion

We have derived explicit expressions of partial wave coefficients for a general complex-k field that includes evanescent and plasmonic fields as two special cases. A stable numerical recurrence algorithm is proposed for their evaluation up to very high order of expansion where conventional projection method is challenged due to the difficulty of convergence in numerical integration. Our approach is corroborated by comparison with previous conventional, but more cumbersome, method, and proves efficient and powerful for the scattering of optical field that can be modelled by a superposition of the complex-k fields which include evanescent waves and simple plasmonic fields as its two special cases. Finally, as examples of applications, our calculation of optical torque exerting on a dissipative spherical particle in a complex-k field shows that the direction of the optical torque may be fixed in the direction of optical spin independent of the particle size, or it may as well changes with particle size, dependent on the polarization of the optical field.

A. Unit dyadic for the expansion coefficients

In this appendix, we present the proof of the expansion Eq. (15) of a unit dyadic in terms of VSWFs. In the following proof, we limit ourselves to the real wave vector k, but the final formula are valid for complex k, as may be understood by analytical continuation, since both sides of Eq. (15) can be regarded as analytic continuation of the case with real wave vector k.

The regular VSWFs form a complete set for regular solution to the vector Helmholtz equation [4, 37], so the product of the unit dyadic Î and exp(i k · r) is expanded in terms of VSWFs as

\hat{I} e^{i k \cdot r} = \sum_{n = 0}^{+ \infty} \sum_{m = - n}^{+ n} [A_{m n}^{(0)} L_{m n} (k, r) + B_{m n}^{(0)} M_{m n} (k, r) + C_{m n}^{(0)} N_{m n} (k, r)],

where the vector coefficients

A_{m n}^{(0)}

,

B_{m n}^{(0)}

and

C_{m n}^{(0)}

depend on k but independent of r. Based on the orthogonality of VSWFs,

\int_{V_{r}} N_{m n} (k, r) \cdot N_{u v}^{*} (p, r) d^{3} r = \frac{2 i^{n} n (n + 1)}{E_{m n}} \frac{π^{2}}{k^{2}} δ (k - p) δ_{m u} δ_{n v},

\int_{V_{r}} M_{m n} (k, r) \cdot M_{u v}^{*} (p, r) d^{3} r = \frac{2 i^{n} n (n + 1)}{E_{m n}} \frac{π^{2}}{k^{2}} δ (k - p) δ_{m u} δ_{n v},

\int_{V_{r}} L_{m n} (k, r) \cdot L_{u v}^{*} (p, r) d^{3} r = \frac{2 i^{n}}{E_{m n}} \frac{π^{2}}{k^{2}} δ (k - p) δ_{m u} δ_{n v},

\int_{V_{r}} M_{m n} (k, r) \cdot N_{u v}^{*} (p, r) d^{3} r = \int_{V_{r}} M_{m n} (k, r) \cdot L_{u v}^{*} (p, r) d^{3} r = \int_{V_{r}} M_{m n} (k, r) \cdot L_{u v}^{*} (p, r) d^{3} r = 0,

where E_mn is given in Eq. (6), it is ready to derive

\int_{V_{r}} [\int_{0}^{\infty} L_{m n}^{*} (p, r) p^{2} d p] e^{i k \cdot r} d^{3} r = \frac{2 i^{n} π^{2}}{E_{m n}} A_{m n}^{(0)},

\int_{V_{r}} [\int_{0}^{\infty} M_{m n}^{*} (p, r) p^{2} d p] e^{i k \cdot r} d^{3} r = \frac{2 i^{n} π^{2} n (n + 1)}{E_{m n}} B_{m n}^{(0)},

\int_{V_{r}} [\int_{0}^{\infty} N_{m n}^{*} (p, r) p^{2} d p] e^{i k \cdot r} d^{3} r = \frac{2 i^{n} π^{2} n (n + 1)}{E_{m n}} C_{m n}^{(0)} .

To work out the integrations on the left sides of Eqs. (27), we first recapitulate the definitions of the scalar and vector spherical harmonics, see, e.g., on P170 of Ref. [2]

Y_{n}^{m} (θ, ϕ) = P_{n}^{m} (\cos θ) e^{i m ϕ}, Y_{n}^{- m} (θ, ϕ) = {(- 1)}^{m} \frac{(n - m)!}{(n + m)!} Y_{n}^{m *},

and

P_{m n} (θ, ϕ) = Y_{n}^{m} (θ, ϕ) e_{r} = P_{n}^{m} (\cos θ) e^{i m ϕ} e_{r},

B_{m n} (θ, ϕ) = \nabla \times [Y_{n}^{m} (θ, ϕ) r] = [i π_{m n} (\cos θ) e_{θ} - τ_{m n} (\cos θ) e_{ϕ}] e^{i m ϕ},

C_{m n} (θ, ϕ) = r \nabla [Y_{n}^{m} (θ, ϕ)] = [τ_{m n} (\cos θ) e_{θ} - i π_{m n} (\cos θ) e_{ϕ}] e^{i m ϕ} .

The expansion of plane wave on P409 of Ref. [4]

e^{i k \cdot r} = \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} i^{n} (2 n + 1) j_{n} (k r) {(- 1)}^{m} Y_{n}^{m} (θ, ϕ) Y_{n}^{- m} (θ^{'}, ϕ^{'}),

implies that

j_{n} (k r) Y_{n}^{m} (θ, ϕ) = \frac{i^{- n}}{4 π} \oint_{4 π} e^{i k \cdot r} Y_{n}^{m} (θ^{'}, ϕ^{'}) d Ω^{'},

where θ′ and ϕ′ denote the polar and azimuthal angles of the wave vector k. Taking the gradient in r-space on Eq. (31) gives

L_{m n} (k, r) = \frac{1}{k} \nabla [j_{n} (k r) Y_{n}^{m} (θ, ϕ)] = \frac{{(- i)}^{n - 1}}{4 π} \oint_{4 π} e^{i k \cdot r} P_{m n} (θ^{'}, ϕ^{'}) d Ω^{'},

and its complex conjugate

L_{m n}^{*} (k, r) = \frac{i^{n - 1} {(- 1)}^{m}}{4 π} \frac{(n + m)!}{(n - m)!} \oint_{4 π} e^{- i k \cdot r} \hat{k} Y_{n}^{- m} (θ^{'}, ϕ^{'}) d Ω^{'},

where one has used Eq. (28) and

P_{m n} (θ^{'}, ϕ^{'}) = P_{n}^{m} (\cos θ^{'}) e^{i m ϕ^{'}} \hat{k}

, which is a counterpart of Eq. (29a) in k-space, with

\hat{k} = k / k

. Multiplying Eq. (33) by k² and integrating over k, one has

\int_{0}^{\infty} L_{m n}^{*} (k, r) k^{2} d k = \frac{i^{n - 1} {(- 1)}^{m}}{4 π} \frac{(n + m)!}{(n - m)!} \int_{V_{r}} e^{- i \hat{k} \cdot r} \hat{k} Y_{n}^{- m} (θ^{'}, ϕ^{'}) d^{3} k,

which, when visualized as an inverse Fourier transform, implies

\int_{V_{r}} [\int_{0}^{\infty} L_{m n}^{*} (p, r) p^{2} d p] e^{i k \cdot r} d^{3} r = 2 π^{2} i^{n - 1} P_{n}^{m} (\cos θ^{'}) e^{i m ϕ^{'}} \hat{k} .

It follows from (27a) that

A_{m n}^{(0)} = - i E_{m n} P_{n}^{m} (\cos θ^{'}) e^{- i m ϕ^{'}} \hat{k},

which is Eq. (16a) if one remembers that θ′ and ϕ′ denote the polar and azimuthal angles of k.

Taking the cross product of Eq. (32) with kr and using the k-space counterpart of Eq. (29a), one derives

\begin{array}{l} M_{m n} (k, r) = - k r \times L_{m n} (k, r) \\ = \frac{i^{- n}}{4 π} \oint_{4 π} {k \times \nabla^{'} [e^{i k \cdot r} Y_{n}^{m} (θ^{'}, ϕ^{'})] - e^{i k \cdot r} k \times [\nabla^{'} Y_{n}^{m} (θ^{'}, ϕ^{'})]} d Ω^{'}, \end{array}

where ∇ operates in k-space. The integration of the first term within the curly brackets vanishes, as may be proved based on the theory of angular momentum, see, e.g., in Sec. 9.6 of Ref. [38], resulting in

M_{m n} (k, r) = \frac{i^{- n}}{4 π} \oint_{4 π} e^{i k \cdot r} B_{m n} (θ^{'}, ϕ^{'}) d Ω^{'},

where B_mn(θ′, ϕ′) is counterpart of Eq. (29b) defined in k-space. Similarly, taking the complex conjugate of Eq. (38), multiplying by k², and integrating over k, one eventually arrives at a form reminiscent of a Fourier transform

\int_{0}^{\infty} M_{m n}^{*} (k, r) k^{2} d k = \frac{i^{n}}{4 π} \frac{{(- 1)}^{m} (n + m)!}{(n - m)!} \int_{V_{k}} e - i^{k \cdot r} B_{- m n} (θ^{'}, ϕ^{'}) d^{3} k,

which suggests

\int_{V_{r}} [\int_{0}^{\infty} M_{m n}^{*} (p, r) p^{2} d p] e^{i k \cdot r} d^{3} r = 2 π^{2} i^{n} \frac{{(- 1)}^{m} (n + m)!}{(n - m)!} B_{- m n} (θ^{'}, ϕ^{'}) .

Eq. (40), together with Eq. (27b) and the definition of B_mn(θ′, ϕ′), leads to Eq. (16b).

Finally, taking the curl of Eq. (38) and noticing $\hat{k} \times B_{m n} (θ^{'}, ϕ^{'}) = C_{m n} (θ^{'}, ϕ^{'})$ , one gets

N_{m n} (k, r) = \frac{1}{k} \nabla \times M_{m n} (k, r) = \frac{i^{- n + 1}}{4 π} \oint_{4 π} e^{i k \cdot r} C_{m n} (θ^{'}, ϕ^{'}) d Ω^{'} .

In a way similar to that from Eq. (38) to Eq. (40), it is easy to reach

\int_{V_{r}} [\int_{0}^{\infty} N_{m n}^{*} (p, r) p^{2} d p] e^{i k \cdot r} d v = 2 π^{2} i^{n - 1} \frac{{(- 1)}^{m} (n + m)!}{(n - m)!} C_{- m n} (θ^{'}, ϕ^{'}),

which, together with Eq. (27c), yields Eq. (16c) in the main text.

Eqs. (32), (38) and (41) correspond to Eqs. (22), (23) and (24) on P172 of Ref. [2]. We observe that both sides of Eq. (15) are analytical functions with the single-valued functions in Eq. (8) and Eq. (9), as a consequence Eq. (15) and Eq. (16) keep effective in the complex plane for the analytic continuation.

B. Integration method for the expansion coefficients

In the previous work Barton and coworkers [39] give the integrals of expansion coefficients. After using the Bessel function’s integral definition and its variant given by Eq. (11.1.16) on P690 of Ref. [40], Almaas and Brevik [15] evaluate the integrals over ϕ with the modified Bessel function I_v (z) = i⁻^v J_v (iz). Then expressions of expansion coefficients only involve integrations of θ. Meanwhile the associated Legendre function ${\tilde{P}}_{n}^{m} (z)$ given in this appendix is different from $P_{n}^{m} (z)$ in the main text, and their relationship is

{\tilde{P}}_{n}^{m} (z) = {(- 1)}^{m} P_{n}^{m} (z) .

The expressions below are in the case of k_y = 0 and p polarization [14]. Consequently, the complex direction cosines and p₁ = 1, p₂ = 0 are determined by Eqs. (13) and (14),

A_{n, m} = Q_{0} (n, m) [\frac{k_{z}}{k} Q_{1} (n, m) - \frac{k_{x}}{k} Q_{2} (n, m)], B_{n, m} = Q_{0} (n, m) Q_{3} (n, m),

Q_{0} (n, m) = {[\frac{(2 n + 1)}{4 π} \frac{(n - m)!}{(n + m)!}]}^{1 / 2} \frac{{(b / R)}^{2}}{n (n + 1) ψ_{n} (k b)},

\begin{array}{l} Q_{1} (n, m) = 2 π \int_{0}^{π / 2} \sin^{2} θ {\begin{cases} \cos (k_{z} b \cos θ) \\ i \sin (k_{z} b \cos θ) \end{cases}} {\tilde{P}}_{n}^{m} (\cos θ) \\ \times [I_{| m - 1 |} (i k_{x} b \sin θ) + I_{| m + 1 |} (i k_{x} b \sin θ)] d θ, \end{array}

Q_{2} (n, m) = 4 π \int_{0}^{π / 2} \sin θ \cos θ {\begin{cases} i \sin (k_{z} b \cos θ) \\ \cos (k_{z} b \cos θ) \end{cases}} {\tilde{P}}_{n}^{m} (\cos θ) I_{| m |} (i k_{x} b \sin θ) d θ,

Q_{3} (n, m) = \frac{- 4 π m}{k_{x} b} \int_{0}^{π / 2} \sin θ {\begin{cases} \cos (k_{z} b \cos θ) \\ i \sin (k_{z} b \cos θ) \end{cases}} {\tilde{P}}_{n}^{m} (\cos θ) I_{| m |} (i k_{x} b \sin θ) d θ,

where integers n, m satisfy n ≥ 1 and −n ≤ m ≤ n, ψ_n (x) = xj_n(x) is the Riccati-Bessel function, and n + m satisfies

{\begin{matrix} e v e n \\ o d d \end{matrix}}

for all braces in Eqs. (44). The surface integrals in Eqs. (44) are carried out over a sphere of arbitrary radius b. The only difference from the previous researches is the exponential decay factor

e^{i k_{z} h}

, which is fixed to

e^{i k_{z} R}

for simplicity [14, 15], dropped in Eq. (44a). In the main text we put

e^{i k_{z} R}

into E₀ in Eqs. (1) and (5).

In the Helmholtz equation, A_n,m and B_n,m are independent of the surface we select [15], which is already confirmed in the non-evanescent waves [41] and the evanescent waves [12]. Considering Q₀ in Eq. (44b) proportion to b/j_n (kb), as a consequence, the integral solutions of Q₁–Q₃ must be proportional to j_n (kb)/b for the completely cancellation of variable b. Observing Q₁ and Q₂ could be related to Q₃ by recursion formulae of associated Legendre functions, we try to work out Q₃ in Eq. (44e) first. Using Euler’s formula e^ix = cos x + i sin x, and replacing I_v (z) with J_v (iz) for Eq. (44e)Q₃ is separated into two integrations’ combination. And each part is Bessel function times associated Legendre functions and exponential term, such as the integration below

T_{0} = \int_{0}^{π} J_{m} (R \sin θ \sin α) e^{(i R \cos θ \cos α)} {\tilde{P}}_{n}^{m} (\cos θ) \sin θ d θ .

This integration is a special case of Gegenbauer’s finite integral given by Eq. (1) on P379 of Ref. [43], and the solution is obtained by using the relationship between associated Legendre functions ${\tilde{P}}_{n}^{m} (z)$ and Gegenbauer polynomials in the complex plane for m ≥ 0, see Eq. (8.936) on P992 of Ref. [44]. Koumandos makes analytic continuation for −n ≤ m < 0 [45]. Finally we obtain

T_{0} = 2 i^{n - m} j_{n} (R) {\tilde{P}}_{n}^{m} (\cos α),

where cos α is determined in Eq. (13). For the real angle 0 ≤ α ≤ π in Eq. (45), Neves proves Eq. (46) by mathematical induction [42] and Cregg proves Eq. (46) by Gegenbauer polynomial [46]. In fact Eq. (46) is presented historically [47] and used in Heisenberg chains [48].

Based on the Bessel functions of integer order and spherical Bessel functions of integer order, see Eqs. (9.1.10) and (10.1.1) of Ref. [23], we obtain

J_{n} (- z) = {(- 1)}^{n} J_{n} (z), j_{n} (- z) = {(- 1)}^{n} j_{n} (z) .

Notating the integrand in Eq. (44e)

f (θ) = \sin θ {\begin{cases} \cos (k_{z} b \cos θ) \\ i \sin (k_{z} b \cos θ) \end{cases}} {\tilde{P}}_{n}^{m} (\cos θ) I_{| m |} (i k_{x} b \sin θ),

and comparing the integrals of f (θ) and f (π − θ) we get

\int_{0}^{π / 2} f (θ) d θ = \frac{1}{2} \int_{0}^{π} f (θ) d θ,

which satisfies n + m both even and odd. Substituting Euler’s formula, Eqs. (46), (47) and (49) into Eq. (44e), after some algebra we finally get

Q_{3} (n, m) = - i^{n} \frac{4 π m}{\sin α} \frac{j_{n} (k b)}{k b} {\tilde{P}}_{n}^{m} (\cos α) .

Emphasizing the range of real angle $0 \leq θ \leq \frac{π}{2}$ and 0 ≤ cos θ ≤ 1 in Eqs. (44c), (44d) and (44e), we use the recurrence relations given by Eqs. (8.733) and (8.735) of Ref. [44], and get the relationship between Q₁ and Q₃ after some algebra,

\begin{array}{l} Q_{1} (n, m) = \frac{(n + m) (n + m - 1)}{2 (2 n + 1) (m - 1)} k_{x} b Q_{3} (n - 1, m - 1) - \frac{k_{x} b Q_{3} (n - 1, m + 1)}{2 (2 n + 1) (m + 1)} \\ - \frac{(n - m + 1) (n - m + 2)}{2 (2 n + 1) (m - 1)} k_{x} b Q_{3} (n + 1, m - 1) + \frac{k_{x} b Q_{3} (n + 1, m + 1)}{2 (2 n + 1) (m + 1)} . \end{array}

Using Eqs. (8.733) of Ref. [44], we get the relation of Q₂ and Q₃,

Q_{2} (n, m) = - \frac{(n + m)}{(2 n + 1)} \frac{k_{x} b}{m} Q_{3} (n - 1, m) - \frac{(n - m + 1)}{(2 n + 1)} \frac{k_{x} b}{m} Q_{3} (n + 1, m) .

Taking Eq. (50) into Eq. (51) and Eq. (52), we obtain the integral solution of Q₁ and Q₂,

\begin{array}{l} Q_{1} (n, m) = \frac{T_{n}}{2} {j_{n - 1} (k b) [{\tilde{P}}_{n - 1}^{m + 1} (\cos α) - (n + m) (n + m - 1) {\tilde{P}}_{n - 1}^{m + 1} (\cos α)] \\ + j_{n + 1} (k b) [{\tilde{P}}_{n + 1}^{m + 1} (\cos α) - (n - m + 1) (n - m + 2) {\tilde{P}}_{n + 1}^{m - 1} (\cos α)]}, \end{array}

Q_{2} (n, m) = T_{n} [(n + m) j_{n - 1} (k b) {\tilde{P}}_{n - 1}^{m} (\cos α) - (n - m + 1) j_{n + 1} (k b) {\tilde{P}}_{n + 1}^{m} (\cos α)],

with the coefficient T_n = 4πiⁿ⁻¹ _(2n + 1). Eq. (50) has the factor j_n (kb)/b and cancels the r-dependent prefactor for B_n,m. However, taking Eqs. (53) into Eq. (44a) could not directly reduce to the factor j_n (kb)/b for A_n,m. We have to transform Eq. (53a) by using the recurrence relations given by Eqs. (8.733) and (8.735) of Ref. [44],

\begin{array}{l} Q_{1} (n, m) = \frac{T_{n}}{\sin α} \frac{j_{n - 1} (k b)}{(2 n - 1)} [(m^{2} + n^{2} - n) {\tilde{P}}_{n}^{m} (\cos α) - (n + m) (n + m - 1) {\tilde{P}}_{n - 2}^{m} (\cos α)] \\ + \frac{T_{n}}{\sin α} \frac{j_{n + 1} (k b)}{(2 n + 3)} [(n - m + 1) (n - m + 2) {\tilde{P}}_{n + 2}^{m} (\cos α) - (m^{2} + n^{2} + 3 n + 2) {\tilde{P}}_{n}^{m} (\cos α)] . \end{array}

Using Eqs. (8.733) of Ref. [44] and the spherical Bessel functions’ recurrence relation given by Eq. (10.1.19) of Ref. [23], we get another form of Q₁ (n, m),

\begin{array}{l} Q_{1} (n, m) = \frac{T_{n} \sin α}{\cos α} [(n + m) j_{n - 1} (k b) {\tilde{P}}_{n - 1}^{m} (\cos α) - (n - m + 1) j_{n + 1} (k b) {\tilde{P}}_{n + 1}^{m} (\cos α)] \\ + \frac{T_{n}}{\cos α \sin α} \frac{j_{n} (k b)}{k b} [n (n - m + 1) {\tilde{P}}_{n + 1}^{m} (\cos α) - (n + 1) (n + m) {\tilde{P}}_{n - 1}^{m} (\cos α)] . \end{array}

Eqs. (55), (53b) and (50) are the explicit expressions of Eqs. (44c), (44d) and (44e) separately. Taking Eqs. (55), (53b) and (50) into Eq. (44a), we work out solutions of A_n,m and B_n,m

\begin{array}{l} A_{n, m} = {[\frac{1}{4 π (2 n + 1)} \frac{(n - m)!}{(n + m)!}]}^{\frac{1}{2}} \frac{4 π i^{n - 1}}{n (n + 1) k^{2} R^{2} \sin α} \\ \times [n (n - m + 1) {\tilde{P}}_{n + 1}^{m} (\cos α) - (n + 1) (n + m) {\tilde{P}}_{n - 1}^{m} (\cos α)], \end{array}

B_{n, m} = - {[\frac{(2 n + 1)}{4 π} \frac{(n - m)!}{(n + m)}]}^{1 / 2} \frac{4 π i^{n} m}{n (n + 1) k^{2} R^{2} \sin α} {\tilde{P}}_{n}^{m} (\cos α) .

The relations between p_mn, q_mn and A_n,m, B_n,m are found after comparing incident fields,

\begin{array}{l} p_{m n} = {(- 1)}^{m} i^{- n + 1} {[\frac{1}{4 π (2 n + 1)} \frac{(n + m)!}{(n - m)}]}^{1 / 2} k^{2} R^{2} A_{n, m}, \\ q_{m n} = {(- 1)}^{m} i^{- n + 2} {[\frac{1}{4 π (2 n + 1)} \frac{(n + m)!}{(n - m)!}]}^{1 / 2} k^{2} R^{2} B_{n, m} . \end{array}

Taking Eqs. (56) into Eqs. (57), and using Eqs. (8.733) of Ref. [44], we finally obtain

p_{m n p} = \frac{1}{n (n + 1)} \times [{{(- 1)}^{m + 1} \sqrt{1 - z^{2}} \frac{d}{d z} {\tilde{P}}_{n}^{m} (z) |}_{z \to \cos α}] = \frac{1}{n (n + 1)} τ_{m n} (\cos α),

q_{m n p} = \frac{1}{n (n + 1)} \times [{(- 1)}^{m} \frac{m}{\sin α} {\tilde{P}}_{n}^{m} (\cos α)] = \frac{1}{n (n + 1)} π_{m n} (\cos α) .

The last equations of Eq. (58a) and Eq. (58b) are due to auxiliary functions in Eq. (8) and the relation in Eq. (43). p_mnp and q_mnp correspond to the Eq. (17) for p polarization (p₁ = 1 and p₂ = 0).

Funding

National Natural Science Foundation of China (NSFC) (Grant No. 11574055); China 973 Projects (Grant No. 2013CB632701 and 2016YFA0301103); Open Project of State Key Laboratory of Surface Physics in Fudan University (KF2016_3). S. Y. Liu is supported by National Natural Science Foundation of China (Grant Nos. 11274277, 11574275), and Zhejiang Provincial Natural Science Foundation of China (LR16A040001).

References and links

1. P. C. Waterman, “Symmetry, unitarity, and geometry in electromagnetic scattering,” Phys. Rev. D 3(4), 825–839 (1971). [CrossRef]

2. L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing (John Wiley & Sons, 1985).

3. M. I. Mishchenko, N. T. Zakharov, N. G. Khlebtsov, T. Wriedt, and G. Videen, “Comprehensive thematic T-matrix reference database: A 2013–2014 update,” J. Quant. Spectrosc. Radiat. Transfer 146, 349–354 (2014), and references therein. [CrossRef]

4. J. A. Stratton, Electromagnetic Theory (McGraw-Hill, 1941).

5. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, 1983).

6. D. W. Mackowski, “Calculation of total cross sections of multiple-sphere clusters,” J. Opt. Soc. Am. A 11(11), 2851–2861 (1994). [CrossRef]

7. Y. L. Xu, “Electromagnetic scattering by an aggregate of spheres,” Appl. Opt. 34(21), 4573–4588 (1995). [CrossRef] [PubMed]

8. Z. F. Lin and S. T. Chui, “Scattering by optically anisotropic magnetic particle,” Phys. Rev. E 69(5), 056614 (2004). [CrossRef]

9. G. Gouesbet and G. Gréhan, Generalized Lorenz-Mie Theories (Springer, 2011). [CrossRef]

10. T. Wriedt and A. Doicu, “Light scattering from a particle on or near a surface,” Opt. Commun. 152(4–6), 376–384 (1998). [CrossRef]

11. A. Doicu, T. Wriedt, and Y. A. Eremin, Light Scattering by Systems of Particles (Springer, 2006). [CrossRef]

12. G. Gouesbet and J. A. Lock, “On the description of electromagnetic arbitrary shaped beams: The relationship between beam shape coefficients and plane wave spectra,” J. Quant. Spectrosc. Radiat. Transfer 162, 18–30 (2015). [CrossRef]

13. J. A. Lock, “Scattering of the evanescent components in the angular spectrum of a tightly focused electromagnetic beam by a spherical particle,” J. Quant. Spectrosc. Radiat. Transfer 162, 95–102 (2015). [CrossRef]

14. A. C. Durand and C. Genet, “Transverse spinning of a sphere in a plasmonic field,” Phys. Rev. A 89(3), 033841 (2014). [CrossRef]

15. E. Almaas and I. Brevik, “Radiation forces on a micrometer-sized sphere in an evanescent field,” J. Opt. Soc. Am. B 12(12), 2429–2438 (1995). [CrossRef]

16. E. Almaas and I. Brevik, “Possible sorting mechanism for microparticles in an evanescent field,” Phys. Rev. A 87(6), 063826 (2013). [CrossRef]

17. J. M. Taylor and G. D. Love, “Multipole expansion of Bessel and Gaussian beams for Mie scattering calculations,” J. Opt. Soc. Am. A 26(2), 278–282 (2009). [CrossRef]

18. J. Chen, J. Ng, P. Wang, and Z. F. Lin, “Analytical partial wave expansion of vector Bessel beam and its application to optical binding”, Opt. Lett. 35(10), 1674–1676 (2010). [CrossRef] [PubMed]

19. A. Y. Bekshaev, K. Y. Bliokh, and F. Nori, “Mie scattering and optical forces from evanescent fields: A complex-angle approach,” Opt. Express 21(6), 7082–7095 (2013). [CrossRef] [PubMed]

20. K. Y. Bliokh, A. Y. Bekshaev, and F. Nori, “Extraordinary momentum and spin in evanescent waves,” Nat. Commun. 5, 3300 (2014). [CrossRef] [PubMed]

21. H. Raether, Surface Plasmons (Springer, 1986).

22. See, http://functions.wolfram.com/PDF/LegendreP2General.pdf. It is noted that Eq. (9) in this paper differs from that given in http://functions.wolfram.com/07.08.26.0005.01 by a factor of (−1)^m.

23. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover Publications, 1972).

24. J. P. Barton, D. R. Alexander, and S. A. Schaub, “Theoretical determination of net radiation force and torque for a spherical particle illuminated by a focused laser beam,” J. Appl. Phys. 66(10), 4594–4602 (1989). [CrossRef]

25. Ø. Farsund and B. U. Felderhof, “Force, torque, and absorbed energy for a body of arbitrary shape and constitution in an electromagnatic radiation field,” Physica A 227(1–2), 108–130 (1996). [CrossRef]

26. J. Ng, Z. F. Lin, C. T. Chan, and P. Sheng, “Photonic clusters formed by dielectric microspheres: Numerical simulations,” Phys. Rev. B 72(8), 085130 (2005). [CrossRef]

27. M. K. Liu, N. Ji, Z. F. Lin, and S. T. Chui, “Radiation torque on a birefringent sphere caused by an electromagnetic wave,” Phys. Rev. E 72(5), 056610 (2005). [CrossRef]

28. A. Salandrino, S. Fardad, and D. N. Christodoulides, “Generalized Mie theory of optical forces,” J. Opt. Soc. Am. B 29(4), 855–866 (2012). [CrossRef]

29. N. Wang, J. Chen, S. Y. Liu, and Z. F. Lin, “Dynamical and phase-diagram study on stable optical pulling force in Bessel beams,” Phys. Rev. A 87(6), 063812 (2013). [CrossRef]

30. A. Cuche, A. C. Durand, E. Devaux, J. A. Hutchison, C. Genet, and T. W. Ebbesen, “Sorting nanoparticles with intertwined plasmonic and thermo-hydrodynamical forces,” Nano Lett. 13, 4230–4235 (2013). [CrossRef] [PubMed]

31. W. J. Wiscombe, “Improved Mie scattering algorithms,” Appl. Opt. 19(9), 1505–1509 (1980). [CrossRef] [PubMed]

32. J. Ng, C. T. Chan, P. Sheng, and Z. F. Lin, “Strong optical force induced by morphology-dependent resonances,” Opt. Lett. 30(15), 1956–1958 (2005). [CrossRef] [PubMed]

33. E. Xifré-Pérez, F. J. García de Abajo, R. Fenollosa, and F. Meseguer, “Photonic binding in silicon-colloid micro-cavities,” Phys. Rev. Lett. 103(10), 103902 (2009). [CrossRef]

34. T. V. Mechelen and Z. Jacob, “Universal spin-momentum locking of evanescent waves,” Optica 3(2), 118–126 (2016). [CrossRef]

35. S. Chang and S. S. Lee, “Optical torque exerted on a sphere in the evanescent field of a circularly-polarized Gaussian laser beam,” Opt. Commun. 151(4–6), 286–296 (1998). [CrossRef]

36. Y. G. Song, S. Chang, and J. H. Jo, “Optically induced rotation of combined Mie particles within an evanescent field of a Gaussian beam,” Jpn. J. Appl. Phys. 38, L380–L383 (1999). [CrossRef]

37. D. Sarkar and N. J. Halas, “General vector basis function solution of Maxwell’s equations,” Phys. Rev. E 56(1) 1102–1112 (1997). [CrossRef]

38. J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, 1999).

39. J. P. Barton, D. R. Alexander, and S. A. Schaub, “Internal and near-surface electromagnetic fields for a spherical particle irradiated by a focused laser beam,” J. Appl. Phys. 64(4), 1632–1639 (1988). [CrossRef]

40. G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical Methods for Physicists, 6th ed. (Elsevier Academic, 2005).

41. A. A. R. Neves, A. Fontes, L. A. Padilha, E. Rodriguez, C. H. B. Cruz, L. C. Barbosa, and C. L. Cesar, “Exact partial wave expansion of optical beams with respect to an arbitrary origin,” Opt. Lett. 31(16), 2477–2479 (2006). [CrossRef] [PubMed]

42. A. A. R. Neves, L. A. Padilha, A. Fontes, E. Rodriguez, C. H. B. Cruz, L. C. Barbosa, and C. L. Cesar, “Analytical results for a Bessel function times Legendre polynomials class integrals,” J. Phys. A: Math. Gen. 39(18), L293–L296 (2006). [CrossRef]

43. G. N. Watson, A Treatise on the Theory of Bessel Functions, 2nd ed. (Cambridge University, 1944).

44. I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products, 7th ed. (Elsevier Academic, 2007).

45. S. Koumandos, “On a class of integrals involving a Bessel function times Gegenbauer polynomials,” Internat. J. Math. Math. Sci. 2007, 1–5 (2007). [CrossRef]

46. P. J. Cregg and P. Svedlindh, “Comment on ‘Analytical results for a Bessel function times Legendre polynomials class integrals’,” J. Phys. A: Math. Theor. 40, 14029 (2007). [CrossRef]

47. B. Podolsky and L. Pauling, “The momentum distribution in hydrogen-like atoms,” Phys. Rev. 34(1), 109–116 (1929). [CrossRef]

48. G. S. Joyce, “Classical Heisenberg model,” Phys. Rev. 155(2), 478–491 (1967). [CrossRef]

Simple algorithm for partial wave expansion of plasmonic and evanescent fields

Abstract

1. Introduction

2. Explicit partial wave expansion

3. Recurrence algorithm for partial wave expansion coefficients

4. Numerical simulation

5. Conclusion

A. Unit dyadic for the expansion coefficients

B. Integration method for the expansion coefficients

Funding

References and links

Cited By

Figures (2)

Equations (80)

Optics Express