Based on an expansion formula for unit dyadic in terms of the vector spherical wave functions, we derive explicit partial wave coefficients for a complex wave vector field that is characterized by a single wave vector with three Cartesian components being arbitrarily constant complex except subject to lossless background constraint and thus includes evanescent waves and simple plasmonic fields as its two special cases. A recurrence method is then proposed to evaluate the partial wave expansion coefficients numerically up to arbitrary order of expansion, offering an efficient tool for the scattering of generic electromagnetic fields that can be modelled by a superposition of the complex wave vector fields such as the evanescent and plasmonic waves. Our approach is validated by analytically working out the integration in the conventional, more cumbersome, projection approach. Comparison of optical forces on a particle in evanescent and plasmonic fields with previous results shows perfect agreement, thereby further corroborating our approach. As examples of its application, we calculate optical force and torque exerting on particles residing in a plasmonic field, with large particle size where the conventional projection method based on the direct numerical integration is unadapted due to the difficulty in convergence. It is found that the direction of optical torque stays parallel to the direction of spin of optical field for some field polarizations and changes for some other polarizations, as the particle radius R varies.
© 2017 Optical Society of America
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  and arbitrary order , 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 kx, ky, kz 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, , 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 nc 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
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 reads22], defined by 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 . 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) ≤ π , 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 [4, 5] that follows the Rodrigues formula
To derive the PWECs pmn and qmn, we first rewrite the electric field given by Eq. (1) into a form showing explicit the polarization
The complex direction cosines of , viz, (sin α cos β, sin α sin β, cos α), are governed byEq. (3), and the argument of the complex kρ lies between and , 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 p1 and p2 that depict the complex amplitudes of the p and s polarizations,
Next we take advantage of the following expansion of a unit dyadic in terms of VSWFs ,Eqs. (5), (12), and (16), leads to the explicit expressions for the PWECs pmn and qmn,
As a special case, for the wave propagating in the xz plane (ky = 0, sin β = 0 and cos β = 1), the p polarization of incident field corresponds to p1 = 1 and p2 = 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  for computing πmn and τmn with real cos α, see also appendix A of Ref. . 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 , defined by
The renormalized auxiliary functions satisfy 
The recurrence algorithm proceeds as follows. Suppose one has got the values of and for n ≤ l and 0 ≤ m ≤ n for each n. Then one can evaluate and based on Eqs. (20a). With for 1 ≤ m ≤ l − 1 obtained by Eq. (20b) and by Eq. (8), one has all for 0 ≤ m ≤ l + 1. Next, the values of for m = 0 and m > 0 are calculated based on Eqs. (20c) and (20d), respectively. The algorithm starts with the initial values
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 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 . 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.  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 Ex = kz/k, Ey = 0, and Ez = −kx/k, with three Cartesian components of the wave vector given by Eq. (3), and p1 = 1, p2 = 0 are determined by Eq. (14). So the PWECs given by Eq. (17) are readily evaluated with Eq. (13). We adopt nd = 1.33 for water as in . The operating wavelength is fixed at λ = 594 nm so εm = −8.767+1.535i  for gold particle and plate. Our results for non-vanishing components of optical force and torque, normalized by the intensity , as a function of particle size are shown in Fig. 1, together with results from Ref. . 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.
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 nc, with an empirical formula for nc given by31]. Here x = 2πnd R/λ is the size parameter [5, 31]. It is noted that nc 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 nc given in Eq. (23), to ensure an accurate description of complex-k field and convergence in simulation results. It is this large nc 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 . 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 Ey = 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 Fig. 1, the spin direction is determined by the wave vector given in Eq. (22), which means s = −ey. In Figs. 1(b) and 1(d), the optical torque is seen to point towards −ey, 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 nd = 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 kx/k0 = −0.40 + 0.30i, ky/k0 = 0.90 − 0.20i and kz/k0 determined by Eq. (3), so that the spin direction s = 0.71ex + 0.64ey − 0.29ez is based on Eq. (24). The only difference between Figs. 2(a) and 2(b) is the polarizations of electric fields: Ex = −0.68 − 0.10i, Ey = 0.45 + 0.22i for Fig. 2(a) and Ex = 0.87 − 0.36i, Ey = 0.19 − 0.53i for Fig. 2(b), the corresponding Ez 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 , i.e. , and , 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 ky = 0 by a rotation of coordinate system. If, in this new coordinate system, we have either Ey = 0 or Ex = Ez = 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.
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.Eq. (6), it is ready to derive
The expansion of plane wave on P409 of Ref. Eq. (31) gives Eq. (28) and , which is a counterpart of Eq. (29a) in k-space, with . Multiplying Eq. (33) by k2 and integrating over k, one has
It follows from (27a) thatEq. (16a) if one remembers that θ′ and ϕ′ denote the polar and azimuthal angles of k.38], resulting in Eq. (29b) defined in k-space. Similarly, taking the complex conjugate of Eq. (38), multiplying by k2, and integrating over k, one eventually arrives at a form reminiscent of a Fourier transform
Finally, taking the curl of Eq. (38) and noticing , one getsEq. (27c), yields Eq. (16c) in the main text.
Eqs. (32), (38) and (41) correspond to Eqs. (22), (23) and (24) on P172 of Ref. . 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  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. , Almaas and Brevik  evaluate the integrals over ϕ with the modified Bessel function Iv (z) = i−v Jv (iz). Then expressions of expansion coefficients only involve integrations of θ. Meanwhile the associated Legendre function given in this appendix is different from in the main text, and their relationship isEqs. (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 , which is fixed to for simplicity [14, 15], dropped in Eq. (44a). In the main text we put into E0 in Eqs. (1) and (5).
In the Helmholtz equation, An,m and Bn,m are independent of the surface we select , which is already confirmed in the non-evanescent waves  and the evanescent waves . Considering Q0 in Eq. (44b) proportion to b/jn (kb), as a consequence, the integral solutions of Q1–Q3 must be proportional to jn (kb)/b for the completely cancellation of variable b. Observing Q1 and Q2 could be related to Q3 by recursion formulae of associated Legendre functions, we try to work out Q3 in Eq. (44e) first. Using Euler’s formula eix = cos x + i sin x, and replacing Iv (z) with Jv (iz) for Eq. (44e)Q3 is separated into two integrations’ combination. And each part is Bessel function times associated Legendre functions and exponential term, such as the integration below
This integration is a special case of Gegenbauer’s finite integral given by Eq. (1) on P379 of Ref. , and the solution is obtained by using the relationship between associated Legendre functions and Gegenbauer polynomials in the complex plane for m ≥ 0, see Eq. (8.936) on P992 of Ref. . Koumandos makes analytic continuation for −n ≤ m < 0 . Finally we obtainEq. (13). For the real angle 0 ≤ α ≤ π in Eq. (45), Neves proves Eq. (46) by mathematical induction  and Cregg proves Eq. (46) by Gegenbauer polynomial . In fact Eq. (46) is presented historically  and used in Heisenberg chains .
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. , we obtain
Notating the integrand in Eq. (44e)Eqs. (46), (47) and (49) into Eq. (44e), after some algebra we finally get
Emphasizing the range of real angle 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. , and get the relationship between Q1 and Q3 after some algebra,
Using Eqs. (8.733) of Ref. , we get the relation of Q2 and Q3,Eq. (50) has the factor jn (kb)/b and cancels the r-dependent prefactor for Bn,m. However, taking Eqs. (53) into Eq. (44a) could not directly reduce to the factor jn (kb)/b for An,m. We have to transform Eq. (53a) by using the recurrence relations given by Eqs. (8.733) and (8.735) of Ref. ,
The relations between pmn, qmn and An,m, Bn,m are found after comparing incident fields,
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]
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]
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]
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]