Theoretical Analysis of Airy–Gauss Bullets Obtained by Means of High Aperture Binary Micro Zonal Plate

We theoretically analyze the methodology for obtaining vectorial three-dimensional bullets, concretely Airy–Gauss bullets. To do this, binary micro zonal plates (BZP) were designed in order to obtain different Airy–Gauss bullets with sub-diffraction main lobe width. Following the vectorial diffraction theory, among the electrical field, we extend the theory to the magnetic field, and thus we analyze several properties such as the Poynting vector and the energy of Airy–Gauss vectorial bullets generated by illuminating the designed BZP with a temporal Gaussian circular polarized pulses.


Introduction
In recent years, due to their nondiffracting, self-accelerating, and self-healing properties, airy beams have attracted interest in linear and non-linear optics fields [1]. Currently, Airy beams are experimentally available, and, as a result, a wide interest has been paid in optical beam manipulation applications. Thus, by using the particular properties of Airy beams, several applications have been developed such as Airy plasma guiding [2], routing surface plasmon polaritons [3], image signal transmission [4], optical micromanipulation [5], optical bullets [1,6], optical trapping [7,8], and biomedical applications [9][10][11].
Experimentally, Airy beams have been generated by imposing a cubic phase on a Gaussian beam and then taking its Fourier transform by using a lens [1,12,13]. Concretely, several methods have been used to do this; thus, spatial light modulators (SLM) have been employed for obtaining a cubic phase profile. However, the beam quality is limited due to the discrete behavior of the device resulting in low efficiency. In a similar manner, liquid crystals and polymeric liquid crystals have also been used, but two Airy beams are simultaneously generated in this case [14,15]. Among these, cubic phase zonal phase plates [16][17][18] have been designed for obtaining high focusing capability.
Lastly, plasmonic and dielectric metasurfaces have transformed photonics research. Metasurfaces are planar versions of metamaterials and have shown unprecedented ability to arbitrarily manipulate light wavefronts, providing efficient solutions to generate desired accelerating beams [19,20]. These artificially designed 2D structured can be made of arrays of sub-wavelength scatterers for arbitrary tuning the characteristics of electromagnetic waves such as the control of Airy beams. Due to the size of the metasurface unit being much less than the wavelength, the limitations of conventional optical devices have relieved. As a result these new devices produce an increase in the bandwidth, generation efficiency, and integration capacity [21][22][23].
In many applications, it is vital that the light beams sharply autofocus their intensity right before a focal point, while keeping a low intensity until that very moment. To this end, the Airy beams were evolved into autofocusing beams [24,25] using a radially Airy symmetric intensity distribution. Generation and propagation of the these circular Airy beams with the abruptly autofocusing property were widely investigated [26][27][28] due to its potential value in optical micro-manipulation [29][30][31][32], and the excellent property makes them especially beneficial for biomedical treatment, laser ablation, and the generation of high-intensity lasers [11,33,34].
Moreover, among the applications of the Airy beams, the interest in the generation of temporal and spatial localized light pulses (optical bullets) has raised in recent years. In this sense, the Airy functions bring the possibility, in a linear regime, to separate spatial and temporal variables to obtain the 3D spatio-temporal beam [13,35,36]. Following, the theoretical procedure proposed by Siviloglou [13], Chong experimentally demonstrated the generation of light bullets in a linear regime, to be exact a spatial Bessel and a temporal Airy [37][38][39], and Abdollahpour obtained the Airy-Airy bullets [6]. Furthermore, theoretically, by solving in paraxial conditions the Schrödinger equation, different nondiffracting beams systems have been described such as Airy-parabolic-cylindrical [40], Airy-Hermite-Gauss [41], and Airy-Tricomi-Gauss [42].
Recently, based on the diffraction theory, we analyzed the generation of vectorial Airy-Airy light bullets by means of binary micro zonal plates (BZP). Thus, we developed a theoretical methodology for obtaining and analyzing Airy-Airy vectorial bullets [43]. In this work, based on the previous results, we studied Airy-Gauss light bullets. To do this, we extended the theoretical analysis to the magnetic field, which gives the possibility of studying additional properties such as the Poynting vector, energy, and the helicity of the resulting Airy-Gauss bullets.

Theoretical Description
As we have stated, Chong experimentally demonstrated the generation of light bullets in a linear regime, to be exact a spatial Bessel and temporal Airy [37][38][39]. The experimental setup, basically, can be divided in two parts (the temporal profile generation), and the second one is the spatial generation. In our case, instead of the case analyzed by Chong, in the first part (temporal profile generation), Gaussian temporal pulses were obtained from a mode-locked laser. In the second part, instead of an Axicon (to obtain a Bessel profile), the spatial Gaussian beam was converted into a spatial Airy beam by using a Binary Zone Plate. In this work we, theoretically, propose a BZP that could be used to obtain Airy-Gauss bullets. Applying the same methodology that we previously described [43], we design a binary zonal plate to focus real or complex transmittance T(x o , y o ) defined at the z = 0 plane in a focal plane located at z = f (see Figure 1). It is important to mention that in this work, the BZP is not fabricated; only the theoretical function is proposed, and to encode it in a material, a directing laser writing technique could be used [17,44]. Furthermore, according to Figure 1, it is assumed that the medium in z > 0 region is water, but air or other linear media can be used. The construction of zonal micro-plates with binary transmittance T ZM (x o , y o ) is based on the following function: where I() and Mod() are the round function and module 1 operator, respectively. According to this definition, T ZM will be a binary amplitude real function. Furthermore, the function ψ(x o , y o ) (Equation (1)) will simultaneously contain information about the transmittance of the binary zonal plate, T(x o , y o ), and a spherical convergent wave of wave number k 0 = ω 0 n/c that focuses at plane z = f according to [45]: In this work, to generate a spatial Airy beam, the designed zonal plates must fulfill that T(x o , y o ) = exp (iβ (x 3 o + y 3 o )). Thus, taking into account Equations (1) and (2), as an example, the transmittance of the generated microplate, with the design parameters summarized in Table 1, can be observed in Figure 1.  Table 1. (b) Diagram for the generation of Airy-Gauss bullets, comprising intensity iso-surfaces of the incident Gaussian beam and the generated spatio-temporal (Airy-Gauss bullet) by the designed binary microplate. Table 1. Design parameters of the binary micro zonal plate shown in Figure 1a.

Parameters Values
The electric field in the Cartesian coordinates solution to the Maxwell's equations in the frequency domain for the region z > 0 generated by the diffraction of an incident field on the designed microplate ( Figure 1b) is given by: where F(ω) is obtained from the Fourier transform of the incident temporal pulse, and the field components, according to Luneburg's vectorial diffraction theory e(x, y, z, ω) [46], are: to the fields generated by the zonal plate, which are related to the polarization state of the incident field (E xo ,E yo ) and also characterized by the designed binary transmittance T ZM : In the same way, the corresponding magnetic field in the Cartesian coordinates solution to the Maxwell's equations in the frequency domain for the region z > 0 is given by: h(x, y, z, ω) can be obtained through the Faraday's law (ikη h(x, y, z, w) = ∇ × e(x, y, z, w)) where η is the impedance in the propagation medium. Thus, the magnetic field components can obtained by introducing Equation (4) in Faraday's law: In this work, we analyzed temporal Gaussian profiles, so we used a F(ω) function given by: where ω 0 corresponds to the design frequency of the microplate Equation (2), and we used the parameter value b = 7.21 × 10 −31 ((rad/s) −2 ), which corresponds to an initial temporal Gaussian pulse of δt = 20 fs. Thus, the Airy-Gauss bullets were obtained by illuminating the microplate by an electric field with a temporal Gaussian profile. The numerical procedure is detailed as follows: 1. The microplate is illuminated by an incident circular electrical field given by a polar- The spectral bandwidth of the incident pulse, ∆ω, is discretized At each z-plane the spatial distributions of the non-diffracting beam ( e(x, y, z, ω i ) and h(x, y, z, ω i )) is numerically obtained by introducing Equations (1) and (2) with the above-mentioned incident electrical field in Equation (4) (for the electrical field) and in Equation (7) (for the magnetic field).

4.
For computational purposes, by taking into account the Convolution theorem, Equations (4) and (7) are numerically solved by means of fast Fourier transform operations by using [47,48]: 5.
Finally, the Airy-Gauss bullets at each z plane are obtained by the discrete inverse Fourier transform by using the spatial distributions of the non-diffracting beam at each frequency ( e(x, y, z, ω i ) and h(x, y, z, ω i )): It is important to note that the validity of the numerical treatment used in this work is widely accepted as it has been compared to commercially available FDTD software [49][50][51]. Moreover, due to the dimensions of the proposed systems, the methodology used in this work presents significant advantages over 3D FDTD (accuracy, reliability, and computational cost).

Results and Discussion
First of all, we are analyzed the averaged total intensity at the studied temporal range at each plane (80 fs). To do this, in order to characterize the electromagnetic field distribution, we defined a normalized temporal averaged electromagnetic intensity to the maximum at focal plane as: where F , E, or H are for electrical and magnetic fields, respectively. Moreover, each temporal averaged vectorial component can be obtained by using in the numerator of Equation (11) only the corresponding component.
In Figure 2, we analyzed the electrical and magnetical temporal average intensity at z = f − 2λ 0 . As in all the cases, the averaged intensity was lower than the maximum value reached at the focal. At this plane, the Airy beam was not completely formed, but the total intensity of the main lobe was close to 15% of the maximum intensity reached. The transversal components (Figure 2c,e) are identical, (as can be deduced from Equation (4)), and both contribute to 80% of the total intensity (Figure 2a). The longitudinal component (Figure 2g) reached the 20% of the total intensity at this plane, presenting high deformation at the main lobe and at the side lobes as well. Moreover, by comparing to the transversal components, an offset in the positions of maximums of all the lobes appeared.
In a similar manner, following the same methodology, the corresponding magnetic field was analyzed at this figure. As it occurs with the electric field, the Airy beam is not totally formed at this plane, presenting lower-intensity values with respect to the maximum reached by the electrical field (Figure 2b). The main difference in the behavior of the electrical and magnetic field is given by the transversal components. As can be seen in Figure 2d,f, H x and H y are quite different, (as can be deduced from Equation (7)), mainly in the side lobes of both perpendicular directions; the H x component presents higher intensity at the x = 0 lobes whereas H y does at the y = 0 axis. The H z component is different to the others, but, in contrast, it is quite similar to E z (Figure 2g,h), with its contribution to the total intensity in this plane near 12%.
At the focal plane (Figure 3), the Airy beam is completely formed. The total electromagnetic intensity (Figure 3a,b) of the main lobe was raised, and the main corner lobes contained a large percentage of the beam's total intensity. It is important to note that at the focal the main lobe of the temporal averaged electrical and magnetic Airy beam presented a subwavelength FWHM value. It is also worth mentioning that the raise (double with respect to z = f − 2λ 0 plane) of the total electromagnetic intensity was mainly given by the corresponding amount of the transversal components one (Figure 3c-f). As previously described, the electric transversal components were quite similar among them, but the magnetic ones were different. At the focal plane, it can be clearly seen that the H x component presents higher intensity at the x = 0 lobes, whereas H y does so at the y = 0 axis. Finally, in relation to the intensity of E z and H z , the values of the intensities are similar to the obtained at the z = f − 2λ 0 plane (Figure 2g,h), but the intensity structure is too different to the other plane, because in this case both components resemble a Vortex-Airy beam [52].  This Vortex-Airy beam behavior can be explained by transforming the electromagnetic field given by Equations (4) and (7) to polar coordinates. Thus, taking into account that the incident field is circularly polarized and at plane z = 0 (x o = ρ cos(θ), y o = ρ sin(θ)) and at any z plane (x = r cos(ψ), y = r, sin(ψ)), the electric field given by Equation (4) is written as: being in this case R = [r 2 + ρ 2 + z 2 − 2rρ cos(θ − ψ)] 1/2 , and the function Ψ is given by: As it can be observed, the second term of function Ψ(E x (ρ, θ)exp(±iθ)) introduces a topological charge term ±1 to the generated fields by the zonal plate E x (ρ, θ), which results in a vortex beam behavior of e z component.
In the same way, the axial component of the magnetic field can be written as: Therefore, h z has the same topological charge as the axial component of the electric field and will present a vortex beam behavior.
Finally, the temporal average intensity of the electric and magnetic vectorial components are analyzed in Figure 4 at z = f + 2λ 0 plane. As can be seen, the behavior is not symmetrical with respect to the focal plane, with the structure of the Airy beam being maintained, but, the main lobe began to deform. The value of the total intensity decreased with respect to the focal, as before, mainly given by the transversal components. Regarding, transversal field components, the behavior was the same as above; the electrical components were quite similar among them, and the H x component presented higher intensity at the x = 0 lobes whereas H y did at the y = 0 axis. It is important to remark that, in contrast to the observed behavior at z = f − 2λ 0 , the structure form of these components corresponds to a lightly deformed Airy beam. However, the longitudinal components are quite similar between them and also to the obtained at z = f − 2λ 0 , in both intensity and structure form.
In relation to this analysis, from the normalized intensities of the electrical and magnetic fields, in Figure 5 we show the deflection along the propagation distance for all the field components. A parabolic trajectory can be clearly seen for both field and each component, which agree to the previous theoretical descriptions [17,53]. Moreover, due to the Airy vortex beam's similar behavior of the axial components (E z and H z ), the deflection is larger.  Once we performed the analysis of the average intensities, we analyzed the corresponding Airy-Gauss vectorial bullets obtained by the proposed methodology described above. At this point, the magnetic intensity was studied, but similar results were obtained with the electrical one. In Figure 6a, we present the temporal variation of the maximum magnetic intensity at each propagation distance. As can be seen, the maximum intensity of the temporal Gaussian pulse was approximately located at a range of ±λ 0 from the focal plane and obtained for times t ∼ 40 fs, a value that was delayed by diffraction effects with respect to the expected time of f n c = 35 fs. Moreover, 50% of the peak intensity can be obtained in the range between −1 to 2 λ 0 .
The temporal evolution of the intensity distributions of the Airy-Gauss bullets at several z planes is shown in Figure 6b-f. As mentioned above, the asymmetrical behavior with respect to the focal plane is clearly demonstrated. At a distance z = f − 2λ 0 (Figure 6b), the spatial structure of the Airy beam appeared with deformed low-intensity secondary lobes, and the main temporal-spatial lobe reached intensity values of 40%. Near the focal plane (Figure 6c), the spatial Airy beam was formed (presenting a well-defined ellipsoid form for all the lobes) with intensity close to 80% at the main spatial lobe. At the focal plane (Figure 6d), the maximum intensity value was reached at the main temporal-spatial lobe. Both secondary spatio-temporal peaks contained around 50% of the power, the spatial Airy being maintained up to 50 fs. Finally, at longer distances from the focal plane, (Figure 6f), the intensity of the central and secondary lobes decreased, but their values were higher than the corresponding values at z < f . Moreover, at z = f + 2λ 0 (Figure 6f), deformations in the spatio-temporal ellipsoidal structures of all the lobes began. As we stated above, the transverse components of the electromagnetic field present the same spatio-temporal structure as the total intensity shown in Figure 6. Moreover, the intensity change along the propagation direction is mainly produced by these components. However, the behavior of the electromagnetic longitudinal components is different to the transversal ones. Due to the fact that H z component response is similar to the corresponding electrical one, we analyzed in Figure 7 the longitudinal magnetic intensity of the Airy-Gauss bullets. In Figure 7a, the temporal variation in the maximum magnetic intensity at each propagation distance is shown. As can be seen, the maximum intensity of the temporal Gaussian pulse was approximately located at all the analyzed range being the higher values obtained at z < f .
The corresponding temporal evolution of the intensity distribution of the longitudinal magnetic component of the Airy-Gauss bullets at several z planes is shown in Figure 7b-f. At a distance z = f − 2λ 0 (Figure 7b), the spatial structure of the Airy beam was not formed (both main and secondary lobes). Near to the focal plane (Figure 7c), the spatial Airy beam was formed, presenting the secondary lobes a well-defined ellipsoid form with a deformed main lobe that presents null intensity similar to an Airy vortex beam. Finally, at longer distances from the focal plane, (Figure 7f), the ellipsoidal form of the secondary lobes was maintained, but the main lobe was deformed.
Finally, taking the advantage of this vectorial analysis, we analyzed several magnitudes such as the Poynting vector, the energy, and the Helicity of the Airy-Gauss bullets at each plane. In this sense, the temporal averaged values of these magnitudes were calculated by: G (x, y, z, ω i )δω (15) where G corresponds to S (Poynting vector); U (total energy) and H (helicity). Thus, S(x, y, z, ω i ) was calculated from the electromagnetic fields e(x, y, z, ω i ) and h(x, y, z, ω i ) obtained from Equations (4) and (7) as: The total energy was calculated according to U(x, y, z, ω i ) = 1 4 ( e(x, y, z, ω i ) · e * (x, y, z, ω i ) + µ h(x, y, z, ω i ) · h * (x, y, z, ω i )) (17) where and µ correspond to the electrical and magnetic permittivity, respectively, of the propagation media. For the analysis of the Poynting vector in Airy-Gauss bullets, in Figure 8 we show the temporal averaged transversal components at different propagation distances. It is important to note that the longitudinal component is an order of magnitude higher than these ones, but, we centered our study on the transversal components. As it can be seen, the magnitude of the transversal Poynting vector was higher at the lobes, where a circulation related to polarization state was observed. Among the circulation at the main lobes, at z = f − 2λ 0 , there was a preferential direction of 45 • pointing from the main lobe with positive slope (sense where the auto-acceleration is produced). At the focal, there was no preferential direction of the net energy flow, but the magnitude was the highest one, and a significant circulation was observed around the main and secondary lobes of the Airy structure. At z > f , the energetic circulation around the main lobes was maintained, but a preferential direction of 45 • appeared pointed in the opposite direction than at z < f (towards the main lobe), being the magnitude of the transversal Poynting vector higher at shorter distances from the focal.  Regarding to the energy, in Figure 9, the temporal averaged normalized energy with respect to the maximum value at the focal plane was analyzed. The value of the energy was maximum at the focal plane, as can be seen in Figure 9a, at z = f − 2λ 0 , the difference among the secondary to the main corners lobes was around the 50%. At the focal, the highest energy peak was at the main lobe, the difference between the main corners lobes being higher than at z < f . Finally, at longer distances from the focal, the energy decreased as well the difference between the main and secondary lobes.
Finally, together with energy, helicity is an important property of light, although it is closely related to its spin, and the polarization of light. In this sense, we obtained that this magnitude conserves along the propagation direction.

Conclusions
In conclusion, based on the vectorial diffraction theory we presented a methodology for obtaining and analyzing vectorial non-diffracting bullets. Thus, we proposed a design methodology to obtain binary zonal microplates (BZP) for generating Airy-Gauss bullets. For near field applications, based on an exact solution to Maxwell's equations, vectorial Airy-Gauss bullets were analyzed, developing a methodology for the study of the corresponding magnetic field. In this article, we presented circular polarized Airy-Gauss bullets generated from Gaussian temporal pulses of 20 fs obtained by the designed microplate of high numerical aperture, resulting in Airy-Gauss bullets with a FWHM close to 0.75 λ 0 . A complete analysis of the electromagnetic propagation was performed, paying special attention to the magnetic components (both transversal and longitudinal). It was observed that the intensity variation in the Airy-Gauss bullets was mainly given by the transversal components and that the contribution of the longitudinal one was maintained along the propagation distance. Taking into account the complete electromagnetic study, several magnitudes were obtained, such as the temporal averaged Poynting vector and the energy of the generated bullets. The transverse temporal averaged Poynting vector presents a circulation in the Airy lobes, the highest magnitude being obtained at the focal.