- Email: [email protected]
Numerical simulation of Bessel beams by FDTD employing the superposition principle
ARTICLE IN PRESS
Optik
Optics
Optik 118 (2007) 315–318 www.elsevier.de/ijleo
Numerical simulation of Bessel beams by FDTD employing the superposition principle Jiefeng Xi, Qing Li, Jia Wang State Key Laboratory of Precision Measurement Technology and Instruments, Department of Precision Instruments, Tsinghua University, Beijing 100084, PR China Received 10 October 2005; received in revised form 8 January 2006; accepted 23 March 2006
Abstract Finite difference time domain (FDTD) method is adopted to build a Bessel beams simulation model according to homogeneousness and linearity of the Maxwell equations in source-free region. Validation for this model is confirmed by comparing the simulation results with the theoretical results solved with a vector Helmholtz equation in free space and good agreement with maximum error 2% has been demonstrated. It is indicated that FDTD could be an effective approach to analyze other complicated models of Bessel beams in source-free region by means of superposition principle. r 2006 Elsevier GmbH. All rights reserved. Keywords: Bessel beams; Finite difference time domain; Near-field optics
1. Introduction Bessel beams, as one of cylindrically symmetric waves, is a diffractionless electromagnetic wave, the scalar theory of which was proposed for the first time by Durnin [1] and Durnin et al. [2]. It plays an important role in many optical aspects [3–6] as it provides convenient techniques for avoiding or reducing an inevitable diffractive spatial spreading. Vector solution of Bessel beams was further calculated by solving a vector potential and vector Helmlotz equation [7,8]. Since the nondiffracting property of Bessel beams, several research works focus on its application on near-field optics [9–11]. Since, however, most numerical algorithms in far-field optics cannot be applied in near-field optics, finite Corresponding author.
E-mail address: [email protected] (J. Wang). 0030-4026/$ - see front matter r 2006 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2006.03.026
difference time domain (FDTD) method is one of the most popular algorithms adapted in near-field optics as it is well-known for solving the electromagnetic problems with arbitrary boundary conditions and inhomogeneous materials [12]. The continuous space is discretized into small cubes called Yee cells with the size less than size of relevant features and the time is discretized into small steps much less than the period of the interested electromagnetic wave so that two curl equations of Maxwell equations become difference equations. Generally, thep result ffiffiffi of this algorithm would be stable if DtpDu= 3v according to numerical stability condition in 3D FDTD [12], where Du is the side length of a cubic cell, Dt is the time step, v is the optical velocity in the considered medium. Mur [13] and Liao et al. [7] absorbing boundary conditions are used practically in the algorithm in order to make actual computation possible. In this paper, exact vector solution of zero-order transversal magnetic (TM) mode Bessel beams has been
ARTICLE IN PRESS J. Xi et al. / Optik 118 (2007) 315–318
concluded by solving vector Helmlotz equation. And then, a zero-order TM mode Bessel beams simulation model is built by means of superposition of plane wave with cylindrical linear polarization. At last, comparison of the simulation result and the theoretical solution is carried out as to the transversal and longitudinal intensity distribution at any cross-section.
Transversal Intensity Longitudinal Intensity
2
Intensity (a.u)
316
1.5
1
0.5
2. Exact solution of Bessel beams The vector monochromatic wave equation can be represented as the following vector Helmholz equation: r2 C þ k2 C ¼ 0,
(1)
where k is wave vector of a monochromatic wave and C could be any vector, such as electrical field, magnetic field and vector potential. According to electromagnetic theory and Bouchal and Olivik’s work [5], a vector potential could be defined as the follows: iX A¼ ðan Mn þ bn Nn þ cn Ln Þ; (2) o n where Mn, Nn and Ln are derived from the solution of scalar Helmholz equation cn [8]. So magnetic induction and electrical vector could be written in this form: iX ð bn M n þ an N n Þ B¼rA¼ o n X E¼ ðan Mn þ bn Nn Þ. ð3Þ n
A, B, E will satisfy vector Helmholz equation, respectively, if cn is solutions of scalar Helmhotz equation so that a generalized vector solution of Bessel beams could be solved easily [5]. Then we could derive zero-order TM mode Bessel beams (shorten as Bessel beams later), or namely zero-order radially polarized wave, in the free space from the generalized solution (as shown in Fig. 1): E r ¼ i2aJ 1 ðarÞ expðibzÞ, E j ¼ 0, 2
E z ¼ a J 1 ðarÞ expðibzÞ=b,
ð4Þ
where r, z are two cylindrical coordinates, a þ b ¼ k2 and k is wave vector of a monochromatic wave. Furthermore, magnetic induction of Bessel beams can be written as 2
2
Br ¼ Bz ¼ 0, Bj ¼ i2J 1 ðarÞ expðibzÞ=v,
ð5Þ
where v is the optical velocity in the considered medium. Hence, it is clearly shown that the transversal distribution of Bessel beams keeps unchanged along the z direction, which is nondiffracting property of Bessel beams in free space. And average energy flow Saverage of
0 -500
-400
-300
-200
-100
0
100
200
300
400
500
Radius (nm)
Fig. 1. Theoretical result of intensity distribution of Bessel beams in transversal and longitudinal component ða ¼ k sin 601Þ.
Bessel beams is only along the z direction, which indicates that energy propagating direction of Bessel beams is z direction and is perpendicular to its equiphase surface.
3. Numerical model of Bessel beams in FDTD Using a commercial FDTD package (Remcom XFDTD 6.1), either plane wave or electric dipole can only be set as an excitation source. Hence, this commercial package cannot be used directly to simulate Bessel beams. However, we can build a numerical model to achieve this goal indirectly in XFDTD based on the theoretical analysis in the last section. According to the Poisson integral formula of Bessel function, Eq. (4) can be written as Z a expðibzÞ 2p iar cos y Er ¼ e cos y dy, 2p 0 Z a2 expðibzÞ 2p iar cos y Ez ¼ e dy. ð6Þ b2p 0 By defining a ¼ k sin j and b ¼ k cos j, Eq. (5) can be expressed as plane waves with cylindrical vector symmetry are interfered to generate a Bessel beams (as shown in Fig. 2). So incident plane waves with cylindrical vector symmetry are set with an identical source-free FDTD model, which is used to calculate time and again. Finally, we develop a Matlab program to superpose the calculated results of complex amplitude to form the complex amplitude of Bessel beams after compensating phase difference among plane waves as post-processing of XFDTD calculation. Validation for this model based on XFDTD program could be confirmed by comparing the simulation result with theoretical one in the free space as analyzed previously. The incident optical waves are linear
ARTICLE IN PRESS J. Xi et al. / Optik 118 (2007) 315–318
317
in Fig. 3(d)). Further study shows, as same as the theoretical predication, that the equiphase surface of simulation result is a plane and the average transversal energy flow of simulation result is almost zero with instantaneous fluctuation.
4. Conclusions and discussions
Fig. 2. Scheme of Bessel beams superposed by linear polarized plane wave.
(a)
At last, hence, it is clearly proven that the model based on a commercial FDTD package is a sound and effective approach to simulate Bessel beams in free space. The error between simulation result and theoretical one are mainly caused by two reasons, the first of which is elevation angle y cannot span continuously and the second of which is numerical dispersion. We estimate the error caused by the two reasons above (as shown in Fig. 4) and the result shows that the total error is dominated by the two reasons. Generally, this model is effective for any source-free simulation region due to the linearity and superposition of Maxwell equations in source-free region. So it could
(b) 0.015
(c)
(d)
Fig. 3. Simulation result on XY cross-section of Bessel beams intensity of: (a) transversal component, (b) longitudinal component, (c) total component ða ¼ k sin 601Þ and (d) the relative error between simulation result and theoretical one.
Estimated Relative Error
0.01 0.005 0 -0.005 -0.01 -0.015 Longitudinal Error Transversal Error
-0.02
-0.025 -500 -400 -300 -200 -100
0
100
200
300
400
500
Radium (nm)
Fig. 4. Estimated error caused by discontinuousness and numerical dispersion.
500 450 400 Total Intensity
FWHM (nm)
polarized with cylindrical vector symmetry, the polarization parameters of which are j ¼ 301 and y spans 0–3601 every 51 (j is the azimuthal angle referenced with x direction and y is the elevation angle referenced with z direction). And the amplitude of the incident wave is 1 V/m and the wavelength is l ¼ 632:8 nm. The model, where the whole space is 1 1 1 mm formed by unit cubic cell with 1000 nm3, is a free space without any dielectric or metal. The time step is set to Dt ¼ 19.2581 as and the program was run for 800 time steps, sufficient for reaching steady state. The FDTD program is run in the existence of a second-order Liao outer radiation boundary condition, which is numerical absorber designed to allow electromagnetic fields to be radiated or scattered by the FDTD geometry to be absorbed with very little reflection from the boundary. The simulation result is shown in Fig. 3(a)–(c), which indicates a very good agreement between them as well as an average relative error 2% in longitudinal component and 0.5% in transversal component (as shown
Transversal Intensity
350
Longitude intensity 300 250 200 150
0
100
200
300
400
500
Distance between aperture and observation plane (nm)
Fig. 5. Simulation result of FWHM of Bessel beams passing a sub-wavelength aperture (diameter of the aperture d ¼ 300 nm, a ¼ k sin 601).
ARTICLE IN PRESS 318
J. Xi et al. / Optik 118 (2007) 315–318
be general to analyze the near-field property of Bessel beams, such as evanescent Bessel beams and the response of Bessel beams to a sub-wavelength aperture. The relation between full-width at half-maximum (FWHM) and the distance, for example, are shown in Fig. 5. Even other beams with different transversal distribution could be calculated under the superposition principle if the exact vector solutions of these beams are known.
References [1] J. Durnin, Exact solution for nondiffracting beams. I. The scalar theory, J. Opt. Soc. Am. A 4 (1987) 651–654. [2] J. Durnin, J.J. Miceli Jr., J.H. Eberly, Diffraction-free beams, Phys. Rev. Lett. 58 (15) (1987) 1499–1501. [3] Q. Zhan, Trapping metallic Rayleigh particles with radial polarization, Opt. Exp. 12 (15) (2004) 3377–3383. [4] S.R. Mishra, A vector wave analysis of a Bessel beam, Opt. Commun. 85 (1991) 159–161. [5] Z. Bouchal, M. Olivik, Non-diffractive vector Bessel beams, J. Mod. Opt. 42 (8) (1995) 1551–1566.
[6] S. Rushin, A. Leizer, Evanescent Bessel beams, J. Opt. Soc. Am. A 15 (5) (1998) 1139–1143. [7] Z.P. Liao, H.L. Wong, G.P. Yang, Y.F. Yuan, A transmitting boundary for transient wave analysis, Sci. Sinica 28 (1984) 1063C1076. [8] J.A. Stratton, Electromagnetic Theory, McGraw-Hill Press, New York, 1941 354pp. [9] J. Arlt, K. Dholakia, J. Soneson, E.M. Wright, Optical dipole traps and atomic waveguides based Bessel light beams, Phys. Rev. A 63 (2001) 063602. [10] D. McGloin, V. Garces-Chavez, K. Dholakia, Interfering Bessel beams for optical micromanipulation, Opt. Lett. 28 (8) (2003) 657–659. [11] T. Grosjean, D. Courjon, Van Labeke, Bessel beams as virtual tips for near-field optics, J. Microsc. (Oxford) 210 (3) (2003) 319–323. [12] K.S. Yee, Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media, IEEE Trans. Antennas Propag. 14 (1966) 302C307. [13] G. Mur, Absorbing boundary conditions for the finite difference approximation of the time-domain electromagnetic-field equations, IEEE Trans. EMC 23 (1981) 1073C1077.
Optik
Optics
Optik 118 (2007) 315–318 www.elsevier.de/ijleo
Numerical simulation of Bessel beams by FDTD employing the superposition principle Jiefeng Xi, Qing Li, Jia Wang State Key Laboratory of Precision Measurement Technology and Instruments, Department of Precision Instruments, Tsinghua University, Beijing 100084, PR China Received 10 October 2005; received in revised form 8 January 2006; accepted 23 March 2006
Abstract Finite difference time domain (FDTD) method is adopted to build a Bessel beams simulation model according to homogeneousness and linearity of the Maxwell equations in source-free region. Validation for this model is confirmed by comparing the simulation results with the theoretical results solved with a vector Helmholtz equation in free space and good agreement with maximum error 2% has been demonstrated. It is indicated that FDTD could be an effective approach to analyze other complicated models of Bessel beams in source-free region by means of superposition principle. r 2006 Elsevier GmbH. All rights reserved. Keywords: Bessel beams; Finite difference time domain; Near-field optics
1. Introduction Bessel beams, as one of cylindrically symmetric waves, is a diffractionless electromagnetic wave, the scalar theory of which was proposed for the first time by Durnin [1] and Durnin et al. [2]. It plays an important role in many optical aspects [3–6] as it provides convenient techniques for avoiding or reducing an inevitable diffractive spatial spreading. Vector solution of Bessel beams was further calculated by solving a vector potential and vector Helmlotz equation [7,8]. Since the nondiffracting property of Bessel beams, several research works focus on its application on near-field optics [9–11]. Since, however, most numerical algorithms in far-field optics cannot be applied in near-field optics, finite Corresponding author.
E-mail address: [email protected] (J. Wang). 0030-4026/$ - see front matter r 2006 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2006.03.026
difference time domain (FDTD) method is one of the most popular algorithms adapted in near-field optics as it is well-known for solving the electromagnetic problems with arbitrary boundary conditions and inhomogeneous materials [12]. The continuous space is discretized into small cubes called Yee cells with the size less than size of relevant features and the time is discretized into small steps much less than the period of the interested electromagnetic wave so that two curl equations of Maxwell equations become difference equations. Generally, thep result ffiffiffi of this algorithm would be stable if DtpDu= 3v according to numerical stability condition in 3D FDTD [12], where Du is the side length of a cubic cell, Dt is the time step, v is the optical velocity in the considered medium. Mur [13] and Liao et al. [7] absorbing boundary conditions are used practically in the algorithm in order to make actual computation possible. In this paper, exact vector solution of zero-order transversal magnetic (TM) mode Bessel beams has been
ARTICLE IN PRESS J. Xi et al. / Optik 118 (2007) 315–318
concluded by solving vector Helmlotz equation. And then, a zero-order TM mode Bessel beams simulation model is built by means of superposition of plane wave with cylindrical linear polarization. At last, comparison of the simulation result and the theoretical solution is carried out as to the transversal and longitudinal intensity distribution at any cross-section.
Transversal Intensity Longitudinal Intensity
2
Intensity (a.u)
316
1.5
1
0.5
2. Exact solution of Bessel beams The vector monochromatic wave equation can be represented as the following vector Helmholz equation: r2 C þ k2 C ¼ 0,
(1)
where k is wave vector of a monochromatic wave and C could be any vector, such as electrical field, magnetic field and vector potential. According to electromagnetic theory and Bouchal and Olivik’s work [5], a vector potential could be defined as the follows: iX A¼ ðan Mn þ bn Nn þ cn Ln Þ; (2) o n where Mn, Nn and Ln are derived from the solution of scalar Helmholz equation cn [8]. So magnetic induction and electrical vector could be written in this form: iX ð bn M n þ an N n Þ B¼rA¼ o n X E¼ ðan Mn þ bn Nn Þ. ð3Þ n
A, B, E will satisfy vector Helmholz equation, respectively, if cn is solutions of scalar Helmhotz equation so that a generalized vector solution of Bessel beams could be solved easily [5]. Then we could derive zero-order TM mode Bessel beams (shorten as Bessel beams later), or namely zero-order radially polarized wave, in the free space from the generalized solution (as shown in Fig. 1): E r ¼ i2aJ 1 ðarÞ expðibzÞ, E j ¼ 0, 2
E z ¼ a J 1 ðarÞ expðibzÞ=b,
ð4Þ
where r, z are two cylindrical coordinates, a þ b ¼ k2 and k is wave vector of a monochromatic wave. Furthermore, magnetic induction of Bessel beams can be written as 2
2
Br ¼ Bz ¼ 0, Bj ¼ i2J 1 ðarÞ expðibzÞ=v,
ð5Þ
where v is the optical velocity in the considered medium. Hence, it is clearly shown that the transversal distribution of Bessel beams keeps unchanged along the z direction, which is nondiffracting property of Bessel beams in free space. And average energy flow Saverage of
0 -500
-400
-300
-200
-100
0
100
200
300
400
500
Radius (nm)
Fig. 1. Theoretical result of intensity distribution of Bessel beams in transversal and longitudinal component ða ¼ k sin 601Þ.
Bessel beams is only along the z direction, which indicates that energy propagating direction of Bessel beams is z direction and is perpendicular to its equiphase surface.
3. Numerical model of Bessel beams in FDTD Using a commercial FDTD package (Remcom XFDTD 6.1), either plane wave or electric dipole can only be set as an excitation source. Hence, this commercial package cannot be used directly to simulate Bessel beams. However, we can build a numerical model to achieve this goal indirectly in XFDTD based on the theoretical analysis in the last section. According to the Poisson integral formula of Bessel function, Eq. (4) can be written as Z a expðibzÞ 2p iar cos y Er ¼ e cos y dy, 2p 0 Z a2 expðibzÞ 2p iar cos y Ez ¼ e dy. ð6Þ b2p 0 By defining a ¼ k sin j and b ¼ k cos j, Eq. (5) can be expressed as plane waves with cylindrical vector symmetry are interfered to generate a Bessel beams (as shown in Fig. 2). So incident plane waves with cylindrical vector symmetry are set with an identical source-free FDTD model, which is used to calculate time and again. Finally, we develop a Matlab program to superpose the calculated results of complex amplitude to form the complex amplitude of Bessel beams after compensating phase difference among plane waves as post-processing of XFDTD calculation. Validation for this model based on XFDTD program could be confirmed by comparing the simulation result with theoretical one in the free space as analyzed previously. The incident optical waves are linear
ARTICLE IN PRESS J. Xi et al. / Optik 118 (2007) 315–318
317
in Fig. 3(d)). Further study shows, as same as the theoretical predication, that the equiphase surface of simulation result is a plane and the average transversal energy flow of simulation result is almost zero with instantaneous fluctuation.
4. Conclusions and discussions
Fig. 2. Scheme of Bessel beams superposed by linear polarized plane wave.
(a)
At last, hence, it is clearly proven that the model based on a commercial FDTD package is a sound and effective approach to simulate Bessel beams in free space. The error between simulation result and theoretical one are mainly caused by two reasons, the first of which is elevation angle y cannot span continuously and the second of which is numerical dispersion. We estimate the error caused by the two reasons above (as shown in Fig. 4) and the result shows that the total error is dominated by the two reasons. Generally, this model is effective for any source-free simulation region due to the linearity and superposition of Maxwell equations in source-free region. So it could
(b) 0.015
(c)
(d)
Fig. 3. Simulation result on XY cross-section of Bessel beams intensity of: (a) transversal component, (b) longitudinal component, (c) total component ða ¼ k sin 601Þ and (d) the relative error between simulation result and theoretical one.
Estimated Relative Error
0.01 0.005 0 -0.005 -0.01 -0.015 Longitudinal Error Transversal Error
-0.02
-0.025 -500 -400 -300 -200 -100
0
100
200
300
400
500
Radium (nm)
Fig. 4. Estimated error caused by discontinuousness and numerical dispersion.
500 450 400 Total Intensity
FWHM (nm)
polarized with cylindrical vector symmetry, the polarization parameters of which are j ¼ 301 and y spans 0–3601 every 51 (j is the azimuthal angle referenced with x direction and y is the elevation angle referenced with z direction). And the amplitude of the incident wave is 1 V/m and the wavelength is l ¼ 632:8 nm. The model, where the whole space is 1 1 1 mm formed by unit cubic cell with 1000 nm3, is a free space without any dielectric or metal. The time step is set to Dt ¼ 19.2581 as and the program was run for 800 time steps, sufficient for reaching steady state. The FDTD program is run in the existence of a second-order Liao outer radiation boundary condition, which is numerical absorber designed to allow electromagnetic fields to be radiated or scattered by the FDTD geometry to be absorbed with very little reflection from the boundary. The simulation result is shown in Fig. 3(a)–(c), which indicates a very good agreement between them as well as an average relative error 2% in longitudinal component and 0.5% in transversal component (as shown
Transversal Intensity
350
Longitude intensity 300 250 200 150
0
100
200
300
400
500
Distance between aperture and observation plane (nm)
Fig. 5. Simulation result of FWHM of Bessel beams passing a sub-wavelength aperture (diameter of the aperture d ¼ 300 nm, a ¼ k sin 601).
ARTICLE IN PRESS 318
J. Xi et al. / Optik 118 (2007) 315–318
be general to analyze the near-field property of Bessel beams, such as evanescent Bessel beams and the response of Bessel beams to a sub-wavelength aperture. The relation between full-width at half-maximum (FWHM) and the distance, for example, are shown in Fig. 5. Even other beams with different transversal distribution could be calculated under the superposition principle if the exact vector solutions of these beams are known.
References [1] J. Durnin, Exact solution for nondiffracting beams. I. The scalar theory, J. Opt. Soc. Am. A 4 (1987) 651–654. [2] J. Durnin, J.J. Miceli Jr., J.H. Eberly, Diffraction-free beams, Phys. Rev. Lett. 58 (15) (1987) 1499–1501. [3] Q. Zhan, Trapping metallic Rayleigh particles with radial polarization, Opt. Exp. 12 (15) (2004) 3377–3383. [4] S.R. Mishra, A vector wave analysis of a Bessel beam, Opt. Commun. 85 (1991) 159–161. [5] Z. Bouchal, M. Olivik, Non-diffractive vector Bessel beams, J. Mod. Opt. 42 (8) (1995) 1551–1566.
[6] S. Rushin, A. Leizer, Evanescent Bessel beams, J. Opt. Soc. Am. A 15 (5) (1998) 1139–1143. [7] Z.P. Liao, H.L. Wong, G.P. Yang, Y.F. Yuan, A transmitting boundary for transient wave analysis, Sci. Sinica 28 (1984) 1063C1076. [8] J.A. Stratton, Electromagnetic Theory, McGraw-Hill Press, New York, 1941 354pp. [9] J. Arlt, K. Dholakia, J. Soneson, E.M. Wright, Optical dipole traps and atomic waveguides based Bessel light beams, Phys. Rev. A 63 (2001) 063602. [10] D. McGloin, V. Garces-Chavez, K. Dholakia, Interfering Bessel beams for optical micromanipulation, Opt. Lett. 28 (8) (2003) 657–659. [11] T. Grosjean, D. Courjon, Van Labeke, Bessel beams as virtual tips for near-field optics, J. Microsc. (Oxford) 210 (3) (2003) 319–323. [12] K.S. Yee, Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media, IEEE Trans. Antennas Propag. 14 (1966) 302C307. [13] G. Mur, Absorbing boundary conditions for the finite difference approximation of the time-domain electromagnetic-field equations, IEEE Trans. EMC 23 (1981) 1073C1077.