Generation of a Bessel beam in FDTD using a cylindrical antenna
Kazem Ardaneh, Remo Giust, Benoit Morel, Francois Courvoisier
GGeneration of a Bessel beam in FDTD using a cylindrical antenna
Kazem Ardaneh(*), Remo Giust, Benoit Morel and Francois Courvoisier( † )FEMTO-ST institute, Univ. Bourgogne Franche-Comt´e, CNRS,15B avenue des Montboucons, 25030, Besanc¸on Cedex, France(*)[email protected]( † )[email protected] is a post-peer-review, pre-copyedit version of an article publishedin Optics Express. The final authenticated version is available online at: https://doi.org/10.1364/OE.385413 Bessel beams are becoming a very useful tool in many areas of optics and photonics, because of the invarianceof their intensity profile over an extended propagation range. Finite-Difference-Time-Domain (FDTD) approachis widely used for the modeling of the beam interaction with nanostructures. However, the generation of theBessel beam in this approach is a computationally challenging problem. In this work, we report an approach forthe generation of the infinite Bessel beams in three-dimensional FDTD. It is based on the injection of the Besselsolutions of Maxwell’s equations from a cylindrical hollow annulus. This configuration is compatible withParticle In Cell simulations of laser plasma interactions. This configuration allows using a smaller computationbox and is therefore computationally more efficient than the creation of a Bessel-Gauss beam from a wall andmodels more precisely the analytical infinite Bessel beam. Zeroth and higher-order Bessel beams with differentcone angles are successfully produced. We investigate the effects of the injector parameters on the error withrespect to the analytical solution. In all cases, the relative deviation is in the range of 0.01-7.0 percent. a r X i v : . [ phy s i c s . op ti c s ] F e b INTRODUCTION
A Bessel beam refers to a family of solutions for the wave equation in which the amplitude of field is expressed by the Besselfunction of the first kind [1]. It has attracted great interest in various branches in optics because its intensity profile is invariantas it propagates. Bessel beams are mainly utilized as optical traps [2–4], for optical manipulation [5], optical acceleration [6–8],light-sheet microscopy [9, 10], nonlinear optics, ultrashort pulse filamentation [11–13], and laser-material processing [14–18].Our main interest here is that Bessel beams can be used to generate long plasma rods with quasi-uniform density [19, 20].Particle-In-Cell (PIC) simulation [21, 22] is a conventional pathway to self-consistently model the interaction between laserbeams with plasmas. It uses the Finite-Difference-Time-Domain (FDTD) algorithm to advance the electric and magnetic fields.A specific challenge arises for the simulation of Bessel beams in FDTD because of their high aspect ratio. Bessel beams arean interference field in which the propagation length Z B is defined by the transverse extent W and the cone angle θ , i.e. theangle made by the interfering waves with the optical axis: Z B ∼ W/ tan θ [23, 24]. Hence, long propagation distances require awide transverse extent. However, the central spot radius ( r = 0 . λ/ sin θ for a zeroth-order Bessel beam) is generally muchsmaller than the beam transverse extent. It is therefore computationally extremely demanding to investigate the laser-plasmainteraction within the central lobe with high spatial resolution. This problem is illustrated in Fig. 1, where we show the time-integrated fluence distribution for a finite-energy zeroth-order Bessel-Gauss femtosecond pulse over a distance of only ≈ µ m ,injected from the left wall of the simulation box. The detail of this simulation is provided in Section .Our objective is to replace the simulation box with a smaller one, as shown by transparent white color box in Fig. 1. Thisis possible because in a Bessel beam, energy flows from the sides with a conical structure and because of the longitudinalinvariance of both the diffraction-free Bessel beams and the uniform plasma distribution. Therefore, we will use a cylindricalannulus injecting the electromagnetic fields and use periodic boundary condition for the surfaces parallel to the optical axis toreproduce the longitudinal invariance.Bessel beams have been previously generated in FDTD using two different approaches. In [25], Wu et al have simulated thegeneration of Bessel beam sources in FDTD, using total-field/scattered-field method [26]. In this approach, the computationaldomain is split into total-field and scattered-field regions. The electromagnetic wave is injected using the surface which separatesthe two regions. In [27], arbitrary order Bessel beams were generated using the scattered-field approach [26]. In this case, thetotal fields are decomposed into known incident fields and unknown scattered fields. The incident fields at all grid points areevaluated using the analytical expression at each time step.Previous techniques of Bessel beam generation have inherent limitations. Indeed, both total-field/scattered-field and scattered-field methods are inappropriate for implementation in PIC codes. In the scattered-field approach, there is no direct access to thetotal field which is required in PIC codes. Moreover, the incident field is needed at any point in the grid. The total-field/scattered-field needs an extra computational region for the scattered field. In this technique particle and field boundaries would be definedin different places because particles are in interaction with the total fields. In contrast with these approaches, injecting theelectromagnetic fields using an antenna has several advantages: (1) there is direct access to total field, (2) incident field is neededonly at the antenna points, (3) no extra computational region is needed for the scattered fields (because PIC approach entirelyworks with total field), and (4) particle and field boundaries are set in same place.In the present work, a cylindrical annulus which we call Bessel antenna is inscribed in the FDTD box and emits the Besselsolutions of Maxwell’s equations. In comparison with Bessel-Gauss FDTD simulation, it significantly reduces the size of thethree-dimensional (3D) computational box. It has been successfully tested for different orders of Bessel beam, different coneangles, different antenna thicknesses, and different antenna radii. The relative deviation between the fields from the Besselantenna and those from the analytical solution is in the range of 0.01-7.0 percent. We readily note that the more straightforwardapproach consisting of injecting the Bessel solutions from the computational walls (square symmetry) does not yield satisfactoryresults in terms of beam symmetry.The paper is organized as follows. We first generate Bessel-Gauss beam in FDTD simulation for reference. Then, we derivein Section the Bessel solutions of the Hertz vector potential. We also provide in this section the analytical model of Bessel pulseto which we will compare our numerical simulations. The Bessel antenna implementation is detailed in Section . The results ofthe simulations are discussed in Section , where we investigate the influence of the antenna parameters. BESSEL-GAUSS BEAM
In this section, we implement the generation of a Bessel-Gauss beam in FDTD. Experimentally, a Bessel-Gauss beam can becreated by focusing a Gaussian beam using an axicon lens [23, 24, 28, 29]. The axicon applies a phase Φ( r ) = − kr sin θ ontothe Gaussian beam, where r is the radial distance in the cylindrical coordinate, k = 2 π/λ is the wavenumber with λ the lasercentral wavelength and θ is the cone angle.For this example, we utilize a computational box of × × µ m . Using a normalization in terms of wavelength, itcan be expressed as λ xy × λ xy × λ z where λ xy = λ/ sin θ , and λ z = λ/ cos θ . We inject from the z min wall a Gaussianbeam polarized along x direction and propagating in the the positive z direction, on which we apply the cylindrically-symmetricphase Φ( r ) . We set the wavelength and cone angle . µ m , ◦ , respectively. The temporal amplitude profile of the beam is aGaussian function of width
100 fs
Full Width at Half Maximum (FWHM). The Gaussian beam waist is w = 6 µ m . We run thesimulation up to t run = 266 fs .The fluence distribution ( (cid:82) t run S d t where S is the magnitude of the Poynting vector) for the resulting Bessel-Gauss beam atthe end of the simulation is shown in Fig. 1. As one can see, a Bessel-Gauss beam with a propagation distance of only ≈ µ m already requires a window width of ≈ µ m . In this example, the ratio between the central lobe radius and the transversedimension is about 0.1. The interaction of this beam with a nanoscale plasma rod will typically drive electron plasma waves ona spatial scale of ∼ . λ . If we resolve the wavelength of the excited plasma waves with 20 grid cells, a PIC simulation willneed a grid of × × which is computationally very expensive. We have implemented the Bessel antenna as analternative approach that can create a Bessel beam in a smaller box (transparent white color box in Fig. 1). FIG. 1. A Bessel-Gauss beam from FDTD simulation. Fluence distribution in zx plane at y = 0 (left), and yx plane at location of themaximum for I (0 , , z ) (right). We aim to create a Bessel beam in the white box with a computational box of λ xy × λ xy × λ z which issmaller than Bessel-Gauss one by a factor of 64. MATHEMATICAL MODELBessel’s solutions
Using the antenna, we will generate ideal Bessel beams by injecting the solution of the electromagnetic wave equation incylindrical coordinates. Several authors have derived the Bessel solutions for Maxwell’s equations in a homogeneous, isotropic,and unmagnetized medium [30–32]. But their solutions were not fully general in terms of polarization states. We briefly recallhere the general formulas and then provide handy normalized forms. We start with the Hertz vector potential Π which reads[33]: ∇ Π − µ (cid:15) ∂ t Π = (1)There are generally two solutions (related to TE and TM modes in cylindrical symmetry) for the electromagnetic fieldssatisfying Eq. (1). These solutions take the following forms for a harmonic function as exp( − iωt ) : ( E , E ) = ( ∇ × ∇ × Π , i ω µ ∇ × Π ) (2a) ( H , H ) = ( − i ω (cid:15) ∇ × Π , ∇ × ∇ × Π ) (2b)The solutions can be determined using two functions given by: ( Q , R ) = ( ∇ × Π , ∇ × ∇ × Π ) (3)Hence: ( E , E ) = ( R , i k Q ) (4a) ( η H , η H ) = ( − i k (cid:15) r Q , R ) (4b)here k = ω/c denotes the wavenumber in the vacuum, η = µ c the impedance of free space, and (cid:15) r the relative permittivityof the propagating medium. For Bessel’s solutions of Eq. (1), one of the following Hertz potentials is usually employed [30–33]. Π x = ˆx J m ( K r ) e i m φ e i ( ω t − β z ) (5a) Π y = ˆy J m ( K r ) e i m φ e i ( ω t − β z ) (5b) Π z = ˆz J m ( K r ) e i m φ e i ( ω t − β z ) (5c)where K is the transverse component (in xy plane) of the wave vector, β the axial component (along the z axis) of thewave vector, m order of Bessel function, φ the azimuth angle, and ( ˆx , ˆy , ˆz ) are unit vectors of cartesian coordinates. We have,therefore, β + K = (cid:15) r k = k . To link with the cone angle θ described in the previous section, we have K = k sin θ and β = k cos θ .Substituting Eq. (5) into Eq. (3) and after some algebra, we can obtain ( Q , R ) vectors associated with the three potentials ( Π x , Π y , Π z ) . In Cartesian coordinates, they are given by:1. Π x Q = e i ( ω t − β z ) e i m φ ˆx ˆy i β J m ( K r ) ˆz K i (cid:2) J m+1 ( K r ) e i φ + J m − ( K r ) e − i φ (cid:3) (6a) R = e i ( ω t − β z ) e i m φ ˆx k + β J m ( K r ) + K (cid:2) J m+2 ( K r ) e i φ + J m − ( K r ) e − i φ (cid:3) ˆy i K (cid:2) J m − ( K r ) e − i φ − J m+2 ( K r ) e i φ (cid:3) ˆz i β K (cid:2) J m − ( K r ) e − i φ − J m+1 ( K r ) e i φ (cid:3) (6b)2. Π y Q = e i ( ω t − β z ) e i m φ − ˆx i β J m ( K r ) ˆy − ˆz K (cid:2) J m+1 ( K r ) e i φ − J m − ( K r ) e − i φ (cid:3) (7a) R = e i ( ω t − β z ) e i m φ ˆx i K (cid:2) J m − ( K r ) e − i φ − J m+2 ( K r ) e i φ (cid:3) ˆy k + β J m ( K r ) − K (cid:2) J m+2 ( K r ) e i φ + J m − ( K r ) e − i φ (cid:3) − ˆz β K (cid:2) J m+1 ( K r ) e i φ + J m − ( K r ) e − i φ (cid:3) (7b)3. Π z Q = e i ( ω t − β z ) e i m φ ˆx i K (cid:2) J m+1 ( K r ) e i φ + J m − ( K r ) e − i φ (cid:3) ˆy K (cid:2) J m+1 ( K r ) e i φ − J m − ( K r ) e − i φ (cid:3) ˆz (8a) R = e i ( ω t − β z ) e i m φ − ˆx i β K (cid:2) J m+1 ( K r ) e i φ − J m − ( K r ) e − i φ (cid:3) − ˆy β K (cid:2) J m+1 ( K r ) e i φ + J m − ( K r ) e − i φ (cid:3) ˆz K J m ( K r ) (8b)The first two solutions, i.e., Π x and Π y , specify Bessel beams with linear, and more generally elliptical polarizations (becausethe field has in general a z component), while the solution Π z generate the purely radial or azimuthal polarizations for m = 0 (using respectively E and E solutions of Eq. (4a)). One can exactly reproduce Eqs. (7), and (8) in [32] from Π y in Eqs. (7)and Eqs. (12a) and (12b) in [31] from Π z in Eqs. (8).For implementation, we modify the above solutions based on cone angle θ and the amplitude of the electric field E . We notethat E corresponds, in the case m = 0 , to the modulus of the electric field on the optical axis ( r = 0 ). The corresponding ( E , H ) fields for TE, and TM modes are given in Eq. (9a), and Eq. (9b), respectively. ( E , H ) = k − ( E k Q , H R ) e i ( ω t − β z ) e i m φ (9a) ( E , H ) = k − ( E R , H k Q ) e i ( ω t − β z ) e i m φ (9b)where H = n/η E (in terms of the energy density µ H = (cid:15) (cid:15) r E ), and n = (cid:15) / is the optical refractive index. FIG. 2. Bessel pulse from Fourier spectrum integration. Eq. (11) integration is calculated for cone angle of θ = 17 ◦ (left panel), θ = 25 ◦ (middle panel), and θ = 30 ◦ (right panel). In all cases, m = 0 , T = 60 fs , T = 180 fs , and ω = 2 . . The red dashed lines show ± c/ sin θ velocities. We show the amplitude of E y field in log-scale to enhance the X-shape of the pulse propagation. Bessel pulse
For the validation of the antenna-generated Bessel pulse, we will need a reference based on an analytical expression. Here, weexpress the Bessel pulse using the monochromatic field derived in section . We work hereafter with the vector potential polarizedalong the x axis ( Π x ) to generate a linearly y polarized Bessel pulse of arbitrary order. A pulse in time-domain is obtained byperforming the integral over the pulse spectrum F ( ω ) . We choose F ( ω ) as a Gaussian function that inverse Fourier transformgives F ( t ) = exp( − iω t ) exp( − t /T ) function in time domain. Hence: F ( ω ) = T √ π e − T ω − ω (10)The y component of the electric field then reads: E y ( x , t ) = − cos θ e i m φ (cid:90) ∞−∞ dωF ( ω ) J m ( ω sin θc r ) e iω ( t − T − z cos θ/c ) (11)Where T corresponds to the time of the peak intensity.We calculate the Eq. (11) integration for m = 0 , , , T = 60 fs ( FWHM = 2 T √ ln 2 = 100 fs ), T = 180 fs , and ω = 2 . . Three different cone angles of θ = 17 ◦ , θ = 25 ◦ , and θ = 30 ◦ are considered. Shown in Fig. 2 is the absolutevalue of E y (normalized to the maximum value) extracted from Eq. (11) for m = 0 . The y coordinate is scaled with thetransverse wavelength λ xy = λ/ sin θ . It, therefore, increases as the cone angle decreases from ◦ to ◦ . In all cases, onesees two propagating waves coming from the bottom ( t = 0 ), interfere around y = 0 , and then propagate away. The interferingwaves propagate with velocities of ± c/ sin θ (red dashed lines). In Section , the E y of the antenna generated Bessel pulse willbe compared with the corresponding Eq. (11) integration. We note again that the width over which we see the pulse propagationis much larger than the width of the central lobe. FIG. 3. Concept of the Bessel antenna in the FDTD simulation box. The yellow annulus with a thickness of δs shows the cylindrical antennawhich has an outer radius R B = N ⊥ λ xy , and a length N (cid:107) λ z . There are N PML grid cells in each transverse direction where the CPMLboundary condition is applied. TABLE I. Sets of simulation parameters for Bessel antenna.Run A B C D E F G H I m θ ◦
30 30 30 25 17 25 25 25 25 δs [ λ ] N (cid:107) N ⊥ BESSEL ANTENNA IMPLEMENTATION
To model the propagation of a plane wave in FDTD, we inject E and B solutions of the wave equation from one wall into thecomputational box. We use a similar methodology for the propagation of the Bessel beam. Since a Bessel beam has a cylindricalsymmetry, we inject the Bessel solution presented in Section from a cylindrical antenna into the simulation box.We implement the antenna in EPOCH code, which is a Particle-In-Cell (PIC) code for laser-plasma simulation [34]. TheFDTD algorithm in EPOCH uses a modified version of the leapfrog scheme in which the fields are updated at both the halftime-step and the full time-step. The full detail of the EPOCH code is presented in [34]. The standard second order Yee’s FDTDscheme is used in the present work, although higher-order schemes (4th- and 6th-order) are also available.The Bessel antenna is inscribed in the simulation box as illustrated in Fig. 3. The cylinder has an outer radius R B and athickness δs . To avoid reflection of the outgoing waves from the ( x min , x max ) and ( y min , y max ) boundaries, perfectly matchedlayers (PML) are used in these surfaces as shown in Fig. 3. A PML is an artificial absorbing layer. It is commonly used totruncate computational regions in numerical methods to simulate problems with open boundaries, particularly in the FDTD [26].In practice, we used the convolutional PML (CPML) method, as presented in [26, 35]. The current CPML implementation inEPOCH uses a number of PML layers ranging between 6 and 16. We use N PML = 6 in the current work.We define the sampling of our computational box with the following procedure. We sample the laser wavelength with n s grid cells. Here, we use n s = 20 . We set the transverse width of the box, excluding the PML boundary layers, as an evennumber N ⊥ of the transverse period of the field λ xy = 2 π/K . The number of points of the computational box is then definedas N xy = int(2 n s N ⊥ / sin θ ) + 2 N PML in the transverse plane. The longitudinal size of the box, including the PML boundarylayers, is an even number of the period λ z = 2 π/β and the corresponding number of points is N z = int(2 n s N (cid:107) / cos θ )+2 N PML .This ensures that the wave propagation is accurately sampled in both longitudinal and transverse directions. The cell size is alsosame in both directions so that the mesh is cubic. (We note, for EPOCH users, that in EPOCH input parameter file, the numberof grid points and length of the computational box have to be reduced by the number of CPML layers and corresponding widthof the CPML conditions.)We will focus on the generation of linearly polarized Bessel beams of order m with a Gaussian temporal profile. Hence,we use the vector potential polarized along the x axis. As shown in section , this will generate a linearly y polarized Besselbeam. Since the Fourier transformation operation cannot be straightforwardly implemented in EPOCH, we use the zeroth-orderapproximation to the integration of Eq. (11), which corresponds to multiplying the ( E , B ) solutions of Eq. (9a) or Eq. (9b) aGaussian time profile peaked at time T : f ( t ) = exp[ − ( t − T ) /T ] . Therefore, our input fields are defined by: E B = k − (cid:60){ E f ( t ) k Q e i [ ω ( t − T ) − β z ] e i m φ } (12a) H B = k − (cid:60){ H f ( t ) R e i [ ω ( t − T ) − β z ] e i m φ } (12b)Here we note that the zeroth-order approximation we performed is equivalent to neglecting the time needed by an illuminatingtemporally Gaussian pulse to reach the extremity of the cylinder in front of the pulse duration. In other words, we consider thatthe Gaussian pulse illuminates at the same time the front side ( z min ) and the rear side ( z max ). This is necessary to obtain theperiodic boundary conditions between these two planes. This approximation is fully justified if the computational box size ismuch smaller than the pulse duration divided by the speed of light. FIG. 4. Different-order Bessel beams. The amplitude of the y component electric field and the intensity of generated Bessel beams are plottedon xy plane at z = 0 : Run A (top row), Run B (middle row), and Run C (bottom row). The first three columns show the y component of theelectric field at different times, t = 26 . , . , and . from left to right. The rightmost column shows the intensity distributionat t = 213 . , the corresponding time of the peak intensity ( T + R B /c sin θ ). The white circle in the right column shows the inner radiusof the Bessel antenna . Finally, we need to define the procedure with which the fields will be updated at the antenna points. In the common sourceimplementation, which is referred to as hard source [26], the fields are defined in source points from predefined functions (likeEqs. 12) and the FDTD update equations do not apply to them. This kind of source implementation, however, scatters all incidentwaves on the source. In our case, however, the fields emitted out of a plasma placed inside the cylindrical antenna should notbounce back.Therefore, we use the soft source scheme [36] which drastically reduces the scattering of the incoming waves. We update thepoints in the antenna by the sum of the values from the FDTD Maxwell’s curl equations and the values dictated by the Eqs. 12.Hence E t+ δ ti ( i + δ i / , j, k ) = E t+ δ ti ( i + δ i / , j, k ) + E t+ δ tBi ( i + δ i / , j, k ) (13a) H t+ δ ti ( i, j + δ j / , k + δ k /
2) = H t+ δ ti ( i, j + δ j / , k + δ k / H t+ δ tBi ( i, j + δ j / , k + δ k / (13b)with cyclic permutation of the indices i , j , and k corresponding to x , y , and z , respectively. Here, δ i is the cell size in direction i . We note that we have independently investigated the opportunity of implementing a transparent source [36], which is basedon subtracting impulse response of the grid. However, given the high sampling of our grid ( n s = 20 grid cells per wavelength)and introducing the source with a smoothly increasing envelope, it appeared that the difference between the soft source andthe transparent one would be less than 5 percent. In contrast, the calculation of the convolution of the injected fields with theimpulse response at all times would be computationally very expensive, this is why we discarded this technical solution. Allresults presented afterwards were obtained using the soft source defined in Eq. (13). RESULTS
In this section, we present our results and compare them to the analytical one from Eq. (11). The antenna is tested for Besselorder of m = 0 , , and , and the cone angle of θ = 17 ◦ , θ = 25 ◦ , and θ = 30 ◦ . The different sets of parameters are summarizedin Table I. For all cases, T = 60 fs , T = 180 fs , and λ = 0 . µ m .We show in Fig. 4 in the three leftmost columns the amplitude of E y in xy plane at different times in the pulse for the FDTDgenerated Bessel beams of order m = 0 , , and . These are calculated for a cone angle of θ = 30 ◦ . The rows correspondrespectively to runs A, B, C where the order m was varied. We observe that the cylindrical symmetry is well preserved for m = 0 and that during the pulse build-up, the fields propagate inwards. The interference progressively takes place at the centerto create the Bessel beams of order m . For the cases of m = 1 , , the pattern effectively rotates with time. (We remark thatmaxima have the same orientation in the 3 times because they are separated by integer number of temporal periods.)In the rightmost column of Fig.4, we show the intensity distribution at a time t = 213 . , i.e. which corresponds to thetime of the peak intensity. The intensity is calculated using /τ (cid:82) t + τt dt S where τ is the period and S is the magnitude ofthe Poynting vector. The intensity distribution shows the expected transverse map for the different order Bessel beams. Thereis a high-intensity core in the zeroth-order Bessel beam while one can see the low-intensity core for the first-order and thesecond-order beams. The cylindrical symmetry is highly apparent, confirming the quality of the injection. FIG. 5. Different-order Bessel beams. The intensity map of generated Bessel beams in xz plane at y = 0 : Run A (left), Run B (middle), andRun C (right). All snapshots are taken at . . In Fig. 5 we show the intensity distribution in xz plane at a time t = 213 . . The intensity is actually invariant with z within ± percent, confirming the ”non-diffracting” behaviour. FIG. 6. Different-order Bessel beams. Comparisons between the y component electric field from the Bessel antenna (blue solid line) and fromthe Eq. (11) (red dashed line). Run A (left column), Run B (middle column), and Run C (right column). The bottom row shows the relativedeviation. In Fig.6, we quantitatively compare the FDTD generated E y field with the analytical expression of Eq. (11). The top rowcompares the absolute values of the fields at time t = 213 . for the different orders m = 0 , , . The bottom row shows therelative deviation defined as χ = | E y − E anay | /E anay in log-scale. In this expression, E y is the FDTD generated field and E anay is from the analytical expression. The relative deviation lies in the range of 0.01-7.0 percent over the full computation window.Interestingly, the error is reduced to less than 3 percent in the central region of the three main lobes which is obviously of higherinterest. Since the deviation increases toward the antenna location, we expect that better agreement might be achieved using alarger annulus radius.We have also investigated the quality of the beam generation for different cone angles. The results are shown for zeroth-orderBessel beam in Fig. 7. The radius of the antenna was scaled according to the transverse period (see Table 1), while using anidentical antenna thickness of λ . The increase in radius required longer computation time to let the pulse reach the centralregion. In all cases, the beams remain of high quality. Fig. 8 compares quantitatively the FDTD generated beam to the analyticalone. We observe that while the error remains small - and even decreases - for lower cone angles at the center, it increases onthe edges. More precisely, the number of constructed lobes decreases as the cone angle is reduced. This issue is still underinvestigation.We have investigated the effect of antenna thickness δs for a fixed cone angle of θ = 25 ◦ . We show the results in Fig. 9 (toprow) for three different thicknesses. As one can observe on the figure, the deviation is reduced to about 1.5 percent when theantenna thickness increases.Finally, the radius of the Bessel antenna R B is another important parameter since this is directly related to the computationaleffort. The bottom row of Fig.9 shows the deviation when the radius of the Bessel antenna is varied from 2 to 4 transversewavelengths λ xy , while maintaining the antenna thickness fixed. More lobes are obviously generated inside the antenna (withthe same transverse period) for the larger box, but there is not a major difference for the three central lobes. In fact, the deviationfor the three central lobes is similar in the three cases, but we note that the deviation increases for the outer lobes for the largestbox. Therefore, the case of λ xy offers an excellent compromise between accuracy in the central lobes and computational effort.We have investigated other components of the electric fields ( E x , E y ) as well. The relative deviation for the E z component issimilar as the deviation for E y . The E x component which is zero in the injected solutions is four orders of magnitude smaller0 FIG. 7. Bessel beam with different cone angles. The amplitude of the y component electric field and the intensity of generated Bessel beamsare plotted on xy plane at z = 0 : Run E (top row), Run D (middle row), and Run A (bottom row). The first three columns show the timeevolution of the y component electric field at three different times from left to right. As the simulation box is different for each row, thecorresponding time for each column is different. The rightmost column shows the intensity distribution. The white circle in the right columnshows the inner radius of the Bessel antenna. than the main component E y . Therefore, the electromagnetic fields generated by the Bessel antenna well satisfy Maxwell’s curlequations.1 FIG. 8. Bessel beam with different cone angles. Comparisons between the y component electric field from the Bessel antenna (blue solid line)and from the Eq. (11) (red dashed line). Run E (left column), Run D (middle column), and Run A (right column). The bottom row shows therelative deviation.FIG. 9. Effects of the antenna thickness and radius on the Bessel beam. The top row shows the relative deviation for different antennathicknesses: Run D (top left), Run F (top middle), and Run G (top right). The bottom row shows the relative deviation for different antennaradii: Run I (bottom left), Run H (bottom middle), and Run D (bottom right). CONCLUSIONS
The generation of Bessel beams in FDTD is an essential and challenging task to enable the investigation of the Bessel beam-plasma interaction with PIC codes. A common implementation is using the so-called Bessel-Gauss beam in which a Bessel beamis generated by focusing a Gaussian beam using an axicon lens, which has been demonstrated here for reference. However, dueto its computational cost, we have developed an alternative approach that is computationally more efficient. This solution isbased on using a cylindrical antenna emitting electromagnetic fields inwards and using the periodic boundary conditions alongthe optical axis.We have first provided the full set of Bessel’s solutions of Maxwell’s equations. The solutions have been then injected into theFDTD grid through the cylindrical antenna. We have validated our approach by comparing the field from the FDTD simulationwith that from the analytical solution. We have shown that different orders of Bessel beams can be successfully generated,with small deviations with regard to the analytical solution, for different beam cone angles. Better agreement was found for thehigher angles. We have investigated the effect of the different antenna parameters on the deviation and shown that thicker antennaprovide a better result, while the radius of the antenna impacts negligibly on the accuracy of central lobes. The electromagneticfields generated by the Bessel antenna well satisfy Maxwell’s curl equations. In conclusion, using a computational box with afull length of two longitudinal periods and a full width of four transverse periods provides good accuracy. In comparison withthe already short Bessel-Gauss beam that was generated, the computational effort has been decreased by a factor 64. Despiterestricted to study longitudinally periodic objects or longitudinally periodic particle distributions, we believe that our work willfind applications in the field of laser-particle scattering [37], particle trapping, nonlinear plasmonics [38] and specifically inlaser-particle acceleration and laser-plasma interaction [6, 8, 13, 19, 20].
FUNDING
European Research Council (ERC) (682032-PULSAR), R´egion Bourgogne Franche-Comt´e and the EIPHI Graduate School(ANR-17-EURE-0002).
ACKNOWLEDGMENTS
We thank the EPOCH support team for help, https://cfsa-pmw.warwick.ac.uk . Numerical simulations have beenperformed using the M´esocentre de Calcul de Franche-Comt´e, PRACE Research Infrastructure resource MARCONI-KNL basedat CINECA, Casalecchio di Reno, Italy, within the Project ”PULSARPIC” (PRA19 4980), and KDK computer system at Re-search Institute for Sustainable Humanosphere, Kyoto University.
DISCLOSURES
The authors declare no conflicts of interest. [1] J. Durnin, J. J. Miceli, and J. H. Eberly, Phys. Rev. Lett. , 1499 (1987).[2] H. Little, C. Brown, V. Garc´es-Ch´avez, W. Sibbett, and K. Dholakia, Opt. Express , 2560 (2004).[3] E. Mcleod and C. Arnold, Nature Nanotechnology , 413 (2008).[4] C. Vetter, R. Steinkopf, K. Bergner, M. Ornigotti, S. Nolte, H. Gross, and A. Szameit, Laser & Photonics Reviews , 1900103 (2019),https://onlinelibrary.wiley.com/doi/pdf/10.1002/lpor.201900103.[5] V. Garc´es-Ch´avez, D. McGloin, H. Melville, W. Sibbett, and K. Dholakia, Nature , 145 (2002).[6] B. Hafizi, E. Esarey, and P. Sprangle, Phys. Rev. E , 3539 (1997).[7] D. Li and K. Imasaki, Applied Physics Letters , 031110 (2005).[8] S. Kumar, A. Parola, P. Di Trapani, and O. Jedrkiewicz, Applied Physics B , 185 (2017).[9] F. O. Fahrbach, P. Simon, and A. Rohrbach, Nature Photonics , 780 (2010).[10] T. A. Planchon, L. Gao, D. E. Milkie, M. W. Davidson, J. A. Galbraith, C. G. Galbraith, and E. Betzig, Nature Methods , 417 (2011).[11] E. Gaizauskas, E. Vanagas, V. Jarutis, S. Juodkazis, V. Mizeikis, and H. Misawa, Opt. Lett. , 80 (2006).[12] D. Roskey, M. Kolesik, J. Moloney, and E. Wright, Opt. Express , 9893 (2007).[13] D. Faccio, E. Rubino, A. Lotti, A. Couairon, A. Dubietis, G. Tamoˇsauskas, D. G. Papazoglou, and S. Tzortzakis, Phys. Rev. A ,033829 (2012). [14] M. Duocastella and C. Arnold, Laser & Photonics Reviews , 607 (2012), https://onlinelibrary.wiley.com/doi/pdf/10.1002/lpor.201100031.[15] C. Xie, V. Jukna, C. Mili´an, R. Giust, I. Ouadghiri-Idrissi, T. Itina, J. M. Dudley, A. Couairon, and F. Courvoisier, Scientific Reports ,8914 (2015).[16] F. Courvoisier, R. Stoian, and A. Couairon, Optics & Laser Technology , 125 (2016).[17] R. Stoian, M. Bhuyan, G. Zhang, G. Cheng, R. Meyer, and F. Courvoisier, Advanced Optical Technologies , 165 (2018).[18] K. Bergner, M. M¨uller, R. Klas, J. Limpert, S. Nolte, and A. T¨unnerman, Appl. Opt. , 5941 (2018).[19] C. G. Durfee and H. M. Milchberg, Phys. Rev. Lett. , 2409 (1993).[20] G. A. Hine, A. J. Goers, L. Feder, J. A. Elle, S. J. Yoon, and H. M. Milchberg, Opt. Lett. , 3427 (2016).[21] J. M. Dawson, Rev. Mod. Phys. , 403 (1983).[22] C. Birdsall and A. Langdon, “Plasma physics via computer simulation,” (Taylor & Francis, 2004) pp. 351–385.[23] J. H. McLeod, J. Opt. Soc. Am. , 592 (1954).[24] D. McGloin and K. Dholakia, Contemporary Physics , 15 (2005).[25] Z. Wu, Y. Han, J. Wang, and Z. Cui, Opt. Express , 28727 (2018).[26] A. Taflove and S. Hagness, “Computational electrodynamics: The finite-difference time-domain method,” (Artech House, 2005) pp.169–327.[27] A. Chen, J. Wang, Y. Han, Z. Cui, and M. Yu, Opt. Express , 26766 (2018).[28] F. Gori, G. Guattari, and C. Padovani, Optics Communications , 491 (1987).[29] V. Jarutis, R. Paˇskauskas, and A. Stabinis, Optics Communications , 105 (2000).[30] S. Mishra, Optics Communications , 159 (1991).[31] Y. Yu and W.-B. Dou, Progress in Electromagnetics Research Letters , 57 (2008).[32] J. J. Wang, T. Wriedt, J. A. Lock, and L. M¨adler, Journal of Quantitative Spectroscopy and Radiative Transfer , 218 (2016).[33] J. Stratton, “Electromagnetic theory,” (McGraw-Hill, 2007) pp. 349–387.[34] T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R.Bell, and C. P. Ridgers, Plasma Physics and Controlled Fusion , 113001 (2015).[35] J. A. Roden and S. D. Gedney, Microwave and Optical Technology Letters , 334 (2000).[36] J. B. Schneider, C. L. Wagner, and O. M. Ramahi, IEEE Transactions on Antennas and Propagation , 1159 (1998).[37] Z. Cui, Y. Han, and L. Han, J. Opt. Soc. Am. A , 1913 (2013).[38] M. Kauranen and A. V. Zayats, Nature Photonics6