1. Introduction

The generalized Lorenz-Mie theory (GLMT) is an extension of the Lorenz-Mie theory [1] for describing the electromagnetic field components of an arbitrary laser beam in terms of spherical harmonic functions [2,3], the coefficients of which being called the beam shape coefficients (BSCs), responsible for correctly modeling the intensity profile of the beam [4]. Their numerical evaluation using quadratures [2] or finite series [4], however, can be a pretty time-consuming, lengthy or awkward task, first because of numerical integrations over the spherical coordinates of the adopted coordinate system, and second because of the inexistence of a single expression, in the latter case. Thus, efficient techniques such as the localized interpretation and, as its natural consequence, the integral localized approximation (ILA), have been developed for calculating on and off-axis BSCs with insignificant difference and much faster than the above mentioned numerical integrations [5–7]. Anyway, regardless of the scheme adopted for evaluating the BSCs, since the development of the GLMT plenty of applications have been benefited by this theoretical methodology (for a resume of the updated subject see, e.g., Refs. [8–10]. and references therein).

The localized approximation developed by Gouesbet et al. is based on the localization principle of van de Hulst [11] and was rigorously justified for on- and off-axis focused Gaussian beams and laser sheets [12,13]. Its improved version, viz., the integral localized approximation, has been used during the years for biomedical optics research such as in optical force and torque calculations exerted on both homogeneous and multilayered spherical dielectric particles in optical trapping systems [14,15]. Even optical forces exerted on hypothetical negative refractive index metamaterials by Gaussian beams have been analyzed with the ILA in the context of the GLMT [16].

In optical trapping experiments, however, there are actually other laser beams of relevant importance, due to their increasing interest because of their properties such as multiple trapping and angular momentum transfer [17–20]. Among those we found the Bessel beams.

Single Bessel beams furnish an easy way for multiple trapping, although their extended focus makes a full three-dimensional trap unachievable. Several methods have been developed for theoretically predicting the behavior of dielectric spherical particles, under the influence of a Bessel beam, in the ray optics regime, where the wavelength λ is much smaller than the diameter d of the particle (λ << d), and in the case where λ >> d, for which a dipole model can be assumed. Also, in the Rayleigh-Gans regime, simple formulas can be used to predict the optical trapping forces with relative success [21]. Even recent works that use the GLMT for studying transverse dynamics of silicon particles under the influence of Bessel beams still formulate their theory based on slightly modified versions of quadratures [22], explicitly emphasizing the time-demanding character of such approach as a function of the particle radius. The inclusion of the ILA for multi-ringed beams in the GLMT may fulfill a bottleneck in this theory and, therefore, really seems to deserve some special attention.

In this way, this paper is devoted to the analysis of the beam shape coefficients in order to describe an ordinary (zero-order) Bessel beam using the GLMT with the integral localized approximation and is organized as follows. In Section 2 we theoretically deduce closed-form solutions of the BSC’s for linearly polarized Bessel beams, which can be automatically extended to circularly polarized beams by applying simple symmetry considerations. In Section 3, numerical results are presented and compared with quadrature methods, which are exact numerical integrations directly derived from the GLMT based on orthogonality conditions. This reveals the efficient and fast character of the ILA algorithm, enabling us to describe, with enough accuracy, the electromagnetic fields of an ordinary Bessel beam. Section 4 shows that the ILA, applied to Bessel beams, makes optical trapping force calculations over spherical micro particles in the Mie regime very simple and accurate when compared to previous works in ray optics, Rayleigh-Gans and Mie regimes. In section 5, our conclusions are presented.

2. The GLMT and the ILA applied to ordinary Bessel beams

Consider an ideal monochromatic zero-order Bessel beam (BB) propagating with speed c and angular frequency ω ₀ parallel to +z, with its optical axis displaced ρ ₀ = (x ₀,y ₀) from the z-axis and making an angle ϕ ₀ relative to the x-axis, according to Fig. 1 . Although we could have shifted the beam in a three-dimensional fashion, we must take into account that, as long as the ideal definition is considered and that the corresponding shift in z be constricted to the focal length (thus keeping the radial Bessel profile still valid), this implies in a more complete but unnecessary mathematical formulation. It is well known that the longitudinal and transverse wave numbers for a BB are given by k_z = kcosθ_a and k_ρ = ksinθ_a, respectively, θ_a being the associated axicon angle and k = ω ₀/c the wavenumber [23]. If this beam is linearly polarized in x or in y, we can write its electric field as

 

Fig. 1 Geometrical description of an ordinary Bessel beam propagating parallel to + z (out of the page). The optical axis makes an angle ϕ ₀ relative to the x-axis and is displaced ρ ₀ from the origin O of the coordinate system.

Download Full Size | PDF

where J ₀(.) is the ordinary Bessel function and E ₀ is the electric field strength. A factor exp(iω0t) has been omitted for convenience. The field profile in Eq. (1) may be considered to be valid as long as β = sinθa/(1 + cosθa) << 1 [22]; i.e., the results presented here will represent good descriptions of real experiments if the axicon angle is small enough that a paraxial approximation can be assumed.

In the framework of the GLMT, there are several ways in which the beam-shape coefficients associated with Eq. (1) can be computed, all of them demanding solely a previous knowledge of the radial components E_r and H_r of the beam as functions of spherical coordinates (r,θ,ϕ) relative to the origin O in Fig. 1 (in the case of optical trapping, for instance, O would coincide with the center of a spherical scatterer, as will be the case in Section 4). The radial component E_r of the electric field for this BB is readily found by putting a multiplicative factor sinθcosϕ (sinθsinϕ) to the right of Eq. (1) in the case of an x-polarized (y-polarized) beam, and replacing ρ by rsinθ and z by rcosθ. Then, for the radial component H_r of the magnetic field, we recall Faraday’s law in the frequency domain and in its differential form, so that we can write E_r and H_r, after some simple algebra, as

where $H_0=E_0/\eta =E_0\sqrt{\epsilon /\mu }$ and η is the intrinsic impedance of the propagating medium, ε and μ being, respectively, the permittivity and permeability of this medium. The slashes represent either an x- or a y-polarized beam, with its corresponding sinϕ or cosϕ factor.

From Eqs. (2) and (3), the beam shape coefficients $g_{n,TE}^m$ and $g_{n,TM}^m$ in the integral localized approximation are readily found from the following relations [7]:

where n and m are integer numbers ranging from 1 ≤ n < ∞ and –n < m < n, respectively, corresponding to a particular spherical harmonic function $Y_n^m\left(\theta ,\varphi \right)=P_n^m\left(\cos \theta \right)\exp \left(im\varphi \right)$, $P_n^m\left(\cos \theta \right)$ being associated Legendre functions. In the integral localized approximation description of $g_{n,TE}^m$ and $g_{n,TM}^m$, the localization operator $\widehat{G}$changes the factor kr to (n + 1/2) and θ to π/2, while $Z_n^m$are normalization factors that read as [7]

By substituting Eq. (2) into Eq. (5) for m = 0, we find

where we have used the normalization factor $Z_n^0$ from Eq. (6). To analytically solve the above integral, we perform the following expansion for a zero-order Bessel function [24]

ε_p being 1 if p = 0 and 2 otherwise. Thus, letting $\varpi =\sin {\theta }_a\left(n+1/2\right)$ and $\xi ={\rho }_{0}k\sin {\theta }_a$, we have from Eq. (8)

which will always be zero whenever ρ ₀ = 0 (on-axis case), because $J_1\left({\rho }_{0}k\sin {\theta }_a\right)=0$.

For m ≠ 0, the analysis is similar to the previous one and we have, from Eqs. (5), (7) and (9),

Notice that, in deriving Eqs. (10) and (11), we have made use of some basic orthogonality conditions for trigonometric functions. Also, because J _{| m | ± 1}(ρ ₀ ksinθ_a) = 0 whenever ρ ₀ = 0 and m ≠ 1, as expected for an on-axis (z-axis) beam, the only non-zero beam-shape coefficients for ρ ₀ = 0 are $g_{n,TE}^m$ and $g_{n,TM}^m$, both of which with a magnitude of $\left(1/2\right)J_0\left(\sin {\theta }_a\left(n+1/2\right)\right)$ that decreases with increasing n.

Similarly, the evaluation of the BSCs $g_{n,TE}^m$ is analogous to the previous derivation, Eq. (4) being used instead of Eq. (5). Alternatively, one may also realize the similitude between Eqs. (2) and (3) with only an interchange between sines and cosines inside the slashes, thus leading to an immediate determination of the $g_{n,TE}^m$:

These are the beam-shape coefficients necessary to fully describe an x- or y- polarized (on- or off-axis) ordinary Bessel beam in the GLMT using the integral localized approximation. They may be greatly simplified if, for example, ϕ ₀ = 0 for x-polarization, or ϕ ₀ = π/2 together with an y-polarized BB, or other analogous particular cases. Notice that we can readily extend the set of equations (10)–(13) to circularly polarized beams by means of simple symmetry relations and writing each of the new TM and TE BSC’s as functions of both TM and TE BSC’s in Eqs. (10)–(13) [15].

3. Numerical results for the BSC’s of an ordinary Bessel beam

In the framework of the generalized Lorenz-Mie theory, the beam shape coefficients can be primarily determined by means of finite series or using quadratures. In the latter case, double and triple integrations over the spherical coordinates (r,θ,ϕ) can be a really time-consuming task, as already pointed out by previous authors [5–9,12,13]. For double integration, if the radial components of the electromagnetic fields of a laser beam are given, we can numerically evaluate the BSC’s $g_{n,TE}^m$ and $g_{n,TM}^m$ by solving the relations

$P_n^m\left(\cos \theta \right)$ being the associated Legendre polynomials, R = kr an auxiliary coordinate and j_n(R) are spherical Bessel functions of integer order n. If E_r and H_r represent an exact solution to Maxwell’s equations, then R is arbitrary because, in such a case, the integrals over the spherical angles θ and ϕ are proportional to j_n(R)/R [12,13,25,26]. But some care must be exercised in choosing a specific R when numerically evaluating Eqs. (14) and (15). It is recommended, because of the van de Hulst principle, and we have adopted it here, to impose R = n + 1/2 [11]. This holds because the most significant contribution of the integrands in Eqs. (14) and (15) in the case of a zero-order Bessel beam also happens around R = n + 1/2, as one can verify by plotting several curves of them and performing the same analysis as those for a focused Gaussian beam [25]. By recalling orthogonality relations of the spherical Bessel functions, the above equations may be written as triple integrals, thus getting rid of the R dependence, with a corresponding undesirable computational cost [12,13].

To compare the numerical calculation performance of the BSCs obtained in section 2 by means of the integral localized approximation with double (Eqs. (14), (15)) and triple integration, a Fortran code was developed and is available under request. The code was then run on a PC (AMD 64 Athlon processor 3500+, 805MHz, with 768Mb of RAM) with windows XP. Table 1 compares the coefficients $g_{n,TM}^1$ calculated via integral localized approximation (called ILA) and quadratures with double (F1) and triple (F2) integration for an x-polarized Bessel beam with wavelength λ = 1064 nm and an axicon angle θ_a = 0.0141 rad (thus, with a transverse spot of Δρ = 2.405v/(ω ₀sinθ_a) = 28.89 μm, v being the velocity of the wave in the propagating medium), with its optical axis coincident with the z-axis, i.e., ρ ₀ = 0, in air. Only the non-zero terms calculated by the ILA method are shown (|m| = 1). Analogous results can be achieved for off-axis beams. For m < 0, the coefficients are immediately calculated by considering symmetry relations [27], although in our code we have preferred to calculate them directly by applying the set of Eqs. (10)-(13). This obviously requires more computational time for applications where these (m < 0) coefficients are needed, such as in optical force calculations for optical trapping.

Table 1. Beam-shape coefficients $g_{n,TM}^1$ evaluated by quadratures (methods F1 and F2) and using the integral localized approximation ILA for an on-axis (ρ₀ = ϕ₀ = 0) zero-order Bessel beam with λ = 1064 nm and θ_a = 0.0141 rad^a

View Table | View all tables in this article

The number of steps for numerically evaluating the BSCs of Table 1, using methods F1 and F2, was set to 200 for integrating both in ϕ and θ. Furthermore, for method F2, the integral over R was performed from 0 to 200 for n = 1 and up to 10, from 0 to 500 for n = 15 and 20, and from 0 to 1500 for n = 50 and 100, all of them with 500 integration points. This ensures a reasonable convergence of both quadrature methods [25].

As expected, $g_{n,TM}^m$ = 0 whenever m ≠ 1 for ρ ₀ = 0.This happens because, as already mentioned in the previous section, J _{| m | ± 1}(ρ ₀ ksinθ_a) = 0. It is easily verified that, for off-axis BBs, significant values for the BSC’s appear for m ≠ 0, as also expected [12,25,27].

From Table 1, one can see that the magnitudes of $g_{n,TM}^m$ (or, equivalently, $g_{n,TE}^m$, as $g_{n,TE}^m=-i{g}_{n,TM}^m$ for ρ ₀ = 0) from the ILA method are very close to those obtained by means of quadratures, when we take into account only the most relevant BSCs in magnitude (low n). In fact, for all the BSCs presented, whenever the relative magnitude $g_{n,TM}^1/g_{1,TM}^1$is small there is a dissimilitude between ILA, F1 and F2 methods (for example, ${g_{100,TM}^1|}_{ILA}$ = 0.280726, ${g_{100,TM}^1|}_{F1}$ ≈0.279240 and ${g_{100,TM}^1|}_{F2}$ ≈0.266963, representing a relative difference of 0.53% and 5.16%, respectively) which, overall, does not contribute significantly to the description of the field intensity profile of the Bessel beam with the parameters chosen. But we must point out that this does not represent a failure of the ILA method for high n, being solely a matter of raising the number of numerical points in θ, ϕ (for F1 and F2) and R (for F2) to ensure an adequate numerical convergence. For example, using 500 points in θ and ϕ leads to ${g_{100,TM}^1|}_{F1}$ ≈0.280659 and ${g_{100,TM}^1|}_{F2}$ ≈0.268462, and increasing n necessarily imposes additional refinements in our numerical integration for F1 and F2.

As we have shown that the ILA can furnish the values of the most significant BSCs with great accuracy, we now turn our attention to the most outstanding advantage of using the integral localized approximation, viz., the computation time of these coefficients. Due to the absence of numerical integrations, calculating $g_{n,TE}^m$ and $g_{n,TM}^m$ using the ILA is much faster than using the quadrature methods F1 and F2. Table 2 shows the average time taken for numerically evaluating the BSCs presented in Table 1. These values were obtained by calculating each ten times and then taking the time average. Notice that, as n increases, so does the average time for F1 and F2. This happens in our simulations because our subroutine uses recursive relations for the spherical Bessel functions, so that, for example, for n = 200, 200 iterations are performed to compute j ₂₀₀(R) in these methods.

Table 2. Elapsed time (in seconds) for computing the BSC’s $g_{n,TM}^1$ of an on-axis (ρ₀ = ϕ₀ = 0) zero-order Bessel beam with λ = 1064 nm and θ_a = 0.0141 rad^a

View Table | View all tables in this article

Thus, the integral localized approximation proves to be a computationally efficient method for computing the beam-shape coefficients necessary to describe an ordinary BB in the framework of the generalized Lorenz-Mie theory. Although this approximation is known since de 90’s, it is noteworthy and remarkable that so few laser beams have ever been described with it [28–30]. It can be shown that good results are also achieved if we replace the x-polarized zero-order Bessel beam by a circularly polarized (xy plane) with the same parameters.

Now, suppose that we are given a set of beam-shape coefficients $g_{n,TE}^m$ and $g_{n,TM}^m$, regardless of the numerical method used to evaluate them. The electric field components in the GLMT can be written as [13]:

with similar expressions for the magnetic field components. In Eq. (16), ${\pi }_n^{\left|m\right|}\left(\theta \right)=\sin \theta P_n^m\left(\cos \theta \right)$ and ${\tau }_n^{\left|m\right|}\left(\theta \right)=\left(1/\sin \theta \right)d{P}_n^m\left(\cos \theta \right)/d\theta$. Furthermore, consider an x-polarized zero-order Bessel beam propagating along + z. The dominant electric field component E_x is written, based on the above equation, as $E_x\left(x,y,z\right)=E_r\sin \theta \cos \varphi +E_{\theta }\cos \theta \cos \varphi -E_{\varphi }\sin \varphi$ but, but, if we further assume an observation point ρ ₀ along the x-axis, we can express the electric field E_x in the same way as Ref. [13]:

Figures 2(a)-(d) are plots of the intensity profiles of Eq. (17) for ϕ = 0 and ρ ₀ = 0, 30, 60 and 90 μm for different ranges of n and with m_max = 15. There is a clear compromise between x and the maximum n entering superposition Eq. (16) for E_r from which a good prediction of the Bessel profile is achieved. For the parameters chosen, a maximum n = 700 is sufficient, as long as the physical analysis is limited to 0 < x < 120 μm. If x is to be kept below approximately 60 μm (i.e., in cases where the physical analysis does not demand the knowledge of the field for distances above x = 60 μm), then 1 ≤ n < 500 could be used to make field calculations faster. Note that this is in contradiction with some previous analysis for focused Gaussian beams, where off-axis descriptions using the ILA may possess serious limitations due to the approximation of the beam itself, i.e., the truncation of the Davis’ series description of this type of beam [13]. In fact, the integral localized approximation, as originally developed, is based on a Gaussian description of order L– or L [6]. For higher ρ ₀’s, m_max may also have to assume higher values.

 

Fig. 2 E_x-field intensity profile for a Bessel beam with λ = 1064 nm, θ_a = 0.0141 rad and displaced (a) ρ ₀ = 0, (b) ρ ₀ = 30 μm, (c) ρ ₀ = 60 μm and (d) ρ ₀ = 90 μm along x. The accuracy of the ILA for increasing x depends on increasing the number of BSCs entering superposition (18) and (19). In all cases, m_max = 15.

Download Full Size | PDF

Thus, for ordinary Bessel beams, the ILA is capable of reproducing the original electric field with great accuracy and reliability. We have also analyzed the intensity profile and relative difference for all the other electromagnetic field components using the ILA formulation. All of them can be reproduced by an adequate choice of n_max and m_max. Similar results were then obtained.

Finally, Fig. 3 shows the three-dimensional electric field intensity profiles for an x-polarized BB with λ = 500 nm and spot Δρ = 10.0 μm for the on-axis case, whereas Fig. 4 represents an x-polarized BB with λ = 532 nm and Δρ = 2.336 μm for an off-axis (ρ ₀ = 20 μm and ϕ ₀ = 0) situation. Excellent agreement is achieved with the ideal situation (paraxial approximation) for the transverse fields and also for the longitudinal fields when compared with vector Bessel beams presented by previous authors [22].

 

Fig. 3 (a) 3D and (b) 2D Ey-field intensity profile for an on-axis x-polarized Bessel beam with λ = 500 nm, Δρ = 10.0 μm (θa ≈ 1.1°). The Fortran code was run by imposing nmax = 1800 (non-zero BSCs occur only for |m| = 1). (c) and (d) are the corresponding Ez-field.

Download Full Size | PDF

 

Fig. 4 .a) 3D and (b) 2D Ey-field intensity profile for an off-axis (ρ0 = 5 μm and ϕ0 = π/2) x-polarized Bessel beam with λ = 532 nm, Δρ = 2.336 μm (θa = 5°, representing a limiting angle for the paraxial approximation). The Fortran code was run by imposing nmax = and mmax = 15. (c) and (d) show the corresponding Ez-field.

Download Full Size | PDF

4. Application to optical force calculations

The success of the integral localized approximation for describing a BB can be used advantageously in biomedical applications, such as in the determination of the optical forces exerted over biological particles under the influence of such beams.

As a first example, let a dielectric spherical particle with radius a and refractive index n_p be immersed in a liquid medium with refractive index n_m, so that the relative refractive index is given by n_rel = n_p/ n_m. Consider that the product k_ρa = 3.5, i.e., for a Bessel beam with λ = 1064 nm and θ_a = 0.0141 rad (thus with a transverse wave number k_ρ = ksinθ_a = (2π/λ)sinθ_a = 8.33x10⁻⁴ m⁻¹) propagating along the liquid medium, our particle has a radius a ≈42.04 μm. The beam is again assumed to be x-polarized. These parameters were chosen in such a manner as to fit those of Ref. [21]. In the Mie regime, expressions for the x-component of the radiation pressure cross-section C_{pr , x} (which is intrinsically related and directly proportional to the x-component of the optical force) can be found elsewhere and will not be reproduced here [25].

Figure 5 reveals the normalized axial force exerted on the particle when the BB is shifted along x for n_rel = 1.1. The force intensity profile for k_ρa = 0.1, i.e., a ≈1.20 μm, is also shown for comparison. In this normalized plot, F_{x ,max| kρa = 3.5}/ F_{x ,max| kρa = 0.1} = 2858.9. To reproduce these slopes, we have set n_max = 700 and m_max = 15, thus ensuring the adequate description of the Bessel beam profile, as already pointed out in previous sections.

 

Fig. 5 Normalized force profile exerted on a dielectric simple particle of relative refractive index nrel = 1.1. The field intensity of the zero-order Bessel beam is also shown as a solid line. All parameters of the incident beam are the same as those of Fig. 2. Positive Fx/Fmax means an attractive force towards the optical axis of the beam.

Download Full Size | PDF

By comparing Fig. 5 with Figs. 2(a) and (b) of Ref. [21], we see that the ILA can be successfully applied in optical force calculations, predicting the same stable equilibrium positions and magnitudes of other previous methods.

Notice, however, that according to Fig. 2, as long as the particle is assumed to be fixed at x = 0, n_max can be significantly reduced because at this point few n-iterations are necessary in calculating C_{pr , x}, to achieve an adequate description of the fields and, consequently, of the optical forces. In this way, for plotting Fig. 5, for instance, we could have safely chosen n_max = 200 and still get the same force profiles as shown.

As a second example, Fig. 6 shows the x component of the radiation pressure cross-section profile C_{pr , x} for silicon spheres (n_p = 1.4496) of four different radii immersed in water (n_m = 1.33). The incident beam is a Bessel beam with Δρ ≈2.35 μm (θ_a ≈7.51°) and a wavelength λ = 802.7 nm, both in water. These parameters were chosen to best fit the theoretical data of Ref. [22], so that negative C_{pr , x} means an attractive force towards the optical axis of the beam.

 

Fig. 6 Radiation pressure cross-section C_pr,x (solid) for an x-polarized Bessel beam displaced along x using the ILA. The beam has λ = 802.7 nm and Δρ ≈2.35 μm in water (n_m = 1.33). The beam intensity is shown as a dotted line. The silicon spheres have a refractive index n_p = 1.4496 and radii a = 1.15 μm (a), 2.15 μm (b), 2.50 μm (c) and 3.42 μm (d). Points of stable equilibrium are close to those predicted in Ref. [21], where a quadrature scheme was adopted for numerically implementing the GLMT.

Download Full Size | PDF

Comparing our results with the above mentioned reference, we see approximately the same points of stable equilibrium. As for the simulation time, each C_{pr , x} curve took, in average, 512 s for n_max = 1300, m_max = 20 and 100 radial positions (we point out that, if we had implemented symmetry relations for calculating the BSCs for m < 0 in our code, or adopted a lower n_max, as pointed out above in deriving Fig. 5, this time would have been reduced further). These upper limits for n and m are again necessary to ensure that the set of BSCs reproduce the electromagnetic fields for this particular BB within the range 0 ≤ x ≤ 20 μm with negligible difference, thus providing the correct picture for the radiation pressure cross-section. This simulation time is radically in contrast to the elapsed time observed in Ref. [22], where it explicitly depended on the radius of the silicon spheres and computational calculations could last more than 3.5 hours for a silicon sphere of radius a = 10 μm (for the same PC configuration, this would represent a multiplicative factor of 25.2. In practice, however, this factor is even higher, as numerical simulations of Ref. [22]. were performed on a higher performance PC). Obviously, here it is the number of BSCs which ultimately establishes the average time, as increasing n_max and m_max turn the evaluation of Eq. (16) and C_{pr , x} more time consuming and, furthermore, more and more BSCs have to be calculated.

In fact, the use of the ILA in the generalized Lorenz-Mie theory provides an efficient tool not only for force calculations, but also for other physical entities such as the electromagnetic field components themselves, as already shown in this paper, or angular momentum, scattering and absorbing cross-sections.

5. Conclusions

For the first time in the literature, the beam-shape coefficients for ordinary Bessel beams were numerically evaluated using the integral localized approximation in the context of the generalized Lorenz-Mie theory. Closed-form solutions were presented, extending the range of applicability of the ILA to beams other than Gaussian and laser sheets.

The method presented here furnishes almost the same results of those based on quadratures for the magnitude of the BSCs necessary to adequately describe the electromagnetic fields of zero-order Bessel beams. But the fundamental advantage lies on the incredibly reduced elapsed time to compute each of these coefficients. This ultimately justifies the adoption of this method for subsequent works.

We have also applied our formalism to compute optical forces in optical trapping systems with incident ordinary Bessel beams, comparing their profiles with those already obtained by previous authors. The similitude of results, together with the lower computational time, reinforce the reliability of our approach and strongly suggest the extension of the integral localized approximation to a wider variety of laser beams, including BB of higher orders. For this particular class of multi-ringed beams, however, it is expectable that some of the symmetry relations in the original development of the ILA method will not hold true. Of course, this can certainly pose difficulties in our attempt to eliminate the angular integration presented in the formulation of the ILA for the BSCs, thus making closed-form solutions more challengeable.

Extending the range of applicability of the ILA certainly would represent a significant gain in the study of optical properties other than forces. We could use this method, for example, in optical torque calculations or simply to find absorbing, scattering and extinction cross-sections for a particle under the influence of a particular laser beam such as Laguerre-Gaussian, for example.