Generation of a Bessel beam in FDTD using a cylindrical antenna

Bessel beams are becoming a very useful tool in many areas of optics and photonics, because of the invariance of their intensity profile over an extended propagation range. Finite-Difference-Time-Domain (FDTD) approach is widely used for the modeling of the beam interaction with nanostructures. However, the generation of the Bessel beam in this approach is a computationally challenging problem. In this work, we report an approach for the generation of the infinite Bessel beams in three-dimensional FDTD. It is based on the injection of the Bessel solutions of Maxwell's equations from a cylindrical hollow annulus. This configuration is compatible with Particle In Cell simulations of laser plasma interactions. This configuration allows using a smaller computation box and is therefore computationally more efficient than the creation of a Bessel-Gauss beam from a wall and models more precisely the analytical infinite Bessel beam. Zeroth and higher-order Bessel beams with different cone angles are successfully produced. We investigate the effects of the injector parameters on the error with respect to the analytical solution. In all cases, the relative deviation is in the range of 0.01-7.0 percent.


INTRODUCTION
A Bessel beam refers to a family of solutions for the wave equation in which the amplitude of field is expressed by the Bessel function of the first kind [1]. It has attracted great interest in various branches in optics because its intensity profile is invariant as it propagates. Bessel beams are mainly utilized as optical traps [2][3][4], for optical manipulation [5], optical acceleration [6][7][8], light-sheet microscopy [9,10], nonlinear optics, ultrashort pulse filamentation [11][12][13], and laser-material processing [14][15][16][17][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 laser beams 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 are an interference field in which the propagation length Z B is defined by the transverse extent W and the cone angle θ, i.e. the angle made by the interfering waves with the optical axis: Z B ∼ W/ tan θ [23,24]. Hence, long propagation distances require a wide transverse extent. However, the central spot radius (r 0 = 0.383λ/ sin θ for a zeroth-order Bessel beam) is generally much smaller than the beam transverse extent. It is therefore computationally extremely demanding to investigate the laser-plasma interaction within the central lobe with high spatial resolution. This problem is illustrated in Fig. 1, where we show the timeintegrated fluence distribution for a finite-energy zeroth-order Bessel-Gauss femtosecond pulse over a distance of only ≈ 20 µ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. This is possible because in a Bessel beam, energy flows from the sides with a conical structure and because of the longitudinal invariance of both the diffraction-free Bessel beams and the uniform plasma distribution. Therefore, we will use a cylindrical annulus injecting the electromagnetic fields and use periodic boundary condition for the surfaces parallel to the optical axis to reproduce the longitudinal invariance.
Bessel beams have been previously generated in FDTD using two different approaches. In [25], Wu et al have simulated the generation of Bessel beam sources in FDTD, using total-field/scattered-field method [26]. In this approach, the computational domain is split into total-field and scattered-field regions. The electromagnetic wave is injected using the surface which separates the two regions. In [27], arbitrary order Bessel beams were generated using the scattered-field approach [26]. In this case, the total fields are decomposed into known incident fields and unknown scattered fields. The incident fields at all grid points are evaluated using the analytical expression at each time step.
Previous techniques of Bessel beam generation have inherent limitations. Indeed, both total-field/scattered-field and scatteredfield methods are inappropriate for implementation in PIC codes. In the scattered-field approach, there is no direct access to the total field which is required in PIC codes. Moreover, the incident field is needed at any point in the grid. The total-field/scatteredfield needs an extra computational region for the scattered field. In this technique particle and field boundaries would be defined in different places because particles are in interaction with the total fields. In contrast with these approaches, injecting the electromagnetic fields using an antenna has several advantages: (1) there is direct access to total field, (2) incident field is needed only at the antenna points, (3) no extra computational region is needed for the scattered fields (because PIC approach entirely works 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 Bessel solutions of Maxwell's equations. In comparison with Bessel-Gauss FDTD simulation, it significantly reduces the size of the three-dimensional (3D) computational box. It has been successfully tested for different orders of Bessel beam, different cone angles, different antenna thicknesses, and different antenna radii. The relative deviation between the fields from the Bessel antenna and those from the analytical solution is in the range of 0.01-7.0 percent. We readily note that the more straightforward approach consisting of injecting the Bessel solutions from the computational walls (square symmetry) does not yield satisfactory results in terms of beam symmetry.
The paper is organized as follows. We first generate Bessel-Gauss beam in FDTD simulation for reference. Then, we derive in Section the Bessel solutions of the Hertz vector potential. We also provide in this section the analytical model of Bessel pulse to which we will compare our numerical simulations. The Bessel antenna implementation is detailed in Section . The results of the 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 be created by focusing a Gaussian beam using an axicon lens [23,24,28,29]. The axicon applies a phase Φ(r) = −kr sin θ onto the Gaussian beam, where r is the radial distance in the cylindrical coordinate, k = 2π/λ is the wavenumber with λ the laser central wavelength and θ is the cone angle.
For this example, we utilize a computational box of 15 × 15 × 28 µm 3 . Using a normalization in terms of wavelength, it can be expressed as 8λ xy × 8λ xy × 32λ z where λ xy = λ/ sin θ, and λ z = λ/ cos θ. We inject from the z min wall a Gaussian beam polarized along x direction and propagating in the the positive z direction, on which we apply the cylindrically-symmetric phase Φ(r). We set the wavelength and cone angle 0.8 µm, 25 • , respectively. The temporal amplitude profile of the beam is a Gaussian function of width 100 fs Full Width at Half Maximum (FWHM). The Gaussian beam waist is w 0 = 6 µm. We run the simulation up to t run = 266 fs.
The fluence distribution ( trun 0 S dt where S is the magnitude of the Poynting vector) for the resulting Bessel-Gauss beam at the end of the simulation is shown in Fig. 1. As one can see, a Bessel-Gauss beam with a propagation distance of only ≈ 20 µm already requires a window width of ≈ 15 µm. In this example, the ratio between the central lobe radius and the transverse dimension is about 0.1. The interaction of this beam with a nanoscale plasma rod will typically drive electron plasma waves on a spatial scale of ∼ 0.05λ. If we resolve the wavelength of the excited plasma waves with 20 grid cells, a PIC simulation will need a grid of 7500 × 7500 × 14100 which is computationally very expensive. We have implemented the Bessel antenna as an alternative 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 the maximum for I(0, 0, z) (right). We aim to create a Bessel beam in the white box with a computational box of 4λxy × 4λxy × 2λz which is smaller than Bessel-Gauss one by a factor of 64.

Bessel's solutions
Using the antenna, we will generate ideal Bessel beams by injecting the solution of the electromagnetic wave equation in cylindrical coordinates. Several authors have derived the Bessel solutions for Maxwell's equations in a homogeneous, isotropic, and unmagnetized medium [30][31][32]. But their solutions were not fully general in terms of polarization states. We briefly recall here the general formulas and then provide handy normalized forms. We start with the Hertz vector potential Π which reads [33]: There are generally two solutions (related to TE and TM modes in cylindrical symmetry) for the electromagnetic fields satisfying Eq. (1). These solutions take the following forms for a harmonic function as exp(−iωt): The solutions can be determined using two functions given by: Hence: here k 0 = ω/c denotes the wavenumber in the vacuum, η = µ 0 c the impedance of free space, and r the relative permittivity of the propagating medium. For Bessel's solutions of Eq. (1), one of the following Hertz potentials is usually employed [30][31][32][33].
where K is the transverse component (in xy plane) of the wave vector, β the axial component (along the z axis) of the wave vector, m order of Bessel function, φ the azimuth angle, and (x,ŷ,ẑ) are unit vectors of cartesian coordinates. We have, therefore, 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: The first two solutions, i.e.,Π x and Π y , specify Bessel beams with linear, and more generally elliptical polarizations (because the field has in general a z component), while the solution Π z generate the purely radial or azimuthal polarizations for m = 0 (using respectively E 1 and E 2 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 0 . We note that E 0 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.

Bessel pulse
For the validation of the antenna-generated Bessel pulse, we will need a reference based on an analytical expression. Here, we express the Bessel pulse using the monochromatic field derived in section . We work hereafter with the vector potential polarized along the x axis (Π x ) to generate a linearly y polarized Bessel pulse of arbitrary order. A pulse in time-domain is obtained by performing the integral over the pulse spectrum F (ω). We choose F (ω) as a Gaussian function that inverse Fourier transform gives F (t) = exp(−iω 0 t) exp(−t 2 /T 2 ) function in time domain. Hence: The y component of the electric field then reads: Where T 0 corresponds to the time of the peak intensity. We calculate the Eq.  Fig. 2 is the absolute value of E y (normalized to the maximum value) extracted from Eq. (11) for m = 0. The y coordinate is scaled with the transverse wavelength λ xy = λ/ sin θ. It, therefore, increases as the cone angle decreases from 30 • to 17 • . In all cases, one sees two propagating waves coming from the bottom (t = 0), interfere around y = 0, and then propagate away. The interfering waves propagate with velocities of ±c/ sin θ (red dashed lines). In Section , the E y of the antenna generated Bessel pulse will be compared with the corresponding Eq. (11) integration. We note again that the width over which we see the pulse propagation is much larger than the width of the central lobe.

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 the computational box. We use a similar methodology for the propagation of the Bessel beam. Since a Bessel beam has a cylindrical symmetry, 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]. The FDTD algorithm in EPOCH uses a modified version of the leapfrog scheme in which the fields are updated at both the half time-step and the full time-step. The full detail of the EPOCH code is presented in [34]. The standard second order Yee's FDTD scheme 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 a thickness δs. To avoid reflection of the outgoing waves from the (x min , x max ) and (y min , y max ) boundaries, perfectly matched layers (PML) are used in these surfaces as shown in Fig. 3. A PML is an artificial absorbing layer. It is commonly used to truncate 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 in EPOCH 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 even number 2N ⊥ of the transverse period of the field λ xy = 2π/K. The number of points of the computational box is then defined as N xy = int(2n s N ⊥ / sin θ) + 2N PML in the transverse plane. The longitudinal size of the box, including the PML boundary layers, is an even number of the period λ z = 2π/β and the corresponding number of points is N z = int(2n s N / cos θ)+2N PML . This ensures that the wave propagation is accurately sampled in both longitudinal and transverse directions. The cell size is also same in both directions so that the mesh is cubic. (We note, for EPOCH users, that in EPOCH input parameter file, the number of grid points and length of the computational box have to be reduced by the number of CPML layers and corresponding width of 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 Bessel beam. Since the Fourier transformation operation cannot be straightforwardly implemented in EPOCH, we use the zeroth-order approximation to the integration of Eq. (11), which corresponds to multiplying the (E, B) solutions of Eq. (9a) or Eq. (9b) a Gaussian time profile peaked at time T 0 : f (t) = exp[−(t − T 0 ) 2 /T 2 ]. Therefore, our input fields are defined by: Here we note that the zeroth-order approximation we performed is equivalent to neglecting the time needed by an illuminating temporally Gaussian pulse to reach the extremity of the cylinder in front of the pulse duration. In other words, we consider that the Gaussian pulse illuminates at the same time the front side (z min ) and the rear side (z max ). This is necessary to obtain the periodic boundary conditions between these two planes. This approximation is fully justified if the computational box size is much smaller than the pulse duration divided by the speed of light. Finally, we need to define the procedure with which the fields will be updated at the antenna points. In the common source implementation, which is referred to as hard source [26], the fields are defined in source points from predefined functions (like Eqs. 12) and the FDTD update equations do not apply to them. This kind of source implementation, however, scatters all incident waves on the source. In our case, however, the fields emitted out of a plasma placed inside the cylindrical antenna should not bounce back.
Therefore, we use the soft source scheme [36] which drastically reduces the scattering of the incoming waves. We update the points 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 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 based on 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 and the transparent one would be less than 5 percent. In contrast, the calculation of the convolution of the injected fields with the impulse response at all times would be computationally very expensive, this is why we discarded this technical solution. All results 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 Bessel order of m = 0, 1, and 2, and the cone angle of θ = 17 • , θ = 25 • , and θ = 30 • . The different sets of parameters are summarized in Table I. For all cases, T = 60 fs, T 0 = 180 fs, and λ = 0.8 µ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 FDTD generated Bessel beams of order m = 0, 1, and 2. These are calculated for a cone angle of θ = 30 • . The rows correspond respectively 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 center to create the Bessel beams of order m. For the cases of m = 1, 2, the pattern effectively rotates with time. (We remark that maxima 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.5 fs, i.e. which corresponds to the time of the peak intensity. The intensity is calculated using 1/τ t+τ t dt S where τ is the period and S is the magnitude of the Poynting vector. The intensity distribution shows the expected transverse map for the different order Bessel beams. There is a high-intensity core in the zeroth-order Bessel beam while one can see the low-intensity core for the first-order and the second-order beams. The cylindrical symmetry is highly apparent, confirming the quality of the injection. In Fig. 5 we show the intensity distribution in xz plane at a time t = 213.5 fs. The intensity is actually invariant with z within ±2 percent, confirming the "non-diffracting" behaviour. Interestingly, the error is reduced to less than 3 percent in the central region of the three main lobes which is obviously of higher interest. Since the deviation increases toward the antenna location, we expect that better agreement might be achieved using a larger annulus radius. We have also investigated the quality of the beam generation for different cone angles. The results are shown for zeroth-order Bessel beam in Fig. 7. The radius of the antenna was scaled according to the transverse period (see Table 1), while using an identical antenna thickness of 1 λ. The increase in radius required longer computation time to let the pulse reach the central region. In all cases, the beams remain of high quality. Fig. 8 compares quantitatively the FDTD generated beam to the analytical one. We observe that while the error remains small -and even decreases -for lower cone angles at the center, it increases on the edges. More precisely, the number of constructed lobes decreases as the cone angle is reduced. This issue is still under investigation.
We have investigated the effect of antenna thickness δs for a fixed cone angle of θ = 25 • . We show the results in Fig. 9 (top row) for three different thicknesses. As one can observe on the figure, the deviation is reduced to about 1.5 percent when the antenna thickness increases.
Finally, the radius of the Bessel antenna R B is another important parameter since this is directly related to the computational effort. The bottom row of Fig.9 shows the deviation when the radius of the Bessel antenna is varied from 2 to 4 transverse wavelengths λ xy , while maintaining the antenna thickness fixed. More lobes are obviously generated inside the antenna (with the same transverse period) for the larger box, but there is not a major difference for the three central lobes. In fact, the deviation for the three central lobes is similar in the three cases, but we note that the deviation increases for the outer lobes for the largest box. Therefore, the case of 2λ 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 is similar as the deviation for E y . The E x component which is zero in the injected solutions is four orders of magnitude smaller

CONCLUSIONS
The generation of Bessel beams in FDTD is an essential and challenging task to enable the investigation of the Bessel beamplasma interaction with PIC codes. A common implementation is using the so-called Bessel-Gauss beam in which a Bessel beam is generated by focusing a Gaussian beam using an axicon lens, which has been demonstrated here for reference. However, due to its computational cost, we have developed an alternative approach that is computationally more efficient. This solution is based on using a cylindrical antenna emitting electromagnetic fields inwards and using the periodic boundary conditions along the optical axis.
We have first provided the full set of Bessel's solutions of Maxwell's equations. The solutions have been then injected into the FDTD grid through the cylindrical antenna. We have validated our approach by comparing the field from the FDTD simulation with 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 the higher angles. We have investigated the effect of the different antenna parameters on the deviation and shown that thicker antenna provide a better result, while the radius of the antenna impacts negligibly on the accuracy of central lobes. The electromagnetic fields generated by the Bessel antenna well satisfy Maxwell's curl equations. In conclusion, using a computational box with a full length of two longitudinal periods and a full width of four transverse periods provides good accuracy. In comparison with the already short Bessel-Gauss beam that was generated, the computational effort has been decreased by a factor 64. Despite restricted to study longitudinally periodic objects or longitudinally periodic particle distributions, we believe that our work will find applications in the field of laser-particle scattering [37], particle trapping, nonlinear plasmonics [38] and specifically in laser-particle acceleration and laser-plasma interaction [6,8,13,19,20].