Generation of Bessel beam sources in FDTD

Zhefeng Wu; Yiping Han; Jiajie Wang; Zhiwei Cui

doi:10.1364/OE.26.028727

1. Introduction

Bessel beams were firstly introduced by Durnin [1] as a non-diffracting electromagnetic wave described by Bessel function of the first kind. Because of its unique properties, such as nondiffraction [2], self-reconstruction or self-healing [3,4], Bessel beams have attracted widespread attention, and promising applications can be found in various fields, including laser materials processing [5], optical manipulation [6], optical acceleration [7], nanolithography [8], imaging [9], and others. Although ideal Bessel beams cannot be generated due to its infinite lateral extent and energy, quasi-Bessel beams can be experimentally generated in microwave [10], terahertz wave [11], and optical regime [12].

Interactions between Bessel beams and complex objects are basic and common topics to be studied in various field. Among the theories for dealing with scattering of objects by structured beams, the generalized Lorenz–Mie theory (GLMT) [13] is a powerful tool which has been widely used. Electromagnetic non-resonance and resonance scattering of zero-order and high-order Bessel beams by a dielectric sphere has been investigated by Mitri et al. [14–16]. Scattering of a homogeneous spheroidal particle was studied by Han [17]. Expanding the GLMT into Debye series, Li et al. [18] analyzed the radiation pressure force exerted on a biological cell which was modeled by a multilayered sphere. Light scattering of an eccentric sphere illuminated by a Bessel beam was analyzed by Wang et al. [19]. Recently, radiation force cross-sections and optical spin torque induced by vector Bessel beams have also been studied [20–22]. Nevertheless, within the framework of GLMT, only regularly shaped particles, such as sphere, cylinder, spheroid, and their composite shapes, can be dealt with. This is because the GLMT is based on the use of separation of variables method, which is of limited use for irregularly shaped particles. Thus various semi-analytical methods and numerical methods have been developed. Yang et al. [23] studied the scattering properties of nonspherical particles using the multilevel fast multipole algorithm (MLFMA). Cui et al. [24] investigated scattering problems involving core-shell structure composite particles by applying the surface integral equations (SIEs). T-matrix method is developed to investigate the acoustical Bessel beam scattering from a rigid finite cylinder and a rigid spheroid by Gong et al. [25,26]. So far, although the finite-difference time-domain (FDTD) method has been extensively used in the prediction of scattering properties of complex particles in the case of plane wave illumination, few research has been reported on the cases of structured beams illumination. The FDTD method [27] discretizes time-dependent Maxwell’s equations by using central-difference approximations on the spatial mesh and solves the resulting finite-difference equation in a leapfrog manner. As a time domain method, FDTD can obtain a broadband result within a single simulation. Besides, specifying the material at every point within the computational domain is allowed in FDTD simulation, leading to easily modeling of various material. To provide an access to solve the problem of Bessel beams interact with complex material, such as plasmonic, nanostructure and metamaterial, this paper focuses on the study of generating Bessel beams sources in FDTD, which is one of key issues required to be addressed as one intends to apply FDTD to deal with Bessel beam scattering of complex objects.

To introduce an incident wave into the simulation domain in FDTD, two common ways, i.e. total-field/scattered-field (TF/ST) and pure scattered-field method [27], can be applied. TF/ST method is an application of the electromagnetic field equivalence principle. In this method, a fictitious closed surface divides computation domain into two regions: total field region and a scattered field region. The incident field is generated and confined in total field region by applying equivalence electric and magnetic current sources on this surface. Using this method, Çapoğlu [28,29] introduced a focused light pulse into FDTD using the superposition of plane wave spectrum. The focused beam can also be excited using a single plane source, which can be viewed as a variation of TF/SF method [30]. On the other hand, pure scattered field method, which based on the linearity of Maxwell’s equation, decomposes total fields into known incident fields and unknown scattered fields. Incident fields at all grid points where there are not free-space are evaluated using analytical expression at each time step. Thus, a variety of shaped beams can be generated using pure scattered field method. It has been applied to the investigation of far-field characteristics of a laser beam scattered by dielectric particles [31], and calculating of optical forces of tightly focused beams on microscopic particles [32]. Compared to the TF/SF method, pure scattered field method is a more straightforward approach to generate shaped beams sources. However, as has been outlined in [27], TF/SF method has advantages in computational efficiency, relatively simple programming of targets, and wide dynamic range. So far, little efforts were devoted to the generation of Bessel beams in FDTD. Using scattered-field FDTD, Chen et al. [33] studied the scattering of complex particles by Bessel beams. In [34] and [35], based on the superposition principle, the scattering results of Bessel beams by objects were obtained in the post-processing by superposing each partial plane wave result, which is much time-consuming. In this paper, a straightforward approach is presented for the generation of Bessel beam sources in FDTD by using the TF/ST method based on the angular spectrum representation (ASR) of the incident Bessel beam.

The rest of the paper is organized as follows. In Section 2, we briefly describe the ASR of Bessel beams which is the theoretical basis of our approach. In Section 3, details concerning the generation of Bessel beam sources in FDTD are presented. In Section 4, the validity and accuracy of the present method are verified. Exemplify results concerning the interactions of Bessel beams with small particles, including far-field, internal and near field, are presented. Conclusions are given in Section 5.

2. Angular spectrum representation of Bessel beams

Scalar description of Bessel beam was first used to analyze the beam properties. It can be obtained by solving the scalar wave equation. For applications under the paraxial condition, scalar Bessel beam gives satisfactory results. But for nonparaxial cases or applications that polarized electromagnetic should be considered, a vectorial treatment is required. One approach to establishing a vector Bessel beam is to solve the vector wave equation directly. With the introduction of magnetic vector potential or Hertz vector potential [36,37], researchers derived a series of vector Bessel beams of arbitrary order in different polarization states, including linear, circular, radial, azimuthal polarizations, transverse electric (TE) and transverse magnetic (TM) modes. Another approach is to represent the vector Bessel beam as a superposition of the beam’s angular spectrum which consists of plane waves with wavevectors on the conical surface [38]. This approach is named by angular spectrum representation (ASR) method and the Bessel beam with arbitrary order and polarization can be described using an integral formalism. Recently, Bessel beams that derived by these two approaches mentioned above was found to have a same functional dependence which led to a general mathematical description [37,39]. The ASR, which is the basis of our approach to generating a Bessel beam in FDTD, is used in this paper.

According to the ASR, assuming the origin of a spherical coordinate system Oxyz is located at the focal center of a lens, the field of an incident laser beam propagating in z-direction can be expressed as [40]:

\begin{matrix} E (r) = \frac{i k f_{0} e^{- i k f_{0}}}{2 π} \int_{θ}^{θ_{m a x}} \int_{ϕ = 0}^{2 π} E_{p w} e^{- i k \cdot r} \sin θ d θ d ϕ, \\ B (r) = \frac{1}{- i ω} \nabla \times E (r) \end{matrix}

where f₀ is the focal length of the lens. Wave vector is k = (k sinθ cosϕ, k sinθ sinϕ, k cosθ), and ϕ and θ are the azimuthal angle and the polar angle, respectively. Considering an ideal L^th-order Bessel beam, the angular spectrum function E_pw in Eq. (1) can be expressed as:

E_{p w} = Q (θ_{0}, ϕ) E_{p w 0} (θ_{0}, ϕ) e^{- i L ϕ} \frac{δ (θ - θ_{0})}{\sin θ},

in which δ(∙) is the Dirac delta function, θ₀ is the half cone angle of the Bessel beam, and E_pw0 is the amplitude of the electric field of a plane wave propagating along the direction determined by ϕ and θ₀. The vector Q, which describes the polarization, is expressed as [41]:

Q (θ, ϕ) = [\begin{matrix} p_{x} (\cos θ \cos^{2} ϕ + \sin^{2} ϕ) - p_{y} (1 - \cos θ) \sin ϕ \cos ϕ \\ - p_{x} (1- \cos θ) \sin ϕ \cos ϕ + p_{y} (\cos θ \sin^{2} ϕ + \cos^{2} ϕ) \\ - p_{x} \sin θ \cos ϕ - p_{y} \sin θ \sin ϕ \end{matrix}],

in which (p_x, p_y) are the parameters that determine the polarization of the Bessel beam. Specifically, (p_x, p_y) = (1, 0), (0, 1), (1, i), (1, -i), (cosϕ, sinϕ), and (-sinϕ, cosϕ) correspond to the polarization that is x-, y-, left circular, right circular, radially, and azimuthally polarized, respectively [16,41]. Substituting Eq. (2) into Eq. (1), omitting the constant term, the electric field of the Bessel beam reduces to a one-dimensional integral:

E (r) = \int_{ϕ = 0}^{2 π} E_{p w 0} Q (θ_{0}, ϕ) e^{- i L ϕ} e^{- i k \cdot r} d ϕ

Before applying Eq. (4) in FDTD simulation, the integral needs to be discrete and expressed in the time domain. A detail discussion is given in Section 3.

3. FDTD implementation of Bessel beam sources

The method we used to create an incident Bessel beam in FDTD is total-field/scattered-field boundary (TF/ST) method. It is an efficient technique to compute scattered field in scattering problem and commonly used to inject a plane wave source by applying a one-dimension (1D) to three-dimension (3D) field projection. Since only a finite number of plane waves can be calculated in practical FDTD simulation, it is necessary to discrete the integral in Eq. (4) as:

E (r, t) = \sum_{n} C_{n} E_{p w 0} Q (θ_{0}, ϕ_{n}) e^{- i L ϕ_{n}} e^{i (ω t - k \cdot r)}, n = 0, 1 ... N,

in which ϕ_n and C_n are numerical integration points of azimuthal angles and weights, e^iωt is the time factor. Many methods can be used to get an approximation to the one-dimensional integral, such as rectangle rule, trapezoidal rule, Simpson's rule, and so on. For simplicity, we choose trapezoidal rule [42] as our one-dimensional quadrature scheme. In the trapezoidal rule, subintervals have the same length, and region under the graph of the integrand is approximated by trapezoids. We can determine ϕ_n and C_n as:

\begin{matrix} ϕ_{n} = 2 π n / N, n = 0, 1 ... N, \\ C_{n} = {\begin{matrix} π / N, n = 0, N \\ 2 π / N, o t h e r s \end{matrix} \end{matrix}

For a zero-order Bessel beam with L = 0, e^-iLϕn equals unity. In this case, it is obvious that the exponential term e^i(ωt-^k^∙^r⁾ in Eq. (5) is a harmonic plane wave which can be represented by 1D FDTD simulation. If vector Q is real, e.g. (p_x, p_y) = (1, 0), (0, 1), (cosϕ, sinϕ), and (-sinϕ, cosϕ), each plane wave with different incident ϕ_n would have the same phase shift at reference sphere surface of the aplanatic system. Under this circumstance, a single 1D FDTD grid is able to represent all the element plane waves in practical simulations. For time domain simulation, this plane wave in the frequency domain can be directly replaced by a broadband incident waveform.

However, for circular polarized Bessel beams whose vector Q is complex, and for higher-order Bessel beams where the term e^-iLϕn is no longer equal to unity, the FDTD implementation becomes complicated. This is because an equivalent additional phase factor should be added to the partial plane waves in these cases. According to the time shifting properties of the Fourier transform h(t-tꞌ) ↔ e^-i^2π^ftꞌF(f), the additional phase factor turns out to be time shift tꞌ in the time domain. That means each wave front of the plane wave in Eq. (5) has different time delay with different wave vectors. An effective approach to deal with this is to allocate a 1D FDTD grid for each plane wave. Although it slightly increases the required memory, it provides an easy access to separately control the retarded time of each angular spectrum element.

Considering a linear polarized L^th order Bessel beam, the time shift tꞌ_n of plane wave spectrum with incident azimuthal angle ϕ_n can be calculated as tꞌ_n = Lϕ_n/2πf, which shows tꞌ_n varies with frequency f. The Eq. (5) can be written as:

E (r, t) = \sum_{n} C_{n} E_{i} (θ_{0}, ϕ_{n}, t - t'_{n} - k \cdot r / ω) Q (θ_{0}, ϕ_{n}), n = 0, 1 ... N .,

The Eq. (7) is the finial expression used in our FDTD program to generate a Bessel beam source. For time domain simulation, the harmonic plane wave and phase term in Eq. (5) are replaced by a modulated Gaussian pulse E_i, which is expressed as:

E_{i} (t) = - \cos (2 π f_{0} t) \exp [- 4 π {(t - t_{0})}^{2} / τ^{2}],

with τ = 6.68fs, f₀ = 5THz, and t₀ = 1.5τ. The full width at a tenth of maximum of the spectrum amplitude of this broadband pulse covers the wavelength ranging from 400nm to 1200nm. It should be emphasized that the circular polarized or higher-order Bessel beams generated in this way cannot be used in broadband simulations since the time shift tꞌ_n of plane wave spectrum varies with frequency f. Although the approach is not suitable for broadband calculation, it still gives a solution to generate higher-order Bessel beams in FDTD.

In practice, numerical dispersion and anisotropy is a crucial drawback of FDTD method which should be taken into consideration. Traditional method FDTD(2,2), which use second-order center difference to approximate partial space and time derivative, introduces error due to its significant numerical dispersion and anisotropy, especially in long-term or large-scale simulation. Using a finer grid can reduce this undesired effect but consumes more computation resources. Higher order method FDTD(2,4) [27], which uses fourth-order center difference to approximate space derivative, has a better performance. Compared to the second-order difference, the fourth-order difference achieves lower numerical dispersion for a coarse grid case, showing great potential in large-scale electromagnetic simulation. In addition, from FDTD(2,2) to FDTD(2,4), without changing the leap-frog time stepping approach, only the iteration part of the code needs to be modified, leading to an easy implementation. As a consequence, higher order method FDTD(2,4) is applied in our code. Assuming electric conductivity σ = 0, the updated equation of the x-component of the electric field at n + 1 time step is expressed as follows:

\begin{array}{l} E_{x} |_{i + 1 / 2, j, k}^{n + 1} = E_{x} |_{i + 1 / 2, j, k}^{n} + \frac{Δ t}{ε (m)} (\frac{27}{24} \frac{H_{z} |_{i + 1 / 2, j + 1 / 2, k}^{n + 1 / 2} - H_{z} |_{i + 1 / 2, j - 1 / 2, k}^{n + 1 / 2}}{Δ y} - \frac{1}{24} \frac{H_{z} |_{i + 1 / 2, j + 3 / 2, k}^{n + 1 / 2} - H_{z} |_{i + 1 / 2, j - 3 / 2, k}^{n + 1 / 2}}{Δ y}, \\ - \frac{27}{24} \frac{H_{y} |_{i + 1 / 2, j, k + 1 / 2}^{n + 1 / 2} - H_{y} |_{i + 1 / 2, j, k - 1 / 2}^{n + 1 / 2}}{Δ z} + \frac{1}{24} \frac{H_{y} |_{i + 1 / 2, j, k + 3 / 2}^{n + 1 / 2} - H_{y} |_{i + 1 / 2, j, k - 3 / 2}^{n + 1 / 2}}{Δ z}) \end{array}

in which ε(m) denotes the average permittivity at grid point (i + 1/2, j, k). In our simulations, grid spacing is ∆x = ∆y = ∆z = 25nm and Cartesian FDTD mesh is used. Time step is chosen as ∆t = 0.45∆x/c in order to fulfill the Courant stability criterion in higher order FDTD, where c is the velocity of light in vacuum. Computation domain is truncated by a uniaxial perfectly matched layer (UPML).

4. Simulation results and discussions

In this section, the validity and accuracy of the present method in this paper are verified by comparing the reconstructed field in FDTD with the original field. Far-field scattering results concerning the interactions of Bessel beams with small particles are then presented. Internal and near-surface field distributions for a two-layer hemisphere particle illuminated by Bessel beams are also displayed to show the potential in solving scattering problems of complex particles.

Amplitude distributions of the x-component electric field of reconstructed Bessel beams (polarization (p_x, p_y) = (1, 0), beam order are L = 0,1,2) on yoz plane in FDTD are displayed in Fig. 1. These beams are assumed propagating in the + z direction and half cone angles are θ₀ = 30°. The incident field is reconstructed by the superposition of 60 plane waves in FDTD. Each plane wave is calculated on the TF/SF boundary and propagates inside the total field region. As shown in Fig. 1, the first three columns display time domain results at three time instants, including 500Δt, 700Δt, and 900Δt. Using Fourier transform, results in the frequency domain at λ = 600nm can be obtained and are shown in the rightmost column. There is a bright region in the center of the zero-order Bessel beam while dark region can be found in the center of the first-order and the second-order beams.

Fig. 1 Amplitude distributions of the x-component electric field of reconstructed Bessel beams with θ₀ = 30° on yoz plane in FDTD simulations. (a) L = 0, (b) L = 1, (c) L = 2. The beams are traveling along the + z direction. The first three columns show the time domain results at different time instants: 500Δt, 700Δt, 900Δt from left to right. The rightmost column shows the results in the frequency domain at λ = 600nm after Fourier transform.

Download Full Size | PDF

To verify the accuracy of the present method, comparisons of the reconstructed field in FDTD and the original field are implemented. The x-component electric fields amplitude along the y-axis are shown in Fig. 2(a), in which EFDTD x,L = 0 and Eana x,L = 0 denote the results reconstructed in FDTD and that calculated using the analytical expression of Eq. (4) in [15], respectively. Relative deviation |EFDTD x− Eana x|/max(|Eana x|) are displayed in Fig. 2 (b), where the maximum value is about 4 × 10⁻³. As shown in Fig. 2, very good agreements have been achieved between Bessel beams generated in FDTD and analytical results.

Fig. 2 Comparisons between Bessel beams with θ₀ = 30°, λ = 600nm in FDTD and analytical results along the y-axis. (a) Amplitude of the x-component electric field. (b) Relative deviation |EFDTD x− Eana x|/max(|Eana x|).

Download Full Size | PDF

The determination of the number of plane wave N used in the reconstruction procedure is a key problem which should be carefully addressed. To analyze the influence of plane wave number N on the reconstruction of Bessel beams in FDTD, Bessel beams with certain wavelength and half cone angle are reconstructed with different numbers of plane waves. The x-component electric field amplitude on xoy plane for reconstructed zero-order Bessel beams in FDTD are displayed in Fig. 3. The original fields are also displayed for the purpose of comparison. As shown in Fig. 3, the less number of plane waves we use, the smaller valid area the reconstructed beam occupies. To reconstruct a Bessel beam accurately in a certain range of area, more number of plane waves are required for larger half cone angle or for shorter wavelength beams. In a further investigation, as Fig. 4 shows, we found that the same conclusion holds for higher-order beams. Hence, N should be carefully chosen according to the computation area, half-cone angle, as well as the wavelength of the incident beam. The standard way to determine the number of plane wave N used in FDTD should be comparing the reconstructed field with the original field. In the following simulations of this paper, N = 60 is used.

Fig. 3 Comparisons of the x-component electric field amplitude on xoy plane for zero-order Bessel beams between FDTD and original field. Different numbers of plane waves are chosen N = 15,30,60 to represent Bessel beams. Parameters of Bessel beams are: (a) θ₀ = 30°, λ = 600nm, (b) θ₀ = 80°, λ = 600nm and (c) θ₀ = 80°, λ = 800nm.

Download Full Size | PDF

Fig. 4 Comparisons of the x-component electric field amplitude on xoy plane for first-order Bessel beams between FDTD and original field. Different numbers of plane waves are chosen N = 15,30,60 to represent Bessel beams. Parameters of Bessel beams are: (a) θ₀ = 30°, λ = 600nm, (b) θ₀ = 80°, λ = 600nm.

Download Full Size | PDF

To further illustrate the validity of our approach in dealing with the scattering problem, the scattering of a homogeneous dielectric particle illuminated by a zero-order Bessel beam is analyzed. Figure 5 shows the scattering intensity distribution in the far zone of a spherical particle in vacuum illuminated by a zero-order Bessel beam (half cone angle is θ₀ = 30°, wavelengths are λ = 532nm and λ = 633nm). Sphere radius is r = 2μm and reflective index is m = 2.0. Note that the scattering angle ϕ_sc = 0° and ϕ_sc = 90° correspond to the xoz-plane and yoz-plane, respectively. For the purpose of comparison, FDTD(2,2), FDTD(2,4) and GLMT are used to solve the same problem. As shown in Fig. 5, the results of FDTD(2,4) achieved higher accuracy than FDTD(2,2). Normally, a grid spacing of at least λ/20 is needed for FDTD calculation [27]. In this example, the incident wavelength is λ = 633nm, grid spacing is ∆x = 25nm≈λ/25 in vacuum. But it decreases to about 1/12.5 of wavelength inside the particle with reflective index m = 2.0. The average relative error of the results from FDTD(2,2) is about 15% while it is less than 4% from FDTD(2,4). It is the numerical dispersion and anisotropy under such large grid spacing that causes the significant error when using FDTD(2,2) while higher order method FDTD(2,4) can get more accurate results. It worth mentioning that the numerical error varies with size parameter, refractive index of the object and it can be reduced by decreasing the grid size.

Fig. 5 Comparisons of scattering intensity for a dielectric sphere (r = 2.0μm, m = 2) illuminated by a zero-order Bessel beam (half cone angle is θ₀ = 30°) obtained from FDTD(2,2), FDTD(2,4) and GLMT. (a) ϕ_sc = 0°, (b) ϕ_sc = 90° at λ = 532nm and (c) ϕ_sc = 0°, (d) ϕ_sc = 90° at λ = 633nm.

Download Full Size | PDF

Far-field results for higher-order Bessel beams are shown in Fig. 6. Scattering intensity of a spherical particle (r = 2μm, m = 1.33) illuminated by a first-order Bessel beam (θ₀ = 20°, λ = 628.3nm) is shown in Fig. 6(a). It is interesting to find that a very low scattering intensity can be observed in forward θ_sc = 0° and backward θ_sc = 180° direction. Figure 6(b) shows the scattering intensity of a prolate spheroid particle with m = 1.65 illuminated by a second-order Bessel beam (θ₀ = 15°, λ = 633nm). Geometric size of the spheroid is determined by semi-major axes a = λ, semi-minor axes b = 0.5λ, and major axes coincide with the z-axis. Surface integral equation method (SIEM) is used to compared with FDTD(2,4). Both results in Fig. 6 show good agreements between different methods.

Fig. 6 Scattering results for higher-order Bessel beams. (a) Spherical particle (r = 2μm, m = 1.33) under a first-order Bessel beam (θ₀ = 20°, λ = 628.3nm) illumination, results are obtained from FDTD(2,4) and GLMT. (b) Spheroid particle (semi-major axes a = λ, semi-minor axes b = 0.5λ) under a second-order Bessel beam (θ₀ = 15°, λ = 633nm) illumination, results are obtained from FDTD(2,4) and SIEM.

Download Full Size | PDF

To show the potential in solving the scattering problems that involving complex particles, the internal and near-surface field distributions for a two-layer hemisphere particle illuminated by Bessel beams (half cone angle is θ₀ = 20°, wavelength is λ = 633nm) are displayed in Fig. 7. Radius and refractive index are r_s = 2μm and m_s = 1.6 for shell region, r_c = 1μm and m_c = 1.33 for core region. As shown in Fig. 7, when beam order is L = 0, a narrow focus is formed at the shadow side surface of the particle. This phenomenon is named by photonic nanojet (PNJ) and has attracted lots of attention due to its potential applications in super-resolution imaging, lithography, and others [43]. Nevertheless, when L = 1 and 2, the PNJ no longer exit while a curved focus can be observed. With the help of modeling technique in FDTD, a wide variety of objects interact with Bessel beams can be calculated in our approach.

Fig. 7 Intensity distributions of electric field on yoz plane when a two-layer hemisphere is illuminated by Bessel beams (θ₀ = 20°, λ = 633nm) with different beam order: (a) L = 0, (b) L = 1, (c) L = 2. Radius and the refractive index of the particle are r_s = 2μm and m_s = 1.6 for shell region, r_c = 1μm and m_c = 1.33 for core region.

Download Full Size | PDF

It worth mentioning that more accurate results can be obtained by using a finer grid mesh in FDTD simulation. When the grid size reduces to 1/40 of wavelength, traditional FDTD(2,2) also shows good performance, but it costs huge computation resources.

5. Conclusions

Generation of structured beam sources in finite-difference time-domain (FDTD) method is a key and difficult task as one intends to apply FDTD to the analysis of the electromagnetic or light scattering by complex objects illuminated by structured beams. In this paper, an approach to generate Bessel beam sources in FDTD is presented. Angular spectrum representation (ASR) is used to describe Bessel beams as a superposition of plane waves, which are injected into FDTD by total-field/scattered-field (TF/ST) method. The validity and accuracy of the present method are carefully verified by comparing the reconstructed field in FDTD with the original field. Higher-order method FDTD(2,4) is applied to reduce error due to numerical dispersion and anisotropy when using large grid spacing in the simulation and the average relative error is less than 4%. The far-field scattering of dielectric spherical or spheroid particle illuminated by zero or higher order Bessel beams are calculated using FDTD, the results are compared with those calculated using GLMT and SIEM. Very good agreements have been achieved which partially indicate the correctness of our method. To show the capability of our method in solving the EM scattering problems of complex particles illuminated by Bessel beams, the internal and near-surface field distributions for a two-layer hemisphere particle illuminated by Bessel beams are also presented. It is worth mentioning that by changing the parameters in vector Q, Bessel beams with different polarizations can be generated directly, and this approach can be applied to the generation of other structured beams in FDTD by using its angular spectrum representation.

Funding

National Natural Science Foundation of China (61431010, 61501350); National Key Basic Research Program of China (2014CB340203); Natural Science Foundation of Shaanxi Province (2018JM6016).

References

1. J. Durnin, “Exact solutions for nondiffracting beams. I. The scalar theory,” J. Opt. Soc. Am. A 4(4), 651–654 (1987). [CrossRef]

2. D. McGloin and K. Dholakia, “Bessel beams: diffraction in a new light,” Contemp. Phys. 46(1), 15–28 (2005). [CrossRef]

3. Z. Bouchal, J. Wagner, and M. Chlup, “Self-reconstruction of a distorted nondiffracting beam,” Opt. Commun. 151(4–6), 207–211 (1998). [CrossRef]

4. N. Mphuthi, R. Botha, and A. Forbes, “Are Bessel beams resilient to aberrations and turbulence?” J. Opt. Soc. Am. A 35(6), 1021–1027 (2018). [CrossRef] [PubMed]

5. F. Courvoisier, R. Stoian, and A. Couairon, "[INVITED] Ultrafast laser micro- and nano-processing with nondiffracting and curved beams,” Opt. Laser Technol. 80, 125–137 (2016). [CrossRef]

6. V. Garcés-Chávez, D. McGloin, H. Melville, W. Sibbett, and K. Dholakia, “Simultaneous micromanipulation in multiple planes using a self-reconstructing light beam,” Nature 419(6903), 145–147 (2002). [CrossRef] [PubMed]

7. D. Li and K. Imasaki, “Vacuum laser-driven acceleration by a slits-truncated Bessel beam,” Appl. Phys. Lett. 86(3), 031110 (2005). [CrossRef]

8. D. Fan, L. Wang, and Y. Ekinci, “Nanolithography using Bessel Beams of Extreme Ultraviolet Wavelength,” Sci. Rep. 6(1), 31301 (2016). [CrossRef] [PubMed]

9. W. Yu, Z. Ji, D. Dong, X. Yang, Y. Xiao, Q. Gong, P. Xi, and K. Shi, “Super-resolution deep imaging with hollow Bessel beam STED microscopy,” Laser Photonics Rev. 10(1), 147–152 (2016). [CrossRef]

10. M. Ettorre, S. M. Rudolph, and A. Grbic, “Generation of Propagating Bessel Beams Using Leaky-Wave Modes: Experimental Validation,” IEEE Trans. Antenn. Propag. 60(6), 2645–2653 (2012). [CrossRef]

11. X. Wei, C. Liu, L. Niu, Z. Zhang, K. Wang, Z. Yang, and J. Liu, “Generation of arbitrary order Bessel beams via 3D printed axicons at the terahertz frequency range,” Appl. Opt. 54(36), 10641–10649 (2015). [CrossRef] [PubMed]

12. J. Arlt and K. Dholakia, “Generation of high-order Bessel beams by use of an axicon,” Opt. Commun. 177(1-6), 297–301 (2000). [CrossRef]

13. G. Gouesbet and G. Gréhan, Generalized lorenz-mie theories (Springer, 2011).

14. F. G. Mitri, “Arbitrary scattering of an electromagnetic zero-order Bessel beam by a dielectric sphere,” Opt. Lett. 36(5), 766–768 (2011). [CrossRef] [PubMed]

15. F. G. Mitri, “Electromagnetic Wave Scattering of a High-Order Bessel Vortex Beam by a Dielectric Sphere,” IEEE Trans. Antenn. Propag. 59(11), 4375–4379 (2011). [CrossRef]

16. F. G. Mitri, R. X. Li, L. X. Guo, and C. Y. Ding, “Resonance scattering of a dielectric sphere illuminated by electromagnetic Bessel non-diffracting (vortex) beams with arbitrary incidence and selective polarizations,” Ann. Phys. 361, 120–147 (2015). [CrossRef]

17. L. Han, “Scattering of a high-order Bessel beam by a spheroidal particle,” J. Quant. Spectrosc. Radiat. Transf. 211, 129–137 (2018). [CrossRef]

18. R. Li, K. F. Ren, X. e. Han, Z. Wu, L. Guo, and S. Gong, “Analysis of radiation pressure force exerted on a biological cell induced by high-order Bessel beams using Debye series,” J. Quant. Spectrosc. Radiat. Transf. 126, 69–77 (2013). [CrossRef]

19. J. J. Wang, Y. P. Han, J. Y. Chang, and Z. Y. Chen, “Light scattering of a Bessel beam by a nucleated biological cell: An eccentric sphere model,” J. Quant. Spectrosc. Radiat. Transf. 206, 22–30 (2018). [CrossRef]

20. F. G. Mitri, R. X. Li, L. X. Guo, and C. Y. Ding, “Optical tractor Bessel polarized beams,” J. Quant. Spectrosc. Radiat. Transf. 187, 97–115 (2017). [CrossRef]

21. R. Li, C. Ding, and F. G. Mitri, “Optical spin torque induced by vector Bessel (vortex) beams with selective polarizations on a light-absorptive sphere of arbitrary size,” J. Quant. Spectrosc. Radiat. Transf. 196, 53–68 (2017). [CrossRef]

22. F. G. Mitri, “Reverse orbiting and spinning of a Rayleigh dielectric spheroid in a J0 Bessel optical beam,” J. Opt. Soc. Am. B 34(10), 2169–2178 (2017). [CrossRef]

23. M. Yang, Y. Wu, X. Sheng, and K. F. Ren, “Computational study of scattering of a zero-order Bessel beam by large nonspherical homogeneous particles with the multilevel fast multipole algorithm,” J. Opt. 19(12), 125606 (2017). [CrossRef]

24. Z. Cui, Y. Han, Z. Chen, and L. Han, “Scattering of Bessel beam by arbitrarily shaped composite particles with core–shell structure,” J. Quant. Spectrosc. Radiat. Transf. 144, 108–116 (2014). [CrossRef]

25. Z. Gong, W. Li, Y. Chai, Y. Zhao, and F. G. Mitri, “T-matrix method for acoustical Bessel beam scattering from a rigid finite cylinder with spheroidal endcaps,” Ocean Eng. 129, 507–519 (2017). [CrossRef]

26. Z. Gong, W. Li, F. G. Mitri, Y. Chai, and Y. Zhao, “Arbitrary scattering of an acoustical Bessel beam by a rigid spheroid with large aspect-ratio,” J. Sound Vibrat. 383, 233–247 (2016). [CrossRef]

27. A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method, 3rd ed. (Artech house, 2005).

28. İ. R. Çapoğlu, A. Taflove, and V. Backman, “Generation of an incident focused light pulse in FDTD,” Opt. Express 16(23), 19208–19220 (2008). [CrossRef] [PubMed]

29. İ. R. Çapoğlu, A. Taflove, and V. Backman, “Computation of tightly-focused laser beams in the FDTD method,” Opt. Express 21(1), 87–101 (2013). [CrossRef] [PubMed]

30. P. Török, P. R. T. Munro, and E. E. Kriezis, “High numerical aperture vectorial imaging in coherent optical microscopes,” Opt. Express 16(2), 507–523 (2008). [CrossRef] [PubMed]

31. W. Sun, Y. Hu, C. Weimer, K. Ayers, R. R. Baize, and T. Lee, “A FDTD solution of scattering of laser beam with orbital angular momentum by dielectric particles: Far-field characteristics,” J. Quant. Spectrosc. Radiat. Transf. 188, 200–213 (2017). [CrossRef]

32. S.-Y. Sung and Y.-G. Lee, “Trapping of a micro-bubble by non-paraxial Gaussian beam: computation using the FDTD method,” Opt. Express 16(5), 3463–3473 (2008). [CrossRef] [PubMed]

33. A. Chen, J. Wang, Y. Han, Z. Cui, and M. Yu, “Implementation of nondiffracting Bessel beam sources in FDTD for scattering by complex particles,” Opt. Express 26(20), 26766–26775 (2018). [CrossRef]

34. J. Xi, Q. Li, and J. Wang, “Numerical simulation of Bessel beams by FDTD employing the superposition principle,” Optik (Stuttg.) 118(7), 315–318 (2007). [CrossRef]

35. A. Elmaklizi, D. Reitzle, A. Brandes, and A. Kienle, “Penetration depth of focused beams in highly scattering media investigated with a numerical solution of Maxwell’s equations in two dimensions,” J. Biomed. Opt. 20(6), 065007 (2015). [CrossRef] [PubMed]

36. S. R. Mishra, “A vector wave analysis of a Bessel beam,” Opt. Commun. 85(2-3), 159–161 (1991). [CrossRef]

37. J. J. Wang, T. Wriedt, J. A. Lock, and L. Mädler, “General description of circularly symmetric Bessel beams of arbitrary order,” J. Quant. Spectrosc. Radiat. Transf. 184, 218–232 (2016). [CrossRef]

38. areT. Čižmár, V. Kollárová, Z. Bouchal, and P. Zemánek, “Sub-micron particle organization by self-imaging of non-diffracting beams,” New J. Phys. 8(3), 43 (2006). [CrossRef]

39. J. A. Lock, “Angular spectrum and localized model of Davis-type beam,” J. Opt. Soc. Am. A 30(3), 489–500 (2013). [CrossRef] [PubMed]

40. L. Novotny and B. Hecht, Principles of nano-optics (Cambridge university, 2006).

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

42. J. Stoer and R. Bulirsch, Introduction to numerical analysis (Springer Science & Business Media, 2013).

43. B. S. Luk’yanchuk, R. Paniagua-Domínguez, I. Minin, O. Minin, and Z. Wang, “Refractive index less than two: photonic nanojets yesterday, today and tomorrow [Invited],” Opt. Mater. Express 7(6), 1820–1847 (2017). [CrossRef]

Generation of Bessel beam sources in FDTD

Abstract

1. Introduction

2. Angular spectrum representation of Bessel beams

3. FDTD implementation of Bessel beam sources

4. Simulation results and discussions

5. Conclusions

Funding

References

Cited By

Figures (7)

Equations (9)

Optics Express