## Abstract

We put forward a new optical system, which is composed of an existing axicon doublet and a newly proposed amplitude filter. The axicon doublet consists of a positive axicon and a negative axicon with high and low refractive indices, respectively. The Bessel beam generated by the axicon doublet propagates as far as more than 200 meters, owing to a small refractive index difference between the double axicons. The newly proposed amplitude filter is used to flatten the axial intensity distribution. Numerical results calculated by the complete Rayleigh-Sommerfeld method demonstrate that the generated Bessel beam propagates stably within a very long axial range. The proposed optical system is expected to have practical applications in tracking far-distance moving targets.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The Bessel beam, as a kind of nondiffraction beams, possesses propagating invariance [1] as well as self-healing [2] and self-accelerating [3] properties along the propagation direction. Therefore, it has broad applications in precise laser machining [4,5], atom trapping and particle manipulation [6,7], optical imaging [8] and microscopy [9–11], optical bullets [12], quantum information [13], volume hologram recording [14], optical communications [15], nonlinear optics [16] and so on. Many previous research works were focused on its generation [17,18] and propagation [19,20]. There typically include the following techniques for generating a Bessel beam: ring-slit method [21], computer generated hologram [22,23], volume holographic method [24], axicon [25,26], spherical-aberration lens [27], resonant cavity [28], and meta-surfaces [17,29].

Among all the above generation methods, the axicon is the most commonly used way due to its simple structure and high conversion efficiency. However, the Bessel beam generated by a single positive axicon can only propagate a few centimeters [30], because the transmitted conical waves have a relatively large conical angle due to fabrication limitation. Wu et al. proposed an axicon doublet, which was composed of a positive axicon and a negative axicon with the same base angles but different materials [31]. Numerical simulations in Ref. [31] reported that the generated Bessel beam could propagate tens of meters. Koronkevich et al. fabricated a diffractive optical element, which transformed a point light source into an axial focused line with propagation distance more than 100 meters [32]. Axial intensities of the above generated Bessel beam are linearly increased with prominent oscillations [22,30–33]. Until now, there is rare research on the generation of uniform Bessel beams, which hampers its applications to stable imaging and detection. Therefore, we propose a new amplitude filter for flattening the axial intensity distribution. On this basis, we put forward a new optical system, which is composed of the axicon doublet and the newly proposed amplitude filter. It is expected that the generated Bessel beam should propagate stably within a long axial range.

This paper is organized as follows. In Section 2, we describe the design principle and structure of the proposed new optical system, including the axicon doublet and the amplitude filter. Formulas for calculating the Bessel beam intensities and evaluating stable performance indices of the Bessel beam are written in detail. In Section 3, we simulate the diffracted Bessel beam intensities by the complete Rayleigh-Sommerfeld method. Numerical results are physically explained. A brief conclusion is drawn in Section 4 together with some discussions.

## 2. Design principles of the axicon doublet and the amplitude filter together with calculating formulas of the generated Bessel beam

#### 2.1 Design of the axicon doublet

For the traditional case of a single positive axicon, as shown in Fig. 1(a), the incident plane wave is refracted by the axicon to form a conical wave. According to Snell’s law, the conical angle ${\theta _1}$ is given by

where $\alpha $ is the base angle of the axicon; ${n_1}$ and ${n_0}$ are the refractive indices of the axicon and the surrounding region, respectively, as shown in Fig. 1(a). Coherent superposition of the conical waves produces a Bessel beam in the transmitted region behind the axicon. The longest propagating distance can be approximately obtained by [1,21,34] where*R*is the radius of the axicon, as shown in Fig. 1(a).

We can conclude from Eq. (2) that the longest propagation distance relies on the axicon radius *R* and the conical angle ${\theta _1}$. We make a typical approximation. If ${n_1} = 1.5$, ${n_0} = 1.0$, $R = 5\; \textrm{cm}$, and $\alpha = 25^\circ $, we can obtain from Eq. (2) that $z_{\textrm{max}}^\textrm{s} = 19.56\; \textrm{cm}$. However, if the positive axicon is applied to far distance focusing, the base angle $\alpha $ should be very small. For instance, if $z_{\textrm{max}}^\textrm{s} \ge 100\; \textrm{m}$, from Eqs. (1) and (2), we can get $\alpha \le 0.05^\circ $, which is much less than the experimental fabrication error.

To reduce the conical angle of the transmitted waves and enlarge the propagation distance, the axicon doublet is used [31], as shown in Fig. 1(b). The axicon doublet consists of a positive axicon and a negative axicon with the same base angles. The double axicons are perfectly attached together and no space is left between the two conical surfaces. When a plane wave normally impinges on the front surface of the axicon doublet, it will go straightly into the positive axicon. On the conical surface of the positive axicon, the refractive wave will deviate from the original propagation direction owing to a refractive index difference between the double axicons. By analogy with Eq. (1), the deviation angle $\delta $ is calculated as

where ${n_1}$ and ${n_2}$ are the refractive indices of the positive axicon and the negative axicon, respectively; $\alpha $ is the base angle of the positive axicon, as shown in Fig. 1(b). Then, the propagating wave will go to the exit surface of the negative axicon and refract twice to the outside space. The outgoing wave behind the axicon doublet has a conical angle ${\theta _2}$ of where ${n_0}$ is the refractive index of the outside space, as seen in Fig. 1(b). If the base angle $\alpha $ is small, Eqs. (3) and (4) can be simplified. Then the conical angle ${\theta _2}$ is approximated asIt is seen from Eq. (5) that a very tiny conical angle ${\theta _2}$ is obtained as long as the refractive index difference between the double axicons is small enough. Consequently, the maximum propagating distance of the Bessel beam and the radius of the axicon doublet are related by

where*R*is the radius of the axicon doublet. If ${\theta _2}$ is very small, the maximum propagation distance $z_{\textrm{max}}^{\textrm{ad}}$ in Eq. (6a) will be a considerably large value in the case of the axicon doublet. We use the same parameters chosen above to make an estimation. When $R = 5\; \textrm{cm}$, ${n_1} = 1.5$, ${n_2} = 1.499$, ${n_0} = 1.0$, and $\alpha = 25^\circ $, we can get $z_{\textrm{max}}^{\textrm{ad}} = 107.22\; \textrm{m}.$

In Fig. 1(b), since the outgoing waves are conical waves with deviation angle ${\theta _2}$ from the optical axis, the transmitted field in the $z = 0$ plane is

where $A(\rho )$ represents the amplitude distribution; $\rho $ is the radial coordinate in the $z = 0$ plane, as shown in Fig. 1(b); $j = \sqrt { - 1} $ is the complex unit. If the incident light is a plane wave with unit amplitude and we ignore the reflection losses on all the three interfaces, then $A(\rho )= 1$ is valid.#### 2.2 Design of the amplitude filter

In this subsection, we propose and design a new amplitude filter to flatten the axial intensity distribution of the generated Bessel beam. It can be divided into two steps, in which we suppress the axial intensity increasing trend and oscillating effect successively.

To examine the reason for the axial intensity increasing trend, we draw a ray tracing picture of the axicon doublet in Fig. 2(a), which is similar to the case of a single positive axicon. In Fig. 2(a), the plane waves are incident upon the front surface of the axicon doublet, after twice refraction, the transmitted waves have the same conical angles. Accordingly, the intensity at a definite axial position *z* depends on the field superposition of all the source points on a circumference with radius $\rho = z \times \tan ({{\theta_2}} )$ in the $z = 0$ plane. By using limit principle, the incident power within an infinitesimal annular region $[{\rho ,\; \rho + \textrm{d}\rho } ]$ is assumed to be concentrated in an infinitesimal axial region $[{z,\; z + \textrm{d}z} ]$, as shown in Fig. 2(a). When the incident light is a plane wave with unit amplitude, the incident power within an infinitesimal annular region $[{\rho ,\; \rho + \textrm{d}\rho } ]$ is calculated as $2\mathrm{\pi }\rho \textrm{d}\rho $, proportional to the annular region radius $\rho $. Therefore, the axial intensity $I(z )$ is proportional to the corresponding annular region radius $\rho $. Taking the relation $\rho = z \times \tan ({{\theta_2}} )$ into account, we now understand why the axial intensity $I(z )$ is almost linearly enlarged with the increase of the axial position *z*.

To suppress the increasing axial intensity into a stable one, a simple way is to impose an amplitude filter in the $z = 0$ plane whose transmission coefficient is proportional to $1/\sqrt \rho .$ For avoiding singularity at the origin point $\rho = 0,$ we assume that the incident light is completely blocked for $\rho < {R_1}$. If we hope to generate a Bessel beam within the axial range $[{{z_1},\; {z_2}} ]$, from Eq. (6b) the amplitude filter must be transparent in the annular regime $[{{R_1},\; {R_2}} ]$ with ${R_1} = {z_1} \times \tan ({{\theta_2}} )$ and ${R_2} = {z_2} \times \tan ({{\theta_2}} )$. The transmission coefficient of the annular amplitude filter is given by

*R*is the radius of the axicon doublet, as shown in Fig. 2(b); ${R_1}$ and ${R_2}$ satisfy $0 < {R_1} < {R_2} \le R$. In contrast, for the case of the axicon doublet without amplitude filter, the transmission coefficient is ${T_0}(\rho )= 1$ for $\rho \le R$ and ${T_0}(\rho )= 0$ for $\rho > R$. It is noted that in Ref. [35] Davidson et al. used a similar formula like Eq. (8) to inhibit the axial intensity increasing trend.

In the second step, we need to suppress the axial intensity oscillations, which are originated from field diffraction effect at the aperture edge where an abrupt cutoff of the incident field is occurred. A common way to relieve the diffraction effect is introducing a gradually absorbing layer near the aperture edge [20,26,34]. Axial intensity oscillations can also be suppressed with a heart-shaped binary amplitude mask, in which the filling ratio of the transparent region is exactly the same as the amplitude distribution of the gradually absorbing layer [36–39]. Motivated by their works, here we apply a gradually varying amplitude mask to stifling axial intensity oscillations of the Bessel beam generated by the axicon doublet. The transmission coefficient of the gradually varying amplitude mask is written as

Through synthesizing the two transmission coefficient functions in Eqs. (8) and (9), we can expect a flat axial intensity of the generated Bessel beam. Hence, we propose such a smoothing amplitude filter whose transmission coefficient is

#### 2.3 Formulas for calculating the Bessel beam intensity and evaluating the stable propagation performance indices

Once the axicon doublet and the smoothing amplitude filter are designed, the diffracted field ${U_2}({x,y,z} )$ of the generated Bessel beam is calculated by the complete Rayleigh-Sommerfeld method as [40]

*z*is the longitudinal coordinate; $T({{x_0},{y_0}} )$ is the transmission coefficient in the $z = 0$ plane, which can be calculated from Eq. (10) with $\rho = \sqrt {x_0^2 + y_0^2} $; ${x_0}$ and ${y_0}$ respectively represent the transverse coordinates along the $x$-axis and the $y$-axis in the $z = 0$ plane; ${U_1}({{x_0},{y_0},z = 0} )$ is the incident field in the $z = 0$ plane, whose formula is given by Eq. (7) with $\rho = \sqrt {x_0^2 + y_0^2} $; $r = \sqrt {{{({x - {x_0}} )}^2} + {{({y - {y_0}} )}^2} + {z^2}} $ represents the distance between the source point $({{x_0},{y_0},z = 0} )$ and an arbitrary observation point $({x,y,z} )$. The diffracted Bessel beam intensity is ${I_2}({x,y,z} )= {|{{U_2}({x,y,z} )} |^2}$, where $|{{U_2}} |$ represents the magnitude of the complex field ${U_2}$ calculated from Eq. (11).

To quantitatively characterize the propagation invariance property of the generated Bessel beam, the axial intensity uniformity and coincidence of the cross-sectional intensity profiles are two important performance indices. For appraising the axial intensity uniformity within the longitudinal range $[{{z_\textrm{a}}\; {z_\textrm{b}}} ]$, we define an axial intensity percentage error (AIPE) function as

For judging the coincidence of the cross-sectional intensity profiles within the central lobe, we define the transverse intensity percentage error (TIPE) function as

*M*indicates the number of the calculated cross-sectional planes; $\overline {{I_x}} ({x,0} )$ is the average transverse intensity along the $x$-axis for all the above

*M*planes, which is given by $\overline {{I_x}} ({x,0} )= \mathop \sum \limits_{m = 1}^M {I_x}({x,0,{z_m}} )/M$; the integral lower limit $- {d_m}/2$ and upper limit ${d_m}/2$ respectively assign the first minimum intensity position on both sides of the central spot in the cross-sectional plane at $z = {z_m}$; ${\delta _\textrm{t}}$ represents the transverse intensity percentage error for the above

*M*planes.

## 3. Performance analysis of the Bessel beam generated by the optical system composed of the axicon doublet and the smoothing amplitude filter

In this section, we choose a set of parameters to simulate the propagating properties of the Bessel beam generated by the proposed optical system consisting of the axicon doublet and the smoothing amplitude filter. The propagating properties include the axial intensity distribution along the optical axis and the transverse intensity distribution on several cross-sectional planes. Parameters are selected as follows. The radius and base angle of the axicon doublet are respectively assigned as $R = 10\; \textrm{cm}$ and $\alpha = 35^\circ $, as shown in Fig. 2(b); the incident plane wave has a unit amplitude and a wavelength of $\lambda = 532\; \textrm{nm}$ in vacuum; the positive and negative axicons are made of H51 and H50 glasses from CDGM company, with refractive indices being 1.5260 and 1.5254 for $\lambda = 532\; \textrm{nm}$, respectively. On setting the axial propagation range of the Bessel beam to [30 210] m, we can calculate the annular transparent region radii as [12.60 88.23] mm from Eq. (6b). The smoothing parameter is chosen as $\varepsilon = 0.15$ in Eq. (10).

#### 3.1 Axial intensity of the long-distance stably propagating Bessel beam

With all the parameters selected above, we can calculate the axial intensity distribution of the Bessel beam from Eq. (11) on setting $x = y = 0$, as shown in Fig. 3. Figures 3(a), 3(b), and 3(c) correspond to the axial intensity distributions of the Bessel beams generated by the axicon doublets without amplitude filter, with the annular amplitude filter, and with the smoothing amplitude filter, respectively. Their transmission coefficients are plotted in Figs. 3(d), 3(e), and 3(f), respectively. For the axicon doublet without amplitude filter, the axial intensity is linearly increased with some oscillations, as seen in Fig. 3(a). Through putting an annular amplitude filter behind the axicon doublet, the axial intensity climbing effect is substantially suppressed while prominent oscillations still exist in the axial intensity profile, as shown in Fig. 3(b). Furthermore, if we attach a smoothing amplitude filter to the axicon doublet, we have obtained a considerably stable axial intensity distribution of the generated Bessel beam during a long axial range, as shown in Fig. 3(c). Numerical results reveal that within the axial range [80 140] m the maximum AIPEs are 50.50%, 12.52%, and 0.22% for Figs. 3(a), 3(b), and 3(c), respectively. If we expand the axial range to [70 150] m, the maximum AIPEs are 65.40%, 14.51%, and 0.45% for Figs. 3(a), 3(b), and 3(c), respectively.

To demonstrate the long-distance propagating property of the Bessel beam generated by the axicon doublet, the axial intensity distribution of a single positive axicon is also investigated for comparison. The single positive axicon is made of H51 glass. Other parameters are the same as above. Numerical results show that the maximum propagation distance for the single positive axicon is only 20.37 cm; in contrast, the Bessel beam generated by the axicon doublet can propagate as far as more than 200 meters. So far we have clearly proved both long distance and axial stable propagation properties of the Bessel beam generated by the proposed optical system.

#### 3.2 Transverse intensity of the long-distance stably propagating Bessel beam

In addition, the cross-sectional intensity patterns are studied to see the central lobe stableness of the generated Bessel beam. To this end, we select four cross-sectional planes with longitudinal positions at ${z_1} = 80\; \textrm{m}$, ${z_2} = 100\; \textrm{m}$, ${z_3} = 120\; \textrm{m}$, and ${z_4} = 140\; \textrm{m}$. For the axicon doublet without amplitude filter, the cross-sectional intensity patterns at ${z_1} = 80\; \textrm{m}$, ${z_2} = 100\; \textrm{m}$, ${z_3} = 120\; \textrm{m}$, and ${z_4} = 140\; \textrm{m}$ are displayed in Figs. 4(a), 4(b), 4(c), and 4(d), respectively. It is seen from Figs. 4(a) to 4(d) that the central spots almost have identical sizes, but the peak intensity increases about twice, as shown by the colorbars. For the axicon doublet with the annular amplitude filter, the intensity patterns in above four cross-sectional planes are drawn in Figs. 4(e), 4(f), 4(g), and 4(h), respectively. The peak intensity differs much less from the colorbars of Figs. 4(e) to 4(h), and the four central spot sizes still keep unchanged. For the axicon doublet with the smoothing amplitude filter, the corresponding four intensity patterns are shown in Figs. 4(i), 4(j), 4(k), and 4(l), respectively. From Figs. 4(i) to 4(l), not only the four central spots have identical sizes, but also their peak intensities almost have no difference.

To see the transverse intensity profiles more accurately, we extract the line-scan intensities along the $x$-axis from Figs. 4(a) to 4(d) and plot them together in Fig. 5(a). It is exhibited in Fig. 5(a) that the peak intensity steadily increases with the increase of the propagation distance. Numerical results reveal that the central spot radii in Fig. 5(a) are 481, 485, 481, and 476 µm for the four transverse planes located at ${z_1} = 80\; \textrm{m}$, ${z_2} = 100\; \textrm{m}$, ${z_3} = 120\; \textrm{m}$, and ${z_4} = 140\; \textrm{m}$, respectively. The TIPE is as large as 83.66%. Figures 5(b) and 5(c) are the same as Fig. 5(a) except for the axicon doublets with the annular amplitude filter and with the smoothing amplitude filter, respectively. A main difference between Fig. 5(b) and 5(a) is that the axial intensity increasing effect is significantly depressed and the peak intensity difference is now confined in a relatively small range. Numerical results indicate that the central spot sizes in Fig. 5(b) are 486, 486, 485, and 498 µm for the four consecutive planes. The TIPE is decreased to 12.44%. In Fig. 5(c), all the four intensity profiles overlap with one another, which suggests that the generated Bessel beam propagates stably along the axial direction. Numerical results reveal that the four central spot sizes in Fig. 5(c) are all identical as 485 µm, and the TIPE is as small as 0.33%. It is concluded that the propagating Bessel beam maintains a very stable transverse intensity profile with the use of the smoothing amplitude filter, due to a constant central spot size and high degree of transverse intensity coincidence.

## 4. Conclusion and discussions

In this paper, we propose a new optical system consisting of an existing axicon doublet and a newly proposed amplitude filter. The axicon doublet extends the propagating distance of the generated Bessel beam up to around three orders of magnitude in comparison with the single positive axicon, owing to a very small refractive index difference between the double axicons. The proposed smoothing amplitude filter precisely manipulates the axial intensity profile so that a stable Bessel beam is produced during the whole propagation range. The generated long-distance stably propagating Bessel beam is expected to have practical applications for far-distance focusing and imaging. Higher order stable Bessel beams with long propagation distance can be generated through inserting a vortex phase plate in the input plane. It is worth to mention that the Bessel beam propagation distance may be furtherly extended to several kilometers through increasing the axicon doublet diameter, decreasing the base angle of the axicon doublet, or reducing the refractive index difference between the double axicons.

Since the proposed axicon doublet can produce very small angle conical waves, it can also be used to analyze parallelism of light. In addition, the maximum propagation distance of the Bessel beam is sensitive to the refractive index difference between the double axicons, so the axicon doublet may also be applied for measuring material refractive index with high precision. Dispersive property investigation of the axicon doublet may open up new application opportunities. For instance, if the double axicon materials are specially chosen to enlarge optical dispersion, the axicon doublet can serve as a color separating and focusing element for optical display.

## Funding

National Natural Science Foundation of China (61735002, 11774243, 11774246, 11474206, 11404224); Youth Innovative Research Team of Capital Normal University (008/20530290053, 008/19530050146, 008/18530500155); Connotative Development Foundation for Distinguished Young Talents in Capital Normal University (2055105); Capacity Building for Sci-Tech Innovation-Fundamental Scientific Research Funds (008/20530290072, 008/19530050170, 008/19530050180, 008/18530500186, 025185305000/142).

## Disclosures

The authors declare no conflicts of interest.

## References

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

**2. **C. Vetter, R. Steinkopf, K. Bergner, M. Ornigotti, S. Nolte, H. Gross, and A. Szameit, “Realization of free-space long-distance self-healing Bessel beams,” Laser Photonics Rev. **13**(10), 1900103 (2019). [CrossRef]

**3. **J. Zhao, P. Zhang, D. Deng, J. Liu, Y. Gao, I. D. Chremmos, N. K. Efremidis, D. N. Christodoulides, and Z. Chen, “Observation of self-accelerating Bessel-like optical beams along arbitrary trajectories,” Opt. Lett. **38**(4), 498–500 (2013). [CrossRef]

**4. **G. Häusler and W. Heckel, “Light sectioning with large depth and high resolution,” Appl. Opt. **27**(24), 5165–5169 (1988). [CrossRef]

**5. **S. Mitra, M. Chanal, R. Clady, A. Mouskeftaras, and D. Grojo, “Millijoule femtosecond micro-Bessel beams for ultra-high aspect ratio machining,” Appl. Opt. **54**(24), 7358–7365 (2015). [CrossRef]

**6. **A. O. Ambriz, S. L. Aguayo, Y. V. Kartashov, V. A. Vysloukh, D. Petrov, H. G. Gracia, J. C. G. Vega, and L. Torner, “Generation of arbitrary complex quasi-non-diffracting optical patterns,” Opt. Express **21**(19), 22221–22231 (2013). [CrossRef]

**7. **K. Dholakia and T. Čižmár, “Shaping the future of manipulation,” Nat. Photonics **5**(6), 335–342 (2011). [CrossRef]

**8. **F. O. Fahrbach and A. Rohrbach, “Propagation stability of self-reconstructing Bessel beams enables contrast-enhanced imaging in thick media,” Nat. Commun. **3**(1), 632 (2012). [CrossRef]

**9. **F. O. Fahrbach, P. Simon, and A. Rohrbach, “Microscopy with self-reconstructing beams,” Nat. Photonics **4**(11), 780–785 (2010). [CrossRef]

**10. **F. O. Fahrbach and A. Rohrbach, “A line scanned light-sheet microscope with phase shaped self-reconstructing beams,” Opt. Express **18**(23), 24229–24244 (2010). [CrossRef]

**11. **F. O. Fahrbach, V. Gurchenkov, K. Alessandri, P. Nassoy, and A. Rohrbach, “Self-reconstructing sectioned Bessel beams offer submicron optical sectioning for large fields of view in light-sheet microscopy,” Opt. Express **21**(9), 11425–11440 (2013). [CrossRef]

**12. **A. Chong, W. H. Renninger, D. N. Christodoulides, and F. W. Wise, “Airy–Bessel wave packets as versatile linear light bullets,” Nat. Photonics **4**(2), 103–106 (2010). [CrossRef]

**13. **M. McLaren, M. Agnew, J. Leach, F. S. Roux, M. J. Padgett, R. W. Boyd, and A. Forbes, “Entangled Bessel-Gaussian beams,” Opt. Express **20**(21), 23589–23597 (2012). [CrossRef]

**14. **A. Badalyan, R. Hovsepyan, P. Mantashyan, V. Mekhitaryan, and R. Drampyana, “Nondestructive readout of holograms recorded by Bessel beam technique in LiNbO_{3}:Fe and LiNbO_{3}:Fe:Cu crystals,” Eur. Phys. J. D **68**(4), 82 (2014). [CrossRef]

**15. **Y. Yuan, T. Lei, Z. Li, Y. Li, S. Gao, Z. Xie, and X. Yuan, “Beam wander relieved orbital angular momentum communication in turbulent atmosphere using Bessel beams,” Sci. Rep. **7**(1), 42276 (2017). [CrossRef]

**16. **C. Chen and S. Huang, “The generation of non-conventional beam in a nonlinear electro-optic photonic crystal,” Jpn. J. Appl. Phys. **59**(9), 092003 (2020). [CrossRef]

**17. **Y. Fan, B. Cluzel, M. Petit, X. L. Roux, A. Lupu, and A. Lustrac, “2D waveguided Bessel beam generated using integrated metasurface-based plasmonic axicon,” ACS Appl. Mater. Interfaces **12**(18), 21114–21119 (2020). [CrossRef]

**18. **R. Chriki, G. Barach, C. Tradosnky, S. Smartsev, V. Pal, A. A. Friesem, and N. Davidson, “Rapid and efficient formation of propagation invariant shaped laser beams,” Opt. Express **26**(4), 4431–4439 (2018). [CrossRef]

**19. **Z. Jiang, Q. Lu, and Z. Liu, “Propagation of apertured Bessel beams,” Appl. Opt. **34**(31), 7183–7185 (1995). [CrossRef]

**20. **R. Borghi, M. Santarsiero, and F. Gori, “Axial intensity of apertured Bessel beams,” J. Opt. Soc. Am. A **14**(1), 23–26 (1997). [CrossRef]

**21. **J. Durnin, J. J. Miceli, and J. H. Eberly, “Diffraction-free beams,” Phys. Rev. Lett. **58**(15), 1499–1501 (1987). [CrossRef]

**22. **A. Vasara, J. Turunen, and A. T. Friberg, “Realization of general nondiffracting beams with computer-generated holograms,” J. Opt. Soc. Am. A **6**(11), 1748–1754 (1989). [CrossRef]

**23. **J. Turunen, A. Vasara, and A. T. Friberg, “Holographic generation of diffraction-free beams,” Appl. Opt. **27**(19), 3959–3962 (1988). [CrossRef]

**24. **A. J. Asuncion and R. A. Guerrero, “Generating superimposed Bessel beams with a volume holographic axicon,” Appl. Opt. **56**(14), 4206–4212 (2017). [CrossRef]

**25. **J. H. Mcleod, “The axicon: a new type of optical element,” J. Opt. Soc. Am. **44**(8), 592–597 (1954). [CrossRef]

**26. **S. Y. Popov and A. T. Friberg, “Apodization of generalized axicons to produce uniform axial line images,” Pure Appl. Opt. **7**(3), 537–548 (1998). [CrossRef]

**27. **R. M. Herman and T. A. Wiggins, “Production and uses of diffractionless beams,” J. Opt. Soc. Am. A **8**(6), 932–942 (1991). [CrossRef]

**28. **A. J. Cox and D. C. Dibble, “Nondiffracting beam from a spatially filtered Fabry–Perot resonator,” J. Opt. Soc. Am. A **9**(2), 282–286 (1992). [CrossRef]

**29. **W. T. Chen, M. Khorasaninejad, A. Y. Zhu, J. Oh, R. C. Devlin, A. Zaidi, and F. Capasso, “Generation of wavelength-independent subwavelength Bessel beams using metasurfaces,” Light: Sci. Appl. **6**, e16259 (2017). [CrossRef]

**30. **R. Dharmavarapu, S. Bhattacharya, and S. Juodkazis, “Diffractive optics for axial intensity shaping of Bessel beams,” J. Opt. **20**(8), 085606 (2018). [CrossRef]

**31. **F. T. Wu, Q. A. Zhang, and W. T. Zheng, “Generating long-distance nondiffracting Bessel beams with equivalent axicon,” Chinese Lasers **38**(12), 1202004 (2011, in Chinese). [CrossRef]

**32. **V. P. Koronkevich, I. A. Mikhaltsova, E. G. Churin, and Y. I. Yurlov, “Lensacon,” Appl. Opt. **34**(25), 5761–5772 (1995). [CrossRef]

**33. **R. Li, X. Yu, T. Peng, Y. Yang, B. Yao, C. Zhang, and T. Ye, “Shaping the on-axis intensity profile of generalized Bessel beams by iterative optimization methods,” J. Opt. **20**(8), 095603 (2018). [CrossRef]

**34. **A. J. Cox and J. D’Anna, “Constant-axial-intensity nondiffracting beam,” Opt. Lett. **17**(4), 232–234 (1992). [CrossRef]

**35. **N. Davidson, A. A. Friesem, and E. Hasman, “Efficient formation of nondiffracting beams with uniform intensity along the propagation direction,” Opt. Commun. **88**(4-6), 326–330 (1992). [CrossRef]

**36. **J. S. Ye, L. J. Xie, X. K. Wang, S. F. Feng, W. F. Sun, and Y. Zhang, “Flattening axial intensity oscillations of a diffracted Bessel beam through a cardioid-like hole,” Opt. Express **26**(2), 1530–1537 (2018). [CrossRef]

**37. **R. Borghi, “Heart diffraction,” Opt. Lett. **42**(11), 2070–2073 (2017). [CrossRef]

**38. **Y. Wang, L. Xie, J. Ye, W. Sun, X. Wang, S. Feng, P. Han, Q. Kan, and Y. Zhang, “Tailoring axial intensity of laser beams with a heart-shaped hole,” Opt. Lett. **42**(23), 4921–4924 (2017). [CrossRef]

**39. **R. Borghi, “‘Tailoring axial intensity of laser beams with a heart-shaped hole,’ by Wang et al.: comment,” Opt. Lett. **43**(14), 3240 (2018). [CrossRef]

**40. **G. D. Gillen and S. Guha, “Modeling and propagation of near-field diffraction patterns: A more complete approach,” Am. J. Phys. **72**(9), 1195–1201 (2004). [CrossRef]